Energy Stability Property of the CPR Method Based on Subcell Second-Order CNNW Limiting in Solving Conservation Laws

This paper studies the energy stability property of the correction procedure via reconstruction (CPR) method with staggered flux points based on second-order subcell limiting. The CPR method with staggered flux points uses the Gauss point as the solution point, dividing flux points based on Gauss weights, with the flux points being one more point than the solution points. For subcell limiting, a shock indicator is used to detect troubled cells where discontinuities may exist. Troubled cells are calculated by the second-order subcell compact nonuniform nonlinear weighted (CNNW2) scheme, which has the same solution points as the CPR method. The smooth cells are calculated by the CPR method. The linear energy stability of the linear CNNW2 scheme is proven theoretically. Through various numerical experiments, we demonstrate that the CNNW2 scheme and CPR method based on subcell linear CNNW2 limiting are energy-stable and that the CPR method based on subcell nonlinear CNNW2 limiting is nonlinearly stable.


Introduction
Computational fluid dynamics (CFD) is widely used in industry, and corresponding CFD commercial softwares have been developed as well. However, most of these software programs use low-order schemes, which have more dissipation errors. Compared with low-order schemes, high-order schemes offer better accuracy and lower computational costs achieving the same error level, especially for problems with complex physics and geometry [1].
In 2007, Huynh [2] introduced a high-order method for solving hyperbolic conservation laws called the flux reconstruction (FR) method. The basic idea of this method is to collocate the solution points and the flux points on the cell, construct the solution polynomial and the flux polynomial by Lagrange interpolations, then use a numerical flux at the cell interface to correct the divergence of the flux polynomial. In the linear case, by selecting the g DG correction function, the FR method can recover the discontinuous Galerkin (DG) [3][4][5][6] method. From the FR method, the spectral difference (SD) [7][8][9][10] method can be recovered by selecting the g SD correction function. In 2009, Wang and Gao [11,12] extended the FR method to two-dimensional (2D) triangle grids and hybrid grids and proposed the lifting collocation penalty (LCP) method. The FR and LCP methods are closely related, and are collectively referred to as the correction procedure via reconstruction (CPR) method.
In 2011, Vincent et al. [13] constructed an energy-stable FR method for a one-dimensional (1D) scalar convection equation. This method is controlled by a scalar parameter c; when c is within a certain range, the method is energy-stable, which is called the energy stable flux reconstruction (ESFR) method. The energy stability of the 1D linear convection equation has been proven. Subsequently, Jameson et al. [14] studied the stability of the ESFR method for a 1D nonlinear equation, and their results showed that aliasing errors are introduced in the discretization of nonlinear fluxes, potentially leading to instability. In 2015, Sheshadri et al. [15] studied the energy stability of the 2D linear convection equation on Cartesian grids, with the results showing that for uniform grids the energy is stable when the parameter c is non-negative. Spiegel et al. [16] studied the de-aliasing strategy of the FR method, which is called the over-integration method. This strategy significantly improves the stability of the FR method.
Because the CPR method is a high-order linear scheme, it tends to produce spurious numerical oscillations when solving strong discontinuities. Therefore, it is necessary to develop a shock capturing strategy suitable for the CPR method. In 2014, Sonntag and Munz [17] proposed a shock capturing algorithm based on the high-order DG method; in this method, the shock regions are divided into several subcells and calculated by a finite volume scheme. This algorithm combines the favourable characteristics of the DG method in smooth regions and the TVD finite volume method in discontinuous regions. Dumbser et al. [18,19] proposed a simple and robust posterior subcell finite volume limiter for the DG method on unstructured grids. The idea of this posterior method is to generate the candidate solution first, then go back to previous time step for correction if elements do not satisfy the posteriori detection criteria. In 2021, Kochi and Ramakrishna [20] proposed a compact subcell weighted essentially non-oscillatory (CSWENO) limiting strategy for the DG method. In 2022, Zhu et al. [21] proposed a CPR method with a subcell shock capturing strategy based on the Gauss solution points and staggered flux points. The troubled cells are detected by a shock indicator and solved by compact nonuniform nonlinear weighted (CNNW) schemes, and the smooth cells are solved by the CPR method with staggered flux points. Shi et al. [22] extended the CPR method with subcell second-order CNNW (CNNW2) limiting to unstructured quadrilateral grids. Liu et al. studied the CPR method with staggered flux points, and found that this method exhibits nonlinear stability with the g DG correction function [23].
In this paper, we investigate the energy stability of the CPR method with staggered flux points based on subcell CNNW2 limiting. The main contributions of this work are as follows: 1. The energy stability of the linear CNNW2 scheme is proved theoretically. By using L 2 energy estimation method [24], the change rate of the energy norm of the numerical solutions constructed by the CNNW2 scheme does not increase over time.
2. Various numerical experiments based on the linear advection equation and Euler equations are conducted. The results show that the CPR method with staggered flux points, the linear scheme of CNNW2, and the CPR method with staggered flux points based on subcell linear CNNW2 limiting are linearly energy stable. In addition, the CPR method with staggered flux points based on subcell nonlinear CNNW2 limiting is nonlinearly stable. This paper is organized as follows. In Section 2, the governing equations and the CPR method with staggered flux points for 1D conservation laws are briefly introduced. In Section 3, the CPR method with staggered flux points based on subcell CNNW2 limiting is provided, and a linear energy stability analysis of the linear CNNW2 method is presented. In Section 4, the numerical results for a series of test cases are presented in detail. The concluding remarks are provided in Section 5. Finally, we prove the energy stability property of the first-order CNNW2 scheme in Appendix A.

Governing Equations
Consider the 2D Euler equations [22] in the conservation form,

∂U ∂t
where U is the vector of the conservative variables, F and G are the inviscid fluxes.
where ρ is the density, u and v are the velocity components, p is the pressure, and E is the total energy. For an ideal gas, the specific heat ratio is γ = 1.4.

CPR Method with Staggered Flux Points
Consider the 1D scalar convection equation ∂u ∂t where u is the variable and f is the flux. Firstly, the equation is transformed from a physical cell E n to a computational cell with the interval I = [−1, 1]. As shown in Figure 1, the red circles represent the solution points. In this paper, we consider a linear transformation, where the linear relationship that maps I onto E n and its inverse are [2] x(ξ) = where, h is the length of the physical cell.
whereû δ = u δ ,f δ = h f δ 2 , h 2 is the Jacobian, and h is the length of the cell. For each cell, the (k − 1)-order approximate solution polynomial is constructed by a Lagrange interpolation asû whereû i δ is the state variable at the solution point ξ i and l i is the (k − 1)-order Lagrange basis function, which has the following form: In addition,f δD is constructed in a similar way: wheref i δD is the flux at the solution point ξ i .
Finally, by introducing the numerical flux and correction function, the spatial semidiscrete scheme of the original CPR [2] is obtained as where L and R represent the left and right interfaces of the cell, respectively, g L and g R are correction functions of order k, andf δI is the numerical flux. The solution points and flux points for the original CPR method are shown in Figure 2. The solution points of this method coincides with the flux points.
To simplify the description, the CPR method with staggered flux points and the original CPR method are named the CPR(Q > P) and CPR(Q = P) methods, respectively. Here, Q represents the number of flux points, and P represents the number of solution points. The difference between the CPR(Q > P) and CPR(Q = P) methods is the selection of flux points, resulting in differentf δD . The solution points and the flux points for the CPR(Q > P) method are shown in Figure 3. It should be noted that the solution points of the CPR(Q > P) method are interleaved with the flux points, and that the lengths of the subcells are determined by the Gauss quadrature weights. The difference is as follows.
The kth−order flux polynomial is constructed by the Lagrange interpolation, and the form is as follows:f wheref i δD is the flux at the flux point ξ f i and where l f i is the kth−order Lagrange basis function, which has the following form: Finally, the spatial semi-discrete scheme of CPR(Q > P) is obtained as  Note that the CPR(Q = P) method using the g DG correction function belongs to the class of ESFR methods, which was originally constructed by Vincent for 1D linear convection equations [13]. Similarly, the CPR(Q > P) method under this correction function is a version of the ESFR(Q > P) method. The nonlinear stability of this method was previously analysed theoretically in [23].

CPR(Q > P) Method with Subcell CNNW2 Limiting
A subcell CNNW2 limiting scheme was proposed by Zhu et al. [21]. Cells in the flow field are divided into troubled cells and smooth cells by shock indicators such as TVB [25] or MDHE [26]. Then, the smooth cells are solved by the CPR(Q > P) method and the troubled cells are solved by the CNNW2 scheme. The CPR (Q > P) method with subcell CNNW2 limiting needs to calculate the common numerical flux at the cell interfaces. If the left side of the interface is a troubled cell and the right side is a smooth cell, then the left value in the numerical flux is provided by the nonuniform nonlinear weighted (NNW) scheme. The right value is provided by the Lagrange interpolation of CPR.
This paper focuses on the energy stability of the CNNW2 scheme. The construction of the CNNW2 scheme for the 1D convection equation is introduced in Section 3.1; see [21] for details. In the following, we take three solution points for each cell as an example.

CNNW2 Scheme
The computational domain [a,b] is divided into N cells, and the length of each cell is h. Each cell is divided into k subcells, and the number of subcells is consistent with the number of solution points. Subcells are obtained based on Gauss quadrature weights. Considering the cell stencil of three solution points, the corresponding state variables are u i , i = 1, 2, 3, as shown in Figure 4. Here, u A and u B are the values at the subcell interfaces and d i , i = 1, 2, 3, 4, 5, 6 is the distance between the solution points and the flux points. The right value of u A and the left value of u B are obtained by the following method.
(1) The intermediate values u B are obtained by an inverse distance weighted interpolation: , , where ω i is the interpolation weight.
(2) The gradient ∂u/∂ξ is calculated by where , (3) With the gradient ∂u/∂ξ and u 2 , u A and u B are recalculated by (4) The gradient is limited in order to control numerical oscillations.
where the limiting function is defined as where m = min(u 1 , Through NNW interpolation, the left and right values of a second-order polynomial with a limiter on the subcell interface can be obtained; thus, the numerical flux at each interface can be obtained. Finally, the second-order difference operator is used to approximate the spatial derivative for the numerical flux.
where ξ f p j is the coordinate of the f p j flux point. It should be pointed out that if φ = 1, the scheme is equivalent to a linear scheme without a limiter. If φ = 0, the scheme is reduced to a first-order scheme.

Proof of the Linear Energy Stability of the CNNW2 Scheme
Consider the 1D scalar convection equation ∂u ∂t The linear energy stability of the scheme is proven by the energy estimation method [24], and periodic boundary conditions are used. The energy norm change rate of the whole physical space and the computational space has the following relationship: According to the above equation, in order to calculate the rate of change of the energy norm over time in the whole physical space, we can consider the energy norm change rate of the nth computational cell. By using the CNNW2 scheme, the spatial semi-discrete scheme of the nth computational cell is ∂u n,1 ∂t The numerical flux adopts the upwind flux [2], and the numerical flux at each interface can be obtained: Using the Gaussian integral formula, we obtain the energy norm change rate of the nth computational cell: ∂u n,1 ∂t w 1 + u n,2 ∂u n,2 ∂t w 2 + u n,3 ∂u n,3 ∂t u n,2 u n,3 where u n,i , i = 1, 2, 3 denotes the state variables at the ith solution point on the nth cell. In addition, w 1 , w 2 , w 3 are the Gauss quadrature weights corresponding to the solution points. The energy norm change rate of the whole computational space is u n,1 u n,3 . where Since According to the average value inequality   (25), the energy norm change rate of the whole physical space is obtained as Since h is positive, the above result shows that the CNNW2 scheme is linear energy stable with φ = 1. In Appendix A, we prove the energy stability property of the first-order CNNW2 scheme. It is worth mentioning that, the linear energy stability of the sucell linear CNNW2 scheme for fifth-order CPR can be proved in a similar deductive procedure.

Numerical Experiments
To simplify the description, the linear CNNW2 scheme is named CNNW2−L, the nonlinear CNNW2 scheme is named CNNW2−N. Since the CPR based on subcell limiting is in fact a hybrid CPR and CNNW scheme, we take the same abbreviation HCCS as in [21]. Then, the CPR(Q > P) method with subcell linear CNNW2 limiting is named HCCS−L2, and the CPR(Q > P) method with subcell nonlinear CNNW2 limiting is named HCCS−N2.
The linear energy stability of the CPR(Q > P), CNNW2−L and HCCS−L2 schemes are studied in the context of 1D [13] and 2D linear convection problems [27]. An isentropic vortex case is used to analyse the errors of the CPR(Q > P), HCCS−L2 and HCCS−N2 schemes. The nonlinear stability of the CPR(Q = P), CPR(Q > P) and HCCS−N2 schemes are studied in the context of 2D subsonic flow over a cylinder [28]. The nonlinear stability of the HCCS−N2 scheme is studied in the context of the 2D Kelvin-Helmholtz (KH) instability problem [29] and transonic flows around the NACA0012 airfoil [30]. The g DG correction function and the explicit third-order TVD Runge-Kutta method [31] are used in this section. The L 2 energy is expressed as where Ω is the computational domain of physical space.
If the calculation does not blow up, then the numerical method is stable. Otherwise, if the density or pressure becomes negative, resulting in calculation blowing up, the numerical method is unstable.

1D Linear Convection Equation
The 1D linear convection Equation (24) is used. The computational domain is defined in [−1, 1]. The polynomial order is 3, and the number of grid cells is 100. Using upwind flux and periodic boundary conditions, the compution time t is 20 and time steps are determined by CFL = 0.1. This case uses the TVB indicator with adjustable parameter M = 1. The initial condition is set as follows: u(x, 0) = e −20x 2 .
(35) Figure 5a shows the numerical solutions of the CPR(Q > P), CNNW2−L and HCCS−L2 schemes. Figure 5b shows the L 2 energy of these three schemes over time. The L 2 energy of the CPR(Q > P) method remains essentially constant over time. This is because the CPR(Q > P) method has high spatial accuracy and low numerical dissipation. The results of the other two schemes decrease gradually over a bounded range, which is because the CNNW2 scheme has more numerical dissipation than the CPR(Q > P) method. Thus, the three schemes are energy-stable. Figure 6 shows the evolution of troubled cells in the HCCS−L2 scheme, where red represents troubled cells.

Isentropic Vortex Test
To analyse the errors of the CPR(Q > P), HCCS−L2 and HCCS−N2 schemes, we solve the isentropic vortex problem. In this case, an isentropic vortex disturbance is added to an uniform flow. The uniform flow is set to (ρ ∞ , u ∞ , v ∞ , p ∞ ) = (1.0, 1.0, 0.0, 1.0), and T ∞ = p ∞ /ρ ∞ . The initial conditions of the vortex are set as follows: where r = (x − x c ) 2 + (y − y c ) 2 , the vortex centre (x c , y c ) = (0.0, 0.0), and the vortex strength = 5.0. The initial conditions of the flow field are as follows: The computational domain is [−10, 10] × [−10, 10]. The LLF flux and periodic boundary conditions are used. The compution time t is 0.1 and time steps are calculated using CFL = 0.1. Table 1 shows the errors of the CPR(Q > P), HCCS−L2 and HCCS−N2 schemes. The errors of the CPR(Q > P) method are the smallest, while the errors of the HCCS−N2 scheme are the largest.

Subsonic Flow Over a Cylinder
The cylinder radius is 0.5, and the far field is ten times the cylinder diameter. The simulation is run at a freestream Mach number of 0.2 with 5th-order accuracy. 178 cells are distributed in the circumferential direction, while 54 cells are distributed in the radial direction. Using the LLF flux, the compution is conducted until t = 30 with time steps determined by CFL = 0.3. The inviscid wall boundary condition is imposed on the wall surface, and the far-field boundary condition is imposed on the far field. This case uses the MDHE indicator with an adjustable parameter a = 0.0005. When solving this problem with straight-sided quadrilateral grids, a pair of wake vortices appears at the rear end of the cylinder. Figure 10 shows the density contours calculated by the CPR(Q = P) and CPR(Q > P) schemes. It can be seen that the CPR(Q = P) scheme blows up at t = 14.25 due to the effect of aliasing errors. Figure 10a shows the flow field before the blow up. The CPR(Q > P) scheme simulates stably to the end. Figure 11 shows the density contour and the distribution of troubled cells calculated by the HCCS−N2 scheme with troubled cells represented by red color. The CPR(Q > P) scheme with subcell nonlinear CNNW2 limiting is shown to be still nonlinearly stable.

2D Kelvin-Helmholtz Instability Problem
The initial conditions are set as follows: ρ(x, y) = 1 2 + 3 4 B, p(x, y) = 1, where B = tanh(15y + 7.5) − tanh(15y − 7.5). The computational domain is defined in The polynomial order is 7, and the number of grid cells is 64 × 64. Using the LLF flux and periodic boundary conditions, the compution time t is 10 with time steps calculated by CFL = 0.3. The Reynolds number goes to infinity, and the Mach number is approximately 0.6. This case uses the MDHE indicator with an adjustable parameter a = 0.5. Figure 12a shows the density contour calculated by the HCCS−N2 scheme. The high resolution of this scheme enables it to capture small-scale features. Figure 12b shows the distribution of troubled cells (in red).

2D Transonic Flow around the NACA0012 Airfoil
The initial condition of transonic flow around the NACA0012 airfoil is a freestream flow condition with Mach number Ma ∞ = 0.8 and angle of attack α = 1.25 • . The inviscid wall boundary condition is imposed on the wall surface, and the far-field boundary condition is imposed on the far field. The numbers of grid cells distributed in the circumferential and radial directions are 120 and 80, respectively, which are generated by solving a partial differential equation, as shown in Figure 13. The polynomial order is 2, and the compution time t is 50 with time steps determined by CFL = 0.5. This case uses the MDHE indicator with an adjustable parameter a = 0.0008. Figure 14 shows the density contour, pressure contour and Mach number contour calculated by the HCCS−N2 scheme. This scheme can better capture the relatively strong shock waves on the upper surface of the airfoil and relatively weak shock waves on the lower surface. Figure 15 shows the pressure coefficient distribution on the airfoil surface calculated by the HCCS−N2 scheme. While there are slight differences in the shock wave regions, the result is in good agreement with the reference solution in the other regions.

Conclusions
This paper addresses the energy stability of the CPR method with staggered flux points (CPR(Q > P)) based on second-order subcell limiting. The linear CNNW2 scheme with φ = 1 is proven to be energy stable. Through numerical tests of 1D and 2D linear convection equations, the results show that linear CNNW2 and CPR(Q > P) method with subcell linear CNNW2 limiting are energy-stable. Through numerical tests of 2D subsonic flow over a cylinder show that the CPR(Q > P) method has better nonlinear stability than the CPR(Q = P) method. The results of 2D Kelvin-Helmholtz (KH) instability problem and transonic flows around the NACA0012 airfoil indicate that the CPR(Q > P) with subcell CNNW2 limiting has good properties in both shock capturing and nonlinear stability.