Strategies to cure numerical shock instability in HLLEM Riemann solver

Strategies to cure numerical shock instability in HLLEM Riemann solver

Sangeeth Simon J. C. Mandal Corresponding author: mandal@iitb.ac.in
Abstract

The HLLEM scheme is a popular contact and shear preserving approximate Riemann solver for cheap and accurate computation of high speed gasdynamical flows. Unfortunately this scheme is known to be plagued by various forms of numerical shock instability. In this paper we present various strategies to save the HLLEM scheme from developing such spurious solutions. A linear scale analysis of its mass and interface-normal momentum flux discretizations reveal that its antidiffusive terms, which are primarily responsible for resolution of linear wavefields, are inadvertently activated along a normal shock front due to numerical perturbations. These erroneously activated terms counteract the favourable damping mechanism provided by its inherent HLL-type diffusive terms and trigger the shock instability. To avoid this, two different strategies are proposed for discretization of these critical flux components in the vicinity of a shock: one that deals with increasing the magnitude of inherent HLL-type dissipation through careful manipulation of specific non-linear wave speed estimates while the other deals with reducing the magnitude of these critical antidiffusive terms. A linear perturbation analysis is performed to gauge the effectiveness of these cures and estimate von-Neumann type stability bounds on the CFL number arising from their use. Results from classic numerical test cases show that both types of modified HLLEM schemes are able to provide excellent shock stable solutions while retaining commendable accuracy on shear dominated viscous flows.

1 Introduction

Computation of high speed gasdynamical flows has benefited tremendously from the development of various approximate Riemann solvers for the Euler system of equations. Although most of these schemes are designed to capture nonlinear waves like shocks and expansion fans, they can be broadly classifed into two groups based on their ability to capture the linear waves like entropy and shear waves accurately. The resulting two classes of approximate Riemann solvers are generally referred to in literature as contact-shear preserving and contact-shear dissipative solvers respectively. Roe scheme [1] is one of the most popular examples of the former category and is widely incorporated in most industry and research codes. However, Roe scheme has several problems like requirement of an entropy fix to resolve sonic expansions, lack of positivity in density and internal energy in near vaccum flows and requirement of the knowledge of full eigen structure of the flux jacobians. On the other hand a popular example of a contact-shear dissipative type is the HLL scheme [2, 3]. Since these are based directly on integral form of conservation laws, they do not require linearizations or complete knowledge of eigenstructures of the Euler system. Consequently, these methods do not suffer from the issues mentioned here that plagues the Roe scheme. However lack of contact and shear ability seriously restricts their applicability to practical problems involving shearing flows, boundary layer flows, flame and material fronts etc [4, 5]. Einfeldt et al [6] proposed an interesting option that combines the advantages of the Roe scheme with that of the HLL scheme called HLLEM scheme. The HLLEM scheme is positivity preserving, contact-shear capturing and does not require an entropy fix. Also it demands knowledge of only the eigenstructure of the linear wavefields of the system.

Unfortunately like other contact-shear preserving Riemann solvers, the HLLEM scheme also suffers from various forms of numerical shock instability during simulations of multidimensional flows which involves regions of strong grid aligned normal shocks. A classic example of such a failure is the Carbuncle phenomenon [7]. A catalogue of such failures can be found in [8]. 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 after shock values, growth in error norms in case of steady state problems and in extreme cases complete breakdown of the solution.

Although there is still no agreement on the exact cause of these spurious solutions, many authors suggest that lack of cross-flow dissipation provided by the contact-shear preserving Riemann solvers along the shock front could be a major trigger [9, 10]. A loss in overall dissipation for a contact-shear preserving scheme in such circumstances are argued to happen due to vanishing of their numerical dissipation terms that corresponds to accuracy on linear waves [10, 11]. On the other hand, a linear wave dissipative scheme like the HLLE does not suffer from this. Many cures for creating a shock stable Roe scheme, for example, are hence based on increasing the component of dissipation corresponding to the linear waves in the overall dissipation matrix [9, 11, 12, 13, 14, 15, 16]. Similarly, in case of the HLLEM scheme, Park et al [17] proposed a control of the antidiffusive terms, which are chiefly responsible for restoring accuracy on linear wavefields, to ensure shock stability. Xie et al [18] showed that merely controlling the antidiffusive terms corresponding to the shear wave could guarantee shock stability. Obayashi et al [19] proposed a Carbuncle free HLLEM scheme based on retaining only its HLLE component in the vicinity of a shock although their scheme could benefit from more rigorous testing.

Recently it was clarified in case of another popular HLL-based contact-shear preserving Riemann solver called HLLC scheme, that the antidiffusive terms, specifically in mass and interface-normal momentum flux discretizations in transverse direction of a shock front are the major triggers for instability [20]. It has been found that the antidiffusive terms present in the HLLC discretization of these critical flux components gets activated in the transverse direction of a numerically perturbed normal shock front and causes a reduction in the magnitude of its inherent HLL-type diffusive terms. The weakened diffusive terms fail to provide requisite damping of perturbations in and variables resulting in unphysical variation of conserved quantity along the shock front and subsequent development of shock unstable solutions. Based on this, the authors suggested an inexpensive shock stable HLLC scheme by simply controlling the antidiffusive terms in these critical flux components using a simple differentiable pressure sensor. An interesting alternative to save the HLLC scheme that does not rely upon controlling the antidiffusive terms was proposed in [21]. It involves a strategy to increase the inherent HLL-type dissipation of this scheme in the vicinity of shock through careful modification of certain nonlinear wavespeed estimates. Since the HLLC scheme and the HLLEM scheme are both derived from the HLL scheme, they share identical HLL-type diffusive terms. Thus, it opens up the possibility of extending these appealing cures developed for the former to the latter provided they share similar instability-trigger mechanism.

In this paper, we first perform a linear scale analysis of the numerical dissipation terms of the HLLEM scheme appearing in the mass and interface-normal momentum flux discretizations in the vicinity of a steady shock front subjected to numerical perturbations with an objective to identify specific terms in its formulation that triggers the instability. Our analysis indicates that the instability in the HLLEM scheme is also promoted due to a weakening of its embedded HLL-type dissipation in the direction parallel to a shock front due to erroneous activation of its antidiffusive terms. Based on these findings, we propose two robust strategies to save the HLLEM scheme from shock instability that are based on methods suggested in [20, 21]. The first method focuses on increasing the inherent HLL-type numerical dissipation existing within the HLLEM scheme. This is achieved by carefully modifying non linear wave speed estimates of this embedded HLLE dissipation. Two multidimensional local solution dependent strategies are discussed to ascertain the quantity of supplementary dissipation to be provided to enhance this HLLE component. The second method focuses on reducing the magnitude of the antidiffusive component of the HLLEM scheme relative to its inherent HLL-type dissipation. A technique of linear perturbation analysis [13] is used to show how these different strategies are able to easily dampen numerical perturbations that may appear in the primitive quantities. Corresponding von Neumann stability bounds on the allowable CFL number for implementation of all of these strategies are also estimated. We also demonstrate through our analysis and numerical examples how simply controlling the antidiffusive terms corresponding to shear waves alone, as suggested in [18], may not suffice to create a shock stable HLLEM scheme. Finally, a suite of first order and second order test cases are used to study the effectiveness of these methods in dealing with shock instabilities.

This paper is organized as follows. In Sec.(2,3) we brief the governing equations and finite volume framework adopted in this paper. Sec.(4) recapitulates the details of the HLLE and the HLLEM schemes that form the basis of our discussion. In Sec.(5) we use a linear scale analysis on the HLLE and HLLEM schemes to identify the specific terms in the HLLEM scheme that triggers the instability. This analysis is supplemented with a linear perturbation analysis that provides insights on the evolution of perturbations of primitive quantities in these schemes. Based on these, we describe strategies to formulate a few shock stable versions of the HLLEM scheme in Sec.(6). The linear perturbation analysis developed earlier is extended to study the behavior of these new shock stable versions in dealing with the perturbations. In Sec.(7) certain classic shock instability test cases are used to demonstrate the robustness of these new formulations. Sec(8) presents some concluding remarks.

2 Governing equations

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

(0)

where are the vector of conserved variables and x and y directional fluxes respectively. These are given by,

(0)

In the above expressions, and stands for density, x-velocity, y-velocity, pressure and specific total energy. The system of equations are closed through the equation of state,

(0)

where is the ratio of specific heats. Present work assumes a calorifically perfect gas with . A particularly useful form of the Eq.(2) is the integral form given by,

(0)

where denotes a control volume over which Eq.(2) describes a Finite Volume balance of the conserved quantities, and denotes the x and y dimensions of the control volume respectively, denotes the boundary surface of the control volume and denotes an infinitesimally small element on . is the outward pointing unit normal vector to the surface .

3 Finite volume discretization

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

Hence i

Figure 1: Typical control volume with its associated interfaces and respective normal vectors .
(0)

where is an appropriate cell averaged conserved state vector, denotes the flux vector at the mid point of each interface while and denotes the unit normal vector and the length of each interface respectively. These are shown in Fig.(1). The interface flux 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 as,

(0)

where indicates the initial conditions of a local Riemann problem across interface shared by the cells and . The matrices and are rotation matrices at the interface given by,

(0)

where denote the components of the normal vector .

In the next section we briefly describe two approximate Riemann solvers named the HLLE scheme and the HLLEM scheme that can be used to estimate the flux at a given interface. To avoid extra notations, henceforth simply represent the local Riemann flux at any interface with outward pointing normal .

4 Recap of HLLE and HLLEM schemes

The original HLL scheme was devised by Harten, Lax and van Leer [2]. Unlike linearized solver of Roe, the HLL scheme does not require the knowledge of the eigenstructure of the Euler system. Instead it assumes a wave structure consisting of two waves that seperates three constant states. Using the integral form of the conservation laws on this wave structure, the HLL Riemann flux can then be written as,

(0)

Here and are the local fluxes at either side of the interface. and are numerical approximations to the speeds of the left most and right most running characteristics that emerge as the solution of the Riemann problem at an interface. It has been shown that under appropriate choice of wavespeeds and , the HLL scheme is both positivity preserving and entropy satisfying [3]. This choice of wavespeeds are given as,

(0)

where are the normal velocity across an interface, are the respective sonic speeds and are the standard Roe averaged quantities at the interface [1]. Using these wavespeed estimates, the HLL scheme is also known as the HLLE scheme [3]. For the purpose of this paper, we note that HLLE scheme can also be rewritten as,

(0)

where the dissipation can be written as [22],

(0)

Although quite accurate in resolution of shocks and expansion fans, the major drawback of HLL scheme is its inability to resolve the contact and shear waves. The loss of accuracy on these waves occur because of the assumption of constant average state between the two waves. Further, it is observed that HLL Riemann solver is free from various forms of numerical shock instability which is generally attributed to its linear wave dissipative nature [13, 23].

Einfeldt et al [6] presented an improved version of the HLLE scheme that can capture the contact and shear waves in addition to retaining the excellent positivity preserving and entropy satisfying property of the original HLLE scheme. Additional accuracy on these waves for HLLE scheme were obtained by explicitly adding necessary antidiffusive terms corresponding to only linear waves. The HLLEM flux can be written as,

(0)

with wave speeds still defined by Eq.(4). Here, terms denote respectively the right eigenvectors of the flux jacobians corresponding to the linearly degenerate waves and the strength of these waves obtained by projecting the total jump onto these eigenvectors. Both these quantities are evaluated at the Roe averaged state. These are given as [24],

(0)
(0)

where operator represents the operation , denotes the tangential component of velocity at an interface and denotes the corresponding Roe averaged quantity. The terms represent coefficients that control the amount of antidissipation being introduced into these waves. An estimate for these that guarantees good accuracy for resolving linear waves are given as [17],

(0)

As in the HLLE scheme, we rewrite the HLLEM scheme as,

(0)

where the dissipation vector is same as in Eq.(4) and represents the embedded HLL-type dissipation while the antidissipation vector unique to the HLLEM scheme is given as,

(0)

In the above expression, the first term denote antidissipation concerning the contact wave while the second term denote that concerning the shear wave. As will be shown in the next section, these terms responsible for accuracy on linear waves will also be responsible for lack of robustness of the HLLEM scheme towards various forms of numerical shock instability.

5 Instability characteristics of HLLE and HLLEM schemes

Studies [10, 20] show that unphysical variation of conserved quantity along a numerically computed shock front could be a plausible reason for manifestation of shock instability. Such unphysical variation could occur independently due to undamped perturbations in density () or flow velocity () variables [20]. In case of other shock unstable approximate Riemann solvers like the Roe scheme and the HLLC scheme, the proliferation of these perturbations are known to occur typically due to reduction in overall numerical dissipation caused by inadvertent activation of their antidiffusive terms on interfaces that are not aligned with the shock front [10, 20]. Specifically it has been shown in [20] that activation of antidiffusive terms in the mass and interface-normal flux component discretizations on these crucial transverse interfaces play a major role in triggering of the instability. In this context we would like to closely study the behaviour of the numerical dissipation of the shock unstable HLLEM scheme in comparison to that of the shock stable HLLE scheme to identify the specific terms of the former that participate in the shock instability mechanism. The study is carried out in a simple setting most commonly known to produce shock instability: in the vicinity of a perturbed strong normal shock. For this purpose, consider a y-directional stencil comprising of three candidate cells namely and located within the upstream coloum of cells of an isolated strong normal shock that exists in a steady supersonic flow as shown in Fig.(2).

Figure 2: Schematic showing the stencil chosen for studying the numerical dissipation characterisitcs of the HLLE and the HLLEM schemes.

The flow is assumed to happen only in positive x direction. Along the transverse direction (ie. along the y-directional interfaces in this case) in the vicinity of the shock, the following assumptions can be made without loss of generality:

(0)

Our particular concern is the evolution of the conserved quantities and in cell () due to fluxes that cross the interfaces (). While the mass fluxes affect the evolution of quantity , the interface-normal fluxes affect the evolution of quantity . The saw-tooth like perturbations in primitive and conserved variables that are characteristic of a shock unstable solution can be thought to occur due to imbalances in these fluxes only [8]. An evolution equation for these quantities can be written as,

(0)
(0)

Notice that since there is no primary flow in the y-direction, subsonic fluxes corresponding to the HLLE or HLLEM scheme will be engaged on these interfaces. Ideally, across these interfaces only information regarding pressure would be transmited through the interface-normal momentum fluxes. However if numerical perturbations in flow quantities does exist, then additional unphysical fluxes would occur through these interfaces. The objective then is to see how the HLLEM scheme differ from the HLLE scheme in its treatment of these unphysical fluxes and the perturbed flow quantities that gets convected because of them. In case of the HLLEM scheme, we pay special attention to the role of its antidiffusive terms on the evolution of and quantities.

5.1 Dissipation analysis of the HLLE scheme

Suppose consider that we employ the HLLE scheme to evaluate the fluxes in Eqs.(5) and (5). Then, we are interested in how the numerical dissipation available in the HLLE discretization of these critical flux components affect the evolution of and respectively. Below, we seperately analyze each of them.

5.1.1 Mass flux

A HLLE discretization for and can be written as,

where we have dropped the from the interface index for convenience. The flux difference is then,

Since our interest lies in contrasting the numerical dissipation behaviours of the HLLE and the HLLEM schemes, it suffices to only consider the dissipation component of the flux difference which will be henceforth denoted as . Using the notation , the dissipation component of the above flux difference can be written as,

(0)

To introduce the effect of small numerical perturbations that are thought to eventually result in shock unstable solutions, consider the existence of a random numerical perturbations of the order (eg. due to round of errors) in the flow variables that exists during the course of computation in the stencil considered. Then, the following approximations are assumed to hold true on the stencil,

(0)

The term can be expanded as . Under the above assumptions is . Assuming and , using Eq.(5) and Eq.(5.1.1), Eq.(5.1.1) can be simplified as,

(0)

Effectively,

(0)

Thus in the presence of perturbations of in flow quantities, the HLLE scheme infuses a net dissipation of the same order into the mass flux discretization.

5.1.2 Interface-normal momentum flux

A HLLE discretization for and can be written as,

The dissipation component of the flux difference in this case denoted as will be,

(0)

The term can be expanded in terms of the mass flux as and under the assumptions made in Eq.(5.1.1) is . Using Eq.(5) and Eq.(5.1.1), Eq.(5.1.2) can be simplified as,

(0)

Which indicates,

(0)

Thus in the presence of perturbations of in flow quantities, the HLLE scheme infuses a net dissipation of the same order into the discretization of interface-normal momentum flux in the y-direction.

5.1.3 Role of the activated diffusive terms in the HLLE scheme

It was seen in Sections.(5.1.1) and (5.1.2) that a numerical perturbation in the flow quantities, in the vicinity of a normal shock, tends to activate the dissipation component of the HLLE scheme in both mass and interface-normal momentum flux components. To confirm whether this amount of dissipation suffices to damp out the perturbations in flow variables like , and which affect the quantities and , we resort to a technique of linear perturbation analysis. Mainly, we seek to develop the linearized evolution equations for perturbation in these variables from the approximate evolution equations given in Eqs.(5) and (5) as per the technique provided in [13]. For developing the evolution equation for we additionally utilize an approximate evolution equation of the total energy in transverse direction given by,

(0)

In the present work, the steady mean flow is chosen to have normalised state values of , , and . The perturbations in each of these variables denoted as and are introduced as a saw tooth profile superimposed on the steady mean flow. Thus for a typical cell ’j’, the perturbed initial conditions are,

(0)

Evolution equations for perturbations in primitive variables , and in case of the HLLE scheme are given in Table 1 where denotes a linearized CFL value. The amplification factors for density, x-velocity and pressure in case of HLLE scheme are respectively. Thus for , any initial perturbation in these flow variables would be effectively damped by the HLLE dissipation. The evolution of these perturbations are plotted in Fig.(3) and is more illustrative of the above discussion. Each plot in Fig.(3) indicate the evolution of all three perturbations: and for a given initial perturbation in only one of those quantities (while setting the other perturbations to be 0). For all the experiments is taken to be 0.2.

Scheme Density Pressure x-velocity
HLLE
HLLEM
HLLEM-SWM-P/E
HLLEM-ADC
HLLEMS
Table 1: Evolution equations for perturbations in primitive variables for various schemes
(a)
(b)
(c)
Figure 3: Comparison of evolution of density, x-velocity and pressure perturbations in the HLLE scheme.

From the behaviour of these evolution equations, it is clear that any perturbation in all three flow variables are attenuated in time by the HLLE scheme. This indicates that the activated dissipative terms of the HLLE scheme as found in Eqs.(5.1.1) and (5.1.2) are capable of damping any unwanted perturbations that may exist in the course of a numerical computation. This in effect ensures that variation of the conserved quantity along the shock front remains spatially and temporally unvarying and the shock structure is preserved through an accurate Rankine-Hugoniot jump across it.

5.2 Dissipation analysis of the HLLEM scheme

Consider that instead of the HLLE scheme, we use the HLLEM scheme to evaluate the mass and interface-normal momentum fluxes in in Eqs.(5) and (5) at these crucial transverse interfaces. Then, we are interested in how the numerical dissipation provided by the HLLEM scheme differ from that already seen in the HLLE scheme. Below, we seperately analyze each of them. Since HLLEM scheme can be thought of as a combination of the HLLE scheme and antidiffusive terms that are responsible for accuracy on linear waves, we specially study the role of these antidiffusive terms in dealing with these perturbations.

5.2.1 Mass flux

A HLLEM discretization for and can be written as,

(0)
(0)

The flux difference is then,

(0)

Since we are interested only in the numerical dissipation provided by the HLLEM scheme, we seek to isolate these terms from the above equation. These terms can be extracted as,

(0)

Notice that we have expressed the numerical dissipation of the HLLEM scheme as a combination of HLL-type dissipation (denoted by ) and dissipation arising from its antidiffusive terms. An order of magnitude estimate for has been already presented in Eq.(5.1.1). We refer to the second term that arises from antidiffusive component of the HLLEM scheme as . Now we seek an estimate for this term. Firstly notice that using the conditions presented in Eq.(5)and Eq.(5.1.1), the wave strength and antidissipation coefficient at each interface can be simplified as,

(0)

Using these the term can be estimated to be,

(0)

Eq.(5.2.1) reveals something interesting. In the presence of numerical perturbations of , it is observed that the antidiffusive terms in the HLLEM massflux discretization is activated and is also of similar to its inherent HLL-type diffusive terms.

5.2.2 Interface-normal momentum flux

A HLLEM discretization for and can be written as,

(0)

and,

(0)

Flux difference in y-direction is then,

The terms corresponding to numerical dissipation can be extracted as,

(0)

We have again expressed the numerical dissipation of the HLLEM scheme as a combination of a HLL-type dissipation (denoted by ) and dissipation arising from the antidiffusive terms. An order of magnitude estimate for has been already presented in Eq.(5.1.2). We refer to the second term that arises from the antidiffusive terms of the HLLEM scheme as . By using the conditions presented in Eq.(5)and Eq.(5.1.1), this term can be further simplified. Note that under the prescribed conditions, the wave strength . Hence it can be shown that,

(0)

Once again, it is observed that numerical perturbations of have activated the antidiffusive terms in the HLLEM interface-normal momentum flux discretization in y-direction. Noticeably, this term is also of similar to its inherent HLL-type diffusive terms.

5.2.3 Role of the activated antidiffusive terms in the HLLEM scheme

To see the effect of these antidiffusive terms more clearly, we develop the corresponding linearized evolution equations for perturbations in , and for the HLLEM scheme similar to the ones described for the HLLE scheme in section(5.1.3). These are given in Table 1. It is seen from these equations that in case of the HLLEM scheme, the amplification factors for perturbations in , and are respectively. This indicates that the HLLEM scheme is incapable of damping out perturbations in and variables. Additionally it is seen that although there exist a mechanism to damp out perturbations in for , any finite is repeatedly fed into which results in the existence of a residual that is growing with time. The behaviour of these evolution equations are given in Fig.(4). To specifically understand the role of antidiffusive terms in triggering shock instability, we rewrite the evolution equations for perturbations in and variables as,