Numerical stability analysis of shock-capturing methods for strong shocks I: second-order MUSCL schemes

We consider a compressible flow governed by the two-dimensional Euler equations written in integral form as

where $\partial\Omega$ denotes boundaries of the control volume $\Omega$ and ${\mathbf{F}}$ is the flux component normal to the boundary $\partial\Omega$. The vectors of conservative variables and flux are given respectively by

where $\rho$, $e$, and $p$ represent density, specific total energy, and pressure respectively, and ${\bf{u}}=\left({u,v}\right)$ is the flow velocity. The directed velocity, $q=un_{x}+vn_{y}$, is the component of velocity acting in the $\mathbf{n}$ direction, where $\mathbf{n}={\left[{{n_{x}},\;{n_{y}}}\right]^{T}}$ is the outward unit vector normal to the surface element ${\rm{d}}S$. The equation of state is in the form

where $\gamma$ is the specific heat ratio.

We consider discretizing the system (1) with a cell-centered finite-volume method over a two-dimensional domain subdivided into some structured quadrilateral cells. A particular control volume $\Omega_{i,j}$ is shown in Fig.1, and the semi-discrete finite volume scheme over it can be written as

where ${\mathbf{U}}_{i,j}$ denotes the average of U on $\Omega_{i,j}$. $\left|\Omega_{i,j}\right|$ is the volume of $\Omega_{i,j}$ and $\mathcal{L}_{k}$ stands for the length of the cell interface. The flux $\mathbf{F}_{k}$ is the calculated numerical flux that is supposed to be constant along the individual cell interface. In the current study, approximate Riemann solvers are implemented to determine the numerical flux $\mathbf{F}_{k}$ at each cell interface.

In our study, we devote our efforts to investigating the shock instability of second-order finite-volume schemes, particularly the MUSCL approach [53, 54]. As seen in Fig.1, considering the interface between $\Omega_{i,j}$ and $\Omega_{i+1,j}$, we can get [55, 56]

where $\mathbf{U}_{i+1/2,j}^{L}$ and $\mathbf{U}_{i+1/2,j}^{R}$ are the variables on the left and right sides of the interface and $\Psi_{i+1/2,j}^{L/R}$ is called limiter. The limiter is a function of the consecutive gradients, which can be written as

where

where $\Delta_{+}$ and $\Delta_{-}$ are the forward and backward difference operators. When $\Delta_{-}\mathbf{U}_{i,j}=0$ and $\Delta_{+}\mathbf{U}_{i+1,j}=0$, $\Psi_{i+1/2,j}^{L}$ and $\Psi_{i+1/2,j}^{R}$ are set to be zero in this paper, respectively. And note that the corresponding scheme becomes first-order accurate if $\Psi_{i+1/2,j}^{L/R}=0$. In the current work, four popular limiters including superbee [57], van Leer [58], van Albada [59], and minmod [60] are employed. Their functions can be written as

As shown in Fig.2, four limiter functions all lie in the shaded region, which is the second-order TVD region [61]. Generally, from the upper boundary of the shadow region to the lower, the limiters become more dissipative [62].

Fig.3 shows the computational grid for the steady shock problem. The initial conditions are prescribed for left ($L:i\leq 5$) and right ($R:i\geq 7$) following the Rankine-Hugoniot conditions across the normal shock as

with

where $M_{0}$ is the upstream Mach number. The intermediate states within the shock ($M:i=6$) are set according to the Hugoniot curve [33]

with

where $\varepsilon\in[0,1]$ is a weighting average that describes the initial state of the internal cell and is called shock position here. The boundary conditions in all directions are imposed through ghost cells. The inflow boundary conditions are set to freestream values, while in order to fix the shock at the same position, the outflow boundary conditions are determined by setting the mass flux at the ghost cells as [35, 63]

Meanwhile, other values are simply extrapolated. The upper and lower boundaries are set as periodic conditions. In order to explore the stability of numerical methods for strong shock capturing, a convergent and stable solution of the 2D planar steady shock problem should be obtained for the linearized stability analysis. Here, the numerical experiment setup presented above is used to carry out numerical tests for the convergent and stable solutions. However, such a convergent and stable flow field is challenging to be obtained, considering that the numerical shock instability may occur in certain cases. In the following section, a simple and effective strategy to enforce a convergent and stable flow field is introduced.

In the current study, the 2D flow field for stability analysis is initialized by projecting the steady flow field from 1D computation onto the 2D domain. This initialization method is also employed in [64, 32]. In this way, when a small initial perturbation is introduced into the 2D flow field, it will be damped if the scheme is stable. Otherwise, the initial perturbation will increase, leading to the instability. As a result, we can investigate the stability of different shock-capturing methods by tracking the evolution of the perturbation error. In this paper, the small initial perturbation is set as $\delta=10^{-7}$.

Unfortunately, it has been demonstrated that the shock instability or the carbuncle phenomenon also occurs in the 1D computation. When it happens, the computation fails to converge and a stable solution cannot be obtained. As a result, it is hard to initialize the 2D flow field. The 1D shock instability or carbuncle usually occurs when the Riemann solver that produces stationary discrete shocks with a single interior point is used, such as Roe under the shock position $\varepsilon={0.1,0.2,0.3}$ [35, 38, 39, 28]. Xie et al.[21] find that the shock can be stabilized by enforcing the consistency of the mass flux across the normal shock, as a result of which, the stable solution can be obtained. In this paper, this method proposed in [21] is used to obtain a stable result in the 1D domain, which can be written as

where the subscripts $L$ and $R$ denote the cells just in front of and after the shock structure, respectively. The modification (14) is used in this paper for Roe Riemann solver when $\varepsilon=0.1,0.2,0.3$, if not mentioned specifically.

The typical carbuncle phenomenon in 2D steady normal shock problem is shown in Fig.4. As shown, the carbuncle phenomenon is characterized by several convex or concave wedge-shaped shock profiles distributed along the transversal direction in a staggered manner, and vortices appear behind the shock profile [21].

An initial perturbation of size $10^{-7}$ is added on each cell at the beginning of the computation, and Fig.5 shows the evolution of the maximal perturbation error. In the current work, we take the norm $\|v\|_{\infty}(t)$ of the transverse velocity as the indicator of perturbation error since it is known that the exact solution for mean flows is $v=0$ [32, 43]. Fig.5 and Fig.6 are computed with the condition of $M_{0}=20$ and $\varepsilon=0.1$, where the computation is more prone to be unstable. The second-order scheme with Roe solver and van Albada limiter is employed. The computational grid is 11$\times$11 Cartesian grid. Third-order TVD Runge Kutta discretization [65] with $\text{CFL}=0.1$ is used throughout the paper, if not mentioned specifically. As shown in Fig.5, according to the evolution of the perturbation error, there are three stages in the development of the shock instability [32]:

• 1.

  In the time interval from 0 to 7.5, which can be called the initial transition region, the perturbation error almost remains static. This might be due to the damping of stable modes [32]. The flow field in this stage is consistent with the initial one and remains steady, as seen in Fig.6(a).

• 2.

  The evolution of $\|v\|_{\infty}(t)$ exhibits an exponential increase in the time interval from 7.5 to 24. Thus, this region is called the exponential growth region. The evolution of the perturbation error at this stage satisfies the exponential law

  $$\|v\|_{\infty}(t)=v_{0}\mathrm{e}^{\lambda_{num}\left(t-t_{0}\right)},$$

  (15)

  where $v_{0}$, $\lambda_{num}$, and $t_{0}$ can be easily obtained from Fig.5, and $\lambda_{num}$ is called the temporal error growth rate in the present work. The flow field becomes unstable during the exponential growth stage, as shown in Fig.6(b) and Fig.6(c), and the shock profile starts to deform.

• 3.

  After time 24, the perturbation error is sufficiently large, nonlinear effects become significant, and the exponential growth is minimized. As seen in Fig.5, the error may fluctuate at a high level during this stage. The flow field, meanwhile, is unphysical. The instability is more severe, as shown in Fig.6(d), and the shock is severely distorted.

As shown in Fig.5, it can be found that the exponential growth stage plays a vital role in the evolution of the perturbation error, determining whether the error will increase or decrease, further stable or unstable, and how quickly it will develop towards instability. Therefore, in the current study we concentrate on the exponential growth stage and relate the shock instability to the temporal error growth rate $\lambda_{num}$.

Fig.7 shows the density contours computed by first and second-order schemes at the same computing time. As shown, the shock profile computed by the second-order scheme distorts more severely than the first-order case. And there are vortices appearing behind the shock in the flow field computed by the second-order scheme. Fig.8 is the comparison of the perturbation error, which also shows that the perturbation error produced by the second-order scheme steps into the exponential growth region earlier and increases more quickly. It is also larger at the nonlinear region than that of the first-order scheme. Both Fig.7 and Fig.8 suggest that the shock instability of the second-order scheme is more severe than the first-order scheme. And in the current study, the shock instability of the second-order scheme will be studied in more detail.

The matrix stability analysis method proposed by Dumbser et al.[32] has been well demonstrated to be a useful tool for assessing the shock instability problem of shock-capturing methods. However, such a stability analysis method is only applicable to first-order schemes, which limits its application in analyzing the stability of second and high-order schemes. In the current study, we extend this method for second-order finite-volume schemes, which is presented as follows.

For the stability analysis of a steady field, we assume that

where $\mathbf{U}_{i,j}^{0}$ represents the steady mean value and $\delta\mathbf{U}_{i,j}$ is small numerical random perturbation. Assume that the random perturbation has no effect on the limiter function. Thus, the limiter is the function of the steady mean value

As a consequence, the matrix stability analysis method can also be applicable for non-differentiable limiters, such as mimod. And the results in section 4 and 5 show that such assumption is appropriate. Substituting (16) into (5), we can get

where $\mathbf{U}_{i+1/2,j}^{L}$ and $\mathbf{U}_{i+1/2,j}^{R}$ are the variables on the left and right sides of the interface between $\Omega_{i,j}$ and $\Omega_{i+1,j}$. E denotes the $4\times 4$ unit matrix. Since the numerical flux $\mathbf{F}_{i+1/2,j}$ is the function of $\mathbf{U}_{i+1/2,j}^{L}$ and $\mathbf{U}_{i+1/2,j}^{R}$, $\mathbf{F}_{i+1/2,j}$ can be linearized around the steady mean value as

To make the expression brief, the following variables are introduced

Then (19) can be written as

Substituting (16) and (21) into (4) and taking into account that the mean-field is steady, finally, we finally get the linear error evolution mode

where

It can be found from (22) that the evolution of perturbation error in $\Omega_{i,j}$ is influenced by the errors in cell $\Omega_{i,j}$ itself, four neighbors, and four sub-adjacent cells. However, for the first-order scheme, only the errors in $\Omega_{i,j}$ and four neighbors will affect the evolution of the error in $\Omega_{i,j}$ [32]. The difference in the linear error evolution model suggests that the shock stability of the first and second-order schemes may be significantly different.

Equation (22) holds for all cells in the computational domain and we finally get the error evolution of all cells in the computational domain

where S is called the stability matrix in the present study. When considering the evolution of the initial error only, (24) can be solved analytically. The solution is

and will remain bounded if the maximal real part of the eigenvalues of S is negative. So the stability criterion is

One should note that equation (19) is accurate only when the numerical flux is differentiable at the mean value, which is not always holding, for example, when the Roe solver is employed and the shock is exactly between two cells ($\varepsilon=0\enspace\text{or}\enspace 1$). In the current work, the cases that have numerical shock structure ($0<\varepsilon<1$) are analyzed. As a result, the numerical flux calculated by the solvers used in this paper is differentiable at the mean value. The gradients of the numerical flux functions such as $\dfrac{\partial\mathbf{F}_{i+1/2,j}}{\partial\mathbf{U}_{i+1/2,j}^{L,0}}$ can be calculated by the centered difference approximation [11]

where $\mathbf{I}_{k}$ is unit vector of which the $kth$ component is 1, and $\delta=10^{-7}$ [32, 11].

It is known that the second-order MUSCL approach (5) can be performed not only in conservative variables but also in primitive variables. The reconstruction of primitive variables is widely used in flow simulation, especially in hypersonic flow simulation, due to its lower cost and broad applicability [66, 67]. Therefore, the matrix stability analysis method in the form of primitive variables is pursued in the current study. The numerical flux $\mathbf{F}_{i+1/2,j}$ is also the function of the primitive variables

where $\mathbf{W}=(\rho,u,v,p)^{T}$ is the vector of primitive variables. Then the flux function $\mathbf{F}_{i+1/2,j}$ can be linearized as

where the variables $\eta_{i+1/2,j}^{L/R}$ and $\beta_{i+1/2,j}^{L/R}$ in (29) are as follows

Substituting (16) and (29) into (4), we can get

where $\xi$ and $\mu$ are computed by (23). And $\frac{\mathrm{d}\mathbf{U}_{i,j}}{\mathrm{d}{\mathbf{W}_{i,j}}}$ is the transformation matrix between the conservative variables and the primitive variables. It can be written as

Equation (31) is the primitive variables form of linear error evolution model in $\Omega_{i,j}$. (24) can also be written in primitive variables form as follows

Considering the evolution of initial errors only, the solution of (33) is

And the stability criterion can be seen in (26).

In section 4.1 and 4.2, we introduce the matrix stability analysis method for the second-order MUSCL scheme. By analyzing (25), we can find the maximal real part of the eigenvalues can describe the exponential growth of the perturbation error, which is $\lambda_{num}$ in (15). Since $\lambda_{num}$ can be obtained easily by numerical experiments, the reliability of the matrix stability analysis method can be quantitatively verified by comparing $\text{max}(\text{Re}(\lambda))$ and $\lambda_{num}$.

The second-order scheme with Roe solver and van Albada limiter is used to perform the quantitative validation. The conditions are $M_{0}=1,2,\cdots,20$ and $\varepsilon=0.1,0.5,0.9$. The computational grid is $11\times 11$ Cartesian grid. Fig.9 provides the comparison between $\text{max}(\text{Re}(\lambda))$ and $\lambda_{num}$, which shows a good agreement between them in both unstable region and stable region. One should note that good agreements can also be obtained by other second-order schemes with different Riemann solvers, limiter functions, and computational conditions. That confirms the reliability of the matrix stability analysis method developed in the current study.

It has been demonstrated that the Riemann solver plays an important role in triggering the shock instability [32, 35, 21, 68, 3]. In this section, we further investigate the influence of the Riemann solver on shock instability for the second-order scheme. Here, several approximate Riemann solvers with different dissipative characteristics are considered, including Roe [69, 70], HLL [7], HLLC [20], AUSM⁺ [42], and van Leer [16]. Note that in this paper, the HLL and HLLC Riemann solvers employ the wave speeds estimate proposed by Davis [71].

The stability analysis is carried out with the condition of $M_{0}=20$ and $\varepsilon=0.1$. All eigenvalues of the stability matrix for the first and second-order schemes with these Riemann solvers are shown in Fig.10. The right column shows the results of second-order schemes with the van Albada limiter, whereas the left column shows the results of first-order ones. According to Fig.10, it can be found that the sign of the maximal real parts of eigenvalues for second-order schemes is closely related to the Riemann solver and is consistent with that of the first-order schemes. It indicates that the Riemann solver still plays a key role in determining the shock stability of second-order schemes as it does for first-order cases. Moreover, As shown, eigenvalues with positive real parts are presented for Roe and HLLC Riemann solvers, which obviously, result in unstable behaviors. It is well-known that both solvers own minimal dissipation on contact and shear waves. On the contrary, the HLL and van Leer Riemann solvers always have eigenvalues with negative real parts, indicating their stability. And both solvers are known as dissipative schemes. What is more interesting is that the maximal real part of the AUSM⁺ Riemann solver is close to zero, indicating that AUSM⁺ Riemann solver is marginally stable.

To study the influence of spatial accuracy, limiter, and Mach number on shock instability, two typical Riemann solvers, Roe and HLL, are employed in this section. The same investigation can also be conducted on other solvers, such as HLLC, van Leer, and AUSM⁺, which is not shown for brevity. As seen in Fig.11, the analysis is carried out for the second-order scheme, and in order to investigate the effect of spatial accuracy, results of the first-order scheme are also presented for comparison. Four limiters (superbee, van Leer, van Albada, and minmod) are employed when considering second-order accuracy. By analyzing the results in Fig.11, the following conclusions can be obtained:

• 1.

  As shown in Fig.11, we can find that the spatial accuracy will significantly affect the shock instability. The second-order scheme with Roe solver has a larger maximal real part of the eigenvalues than the first-order scheme, which suggests that the second-order scheme with Roe solver is more prone to be unstable than the first-order scheme. However, the maximal real part of the eigenvalues of the second-order scheme with HLL solver is less than that of the first-order case, indicating that the second-order scheme with the HLL solver is more stable than the first-order scheme. The same conclusion as HLL solver can also be obtained from AUSM⁺ solver, which is marginally stable as shown in Section 5.1. As a result, we can infer that the stability of the stable and marginally stable schemes will be better as the spatial accuracy is enhanced to second-order, while the stability of the unstable schemes will be poorer.

• 2.

  For the second-order scheme, the limiter is a critical factor affecting the shock instability. The second-order scheme has the largest maximal real part of the eigenvalues when employing the superbee limiter with the lowest dissipation, while has the smallest maximal real part when using the minmod limiter with the highest dissipation. This indicates that the computation is more prone to be unstable when employing the limiter with lower dissipation. One exception is the HLL solver when $M_{0}=2$. As shown in Fig.11(b), the second-order scheme with HLL solver is more unstable when using minmod limiter, but it is still stable. Anyway, if there is a need to alleviate the shock instability problem in supersonic and hypersonic flow simulations, the limiter with higher dissipation, such as minmod, is preferred.

• 3.

  It can be seen from Fig.11 that the Mach number has a significant impact on the shock instability. When the Mach number is modest, the maximal real part of the eigenvalues increases quickly with the increasing Mach number, but when the Mach number is greater than a threshold, the growth will slow down and the maximal real part of the eigenvalues even remains unchanged. The threshold is determined by the Riemann solver, spatial accuracy, and limiter. It implies that when the Mach number is modest, its influence is evident, and the shock instability becomes more severe as the Mach number increases. However, the influence of the Mach number is limited. If the Mach number exceeds a threshold, it will do little to the shock stability. It should also be noted that the scheme with HLL solver is most stable when the Mach number is 2 or 3 (for the first-order scheme and second-order scheme with superbee limiter, the Mach number is 2, and for the others is 3.).

It has been found that the numerical shock structure plays a vital role in shock instability. In this section, we investigate the impact of the numerical shock structure for the second-order scheme using the matrix stability analysis method. Since we are more concerned about the unstable behavior, Roe solver, which is proved to be more prone to be unstable, is employed to perform the remaining analysis. Fig.12 shows the effect of the numerical shock structure. As shown, with $\varepsilon$ increasing, the maximal real part of the eigenvalues of the second-order scheme decreases, and when $\varepsilon$ exceeds a threshold, the maximal real part of the eigenvalues begins to be less than 0. The threshold is related to the Mach number, around 0.4 when $M_{0}=2$ and around 0.6 when $5<M_{0}<20$. Fig.14 is the convergent results of the 1D shock with different $\varepsilon$. According to Fig.12 and Fig.14, we can conclude that the second-order scheme will be more stable when the internal shock conditions are closer to the downstream states. The density contours with different numerical shock positions are shown in Fig.15, demonstrating the validity of the matrix stability analysis.

Fig.12 also indicates that the numerical shock structure has the same impact on both first and second-order schemes. The unstable scheme can be stable by setting the internal shock conditions closer to the downstream states, whether for the first or second-order schemes. Moreover, When $\varepsilon<0.6$, both the first and second-order schemes with Roe Riemann solver are unstable. And the maximal real part of the eigenvalues of the second-order scheme is greater than that of the first-order scheme, which suggests that the stability of the second-order scheme is worse than the first-order scheme. However, as shown in Fig.12 and Fig.13, when $\varepsilon\geq 0.6$ the maximal real parts of the eigenvalues for both first and second-order schemes are less than zero or slightly larger than zero, indicating that the computation is stable or marginally stable. Furthermore, in the region of $\varepsilon\geq 0.6$, the maximal real part of the eigenvalues of the first-order scheme is greater than that of the second-order scheme, suggesting that the stability of the second-order scheme is better than the first-order scheme, which is consistent with the conclusion in Section 5.2.

In addition to the influence of the Mach number, spatial accuracy, limiter function, Riemann solver, and the numerical shock structure, the computational grid also plays a significant role in shock instability [43, 5, 15, 32, 25]. In this section, we study further about how the grid affects the stability of the second-order scheme in capturing strong shocks from the perspectives of aspect ratio and distortion angle.

In the above sections, we devote our efforts to exploring the primary numerical characteristics of shock instability for the second-order finite-volume scheme. In this section, we will study further about an important problem about the mechanism underlying the shock instability: spatial localization of the source of instability. Dumbser et al.[32] and Chen et al.[25] investigate this problem and find that the shock instability origins from the upstream of the shock when the shock is resolved exactly between two cells. However, in actual simulation, the captured shock is rarely exactly between two cells. The effect of numerical shock structure is studied by Xie et al.[21]. Based on a series of well-designed numerical experiments, they find that the perturbations inside the shock structure are the main factors to trigger the instability. These studies concentrate on first-order schemes. In this section, we will investigate the spatial localization of the source of instability for the second-order scheme using the matrix stability analysis method.

As shown in Fig.21, in order to explore the location of the source of instability, three test cases are employed. In the first one (Fig.21(a)), all cells apart from the ghost cells are set to the upstream state $\mathbf{U}_{L}$ while the right boundary conditions are set to intermediate state $\mathbf{U}_{M}$ and downstream state $\mathbf{U}_{R}$. Since the boundaries are supposed to be error-free, so their direct contribution to error evolution vanishes in the stability matrix except from the flux-gradient [32]. This procedure prevents any error from developing in numerical shock structure and downstream region, as a result of which we can analyze the stability of the upstream region separately. Identically, in order to analyze the stability of numerical shock structure and downstream region, similar procedures are shown in Fig.21(b) and (c).

The stabilities of three test cases are investigated using the matrix stability analysis method. Fig.22 shows the distribution of the eigenvalues in the complex plane, from which we can find that the maximal real part of the eigenvalues of the numerical shock structure exceeds 0, while the others are less than 0. As a result, we can infer that in the cases of the shock on the upstream boundary and downstream boundary, the stability analysis predicts a stable behavior, while if the numerical shock structure is contained, the analysis predicts unconditional instability. Moreover, the maximal real part of eigenvalues for the downstream field is larger than that of the upstream field, which may imply that the downstream is more susceptible to being influenced by the unstable information from the numerical shock structure than the upstream. Additionally, the information of the spatial-behavior of the unstable mode is contained in the eigenvectors of the maximal eigenvalue [32]. Fig.23 depicts the unstable mode for the four primitive variables. As shown in Fig.23, the numerical shock structure is the most unstable location, and unstable information will propagate more easily downstream. Thus, according to the analysis in this section, we can conclude that the shock instability originates from the numerical shock structure. And the downstream is more susceptible to being influenced by spurious errors in the shock structure.