A carbuncle cure for the Harten-Lax-van Leer contact (HLLC) scheme using a novel velocity-based sensor

A hybrid numerical flux scheme is proposed by adapting the carbuncle-free modified Harten-Lax-van Leer contact (HLLCM) scheme to smoothly revert to the Harten-Lax-van Leer contact (HLLC) scheme in regions of shear. This hybrid scheme, referred to as the HLLCT scheme, employs a novel, velocity-based shear sensor. In contrast to the non-local pressure-based shock sensors often used in carbuncle cures, the proposed shear sensor can be computed in a localized manner meaning that the HLLCT scheme can be easily introduced into existing codes without having to implement additional data structures. Through numerical experiments, it is shown that the HLLCT scheme is able to resolve shear layers accurately without succumbing to the shock instability.


Introduction
Simulation of high-speed, compressible flows requires proper design of inviscid numerical flux functions to capture shock waves, contact discontinuities, and shear layers in a stable and accurate manner. Grounded in the seminal work of Godunov [1] , approximate Riemann solvers represent the state-of-the-art technique in designing flux functions whereby the inviscid fluxes are obtained by assuming an approximate solution structure for the Riemann problem across a cell interface. One of the simplest and most robust flux schemes is the Harten-Lax-van Leer (HLL) scheme [2] in which an incomplete two-wave solution structure is assumed. However, this incompleteness introduces non-physical dissipation of density and tangential momenta which precludes exact capturing of contact discontinuities and shear waves, resulting in highly diffused material interfaces and shear layers. Hence, it is unsuitable for use in many practical flow scenarios such as boundary layer flows [3] . Toro et al. [4] proposed a modification to the HLL scheme known as the HLL contact (HLLC) scheme [4] , which allowed exact capturing of the contact discontinuities and shear waves. The complete three-wave HLLC scheme has become one of the most widely used numerical flux schemes owing to its simplicity and accuracy. However, being a contact-/shear-preserving flux scheme, the HLLC scheme invariably succumbs to a multi-dimensional shock instability problem.
Quirk [5] presented the first systematic analysis for failings of approximate Riemann solvers, including the various manifestations of the shock instability problem such as carbuncle phenomena, kinked Mach stems, and odd-even decoupling. He hypothesized that strong, grid-aligned shocks provide a 'systematic perturbation' to the pressure field and shock-unstable flux schemes feed the pressure perturbations into the density field causing density perturbations to grow unbounded culminating in the instability. Pandolfi and D'Ambrosio [6] , through odd-even decoupling analyses, confirmed that flux schemes which account for contact surfaces explicitly are susceptible to shock instability. Xu and Li [7] mentioned that the shear-preserving ability of flux schemes provides insufficient dissipation along the shock front which leads to shock instability.
The conventional approach to overcoming the shock instability problem is to introduce additional dissipation in the vicinity of the shocks, usually by switching to a complementary shock-stable scheme. There have been several such attempts to hybridize the HLLC scheme in the recent years [8][9][10][11][12] . With the exception of the work of Huang et al. [11] who used the rotated Riemann solver formulation [13] with a localized velocity-based sensor, the remaining studies relied on some form of pressure-based shock sensor computed in a non-localized manner about the concerned cell interface. Among these studies, the work of Shen et al. [10] is unique and particularly illuminating. They designed a modified HLLC scheme called HLLCM (where 'M' refers to the modification) based on the Rankine-Hugoniot jump conditions which preserves contact waves but smears shear waves. In this respect, one may view it as a 'two-and-a-halfwave' Riemann solver with shear velocity dissipation across the contact wave. Using matrix stability analysis [14] , they showed that the HLLCM scheme is shock-stable. Then, they proposed a hybrid HLLC-HLLCM scheme by using a pressure-based shock sensor as mentioned before and demonstrated the effectiveness of the hybrid scheme for a variety of test cases.
Unlike the more common HLLC-HLL hybrids, the HLLC-HLLCM hybrid scheme uses a contact-preserving complementary scheme. This presents a unique opportunity to explore a radically new idea: to deactivate dissipation across shear layers rather than injecting additional dissipation along shock fronts. In other words, the HLLCM scheme is viewed as the base scheme which is applied everywhere except near shear layers where the HLLC scheme is used to resolve the shear layers accurately. The main objective of the present work is to explore the viability of this approach. In this approach, the role of the sensor now becomes one of determining regions of shear flow. To this end, a velocity-based sensor is introduced in this paper. This sensor can be computed in a localized manner as opposed to the non-local shock sensors used in most of the carbuncle cures reported in past studies [6,[8][9][10]12,[15][16][17][18] . The paper is organized as follows. The governing equations and discretization schemes are introduced in Section 2. The derivations of the HLL, HLLC, and HLLCM approximate Riemann solvers are presented in Section 3 along with suitable wave-speed estimations. The new localized, velocity-based sensor is proposed and discussed in Section 4. The results of a variety of numerical experiments are presented in Section 5. Finally, the concluding remarks are provided in Section 6.

Governing equations
The equations governing the dynamics of inviscid gas flow are given as where the state vector Q and flux vector F (Q) are defined by in which ρ is the gas density, u is the gas velocity, p is the gas pressure, and E is the specific total energy. The total energy comprises of the internal energy e and the kinetic energy k, which can be related to the other state variables through the ideal gas equation of state as follows: where γ is the ratio of specific heat capacities taken to be 1.4 in the present study.

Spatial and temporal discretization
In finite volume methodology (FVM), the domain is discretized into non-overlapping cells, and the solution is obtained in terms of the cell averages. The cell average of a function φ(x) over the ith cell Ω i is given by where |Ω i | refers to the volume of the cell. Applying the averaging operation to Eq. (1) and transforming the volume integral of the divergence term into a surface integral using the Gauss theorem, as is customary in FVM, yields an ordinary differential equation (ODE) for Q i as follows: where ∂Ω i refers to the boundary of cell Ω i . For a polyhedral cell, the surface integral on the right-hand side of Eq. (5) can be approximated as a summation over each of the planar faces that make up ∂Ω i as follows: Note that the physical flux F has been replaced by a numerical flux F in the above approximation. The numerical flux is evaluated at the face centroids using the left-and right-based interpolated state variables Q f,L and Q f,R . The left (L) and right (R) sides of a face are defined with respect to the face unit normal vectors n f which point out of Ω i as shown in Fig. 1. Lastly, S f is the face surface area.
In the present study, we restrict to the first-order upwind scheme for spatial interpolation since higher order schemes tend to be less susceptible to shock instability problem [6] . As for the temporal discretization, the explicit third-order total-variation-diminishing (TVD) Runge-Kutta scheme [19] was used for transient problems while an implicit matrix-free lower-upper symmetric Gauss-Seidel (LU-SGS) algorithm [20] was used for steady-state problems. With the spatial and temporal schemes defined, we turn our attention to the numerical flux function itself which is the main focus of this paper.

Review of HLLC and HLLCM schemes
In the first-order upwind scheme, the variables are assumed to remain constant within each cell, i.e., no spatial variation within the cell. Hence, a discontinuity arises at each cell interface f . Godunov [1] put forth the idea of treating each discontinuity as a Riemann problem and solving this Riemann problem exactly to derive the numerical flux F . Godunov's method is expensive for gas dynamics due to the iterative procedure involved in computing the exact solution structures which may consist of rarefaction fans, contact discontinuities, and shock waves, depending on the initial condition. Moreover, Godunov's method itself is not immune to the shock instability problem. One way to improve the efficiency is to replace the exact solution structure with an approximate one for which it is simpler to compute the intermediate states and flux. These are known as approximate Riemann solvers.
The schematic of a general 3-wave approximate Riemann solver is shown in Fig. 2. The approximate solution structure consists of three waves traveling at speeds s L < s * < s R in the direction normal to the face and two intermediate states Q L * and Q R * . For a clearer understanding of how each flux component is treated by the approximate Riemann solver, the state vector is written in the following form whereby the momentum equation is decomposed into the normal and tangential directions: From hereon, the variables u ≡ u · n f and v = u − un f refer to the normal and tangential velocities, respectively, with respect to the face normal direction.
The approximate solution is required to obey the following Rankine-Hugoniot (RH) jump conditions: Note that the intermediate fluxes F L * and F R * need not necessarily be equal to F (Q L * ) and F (Q R * ), respectively. Denoting F L = F (Q L ) and F R = F (Q R ) for conciseness, the final numerical flux F is chosen as follows: The HLLC and HLLCM schemes can be derived by assuming different forms for the intermediate states Q L * and Q R * . These assumptions are summarized in Fig. 3. The intermediate state for the HLLC scheme [4] is given by where α K ≡ ρ K (s K − u K ). Similarly, the intermediate state for the HLLCM scheme [10] is given as follows: The intermediate wave-speed s * can be obtained through the following closed form expression: Comparing Eqs. (10) and (11), it can be observed that the HLLC and HLLCM schemes behave identically when v L = v R , i.e., when the flow is locally one-dimensional. The two schemes differ in their treatment of the intermediate tangential velocity. The HLLCM introduces artificial dissipation to the intermediate tangential velocity while the HLLC scheme does not. This shear viscosity improves shock stability, but it also results in highly diffused shear layers.

Wave-speeds
In order to complete the definition of the flux schemes, the wave-speeds s L and s R are required. One possible choice proposed by Davis [21] is to compute the wave-speeds based on the minimum and maximum eigenvalues of the left and right states as follows: where c = γp/ρ is the acoustic speed. Another possible choice recommended by Einfeldt [22] is to use the eigenvalues of the Roe-averaged state as follows: where quantities with overhats (ˆ) refer to the respective Roe-averages. Batten et al. [23] recommend Eq. (14) for wave-speed estimation since it is less diffusive than Eq. (13) and it returns the exact shock velocity in the case of an isolated shock. Recently, Liu et al. [24] showed that the choice of wave-speed plays a significant role in the stability of the flux scheme and demonstrated that Eq. (13) provides better shock stability compared with Eq. (14). Therefore, Eq. (13) will be used in this study.

.1 Motivation
In the hybrid HLLC-HLLCM scheme proposed by Shen et al. [10] , the intermediate states are defined as Q Hybrid where w is a pressure-based switching function defined as follows: The sensor w must be evaluated in a non-local manner, considering not only the face f but also all its neighbors. For a structured 2D mesh, the set of neighboring faces to face (i+1/2, j) is shown in Fig. 4. Such non-local operations are unavoidable for determining regions near shocks, but they contribute to additional implementation costs. Therefore, it would be convenient to have a localized sensor instead. This provided the motivation to develop a localized 'shear sensor'. Unlike the original hybrid HLLC-HLLCM scheme in which dissipation is introduced in regions near shocks by switching from the HLLC scheme to the HLLCM scheme, the shear sensor may be used to eliminate dissipation in regions close to shear layers by switching from the HLLCM scheme to the HLLC scheme. The construction of such a shear sensor will be presented next.

A localized velocity-based shear sensor
The objective is to design a sensor that detects shear layers by using state values from the immediate neighboring cells of a face. The distinguishing characteristic of a genuine shear layer is that the streamwise velocity changes more rapidly compared with the cross-stream velocity across a shear layer. With this idea in mind, a localized shear sensor Θ could take the form where '∆' is the difference operator across the face, i.e., ∆φ ≡ φ R − φ L . When a face is oriented parallel to a grid-aligned shear layer exactly as shown in Fig. 5(a), u L = u R = 0. Therefore, the numerator of Θ vanishes and Θ = 0. Even when the shear layer is not perfectly grid-aligned, as long as the face is oriented approximately parallel to a shear layer as shown in Fig. 5(b), the tangential velocity changes much faster than the normal velocity. Hence, it is expected that |∆u| ∆v < ∆u across the face. Under these conditions, the numerator of Θ is much smaller than its denominator and, as a result, Θ is very small. To minimize the dissipation across genuine shear layers, the switching parameter w was taken to be where the constant τ was set to 0.5 based on numerical experiments. A small number 10 −12 was added to the denominator to avoid division by zero. Equation (18) performed well for several benchmark problems such as hypersonic flow past a blunt body [25] and double Mach reflection [26] problems commonly used to assess the shock stability of flux schemes. However, it was still inadequate in suppressing shock instability in problems where the instability manifested as odd-even decoupling, e.g., Quirk's test [5] . To overcome this deficiency, w was modified slightly by introducing a small reference velocity ε u as follows: We require ε u to increase faster than the velocity differences that may arise from numerical perturbations along a shock and yet remain several orders of magnitude smaller than the velocity difference across a genuine shear layer. Furthermore, we also require ε u to vanish at grid-aligned shear layers to preserve them exactly. Based on these considerations, ε u was computed as follows: The rationale behind this modification is as follows. Consider a face oriented parallel to the shock propagation direction. In the absence of any perturbations, all the properties remain the same on both sides of the face. Since ∆v = 0, the value of w is irrelevant. However, numerical perturbations do arise due to various reasons, and it has been found that the shock instability is linked to perturbations in the shock-normal momentum [27] which is the tangential momentum component with respect to the face. Hence, suppose that a slight perturbation is introduced only in the tangential velocity after the flow passes through the shock as shown in Fig. 6(a). Since pressure is computed based on the difference between the total and kinetic energies (see Eq. (3)), a pressure difference occurs between the two cells with a lower pressure in the cell with larger u and a higher pressure in the cell with smaller u . With u L = u R = 0 and w = 1 initially, the HLLC scheme is used. The HLLC scheme responds to the pressure difference (but not the tangential velocity difference which caused it) by inducing a mass flux across the face. However, the transfer of mass also results in the transfer of tangential momentum (ρvu) ∆p from the high-pressure cell to the low-pressure cell, thereby increasing the tangential velocity difference. As a result, the pressure perturbation continues to persist, and the process repeats.
The imbalance of pressure on the top and bottom faces of the cell changes the normal momentum in the cells, causing an increase in ε u . As the above-described cycle continues, ε u continues to grow. Recall that ε u was designed to increase faster than ∆u that arises from small numerical perturbations. As such, |∆u| , ∆u ε u , and w ≈ 0. Consequently, a tangential momentum flux (ρuv) ∆u is activated in the direction opposite to the pressure induced flux (ρuv) ∆p as shown in Fig. 6(b). This tangential momentum flux equalizes the tangential velocity difference between the cells bringing ∆v to approximately zero which, in turn, dampens the pressure perturbations and thereby ends the cycle as depicted in Fig. 6(c). In this manner, the proposed sensor alleviates the shock instability problem.
In contrast to the pressure-based switch, the localized velocity-based shear sensor provides a continuous blending of the HLLC and HLLCM schemes. This blended scheme will be referred to as the HLLCT scheme whereby 'T' indicates the modifications pertaining to the tangential momentum. Fig. 6 Restoring mechanism of the proposed sensor in regions susceptible to shock instability

Matrix stability analysis
Matrix stability analysis is a linear analysis technique introduced by Dumbser et al. [14] to characterize the shock stability of a numerical flux scheme. The analyses were performed for the HLLC, HLLCM, and HLLCT schemes, using wave-speed estimates from Eqs. (13) and (14). For the wave-speed estimates from Eq. (14), the fluxes were linearized about a grid-aligned M = 7 shock initialized on a dimensional 11 × 11 uniform Cartesian grid as follows: On the other hand, as the wave-speed estimates from Eq. (13) do not admit an isolated gridaligned shock, the fluxes were linearized about the one-dimensional numerical shock solution instead, for a truly meaningful analysis. The one-dimensional problem was initialized with an M = 7 shock using Eq. (21) on a uniform 11-cell grid. Ghost cells were used to prescribe the boundary conditions. The ghost cell on the left was fixed at the supersonic state while the pressure and density were updated at the ghost cell on the right using the Dirichlet and Neumann conditions, respectively. The horizontal velocity at the right ghost cell was updated to enforce a mass flow rate of unity. The problem was computed using the HLLC scheme until the density residual dropped below 10 −15 . Then, the one-dimensional solution was mapped to the two-dimensional 11 × 11 uniform Cartesian grid for the stability analysis. Since the HLLC, HLLCM, and HLLCT schemes are identical in one-dimension, the same numerical shock solution was used for all three schemes. The conserved variables Q were perturbed by an amplitude of ±10 −6 to approximate the flux gradients. The eigenvalues λ of the flux Jacobian matrices are plotted in Fig. 7. A flux scheme is deemed unstable if max(Re(λ)) 0, i.e., if the maximum real part of any of the eigenvalues is positive.
It can be noticed from the results that all three schemes are unstable when using the wavespeed estimates from Eq. (14) as evident from the positive values of max(Re(λ)). In comparison, using wave-speed estimates from Eq. (13) reduces max(Re(λ)) significantly indicating an improvement in shock stability. Yet, max(Re(λ)) still remains positive for the HLLC scheme. On the other hand, with max(Re(λ)) < 0, both the HLLCM and HLLCT schemes are stable when using wave-speed estimates from Eq. (13).

Numerical experiments
In this section, the HLLCT scheme is compared with the HLLC scheme, the hybrid HLLC-HLLCM scheme with the pressure-based shock sensor (Eq. (27)) and, where shear layers are concerned, with the HLLCM scheme as well. The comparison is performed only for twodimensional test cases since shock instability is multi-dimensional in nature. Moreover, the behaviors of the HLLC, HLLCM, and HLLCT schemes are identical in one-dimension.

Contact discontinuity
This problem is a test of a scheme's ability to preserve a grid-aligned contact discontinuity. The domain x × y ∈ [0, 1] × [0, 1] was discretized into a uniform 50 × 50 grid. The problem was initialized as follows: (p, u x ) = (1, 0), (ρ, u y ) = (10, 1), x < 0.5, where u x and u y refer to the horizontal and vertical velocity components, respectively. The Neumann boundary conditions were applied on all sides. The simulations were run for 100 iterations using the implicit solver. ρ and u y profiles along the horizontal line y = 0.5 are shown in Fig. 8. It can be easily shown from Eqs. (19) and (20) that the shear sensor w = 1 everywhere for this problem. Therefore, the HLLCT scheme reverts to the HLLC scheme and preserves the isolated contact discontinuity exactly. The HLLCM scheme, on the other hand, diffuses it. Due to the absence of any pressure difference, the HLLC-HLLCM scheme performed identically to the HLLC scheme as well and, hence, it is not included in Fig. 8.
The left boundary was fixed at the initial values while the Neumann boundary condition was applied to the remaining boundaries. The simulations were run for 1 500 iterations which was sufficient for the density residuals to drop to 10 −12 . The final results are plotted in Fig. 9 with sixteen equispaced velocity magnitude contours from u = 2.45 to u = 3.95. The minimum value of the shear sensor w (Eq. (19)) over the faces of each cell is plotted in Fig. 9(c) next to the velocity magnitude contours for the HLLCT scheme. It is clear that the HLLCT scheme performs on par with the HLLC scheme in capturing the shear layer while the HLLCM scheme diffuses it. It is immediately obvious that the velocitybased shear sensor identifies the regions of shear accurately and becomes approximately equal to one. ρ and u profiles sampled along the vertical line x = 1 are plotted in Fig. 10 for a more quantitative comparison. This reveals that the performance of the HLLCT scheme is indistinguishable from that of the HLLC scheme. It is interesting to note that the performance of the HLLCT scheme is not affected adversely in the regions where w ≈ 0. This is most likely due to the slow variations in the tangential velocities in those regions.  Fig. 11. The grid has been refined in the vertical direction close to the wall to include approximately 30 cells within the boundary layer region. A stagnation pressure of p 0 = 10 5 Pa and stagnation temperature of T 0 = 300 K was specified at the inlet boundary while a static pressure of p = 97 250 Pa was specified at the outlet boundary. The flat plate was modeled as an adiabatic, no-slip wall. The dynamic viscosity was taken to be µ = 1.87 × 10 −5 N·s·m −2 , and the specific heat capacity at constant pressure was taken to be c p = 1 004.5 J·kg −1 ·K −1 . Spatial discretization was performed using the second-order upwind scheme for the inviscid fluxes and the second-order central scheme for the viscous fluxes. The problem was computed using the HLLC, HLLCM, and HLLCT schemes until the density residuals dropped to below 10 −6 . The normalized horizontal velocities u/u ∞ at the exit plane at x = 0.304 8 m are compared with the Blasius solution in Fig. 12. It can be seen that, while the HLLCM scheme produces a diffused boundary layer, the HLLCT scheme resolves the boundary layer as accurately as the HLLC scheme which matches well with the analytical solution. This problem involves a two-dimensional steady isolated standing shock. Following the work of Dumbser et al. [14] , the problem was computed on the domain x × y ∈ [0, 1] × [0, 1] discretized uniformly into a 25 × 25 grid. The problem was initialized using Eq. (21) for an M = 20 shock that is initially located at x = 0.5. Perturbations of the relative order of 10 −6 were introduced to the initial conditions. The left and right boundary conditions were the same as those described in Subsection 4.3. The top and bottom boundaries were set to be periodic. The problem was computed using the implicit solver for 10 000 iterations. Fourteen equispaced density contours from ρ = 1.2 to ρ = 6.4 are plotted in Fig. 13.
The results demonstrate that the HLLCT scheme is able to remain stable for very strong shocks. Notice that the shock position at x = 0.5 is captured correctly by the HLLCT scheme. Observing the plot of min(w) in Fig. 13(b), it can be seen that the shear dissipation is the strongest not at the shock itself but slightly downstream of it around x ≈ 0.65. It is likely that  the numerical perturbations cause the value of ε u to become large enough about this location for the shear dissipation to be activated.

Quirk's test: odd-even decoupling problem
This classic test proposed by Quirk [5] involves an M = 6 shock propagating along a perturbed grid. The domain x × y ∈ [0, 800] × [−20, 20] was discretized into a uniform 800 × 40 mesh. The horizontal line y = 0 was then perturbed in a sawtooth fashion as follows: The problem was initialized with the shock positioned at x = 5 as follows: The transient simulation was performed until t = 105 s at which point the shock front would have propagated to x = 750. The final results are shown in Fig. 14 using twenty-five equispaced density contours from ρ = 1 to ρ = 5.8. The pressure-based shock sensor (Eq. (16)) and velocity-based shear sensor are plotted below the density contours of the HLLC-HLLCM and HLLCT schemes, respectively.
While the HLLC scheme breaks down the shock front and produces a prominent carbuncle, both the HLLC-HLLCM and HLLCT schemes keep it intact. It can be observed from Fig. 14(c) that shear dissipation is the most active around y = 0 where the grid is perturbed. It is interesting to note that the region of increased dissipation around y = 0 resembles the shock instability pattern observed for the HLLC scheme in Fig. 14(a).
The computation was repeated for an M = 20 shock to demonstrate that the choice of numerical constants in Eq. (20) applies for a wide range of shock strengths. The computation was performed on the extended domain x × y ∈ [0, 1 500] × [−20, 20] which was discretized into a uniform 1 500 × 40 mesh. The horizontal line y = 0 was perturbed as before. The problem was initialized with the shock positioned at x = 5 as follows: The problem was computed until t = 61 s at which point the shock would have propagated to approximately x = 1 448.5. The final results are shown in Fig. 15 using twenty-six equispaced density contours from ρ = 1 to ρ = 6. Once again, the HLLC scheme breaks down the shock front completely. Interestingly, the HLLC-HLLCM hybrid scheme also becomes susceptible to the shock instability due to insufficient dissipation behind the shock. It can be seen from Fig. 15(b) that the pressure-based sensor does not encompass the entire shock. On the other hand, The HLLCT scheme remains stable and captures the shock front accurately. This test demonstrates that the choice of the numerical constants used for ε u in Eq. (20) is applicable for a wide range of shock strengths.

Hypersonic flow past blunt body
This problem involves an M ∞ = 20 flow past a cylinder. The gridpoints (x i,j , y i,j ) were obtained from the expressions where θ j = 11π 15 (j/J − 0.5), x ci = 1.8(I − i)/I, and r i = 1 + 2.4(I − i)/I. The indices i and j run from 0 to I and 0 to J, respectively. Hence, the grid consists of I × J cells in total. A 10 × 200 grid is shown in Fig. 16(a). The problem was computed on a 20 × 400 grid to trigger a prominent carbuncle. The problem was initialized as follows: (ρ, u x , u y , p) = (1.4, 20, 0, 1).
The initial conditions were applied at the inlet. The reflective boundary condition was prescribed along the cylinder, and the Neumann boundary conditions were applied at the top and bottom boundaries. The simulations were run until the density residuals dropped down to 10 −11 . The final results are shown in Fig. 16 using fifteen equispaced density contours from ρ = 1.4 to ρ = 8.4. The pressure-based shock sensor and velocity-based shear sensor are plotted next to the density contours of the HLLC-HLLCM and HLLCT schemes, respectively.
The HLLC scheme produces a very prominent carbuncle, but the HLLCT scheme is carbunclefree similar to the HLLC-HLLCM scheme. It is interesting to note that min(w) remains at about 0.75 near the shock for the HLLCT scheme and yet, this is sufficient to suppress the shock instability in this problem.

Double Mach reflection problem
The double Mach reflection (DMR) problem proposed by Woodward and Colella [26] involves an M = 10 shock impinging on an inclined ramp to form a system of complex shock reflections. The DMR problem was solved on a modified computational domain proposed by Vevek et al. [28] to prevent the formation of numerical artefacts. A coarse mesh of the modified domain is shown in Fig. 17. The were prescribed at the right and left boundaries, respectively. The shock was initially positioned at x = −0.9 and allowed to propagate for t = 0.08 s at which time it would approximately impinge on the ramp. Next, the region to the left shock is re-initialized with the exact postshock conditions given above to eliminate the numerical noise generated as the shock is smeared on the grid. Then, the simulation is continued for another 0.19 s until t = 0.27 s. The final results are plotted in Fig. 18 using forty-three equispaced density contours from ρ = 1.4 to ρ = 22.4. The pressure-based shock sensor and velocity-based shear sensor are plotted below the density contours of the HLLC-HLLCM and HLLCT schemes, respectively. The HLLC scheme produces a kinked Mach stem near x = 2.5. While the HLLC-HLLCM scheme reduces the extent of the kink, it does not eliminate it completely. On the other hand, the HLLCT scheme removes the kink entirely. Secondly, it is clear from the plot of the shear sensor in Fig. 18(c) that the value of w is approximately one near the slip lines. Consequently, the HLLCT scheme behaves like the HLLC scheme and captures the slip lines with little dissipation. Comparing Figs. 18(b) and 18(d), it can be seen that the HLLCM scheme has destroyed the roll-up of the primary slip line near the bottom wall while the HLLCT scheme has improved the resolution of the slip line significantly.
The left boundary was set to supersonic inflow using the above initial conditions while the Neumann boundary condition was applied at the right boundary. All other boundaries were set to be reflective. The special treatment suggested by Woodward and Colella [26] was applied at the step corner to minimize the effects of the unphysical 'boundary layer'. The problem was computed till t = 3 s, and the final results are plotted in Fig. 19, using the forty-one velocity magnitude contours from u = 0 to u = 4 to show the slip line emanating from the triple point clearly. The pressure-based shock sensor and velocity-based shear sensor are plotted below the velocity magnitude contours of the HLLC-HLLCM and HLLCT schemes, respectively.
The HLLC scheme results in instabilities downstream of the normal shock which are eliminated by the other three schemes. While the HLLCM scheme diffuses the slip lines issuing from the triple point almost completely, both the HLLC-HLLCM and HLLCT schemes are able to capture them well. Observing the plot of min(w) in Fig. 19(d), it is evident that the shear sensor identifies the slip lines accurately and attains a value close to unity, thus reducing the dissipation across the slip lines. For a more quantitative comparison, ρ and u profiles sampled along the vertical line x = 1.25 are plotted in Fig. 20. The results near the slip lines are plotted in the insets. Interestingly, the HLLCT scheme performs slightly better than the HLLC-HLLCM scheme as benchmarked against the HLLC scheme. This is mostly likely because, unlike the pressure-based shock sensor, the velocity-based shear sensor does not vanish completely around the shock. Therefore, the HLLCT scheme introduces lesser dissipation than the HLLC-HLLCM scheme. Nevertheless, this is sufficient in eliminating the shock instability.

Computational efficiency
It has been demonstrated thus far that the HLLCT scheme is able to resolve shear layers well without succumbing to the shock instability problem. However, it is also important to ensure that the scheme is computationally efficient. In Table 1, a comparison of the computation time of the different schemes is provided for the two largest problems considered in this study, namely, the DMR and forward-facing step problems. The computation time has been normalized with respect to that of the HLLC scheme. It is clear that the HLLCT scheme remains competitive with the original HLLC scheme. Though the HLLC-HLLCM scheme takes the longest time amongst the schemes, it is not far behind the others. This shows that, if implemented properly, shock sensor computation is not necessarily a demanding operation. Its main drawback lies in introducing it within existing codes which have been designed to compute fluxes in a localized manner. In this regard, it is straightforward to incorporate the HLLCT scheme in such codes.

Conclusions
In conventional carbuncle cures, the HLLC scheme is coupled with a more diffusive, shockstable scheme, and the complementary scheme is applied in the vicinity of shocks. This requires a suitable shock sensor which needs to be computed in a non-local manner, considering all the neighboring faces to assure the elimination of the shock instability. The contact-preserving shock-stable HLLCM scheme proposed by Shen et al. [10] presents a unique opportunity to explore the viability of a different approach to curing the carbuncle problem. Since the HLLCM scheme introduces spurious dissipation only across shear layers, it was hypothesized that it might be sufficient to restrict the action of the HLLC scheme to regions of shear and use the HLLCM scheme everywhere else. To evaluate this hypothesis, a hybrid HLLC-HLLCM scheme called the HLLCT scheme was designed. The crucial advantage of the HLLCT scheme comes from a velocity-based shear sensor which can be computed in a localized manner based on the left and right states only. Numerical experiments on grid-aligned and non-grid-aligned shear layers show that the HLLCT scheme is able to capture shear layers with as little dissipation as the HLLC scheme itself, attesting to the effectiveness of the shear sensor. Interestingly, the use of the non-shear-preserving HLLCM scheme in regions away from shear layers does not seem to have an adverse effect on the results. More importantly, the HLLCT scheme is also found to be free from any shock instabilities for a number of classic test cases involving strong grid-aligned shocks for which the HLLC scheme fails.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.