A simple cure for numerical shock instability in HLLC Riemann solver

TL;DR

This paper introduces the HLLC-SWM framework, a modification to the HLLC Riemann solver that enhances shock stability by selectively increasing dissipation near shocks, while preserving accuracy on other wavefields.

Contribution

The paper proposes a novel framework, HLLC-SWM, that improves shock stability of the HLLC scheme by adding targeted dissipation without sacrificing accuracy on contact and shear waves.

Findings

HLLC-SWM schemes are shock stable over a wide Mach number range.

They retain positivity and exact contact resolution.

HLLC-SWM-P is highly accurate for viscous flows.

Abstract

The Harten-Lax-van Leer with contact (HLLC) scheme is known to be plagued by various forms of numerical shock instabilities. In this paper, we propose a new framework for developing shock stable, contact and shear preserving approximate Riemann solvers based on the HLLC scheme for the Euler system of equations. The proposed framework termed as HLLC-SWM (\textbf{S}elective \textbf{W}ave \textbf{M}odified) scheme identifies and increases the magnitude of the inherent diffusive HLL component within the HLLC scheme in the vicinity of a shock wave while leaving its antidiffusive component unmodified to retain accuracy on linearly degenerate wavefields. We present two strategies to compute the requisite supplementary dissipation which results in HLLC-SWM-E and HLLC-SWM-P variants. Through a linear perturbation analysis of the HLLC-SWM framework, we clarify how the additional dissipation introduced in this way helps in damping of unphysical perturbations in primitive quantities under a derived CFL constraint. A matrix based stability analysis of a steady two-dimensional normal shock is used to show that both variants of the HLLC-SWM scheme are shock stable over a wide range of inlet Mach numbers. Results from standard test cases demonstrate that the HLLC-SWM schemes are capable of computing shock stable solutions on a variety of problems while retaining positivity and exact inviscid contact ability. On viscous flows, while the HLLC-SWM-P variant is quite accurate, the HLLC-SWM-E variant introduces slight inaccuracy which can be corrected through a simple Mach number based switching function.

Figures

Equations

$$\frac{\partial \acute{\mathbf{U}}}{\partial t}+\frac{\partial \acute{\mathbf{F}}(\mathbf{U})}{\partial x}+\frac{\partial \acute{\mathbf{G}}(\mathbf{U})}{\partial y}=0$$

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

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

$$\frac{\partial }{\partial t}{\int }_{\Omega }\acute{\mathbf{U}}dxdy+{\oint }_{d\Omega }[(\acute{\mathbf{F}},\acute{\mathbf{G}}).\mathbf{n}]dl=0$$

$$\frac{d{\mathbf{U}}_i}{dt}=-\frac{1}{|\Omega {|}_i}\sum _{k=1}^4[(\acute{\mathbf{F}},\acute{\mathbf{G}}{)}_k.{\mathbf{n}}_k]\Delta s_k$$

$$\frac{d{\mathbf{U}}_i}{dt}=-\frac{1}{|\Omega {|}_i}\sum _{k,m=1}^4[{\mathbf{T}}_k^{-1}\mathbf{F}({\mathbf{U}}_{\mathbf{L}},{\mathbf{U}}_{\mathbf{R}})]\Delta s_k$$

$${\mathbf{U}}_{\ast }^{HLL}$$

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

$$S_L=min(0,u_{nL}-a_L,{\tilde{u}}_n-\tilde{a})$$

$$S_R=max(0,u_{nR}+a_R,{\tilde{u}}_n+\tilde{a})$$

$${\mathbf{F}}_{HLLC}=\{\begin{array}{ll}{\mathbf{F}}_L+S_L({\mathbf{U}}_{\ast L}^{HLLC}-{\mathbf{U}}_L),& \text{if }{S}_L\le 0\le S_M\\ {\mathbf{F}}_R+S_R({\mathbf{U}}_{\ast R}^{HLLC}-{\mathbf{U}}_R),& \text{if }{S}_M\le 0\le S_R\end{array}$$

$${\mathbf{U}}_{\ast L/R}^{HLLC}$$

$$S_M=\frac{p_R-p_L+{\rho }_L{u}_L(S_L-u_L)-{\rho }_R{u}_R(S_R-u_R)}{{\rho }_L(S_L-u_L)-{\rho }_R(S_R-u_R)}$$

$${\mathbf{F}}_{HLLC}=\{\begin{array}{ll}{\mathbf{F}}_{HLL}+S_L({\mathbf{U}}_{\ast L}^{HLLC}-{\mathbf{U}}_{\ast }^{HLL}),& \text{if }{S}_L\le 0\le S_M\\ {\mathbf{F}}_{HLL}+S_R({\mathbf{U}}_{\ast R}^{HLLC}-{\mathbf{U}}_{\ast }^{HLL}),& \text{if }{S}_M\le 0\le S_R\end{array}$$

$${\mathbf{F}}_{HLL}=\frac{1}{2}({\mathbf{F}}_L+{\mathbf{F}}_R)+{\mathbf{D}}_{HLL}$$

$${\mathbf{D}}_{HLL}=\frac{S_R+S_L}{2(S_R-S_L)}({\mathbf{F}}_L-{\mathbf{F}}_R)-\frac{S_L{S}_R}{(S_R-S_L)}({\mathbf{U}}_L-{\mathbf{U}}_R)$$

$${\mathbf{D}}_{HLL}=\frac{|S_R|-|S_L|}{2(S_R-S_L)}({\mathbf{F}}_L-{\mathbf{F}}_R)+\frac{|S_L|S_R-|S_R|S_L}{2(S_R-S_L)}({\mathbf{U}}_L-{\mathbf{U}}_R)$$

$${\overline{\mathbf{D}}}_{HLL}=\frac{|\overline{S_R}|-|\overline{S_L}|}{2(S_R-S_L)}({\mathbf{F}}_L-{\mathbf{F}}_R)+\frac{|\overline{S_L}|S_R-|\overline{S_R}|S_L}{2(S_R-S_L)}({\mathbf{U}}_L-{\mathbf{U}}_R)$$

$$\overline{S_L}=S_L-\alpha {ϵ}^L$$

$$\overline{S_R}=S_R+\alpha {ϵ}^R$$

$${\mathbf{F}}_{HLLC}=\underset{HLLFlux}{\underbrace{{\mathbf{F}}_L+S_L({\mathbf{U}}_{HLL}^{\ast }-{\mathbf{U}}_L)}}+\underset{HLLCAntidissipation}{\underbrace{S_L({\mathbf{U}}_{HLLC}^{\ast }-{\mathbf{U}}_{HLL}^{\ast })}}$$

$${\mathbf{D}}_{HLL}=\frac{|S_R|-|S_L|}{2(S_R-S_L)}({\mathbf{F}}_L-{\mathbf{F}}_R)+\frac{|S_L|S_R-|S_R|S_L}{2(S_R-S_L)}({\mathbf{U}}_L-{\mathbf{U}}_R)$$

$$a_0$$

$$a_1$$

$$\overline{a_0}$$

$$\overline{a_1}$$

$$\overline{a_0}$$

$$\overline{a_1}$$

$$\{\begin{array}{ll}{\rho }_j={\rho }_0+\hat{\rho },p_j=p_0+\hat{p},u_j=u_0+\hat{u},v_j=0,& \text{if }jiseven\\ {\rho }_j={\rho }_0-\hat{\rho },p_j=p_0-\hat{p},u_j=u_0-\hat{u},v_j=0,& \text{if }jisodd\end{array}$$

$$\displaystyle\left[\begin{array}[]{c}\hat{\rho}\\ \hat{u}\\ \hat{p}\end{array}\right]^{n+1}=\left[\begin{array}[]{r}1\\ 0\\ 0\end{array}\begin{array}[]{r}0\\ 1\\ 0\end{array}\begin{array}[]{r}\frac{-2\lambda}{\gamma}\\ 0\\ 1-2\lambda\end{array}\right]\left[\begin{array}[]{c}\hat{\rho}\\ \hat{u}\\ \hat{p}\end{array}\right]^{n}$$

$$\displaystyle\left[\begin{array}[]{c}\hat{\rho}\\ \hat{u}\\ \hat{p}\end{array}\right]^{n+1}=\left[\begin{array}[]{r}1-\frac{2\lambda\alpha\epsilon}{\sqrt{\gamma}}\\ 0\\ 0\end{array}\begin{array}[]{r}0\\ 1-\frac{2\lambda\alpha\epsilon}{\sqrt{\gamma}}\\ 0\end{array}\begin{array}[]{r}-\frac{2\lambda}{\gamma}\\ 0\\ 1-2\lambda(1+\frac{\alpha\epsilon}{\sqrt{\gamma}})\end{array}\right]\left[\begin{array}[]{c}\hat{\rho}\\ \hat{u}\\ \hat{p}\end{array}\right]^{n}$$

$$0<\lambda <1$$

$$0<\lambda <\frac{1}{1+\frac{\alpha {ϵ}^{L,R}}{\sqrt{\gamma }}}$$

Full text

††journal: Journal of Computational Physics

1 Introduction

Accurate and robust computation of high-speed gas dynamical flows continues to be a challenge due to the presence of strong shocks, expansion fans, contact surfaces and their mutual interactions. The approximate Riemann solvers are one of the most popular methods available for numerical computations of such flows. However, it is well known that the contact and shear preserving variants of these schemes like the Roe [1], HLLEM [2] and HLLC [3] are susceptible to numerical shock instability during multidimensional computations while the linear wavefield dissipative schemes like the HLL [4] and the HLLE [5] are free from them. These essentially multidimensional phenomena are primarily observed during the simulation of strong grid aligned normal shocks. Benign perturbations in the initial conditions are rapidly amplified in the numerical shock region leading to spurious solutions. Most commonly, these instabilities manifest as oscillations in the shock profile, polluted aftershock values, growth in error norms in case of steady state problems and in extreme cases complete breakdown of the solution. Additionally, the ocurrences of these instabilities are also highly sensitive to parameters like grid aspect ratio [6], inflow Mach number[7], numerical shock structure [7, 8], perturbations in upstream supersonic flow [7] and order of accuracy of the solution [9] among other things. A detailed catalouge of these failures can be found in [7, 10, 11].

Quirk [10] conjectured that instability is caused because of coupling between pressure perturbations and density perturbations along the length of the shock which affects the numerical shock speed. Dumbser et al [7] developed a useful matrix based stability analysis that can predict shock instability for a given scheme under various grid configurations, inlet Mach numbers and boundary conditions. Kun Xu [12] and others [7, 8, 13] located the source of instability to be within the numerical shock region. Some researchers [13, 14, 15, 16] identify that inadequate dissipation of shear wave could result in these instabilities. Sanders et al [17] analyzed the numerical truncation error of upwind schemes to establish that shock instabilities in these schemes arise due to lack of cross-flow dissipation provided by these schemes in the vicinity of a shock wave. Shen et al [18] showed how the Roe scheme suffers an overall loss of numerical dissipation along a numerically perturbed shock front due to the vanishing of terms associated with the contact and shear waves. Further, they demonstrated that the addition of supplementary dissipation on interfaces normal to the shock front mitigated this problem and produced shock stable solutions while ensuring consistency in mass flux across the shock.

Based on these observations, various methods have been proposed in the literature to cure approximate Riemann solvers of this deficiency. All of these cures are based on introducing a certain amount of numerical dissipation that successfully attenuates the instabilities while inadvertently smearing the contact and shear waves. Thus there are hybrid Riemann solvers like Roe-HLLE [19, 20, 21] and HLLC-HLL [15, 22, 23, 24] and rotated Riemann solvers [16, 21, 25] that treat these instabilities with some success. While the hybrid methods require explicit switching sensors and careful choice of complementary Riemann solvers, the rotated methods are computationally at least twice as expensive to implement as their grid aligned counterparts. As yet, one of the simplest methods to deal with shock instabilities in Roe scheme involves increasing its overall numerical dissipation by modifying the eigenvalues of its linearized Jacobian matrix; including the ones corresponding to the linear waves. Peery and Imlay [26] demonstrated a cure for ’Carbuncle phenomenon’ in Roe scheme based on this idea. The modification they proposed used a direction based pressure gradient to tune the magnitude of the dissipation introduced. Lin [27] pointed out the limitations of this pressure sensor and proposed an improved version that works even on viscous flows. Sanders et al [17] suggested using the idea of an entropy fix on linear waves in conjunction with a simple ‘H’shaped multidimensional stencil to estimate the requisite dissipation to cure the Roe scheme. Later Pandolfi et al [11] pointed out a more efficient way to implement the same that avoids unwanted dissipation at shear layers for more anisotropic grids. However, as noted by some authors [7, 9], there is no physical justification for modifying the linear eigenvalues using an entropy fix.

Incidentally, a majority of cures for shock instability have been developed for the Roe scheme. However, Roe scheme suffers several drawbacks like requirement of an additional entropy fix on its nonlinear wave speeds to admit entropy satisfying solutions [28], knowledge of the full eigenstructure of the flux Jacobians making it difficult to be extended to general system of governing equations with complex equations of state [29] and lack of positivity [2]. These shortcomings of Roe scheme has brought forth an increasing interest in another equally accurate and efficient contact and shear preserving Riemann solver - the HLLC Riemann solver [3] which also, unfortunately, succumbs to various numerical shock instabilities. Hence the need for delivering the HLLC solver from these instabilities is paramount.

In this paper, we present a very simple and inexpensive strategy to develop a shock stable HLLC scheme. The strategy is termed HLLC-SWM (Selective Wave Modified). The novelty of the formulation lies in first decomposing the HLLC scheme into its diffusive HLL component and the antidiffusive component that restores the accuracy on linear waves. Once distinguished, we propose a strategy to adaptively increase the numerical dissipation associated with this inherent HLL component alone in the vicinity of shocks through appropriate modification of selected nonlinear wave speeds. In contrast to the cures proposed in [11, 17, 27], the present method only targets nonlinear wave speeds that resides in the diffusive component while leaving the middle wave speed which is crucial to the contact and shear preserving ability of the HLLC solver completely undisturbed. Various theoretical analyses and an array of numerical test problems confirm that this simple technique is found to provide sufficient dissipation in the vicinity of the shock to eliminate shock instabilities.

The outline of the paper is as follows. In Sections 2 and LABEL:sec:_FVM_framework a brief outline of the governing equation and the Finite Volume framework used in the present work is discussed. Section 4 reviews the HLL and the HLLC Riemann solvers. In Section 5 we discuss the formulation of the HLLC-SWM framework and prove the effectiveness of the proposed strategy in increasing the overall dissipation of the HLLC scheme. In Section 6, a linear analysis technique is used to investigate how the HLLC-SWM framework is effective in suppressing perturbations in primitive quantities. In Sections 7 and 8, we describe two methods to determine the amount of dissipation required by the HLLC-SWM formulation to suppress shock instabilities over a wide range of Mach numbers. In Section 9 the HLLC-SWM formulation is extensively tested on a series of numerical test problems. Some concluding remarks are given in Section 10.

2 Governing equations

The governing equations for two dimensional inviscid compressible flow can be expressed in their conservative form as,

[TABLE]

where $\mathbf{\acute{U}},\mathbf{\acute{F}(U)},\mathbf{\acute{G}(U)}$ are the vector of conserved variables and x and y directional fluxes respectively. These are given by,

[TABLE]

In the above expressions, $\rho,u,v,p$ and $E$ stands for global density, x-velocity, y-velocity, pressure and specific total energy. The system of equations are closed through the equation of state,

[TABLE]

where $\gamma$ is the ratio of specific heats. Present work assumes a calorifically perfect gas with $\gamma=1.4$. A particularly useful form of the Eq.(1) is the integral form given by,

[TABLE]

where $\varOmega$ denotes a control volume over which Eq.(16) describes a Finite Volume balance of the conserved quantities, $dx$ and $dy$ denotes the x and y dimensions of the control volume respectively, $d\varOmega$ denotes the boundary surface of the control volume and $dl$ denotes an infinitesimally small element on $d\varOmega$. $\mathbf{n}$ is the outward pointing unit normal vector to the surface $d\varOmega$.

3 Finite volume discretization

In this paper we seek a Finite Volume numerical solution of Eq.(16) by discretizing the equation on a computational mesh consisting of structured quadrilateral cells as shown in Fig.(1). For a typical cell $i$ belonging to this mesh, a semi-discretized version of such a solution can be written as,

[TABLE]

where ${\mathbf{U}_{i}}$ is an appropriate cell averaged conserved state vector, $({\mathbf{\acute{F},\acute{G}}})_{k}$ denotes the flux vector at the mid point of each interface $I_{k}$ while $\mathbf{n}_{k}$ and $\Delta s_{k}$ denotes the unit normal vector and the length of each $I_{k}$ interface respectively. These are shown in Fig.(1). The interface flux $({\mathbf{\acute{F},\acute{G}}})_{k}.\mathbf{n}_{k}$ can be obtained by various methods. One of the most popular class of methods to compute this are the approximate Riemann solvers. A conventional two state approximate Riemann solver uses the rotational invariance property of Euler equations to express the term $({\mathbf{\acute{F},\acute{G}}})_{k}.\mathbf{n}_{k}$ as,

[TABLE]

where $\mathbf{U_{L}}=\mathbf{T}_{k}(\mathbf{U}_{i}),\mathbf{U_{R}}=\mathbf{T}_{k}(\mathbf{U}_{m})$ indicates the initial conditions of a local Riemann problem across $k^{th}$ interface formed by the $i^{th}$ cell and its neighbouring cell $m$. The matrices $\mathbf{T}_{k}$ and $\mathbf{T}_{k}^{-1}$ are rotation matrices at the $k^{th}$ interface and is described in [30].

In the next section we briefly describe two approximate Riemann solvers namely the HLL and the HLLC that can be used to estimate the flux $\mathbf{F}(\mathbf{U_{L}},\mathbf{U_{R}})$ at a given interface. To avoid extra notations, henceforth $\mathbf{F}$ simply represent the local Riemann flux at any interface with outward pointing normal $\mathbf{n}$.

4 A review of two classic Riemann solvers

The following subsections briefs two classic Riemann solvers on which the present work builds upon.

4.1 The HLL Riemann solver

The HLL Riemann solver, devised by Harten, Lax and van Leer [4] assumes a wave structure consisting of two waves that seperates three constant states. These two waves represent the genuinely nonlinear waves like shocks and expansion fans. The average state $\mathbf{U}^{*}_{HLL}$ can be derived through application of integral form of conservation laws as,

[TABLE]

Using Rankine-Hugoniot conditions, the HLL Riemann flux can then be written as,

[TABLE]

Here $S_{L}$ and $S_{R}$ are numerical approximations to the speeds of the left most and right most running characteristics that emerge as the solution of Riemann problem at the cell interface of interest. It has been shown that under appropriate choice of wavespeeds $S_{L}$ and $S_{R}$, the HLL scheme is both positivity preserving and entropy satisfying [5]. In this work, $S_{L}$ and $S_{R}$ are computed as [5],

[TABLE]

where $u_{nL,R}$ are the normal velocity across an interface, $a_{L,R}$ are the respective sonic speeds and $\tilde{u}_{n},\tilde{a}$ are the standard Roe averaged quantities at the interface [1]. Using these wavespeed estimates, the HLL scheme is also known as the HLLE scheme [5]. It has been shown that under these wave speed, the HLLE scheme is both positivity preserving and entropy satisfying. Although the HLLE scheme is very robust for shock capturing applications, one of its major drawback is its inability to resolve the linearly degenerate waves of the Euler system that govern the contact and the shear phenomenon. Such an inability severely restricts the scheme from being used in a full Navier-Stokes solver because linear wave resolution is generally regarded as a prerequisite for accurate resolution of viscous phenomenon [31, 32]. Gressier et al [9] have conjectured that upwind schemes that posess exact contact ability cannot be free from shock instabilities. Hence by corollary, the shock stable behaviour of the HLLE scheme is attributed to its linear wave dissipative nature [11].

4.2 The HLLC Riemann solver

Toro et al [3] mitigated the short coming of the HLL Riemann solver by introducing the contact and shear preserving feature into it. The resulting solver is called the HLLC ( C for Contact) scheme. This improvement was achieved by adding a third wave to the one dimensional HLL wave structure, called a contact wave with speed $S_{M}$. Once again by using the integral form of conservation laws and appropriate flow conditions relating to the middle wave, closed form expressions for the middle state quantities $\mathbf{U}_{*L/R}^{HLLC}$ can be derived. Based on this the HLLC interface flux $\mathbf{F}_{HLLC}$ can be written as,

[TABLE]

The two intermediate states $\mathbf{U}_{*L}^{HLLC}$ and $\mathbf{U}_{*R}^{HLLC}$ that exist across the additional wave $S_{M}$ are obtained as,

[TABLE]

where $u_{tL/R}$ denote the local tangential velocity at an interface. In the above expressions, $S_{L}$ and $S_{R}$ can be obtained using Eq.(21). Batten et al [33] provides a closed form expression for $S_{M}$ as,

[TABLE]

The HLLC Riemann solver is one of the simplest known Riemann solver to be able to resolve all the waves of the Euler system accurately. Like the HLL Riemann solver, it is also an entropy satisfying and positively conservative scheme under appropriate wave speed selection [3, 33]. However, albeit all the attractive features the scheme posesses, the HLLC Riemann solver is known to be highly susceptible to numerical shock instabilities. In the next section, a simple and inexpensive method of protecting the HLLC scheme from the shock instability phenomenon is presented. The idea involves first splitting the standard HLLC scheme into its diffusive HLL component and an antidiffusive component responsible for restoring linear wave accuracy to this HLL component. Then, a multidimensional solution dependent strategy is used to enhance the magnitude of this diffusive component in the vicinity of shocks through simple modifications of certain carefully chosen nonlinear wavespeed estimates. The resulting solver is found to be completely free of numerical shock instabilities yet retaining the accuracy and efficiency of the original HLLC solver.

5 Formulation of a shock stable HLLC scheme

We propose the following strategy to save the HLLC scheme from instabilities. We first decompose the HLLC scheme into its diffusive component and the antidiffusive component responsible for restoration of accuracy on contact and shear waves. To achieve this we propose to rewrite the HLLC scheme as,

[TABLE]

where $\mathbf{F}_{HLL}=\mathbf{F}_{L/R}+S_{L/R}(\mathbf{U}^{*}_{HLL}-\mathbf{U}_{L/R})$. It is seen that $\mathbf{F}_{HLL}$ can be further recast into a second order central plus upwind dissipation form as,

[TABLE]

where the numerical dissipation of the embedded HLL scheme is [34],

[TABLE]

This can also be written in an equivalent form inspired from [35] as,

[TABLE]

Rewriting HLLC scheme in the form presented in Eq.(29) helps in clear identification of the diffusive and antidiffusive components that make up the scheme. While the overall disspation vector $\mathbf{D}_{HLL}$ occurs as a part of the HLL component of the flux, the accuracy on contact and shear waves that makes HLLC more accurate than the HLL solver is provided through the antidissipation vectors $S_{L/R}(\mathbf{U}^{*HLLC}_{L/R}-\mathbf{U}^{*HLL})$ in each case. Once the inherent dissipation vector $\mathbf{D}_{HLL}$ is identified, a strategy to increase the magnitude of it is sought near shock waves. Davis [36] has noted that choice of the wavespeeds $S_{L}$ and $S_{R}$ can affect the amount of dissipation that is introduced into the HLL scheme. Hence it is envisioned that additional dissipation to deal with shock instabilities can be introduced into the HLLC scheme through careful manipulations of the wavespeeds $S_{L}$ and $S_{R}$ present in $\mathbf{D}_{HLL}$. In the present work, a modified dissipation vector $\overline{\mathbf{D}}_{HLL}$ is proposed,

[TABLE]

where,

[TABLE]

$S_{L}$ and $S_{R}$ in the above equations are obtained from Eq.(21). The term $\epsilon^{L,R}$ is an appropriate choice of dissipation of the order $\mathcal{O}(S_{L},S_{R})$ and $\alpha$ is a coefficient through which the quantity of $\epsilon^{L,R}$ introduced into $S_{L},S_{R}$ can be controlled. Both $\epsilon^{L,R}$ and $\alpha$ are by definition nonnegative quantities and will be defined completely later. It can be seen intutitvely that modification of wavespeeds in Eq.(34) amounts to spreading of the wave structure of the inherent HLL scheme thereby increasing the magnitude of $\mathbf{D}_{HLL}$ which will reflect as an increase in the overall dissipation of the HLLC scheme. Eqs.(29,30,33 and 34) describe the proposed modified HLLC scheme.

Note that the contact wave speed $S_{M}$ and the wavespeeds $S_{L},S_{R}$ associated with the antidiffusion component are left completely unaffected by this modification. This is done purposefully to avoid affecting the contact and shear capturing ability of the HLLC solver. It is important to emphasize here that although the technique proposed here aligns with the idea of wave speed modification to achieve shock stability as recommended in [17], it differs from it in that instead of targeting eigenvalues of the linear fields, the dissipation is added through the wave speed estimates of the nonlinear fields. This can be justified by the fact that in the two wave formulation of the HLL scheme, $S_{L}$ and $S_{R}$ are the numerical characteristics associated with the shock waves and hence any dissipation infused through them would directly affect shock capturing ability only. Owing to this selective wave modification strategy the proposed modified version of the HLLC scheme will be henceforth refered to in this work as HLLC-SWM (Selective Wave Modified) scheme. It is emphasized that the modification presented in Eq.(34) must be used in conjuction with the form of dissipation vector $\overline{\mathbf{D}}_{HLL}$ presented in Eq.(33) for it to be effective. In the following subsection, we provide an analytical proof showing how the modification proposed in Eq.(34) manages to improve the overall dissipation of the HLLC scheme.

5.1 Effect of $\epsilon^{L,R}$ on numerical dissipation of the HLLC scheme

In this section we aim to analytically demonstrate how the $\epsilon^{L,R}$ parameter introduced through modified wavespeeds in Eq.(34) into the dissipation vector form in Eq.(33) amounts to overall increase in dissipation of the HLLC scheme. Although the explanation is provided through an expression of $\mathbf{F}_{HLLC}$ corresponding to condition $S_{L}\leq 0\leq S_{M}$, it remains symmetrically valid for the case when $S_{M}\leq 0\leq S_{R}$. As observed previously from Eq.(29), the HLLC flux can be interpreted as a combination of a diffusive HLL flux plus an antidiffusive component responsible for accuracy on contact and shear waves as,

[TABLE]

The numerical dissipation vector $\mathbf{D}_{HLL}$ of the this inherent HLL component can be expressed as in Eq.(35),

[TABLE]

Let $a_{0}$ and $a_{1}$ be coefficients of $\mathbf{D}_{HLL}$ defined as,

[TABLE]

Let $\overline{a_{0}}$ and $\overline{a_{1}}$ be the corresponding coefficients of the modified dissipation vector $\overline{\mathbf{D}}_{HLL}$ given in Eq.(33)

[TABLE]

For the sake of convenience, in the above expressions, the value of $\alpha$ is taken to be 1. Since in subsonic flows, $S_{L}<0$ and $S_{R}>0$, it becomes immediately evident that the denominator of $\overline{a_{0}}$ and $\overline{a_{1}}$ will be always postive. Further, coefficients in Eq.(37) can be simplified as,

[TABLE]

Some interesting observations are forthcoming at this point. It is clear from Eq.(38) that both coefficients $\overline{a_{0}}$ and $\overline{a_{1}}$ turn out to be respectively $a_{0}$ and $a_{1}$ with additional terms that are functions of the dissipation parameter $\epsilon^{L,R}$. If we choose $\epsilon^{R}>\epsilon^{L}$ then these terms are positive additions and we can ensure that $\overline{a_{0}}>a_{0}$ and $\overline{a_{1}}>a_{1}$. Note that such a constraint on $\epsilon^{L,R}$ can be imposed simply by setting say, $\epsilon^{L}\sim|S_{L}|$ and $\epsilon^{R}\sim|S_{R}|$ in a subsonic flow. Consider a simple choice $\epsilon^{L}=\epsilon^{R}=\epsilon$. Then the additional term in $\overline{a_{0}}$ vanish while the additional term in $\overline{a_{1}}$ becomes $\epsilon/2$. Thus for such a choice, it is seen that the modification presented in Eq.(34) effectively affects only the $\mathbf{U}_{L}-\mathbf{U}_{R}$ vector (ie coefficient $a_{1}$) and nowhere else in the entire HLLC formulation. This fact can be exploited to reduce the cost associated with implementation of the present cure to an already existing HLLC Riemann solver code. The above discussion clearly indicates that the modified HLL dissipation vector $\overline{\mathbf{D}}_{HLL}$ provides an excellent framework to introduce adequate dissipation into the original HLLC Riemann solver. However as seen above, the definition of $\epsilon^{L,R}$ has a significant role to play in the effectiveness of $\overline{\mathbf{D}}_{HLL}$.

In the next section we perform a linear perturbation analysis on the HLLC-SWM framework to reveal the intricate role of $\epsilon^{L,R}$ parameter in damping the numerical perturbations in flow field values which are known to be the prime instigators of these instabilities.

6 Linear perturbation analysis of HLLC-SWM scheme

A technique of predicting whether a numerical flux function will produce shock unstable solution based on how the flux function evolves an initially perturbed primitive variable profile was established by Quirk [10] and others [9, 11]. Quirk conjectured that schemes in which pressure perturbations feed into density perturbations are prone to produce shock instabilities. In this regard, it is interesting to use this technique to study the evolution of primitive variable perturbations using the HLLC-SWM scheme. Of specific interest is in discerning the effect of the dissipation parameters $\epsilon^{L,R}$ and $\alpha$ on this evolution. A steady mean flow with normalised state values of $\rho_{0}=1$, $u_{0}\neq 0$, $v_{0}=0$ and $p_{0}=1$ is chosen as the base flow. On a typical y-directional stencil the cells are marked as “even ”and “odd ”[11]. The perturbations are introduced as a saw tooth profile in density, x-velocity and pressure variables. For a typical cell ’j’ belonging to this stencil, these are initialized as,

[TABLE]

Such saw tooth profiles are characteristic of most shock unstable solutions. The crux of the analysis is to develop equations that describe the temporal evolution of these saw tooth perturbations described in Eq.(39). Note that we avoid introducing perturbations $\hat{v}$ in y-directional velocity in the present analysis because it has been reported that effect of this component is to dampen out perturbations in other components eventually stabilizing the scheme [11, 18]. Had it been included, then isolating the effect of parameter $\epsilon^{L,R}$ on physical quantitites would have been difficult.

For the HLLC scheme, the evolution equation for $\hat{\rho},\hat{p}$ and $\hat{u}$ can be derived and arranged in a matrix form as,

[TABLE]

Corresponding evolution equations for the HLLC-SWM scheme can be derived as,

[TABLE]

where $\lambda=\sqrt{\gamma}\frac{\Delta t}{\Delta y}$ denotes a linearized CFL number. The eigenvalues of the coefficient matrices in these cases represent the amplification factors of the respective perturbation quantities. It is seen from these equations that in case of the HLLC scheme, the amplification factors for perturbations in $\hat{\rho},\hat{u},\hat{p}$ are respectively $(1,1,1-2\lambda)$. A von Neumann type stability bound on $\lambda$ can be obtained from these as,

[TABLE]

It becomes evident that the perturbations in $\rho$ and $u$ are not damped by the HLLC scheme while those in $p$ are attenuated for $0<\lambda<1$. These behaviours can be observed in Figs.(2(b),3(b), 3(d)) that form the right panel of Fig.(3). Each plot in these figures indicate the evolution of perturbation in either density, x-velocity or pressure (indicated by $\hat{\rho},\hat{u},\hat{p}$ respectively) for a given initial perturbation in one of those quantities (while setting the other initial perturbations to be 0). These plots correspond to a value of $\lambda=0.2$. However, as seen from Fig.(3(d)), although $\hat{p}$ gets attenuated with time, any finite $\hat{p}$ is repeatedly feeding into $\hat{\rho}$ which results in existence of a residual $\hat{\rho}$ that is growing in time. These undamped perturbations in $\rho$ and $u$ may eventually amplify leading to unphysical mass flux jumps across the shock wave. Such jumps could then force a distortion of the shock profile causing the typical ’bulge’ in shock structure and an eventual shock instability solution [18]. Since the HLLC scheme is known to produce shock instabilites, it can be inferred from these observations that Quirk’s conjecture is perfectly valid in its case. Gressier et al [9] have termed schemes with atleast one undamped perturbation in physical variable (ie. with non-decaying amplifcation factors) as marginally stable schemes. By due consideration the HLLC scheme is clearly a marginally stable scheme.

Examining Eqs.(71) and the left panel of Fig.(3) in the light of above arguments, it is easy to see why the HLLC-SWM framework will result in a shock stable scheme. First of all we notice that all the three evolution equations of the HLLC-SWM scheme essentially comprise of a HLLC component and extra terms that contain the parameters $\epsilon^{L,R}$ and $\alpha$. The amplification factors of the perturbations $\hat{\rho},\hat{u},\hat{p}$ are respectively $\left(1-\frac{2\lambda\alpha\epsilon^{L,R}}{\sqrt{\gamma}},1-\frac{2\lambda\alpha\epsilon^{L,R}}{\sqrt{\gamma}},1-2\lambda(1+\frac{\alpha\epsilon^{L,R}}{\sqrt{\gamma}})\right)$. These lead to the following stability bound on $\lambda$,

[TABLE]

For a $\lambda$ in these bounds, it can be seen that any value for $\epsilon^{L,R},\alpha>0$ will introduce a damping factor of $\frac{2\lambda\alpha\epsilon^{L,R}}{\sqrt{\gamma}}$ into each of these perturbations. Notice that while this additional factor damps $\hat{u}$, it provides an ancillary damping to $\hat{p}$ over the preexisting damping coefficient $(1-2\lambda)$. In case of $\hat{\rho}$, the additional damping factor not only ensures that independent perturbations in density are attenuated, but also that any residual perturbations arising in $\hat{\rho}$ due to action of $\hat{p}$ are also suppressed. Thus in general, the additional factor ensures that the perturbations in all the three primitive quantities die down in time. These behaviours can be clearly observed in left panel of Fig.(3). The plots correspond to $\lambda=0.2$ and $\alpha,\epsilon^{L,R}=1.0$. Given these damping characterisitcs, the HLLC-SWM scheme is very unlikely to produce variation of mass flux in the vicinity of the shock profile. This would indicate that numerical solutions of the HLLC-SWM scheme would remain shock instability free [18]. Such a scheme with decaying amplification factors are termed as *strictly stable *scheme and are known to be shock instability free [9]. Thus modified dissipation $\overline{\mathbf{D}}_{HLL}$ proposed in Eq.(33) thus helps convert the HLLC scheme from its inherent marginally stable behaviour to a strictly stable behaviour.

Although various authors [17, 26, 27] have successfully used eigenvalue modification for the Roe scheme to cure shock instabilities it is not entirely made clear as to how these modifications affected the evolution of errors in physical variables during a computation to eventually secure a shock stable solution. The analysis presented here clarifies this in the context of the HLL based schemes. As seen from Eq.(71), the wavespeed modification suggested in Eq.(34) for the HLLC scheme expressed in the form shown in Eq.(29) directly affects the amplification factors that control the growth or decay of perturbations in physical quantities. A damping of these perturbations is assured through the dissipation parameters $\epsilon^{L,R}$ and $\alpha$ which eventually promises shock stability.

In the next section, we discuss the some strategies to estimate the quantity of the dissipation parameter $\epsilon^{L,R}$ to be used while calculating the HLLC-SWM flux at any interface.

7 Estimation of dissipation parameter $\epsilon^{L,R}$

Several choices for estimation of $\epsilon^{L,R}$ can be designed. Since the purpose of $\epsilon^{L,R}$ is primarily to increment the value of $\overline{S_{L}}$ and $\overline{S_{R}}$ locally in the vicinity of a shock wave and tend to zero elsewhere, any viable shock sensor is a good candidate for it. Since shock instability is primarily a multidimensional phenomenon, it is also advisable to have an estimation of $\epsilon^{L,R}$ that takes this fact into account [17]. Below we discuss some strategies to achieve a satisfactory estimate of $\epsilon^{L,R}$ based on these considerations.

7.1 Eigenvalues based estimation of $\epsilon^{L,R}$

One possible estimate for $\epsilon^{L,R}$ can be based upon the idea of an entropy fix. Several authors [11, 17, 26, 27] have used such entropy fix based estimates to save the Roe scheme from shock instabilities. For a typical cell interface $(i,j+1/2)$, $\epsilon^{L,R}$ may be defined as,

[TABLE]

where each $\eta_{k}$ ($k$=1…4) corresponds to the dissipation quantity calculated at predefined interfaces that form a stencil around $(i,j+1/2)$ as shown in Fig.(4).

Each $\eta_{k}$ can be evaluated using a common entropy fix formula,

[TABLE]

where $p$ iterates over all the characteristic wavespeeds $\lambda_{p}$ of the Riemann problem at $k^{th}$ interface. An $\epsilon^{L,R}$ at any interface defined in the above form represents the largest one dimensional entropy correction of all the associated interfaces considered around it. As pointed out in [11] this definition is designed carefully to avoid extraneous dissipation from being evaluated when computing flow features other than shocks with naturally occuring gradients; like the boundary layers. This is primarily achieved by excluding the $({i,j+1/2})^{th}$ interface itself from the stencil used to calculate $\epsilon^{L,R}$. Further, it is important to state here that the idea of an entropy fix is utilized in the above strategy only to determine an estimate for $\epsilon^{L,R}$ that is $\mathcal{O}(S_{L},S_{R})$. Since the HLLC-SWM framework builds upon the HLLC scheme which by default is entropy satisfying, it does not require fixes to ensure entropy satisfying solutions. It is well known from the works of many authors [12, 17, 16, 18, 21, 37] that one of the reasons for shock instability appearing in solutions of contact and shear preserving schemes has to do with its inability to provide adequate multidimensional dissipation in the vicinity of a shock front. The above definition of $\epsilon^{L,R}$ along with modification proposed in Eq.(34) may be considered as an effective strategy to address this concern in case of the HLLC scheme. The dissipation provided through $\epsilon^{L,R}$ has a multidimensional flavor because the dissipation introduced at a longitudinal interface is a function of the solution variation in the transverse direction and vice versa [17]. An HLLC-SWM framework that uses this estimate for $\epsilon^{L,R}$ is hereby termed HLLC-SWM-E (Eigenvalue based) scheme.

7.2 Pressure based estimation of $\epsilon^{L,R}$

Although quite appealing, the entropy fix based estimation of $\epsilon^{L,R}$ discussed above may cause certain loss of accuray in shear dominated viscous flows. This is due to the fact that $\epsilon^{L,R}$ in this definition is employed both as a shock sensor and as an estimate of the quantity of dissipation to be provided at it. Thus, in flows with strong velocity gradients, the value of $\epsilon^{L,R}$, instead of tending towards zero within the shear layers, could become a non zero quantity thereby introducing erroneous dissipation into the scheme and compromising the accuracy of the solution. An interesting alternative to this would be to use seperate functions for shock detection and estimation of dissipation quantity. This can be done by including a pressure based shock sensor in the definition of $\epsilon^{L,R}$. The advantage of pressure based sensor is that they become active at shocks and remain inactive in viscous dominated flows and inviscid contacts. Hence we define a new $\epsilon^{L,R}$ through a pressure based shock sensor as,

[TABLE]

where $\omega_{i,j+1/2}$ is a pressure based shock sensor that can be constructed as [38],

[TABLE]

where $f_{k}$’s are pressure ratio based functions evaluated on a predefined stencil around $({i,j+1/2})$ interface shown in Fig.(4). At $k^{th}$ interface, $f_{k}$ is defined as,

[TABLE]

Here $p_{R}$ and $p_{L}$ denotes the right and left cell center pressures across $k^{th}$ interface. A value of $\beta=5.0$ is used. Note that this pressure based $\epsilon^{L,R}$ also retains the merits of the shock sensor discussed in Sec.(7.1) ie, a multidimensional stencil and local solution dependence. An HLLEM-SWM scheme that uses this kind of $\epsilon^{L,R}$ will hereby be addressed as HLLC-SWM-P (Pressure based) scheme. In the next section, we complete the description of HLLC-SWM framework by estimating an appropriate value for the parameter $\alpha$.

8 Estimation of dissipation parameter $\alpha$

The dissipation parameter $\alpha$ controls the extent to which the wave speeds $S_{L}$ and $S_{R}$ are modified by the $\epsilon^{L,R}$ parameter. Too large a value of $\alpha$ could result in larger magnitudes of $\overline{S_{L}}$ and $\overline{S_{R}}$ which as seen from Eq.(38) will result in larger magnitude of $\overline{\mathbf{D}}_{HLL}$ and eventually threaten to spoil the accuracy of the solution. On the other hand, too small a value could fail to provide adequate numerical disspation to counter the threat of shock instability. Thus, choice of $\alpha$ plays a crucial role in the overall efficiency of the HLLC-SWM class of schemes and hence must be carefully chosen. We envision that a reasonable lower bound on $\alpha$ that is capable of dealing with shock instabilities can be obtained by studying the effect of varying $\alpha$ on the stability of a shock profile of a simple shock instability test problem: the steady simulation of a two dimensional thin shock. Ideally, the steady state solution of this problem will consist of the original unperturbed shock profile that will be preserved forever. However, presence of shock instability would immediately manifest as unsteady perturbations in physical quantities in the aftershock region. Instead of directly resorting to cumbersome numerical simulations to achieve this objective, we utilize the technique of matrix based stability analysis as proposed by Dumsber et al [7]. The present work chooses a $M=7$ two dimensional thin steady shock to be located between the 6th and 7th cells of a regular Cartesian grid comprising of 11 $\times$ 11 cells on a 1.0$\times$ 1.0 domain. The initial conditions on the supersonic side of the shock are chosen as $(\rho,u,v,p)_{L}=(1.0,1.0,0.0,0.01457)$ while Rankine-Hugoniot conditions are used to determine the corresponding conditions on the subsonic side. Periodic conditions are applied on the top and bottom boundaries. Since the HLLC scheme exactly preserves a thin shock, the initial and boundary conditions can be exactly specified [7]. The eigenvalue problem is solved using the eigenvalue function eig() from linalg package in python version 2.7.6.

8.1 Estimation of dissipation parameter $\alpha$ for the HLLC-SWM-E scheme

In this section we estimate the value of $\alpha$ that ensures a shock stable HLLC-SWM-E scheme. Fig.(5(a)) shows the variation of the max Real part of eigenvalues $\theta$ (indicated in the boxes) as a function of varying $\alpha$ values computed using the HLLC-SWM-E scheme.

Since $max(Re(\theta))$ represents the largest error growth rate, all $max(Re(\theta))$ value in Fig.(5(a)) lying above 0.0 corresponds to unstable (exponentially amplifying) solutions while those lying below represents shock stable (exponentially decaying) solutions on this problem. An $\alpha=0.0$ corresponds to the unmodified HLLC scheme which is known to be unstable on this problem. The positive $max(Re(\theta))$ value of +8.15562 associated with $\alpha=0.0$ confirms that the numerical errors have a propensity to grow and eventualy spoil the solution. It is interesting to note that even $\alpha=1$ does not ensure complete shock stability for the HLLC-SWM-E scheme. This means that unlike suggestions in [11, 17], for the Roe scheme, where an estimate of $\epsilon^{L,R}$ discussed in Sec.(7.1) was claimed to ensure shock instability, an $\alpha>1.0$ is necessary to save the HLLC scheme. Further, it is interesting to note that beyond $\alpha=1.0$, a slight increase in $\alpha$ is seen to stabilize the solution in the linear limit. With increasing $\alpha$, the maximum error growth rate is found to drop indicating a tendency to establish a shock stable solution.

We verify these predictions in case of the HLLC-SWM-E scheme using a numerical simulation of the above problem carried out using a Finite volume first order accurate Euler code. Fig.(6) shows plots of cell centered density values versus the y locations extracted right behind the shock, ie. from the 7th coloum of cells in the computational domain. Since the problem has an exact steady state solution, any deviation of value of the density variable from its exact steady value of 5.44 clearly indicates presence of instabilities. The solutions are obtained at t=20 units.

Fig.(6(a)) clearly shows the saw tooth profiled density variation that is typical of a shock unstable solution [10]. These undulations in physical variables may enhance or decay depending on parameters like grid aspect ratio or inflow Mach number of the specific problem considered. Further, as predicted by the analysis, increasing $\alpha$ values are seen to stabilize the solution. For $\alpha\geq 2.5$ the HLLC-SWM-E scheme is able to infuse sufficient dissipation to render a shock stable solution on this problem. However as reported in [7], inflow Mach number can have a predominant effect on the proliferation of these instabilities. Hence, a choice of $\alpha$ that is capable of dealing with instabilities over a wide range of Mach numbers is preferred. To include the effect of Mach numbers, the above analysis is repeated for three additional Mach numbers; M=3, M=10 and M=20. Fig.(5(b)) shows the variation of $max(Re(\theta))$ with respect to $\alpha$ for these additional Mach numbers.

It is seen from Fig.(5(b)) that in general, higher Mach numbers corresponds to larger error growth rates for any value of $\alpha$. For $\alpha=0.0$ which corresponds to the HLLC scheme, increase in Mach numbers leads to a drastic decrease in theoretical shock stability. This is an expected behaviour from a typical shock instability prone scheme like the HLLC [7]. Further at any given Mach number, an increasing value of $\alpha$ is seen to reduce the error growth rate. In case of the HLLC-SWM-E scheme, each Mach number has a corresponding $\alpha$ value at which the $max(Re(\theta))$ crosses over into the stable region. From Fig.(5(b)) it is seen that choice of $\alpha=3.5$ favors a linearly stable solution over the range of Mach numbers investigated.

8.2 Estimation of dissipation parameter $\alpha$ for the HLLC-SWM-P scheme

Fig.(7) shows the variation of $max(Re(\theta))$ with $\alpha$ for the stationary shock problem computed using the HLLC-SWM-P scheme. It is evident that the behaviour of the stability region for the HLLC-SWM-P scheme is quite similar to that of the HLLC-SWM-E scheme. A value of $\alpha=3.5$ is seen to assure theoretical shock stability. We have confirmed this prediction with numerical experiments similar to that done in Sec.(8.1).

We present the eigenvalue spectrums for M=7 stationary normal shock computed using the HLLC-SWM schemes in Fig.(8). An $\alpha=3.5$ is used for both schemes. Corresponding specturm for the original HLLC is also provided for comparison. The eigenvalue spectrum visibly shows how the proposed HLLC-SWM-E and the HLLC-SWM-P schemes utilize the modified dissipation $\overline{\mathbf{D}}_{HLL}$ to shift the unstable error growth rates to the stable region. While the HLLC scheme has the largest error growth rate of +8.15562, the HLLC-SWM-E and the HLLC-SWM-P schemes with $\alpha=3.5$ are able to bring this down to -0.66053 and -0.59317 respectively. Amongst the HLLC-SWM schemes, the HLLC-SWM-E scheme turns out to be more stable. Further, it is important to note that the shock stability predicted by the matrix analysis for the proposed HLLC-SWM schemes are independant of the time discretization technique adopted or the CFL number employed because only spatial discretization of the governing equations has been used to construct the stability matrix.

9 Numerical verification of the HLLC-SWM schemes

In this section, we demonstrate the efficacy of the proposed HLLC-SWM schemes through a series of strict test problems on which the original HLLC scheme is known to produce shock unstable results. It has been reported in [9] that the intensity of instability manifestation in certain test problems reduces as the order of accuracy of the computation increases; although it does not completely remove it. Based on this, certain test cases are computed as plain first order while second order solutions are sought for others to clearly showcase the improvement in the computed flow field using the proposed schemes. Second order spatial accuracy in primitive variables is achieved by limiting gradients obtained using Green Gauss method [39] with Barth Jersperson limiter [40]. Second order time accuracy is achieved using strong stability preserving variant of the Runge Kutta method [41]. All boundary conditions are set using ghost cells. A value of $\alpha=3.5$ is used to configure both the HLLC-SWM schemes.

9.1 Positivity tests

Batten et al [33] has proved analytically that the original HLLC scheme is a positively conservative scheme and that it passes the tough numerical test cases found in [2]. This section investigates whether the modification of nonlinear wavespeeds proposed in Eq.(34) adversely affects this important property of the HLLC scheme. Computations are done on 100$\times$50 structured Cartesian mesh on a domain of size 1.0$\times$1.0. First order solutions are sought and the CFL is chosen as 0.5. The first test case is a Riemann problem with initial conditions $(\rho,u,v,p)_{L}=(1.0,-2.0,0.0,0.4)$ and $(\rho,u,v,p)_{R}=(1.0,2.0,0.0,0.4)$. The boundaries are maintained as open using simple extrapolation. The solution to this problem describes two receding expansion waves that originate from the center of the domain and move toward the boundary.

Fig.(9) shows the result of this test case where both density and total energy are plotted along the length of the tube after $t=0.1$. Exact solution to this problem can be computed using the exact Riemann solver of Godunov [42] and is also plotted in the same figure for reference. It is clearly seen that both HLLC-SWM schemes are able to compute the solution without failure and admits physically acceptable density and total energy profiles.

The second Riemann problem is initialized using the condition $(\rho,u,v,p)_{L}=(1.0,-1.0,-2.0,5.0)$ and $(\rho,u,v,p)_{R}=(1.0,1.0,2.0,5.0)$. The solution of this problem consists of a pair of receding expansion fans seperated by a shearing flow in y-direction. Fig.(10) shows the solution to this problem computed using the HLLC-SWM schemes shown at $t=0.11$. For reference, the solution computed by the HLLC scheme at the same time is also provided in Fig.(10). It is seen that the HLLC-SWM schemes are able to compute the solution to this problem too without failing. The solutions computed by the HLLC-SWM schemes are atleast as accurate as the HLLC scheme. Solutions to these demanding test problems suggest that the modification proposed in Eq.(34) seems to have no adverse effect on the positivity property of the HLLC solver and that both HLLC-SWM schemes are positively conservative in nature. It is also observed that choice of value of $\alpha$ does not affect this property.

9.2 Moving shock instability problem (Inviscid)

Shock instability that occurs in a M=6 moving shock propagating down a computational tube filled a stationary fluid defined by $(\rho,u,v,p)=(1.4,1.0,0.0,1.0)$ was first reported by Quirk [10]. The instability is triggered by perturbing the centerline grid of the computational domain to an order of 1E-6. First order solution is sought. The CFL for the computations were taken to be 0.5 and simulations were run till shock reached a location x=650. The results showing 50 density contours equally spanning values from 1.4 to 7.34 is shown in Fig.(11).

Fig.(11(a)) clearly shows the deteriorated condition of the shock profile calculated by the original HLLC scheme. Figs.(11(b),11(c)) demonstrate that both the HLLC-SWM-E and the HLLC-SWM-P schemes are quite capable of suppressing the instabilities and maintaining a clean shock profile.

9.3 Inclined stationary shock instability problem (Inviscid)

This problem was reported in [43] and is a variant of the standing shock instability problem discussed previously. The initial shock wave of strength M=7 is set up along the line y=2(x-12) at an angle of $63.43^{o}$ with respect to the global x-direction. The computation mesh consists of 50$\times$30 cells uniformly spanning a domain of size 50.0$\times$30.0. The left and right boundaries are maintained as supersonic inlet and subsonic outlet respectively while periodic conditions are applied to top and bottom walls. CFL of the computation was chosen as 0.5 and first order solution was sought at t=5. Fig.(12) shows the result of this experiment where thirty density contours uniformly spanning 1.0 to 5.44 are shown.

Fig.(12(a)) shows the presence of instability in the HLLC scheme. Fig.(13(a)) shows the density variation along the shock front extracted just behind the shock. The density variation has a typical saw tooth profile. Fig.(12(b)) shows the density contours computed by the HLLC-SWM-E scheme. It was noticed that although till t=2 units there was no trace of instability in the computed field, some sort of perturbations were observed in the simulation later. These perturbations however remained temporaly bounded and did not degrade the shock profile. The respective density variation is shown in Fig.(13(b)). It can be seen that the nature of these perturbations are not similar to the spatially oscillatory one observed for the HLLC scheme. Fig.(12(c)) shows the density contours computed by the HLLC-SWM-P scheme. The corresponding density variation is shown in Fig.(13(c)). It can be seen that introduction of the pressure sensor helps the HLLC-SWM-P scheme to outperform the HLLC-SWM-E scheme on this test case. The instability has been drastically minimized and the computed field is preserved very close to the exact solution.

9.4 Supersonic flow over a blunt body (Inviscid)

The steady state numerical solution of hypersonic flow over cylindrical body is a routine test problem to investigate the susceptibility of numerical schemes to a particularly dreaded form of shock instability: the carbuncle phenomenon. This failure is characterized by visible distortions that appear in the expected smooth bow shock profile in front of the body. In this work, a cylindrical body of unit radius is subjected to a M=20 flow. We use 320$\times$20 body fitted structured quadrilateral cells in circumferential and radial directions respectively. Detailed problem description can be found in [22]. The computation of this problem is carried out with plain first order accuracy. The CFL for the computations were taken to be 0.5 and simulations were run for 30,000 iterations. The results showing twenty density contours equally spanning value from 1.4 to 8.5 is shown in Fig.(14).

Fig.(14(a)) clearly shows the onset of instabilities in the shock profile calculated by the HLLC scheme. Figs.(14(b), 14(c)) on the other hand demonstrates that the HLLC-SWM schemes are quite capable of computing clean shock profiles and post shock subsonic reions which are completely free of any instabilities. A quantitative review of the performance of these schemes is given in Fig.(15(a)) which shows the pressure value extracted along the stagnation line (j=160). Specific locations of interest on the plot are labeled as ‘A ’and ‘B ’and are shown magnified in Figs.(15(b)) and (15(c)).

Referring to Fig.(15(b)), at location ‘A ’, it is seen that the instability has rendered a much diffused shock in case of the HLLC scheme. In comparison, while both HLLC-SWM schemes produce crisper shocks, the HLLC-SWM-E scheme seems to be slightly better than HLLC-SWM-P scheme in predicting the after shock pressure value. Location ‘B ’which is close to the stagnation point clearly shows that both the HLLC-SWM schemes produce much accurate stagnation pressure values (which is closer to the expected inviscid value of 515.60) as compared to the HLLC scheme. This example shows the fidelty of the proposed HLLC-SWM schemes in computing Carbuncle free solutions on such severe grids.

9.5 Double Mach Reflection problem (Inviscid)

Development of kinked Mach stem in a situation of a double Mach reflection has been reported in various literatures [10, 44]. An oblique shock corresponding to $M=10$, making a $60^{o}$ angle with the bottom wall is made to propagate through a domain of size $4.0\times 1.0$ divided into $480\times 120$ structured Cartesian cells containing stationary fluid. Detailed problem description can be found in [44]. First order solution in time and space is sought for this case. The simulation is run till $t=0.02$ with CFL of 0.8. Fig.(16) shows results of the experiment with twenty five equispaced density contours spanning values from 1.4 to 21.

In Fig.(16(a)), the presence of kinked Mach stem and the subsequent triple point produced by the HLLC scheme is clearly visible. However, both the HLLC-SWM schemes produce clean shock profiles which are devoid of the kink and triple point.

9.6 Diffraction of a moving normal shock over a 900 corner (Inviscid)

The problem of a Mach 5.09 normal shock’s sudden expansion around a $90^{o}$ corner was studied in [10]. We use an unit dimensional domain meshed with 400 $\times$ 400 cells with the right angled corner where the initial normal shock stands is located at $x=0.05$, $y=0.45$. The problem is computed to second order accuracy. Simulations are run for $t=0.1561$ units with CFL of 0.4. Thirty density contours equally spanning values of 0 to 7 are shown in Fig.(17).

Fig.(17(a)) clearly demonstrates the extent of instability present in the computed shock profile produced by the HLLC scheme. A part of the normal shock in the right corner is completely distorted in this case. However, both the HLLC-SWM schemes are able to satisfactorily compute the solution without producing any trace of instability in the normal shock. It must also be noted that HLLC-SWM schemes does not produce any unphysical expansion shocks like the Roe scheme [10].

9.7 Supersonic flow over forward facing step (Inviscid)

A M=3 supersonic flow over forward facing step is another standard problem that has been extensively studied earlier in [44]. The computational domain is of size $3\times 1$ and is meshed with $120\times 40$ structured Cartesian cells. A 0.2 units high step is located at a distance of 0.6 units from the inlet. The problem is computed to second order accuracy. The simulations are run for t=4 with CFL as 0.5. Fig.(18) shows forty equispaced density contours spanning 0.2 to 7.0.

It is seen from Fig.(18(a)) that predominant oscillations are visible in the shock profile computed by HLLC scheme. These are most prominent at the normal shock stems of the primary shock near the bottom and top boundaries. Further, HLLC solution also has a severe kinked stem in the reflected shock. In comparison, solution computed by both HLLC-SWM scheme shown in Figs.(18(b),18(c)) shows clean shock profiles with no remnants of instabilities visible. The infamous kink has also improved considerably in these solutions.

9.8 Two dimensional shear flow (Inviscid)

The problem describes two fluids with different densities sliding at different speeds over each other and investigates the inviscid contact capturing ability of a given scheme [45]. The conditions for the top and bottom fluids are chosen as $(\rho,p,M)_{top}$ = (1, 1, 2) and $(\rho,p,M)_{bottom}$ = (10, 1, 1.1) respectively. A coarse grid consisting of 10$\times$10 cells on a domain of 1.0$\times$1.0 is used. CFL for the simulations are taken to be 1.0 and simulations were run for 1000 iterations. All simulations are plain first order accurate. The results showing fifty equispaced density contours spanning 1.0 to 10.0 is shown in Fig.(19). To clearly distinguish between the behaviours of contact preserving and contact dissipative schemes, the solution for the HLLE scheme is also given for comparison.

It can be seen from Fig.(19) that behaviours of the HLLC-SWM schemes closely matches that of the HLLC scheme. All these schemes are able to clearly preserve the inviscid contact. On the other hand, the HLLE scheme produces a very diffused density profile. Fig.(19(e)) shows variation of density along y-location at the center of the domain which confirms the observation made above. From the plot it can be inferred that the modified dissipation vector $\overline{\mathbf{D}}_{HLL}$ that introduces additional dissipation into the HLLC scheme to deal with shock instabilities has not affected the ability to resolve inviscid contact discontinuities.

9.9 Laminar flow over flat plate (Viscous)

The original HLLC scheme is well known for its ability to accurately resolve shear dominated viscous flow features including practically significant boundary layer flows. It is hence interesting to investigate how modified dissipation vector $\overline{\mathbf{D}}_{HLL}$ will affect the performance of the HLLC scheme on these types of flows. To perform this study, the test case of a $M=0.1$ subsonic laminar flow of air over a flat plate of length L=0.3048 m in atmospheric conditions is chosen. The total length of the domain is 0.381 m in x direction and 0.1 m in y direction. The domain is divided into 31$\times$33 Cartesian cells with 15 cells within the boundary layer. Viscous fluxes were discretized using simple averaging. CFL was taken to be 0.7. The flow was considered to have achieved steady state when the horizontal velocity residuals dropped to the order of $1e^{-7}$. The normalized longitudinal velocity profiles $(\frac{u}{u_{\infty}})$ are plotted against the Blasius parameter $\eta=y\sqrt{u_{\infty}/\mu L}$ in Fig.(20).

It is seen from Fig.(20) that the HLLC-SWM-E scheme suffers a slight loss of accuracy within the boundary layer. However close to the wall, the HLLC-SWM-E scheme is as accurate as the HLLC scheme in resolving the velocity gradients. On the other hand, the HLLC-SWM-P scheme is remarkably accurate everywhere within the boundary layer and exactly matches the Blasius profile. As discussed previously, the dissipative behavior of the HLLC-SWM-E scheme could be explained from the fact that since the $\epsilon^{L,R}$ in Eq.(34) computes the maximum difference in wavespeeds at any interface, it could be assuming a non zero value in the shear layer due to existence of gradients $\Delta_{x}(u,v)$ and $\Delta_{y}(u,v)$. This would then enlarge the magnitude of the $\overline{\mathbf{D}}_{HLL}$ term thus making the HLLC-SWM-E scheme more dissipative. To control this unwanted dissipation in the shear dominated low Mach number flows, a simple and inexpensive Mach number based switching function $f$ is defined as,

[TABLE]

where $M$ is the local flow Mach number. The function $f$ can now be integrated into Eq.(34) as,

[TABLE]

It is apparent that the switching function $f$ smoothly transits from its value of 1 for $M>0.3$ to its value of 0 as $M\rightarrow$ 0. This enables reduction of dissipation in the HLLC-SWM-E scheme thus recovering the full accuracy of the original HLLC scheme within the viscous layers. The results showing the effect of $f$ on HLLC-SWM-E is also shown in Fig.(20). It is evident that introduction of switching function clearly restores the accuracy of the solution.

10 Concluding remarks

We have proposed a simple cure for numerical shock instabilities in the HartenLaxvanLeer with contact (HLLC) approximate Riemann solver. Firstly we showed that the HLLC scheme can be split into its respective diffusive HLL component and an antidiffusive component responsible for accuracy on contact and shear waves. Once identified, we suggest a simple modificiation to certain nonlinear wavespeed estimates appearing only in the numerical dissipation of the HLL component that increases its magnitude. Two entropy fix based strategies for obtaining multidimensional dissipation in the vicinity of a shock wave was presented which resulted in the HLLC-SWM-E and HLLC-SWM-P variants of the proposed framework. Through a linear analysis of a mean flow with perturbations and a matrix based stability analysis of an isolated two dimensional stationary normal shock, the theoretical effectiveness of these schemes in computing shock stable solutions over a wide range of Mach numbers was established. Further, using numerical test cases we show that the HLLC-SWM schemes provide instability free, clean shock profiles on a variety of problems while retaining the positivity preserving and exact inviscid contact abilities of the underlying HLLC scheme. When computing the problem of a subsonic laminar flow over flat plate, it was found that the HLLC-SWM-P scheme was as accurate as the HLLC scheme while the HLLC-SWM-E scheme introduced slight extraneous dissipation into the solution. A simple strategy that involves a Mach number based switching function is proposed to recover the accuracy of the HLLC-SWM-E scheme on such problems. The success of these HLLC-SWM schemes in dealing with shock instability phenomenon alludes to the fact that cures for these spurious failures can be designed by introducing adequate multidimensional dissipation in the vicinity of a shock front. Lastly, the development of the HLLC-SWM-P scheme suggests that it is possible to circumvent the conjecture in Gressier et al [9] and create a robust shock stable upwind scheme that is accurate on contact and shear waves.