Abstract

By introducing hybrid technique into high-order CPR (correction procedure via reconstruction) scheme, a novel hybrid WCNS-CPR scheme is developed for efficient supersonic simulations. Firstly, a shock detector based on nonlinear weights is used to identify grid cells with high gradients or discontinuities throughout the whole flow field. 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.

1. 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 [16].

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, 1012]. 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 -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 shock-capturing 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 shock-capturing 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 shock-capturing 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.

2. Brief Review of WCNS and CPR

In this section, a brief review of WCNS and CPR is presented based on two-dimensional hyperbolic conservation lawswhere is a conservative variable vector and and are the flux vectors.

2.1. Brief Review of Third-Order WCNS

Equation (1) is discretized by a finite difference scheme at every node where is the approximation of the spatial derivatives. The following fourth-order flux difference operators are usedwhere 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 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 into two smaller substencils and .

Step 2. Construct an interpolation polynomial in each substencil , which satisfies . We have

Step 3. Calculate the nonlinear weights based on linear weights and a smoothness indicator for each substencil . The linear weights are

The improved Z-type nonlinear weights [24, 25] are defined bywhere is a small number, is a smoothness indicator, , and and are parameters. We use and to improve the performance of the nonlinear weights.

An improved smoothness indicator [25] is employed.

Step 4. Get weighted interpolation value at .The formulations of is symmetric with by .

2.2. Review of Third-Order CPR

For the third-order CPR framework, there are solution points , in the -th cell, as shown in Figure 2. The nodal values of the conservative variable at solution points are updated by the following differential equation

The reconstructed flux and are expressed aswhere and are flux correction polynomials

Here, correction functions , , , and are all cubic polynomials, and and are Riemann fluxes corresponding to and , respectively.

Moreover, the solution polynomial iswhere and 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 for Correction function for scheme is used, which is expressed as where is Legendre polynomials. More specifically, and

3. 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 [2628]. 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 represents the number of candidate stencils in WCNS with 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.

3.2. 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 and the solution points of WCNS by .

3.2.1. 1D Interface Treatment

First, WCNS in the -th cell needs information of ghost points from the neighboring -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 . That is,

Second, we apply provided by -th CPR cell and computed in -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 -th cell is a smooth cell at the -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 -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 -th solution point as an example.

Step 1. Divide the stencil into two substencils .

Step 2. Construct interpolation polynomials and in and , respectively, we havewhere and are linear Lagrange polynomials.

Step 3. Calculate the nonlinear weights based on linear weights and . The linear weights can be obtained by solving the following equations:Then, nonlinear weights can be calculated by equation (6).

Step 4. Get the interpolation value .

3.2.2. 2D Interface Treatment

In this subsection, we consider the interface treatment between -th CPR cell and -th WCNS cell in 2D computational domain in Figure 5(a).

Similarly, WCNS needs information of ghost points from the -th CPR cell. The solution values at ghost points of WCNS in the CPR cell are calculated via CPR interpolation using data at . Using dimension by dimension interpolation strategy, we do the Lagrange interpolation in the direction first, as shown in Figure 5(b), and then in the direction, as shown in Figure 5(c). Using solution polynomial equation (15) of CPR, we haveat the ghost points .

Then, we calculate the numerical fluxes at cell interface. As shown in Figure 6, two types of flux points, and , do not coincide with each other. Thus, and need to be interpolated before calculating numerical fluxes. is obtained by CPR linear interpolation using equation (15) based on values in -th cell, while is obtained through WCNS nonlinear interpolation dimension by dimension. The interpolated values of and are denoted by and correspondingly.

Therefore, for the CPR scheme, the numerical flux at is

Meanwhile, the numerical flux at for WCNS is

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.

4. Numerical Investigation

In this section, we will test spatial accuracy, shock-capturing 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 error and error are computed to measure local solution quality and global solution quality, respectively. The error is estimated bywhere and are the numerical solution and the exact solution at the -th solution point. Here, is the total number of solution points.

4.1. 1D Linear Wave Equation

Consider the 1D linear wave equationwith the initial conditionand periodic boundary conditions. and time step is used to solve the problem till 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 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.

4.2. 1D Euler Equations

Consider 1D Euler equationswherewhere is the density, is the velocity, is the total energy, and is the pressure, which is related to the total energy by with . Three kinds of shock-tube problems are considered in this subsection.

4.2.1. Shock-Tube Problem of Sod

Consider the Sod problem with initial conditions [29]and Dirichlet boundary conditions. Time step , , and are used to solve the problem till by the hybrid WCNS-CPR scheme. The reference solution is obtained by the WCNS3 scheme with . As shown in Figure 8, shocks and contact discontinuities can be well captured. Near the contact discontinuity () and the shock (), 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.

4.2.2. Shu-Osher Problem

Consider the Shu-Osher problem with initial conditions [30]and Dirichlet boundary conditions. Time step , , and are used to solve the problem till by the hybrid WCNS-CPR scheme.

The reference solution is obtained by the WCNS3 scheme with . 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.

4.2.3. Shock-Tube Problem of Lax

Consider the Lax problem with initial condition [31]and Dirichlet boundary conditions. Time step , , and are adopted to solve the problem till by the hybrid WCNS-CPR scheme. The reference solution is obtained by the WCNS3 scheme with . We know from Figure 10 that at and , shocks and discontinuities are all captured without obvious oscillations.

4.3. 2D Euler Equations

Let us consider the 2D Euler equations of gas dynamics: For an ideal gas, . Here, ,,,, and are the density, the direction velocity, the direction velocity, the pressure, and the total energy, respectively.

4.3.1. 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, , with the isotropic vortex perturbations added to the mean flow. There are perturbations in and and no perturbation in entropy :where and the vortex strength . The computational domain is with periodic boundary conditions. The problem is solved till with time step for accuracy test. The parameter is in this test.

Figure 11 shows distribution of different cells, NI in the direction (NI in the 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 , 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 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.

Remark 2. WCNS represents WCNS with characteristic projection.

4.3.2. 2D Riemann Problem

We solve the Riemann problem for equation (33) on the computational domain of with initial conditionsand Dirichlet boundary conditions until . The parameter is in this test.

Compared with the results of WCNS in 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.

4.3.3. 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 on the lower boundary, and its direction is 60° from the -axis. Inflow and slip wall boundary conditions are specified at the bottom boundary for and , respectively. For the upper boundary (), a time-dependent boundary based on the analytical propagation speed of the oblique shock is imposed.

The initial condition is

The simulation is run until using time step .

As shown in Figures 15(a) and 15(b), essential small- and large-scale structures of the flow field are well captured compared with the results of WCNS in 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 and , the hybrid scheme is slightly more efficient than WCNS. However, for and , the hybrid time is about 1.6 times more efficient than the WCNS scheme.

5. 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.

Acknowledgments

This study was supported by the National Natural Science Foundation of China (Grant Nos. 11902344 and 11971481), the Basic Research Foundation of National Numerical Wind Tunnel Project (No. NNW2018-ZT4A08), and the foundation of State Key Laboratory of Aerodynamics (Grant No. SKLA2019010101). The authors would like to thank assistant researcher Chen Li from the State Key Laboratory of Aerodynamics for his valuable suggestions.