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

Authors: Jia Guo; Huajun Zhu; Zhen-Guo Yan; Lingyan Tang; Songhe Song
Journal: International Journal of Aerospace Engineering
 (2020)

---

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


---

## Body

## 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 [1–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–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 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 Section2, 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 laws(1)Ut+FUx+GUy=0,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 xi,yj(2)Uti,j=Ei,j′=−Fxi,j−Gyi,j,where Ei,j′ is the approximation of the spatial derivatives. The following fourth-order flux difference operators are used(3)Fxi,j=43hF̂i+1/2,j−F̂i−1/2,j−16hFi+1,j−Fi−1,j,Gyi,j=43hĜi,j+1/2−Ĝi,j−1/2−16hGi,j+1−Gi,j−1,where h is the grid spacing. Here, F̂i−1/2,j=F̂i−1/2Ui−1/2,jL,Ui−1/2,jR,Ĝi,j−1/2=Ĝj−1/2Ui,j−1/2L,Ui,j−1/2R 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.Figure 1
Distribution of solution points and flux points of 2D WCNS ini,j-th cell. Solution points are marked in blue dot and flux points in orange square.To capture discontinuities, the conservative variables at cell edgesUi−1/2,jR 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 stencilSi=Ui−1,j,Ui,j,Ui+1,j into two smaller substencils Si1=Ui−1,j,Ui,j and Si2=Ui,j,Ui+1,j.Step 2.
Construct an interpolation polynomialpmx in each substencil Sim,m=1,2, which satisfies pmxi=Ui,j. We have(4)p1xi−1/2=12Ui−1,j+12Ui,j,p2xi−1/2=32Ui,j−12Ui+1,j.Step 3.
Calculate the nonlinear weightsωk based on linear weights dm and a smoothness indicator ISk for each substencil Sim,m=1,2. The linear weights are(5)d1=34,d2=14.The improved Z-type nonlinear weights [24, 25] are defined by(6)ωk=βk∑m=12βm,(7)βk=dkCq+τ5ε+ISkpz,(8)k=1,2,where ε=10−6 is a small number, ISk is a smoothness indicator, τ5=IS0−IS1, and Cq and pz are parameters. We use pz=2 and Cq=4 to improve the performance of the nonlinear weights.An improved smoothness indicator [25] is employed.(9)IS0=143Ui,j−4Ui−1,j+Ui−2,j2+Ui,j−2Ui−1,j+Ui−2,j2,IS1=143Ui+2,j−4Ui+1,j+3Ui,j2+Ui+2,j−2Ui+1,j+Ui,j2.Step 4.
Get weighted interpolation value atxi−1/2.(10)Ui−1/2,jR=∑m=12ωmpmxi−1/2.
The formulations ofUi−1/2,jL is symmetric with Ui−1/2,jR by i.


### 2.2. Review of Third-Order CPR

For the third-order CPR framework, there are32 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(11)∂Ui,j∂t+∂F̃i,jξ,η∂ξ+∂G̃i,jξ,η∂η=0,ξ,η∈−1,1×−1,1.Figure 2
Distribution of solution points and flux points of 2D CPR ini,j-th cell. Solution points are marked in green dot and flux points on element boundary in purple square. The flux points inside the CPR cell coincide with the solution points.The reconstructed fluxF̃i,jξ,η and G̃i,jξ,η are expressed as(12)F̃i,jξ,η=Fi,jξ,η+σi,j1ξ,η,G̃i,jξ,η=Gi,jξ,η+σi,j2ξ,η,where σi,j1ξ,η and σi,j2ξ,η are flux correction polynomials(13)σi,j1ξl,ηm=F̂−1,ηm−Fi,j−1,ηmgLξl+F̂1,ηm−Fi,j1,ηmgRξl,σi,j2ξl,ηm=Ĝξl,−1−Gi,jξl,−1gLηm+Ĝξl,1−Gi,jξl,1gRηm.Here, correction functionsgLξ, gRξ, gLη, and gRη are all cubic polynomials, and F̂ and Ĝ are Riemann fluxes corresponding to F and G, respectively.(14)F̂−1,η=F̂Ui−1,j1,η,Ui,j−1,η,F̂1,η=F̂Ui,j1,η,Ui+1,j−1,η.Moreover, the solution polynomial is(15)Ui,jξ,η=∑l=13∑m=13Ui,j,l,mLlξLmηwhere Llξ and Lmη are the 1D Lagrange polynomials in the ξ and η directions.The flux polynomial is(16)Fi,jξ,η=∑l=13∑m=13Fi,j,l,mLlξ·Lmη,Gi,jξ,η=∑l=13∑m=13Gi,j,l,mLlξ·Lmη.Remark 1.
In this paper, solution points are the Gauss points, which are positioned atξ1=−15/5,ξ2=0,ξ3=15/5 for K=3. Correction function for gDG scheme is used, which is expressed as gDG,L=−1k/2LK−LK−1,gDG,R=1/2LK+LK−1, where LK is Legendre polynomials. More specifically, L3ξ=1/25ξ3−3ξ and L2ξ=1/23ξ2−1.


## 2.1. Brief Review of Third-Order WCNS

Equation (1) is discretized by a finite difference scheme at every node xi,yj(2)Uti,j=Ei,j′=−Fxi,j−Gyi,j,where Ei,j′ is the approximation of the spatial derivatives. The following fourth-order flux difference operators are used(3)Fxi,j=43hF̂i+1/2,j−F̂i−1/2,j−16hFi+1,j−Fi−1,j,Gyi,j=43hĜi,j+1/2−Ĝi,j−1/2−16hGi,j+1−Gi,j−1,where h is the grid spacing. Here, F̂i−1/2,j=F̂i−1/2Ui−1/2,jL,Ui−1/2,jR,Ĝi,j−1/2=Ĝj−1/2Ui,j−1/2L,Ui,j−1/2R 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.Figure 1
Distribution of solution points and flux points of 2D WCNS ini,j-th cell. Solution points are marked in blue dot and flux points in orange square.To capture discontinuities, the conservative variables at cell edgesUi−1/2,jR 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 stencilSi=Ui−1,j,Ui,j,Ui+1,j into two smaller substencils Si1=Ui−1,j,Ui,j and Si2=Ui,j,Ui+1,j.Step 2.
Construct an interpolation polynomialpmx in each substencil Sim,m=1,2, which satisfies pmxi=Ui,j. We have(4)p1xi−1/2=12Ui−1,j+12Ui,j,p2xi−1/2=32Ui,j−12Ui+1,j.Step 3.
Calculate the nonlinear weightsωk based on linear weights dm and a smoothness indicator ISk for each substencil Sim,m=1,2. The linear weights are(5)d1=34,d2=14.The improved Z-type nonlinear weights [24, 25] are defined by(6)ωk=βk∑m=12βm,(7)βk=dkCq+τ5ε+ISkpz,(8)k=1,2,where ε=10−6 is a small number, ISk is a smoothness indicator, τ5=IS0−IS1, and Cq and pz are parameters. We use pz=2 and Cq=4 to improve the performance of the nonlinear weights.An improved smoothness indicator [25] is employed.(9)IS0=143Ui,j−4Ui−1,j+Ui−2,j2+Ui,j−2Ui−1,j+Ui−2,j2,IS1=143Ui+2,j−4Ui+1,j+3Ui,j2+Ui+2,j−2Ui+1,j+Ui,j2.Step 4.
Get weighted interpolation value atxi−1/2.(10)Ui−1/2,jR=∑m=12ωmpmxi−1/2.
The formulations ofUi−1/2,jL is symmetric with Ui−1/2,jR by i.


## 2.2. Review of Third-Order CPR

For the third-order CPR framework, there are32 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(11)∂Ui,j∂t+∂F̃i,jξ,η∂ξ+∂G̃i,jξ,η∂η=0,ξ,η∈−1,1×−1,1.Figure 2
Distribution of solution points and flux points of 2D CPR ini,j-th cell. Solution points are marked in green dot and flux points on element boundary in purple square. The flux points inside the CPR cell coincide with the solution points.The reconstructed fluxF̃i,jξ,η and G̃i,jξ,η are expressed as(12)F̃i,jξ,η=Fi,jξ,η+σi,j1ξ,η,G̃i,jξ,η=Gi,jξ,η+σi,j2ξ,η,where σi,j1ξ,η and σi,j2ξ,η are flux correction polynomials(13)σi,j1ξl,ηm=F̂−1,ηm−Fi,j−1,ηmgLξl+F̂1,ηm−Fi,j1,ηmgRξl,σi,j2ξl,ηm=Ĝξl,−1−Gi,jξl,−1gLηm+Ĝξl,1−Gi,jξl,1gRηm.Here, correction functionsgLξ, gRξ, gLη, and gRη are all cubic polynomials, and F̂ and Ĝ are Riemann fluxes corresponding to F and G, respectively.(14)F̂−1,η=F̂Ui−1,j1,η,Ui,j−1,η,F̂1,η=F̂Ui,j1,η,Ui+1,j−1,η.Moreover, the solution polynomial is(15)Ui,jξ,η=∑l=13∑m=13Ui,j,l,mLlξLmηwhere Llξ and Lmη are the 1D Lagrange polynomials in the ξ and η directions.The flux polynomial is(16)Fi,jξ,η=∑l=13∑m=13Fi,j,l,mLlξ·Lmη,Gi,jξ,η=∑l=13∑m=13Gi,j,l,mLlξ·Lmη.Remark 1.
In this paper, solution points are the Gauss points, which are positioned atξ1=−15/5,ξ2=0,ξ3=15/5 for K=3. Correction function for gDG scheme is used, which is expressed as gDG,L=−1k/2LK−LK−1,gDG,R=1/2LK+LK−1, where LK is Legendre polynomials. More specifically, L3ξ=1/25ξ3−3ξ and L2ξ=1/23ξ2−1.


## 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 [26–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:(17)NI=1rr+11/2∑k=0r1/r+1−ωk/dk/∑k=0rωk/dk1/r+121/2,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.Figure 3
Three different types of cells. Troubled cells are marked in dark blue; buffer cells are in light blue; and smooth cells are in green.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 byc and the solution points of WCNS by w.


#### 3.2.1. 1D Interface Treatment

First, WCNS in thei+1-th cell needs information of ghost points gi,k,k=1,2,3 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 ci,k,k=1,2,3. That is,(18)ugi,k=∑k=13uci,kLk.Second, we applyui+1/2,interface−CPR provided by i-th CPR cell and ui+1/2,interface+WCNS computed in i+1-th WCNS cell with the values of ghost points to calculate the numerical fluxes at interface:(19)f̂i+1/2,interface=fui+1/2,interface−CPR,ui+1/2,interface+WCNS.In addition, information exchange of switching strategy at different moments should also be considered. Wheni-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.Figure 4
1D grid, different solution points and ghost points.Step 1.
Divide the stencilS=wi−1,3,wi,1,wi,2 into two substencils S1=wi−1,3,wi,1,S2=wi,1,wi,2.Step 2.
Construct interpolation polynomialsu1ci,1 and u2ci,1 in S1 and S2, respectively, we have(20)u1ci,1=uwi−1,3l1ci,1+uwi,1l2ci,1,u2ci,1=uwi,1l1ci,1+uwi,2l2ci,1,where l1 and l2 are linear Lagrange polynomials.Step 3.
Calculate the nonlinear weights based on linear weightsd1 and d2. The linear weights can be obtained by solving the following equations:(21)d1u1ci,1+d2u2ci,1=uwi−1,3L1ci,1+uwi,1L2ci,1+uwi,2L3ci,1.
Then, nonlinear weights can be calculated by equation (6).Step 4.
Get the interpolation valueuci,1=ω1u1ci,1+ω2u2ci,1.


#### 3.2.2. 2D Interface Treatment

In this subsection, we consider the interface treatment betweeni,j-th CPR cell and i+1,j-th WCNS cell in 2D computational domain in Figure 5(a).Figure 5
2D different points near the interface of hybrid scheme. (a) 2D grid and different solution points. (b) 2Dx direction interpolation. (c) 2D y direction interpolation. (d) 2D ghost points.
(a)(b)(c)(d)Similarly, WCNS needs information of ghost pointsgi,j,i0,j0,i0=1,2,3;j0=1,2,3 from the i,j-th CPR cell. The solution values at ghost points of WCNS in the CPR cell are calculated via CPR interpolation using data at ci,j,i0,j0,i0=1,2,3;j0=1,2,3. 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(22)ugi,j,i1,j1=∑i0=13∑j0=13uci,j,i0,j0Li0ξmLj0ηm,at the ghost points gi,j,i1,j1,i1=1,2,3;j1=1,2,3.Then, we calculate the numerical fluxes at cell interface. As shown in Figure6, two types of flux points, fpkWCNS and fpkCPR, do not coincide with each other. Thus, ufpkWCNS− and ufpkCPR+ need to be interpolated before calculating numerical fluxes. ufpkWCNS− is obtained by CPR linear interpolation using equation (15) based on values in i,j-th cell, while ufpkCPR+ is obtained through WCNS nonlinear interpolation dimension by dimension. The interpolated values of ufpkWCNS− and ufpkCPR+ are denoted by ufpkWCNS−CPR,interp and ufpkCPR+WCNS,interp correspondingly.Therefore, for the CPR scheme, the numerical flux atfpkCPR is(23)f̂fpkCPR=fufpkCPR−,ufpkCPR+WCNS,interp.Meanwhile, the numerical flux atfpkWCNS for WCNS is(24)f̂fpkWCNS=fufpkWCNS−CPR,interp,ufpkWCNS+.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.Figure 6
Distribution of different flux points on cell interfaces. Flux points of CPR are marked in purple square and those of WCNS in orange square.


## 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–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:(17)NI=1rr+11/2∑k=0r1/r+1−ωk/dk/∑k=0rωk/dk1/r+121/2,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.Figure 3
Three different types of cells. Troubled cells are marked in dark blue; buffer cells are in light blue; and smooth cells are in green.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 byc and the solution points of WCNS by w.


### 3.2.1. 1D Interface Treatment

First, WCNS in thei+1-th cell needs information of ghost points gi,k,k=1,2,3 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 ci,k,k=1,2,3. That is,(18)ugi,k=∑k=13uci,kLk.Second, we applyui+1/2,interface−CPR provided by i-th CPR cell and ui+1/2,interface+WCNS computed in i+1-th WCNS cell with the values of ghost points to calculate the numerical fluxes at interface:(19)f̂i+1/2,interface=fui+1/2,interface−CPR,ui+1/2,interface+WCNS.In addition, information exchange of switching strategy at different moments should also be considered. Wheni-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.Figure 4
1D grid, different solution points and ghost points.Step 1.
Divide the stencilS=wi−1,3,wi,1,wi,2 into two substencils S1=wi−1,3,wi,1,S2=wi,1,wi,2.Step 2.
Construct interpolation polynomialsu1ci,1 and u2ci,1 in S1 and S2, respectively, we have(20)u1ci,1=uwi−1,3l1ci,1+uwi,1l2ci,1,u2ci,1=uwi,1l1ci,1+uwi,2l2ci,1,where l1 and l2 are linear Lagrange polynomials.Step 3.
Calculate the nonlinear weights based on linear weightsd1 and d2. The linear weights can be obtained by solving the following equations:(21)d1u1ci,1+d2u2ci,1=uwi−1,3L1ci,1+uwi,1L2ci,1+uwi,2L3ci,1.
Then, nonlinear weights can be calculated by equation (6).Step 4.
Get the interpolation valueuci,1=ω1u1ci,1+ω2u2ci,1.


### 3.2.2. 2D Interface Treatment

In this subsection, we consider the interface treatment betweeni,j-th CPR cell and i+1,j-th WCNS cell in 2D computational domain in Figure 5(a).Figure 5
2D different points near the interface of hybrid scheme. (a) 2D grid and different solution points. (b) 2Dx direction interpolation. (c) 2D y direction interpolation. (d) 2D ghost points.
(a)(b)(c)(d)Similarly, WCNS needs information of ghost pointsgi,j,i0,j0,i0=1,2,3;j0=1,2,3 from the i,j-th CPR cell. The solution values at ghost points of WCNS in the CPR cell are calculated via CPR interpolation using data at ci,j,i0,j0,i0=1,2,3;j0=1,2,3. 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(22)ugi,j,i1,j1=∑i0=13∑j0=13uci,j,i0,j0Li0ξmLj0ηm,at the ghost points gi,j,i1,j1,i1=1,2,3;j1=1,2,3.Then, we calculate the numerical fluxes at cell interface. As shown in Figure6, two types of flux points, fpkWCNS and fpkCPR, do not coincide with each other. Thus, ufpkWCNS− and ufpkCPR+ need to be interpolated before calculating numerical fluxes. ufpkWCNS− is obtained by CPR linear interpolation using equation (15) based on values in i,j-th cell, while ufpkCPR+ is obtained through WCNS nonlinear interpolation dimension by dimension. The interpolated values of ufpkWCNS− and ufpkCPR+ are denoted by ufpkWCNS−CPR,interp and ufpkCPR+WCNS,interp correspondingly.Therefore, for the CPR scheme, the numerical flux atfpkCPR is(23)f̂fpkCPR=fufpkCPR−,ufpkCPR+WCNS,interp.Meanwhile, the numerical flux atfpkWCNS for WCNS is(24)f̂fpkWCNS=fufpkWCNS−CPR,interp,ufpkWCNS+.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.Figure 6
Distribution of different flux points on cell interfaces. Flux points of CPR are marked in purple square and those of WCNS in orange square.


## 3.2.1. 1D Interface Treatment

First, WCNS in thei+1-th cell needs information of ghost points gi,k,k=1,2,3 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 ci,k,k=1,2,3. That is,(18)ugi,k=∑k=13uci,kLk.Second, we applyui+1/2,interface−CPR provided by i-th CPR cell and ui+1/2,interface+WCNS computed in i+1-th WCNS cell with the values of ghost points to calculate the numerical fluxes at interface:(19)f̂i+1/2,interface=fui+1/2,interface−CPR,ui+1/2,interface+WCNS.In addition, information exchange of switching strategy at different moments should also be considered. Wheni-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.Figure 4
1D grid, different solution points and ghost points.Step 1.
Divide the stencilS=wi−1,3,wi,1,wi,2 into two substencils S1=wi−1,3,wi,1,S2=wi,1,wi,2.Step 2.
Construct interpolation polynomialsu1ci,1 and u2ci,1 in S1 and S2, respectively, we have(20)u1ci,1=uwi−1,3l1ci,1+uwi,1l2ci,1,u2ci,1=uwi,1l1ci,1+uwi,2l2ci,1,where l1 and l2 are linear Lagrange polynomials.Step 3.
Calculate the nonlinear weights based on linear weightsd1 and d2. The linear weights can be obtained by solving the following equations:(21)d1u1ci,1+d2u2ci,1=uwi−1,3L1ci,1+uwi,1L2ci,1+uwi,2L3ci,1.
Then, nonlinear weights can be calculated by equation (6).Step 4.
Get the interpolation valueuci,1=ω1u1ci,1+ω2u2ci,1.


## 3.2.2. 2D Interface Treatment

In this subsection, we consider the interface treatment betweeni,j-th CPR cell and i+1,j-th WCNS cell in 2D computational domain in Figure 5(a).Figure 5
2D different points near the interface of hybrid scheme. (a) 2D grid and different solution points. (b) 2Dx direction interpolation. (c) 2D y direction interpolation. (d) 2D ghost points.
(a)(b)(c)(d)Similarly, WCNS needs information of ghost pointsgi,j,i0,j0,i0=1,2,3;j0=1,2,3 from the i,j-th CPR cell. The solution values at ghost points of WCNS in the CPR cell are calculated via CPR interpolation using data at ci,j,i0,j0,i0=1,2,3;j0=1,2,3. 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(22)ugi,j,i1,j1=∑i0=13∑j0=13uci,j,i0,j0Li0ξmLj0ηm,at the ghost points gi,j,i1,j1,i1=1,2,3;j1=1,2,3.Then, we calculate the numerical fluxes at cell interface. As shown in Figure6, two types of flux points, fpkWCNS and fpkCPR, do not coincide with each other. Thus, ufpkWCNS− and ufpkCPR+ need to be interpolated before calculating numerical fluxes. ufpkWCNS− is obtained by CPR linear interpolation using equation (15) based on values in i,j-th cell, while ufpkCPR+ is obtained through WCNS nonlinear interpolation dimension by dimension. The interpolated values of ufpkWCNS− and ufpkCPR+ are denoted by ufpkWCNS−CPR,interp and ufpkCPR+WCNS,interp correspondingly.Therefore, for the CPR scheme, the numerical flux atfpkCPR is(23)f̂fpkCPR=fufpkCPR−,ufpkCPR+WCNS,interp.Meanwhile, the numerical flux atfpkWCNS for WCNS is(24)f̂fpkWCNS=fufpkWCNS−CPR,interp,ufpkWCNS+.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.Figure 6
Distribution of different flux points on cell interfaces. Flux points of CPR are marked in purple square and those of WCNS in orange square.


## 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.TheL∞ error and L2 error are computed to measure local solution quality and global solution quality, respectively. The L2 error is estimated by(25)ErrorL2=∑i=1NUih−Uri2N,where Uih and Uri are the numerical solution and the exact solution at the i-th solution point. Here, N is the total number of solution points.


### 4.1. 1D Linear Wave Equation

Consider the 1D linear wave equation(26)ut+ux=0,with the initial condition(27)ux,0=sinπLx,x∈−L,Land 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 Table1 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).Table 1
Test of single schemes and hybrid WCNS-CPR scheme in−6,6 for 1D linear equation.
NormDOFsWCNSCPRHybrid WCNS-CPRErrorOrderErrorOrderErrorOrderL∞1203.00E−4—1.28E−4—3.00E−4—2403.76E−53.001.61E−52.993.40E−53.144804.70E−63.002.02E−62.992.02E−64.07L21202.12E−4—7.99E−5—2.07E−4—2402.66E−52.999.96E−63.001.60E−53.694803.32E−63.001.24E−63.011.24E−63.69Secondly, 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.Figure 7
Numerical results of third-order single scheme and hybrid WCNS-CPR scheme. (a) WCNS in−6,6 with DOFs=120. (b) CPR in −6,6 with DOFs=120. (c) Hybrid WCNS-CPR scheme in −6,6 with DOFs=120. (d) Hybrid WCNS-CPR scheme in −6,6 with DOFs=240. (e) Hybrid WCNS-CPR scheme in −6,6 with DOFs=480.
(a)(b)(c)(d)(e)


### 4.2. 1D Euler Equations

Consider 1D Euler equations(28)ut+fux=0,where(29)u=ρρvE,fu=ρvρv2+pvE+p,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ρv2 with γ=1.4. 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](30)ρ,u,p=0.125,0,0.1,−5≤x<0,1.0,0,1.0,0≤x≤5,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.Figure 8
Density numerical solution and shock detection result of the hybrid WCNS-CPR scheme for the 1D Sod problem.


#### 4.2.2. Shu-Osher Problem

Consider the Shu-Osher problem with initial conditions [30](31)ρ,u,p=3.857143,2.629369,10.33333,x<4,1.0+0.2sin5x,0,1.0,x≥4,and Dirichlet boundary conditions. Time step ΔT=0.0001, DOFs=600, and δNI=0.2 are used to solve the problem till T=1.8 by the hybrid WCNS-CPR scheme.The reference solution is obtained by the WCNS3 scheme withDOFs=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.Figure 9
Density numerical solution and shock detection result of the hybrid WCNS-CPR scheme for the 1D Shu-Osher problem.


#### 4.2.3. Shock-Tube Problem of Lax

Consider the Lax problem with initial condition [31](32)ρ,u,p=0.445,0.698,3.528,0.5,0,0.571,−5≤x<0,0≤x≤5,and Dirichlet boundary conditions. Time step ΔT=0.0001, DOFs=300, and δNI=0.2 are adopted to solve the problem till T=1.5 by the hybrid WCNS-CPR scheme. The reference solution is obtained by the WCNS3 scheme with DOFs=2100. We know from Figure 10 that at x≈2.3 and x≈3.8, shocks and discontinuities are all captured without obvious oscillations.Figure 10
Density numerical solution and shock detection result of the hybrid WCNS-CPR scheme for the 1D Lax problem.


### 4.3. 2D Euler Equations

Let us consider the 2D Euler equations of gas dynamics:(33)∂∂tρρuρvE+∂∂xρuρu2+pρuvuE+p+∂∂yρvρuvρv2+pvE+p=0,(34)p=γ−1·E−ρ2u2+v2. 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.


#### 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, ρ,u,v,p=1,1,1,1, with the isotropic vortex perturbations added to the mean flow. There are perturbations in u,v and T=p/ρ and no perturbation in entropy S=p/ργ:(35)Δu,Δv=ε2πe0.51−r2−y,x,ΔT=−γ−1ε28γπ2e1−r2,where r=x2+y2 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.Figure11 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.Figure 11
Distribution of different cells, errors, and NI in thex direction at T=1.0 for DOFs=120,240,480 (from left to right).From Table2, 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.Table 2
Accuracy test of single scheme and hybrid WCNS-CPR scheme for the Euler equations.
NormDOFsWCNSCPRHybrid WCNS-CPRErrorOrderErrorOrderErrorOrderL∞1204.86E−3—2.14E−3—4.86E−3—2406.87E−42.823.17E−42.766.86E−42.824808.70E−52.986.48E−52.296.48E−53.40L21202.21E−4—1.04E−4—2.34E−4—2403.18E−52.801.95E−52.424.20E−52.484804.02E−62.983.76E−62.373.76E−63.48From Table3, 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. Ratio=timeofWCNS∗/timeofCPR.Table 3
The comparison of CPU time in isentropic vortex problem.
DOFsWCNS∗CPRTime (min)RatioTime (min)120×120452.3719240×2401812.1385480×4807292.01363


#### 4.3.2. 2D Riemann Problem

We solve the Riemann problem for equation (33) on the computational domain of x,y∈0,1×0,1 with initial conditions(36)ρ,u,v,px,y,0=1.5,0,0,1.5,if0.8<x<1.0and0.8<y<1.0,0.5323,1.206,0,0.3,if0<x<0.8and0.8<y<1.0,0.138,1.206,1.206,0.029,if0<x<0.8and0<y<0.8,0.5323,0,1.206,0.3,if0.8<x<1.0and0<y<0.8,and Dirichlet boundary conditions until T=0.8. The parameter is δNI=0.1 in this test.Compared with the results of WCNS in Figure12, 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.Figure 12
Result calculated by WCNS.DOFs=960×960. Density from 0.13 to 1.78 with 25 contours.Figure 13
Results calculated by the hybrid WCNS-CPR scheme.NI=0.1. DOFs=960×960. (a) Density from 0.13 to 1.78 with 25 contours. (b) Distribution of different cells: red represents WCNS, blue represents CPR. (c) NI in the x direction. (d) NI in the y direction.
(a)(b)(c)(d)


#### 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 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 initial condition is(37)ρ,u,v,p=1.4,0.0,0.0,1.0,8.0,7.145,−4.125,116.5,ify<3x−16,ify≥3x−16.The simulation is run untilT=0.2 using time step ΔT=0.0001.As shown in Figures15(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 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.Figure 14
Density counters from 1.5 to 21.7 with 30 contours of double Mach reflection problem calculated by single WCNS.DOFs=1200×300.Figure 15
Results of the hybrid scheme.δNI=0.9. DOFs=1200×300. Density from 1.5 to 21.7 with 30 contours. (a) Density contours from 1.5 to 21.7 with 30 contours. (b) Distribution of different cells: the red area represents WCNS cells, and the blue area represents CPR cells. (c) Distribution of NI in the x direction. (d) Distribution of NI in the y direction.
(a)(b)(c)(d)Table 4
CPU time comparison on solving the double Mach reflection problem.
DOFsWCNSHybrid WCNS-CPR schemeδNI=0.1δNI=0.5δNI=0.9600×15063 min60 min58 min56 min1200×300314 min217 min207 min202 min


## 4.1. 1D Linear Wave Equation

Consider the 1D linear wave equation(26)ut+ux=0,with the initial condition(27)ux,0=sinπLx,x∈−L,Land 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 Table1 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).Table 1
Test of single schemes and hybrid WCNS-CPR scheme in−6,6 for 1D linear equation.
NormDOFsWCNSCPRHybrid WCNS-CPRErrorOrderErrorOrderErrorOrderL∞1203.00E−4—1.28E−4—3.00E−4—2403.76E−53.001.61E−52.993.40E−53.144804.70E−63.002.02E−62.992.02E−64.07L21202.12E−4—7.99E−5—2.07E−4—2402.66E−52.999.96E−63.001.60E−53.694803.32E−63.001.24E−63.011.24E−63.69Secondly, 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.Figure 7
Numerical results of third-order single scheme and hybrid WCNS-CPR scheme. (a) WCNS in−6,6 with DOFs=120. (b) CPR in −6,6 with DOFs=120. (c) Hybrid WCNS-CPR scheme in −6,6 with DOFs=120. (d) Hybrid WCNS-CPR scheme in −6,6 with DOFs=240. (e) Hybrid WCNS-CPR scheme in −6,6 with DOFs=480.
(a)(b)(c)(d)(e)


## 4.2. 1D Euler Equations

Consider 1D Euler equations(28)ut+fux=0,where(29)u=ρρvE,fu=ρvρv2+pvE+p,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ρv2 with γ=1.4. 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](30)ρ,u,p=0.125,0,0.1,−5≤x<0,1.0,0,1.0,0≤x≤5,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.Figure 8
Density numerical solution and shock detection result of the hybrid WCNS-CPR scheme for the 1D Sod problem.


### 4.2.2. Shu-Osher Problem

Consider the Shu-Osher problem with initial conditions [30](31)ρ,u,p=3.857143,2.629369,10.33333,x<4,1.0+0.2sin5x,0,1.0,x≥4,and Dirichlet boundary conditions. Time step ΔT=0.0001, DOFs=600, and δNI=0.2 are used to solve the problem till T=1.8 by the hybrid WCNS-CPR scheme.The reference solution is obtained by the WCNS3 scheme withDOFs=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.Figure 9
Density numerical solution and shock detection result of the hybrid WCNS-CPR scheme for the 1D Shu-Osher problem.


### 4.2.3. Shock-Tube Problem of Lax

Consider the Lax problem with initial condition [31](32)ρ,u,p=0.445,0.698,3.528,0.5,0,0.571,−5≤x<0,0≤x≤5,and Dirichlet boundary conditions. Time step ΔT=0.0001, DOFs=300, and δNI=0.2 are adopted to solve the problem till T=1.5 by the hybrid WCNS-CPR scheme. The reference solution is obtained by the WCNS3 scheme with DOFs=2100. We know from Figure 10 that at x≈2.3 and x≈3.8, shocks and discontinuities are all captured without obvious oscillations.Figure 10
Density numerical solution and shock detection result of the hybrid WCNS-CPR scheme for the 1D Lax problem.


## 4.2.1. Shock-Tube Problem of Sod

Consider the Sod problem with initial conditions [29](30)ρ,u,p=0.125,0,0.1,−5≤x<0,1.0,0,1.0,0≤x≤5,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.Figure 8
Density numerical solution and shock detection result of the hybrid WCNS-CPR scheme for the 1D Sod problem.


## 4.2.2. Shu-Osher Problem

Consider the Shu-Osher problem with initial conditions [30](31)ρ,u,p=3.857143,2.629369,10.33333,x<4,1.0+0.2sin5x,0,1.0,x≥4,and Dirichlet boundary conditions. Time step ΔT=0.0001, DOFs=600, and δNI=0.2 are used to solve the problem till T=1.8 by the hybrid WCNS-CPR scheme.The reference solution is obtained by the WCNS3 scheme withDOFs=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.Figure 9
Density numerical solution and shock detection result of the hybrid WCNS-CPR scheme for the 1D Shu-Osher problem.


## 4.2.3. Shock-Tube Problem of Lax

Consider the Lax problem with initial condition [31](32)ρ,u,p=0.445,0.698,3.528,0.5,0,0.571,−5≤x<0,0≤x≤5,and Dirichlet boundary conditions. Time step ΔT=0.0001, DOFs=300, and δNI=0.2 are adopted to solve the problem till T=1.5 by the hybrid WCNS-CPR scheme. The reference solution is obtained by the WCNS3 scheme with DOFs=2100. We know from Figure 10 that at x≈2.3 and x≈3.8, shocks and discontinuities are all captured without obvious oscillations.Figure 10
Density numerical solution and shock detection result of the hybrid WCNS-CPR scheme for the 1D Lax problem.


## 4.3. 2D Euler Equations

Let us consider the 2D Euler equations of gas dynamics:(33)∂∂tρρuρvE+∂∂xρuρu2+pρuvuE+p+∂∂yρvρuvρv2+pvE+p=0,(34)p=γ−1·E−ρ2u2+v2. 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.


### 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, ρ,u,v,p=1,1,1,1, with the isotropic vortex perturbations added to the mean flow. There are perturbations in u,v and T=p/ρ and no perturbation in entropy S=p/ργ:(35)Δu,Δv=ε2πe0.51−r2−y,x,ΔT=−γ−1ε28γπ2e1−r2,where r=x2+y2 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.Figure11 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.Figure 11
Distribution of different cells, errors, and NI in thex direction at T=1.0 for DOFs=120,240,480 (from left to right).From Table2, 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.Table 2
Accuracy test of single scheme and hybrid WCNS-CPR scheme for the Euler equations.
NormDOFsWCNSCPRHybrid WCNS-CPRErrorOrderErrorOrderErrorOrderL∞1204.86E−3—2.14E−3—4.86E−3—2406.87E−42.823.17E−42.766.86E−42.824808.70E−52.986.48E−52.296.48E−53.40L21202.21E−4—1.04E−4—2.34E−4—2403.18E−52.801.95E−52.424.20E−52.484804.02E−62.983.76E−62.373.76E−63.48From Table3, 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. Ratio=timeofWCNS∗/timeofCPR.Table 3
The comparison of CPU time in isentropic vortex problem.
DOFsWCNS∗CPRTime (min)RatioTime (min)120×120452.3719240×2401812.1385480×4807292.01363


### 4.3.2. 2D Riemann Problem

We solve the Riemann problem for equation (33) on the computational domain of x,y∈0,1×0,1 with initial conditions(36)ρ,u,v,px,y,0=1.5,0,0,1.5,if0.8<x<1.0and0.8<y<1.0,0.5323,1.206,0,0.3,if0<x<0.8and0.8<y<1.0,0.138,1.206,1.206,0.029,if0<x<0.8and0<y<0.8,0.5323,0,1.206,0.3,if0.8<x<1.0and0<y<0.8,and Dirichlet boundary conditions until T=0.8. The parameter is δNI=0.1 in this test.Compared with the results of WCNS in Figure12, 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.Figure 12
Result calculated by WCNS.DOFs=960×960. Density from 0.13 to 1.78 with 25 contours.Figure 13
Results calculated by the hybrid WCNS-CPR scheme.NI=0.1. DOFs=960×960. (a) Density from 0.13 to 1.78 with 25 contours. (b) Distribution of different cells: red represents WCNS, blue represents CPR. (c) NI in the x direction. (d) NI in the y direction.
(a)(b)(c)(d)


### 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 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 initial condition is(37)ρ,u,v,p=1.4,0.0,0.0,1.0,8.0,7.145,−4.125,116.5,ify<3x−16,ify≥3x−16.The simulation is run untilT=0.2 using time step ΔT=0.0001.As shown in Figures15(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 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.Figure 14
Density counters from 1.5 to 21.7 with 30 contours of double Mach reflection problem calculated by single WCNS.DOFs=1200×300.Figure 15
Results of the hybrid scheme.δNI=0.9. DOFs=1200×300. Density from 1.5 to 21.7 with 30 contours. (a) Density contours from 1.5 to 21.7 with 30 contours. (b) Distribution of different cells: the red area represents WCNS cells, and the blue area represents CPR cells. (c) Distribution of NI in the x direction. (d) Distribution of NI in the y direction.
(a)(b)(c)(d)Table 4
CPU time comparison on solving the double Mach reflection problem.
DOFsWCNSHybrid WCNS-CPR schemeδNI=0.1δNI=0.5δNI=0.9600×15063 min60 min58 min56 min1200×300314 min217 min207 min202 min


## 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, ρ,u,v,p=1,1,1,1, with the isotropic vortex perturbations added to the mean flow. There are perturbations in u,v and T=p/ρ and no perturbation in entropy S=p/ργ:(35)Δu,Δv=ε2πe0.51−r2−y,x,ΔT=−γ−1ε28γπ2e1−r2,where r=x2+y2 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.Figure11 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.Figure 11
Distribution of different cells, errors, and NI in thex direction at T=1.0 for DOFs=120,240,480 (from left to right).From Table2, 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.Table 2
Accuracy test of single scheme and hybrid WCNS-CPR scheme for the Euler equations.
NormDOFsWCNSCPRHybrid WCNS-CPRErrorOrderErrorOrderErrorOrderL∞1204.86E−3—2.14E−3—4.86E−3—2406.87E−42.823.17E−42.766.86E−42.824808.70E−52.986.48E−52.296.48E−53.40L21202.21E−4—1.04E−4—2.34E−4—2403.18E−52.801.95E−52.424.20E−52.484804.02E−62.983.76E−62.373.76E−63.48From Table3, 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. Ratio=timeofWCNS∗/timeofCPR.Table 3
The comparison of CPU time in isentropic vortex problem.
DOFsWCNS∗CPRTime (min)RatioTime (min)120×120452.3719240×2401812.1385480×4807292.01363


## 4.3.2. 2D Riemann Problem

We solve the Riemann problem for equation (33) on the computational domain of x,y∈0,1×0,1 with initial conditions(36)ρ,u,v,px,y,0=1.5,0,0,1.5,if0.8<x<1.0and0.8<y<1.0,0.5323,1.206,0,0.3,if0<x<0.8and0.8<y<1.0,0.138,1.206,1.206,0.029,if0<x<0.8and0<y<0.8,0.5323,0,1.206,0.3,if0.8<x<1.0and0<y<0.8,and Dirichlet boundary conditions until T=0.8. The parameter is δNI=0.1 in this test.Compared with the results of WCNS in Figure12, 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.Figure 12
Result calculated by WCNS.DOFs=960×960. Density from 0.13 to 1.78 with 25 contours.Figure 13
Results calculated by the hybrid WCNS-CPR scheme.NI=0.1. DOFs=960×960. (a) Density from 0.13 to 1.78 with 25 contours. (b) Distribution of different cells: red represents WCNS, blue represents CPR. (c) NI in the x direction. (d) NI in the y direction.
(a)(b)(c)(d)


## 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 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 initial condition is(37)ρ,u,v,p=1.4,0.0,0.0,1.0,8.0,7.145,−4.125,116.5,ify<3x−16,ify≥3x−16.The simulation is run untilT=0.2 using time step ΔT=0.0001.As shown in Figures15(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 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.Figure 14
Density counters from 1.5 to 21.7 with 30 contours of double Mach reflection problem calculated by single WCNS.DOFs=1200×300.Figure 15
Results of the hybrid scheme.δNI=0.9. DOFs=1200×300. Density from 1.5 to 21.7 with 30 contours. (a) Density contours from 1.5 to 21.7 with 30 contours. (b) Distribution of different cells: the red area represents WCNS cells, and the blue area represents CPR cells. (c) Distribution of NI in the x direction. (d) Distribution of NI in the y direction.
(a)(b)(c)(d)Table 4
CPU time comparison on solving the double Mach reflection problem.
DOFsWCNSHybrid WCNS-CPR schemeδNI=0.1δNI=0.5δNI=0.9600×15063 min60 min58 min56 min1200×300314 min217 min207 min202 min


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



---

*Source: 8825445-2020-10-14.xml*