Further studies on numerical instabilities of Godunov-type schemes for strong shocks

TL;DR

This paper investigates the causes of numerical shock instabilities in Godunov-type schemes for strong shocks and proposes an entropy-control technique to enhance their robustness without sacrificing accuracy.

Contribution

It introduces a general entropy-control method that improves the stability of Godunov-type schemes at strong shocks by ensuring sufficient entropy production inside the shock structure.

Findings

The instability is linked to insufficient entropy production.

The proposed entropy-control technique enhances robustness.

Numerical tests confirm improved accuracy and stability.

Abstract

In this paper, continuous research is undertaken to explore the underlying mechanism of numerical shock instabilities of Godunov-type schemes for strong shocks. By conducting dissipation analysis of Godunov-type schemes and a sequence of numerical experiments, we are able to clarify that the instability may be attributed to insufficient entropy production inside the numerical shock structure. As a result, a general entropy-control technique for improving the robustness of various Godunov-type schemes at strong shocks is developed. It plays a part in guaranteeing that enough entropy is produced inside the numerical shock structure. Furthermore, such a modified approach does not introduce any additional numerical dissipation on linear degenerate waves to suppress the shock instability. Numerical results that are obtained for various test cases indicate that the proposed methods have a good performance in terms of accuracy and robustness.

Tables

Table 1. Table 1: 1D test results of Godunov-type flux functions (S: stable and converged, U: unstable).

Riemann solver	Test problem	$\epsilon =0.0$	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
HLLEM	1D	U	U	U	U	S	S	S	S	S	S
HLLE	1D	U	U	U	U	S	S	S	S	S	S

Equations

$$\frac{\partial }{\partial t}{\int }_{\Omega }\mathbf{U}d\Omega +{\oint }_{\partial \Omega }\mathbf{F}dS=0,$$

$$\mathbf{U}=\left[{\begin{array}[]{{c}}\rho\\ {\rho u}\\ {\rho v}\\ {\rho e}\end{array}}\right],\quad\mathbf{F}=\left[{\begin{array}[]{{c}}\rho q\\ {\rho uq+pn_{x}}\\ {\rho vq+pn_{y}}\\ {\left(\rho e+p\right)q}\end{array}}\right],$$

$$p=\left(\gamma -1\right)\rho \left[e-\frac{1}{2}\left(u^2+v^2\right)\right],$$

$$\frac{\mathrm{d}{\mathbf{U}}_i}{\mathrm{d}t}+\frac{1}{\left|{\Omega }_i\right|}\sum _{{\Gamma }_{ij}\subset \partial {\Omega }_i}\left|{\Gamma }_{ij}\right|{\mathbf{F}}_{ij}=0.$$

$$\omega (x/t;{\mathbf{U}}_L,{\mathbf{U}}_R)=\{\begin{array}{ll}{\mathbf{U}}_L& \text{x/t<SL}\\ {\mathbf{U}}^{\ast }& \text{SL<x/t<SR}\\ {\mathbf{U}}_R& \text{SR<x/t}\end{array},$$

$${\int }_{-\Delta x/2}^{\Delta x/2}\omega (x/t;{\mathbf{U}}_L,{\mathbf{U}}_R)dx=\Delta x/2\left({\mathbf{U}}_L+{\mathbf{U}}_R\right)-\Delta t\mathbf{F}\left({\mathbf{U}}_R\right)+\Delta t\mathbf{F}\left({\mathbf{U}}_L\right)$$

$${\mathbf{U}}_i^{n+1}=\frac{1}{\Delta x}{\int }_0^{\Delta x/2}{\omega }_{i-1/2}^n\left(x/\Delta t\right)dx+\frac{1}{\Delta x}{\int }_{-\Delta x/2}^0{\omega }_{i+1/2}^n\left(x/\Delta t\right)dx.$$

$${\mathbf{U}}_i^{n+1}={\mathbf{U}}_i^n-\frac{\Delta t}{\Delta x}\left({\mathbf{F}}_{i+1/2}^n-{\mathbf{F}}_{i-1/2}^n\right)$$

$${\mathbf{F}}_{\mathrm{HLL}}=\{\begin{array}{ll}\mathbf{F}\left({\mathbf{U}}_L\right)& \text{SL>0}\\ \frac{S_R\mathbf{F}\left({\mathbf{U}}_L\right)-S_L\mathbf{F}\left({\mathbf{U}}_R\right)}{S_R-S_L}+\frac{S_L{S}_R}{S_R-S_L}\left({\mathbf{U}}_R-{\mathbf{U}}_L\right)& \text{SL≤0≤SR}\\ \mathbf{F}\left({\mathbf{U}}_R\right)& \text{SR<0}\end{array}.$$

$$S_L=\min (q_L-a_L,\hat{q}-\hat{a}),S_R=\max (q_R+a_R,\hat{q}+\hat{a})$$

$$S_L=\min (q_L-a_L,q_R-a_R),S_R=\max (q_L+a_L,q_R+a_R).$$

$${\mathbf{F}}_{\mathrm{HLL}}=\frac{S_R^{+}\mathbf{F}\left({\mathbf{U}}_L\right)-S_L^{-}\mathbf{F}\left({\mathbf{U}}_R\right)}{S_R^{+}-S_L^{-}}+\frac{S_L^{-}{S}_R^{+}}{S_R^{+}-S_L^{-}}\left({\mathbf{U}}_R-{\mathbf{U}}_L\right)$$

$$\omega (x/t;{\mathbf{U}}_L,{\mathbf{U}}_R)=\{\begin{array}{ll}{\mathbf{U}}_L& \text{x/t<SL}\\ {\mathbf{U}}^{\ast }+(x-\overline{q}t)({\hat{\delta }}_2^{\ast }{\hat{\alpha }}_2{\hat{\mathbf{R}}}_2+{\hat{\delta }}_3^{\ast }{\hat{\alpha }}_3{\hat{\mathbf{R}}}_3)& \text{SL≤x/t≤SR}\\ {\mathbf{U}}_R& \text{SR<x/t}\end{array}$$

$$\overline{q}=\frac{S_L+S_R}{2}.$$

$${\mathbf{U}}_R-{\mathbf{U}}_L=\sum _{k=1}^4{\hat{\alpha }}_k{\hat{\mathbf{R}}}_k$$

$${\hat{\alpha }}_2=\Delta \rho -\frac{\Delta p}{{\hat{a}}^2},{\hat{\alpha }}_3=\hat{\rho }\left(-n_y\Delta u+n_x\Delta v\right),$$

$$\widehat{\bf{R}}_{2}={\left[{\begin{array}[]{*{20}{c}}1&{\widehat{u}}&{\widehat{v}}&{\left({{{\widehat{u}}^{2}}+{{\widehat{v}}^{2}}}\right)/2}\end{array}}\right]^{\rm{T}}},\begin{array}[]{*{20}{c}}{}\hfil&{}\hfil\end{array}\widehat{\bf{R}}_{3}={\left[{\begin{array}[]{*{20}{c}}0&{-{n_{y}}}&{{n_{x}}}&{-\widehat{u}{n_{y}}+\widehat{v}{n_{x}}}\end{array}}\right]^{\rm{T}}},$$

$${\hat{\delta }}_k^{\ast }=\frac{2}{T(S_R-S_L)}{\hat{\delta }}_k\text{for}k=2,3$$

$${\hat{\delta }}_2={\hat{\delta }}_3=\frac{\hat{a}}{\hat{a}+\left|\hat{q}\right|}.$$

$${\mathbf{F}}_{\mathrm{HLLEM}}=\frac{S_R^{+}\mathbf{F}\left({\mathbf{U}}_L\right)-S_L^{-}\mathbf{F}\left({\mathbf{U}}_R\right)}{S_R^{+}-S_L^{-}}+\frac{S_L^{-}{S}_R^{+}}{S_R^{+}-S_L^{-}}\left({\mathbf{U}}_R-{\mathbf{U}}_L-{\hat{\delta }}_2{\hat{\alpha }}_2{\hat{\mathbf{R}}}_2-{\hat{\delta }}_3{\hat{\alpha }}_3{\hat{\mathbf{R}}}_3\right).$$

$${\mathbf{U}}_L=\left(\begin{array}{c}1\\ 1\\ 0\\ \frac{1}{\gamma (\gamma -1)M_0^2}+\frac{1}{2}\end{array}\right),{\mathbf{U}}_R=\left(\begin{array}{c}f(M_0)\\ 1\\ 0\\ \frac{g(M_0)}{\gamma (\gamma -1)M_0^2}+\frac{1}{2f(M_0)}\end{array}\right),$$

$$f(M_0)={\left(\frac{2}{(\gamma +1)M_0^2}+\frac{\gamma -1}{\gamma +1}\right)}^{-1},g(M_0)=\frac{2\gamma M_0^2}{\gamma +1}-\frac{\gamma -1}{\gamma +1},$$

$$(\rho u{)}_{imax+1,j}=(\rho u{)}_0=1$$

$${\rho }_M$$

$$u_M$$

$$p_M$$

$${\alpha }_{\rho }$$

$${\alpha }_u$$

$${\alpha }_p$$

$$\omega (x/t;{\mathbf{U}}_L,{\mathbf{U}}_R)=\{\begin{array}{ll}{\mathbf{U}}_L& \text{x/t<SL}\\ {\mathbf{U}}^{\ast }+(x-\hat{u}t){\hat{\delta }}_2^{\ast }{\hat{\alpha }}_2{\hat{\mathbf{R}}}_2& \text{SL<x/t<SR}\\ {\mathbf{U}}_R& \text{SR<x/t}\end{array}$$

$${\mathbf{F}}_{\mathrm{HLLEM}}=\{\begin{array}{ll}\mathbf{F}\left({\mathbf{U}}_L\right)& \text{SL>0}\\ \frac{S_R\mathbf{F}\left({\mathbf{U}}_L\right)-S_L\mathbf{F}\left({\mathbf{U}}_R\right)}{S_R-S_L}+\frac{S_L{S}_R}{S_R-S_L}\left({\mathbf{U}}_R-{\mathbf{U}}_L-{\hat{\delta }}_2{\hat{\alpha }}_2{\hat{\mathbf{R}}}_2\right)& \text{SL≤0≤SR}\\ \mathbf{F}\left({\mathbf{U}}_R\right)& \text{SR<0}\end{array}$$

$$S_L=\hat{u}-\hat{a},S_R=\hat{u}+\hat{a}.$$

Full text

1 Introduction

Numerical shock prediction is one of the most important issues in computational fluid dynamics. Generally, the simulation of compressible fluid flows containing shocks can be implemented in two alternative approaches: shock-fitting and shock-capturing. The shock-fitting technique is able to provide very accurate solutions and has considerable advantages in efficiency in terms of computational cost. Unfortunately, though considerable progress has been achieved in the shock-fitting technique, it has always been unpopular due to clear limitations in simulating three-dimensional flow fields and complex shock patterns [1]. Nowadays, the shock-capturing prevails over the shock-fitting in practical simulations of fluid flows involving shock waves. Based on the mathematical theory of weak solutions, the shock-capturing method is able to naturally recognize shock without special treatment. The Godunov-type schemes with their clear physical interpretations are one of the most famous shock-capturing methods. However, modern Godunov-type schemes that have minimal dissipation on discontinuities are still limited in the cases where strong shocks exist. The solutions of flows involving strong shock waves are often characterized by the appearance of numerical anomalies. The most commonly observed shock anomaly is the carbuncle phenomenon.

Since its first discovery by Peery and Imlay [2], the carbuncle has been extensively studied by many researchers during the past three decades. Readers are referred to references [3, 4, 5, 6] and references therein for detail literature reviews of this subject. Although the carbuncle commonly refers to a spurious solution of blunt body calculations in which a protuberance grows ahead of the bow shock along the stagnation line [7], it has been generally considered to be a concrete manifestation of numerical shock instabilities. A number of studies have shown that these numerical shock instabilities are more prone to occur in cases where numerical methods own minimal dissipation on discontinuities and the shocks are sufficiently strong. Apart from the compressible Euler equations of gas dynamics, recent studies have shown that these numerical shock instabilities also plague numerical approximations to other hyperbolic systems of conservation laws, such as the shallow water equations [8, 9, 10] and the MHD equations [11]. Actually, the stability of numerical approximations to systems of conservation laws has not been well established [12, 13]. Currently, the entropy stability is the only rigorous theory of the numerical instability for numerical approximations to systems of conservation laws. Godunov scheme [14] and some Godunov-type schemes (e.g. the Roe scheme with Harten’s entropy fix [15], the HLLE/HLLEM schemes [16, 17] and the Osher-Solomon scheme [18]) satisfy a discrete version of the entropy conditions which ensure that the computed entropy will increase across shocks [19, 20]. The entropy conditions rescue Godunov-type schemes from admitting nonphysical expansion shocks. However, in some cases, it may not suffice to ensure the convergence of approximate solutions. It has been theoretically established that Godunov and most Godunov-type methods are still unstable for steady shock waves even for the one-dimensional case [21, 22]. Recent results in studies [23, 24, 25] even demonstrate that entropy solutions may not be unique. Consequently, numerical approximations that are designed to satisfy the entropy conditions may not converge to the entropy solutions we desire. That is the reason why authors in [26, 27, 28] suspect that carbuncles may be the entropy solutions of the Euler system and have their own physical background. Considering the failure of entropy stable schemes on solving a steady shock in one dimension, Ismail et al. [3, 29] are the first to investigate the effects of the entropy production on shock stability. They find that proper entropy should be produced inside the shock structure to avoid excessive smearing on discontinuities and spurious oscillations. With suitable amount of diffusion on entropy conservative flux [30], their entropy-consistent solver successfully eliminates the one-dimensional shock instability even though it does not remove the multidimensional shock instability, particularly the carbuncle phenomenon.

The importance of numerical shock structure to trigger the shock instability has also been emphasized by many other researchers [31, 32, 33, 34, 35]. In our previous work [6], we make efforts to clarify the strong relation between the onset of numerical shock instability and the shock structure. It has been clarified that the spatial location of the shock instability originates from the intermediate states inside the shock structure. Furthermore, it has been demonstrated that if the mass flux across the normal shock is correctly preserved during the computation, then the shock wave can be captured stably. These are useful for developing possible strategies for suppressing numerical shock instabilities. For example, the robustness of numerical schemes can be improved by introducing multidimensional dissipation to damp perturbation errors inside the numerical shock structure and then help maintain the correct mass flux across the normal shock. Such a strategy is commonly used to cure shock instabilities of low diffusion approximate Riemann solvers. For example, the hybrid technique [7, 36, 37, 38, 39] that introduces more numerical dissipation in the vicinity of shocks to suppress instabilities by hybridizing a low dissipative solver and a more dissipative one. In some rotated Riemann solvers [40, 41, 42], artificial dissipation is introduced in the vicinity of shocks to improve the robustness of shock-capturing schemes. Several authors [43, 6] resort to modifying the diffusion term of the complete Riemann solvers to include more numerical dissipation on linear degenerate waves to cure the shock instability. Moreover, some recent studies [44, 5, 45] have demonstrated that the carbuncle instability can be cured by imposing some artificial viscosity analogous to physical diffusion on shock-capturing schemes. Although introducing additional numerical dissipation or artificial viscosity has been well-demonstrated to be effective in eliminating the shock instability, such numerical viscosity will probably smear discontinuities and decrease overall accuracy. Actually, it is still a challenge to determine suitable amount of additive diffusion for low diffusion approximate Riemann solvers in the sense that it is large enough to make the resulting solvers carbuncle free and minimum enough to resolve the discontinuities as crisply as possible. Thus, it is an attractive issue on how to cure the carbuncle problem of shock-capturing schemes without resorting to additional numerical dissipation.

Kitamura et al. [46] develop a new pressure flux for AUSM-family schemes. They find that the dissipation inside the numerical shock structure must be proportional to Mach number. Hence, more numerical diffusion is included to smear the shock profile, but the accuracy on resolving contact and shear waves are maintained. In a recent work of Sangeeth and Mandal [47], a simple cure for numerical shock instability in the HLLC Riemann solver [48] is proposed. By a newly developed selective wave modification strategy, they succeed to eliminate the shock instability problem of the HLLC scheme without compromising on its contact and shear preserving ability. Zaide et al. [49, 35] notice that the shock with the unphysical intermediate states are highly sensitive to perturbations. To circumvent shock anomalies, they introduce the interpolated fluxes [50] that do not depend on intermediate states inside the shock. Numerical results demonstrate that such a technique is effective in eliminating several shock anomalies, including the carbuncle problem in one dimension. Apart from circumventing the intermediate state, another attractive approach is to limit the propagation of perturbations that are generated by the intermediate cell. As discussed in the previous study [6], the instability is triggered by the erroneous mass flux behind shocks. It is caused by the perturbation errors that are generated inside the numerical shock structure and propagate downstream. Such an observation paves the way for suppressing the shock instability from a new perspective. That is, we may control the transportation of the erroneous perturbations from the intermediate states to the downstream cells, thus suppressing the instability. This technique is more attractive because it does not introduce additional numerical dissipation on linear degenerate waves to enhance the stability. In a recent study [51], we have shown that the robustness of the Roe method [52] can be improved by controlling the pressure perturbations that are transported from the intermediate states inside the shock structure to the downstream regions. In the current study, we will continue to explore the mechanism of the numerical shock instability and provide an in-depth explanation for the rationality of this strategy. A focus is on clarifying the connection between the entropy production inside the numerical shock structure and the stability, which is inspired by the pioneer work of Ismail et al. [29]. Combining the numerical dissipation analysis of Godunov-type schemes and a sequence of numerical experiments, we are able to clarify that the numerical instability at strong shocks can be attributed to inappropriate entropy production inside the shock structure. As a result, a general entropy-control technique is developed to improve the robustness of various Godunov-type schemes at strong shocks. Different from common cures for carbuncles which rely on multidimensional dissipation, the current modified approach does not introduce any additional numerical dissipation to improve the robustness, thus maintaining the accuracy of numerical methods on discontinuities.

The outline of the rest of this paper is as follows. In section 2, governing equations of compressible flows and their related finite volume discretization are presented. Two classical HLL-type schemes are also reviewed in the same section. In section 3, we present a mechanism study of the numerical shock instability. In section 4, the shock instability problem is revisited through a linearized perturbation analysis. A general technique for improving the robustness of Godunov-type schemes at strong shocks is proposed in section 5. The accuracy and robustness of the proposed methods are tested in section 6. Section 7 contains conclusions and an outlook to future developments.

2 Governing equations and finite volume discretization

The two-dimensional Euler equations may be written in integral form as

[TABLE]

where $\partial\Omega$ denote boundaries of the control volume $\Omega$. The state vector and flux vector are defined as

[TABLE]

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 $dS$. The equation of state is in the form

[TABLE]

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. The semi-discrete finite volume scheme over a particular control volume ${\Omega_{i}}$ can be written as

[TABLE]

In the above expression, ${{\bf{U}}_{i}}$ is the cell average of ${\bf{U}}$ on ${\Omega_{i}}$ , $\left|{{\Omega_{i}}}\right|$ denotes the volume of ${\Omega_{i}}$. ${\Gamma_{ij}}$ denotes the common edge of two neighboring cells ${\Omega_{i}}$ and ${\Omega_{j}}$, and $|{\Gamma_{ij}}|$ is the length of face ${\Gamma_{ij}}$. The flux ${{\bf{F}}_{ij}}$ is the calculated numerical flux that is supposed to be constant along the individual face ${\Gamma_{ij}}$. In the following section, we give a brief review of two classical approximate Riemann solvers for determining the numerical flux ${{\bf{F}}_{ij}}$ at each cell interface, i.e., HLLE scheme [16, 17] and HLLEM scheme [16].

2.1 HLL approximate Riemann solver

The HLL Riemann solver proposed by Harten, Lax, and van Leer [53] is one of the most reliable approaches to solve the Riemann problem approximately. Different from the exact Riemann solution with a large amount of detail, the HLL solver assumes an average intermediate state $\bf{U}^{*}$ between the fastest and slowest waves, $S_{L}$, $S_{R}$. The resulting approximate Riemann solution can be written as

[TABLE]

where the approximate intermediate state ${\bf{U}}^{*}$ is defined to be consistent with the following integral form of the conservation law, i.e.,

[TABLE]

for $\Delta x/2>\Delta t\cdot\max\left\{{\left|S_{L}\right|,\left|S_{R}\right|}\right\}$. The cell averages at the next time level are then obtained after computing the approximate Riemann solvers at the cell interfaces, i.e.,

[TABLE]

The above equation can be rewritten in the conservative form as

[TABLE]

with the numerical flux function defined by

[TABLE]

To completely determine the interface flux, the wave speeds need to be estimated. Concerning the stability and positivity preserving property, Einfeldt et al. [17] proposed a way of computing the wave speeds,

[TABLE]

where $\widehat{(\cdot)}$ is Roe’s averaged variable with respect to ${\bf U}_{L}$ and ${\bf U}_{R}$, $a_{L}$ and $a_{R}$ are the sound speeds of the left and the right states. The estimate (10) can return the exact shock velocity in the case of an isolated shock [54]. Davis [55] also suggests another simpler and more diffusive estimate,

[TABLE]

Practically, the flux function (9) can be further rewritten into a unified form

[TABLE]

with the wave speeds $S_{L}^{-}=\min(S_{L},0)$ and $S_{R}^{+}=\max(S_{R},0)$.

2.2 HLLEM approximate Riemann solver

Due to the assumption of a two-wave configuration, the HLL scheme introduces excess dissipation on linear waves. Therefore, it provides poor resolution of physical features such as contact surfaces, shear waves and material interfaces. The difficulty with this method can be resolved by modifying the intermediate state in (5) through a linear distribution approach. For the HLLEM scheme, the solution of the Riemann problem is approximated as

[TABLE]

where $\overline{q}$ denotes the numerical approximation of the velocity at the contact discontinuity, it is defined as a simple arithmetic average of wavespeeds $S_{L}$ and $S_{R}$,

[TABLE]

$\widehat{\bf{R}}_{2}$ and $\widehat{\bf{R}}_{3}$ are the second and third right eigenvectors of the flux Jacobian evaluated at a Roe’s averaged state $\widehat{\bf{U}}$. $\widehat{\alpha}_{2}$ and $\widehat{\alpha}_{3}$ are the coefficients of the projection of ${{\bf{U}}_{R}}-{{\bf{U}}_{L}}$ onto these eigenvectors, i.e.,

[TABLE]

with

[TABLE]

and

[TABLE]

The parameters ${\hat{\delta}}^{*}_{2}$ and ${\hat{\delta}}^{*}_{3}$ denote anti-diffusion coefficients which play a role in controlling the amount of anti-diffusion in the linear degenerate fields,

[TABLE]

with

[TABLE]

One should notice that the Roe’s averaged velocity $\widehat{q}$ instead of $\overline{q}$ is used to calculate ${\hat{\delta}}_{2}$ and ${\hat{\delta}}_{3}$, which resolves the stationary contact discontinuity exactly [56]. The numerical scheme corresponding to the Riemann solver (13) can be written again in conservation form as that in (8) with the following flux function,

[TABLE]

One can observe that the HLLEM scheme will reduce to the HLLE scheme if the anti-diffusion coefficients ${\hat{\delta}}_{2}$ and ${\hat{\delta}}_{3}$ are eliminated.

3 Exploring the mechanism of numerical shock instability

The carbuncle phenomenon is conventionally referred to the distorted shock ahead of the blunt body in a supersonic or hypersonic flow. In order to analyze the carbuncle, it is wise to resort to a simplified test case which shares essential characteristics of the blunt body carbuncle and can be analyzed in a simplified manner. The problem considered here is the steady normal shock problems in one and two dimensions [57]. It has been well demonstrated that if a scheme encounters the instability for the steady shock, then it will also suffer from the carbuncle [33, 3, 57]. In the current section, we start by analyzing the instability characteristics of the one-dimensional steady shock problem, i.e., the one-dimensional carbuncle. It is a one-dimensional manifestation of the carbuncle which is often deemed as a multidimensional phenomenon. Attempts are made to explore the connection between the entropy production inside the numerical shock structure and the instability. To this end, we conduct a mechanism study through two interlinked approaches: an analysis of the numerical entropy production and numerical experiments. Then we generalize the approach to the multidimensional case using the dimensional splitting method at the end of this section.

3.1 One-dimensional carbuncle problem

3.1.1 Numerical experiment setup

The one-dimensional carbuncle is a steady shock problem which is a dimensionally-reduced case of its two-dimensional counterpart: the 1.5D steady normal shock test [57]. Here, we repeat the setup of this numerical test for self-containedness. Readers are suggested to refer to references [57, 6] for more detailed descriptions.

In Fig. 1, a schematic diagram is presented to show the computational grid and conditions for the steady shock problem. For the one-dimensional case, only one cell is used in the $j$ direction and $50$ equally spaced cells are used in the $i$ direction. The initial conditions are prescribed for left $(L:i\leq 12)$ and right $(R:i\geq 14)$ following the Rankine-Hugoniot conditions across the normal shock as

[TABLE]

with

[TABLE]

where ${M_{0}}$ and $\gamma$ represent the upstream Mach number and the gas specific heat ratio. The inflow boundary conditions are set to freestream values. The mass flux at the ghost cell is prescribed as

[TABLE]

in order for the mass in the whole computational domain to remain constant and to fix the shock at the same position [29]. Meanwhile, other values are simply extrapolated. The intermediate states within the shock are assumed to lie on the Hugoniot curve [34]. The internal shock conditions $(M:i=13)$ are as follows:

[TABLE]

with

[TABLE]

where $\varepsilon=0.0,\;...\;,\;0.9$ is a discrete weighting average. It describes the initial state of the internal cell and is called shock position here. We apply the first-order finite volume method with HLLE/HLLEM schemes defined in section 2 on this problem for various combinations of Mach numbers and shock positions $\varepsilon$.

3.1.2 The profile of a one-dimensional carbuncle

It was demonstrated both analytically [21, 22] and numerically [57] that Riemann solvers which are able to produce stationary discrete shocks with a single interior point are not always stable. As with the Godunov scheme [14] or the Roe scheme [52], the HLLE and HLLEM solvers with wave speeds (10) are able to resolve shocks with only one intermediate state inside the numerical shock structure. They suffer from the one-dimensional carbuncle which happens for strong enough shock, over a range of shock positions. We repeat the one-dimensional steady shock test in [57] to evaluate the behaviors of the HLLE/HLLEM schemes for a wide Mach number spectrum and different shock positions. The results are also presented again for the convenience of further discussion. In Table 1, stability results of numerical schemes for 1D steady shock problem are summarized. As expected, both schemes produce unstable results for certain shock positions $\varepsilon=0.1\sim 0.3$. Compared with the HLLEM scheme, the HLLE flux omits the contact wave in the solution of the Riemann problem, thus it contains excess numerical dissipation on contact discontinuities. However, it is demonstrated in Table 1 that such numerical dissipation is of little use against the one-dimensional carbuncle.

In Fig. 2, solutions of the HLLEM scheme for a stable shock position are presented. As shown, the computation converges to a weak solution of the Euler equations. An intermediate state exists inside the shock structure which is bounded by the left and right primitive variables. However, the solution fails to converge to a steady state for certain unstable shock positions. Fig. 3 demonstrates essential behaviors of the one-dimensional carbuncle. For the unstable case with setup in section 3.1.1, the intermediate state oscillates inside the shock and enters into an approximate limit cycle involving the shedding of spurious waves and their reflection from the downstream boundary. The HLLE scheme demonstrates similar results that are not shown for clarity. One should note that these unstable behaviors of the one-dimensional carbuncle have also been discussed in previous work by Ismail [3, 29, 58] and Zaide [49, 35] respectively in the context of the Roe scheme [52]. It is a common anomalous phenomenon for Godunov-type schemes, particularly the one that has minimal dissipation on shocks.

3.2 A heuristic discussion of numerical entropy production and shock instability

As we have demonstrated above, for Godunov-type schemes that capture steady shocks with a single interior point, the instability is inevitable. One may argue that adding the dissipation corresponding to nonlinear waves can spread the shock profile, thus suppressing the shock instability. This is the case. In our own experience, shock-capturing schemes that yield a steady shock structure with more than two internal cells are generally stable for the one-dimensional carbuncle. The list includes the Osher scheme [18], van Leer scheme [59], Roe scheme with Harten’s entropy fix [15], to name a few. Whereas, the finite discretization of the Euler equations is expected to capture the shock with minimal smearing. Thus, it is still attractive to improve the stability of Godunov-type schemes that resolve the shock with only one intermediate point at most.

In terms of continuum mechanics on which the Euler equations are established, the shock wave is regarded as a jump discontinuity. Theoretical studies have established that inviscid shock waves for the ideal gas are uniformly stable in one and multiple space dimensions [60, 61]. However, due to the nature of finite discretization, it is impossible to capture the shock with zero thickness. The captured shock occupies several cells that is often greater than the physical width. Hence, the shock instability of numerical approximations to the Euler systems may come from a loss of information between the continuous and discrete levels. Without explicit knowledge of these discretization errors, it is hard to theoretically investigate the instability of corresponding numerical methods.

However, the state of the intermediate point inside the numerical shock structure is defined and controlled by artificial dissipation introduced by the numerical scheme. Hence, the intermediate state is a direct representation of the characteristics of the numerical dissipation. Thus, it paves the way for exploring the mechanism of numerical shock instability through analyzing the characteristics of the shock structure. In our previous study [6], it has been shown that the shock will be stabilized if the intermediate state is fixed to a proper value. Such a value is obtained by the solution of the steady normal shock with an artificial condition that enforces the mass flux consistency (i.e., $(\rho u)_{L}=(\rho u)_{R}$) across the normal shock. However, such a technique cannot be implemented in a practical simulation. There is still an important question that needs to be answered. What the intermediate state should be ? Certain useful criterion, which can be used to guarantee that the intermediate state lies inside the stable spectrum, is needed.

The entropy condition plays a decisive role in the numerical prediction of the system of conservation laws. For the Euler systems, the entropy condition requires that the entropy should increase across the shock front. However, such a condition does not guarantee that the entropy is produced in right amounts. The numerical shock profile can be considered to be an approximation of a viscous shock profile governed by the Euler equations with an additional viscosity introduced by the numerical discretization. Analogous to the physical viscosity, the numerical viscosity will also provide a dissipative mechanism that transforms the energy into heat. The whole system should obey the second thermodynamics law and tends to a stable state with the maximum entropy. Without a guarantee that enough entropy is produced inside the shock, the numerical shock structure will be unstable.

Actually, the viscous shock structure of the Navier-Stokes equations may shed some light on this problem. In Fig. 4, the entropy profile for the NS shock has been plotted, it is generated by an open code provided by Nishikawa [62]. It is demonstrated that an entropy overshoot is observed inside the shock. One should note that theoretical studies [63, 64, 65] have established that planar viscous Navier-Stokes shock waves for the ideal gas are also uniformly stable in one and multiple space dimensions. Thus, the key factor to the shock instability lies on the characteristics of the dissipation introduced in the numerical procedure and the entropy production can be used to determine what the proper dissipation for stabilizing the shock is. In the following sections, we will demonstrate that if enough entropy production is guaranteed inside the shock, then the instability can be successfully eliminated. More entropy production will make the solution converge to a stable state more quickly.

3.3 Dissipative mechanism in Godunov-type schemes

As discussed in the above section, the numerical shock instability may be attributed to insufficient entropy production inside the numerical shock structure. To increase the entropy production around shocks, there are two alternative approaches. One is to impose more numerical dissipation on nonlinear waves to smear the shock profile, thus more entropy production can be obtained. The other is to reduce the numerical dissipation on entropy waves in the vicinity of shocks, then a larger entropy state will be produced inside the numerical shock structure. In the current study, we prefer the latter one and try to avoid adding additional dissipation corresponding to the nonlinear wave to stabilize shocks. There are at least two reasons to do that. Firstly, more numerical dissipation is harmful for simulation that needs high accuracy. Secondly, the dissipation corresponding to nonlinear waves helps little in eliminating the multidimensional shock instability. Actually, the ultimate goal of this study is to stabilize shocks without introducing any additional numerical dissipation corresponding to whether nonlinear waves or linear waves to the target Riemann solver. To facilitate a further investigation of the relation between entropy production and shock stability, we first conduct an analysis of numerical dissipation in the HLL-family schemes. And then, two modified schemes are constructed by modifying the dissipative terms in the target HLLEM scheme to control the entropy production. An essential role is played by the modified equation associated with these discrete schemes, which is found to be relevant even for solutions containing shock waves.

3.3.1 Numerical dissipation analysis of HLL-family schemes

In this section, the dissipative mechanisms of Godunov-type schemes particularly the HLL-family schemes are carefully analyzed. First, we consider the approximate solution of the Riemann problem for the HLLEM flux defined in Eq. (13). It can be written as

[TABLE]

where the dissipative term corresponding to shear waves has been omitted in the one-dimensional case. The resulting flux function in the corresponding conservation form (8) can be obtained,

[TABLE]

where the following simplified wave speed estimates are used for facilitating further analysis,

[TABLE]

It is justified for the one-dimensional steady shock problem considered here. One should note that the HLLEM scheme (27) with the wave speed estimates (28) is identical to the Roe scheme [52], so results of the current analysis can also be applied to the Roe scheme. The above flux function for the subsonic case may also be written in the form

[TABLE]

According to Einfeldt[16], the numerical viscosity matrix ${\bf{Q}}_{i+1/2}$ can be rewritten by

[TABLE]

with

[TABLE]

Here, $\widehat{\bf{A}}$ is Roe Jacobian matrix and $\widehat{\bf{R}}=(\widehat{\bf{R}}_{1},\widehat{\bf{R}}_{2},\widehat{\bf{R}}_{3})$. The last term on the right of Eq. $(\ref{eq3.3.1.5})$ is negative and plays a role in reducing numerical diffusion on entropy waves. One possible way to increase the entropy production is to increase the coefficient $\hat{{\delta}}_{2}$, thus less diffusion will be imposed on entropy waves. To demonstrate this, the modified equation associated with the HLLEM scheme is considered.

In reference [66], the modified equation for the HLL scheme is obtained through a simple approach. Here, we follow closely the work of Wang et al. [66] and extend the strategies therein to the HLLEM method. Furthermore, the characteristic form of the modified equation is also derived. It serves as an essential role in analyzing the artificial viscosity on the individual characteristic waves, particularly the entropy waves. For the steady shock problem in the current study, we focus only on the numerical dissipation introduced by the flux function. Therefore, only the semi-discrete equations are studied. The HLLEM scheme applied to the one-dimensional Euler system can be written in semi-discrete conservative form as

[TABLE]

where ${{\bf{F}}_{i+1/2}}$ is the numerical flux which is the approximation of the flux through the interface $x_{i+1/2}$. $\Delta x$ is the spatial cell size and $x_{i}=i\Delta x$, $i=1,...,N$. For the Godunov-type schemes in (29), (30) and (31), the system (32) can be recast into the equivalent viscosity form

[TABLE]

where ${\bf{Q}}_{i+1/2}$ is the numerical viscosity coefficient matrix. The modified equation of (33) has the form

[TABLE]

where $\widetilde{\bf{U}}$ is assumed to be a polynomial interpolated through the mesh values of $\bf{U}$ in the neighborhood of the mesh point $x_{i}$. The second order and even higher order error terms are suppressed in Eq. (34) in the context of modified equations.

To derive the modified equation associated with the HLLEM scheme, the numerical viscosity coefficient matrix needs to be determined, it can be written by sending ${\bf{U}}_{j}^{n}$ and ${\bf{U}}_{j+1}^{n}$ to ${\bf{U}}$,

[TABLE]

where $\bf{A}\left(\bf{U}\right)$ is the Jacobian of the flux function $\bf{F}\left(\bf{U}\right)$ and ${\bf{A}}\left({\bf{U}}\right)=\frac{{\partial{\bf{F}}\left({\bf{U}}\right)}}{{\partial{\bf{U}}}}$, $\bf{I}$ is the identity matrix, $S_{L}:=S_{L}\left({\bf{U}},{\bf{U}}\right)$ and $S_{R}:=S_{R}\left({\bf{U}},{\bf{U}}\right)$.

We want to analyze the artificial viscosity on the individual characteristic waves, especially the entropy wave. To this end, we derive the characteristic form of the modified equation. The right eigenvectors of the Jacobian ${\bf{A}}\left({\bf{U}}\right)$ are $\bf{R}_{1}$, $\bf{R}_{2}$, $\bf{R}_{3}$ and their corresponding eigenvalues are

[TABLE]

The transformation matrices $\bf{R}$ and its inverse $\bf{R}^{-1}$ are applied to the modified Eq. (34), by sending ${\widetilde{\bf{U}}}$ to ${\bf{U}}$,

[TABLE]

The Eq. (37) can be further simplified by using the relations ${\bf{A}}={\bf{R\Lambda}}{{\bf{R}}^{-1}}$ and ${\rm{d}}{\bf{W}}={{\bf{R}}^{-1}}{\rm{d}}{\bf{U}}$, the modified equation in characteristic variables can be obtained as

[TABLE]

where ${\bf{\Lambda}}={\rm{diag}}\left({u-a,u,u+a}\right)$ is the matrix of eigenvalues, the characteristic variables are defined by

[TABLE]

$\bf{V}$ is the viscosity matrix of the characteristic form that is expressed as

[TABLE]

The Eq. (38) can be written in a form where each individual characteristic wave is explicit presented

[TABLE]

where $\mu_{i}$ can be interpreted as numerical viscosity coefficients for the individual characteristic waves. To assess the numerical entropy change, we should only consider the coefficient of artificial viscosity corresponding to the entropy wave ${\rm{d}}w_{2}$,

[TABLE]

It can be observed that the last term on the right of Eq. (42) is nonpositive and plays a role in reducing numerical diffusion on the entropy wave. For the HLL flux, the coefficient $\delta_{2}$ equals zero. That is the reason why it contains excessive dissipation on the entropy wave. One possible way to increase the entropy production is to increase the coefficient ${\delta}_{2}$, thus less diffusion will be applied on the entropy wave. However, this will not work, because ${\hat{{\delta}}}_{2}$ used in Eq. (30) should be chosen such that a TVD-type condition for the viscosity matrix ${\bf{Q}}_{i+1/2}$ is valid [16],

[TABLE]

where ${\lambda}^{*}_{k}$ represents the eigenvalues of ${\bf{Q}}_{i+1/2}$ and ${\hat{\lambda}_{k}}$ denotes the eigenvalues of $\widehat{\bf{A}}$. It follows from (43) that the anti-diffusion coefficient ${\hat{\delta}}_{2}$ should satisfy the following condition

[TABLE]

Considering the definition of ${\widehat{{\delta}}}_{2}$ presented in Eq. (19), it is not desirable to increase the magnitude of the coefficient ${\hat{{\delta}}}_{2}$. Hence, we resort to modifying the dissipative terms in the wave strength ${\hat{\alpha}}_{2}$ to control the numerical dissipation on the entropy wave.

3.3.2 Two alternative modified schemes

In this section, two alternative modified schemes which contain less numerical dissipation on entropy waves are constructed. To avoid smearing the shock profile, the solutions ${\bf U}_{L}$ and ${\bf U}_{R}$ of the Riemann problem (26) for the supersonic case should not be altered. Here, two alternative approaches are considered. First, the average state ${\bf{U}}_{i+1/2}$ in (26) can be modified as

[TABLE]

with the following modified wave strength as

[TABLE]

To indicate this modification, the modified scheme is denoted as the HLLEM-${\rho}$ scheme. The modification of the average state ${\bf{U}}_{i+1/2}$ does not change the integral (6), therefore the Riemann solver (45) can still be written in the following conservation form

[TABLE]

with the modified wave strength $\widehat{\alpha}^{\ast}_{2}$ defined in $(\ref{eq3.3.2.2})$.

Alternatively, the pressure difference term in $\widehat{\alpha}_{2}$ can also be modified. The new scheme denoted as the HLLEM-$p$ can be expressed by

[TABLE]

with a modified wave strength defined as

[TABLE]

The coefficients $f_{\rho}$ and $f_{p}$ are used to determine numerical diffusion on entropy waves. In order to demonstrate how they will work, the dissipation mechanisms of these Godunov-type schemes are analyzed. To this end, we apply the strategies presented in section 3.3.1 to derive modified equations of these two new schemes. Here, the HLLEM-$p$ scheme is used to demonstrate the main procedure, the derivation of the modified equation of the HLLEM-$\rho$ scheme is described in Appendix A.

As demonstrated in the above section, the modified equation of the HLLEM-$p$ scheme can be derived if its numerical viscosity coefficient matrix is determined. The flux function defined in Eq. (48) and (49) can be rewritten by

[TABLE]

where the numerical viscosity matrix ${\bf{Q}}_{i+1/2}$ is defined by

[TABLE]

with the following relation

[TABLE]

The modified equation corresponding to the HLLEM-$p$ scheme can be written as

[TABLE]

where the viscosity matrix can be obtained,

[TABLE]

The modified equation in characteristic variables can be obtained following the procedure presented in Eq. (37), it can be written as

[TABLE]

where the following relation has been used,

[TABLE]

with

[TABLE]

To access the numerical viscosity corresponding to the entropy wave, we need to consider the following modified equation for the entropy wave ${\rm{d}}w_{2}$ that is extracted from Eq. (55),

[TABLE]

where the numerical viscosity coefficient $\mu_{2}$ is defined in Eq. (42). Considering (56) and (57), Eq. (58) can be splitted into the following subsystems,

[TABLE]

Due to the nonpositivity of the last term at the right side of Eq. (60), the dissipative term plays a role in smearing the wave ${\rm{d}}w_{2}^{2}$, thus a reduced $f_{p}$ will lead to a smaller ${\rm{d}}w_{2}^{2}$. Combining with the relations (56) and (57), an increased entropy will be obtained. Similarly, an increased $f_{\rho}$ can also lead to an increased entropy, which is indicated by the dissipation analysis of HLLEM-$\rho$ flux in Appendix A. However, one should be noted that the stability condition (43) for the contact discontinuity is no longer valid when $f_{\rho}$ is greater than one. Thus, the HLLEM-${\rho}$ flux can not be used for the multidimensional case where contact discontinuities may exist. Whereas, it works well for the mentioned steady shock problem in one dimension.

3.4 Experimental study on entropy production and shock instability

In the previous section, we have demonstrated that the numerical dissipation of HLLEM scheme on entropy waves can be controlled by modifying the wave strength in the dissipative term. Furthermore, it is indicated that the dissipation on entropy waves can be reduced by increasing the magnitude of density difference $\Delta{\rho}$ or reducing the magnitude of the pressure difference $\Delta p$ in the wave strength ${\widehat{\alpha}}_{2}$. If these modified schemes are applied at shocks, then a larger entropy production can be obtained inside the numerical shock structure. The coefficients $f_{\rho}$ and $f_{p}$ play a decisive role in controlling the entropy production. Since it is hard to analytically determine the amount of $f_{\rho}$ and $f_{p}$, we resort to conducting numerical experiments to determine their proper values.

Considering the one-dimensional steady shock problem in section 3.1, the computations are conducted with the first-order HLLEM scheme. A three-stage first-order Runge-Kutta method [67] is used with CFL=0.5 until the time step reaches 40,000. The Mach number is set as ${M_{0}}=6.0$ and the shock position is $\varepsilon=0.3$. The modified schemes HLLEM-${\rho}$ and HLLEM-$p$ are applied at the interfaces of cell M to control the entropy production in it. As demonstrated in section 3.1.2, the HLLEM scheme will produce the one-dimensional carbuncle under this condition. Thus, to obtain a stable and converged solution, we gradually increase the magnitude of the coefficient $f_{\rho}$ or reduce the magnitude the coefficient $f_{p}$ to increase the entropy production inside the intermediate cell M. In Fig. 5, entropy profiles of stable solutions computed with different coefficients $f_{\rho}$ and $f_{p}$ are illustrated. The corresponding residual histories are demonstrated in Fig. 6. As shown in Fig. 5, a larger entropy $S=\log\frac{p}{{\rho}^{\gamma}}-\log\frac{p_{L}}{{\rho}_{L}^{\gamma}}$ will be produced inside the numerical shock structure (i.e., the cell M) in the case that a larger $f_{\rho}$ or a smaller $f_{p}$ is used. Such results are in accordance with the dissipation analysis presented in section 3.3.2 and they show that if enough entropy production is guaranteed inside the shock, then the instability can be successfully eliminated. Residual histories presented in Fig. 6 indicate that more entropy production will make the solution converge to a stable state more quickly. It should be noted that similar results are also observed for various combinations of Mach numbers $M_{0}$ and the shock positions $\varepsilon$. In Fig. 7, we plot the stable thresholds of $f_{{\rho},threshold}$ and $f_{p,threshold}$ for different Mach numbers $M_{0}$ and their corresponding entropy thresholds are illustrated in Fig. 8.

3.5 Extension to the multidimensional case

So far, we have discussed the shock instability for the one-dimensional steady shock problem. In this part, the entropy-control technique developed in the above sections is extended to the two-dimensional case. As shown in Fig. 1, the numerical setup for the two-dimensional steady shock problem is the same as its one-dimensional counterpart except that 25 cells are used in the $j$ direction. This simple numerical test is commonly used to analyze numerical shock instability of Godunov-type schemes. It is believed that numerical methods that produce shock instabilities for this problem also produce carbuncle solutions for the blunt body problem [33, 3, 57].

As shown in Fig. 9, the conservative quantities inside the numerical shock structure may suffer errors when fluids pass through shock with perturbations. Due to discrepancies of physical quantities in the transverse of shock, new Riemann problems may exist at interfaces (see the red solid lines in Fig. 9). In the context of Godunov-type methods, one-dimensional planar waves will be generated at these interfaces. One should be noted that, physically, the conservative quantities in the transverse direction should be uniform and there should be no mass flux appearing in the $j$ direction due to the nature of quasi-one-dimensional flows. Thus, these waves may be not actually present. During the computation of this two-dimensional problem, dissipation would be added to stabilize these waves which may be much more than is actually required. Inside the numerical shock structure, it leads to insufficient entropy production that triggers the instability. To further validate this heuristic consideration, we also resort to numerical experiments.

The computations are conducted with the first-order HLLEM scheme for the two-dimensional steady shock problem. A three-stage first-order Runge-Kutta method [67] is used with CFL=0.5 until the time step reaches 40,000. To increase the entropy production, the HLLEM-$p$ fluxes are applied at all the interfaces of the intermediate cells inside the numerical shock structure. Stable and converged solutions can be obtained by reducing the magnitude of the coefficient $f_{p}$ properly. In Fig. 10, computed solutions by the original HLLEM scheme and the HLLEM scheme with the HLLEM-$p$ fluxes applied at interfaces of intermediate cells are demonstrated. It can be observed that the computed solution of the HLLEM scheme exhibits obvious shock instabilities, the shock front is severely distorted. With the HLLEM-$p$ flux applied at interfaces of intermediate cells inside the shock, the instability is barely observed and the shock is captured stably without any anomalies. The entropy profile extracted from the stable solutions are plotted in Fig. 11. It indicates that a smaller $f_{p}$ leads to larger entropy production inside the numerical shock structure, which is again well accordant with the dissipation analysis presented in section 3.3.2. Residual histories presented in Fig. 12 indicate that more entropy production will make the solution converge to a stable state more quickly. It should be noted that similar results are also observed for various combinations of Mach numbers $M_{0}$ and the shock positions $\varepsilon$. In Fig. 13, we also provide the stable thresholds of $f_{p,threshold}$ for different Mach numbers $M_{0}$ and their corresponding entropy thresholds are also included.

4 Numerical shock instability revisited: A linearized perturbation analysis

In the previous sections, we focus on analyzing the stability of approximate Riemann solvers to solve steady shocks. It has been demonstrated that the root for such numerical instabilities is due to the improper entropy production inside the shock structure. In this section, we revisit the instability of approximate Riemann solvers for the steady shock problem through a linearized perturbation analysis. This analysis method is first used by Quirk [7] and then followed by many other researchers [68, 4, 6, 43] to explore the mechanism of numerical shock instabilities. One advantage of the linearized perturbation analysis is that it is able to provide us an intuitive way to understand the mechanism of shock instability in the view of perturbations. However, one should always be caution that the linear analysis is supposed to only apply to cases where the perturbation errors are small enough. Thus, it can only provide a qualitative description of the instability.

4.1 One-dimensional case

In section 3, the improper entropy production inside the numerical shock structure is found to be the root cause of the one-dimensional carbuncle. By modifying the wave strength of dissipative term in Godunov-type schemes, enough entropy production can be guaranteed inside shocks. In order to further understand the relation between the entropy production and the instability, the instability problem is revisited in the context of the linearized perturbation analysis. In our previous work [6], it has been found that if the mass flux across the normal shock is accurately preserved (i.e., $(\rho u)^{n}_{L}=(\rho u)^{n}_{R}$), then the shock could be stabilized. The erroneous mass flux is originated from the intermediate states inside the shock structure. Thus, we need to clarify how the entropy-control technique could influence the mass flux perturbation behind the shock. Here, the linearized perturbation analysis is conducted again for the numerical flux functions HLLEM-$\rho$ and HLLEM-$p$. Readers are referred to reference [6] for more details.

At time $t^{n}$, the instability happens. It is assumed that there are some perturbation errors being generated in the cell M, they are expressed as

[TABLE]

where $\delta(\cdot)$ denote the perturbation errors which represent small discrepancies from the stable steady states. ${\left(\cdot\right)^{*}}$ denote the stable steady solutions that are assumed to locate at a Hugoniot curve. They have the following relationship with states in the cell L and the cell R,

[TABLE]

where the coefficients $\alpha_{\rho}$, $\alpha_{u}$ and $\alpha_{p}$ are defined in Eq. (25).

To explore how the perturbation errors generated in cell M influence the mass flux perturbation in cell R, we need to consider the following conservative scheme

[TABLE]

The subscript $R,R^{\prime}$ denotes the interface between the cell $R$ and the cell $R^{\prime}$ , the subscript $M,R$ denotes the interface between the cell $M$ and the cell $R$. For the HLLEM-$p$ scheme, the numerical fluxes at the interfaces in Eq. (63) can be written as

[TABLE]

Here, the subscript M denotes the cell M, the subscript R denotes the cell R that is at the right side of M. $\widehat{(\cdot)}$ are Roe’s averaged variables between states in cells M and R. Inserting (64) into (63), after some calculations we will obtain that

[TABLE]

with

[TABLE]

where the coefficients $\omega$ and $\sigma$ are defined as

[TABLE]

In Eq. (66), $\nu$ denotes the Courant number and the coefficients ${\theta}_{\rho}$, ${\theta}_{u}$ and ${\theta}_{p}$ are functions of the freestream Mach number $M_{0}$ and the shock position parameter $\varepsilon$. A similar result can also be obtained for the HLLEM-${\rho}$ flux. The recursive formulas of perturbed quantities can still be written as Eq. (65), but different coefficients are obtained as

[TABLE]

where the coefficients $\omega$ and $\sigma$ are still defined in Eqs. (67).

Due to the nonpositivity of the wave speed $S_{M}$, it can be observed from (66) that a smaller $f_{p}$ can lead to a reduced ${\theta}_{p}$, then a reduced mass flux $\delta(\rho u)_{R}^{n+1}$ will be obtained. Thus, the instability can be suppressed. Similarly, a larger $f_{\rho}$ leads to a reduced ${\theta}_{\rho}$, then a reduced mass flux $\delta(\rho u)_{R}^{n+1}$ will also be obtained, which plays a role in eliminating the instability. In conclusion, the entropy-control technique is equal to reducing the perturbation errors of mass flux which contribute to the instability in the view of perturbations analysis.

4.2 Two-dimensional case

The magnitude of the velocity in the transverse direction of the shock wave has been well recognized to be a proper parameter to use to show the magnitude of the multidimensional carbuncle phenomenon [33, 69]. For the steady normal shock problem considered in section 3, physically, there should be no mass flux appearing in the transverse direction and the transverse velocity should be zero. Hence, any erroneous mass flux developed in this direction results from the instability.

The dimensional splitting method [70] is used to update the solutions of the two-dimensional normal shock problem. For a particular cell $(i,j)$ inside the shock structure, the updated cell average ${\bf{U}}_{i,j}^{n+1/2}$ is obtained by solving the one-dimensional conservative scheme in the $x$-direction,

[TABLE]

where ${\bf{F}}_{i+1/2,j}^{n}$ is the numerical flux function for the one-dimensional Riemann problem between cells $(i,j)$ and $(i+1,j)$. Then the state ${\bf{U}}_{i,j}^{n+1/2}$ is used as data to update the solution in the $y$-direction,

[TABLE]

where ${\bf{G}}_{i,j+1/2}^{n+1/2}$ is the numerical flux function between cells $(i,j)$ and $(i,j+1)$. At the beginning of the instability, perturbations are generated inside the shock structure. As shown in Fig. 9, the perturbed states are distributed along the $j$ direction in a saw-toothed manner. To facilitate further analysis, it is assumed that the states along the $y$ direction inside the shock structure are initialized as follows,

[TABLE]

and

[TABLE]

where ${()}^{*}$ represent the stable steady solutions that are assumed to be uniform along the transverse direction. In what follows, we omit the subscript $i$ for clarity. In the two-dimensional case, we need to clarify how the perturbations will promote the perturbed mass flux in the transverse direction. Hence, the following conservative scheme is considered,

[TABLE]

For the flux function HLLEM-$p$, the numerical fluxes at the interfaces in (73) can be written as

[TABLE]

where ${\widehat{\left(\cdot\right)}_{j+1/2}}$ denote Roe’s averaged variables between states in cells $(i,j)$ and $(i,j+1)$. Inserting (74) into (73), we can obtain the following relationship (see Appendix B for the detailed derivation)

[TABLE]

with

[TABLE]

It can be observed from (76) that a smaller $f_{p}$ can lead to a reduced ${\theta}_{p}$, then a reduced mass flux $\delta(\rho v)_{j}^{n+1}$ will be obtained. Thus, the instability can be suppressed. Similar to its one-dimensional case, the entropy-control technique is equal to reducing the mass flux perturbation errors which contribute to the instability in the view of perturbations analysis.

5 A general technique for suppressing the shock instability by an entropy control

In above sections, efforts have been devoted to investigating the underlying mechanism of numerical shock instability for the HLL-type schemes. It has been demonstrated that the shock instability can be attributed to improper entropy production inside the numerical shock structure. In this section, we extend the entropy-control technique discussed in section 3.3.2 to other Godunov-type schemes to demonstrate its generality.

The proposed numerical flux function with an Entropy Control can be expressed as

[TABLE]

where $``\rm{X}"$ denotes certain Godunov-type schemes. The entropy-control term ${\bf{F}}_{\rm{EC}}$ can be derived directly from the mentioned HLLEM-$p$ flux . It is defined by

[TABLE]

where ${\widehat{\delta}}_{2}$ represents the anti-diffusion coefficient defined in Eq. (19), ${\widehat{{\bf{R}}}}_{2}$ denotes the second right eigenvector defined in Eq. (17), $\widehat{a}$ is the sound speed evaluated by Roe’s averaged variables, $\Delta p$ is the pressure difference $\Delta p=p_{R}-p_{L}$. The wave speeds $S^{-}_{L}:=\min(0,S_{L})$ and $S^{+}_{R}:=\max(0,S_{R})$ are evaluated by (10), which provides minimal diffusion on isolated shocks. Alternatively, the estimate (11) can also be used to compute wave speeds $S^{-}_{L}$ and $S^{+}_{R}$, they are more diffusive and useful for suppressing possible post-shock oscillations. To complete the evaluation of the resulting flux function ${\bf{F}}_{\rm{X-EC}}$, $f_{p}$ needs to be determined, it is written by

[TABLE]

with

[TABLE]

In the vicinity of shock waves, the pressure difference is increased and the function $f_{p}$ becomes small, the entropy control term ${\bf{F}}_{\rm{EC}}$ is activated to increase the entropy production inside the numerical shock structure, thus enhancing the stability. In the smooth region, the function $f_{p}$ approaches unit, the modified flux turns into its original form. It should be noted that the entropy-control term ${\bf{F}}_{\rm{EC}}$ does not involve any numerical dissipation on linear degenerate waves, thus it avoids the risk of introducing unnecessary numerical dissipation to contaminate flow structures. In the current study, the entropy-control technique is applied to four classical Godunov-type schemes, i.e., “X” denote Godunov’s exact Riemann solver [14], Roe [52], HLLEM [16] and HLLC [48] schemes, which are known to resolve discontinuities with minimal diffusion but are plagued by numerical shock instabilities.

6 Numerical experiments

In what follows, we apply our approach to a series of numerical experiments already used for the validation of several improved Godunov-type schemes[6, 51]. The former four problems are selective test cases to access the robustness of the proposed method against strong shock waves. A focus of other problems is on accuracy for resolving shear layers.

6.1 Quirk’s odd-even grid perturbation problem

Quirk’s odd-even decoupling problem [7] is a well-known test used to access the resistance of numerical schemes against shock instabilities. A plane normal shock wave moving from the left to the right at Mach 6 is computed on a regular grid system with a slight perturbation in one of the grid lines parallel to the direction of the flow. The computational domain is covered by an $800\times 20$ structured grid involving the following grid perturbations:

[TABLE]

where $y(i,j_{mid})$ is the $y$ coordinate of a vertex $(i,11)$, $y_{mid}$ is the $y$ coordinate of the halfway line. The initial conditions are $\left({\rho,u,v,p}\right)=\left({1.4,0,0,1}\right)$, the inlet boundary condition is set to post shock values, nonreflecting simple wave boundary conditions are imposed at the outlet. The upper and lower boundaries are considered as solid walls. This problem is solved by using different flux functions, along with the two-stage second-order Runge-Kutta method [71]. The CFL number is set to 0.5 for all the schemes to compute the solutions. Fig. 14 shows the contour plots of density at the time levels $t=0$, $t=25$ and $t=50$, where there are $20$ contour levels varying from $1.5$ to $6.0$.

As shown in Fig. 14, the shocks were destroyed by all the original low diffusion approximate Riemann solvers. Compared with the Roe and HLLEM solvers, the Godunov’s exact Riemann solver and HLLC scheme are less prone to the irregularities. Whereas, typical carbuncles were developed and could be observed clearly. With the entropy-control technique, all the schemes produce no odd-even decoupling along the shock. The shock waves were kept all the way through.

6.2 Shock diffraction

The problem consists of the diffraction of a plane normal shock wave moving at the Mach number 5.09 around a $90^{\circ}$ corner. The computational domain is a unit square $\left[{0,1}\right]\times\left[{0,1}\right]$ that is discretized into $400\times 400$ uniform cells. The corner is located at the midpoint of the left boundary: the lower half is treated as a wall; the top half is taken as an inflow. The shock is initially located at the midpoint of the left boundary. To the right of the shock, the domain is initialized to pre shock conditions of $\rho=1.4$, $p=1$, $u=0$ and $v=0$. The domain to the left of the shock is initialized to post shock conditions. The inflow boundary condition is imposed at the left boundary and the slip boundary condition is imposed at the corner and the top boundary. The right and bottom boundaries are taken as outflow. The simulations are conducted by 2nd order accurate schemes with the CFL value of 0.5. The plots shown in Fig. 15 have been generated at time $t=0.15$. A total of 20 density contour levels varying from 0.5 to 5.0 have been illustrated. As shown, solutions computed by different original flux functions all exhibit obvious post-shock oscillations. Whereas, with the proposed entropy-control technique, these post-shock oscillations are completely eliminated and the instabilities are not observed. We remark that the first-order accurate numerical results show the similar phenomena.

6.3 Double Mach reflection problem

It is well-known that many low diffusion Riemann solvers produce kinked Mach stems in the double Mach reflection problem. This test case was first studied by Woodward and Colella [72], and later by other scholars to test the performance of numerical schemes. The computational domain is $\left[{0,\;4}\right]\times\left[{0,\;1}\right]$ which has been divided into 960 cells along the length and 240 cells along the width. The shock with a Mach number of 10 is initially set up to be inclined at an angle of 60 with the bottom reflecting wall. The domain in front of the shock is initialized with pre shock values given as $\rho=1.4$, $u=0$, $v=0$, $p=1$. The domain behind the shock is initialized to post shock values. The computations are performed by first-order numerical schemes and the third-order TVD Runge-Kutta time discretization [71] with $\rm{CFL}=0.5$ up to $t=0.2$. The density contours computed by different schemes are shown in Fig. 16, where 20 contour levels varying from 2.0 to 20.0 are used.

As shown, all the original methods produce the unphysical triple points, i.e., the kinked Mach stems. With the proposed entropy-control technique, these modified schemes are all able to resolve shocks without any irregularities and the kinked Mach stems are barely noticeable.

6.4 Hypersonic flow over a cylinder

The simulation of hypersonic inviscid flow over a cylinder is a typical problem to assess the performance of a scheme with respect to the carbuncle. This test problem is computed on a quadrilateral grid with 60 cells in the radial direction, 280 cells in the circumferential direction. The freestream Mach number is 20. The computational domain has been initialized with values $\rho=1.4$, $p=1.4$, $u=20$ and $v=0$. The simulations are conducted by using different flux functions, along with the LU-SGS approach [73]. Computations are conducted for 50,000 time steps with $\rm{CFL}=0.5$.

In Fig. 17, the density contours by different approximate Riemann solvers are demonstrated, where 20 contour levels varying from 1.5 to 8.5 are used. As shown in Fig. 17, all the unmodified numerical schemes produce the carbuncle phenomena. The pressure profiles along the stagnation line for different schemes have been plotted in Fig. 18. As shown, the Roe scheme does not produce a correct solution due to the carbuncle. Whereas, the Roe-EC solver is able to produce a correct pressure profile with only one interior point. For other approximate Riemann solvers, similar results are obtained but not shown for clarity.

6.5 Flat plate boundary layer

A laminar boundary layer problem is used here to assess the ability of the proposed scheme to resolve shear layers. The Mach number of the freestream is $Ma=0.3$ and the Reynolds number is ${\mathop{\rm Re}\nolimits}={10^{5}}$. The problem is computed by a second-order Navier-Stokes code. In Fig. 19, the computational domain is discretized into $120\times 30$ non-uniform cells. The non-slip adiabatic boundary condition is imposed at the plate and the bottom boundary before the flat plate is treated as a symmetry plane. The non-reflecting boundary condition based on the Riemann invariants is adopted for the other boundaries. The computations were conducted for 50,000 steps with $\rm{CFL}=0.5$ , and all the computations achieved at least three orders of magnitude reductions of the density residuals. In Fig. 20, velocity profiles computed by different schemes are compared. As shown, the Roe scheme with an entropy-control technique gives an almost identical result with that of the original Roe flux. The solutions match with the exact Blasius solutions very well. The entropy-control term does not decrease the accuracy of numerical schemes for resolving shear layers. Other flux functions give the same results, but not shown for clarity. Results computed by the HLLE scheme are also included for comparison. It can be observed that it gives a dissipative and inaccurate solution.

6.6 Elling’s physical carbuncle problem

The last test case considered is a challenging problem. It is first proposed by Elling [27] and later studied by other scholars [74] to test numerical schemes to resolve physical carbuncle problem. In references [26, 27], the author proposes that carbuncles may be a special class of entropy solutions which can be physically correct in some circumstances. The description of this problem is presented in Fig. 21. It is set up to model the interaction of a vortex filament with a strong shock. Numerical schemes that are able to avoid unphysical carbuncles must also allow this physical carbuncle-like structure. Thus, it is a suitable test case to access whether the improvement for carbuncles will reduce the accuracy of numerical methods on resolving shear layers.

Following Kemm [74], this case is tested in the computational domain $[0,100]\times[0,4]$ which has been divided into 100 cells along the length and 40 cells along the width. The shock of a Mach number $M_{0}=20$ is located at a cell interface. The domain pre to the shock is initialized with $\rho=1$, $u=1$ , $v=0$, $p=\frac{1}{\gamma M_{0}^{2}}$. The post shock condition is calculated by the Rankine-Hugoniot relations. In the supersonic inflow region, only the middle slice of the computational grid is changed and the velocity is artificially set to zero. Computations are performed by second-order accurate numerical schemes using the third-order TVD Runge-Kutta time discretization [71] with $\rm{CFL}=0.2$ up to $t=100$. In Fig. 22, the entropy contours $s=\log(p/{\rho}^{\gamma})$ computed by different approximate Riemann solvers are demonstrated, where 15 contour levels varying from -6.0 to -2.5 are used. As shown, the numerical results computed by Godunov’s exact Riemann solver, Roe, HLLEM and HLLC are rather similar. With their modified versions, the physical carbuncles are properly reproduced. This indicates that the proposed entropy-control technique involves minimal diffusion on shear layers. Whereas, the HLLE scheme smears the physical carbuncle severely, demonstrating its excessive amount of shear viscosity. The last picture presents the result computed by the hybrid solver Roe-HLLE suggested in [7]. This type of numerical methods usually apply a switching function to combine a complete Riemann solver (e.g. Roe scheme, HLLC scheme, etc.) with an incomplete Riemann solver (e.g. HLLE scheme, van Leer scheme, etc.) to eliminate shock instabilities while maintaining high resolution in smooth regions. As shown, the physical carbuncle is also destroyed. This indicates that hybrid approximate Riemann solvers may lose their high resolution in the vicinity of shock-vortex interactions.

7 Conclusion

In the current study, we devote our efforts to developing a general approach to improve the robustness of Godunov-type schemes at strong shocks without compromising their accuracy on discontinuities. To this end, a shock instability analysis has been carried out by combining two interlinked techniques: dissipation analysis of Godunov-type schemes and numerical experiments. Results of the stability analysis demonstrate that the numerical shock instability of Godunov-type schemes is due to inappropriate entropy production inside the numerical shock structure. Moreover, if enough entropy production is guaranteed inside shocks, then the instability problem can be successfully eliminated. Such results allow to propose a general approach that can be applied to suppressing numerical shock instabilities of various Godunov-type schemes. The essential component of this approach is an entropy-control term that plays a role in increasing the entropy production in the numerical shock structure. A further linearized analysis illustrates that this entropy-control term is equivalent to reducing mass flux perturbations behind shocks, which is found to trigger the shock instability in the view of perturbations. Several stringent shock wave problems are conducted to indicate that the proposed entropy-control approach can effectively improve the robustness of low diffusion Godunov-type schemes at strong shocks. Meanwhile, comparisons of computed solutions from a boundary layer problem and a shock-vortex interaction problem indicate that these modified methods preserve the accuracy of the original solvers on resolving shear layers.

The proposed entropy-control technique is expected to also applied to approximations to other hyperbolic systems of conservation laws, such as the shallow water equations and the MHD equations. Extensions to cases of high-order schemes will also be considered in further investigations.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (11472004), the Foundation of Innovation of NUDT (B150106).

Appendix A

Here the modified equation for the HLLEM-$\rho$ scheme is presented. As illustrated in section 3.3.1, the modified equation of the HLLEM-${\rho}$ scheme can be derived if its numerical viscosity coefficient matrix is determined. The flux function of the HLLEM-${\rho}$ can be expressed as

[TABLE]

where the numerical viscosity coefficient matrix ${\bf{Q}}_{i+1/2}$ is defined by

[TABLE]

with the following relation

[TABLE]

In Eq. (84), $\widehat{\alpha}_{2}^{*}$ denotes the modified wave strength corresponding to entropy waves and is defined in Eq. (46), $\widehat{\bf{R}}^{*,-1}$ is determined by Eq. (84).

The modified equation corresponding to the HLLEM-$\rho$ scheme can be written as

[TABLE]

where the viscosity matrix can be obtained as

[TABLE]

The modified equation in characteristic variables can be obtained following the procedure presented in Eq. (37), it can be written as

[TABLE]

where the viscosity coefficient matrix $\bf{V}$ is defined in (40), the characteristic variables $\rm{d}{\bf{W}}^{1}$ and $\rm{d}{\bf{W}}^{2}$ are defined in Eq. (56) and Eq. (57).

To access the numerical viscosity corresponding to the entropy wave, we need to consider the following modified equation for the entropy wave ${\rm{d}}w_{2}$,

[TABLE]

where the numerical viscosity coefficient $\mu_{2}$ is defined in Eq. (42). Following Eq. (56) and (57), the entropy wave ${\rm{d}}w_{2}$ can be splitted into the following form,

[TABLE]

with

[TABLE]

Considering (89) and (90), Eq. (88) can be splitted into the following subsystems,

[TABLE]

Due to the nonpositivity of the last term at the right side of Eq. (91), a larger $f_{\rho}$ will lead to less diffusion on ${\rm{d}}w_{2}^{1}$. Combing with the relation (89), an increased entropy will be obtained.

Appendix B

Here we give a detailed description of how to obtain the formulas presented in Eq. (75) and Eq. (76). At time $t^{n}$, it is assumed that the states along the $y$ direction inside the shock structure are distributed in a saw-toothed manner. These solutions are given by

[TABLE]

and

[TABLE]

where ${()}^{*}$ represent the stable steady solutions. In what follows, we omit the subscript $i$ for clarity. In the two-dimensional case, we need to clarify how the perturbations will promote the perturbed mass flux in the transverse direction. Hence, the following momentum equation is considered,

[TABLE]

For the flux function HLLEM-${\rho}$, the numerical fluxes at the interfaces in (95) can be written as

[TABLE]

where ${\widehat{\left(\cdot\right)}_{j+1/2}}$ denote Roe’s averaged variables between states in cells $(i,j)$ and $(i,j+1)$. The Roe’s averaged varibles used can be approximated as

[TABLE]

and

[TABLE]

where the second order and even higher order perturbation terms are suppressed in the context of linearized analysis. Similarly, the Roe’s averaged sound speed can be written by

[TABLE]

Inserting Eqs. (93), (94) and (97 99) into Eqs. (96), one can obtain,

[TABLE]

Inserting Eqs. (100) into Eq. (95), it can be obtained that

[TABLE]

It should be noted that

[TABLE]

Subtracting Eq. (102) from Eq. (101), it can be obtained that

[TABLE]

where the term $\frac{\Delta t}{\Delta y}$ can be approximated as

[TABLE]

In Eq. (104), $\nu$ denotes the Courant number.

Inserting (104) into (103), we can obtain the following relationship

[TABLE]

with

[TABLE]