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

Abstract

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.

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

2. Governing Equations and Discretization Methods

2.1. Governing Equations

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

(1)${\displaystyle \frac{\partial \mathit{U}}{\partial t}}+{\displaystyle \frac{\partial \mathit{F}}{\partial x}}+{\displaystyle \frac{\partial \mathit{G}}{\partial y}}=0,$

where $\mathit{U}$ is the vector of the conservative variables, $\mathit{F}$ and $\mathit{G}$ are the inviscid fluxes.

(2)$\mathit{U}=\left[\begin{array}{c}\rho \\ \rho u\\ \rho v\\ E\end{array}\right],\mathit{F}=\left[\begin{array}{c}\rho u\\ \rho u^2+p\\ \rho uv\\ u(E+p)\end{array}\right],\mathit{G}=\left[\begin{array}{c}\rho v\\ \rho uv\\ \rho v^2+p\\ v(E+p)\end{array}\right],$

(3)$E=\frac{p}{\gamma -1}+\frac{1}{2}\rho (u^2+v^2),$

where $\rho$ 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 $\gamma$ = 1.4.

2.2. CPR Method with Staggered Flux Points

Consider the 1D scalar convection equation

(4)${\displaystyle \frac{\partial u}{\partial t}}+{\displaystyle \frac{\partial f}{\partial x}}=0,$

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]

(5)$x\left(\xi \right)=\frac{x_n+x_{n+1}}{2}+\xi h/2\mathrm{and}\xi \left(x\right)=2\left(x-\frac{x_n+x_{n+1}}{2}\right)/h$

where, h is the length of the physical cell.

Upon transformation to computational cells. The CPR method for Equation (4) becomes

(6)${\displaystyle \frac{\partial {\widehat{u}}^{\delta }}{\partial t}}+{\displaystyle \frac{\partial {\widehat{f}}^{\delta }}{\partial \xi }}=0,$

where ${\widehat{u}}^{\delta }=u^{\delta }$, ${\widehat{f}}^{\delta }=\frac{h{f}^{\delta }}{2}$, $\frac{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

(7)${\widehat{u}}^{\delta }=\sum _{i=1}^k{\widehat{u_i}}^{\delta }{l}_i,$

where ${\widehat{u_i}}^{\delta }$ is the state variable at the solution point ${\xi }_i$ and $l_i$ is the $(k-1)$-order Lagrange basis function, which has the following form:

(8)$l_i=\prod _{j=1,j\ne i}^k\left(\frac{\xi -{\xi }_j}{{\xi }_i-{\xi }_j}\right).$

In addition, ${\widehat{f}}^{\delta D}$ is constructed in a similar way:

(9)${\widehat{f}}^{\delta D}=\sum _{i=1}^k{\widehat{f_i}}^{\delta D}{l}_i,$

where ${\widehat{f_i}}^{\delta D}$ is the flux at the solution point ${\xi }_i$.

Finally, by introducing the numerical flux and correction function, the spatial semi-discrete scheme of the original CPR [2] is obtained as

(10)${\displaystyle \frac{\partial {\widehat{u_i}}^{\delta }}{\partial t}}=-\left(\sum _{i=1}^k{\widehat{f_i}}^{\delta D}\frac{d{l}_i}{d\xi }\left({\xi }_i\right)+({\widehat{f_L}}^{\delta I}-{\widehat{f_L}}^{\delta D})\frac{d{g}_L}{d\xi }\left({\xi }_i\right)+({\widehat{f_R}}^{\delta I}-{\widehat{f_R}}^{\delta D})\frac{d{g}_R}{d\xi }\left({\xi }_i\right)\right),$

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, and ${\widehat{f}}^{\delta 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 different ${\widehat{f}}^{\delta 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:

(11)${\widehat{f}}^{\delta D}=\sum _{i=1}^{k+1}{\widehat{f_i}}^{\delta D}{l}_i^f,$

where ${\widehat{f_i}}^{\delta D}$ is the flux at the flux point ${\xi }_i^f$ and where $l_i^f$ is the kth−order Lagrange basis function, which has the following form:

(12)$l_i^f=\prod _{j=1,j\ne i}^{k+1}\left(\frac{\xi -{\xi }_j^f}{{\xi }_i^f-{\xi }_j^f}\right).$

Finally, the spatial semi-discrete scheme of CPR $(Q>P)$ is obtained as

(13)${\displaystyle \frac{\partial {\widehat{u_i}}^{\delta }}{\partial t}}=-\left(\sum _{i=1}^{k+1}{\widehat{f_i}}^{\delta D}\frac{d{l}_i^f}{d\xi }\left({\xi }_i\right)+({\widehat{f_L}}^{\delta I}-{\widehat{f_L}}^{\delta D})\frac{d{g}_L}{d\xi }\left({\xi }_i\right)+({\widehat{f_R}}^{\delta I}-{\widehat{f_R}}^{\delta D})\frac{d{g}_R}{d\xi }\left({\xi }_i\right)\right).$

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

The $g_{DG}$ correction function [13] is provided by

(14)$g_{DG,L}=\frac{{(-1)}^k}{2}(L_k-L_{k+1}),$

(15)$g_{DG,R}=\frac{1}{2}(L_k+L_{k+1}),$

where $L_k$ is the Legendre polynomial of order k. When k is 4, the expressions of the Legendre polynomials are $L_4=\left(35x^4-30x^2+3\right)/8$ and $L_5=\left(63x^5-70x^3+15x\right)/8$.

3. CPR $(\mathit{Q}>\mathit{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.

3.1. 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_A^{\left(1\right)}$, $u_B^{\left(1\right)}$ are obtained by an inverse distance weighted interpolation:

(16)${\omega }_1=\frac{\left(1/d_2\right)}{\left(1/d_2\right)+\left(1/d_3\right)},{\omega }_2=\frac{\left(1/d_3\right)}{\left(1/d_2\right)+\left(1/d_3\right)},u_A^{\left(1\right)}={\omega }_1{u}_1+{\omega }_2{u}_2,$

(17)${\omega }_3=\frac{\left(1/d_4\right)}{\left(1/d_4\right)+\left(1/d_5\right)},{\omega }_4=\frac{\left(1/d_5\right)}{\left(1/d_4\right)+\left(1/d_5\right)},u_B^{\left(1\right)}={\omega }_3{u}_2+{\omega }_4{u}_3,$

where ${\omega }_i$ is the interpolation weight.

(2) The gradient $\partial u/\partial \xi$ is calculated by

(18)$\frac{\partial u}{\partial \xi }={\omega }_5{\left(\frac{\partial u}{\partial \xi }\right)}^{\left(1\right)}+{\omega }_6{\left(\frac{\partial u}{\partial \xi }\right)}^{\left(2\right)}={\omega }_5\frac{u_2-u_A^{\left(1\right)}}{d_3}+{\omega }_6\frac{u_B^{\left(1\right)}-u_2}{d_4},$

where

(19)${\omega }_5=\frac{\left(1/d_3\right)}{\left(1/d_3\right)+\left(1/d_4\right)},{\omega }_6=\frac{\left(1/d_4\right)}{\left(1/d_3\right)+\left(1/d_4\right)}.$

(3) With the gradient $\partial u/\partial \xi$ and $u_2$, $u_A$ and $u_B$ are recalculated by

(20)$u_A^{\left(2\right)}=u_2-\frac{\partial u}{\partial \xi }{d}_3,u_B^{\left(2\right)}=u_2+\frac{\partial u}{\partial \xi }{d}_4,$

(4) The gradient is limited in order to control numerical oscillations.

(21)$u_A^R=u_2-\varphi \frac{\partial u}{\partial \xi }{d}_3,u_B^L=u_2+\varphi \frac{\partial u}{\partial \xi }{d}_4,\varphi =min\left(lim\left(u_A^{\left(2\right)}\right),lim\left(u_B^{\left(2\right)}\right)\right),$

where the limiting function is defined as

(22)$lim\left(u\right)=\left\{\begin{array}{l}min\left(1,\frac{M-u_2}{u-u_2}\right),\mathrm{if}u>u_2,\\ min\left(1,\frac{m-u_2}{u-u_2}\right),\mathrm{if}u<u_2,\\ 1,\mathrm{if}u=u_2,\end{array}\right.,$

where $m=min\left(u_1,u_2,u_3\right),M=max\left(u_1,u_2,u_3\right)$.

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.

(23)$\frac{\partial f}{\partial \xi }=\frac{{\widehat{f}}_{n,f{p}_{j+1}}-{\widehat{f}}_{n,f{p}_j}}{{\xi }_{f{p}_{j+1}}-{\xi }_{f{p}_j}},$

where ${\xi }_{f{p}_j}$ is the coordinate of the $f{p}_j$ flux point. It should be pointed out that if $\varphi =1$, the scheme is equivalent to a linear scheme without a limiter. If $\varphi =0$, the scheme is reduced to a first-order scheme.

3.2. Proof of the Linear Energy Stability of the CNNW2 Scheme

Consider the 1D scalar convection equation

(24)$\frac{\partial u}{\partial t}+\frac{\partial u}{\partial x}=0.$

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:

(25)$\frac{\partial }{\partial t}{\int }_a^b{u}^{2}dx=\frac{\partial }{\partial t}\sum _{n=1}^N{\int }_{x_n}^{x_{n+1}}{u}^{2}dx=\frac{h}{2}\frac{\partial }{\partial t}\sum _{n=1}^N{\int }_{-1}^1{u}^{2}d\xi .$

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

(26)$\begin{array}{c}\frac{\partial u_{n,1}}{\partial t}=\frac{{\widehat{f}}_{n,fp1}-{\widehat{f}}_{n,fp2}}{d_1+d_2},\\ \frac{\partial u_{n,2}}{\partial t}=\frac{{\widehat{f}}_{n,fp2}-{\widehat{f}}_{n,fp3}}{d_3+d_4},\\ \frac{\partial u_{n,3}}{\partial t}=\frac{{\widehat{f}}_{n,fp3}-{\widehat{f}}_{n,fp4}}{d_5+d_6}.\end{array}$

The numerical flux adopts the upwind flux [2], and the numerical flux at each interface can be obtained:

(27)$\begin{array}{cc}\hfill & {\widehat{f}}_{n,fp1}=u_{n-1,3}+\frac{d_6^2}{d_5+d_6}\frac{u_{n-1,3}-u_{n-1,2}}{d_4+d_5}+\frac{d_5{d}_6}{d_5+d_6}\frac{u_{n,1}-u_{n-1,3}}{d_1+d_6},\hfill \\ \hfill & {\widehat{f}}_{n,fp2}=u_{n,1}+\frac{d_2^2}{d_1+d_2}\frac{u_{n,1}-u_{n-1,3}}{d_1+d_6}+\frac{d_1{d}_2}{d_1+d_2}\frac{u_{n,2}-u_{n,1}}{d_2+d_3},\hfill \\ \hfill & {\widehat{f}}_{n,fp3}=u_{n,2}+\frac{d_4^2}{d_3+d_4}\frac{u_{n,2}-u_{n,1}}{d_2+d_3}+\frac{d_3{d}_4}{d_3+d_4}\frac{u_{n,3}-u_{n,2}}{d_4+d_5},\hfill \\ \hfill & {\widehat{f}}_{n,fp4}=u_{n,3}+\frac{d_6^2}{d_5+d_6}\frac{u_{n,3}-u_{n,2}}{d_4+d_5}+\frac{d_5{d}_6}{d_5+d_6}\frac{u_{n+1,1}-u_{n,3}}{d_1+d_6}.\hfill \end{array}$

Using the Gaussian integral formula, we obtain the energy norm change rate of the nth computational cell:

(28)$\begin{array}{cc}\hfill & \frac{\frac{1}{2}\partial {\int }_{-1}^1{u}^{2}d\xi }{\partial t}\hfill \\ \hfill & =u_{n,1}\frac{\partial u_{n,1}}{\partial t}{w}_1+u_{n,2}\frac{\partial u_{n,2}}{\partial t}{w}_2+u_{n,3}\frac{\partial u_{n,3}}{\partial t}{w}_3\hfill \\ \hfill & =\left(\frac{d_5{d}_6}{d_5+d_6}\frac{1}{d_1+d_6}-1-\frac{d_2^2}{d_1+d_2}\frac{1}{d_1+d_6}+\frac{d_1{d}_2}{d_1+d_2}\frac{1}{d_2+d_3}\right)u_{n,1}^2\hfill \\ \hfill & +\left(1+\frac{d_6^2}{d_5+d_6}\frac{1}{d_4+d_5}-\frac{d_5{d}_6}{d_5+d_6}\frac{1}{d_1+d_6}+\frac{d_2^2}{d_1+d_2}\frac{1}{d_1+d_6}\right)u_{n-1,3}{u}_{n,1}\hfill \\ \hfill & +\left(-\frac{d_6^2}{d_5+d_6}\frac{1}{d_4+d_5}\right)u_{n-1,2}{u}_{n,1}+\left(-\frac{d_1{d}_2}{d_1+d_2}\frac{1}{d_2+d_3}\right)u_{n,2}{u}_{n,1}\hfill \\ \hfill & +\left(-1-\frac{d_4^2}{d_3+d_4}\frac{1}{d_2+d_3}+\frac{d_3{d}_4}{d_3+d_4}\frac{1}{d_4+d_5}+\frac{d_1{d}_2}{d_1+d_2}\frac{1}{d_2+d_3}\right)u_{n,2}^2\hfill \\ \hfill & +\left(\frac{d_4^2}{d_3+d_4}\frac{1}{d_2+d_3}+1+\frac{d_2^2}{d_1+d_2}\frac{1}{d_1+d_6}-\frac{d_1{d}_2}{d_1+d_2}\frac{1}{d_2+d_3}\right)u_{n,1}{u}_{n,2}\hfill \\ \hfill & +\left(-\frac{d_3{d}_4}{d_3+d_4}\frac{1}{d_4+d_5}\right)u_{n,3}{u}_{n,2}+\left(-\frac{d_2^2}{d_1+d_2}\frac{1}{d_1+d_6}\right)u_{n-1,3}{u}_{n,2}\hfill \\ \hfill & +\left(\frac{d_3{d}_4}{d_3+d_4}\frac{1}{d_4+d_5}-1-\frac{d_6^2}{d_5+d_6}\frac{1}{d_4+d_5}+\frac{d_5{d}_6}{d_5+d_6}\frac{1}{d_1+d_6}\right)u_{n,3}^2\hfill \\ \hfill & +\left(1+\frac{d_4^2}{d_3+d_4}\frac{1}{d_2+d_3}-\frac{d_3{d}_4}{d_3+d_4}\frac{1}{d_4+d_5}+\frac{d_6^2}{d_5+d_6}\frac{1}{d_4+d_5}\right)u_{n,2}{u}_{n,3}\hfill \\ \hfill & +\left(-\frac{d_4^2}{d_3+d_4}\frac{1}{d_2+d_3}\right)u_{n,1}{u}_{n,3}+\left(-\frac{d_5{d}_6}{d_5+d_6}\frac{1}{d_1+d_6}\right)u_{n+1,1}{u}_{n,3},\hfill \end{array}$

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

(29)$\begin{array}{cc}\hfill & \frac{\frac{1}{2}{\sum }_{n=1}^N\partial {\int }_{-1}^1{u}^{2}d\xi }{\partial t}\hfill \\ \hfill & =\left(\frac{d_5{d}_6}{d_5+d_6}\frac{1}{d_1+d_6}-1-\frac{d_2^2}{d_1+d_2}\frac{1}{d_1+d_6}+\frac{d_1{d}_2}{d_1+d_2}\frac{1}{d_2+d_3}\right)\sum _{n=1}^N{u}_{n,1}^2\hfill \\ \hfill & +\left(-1-\frac{d_4^2}{d_3+d_4}\frac{1}{d_2+d_3}+\frac{d_3{d}_4}{d_3+d_4}\frac{1}{d_4+d_5}+\frac{d_1{d}_2}{d_1+d_2}\frac{1}{d_2+d_3}\right)\sum _{n=1}^N{u}_{n,2}^2\hfill \\ \hfill & +\left(\frac{d_3{d}_4}{d_3+d_4}\frac{1}{d_4+d_5}-1-\frac{d_6^2}{d_5+d_6}\frac{1}{d_4+d_5}+\frac{d_5{d}_6}{d_5+d_6}\frac{1}{d_1+d_6}\right)\sum _{n=1}^N{u}_{n,3}^2\hfill \\ \hfill & +\left(1+\frac{d_6^2}{d_5+d_6}\frac{1}{d_4+d_5}-2\frac{d_5{d}_6}{d_5+d_6}\frac{1}{d_1+d_6}+\frac{d_2^2}{d_1+d_2}\frac{1}{d_1+d_6}\right)\sum _{n=1}^N{u}_{n-1,3}{u}_{n,1}\hfill \\ \hfill & +\left(\frac{d_4^2}{d_3+d_4}\frac{1}{d_2+d_3}+1+\frac{d_2^2}{d_1+d_2}\frac{1}{d_1+d_6}-2\frac{d_1{d}_2}{d_1+d_2}\frac{1}{d_2+d_3}\right)\sum _{n=1}^N{u}_{n,1}{u}_{n,2}\hfill \\ \hfill & +\left(1+\frac{d_4^2}{d_3+d_4}\frac{1}{d_2+d_3}-2\frac{d_3{d}_4}{d_3+d_4}\frac{1}{d_4+d_5}+\frac{d_6^2}{d_5+d_6}\frac{1}{d_4+d_5}\right)\sum _{n=1}^N{u}_{n,2}{u}_{n,3}\hfill \\ \hfill & +\left(-\frac{d_2^2}{d_1+d_2}\frac{1}{d_1+d_6}\right)\sum _{n=1}^N{u}_{n-1,3}{u}_{n,2}+\left(-\frac{d_6^2}{d_5+d_6}\frac{1}{d_4+d_5}\right)\sum _{n=1}^N{u}_{n-1,2}{u}_{n,1}\hfill \\ \hfill & +\left(-\frac{d_4^2}{d_3+d_4}\frac{1}{d_2+d_3}\right)\sum _{n=1}^N{u}_{n,1}{u}_{n,3}.\hfill \end{array}$

Taking ${\sum }_{n=1}^N{\int }_{-1}^1{u}^{2}d\xi =E$, the above equation can be converted into

(30)$\begin{array}{cc}\hfill \frac{1}{2}\frac{\partial E}{\partial t}& =\left(a_1+b_1\right)\sum _{n=1}^N\left(u_{n,1}^2+u_{n,2}^2\right)+\left(-2a_1-2b_1\right)\sum _{n=1}^N\left(u_{n,1}{u}_{n,2}\right)\hfill \\ \hfill & +\left(c_1\right)\sum _{n=1}^N\left(u_{n,1}^2+u_{n,2}^2\right)+\left(-2c_1\right)\sum _{n=1}^N\left(u_{n,1}{u}_{n-1,2}\right)\hfill \\ \hfill & +\left(a_2+b_2\right)\sum _{n=1}^N\left(u_{n,1}^2+u_{n,3}^2\right)+\left(-2a_2-2b_2\right)\sum _{n=1}^N\left(u_{n,1}{u}_{n-1,3}\right)\hfill \\ \hfill & +\left(c_2\right)\sum _{n=1}^N\left(u_{n,1}^2+u_{n,3}^2\right)+\left(-2c_2\right)\sum _{n=1}^N\left(u_{n,1}{u}_{n,3}\right)\hfill \\ \hfill & +\left(a_3+b_3\right)\sum _{n=1}^N\left(u_{n,2}^2+u_{n,3}^2\right)+\left(-2a_3-2b_3\right)\sum _{n=1}^N\left(u_{n,2}{u}_{n,3}\right)\hfill \\ \hfill & +\left(c_3\right)\sum _{n=1}^N\left(u_{n,2}^2+u_{n,3}^2\right)+\left(-2c_3\right)\sum _{n=1}^N\left(u_{n,2}{u}_{n-1,3}\right),\hfill \end{array}$

where

(31)$\begin{array}{cc}\hfill & a_1=-\frac{1}{2}-\frac{1}{2}\frac{d_4^2}{d_3+d_4}\frac{1}{d_2+d_3}-\frac{1}{2}\frac{d_2^2}{d_1+d_2}\frac{1}{d_1+d_6},\hfill \\ \hfill & b_1=\frac{d_1{d}_2}{d_1+d_2}\frac{1}{d_2+d_3},\hfill \\ \hfill & c_1=\frac{1}{2}\frac{d_6^2}{d_5+d_6}\frac{1}{d_4+d_5},\hfill \\ \hfill & a_2=-\frac{1}{2}-\frac{1}{2}\frac{d_6^2}{d_5+d_6}\frac{1}{d_4+d_5}-\frac{1}{2}\frac{d_2^2}{d_1+d_2}\frac{1}{d_1+d_6},\hfill \\ \hfill & b_2=\frac{d_5{d}_6}{d_5+d_6}\frac{1}{d_1+d_6},\hfill \\ \hfill & c_2=\frac{1}{2}\frac{d_4^2}{d_3+d_4}\frac{1}{d_2+d_3},\hfill \\ \hfill & a_3=-\frac{1}{2}-\frac{1}{2}\frac{d_4^2}{d_3+d_4}\frac{1}{d_2+d_3}-\frac{1}{2}\frac{d_6^2}{d_5+d_6}\frac{1}{d_4+d_5},\hfill \\ \hfill & b_3=\frac{d_3{d}_4}{d_3+d_4}\frac{1}{d_4+d_5},\hfill \\ \hfill & c_3=\frac{1}{2}\frac{d_2^2}{d_1+d_2}\frac{1}{d_1+d_6},\hfill \end{array}$

Since $d_1=d_6\approx 0.2254$, $d_2=d_5\approx 0.3302$ and $d_3=d_4\approx 0.4444$, $\left|a_i\right|\ge \left|b_i+c_i\right|$ and $a_i\le 0,b_i\ge 0,c_i\ge 0,i=1,2,3$ are satisfied. According to the average value inequality $a^2+b^2\ge 2ab$,

(32)$\begin{array}{cc}\hfill \frac{1}{2}\frac{\partial E}{\partial t}& =\left(a_1+b_1\right)\sum _{n=1}^N\left(u_{n,1}^2+u_{n,2}^2\right)+\left(-2a_1-2b_1\right)\sum _{n=1}^N\left(u_{n,1}{u}_{n,2}\right)\hfill \\ \hfill & +\left(c_1\right)\sum _{n=1}^N\left(u_{n,1}^2+u_{n,2}^2\right)+\left(-2c_1\right)\sum _{n=1}^N\left(u_{n,1}{u}_{n-1,2}\right)\hfill \\ \hfill & +\left(a_2+b_2\right)\sum _{n=1}^N\left(u_{n,1}^2+u_{n,3}^2\right)+\left(-2a_2-2b_2\right)\sum _{n=1}^N\left(u_{n,1}{u}_{n-1,3}\right)\hfill \\ \hfill & +\left(c_2\right)\sum _{n=1}^N\left(u_{n,1}^2+u_{n,3}^2\right)+\left(-2c_2\right)\sum _{n=1}^N\left(u_{n,1}{u}_{n,3}\right)\hfill \\ \hfill & +\left(a_3+b_3\right)\sum _{n=1}^N\left(u_{n,2}^2+u_{n,3}^2\right)+\left(-2a_3-2b_3\right)\sum _{n=1}^N\left(u_{n,2}{u}_{n,3}\right)\hfill \\ \hfill & +\left(c_3\right)\sum _{n=1}^N\left(u_{n,2}^2+u_{n,3}^2\right)+\left(-2c_3\right)\sum _{n=1}^N\left(u_{n,2}{u}_{n-1,3}\right)\hfill \\ \hfill & \le 2\left(a_1+b_1\right)\sum _{n=1}^N\left(u_{n,1}{u}_{n,2}\right)+\left(2c_1\right)\sum _{n=1}^N\left(u_{n,1}{u}_{n-1,2}\right)\hfill \\ \hfill & +\left(-2a_1-2b_1\right)\sum _{n=1}^N\left(u_{n,1}{u}_{n,2}\right)+\left(-2c_1\right)\sum _{n=1}^N\left(u_{n,1}{u}_{n-1,2}\right)\hfill \\ \hfill & +2\left(a_2+b_2\right)\sum _{n=1}^N\left(u_{n,1}{u}_{n-1,3}\right)+\left(2c_2\right)\sum _{n=1}^N\left(u_{n,1}{u}_{n,3}\right)\hfill \\ \hfill & +\left(-2a_2-2b_2\right)\sum _{n=1}^N\left(u_{n,1}{u}_{n-1,3}\right)+\left(-2c_2\right)\sum _{n=1}^N\left(u_{n,1}{u}_{n,3}\right)\hfill \\ \hfill & +2\left(a_3+b_3\right)\sum _{n=1}^N\left(u_{n,2}{u}_{n,3}\right)+\left(2c_3\right)\sum _{n=1}^N\left(u_{n,2}{u}_{n-1,3}\right)\hfill \\ \hfill & +\left(-2a_3-2b_3\right)\sum _{n=1}^N\left(u_{n,2}{u}_{n,3}\right)+\left(-2c_3\right)\sum _{n=1}^N\left(u_{n,2}{u}_{n-1,3}\right)\hfill \\ \hfill & =0,\hfill \end{array}$

Substituting the above equation into Equation (25), the energy norm change rate of the whole physical space is obtained as

(33)$\begin{array}{cc}\hfill & \frac{\partial }{\partial t}{\int }_a^b{u}^{2}dx=\frac{h}{2}\frac{{\sum }_{n=1}^N\partial {\int }_{-1}^1{u}^{2}d\xi }{\partial t}\le 0.\hfill \end{array}$

Since h is positive, the above result shows that the CNNW2 scheme is linear energy stable with $\varphi =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.

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

(34)${\left({\int }_{\Omega }{u}^{2}dx\right)}^{\frac{1}{2}}$

where $\Omega$ 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.

4.1. Linear Energy Stability Test

4.1.1. 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:

(35)$u(x,0)=e^{-20x^2}.$

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.

4.1.2. 2D Linear Convection Equation

Consider the 2D linear convection equation of

(36)$\frac{\partial u}{\partial t}+\frac{\partial u}{\partial x}+\frac{\partial u}{\partial y}=0.$

The computational domain is defined in $[-5,5]\times [-5,5]$. The polynomial order is 3. $60\times 60$ uniform quadrilateral grids are used for the computations. Using upwind flux and periodic boundary conditions, the compution time t is 20 with time steps calculated by CFL = 0.1. This case uses the TVB indicator with adjustable parameter $M=0.6$. The initial condition is set as follows:

(37)$u(x,y,0)=e^{-\left(x^2+y^2\right)}.$

Figure 7, Figure 8 and Figure 9 respectively show the numerical solutions of the CPR $(Q>P)$, CNNW2−L and HCCS−L2 schemes along with the changes in the $L^2$ energy norm over time. The results are consistent with the analysis based on the 1D convection equation. The three schemes are energy-stable, and the CPR $(Q>P)$ method is the highest in accuracy. Figure 9c shows the distribution of troubled cells for the HCCS−L2 scheme, where troubled cells are marked with the red color.

4.2. Nonlinear Stability Test

4.2.1. 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 $({\rho }_{\infty },u_{\infty },v_{\infty },p_{\infty })=(1.0,1.0,0.0,1.0)$, and $T_{\infty }=p_{\infty }/{\rho }_{\infty }$. The initial conditions of the vortex are set as follows:

(38)$\begin{array}{cc}\hfill & \Delta u=-\left(y-y_c\right)\frac{\epsilon }{2\pi }exp\left(\frac{1-r^2}{2}\right),\hfill \\ \hfill & \Delta v=\left(x-x_c\right)\frac{\epsilon }{2\pi }exp\left(\frac{1-r^2}{2}\right),\hfill \\ \hfill & \Delta T=-\frac{(\gamma -1)\epsilon }{8\gamma {\pi }^2}exp\left(1-r^2\right),\hfill \end{array}$

where $r=\sqrt{{(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:

(39)$(\rho ,u,v,p)=[{(T_{\infty }+\Delta T)}^{\frac{\gamma }{\gamma -1}},u_{\infty }+\Delta u,v_{\infty }+\Delta v,{(T_{\infty }+\Delta T)}^{\frac{1}{\gamma -1}}]$

The computational domain is $[-10,10]\times [-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.

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

4.2.3. 2D Kelvin-Helmholtz Instability Problem

The initial conditions are set as follows:

(40)$\begin{array}{cc}\hfill \rho (x,y)& =\frac{1}{2}+\frac{3}{4}B,p(x,y)=1,\hfill \\ \hfill u(x,y)& =\frac{1}{2}(B-1),v(x,y)=\frac{1}{10}sin\left(2\pi x\right),\hfill \end{array}$

where $B=tanh(15y+7.5)-tanh(15y-7.5)$. The computational domain is defined in $[-1,1]\times [-1,1]$. The polynomial order is 7, and the number of grid cells is $64\times 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).

4.2.4. 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 $M{a}_{\infty }=0.8$ and angle of attack $\alpha =1.{25}^{\circ }$. 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.

5. 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 $\varphi$ = 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. It also shows that the CPR $(Q>P)$ method with subcell nonlinear CNNW2 limiting is nonlinearly stable.

Author Contributions

Methodology, writing Z.-G.Y.; Validation, F.J.; Methodology, writing, validation, R.L.; Supervision, methodology, writing, H.Z.; Supervision, X.F. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

The authors thank the editor and referees for their valuable comments and suggestions that helped us improve the results of this paper.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Table

Figure 1. A linear transformation between a physical cell and its computational cell.

Figure 2. Location of solution points and flux points of the original CPR method (p4).

Figure 3. Location of solution points and flux points of the CPR [Forumla omitted. See PDF.] method (p4).

Figure 4. The position relationship between the solution points (red circles) and the flux points (blue squares) in the CNNW2 scheme. Here, [Forumla omitted. See PDF.] denotes the jth flux point, j = 1, 2, 3, 4.

Figure 5. Numerical solutions and [Forumla omitted. See PDF.] energy of the CPR [Forumla omitted. See PDF.], CNNW2−L and HCCS−L2 schemes.

Figure 6. Distributions of troubled cells corresponding to the HCCS−L2 scheme.

Figure 7. Numerical solution and [Forumla omitted. See PDF.] energy of the CPR [Forumla omitted. See PDF.] scheme.

Figure 8. Numerical solution and [Forumla omitted. See PDF.] energy of the CNNW2−L scheme.

Figure 9. Numerical solution, [Forumla omitted. See PDF.] energy and troubled cell distribution of the HCCS−L2 scheme.

Figure 9. Numerical solution, [Forumla omitted. See PDF.] energy and troubled cell distribution of the HCCS−L2 scheme.

Figure 10. Density contours of the CPR [Forumla omitted. See PDF.] and CPR [Forumla omitted. See PDF.] schemes.

Figure 11. Density contour and troubled cell distributions in the HCCS−N2 scheme.

Figure 12. Density contour and troubled cell distributions in the HCCS−N2 scheme.

Figure 13. NACA0012 airfoil grid.

Figure 14. Density contour, pressure contour, and Mach number contour of the NACA0012 airfoil.

Figure 15. Pressure coefficient distribution on the airfoil surface.

Table 1

Errors of three schemes.

Error	Mesh	CPR (Q > P)	HCCS−L2	HCCS−N2
$L_1$ error	20 × 20	1.01 × 10${}^{-5}$	4.87 × 10${}^{-5}$	6.73 × 10${}^{-5}$
	40 × 40	7.07 × 10${}^{-7}$	1.21 × 10${}^{-5}$	1.54 × 10${}^{-5}$
	80 × 80	4.05 × 10${}^{-8}$	1.80 × 10${}^{-5}$	2.50 × 10${}^{-6}$
$L_{\infty }$ error	20 × 20	6.52 × 10${}^{-4}$	5.71 × 10${}^{-3}$	9.07 × 10${}^{-3}$
	40 × 40	7.55 × 10${}^{-5}$	2.10 × 10${}^{-3}$	3.67 × 10${}^{-3}$
	80 × 80	4.21 × 10${}^{-6}$	6.32 × 10${}^{-4}$	1.44 × 10${}^{-3}$

Appendix A

We prove the energy stability property of CNNW2 with $\varphi$ = 0 in Equation (21). The energy norm change rate of the nth computing cell is (A1)$$\begin{array}{cc}\hfill \frac{\frac{1}{2}\partial {\int }_{-1}^1{u}^{2}d\xi }{\partial t}& =u_{n,1}\frac{\partial u_{n,1}}{\partial t}{w}_1+u_{n,2}\frac{\partial u_{n,2}}{\partial t}{w}_2+u_{n,3}\frac{\partial u_{n,3}}{\partial t}{w}_3\hfill \\ \hfill & =(-1)u_{n,1}^2+(-1)u_{n,2}^2+(-1)u_{n,3}^2\hfill \\ \hfill & +u_{n-1,3}{u}_{n,1}+u_{n,1}{u}_{n,2}+u_{n,2}{u}_{n,3}.\hfill \end{array}$$ The energy norm change rate of the whole computational space is (A2)$$\begin{array}{cc}\hfill \frac{\frac{1}{2}{\sum }_{n=1}^N\partial {\int }_{-1}^1{u}^{2}d\xi }{\partial t}& =(-1)\sum _{n=1}^N{u}_{n,1}^2+(-1)\sum _{n=1}^N{u}_{n,2}^2+(-1)\sum _{n=1}^N{u}_{n,3}^2\hfill \\ \hfill & +\sum _{n=1}^N{u}_{n-1,3}{u}_{n,1}+\sum _{n=1}^N{u}_{n,1}{u}_{n,2}+\sum _{n=1}^N{u}_{n,2}{u}_{n,3}\hfill \\ \hfill & =(-\frac{1}{2})\sum _{n=1}^N\left(u_{n,1}^2+u_{n,2}^2\right)+\sum _{n=1}^N\left(u_{n,1}{u}_{n,2}\right)\hfill \\ \hfill & +(-\frac{1}{2})\sum _{n=1}^N\left(u_{n,1}^2+u_{n,3}^2\right)+\sum _{n=1}^N\left(u_{n,1}{u}_{n-1,3}\right)\hfill \\ \hfill & +(-\frac{1}{2})\sum _{n=1}^N\left(u_{n,2}^2+u_{n,3}^2\right)+\sum _{n=1}^N\left(u_{n,2}{u}_{n,3}\right)\hfill \\ \hfill & \le \left(-1\right)\sum _{n=1}^N\left(u_{n,1}{u}_{n,2}\right)+\sum _{n=1}^N\left(u_{n,1}{u}_{n,2}\right)\hfill \\ \hfill & +\left(-1\right)\sum _{n=1}^N\left(u_{n,1}{u}_{n-1,3}\right)+\sum _{n=1}^N\left(u_{n,1}{u}_{n-1,3}\right)\hfill \\ \hfill & +\left(-1\right)\sum _{n=1}^N\left(u_{n,2}{u}_{n,3}\right)+\sum _{n=1}^N\left(u_{n,2}{u}_{n,3}\right)\hfill \\ \hfill & =0.\hfill \end{array}$$ By substituting the above equation into Equation (25), the energy norm change rate of the whole physical space can be obtained as (A3)$$\begin{array}{cc}\hfill \frac{\partial }{\partial t}{\int }_a^b{u}^{2}dx& =\frac{h}{2}\frac{{\sum }_{n=1}^N\partial {\int }_{-1}^1{u}^{2}d\xi }{\partial t}\le 0.\hfill \end{array}$$ Therefore, the CNNW2 scheme is linear energy stable when $\varphi$ = 0.

References

1. Castonguay, P.; Vincent, P.E.; Jameson, A. A new class of high-order energy stable flux reconstruction schemes for triangular elements. J. Sci. Comput.; 2012; 51, pp. 224-256. [DOI: https://dx.doi.org/10.1007/s10915-011-9505-3]

2. Huynh, H.T. A Flux Reconstruction Approach to High-Order Schemes Including Discontinuous Galerkin Methods. Proceedings of the 18th AIAA Computational Fluid Dynamics Conference; Miami, FL, USA, 25–28 June 2007.

3. Cockburn, B.; Karniadakis, G.E.; Shu, C.W. The Development of Discontinuous Galerkin Methods; Springer: Berlin/Heidelberg, Germany, 2000.

4. Cockburn, B.; Shu, C.W. TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws. II. General framework. Math. Comput.; 1989; 52, pp. 411-435.

5. Johnson, C.; Pitkranta, J. An analysis of the discontinuous Galerkin method for a scalar hyperbolic equation. Math. Comput.; 1986; 46, pp. 1-26. [DOI: https://dx.doi.org/10.1090/S0025-5718-1986-0815828-4]

6. Cockburn, B.; Shu, C.W. The Runge-Kutta local projection-discontinuous-Galerkin finite element method for scalar conservation laws. ESAIM Math. Model. Numer. Anal.; 1991; 25, pp. 337-361. [DOI: https://dx.doi.org/10.1051/m2an/1991250303371]

7. Liu, Y.; Vinokur, M.; Wang, Z.J. Discontinuous spectral difference method for conservation laws on unstructured grids. Computational Fluid Dynamics 2004, Proceedings of the Third International Conference on Computational Fluid Dynamics, ICCFD3, Toronto, ON, Canada, 12–16 July 2004; Springer: Berlin/Heidelberg, Germany, 2006; pp. 449-454.

8. Nakatsu, T.; Ichiyama, S.; Hiratake, J.; Saldanha, A.; Kobashi, N.; Sakata, K.; Kato, H. Structural basis for the spectral difference in luciferase bioluminescence. Nature; 2006; 440, pp. 372-376. [DOI: https://dx.doi.org/10.1038/nature04542] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/16541080]

9. Kopriva, D.A. A Conservative Staggered-Grid Chebyshev Multidomain Method for Compressible Flows. II. A Semi-Structured Method. J. Comput. Phys.; 1996; 128, pp. 475-488. [DOI: https://dx.doi.org/10.1006/jcph.1996.0225]

10. Kopriva, D.A.; Kolias, J.H. A Conservative Staggered-Grid Chebyshev Multidomain Method for Compressible Flows. J. Comput. Phys.; 1996; 125, pp. 244-261. [DOI: https://dx.doi.org/10.1006/jcph.1996.0091]

11. Wang, Z.J.; Gao, H. A unifying lifting collocation penalty formulation including the discontinuous Galerkin, spectral volume/difference methods for conservation laws on mixed grids. J. Comput. Phys.; 2009; 228, pp. 8161-8186. [DOI: https://dx.doi.org/10.1016/j.jcp.2009.07.036]

12. Wang, Z.J.; Gao, H. A unifying lifting collocation penalty formulation for the Euler equations on mixed grids. Proceedings of the 47th AIAA Aerospace Sciences Meeting Including The New Horizons Forum and Aerospace Exposition; Orlando, FL, USA, 5–8 January 2009; 401.

13. Vincent, P.E.; Castonguay, P.; Jameson, A. A New Class of High-Order Energy Stable Flux Reconstruction Schemes. J. Sci. Comput.; 2011; 47, pp. 50-72. [DOI: https://dx.doi.org/10.1007/s10915-010-9420-z]

14. Jameson, A.; Vincent, P.E.; Castonguay, P. On the Non-linear Stability of Flux Reconstruction Schemes. J. Sci. Comput.; 2012; 50, pp. 434-445. [DOI: https://dx.doi.org/10.1007/s10915-011-9490-6]

15. Sheshadri, A.; Jameson, A. On the stability of the flux reconstruction schemes on quadrilateral elements for the linear advection equation. J. Sci. Comput.; 2016; 67, pp. 769-790. [DOI: https://dx.doi.org/10.1007/s10915-015-0102-8]

16. Spiegel, S.C.; Huynh, H.; Debonis, J.R. De-Aliasing through Over-Integration Applied to the Flux Reconstruction and Discontinuous Galerkin Methods. Proceedings of the 22nd AIAA Computational Fluid Dynamics Conference; Dallas, TX, USA, 22–26 June 2015.

17. Sonntag, M.; Munz, C.D. Shock capturing for discontinuous Galerkin methods using finite volume subcells. Proceedings of the Finite Volumes for Complex Applications VII-Elliptic, Parabolic and Hyperbolic Problems: FVCA 7; Berlin, Germany, 15–20 June 2014; Springer International Publishing: Berlin/Heidelberg, Germany, 2014; pp. 945-953.

18. Dumbser, M.; Zanotti, O.; Loubère, R.; Diot, S. A posteriori subcell limiting of the discontinuous Galerkin finite element method for hyperbolic conservation laws. J. Comput. Phys.; 2014; 278, pp. 47-75. [DOI: https://dx.doi.org/10.1016/j.jcp.2014.08.009]

19. Dumbser, M.; Loubère, R. A simple robust and accurate a posteriori sub-cell finite volume limiter for the discontinuous Galerkin method on unstructured meshes. J. Comput. Phys.; 2016; 319, pp. 163-199. [DOI: https://dx.doi.org/10.1016/j.jcp.2016.05.002]

20. Kochi, S.R.S.P.; Ramakrishna, M. A compact subcell WENO limiting strategy using immediate neighbours for Runge-Kutta discontinuous Galerkin methods. Int. J. Comput. Math.; 2021; 98, pp. 608-626. [DOI: https://dx.doi.org/10.1080/00207160.2020.1770234]

21. Zhu, H.; Liu, H.; Yan, Z.-G.; Shi, G.; Deng, X. A priori subcell limiting based on compact nonuniform nonlinear weighted schemes of high-order CPR method for hyperbolic conservation laws. Comput. Fluids; 2022; 241, 105456. [DOI: https://dx.doi.org/10.1016/j.compfluid.2022.105456]

22. Shi, G.Q.; Yan, Z.G.; Zhu, H.J. Troubled cell detection in subcell limiting for high-order CPR method on unstructured meshes. Acta Aerodyn. Sin.; 2023; 41, pp. 38-52.

23. Liu, R.; Zhu, H.J.; Jia, F.R. An energy stable flux reconstruction method based on over-integration. Math. Numer. Sin.; 2023; in press

24. Yamaleev, N.K.; Carpenter, M.H. Third-order Energy Stable WENO scheme. J. Comput. Phys.; 2009; 228, pp. 3025-3047. [DOI: https://dx.doi.org/10.1016/j.jcp.2009.01.011]

25. Cockburn, B.; Shu, C.W. The Runge–Kutta discontinuous Galerkin method for conservation laws V: Multidimensional systems. J. Comput. Phys.; 1998; 141, pp. 199-224. [DOI: https://dx.doi.org/10.1006/jcph.1998.5892]

26. Shi, G.Q.; Zhu, H.; Yan, Z.G. A priori subcell limiting approach for the FR/CPR method on unstructured meshes. arXiv; 2021; arXiv: 2108.03446[DOI: https://dx.doi.org/10.4208/cicp.OA-2021-0165]

27. Sheshadri, A.; Jameson, A. An analysis of stability of the flux reconstruction formulation on quadrilateral elements for the linear advection–diffusion equation. J. Sci. Comput.; 2018; 74, pp. 1757-1785. [DOI: https://dx.doi.org/10.1007/s10915-017-0513-9]

28. Cox, C.; Trojak, W.; Dzanic, T.; Witherden, F.D.; Jameson, A. Accuracy, stability, and performance comparison between the spectral difference and flux reconstruction schemes. Comput. Fluids; 2021; 221, 104922. [DOI: https://dx.doi.org/10.1016/j.compfluid.2021.104922]

29. Rueda-Ramírez, A.M.; Gassner, G.J. A subcell finite volume positivity-preserving limiter for DGSEM discretizations of the Euler equations. arXiv; 2021; arXiv: 2102.06017

30. Park, J.S.; Yoon, S.H.; Kim, C. Multi-dimensional limiting process for hyperbolic conservation laws on unstructured grids. J. Comput. Phys.; 2010; 229, pp. 788-812. [DOI: https://dx.doi.org/10.1016/j.jcp.2009.10.011]

31. Gottlieb, S.; Shu, C.W.; Tadmor, E. Strong stability-preserving high-order time discretization methods. SIAM Rev.; 2001; 43, pp. 89-112. [DOI: https://dx.doi.org/10.1137/S003614450036757X]