An all speed second order IMEX relaxation scheme for the Euler equations

We present an implicit-explicit finite volume scheme for the Euler equations. We start from the non-dimensionalised Euler equations where we split the pressure in a slow and a fast acoustic part. We use a Suliciu type relaxation model which we split in an explicit part, solved using a Godunov-type scheme based on an approximate Riemann solver, and an implicit part where we solve an elliptic equation for the fast pressure. The relaxation source terms are treated projecting the solution on the equilibrium manifold. The proposed scheme is positivity preserving with respect to the density and internal energy and asymptotic preserving towards the incompressible Euler equations. For this first order scheme we give a second order extension which maintains the positivity property. We perform numerical experiments in 1D and 2D to show the applicability of the proposed splitting and give convergence results for the second order extension.

Authors

• 2 publications
• 1 publication
• 3 publications
• 8 publications
• An all speed second order well-balanced IMEX relaxation scheme for the Euler equations with gravity

We present an implicit-explicit well-balanced finite volume scheme for t...
12/19/2019 ∙ by Andrea Thomann, et al. ∙ 0

• Asymptotic Preserving and Low Mach Number Accurate IMEX Finite Volume Schemes for the Euler Equations

In this paper the design and analysis of a class of second order accurat...
07/03/2019 ∙ by Dr. K. R. Arun, et al. ∙ 0

• A second-order well-balanced finite volume scheme for Euler equations with gravity

We present a well-balanced, second order, Godunov-type finite volume sch...
03/12/2018 ∙ by Deepak Varma, et al. ∙ 0

• On Fake Accuracy Verification

In this paper, we reveal a mechanism behind a fake accuracy verification...
08/27/2020 ∙ by Hiroaki Nishikawa, et al. ∙ 0

• Efficient parallel 3D computation of the compressible Euler equations with an invariant-domain preserving second-order finite-element scheme

We discuss the efficient implementation of a high-performance second-ord...
06/30/2020 ∙ by Matthias Maier, et al. ∙ 0

• Consistent Internal Energy Based Schemes for the Compressible Euler Equations

Numerical schemes for the solution of the Euler equations have recently ...
06/26/2019 ∙ by R. Herbin, et al. ∙ 0

• A time adaptive multirate Dirichlet-Neumann waveform relaxation method for heterogeneous coupled heat equations

We consider partitioned time integration for heterogeneous coupled heat ...
07/01/2020 ∙ by Peter Meisrimel, et al. ∙ 0

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

Abstract

We present an implicit-explicit finite volume scheme for the Euler equations. We start from the non-dimensionalised Euler equations where we split the pressure in a slow and a fast acoustic part. We use a Suliciu type relaxation model which we split in an explicit part, solved using a Godunov-type scheme based on an approximate Riemann solver, and an implicit part where we solve an elliptic equation for the fast pressure. The relaxation source terms are treated projecting the solution on the equilibrium manifold. The proposed scheme is positivity preserving with respect to the density and internal energy and asymptotic preserving towards the incompressible Euler equations. For this first order scheme we give a second order extension which maintains the positivity property. We perform numerical experiments in 1D and 2D to show the applicability of the proposed splitting and give convergence results for the second order extension.

Keywords

finite volume methods, Euler equations, positivity preserving, asymptotic preserving, relaxation, low Mach scheme, IMEX schemes

1 Introduction

We consider the non-dimensional Euler equations in -space dimensions which are given by the following set of equations [1]

 ρt+∇⋅(ρu)=0(ρu)t+∇⋅(ρu⊗u+1M2pI)=0Et+∇⋅(u(E+p))=0, (1.1)

where the total energy is given by

 E=ρe+12M2ρ|u|2 (1.2)

and denotes the internal energy. The density is denoted by ,

is the velocity vector and

is a given Mach number which controls the ratio between the velocity of the gas and the sound speed. Depending on the magnitude of the Mach number the characteristic nature of the flow changes. This makes the numerical simulation of these flows very challenging but also a very interesting research subject with a wide range of applications, for example in astrophysical stellar evolution or multiphase flows [2, 3]. For large Mach numbers, the flow is governed by compressible effects whereas in the low Mach limit the compressible equations converge towards the incompressible regime. This behaviour was studied for example in [4, 1, 5]. We refer to [6] for a study of the full Euler equations.

Standard schemes designed for compressible flows like the Roe scheme [7] or Godunov type schemes fail due to exessive diffusion when applied in the low Mach regime. A lot of work is dedicated to cure this defect, see for instance [8, 9, 1, 2]. Another way to ensure accurate solutions in the low Mach regime is the development of asymptotic preserving schemes which are consistent with its limit behaviour as tends to zero, see for example [10, 11, 12] and references therein.

Due to the hyperbolic nature of (1.1) the time step for an explicit scheme is restricted by a stability CFL condition that depends on the inverse of the fastest wave speed. In the case of (1.1) the acoustic wave speeds tend to infinity as tends to zero which leads to very small time steps to guarantee the stability of the explicit scheme. As a side effect all waves will be resolved by the explicit scheme although the fast waves are not necessary to capture the motion of the fluid as they carry a negligible amount of energy. Implicit methods on the other hand allow larger time steps but introduce diffusion on the slow wave which leads to a loss of accuracy. In addition at each iteration a non-linear often ill conditioned algebraic system has to be solved. Implicit-explicit (IMEX) methods try to overcome those disadvantages by treating the stiff parts implicitly and thus allow for a Mach number independent time step. Many of those schemes are based on a splitting of the pressure in the spirit of Klein [13] since the stiffness of the system is closely related with the pressure, see for example [10, 14, 13].

Another way to avoid solving non-linear implicit systems is using a relaxation approach. See [15, 16, 17, 18, 19] for references on relaxation. The idea of using relaxation is to linearise the equations which makes it easier to solve them implicitly as done in [20] through Jin Xin relaxation [15]. The linear degenerate structure of the resulting relaxation models allows to use accurate Godunov-type finite volume methods due to the knowledge about the Riemann solution also with implicit schemes as done in [21]. Here, we want to combine a Suliciu type relaxation with an IMEX scheme and the splitting of the pressure. The relaxation model we use is based on relaxing the slow and the fast pressure separately. The flux is then split such that only some relaxation variables are treated implicitly which leads to solving only one equation implicitly. This results in a computationally cheap scheme. We do not split inside the flux of the Euler equations but treat all terms explicitly. This leads to a conservative explicit part for which we use a Godunov-type approximate Riemann solver. The use of a Riemann solver allows us to prove with little effort the positivity preserving property of the scheme which is important in physical applications. In addition the way the flux is splitted leads to an asymptotic preserving scheme.

To have a relevant scheme for applications an extension to second or higher order is necessary. Higher order schemes in time can be achieved by using IMEX Runge Kutta (RK) methods as in [12, 22]. As standard in finite volume schemes higher order in space is achieved by reconstructing the cell interface values using WENO schemes for example [23, 24]. We use a MUSCL approach [25] to achieve a second order extension to the first order scheme which preserves the positivity of density and internal energy.

The paper is organized as follows. In the next section we describe the relaxation model that is used to derive the IMEX scheme. Section 3 is dedicated to the splitting of the flux into implicit and explicit terms and the structure of the first order time semi-discrete scheme. The asymptotic preserving property is proved in Section 4. Next, we give the derivation of the fully discrete scheme in Section 5. Therein the Godunov type scheme for the explicit part is given as well as the proof of the positivity of the resulting numerical scheme. In addition, we show that the diffusion introduced by the Riemann solver is independent of the Mach number due to the splitting described in Section 3. In Section 6 we give a second order extension for the first order scheme for which we show that it preserves the positivity property. It is followed by a section of numerical results to validate the theoretical results. A section of conclusion completes this paper.

2 Suliciu Relaxation model

To simplify the non-linear structure of the Euler equations (1.1) we make use of the Suliciu relaxation approach [16, 26, 27] and references therein. Compared to the Jin Xin relaxation [15] the original system of equations remains unchanged. Whereas in the Jin Xin relaxation approach the linearisation of the equations is achieved by relaxing every component of the flux function, in the Suliciu type relaxation single variables are relaxed in a fashion that is tailored to the problem. This leads to a reduced diffusion introduced by the relaxation process compared to the Jin Xin relaxation.

Following the usual Suliciu relaxation procedure, the key element is the relaxation of the pressure by introducing a new variable ,

 (ρπ)t+∇⋅(ρπu)+a2∇⋅u=ρϵ(p−π). (2.1)

Here denotes the relaxation time. Equation (2.1) is derived by multiplying the density equation by . The non-linearity that thereby arises is replaced by a constant parameter , in the following called the relaxation parameter. To guarantee a stable diffusive approximation of the original Euler equations the relaxation parameter must meet the sub-characteristic condition [17]. We want to profit from the properties of the Suliciu relaxation also for the non-dimensional Euler equations. In the following, we describe the relaxation model that was first introduced in [21]. Following [28], in a first step the pressure is decomposed into a slow dynamics and a fast acoustic component

 pM2=p+1−M2M2p.

Approximating the slow and the fast pressure in the momentum equation in (1.1) by the variables and respectively, momentum equation becomes

 (ρu)t+∇⋅(ρu⊗u+π+1−M2M2ψ)=0

The evolution of the new variables and are then developed in the spirit of a Suliciu relaxation approach. The evolution for is given by the Suliciu relaxation equation (2.1). However, applying the standard Suliciu relaxation method also on the pressure , would only lead to non relevant diffusion terms and not to a low Mach scheme. To overcome this, the authors of [21] introduced a new velocity variable which is relaxed to the system velocity and couples to the pressure . The form of the evolution equations for and is chosen, such that

• the resulting model is in conservation form

• the resulting model has ordered eigenvalues, which results in a clear wave structure

• the resulting model is a stable diffusive approximation of the non-dimensional Euler equations (1.1)

• the resulting numerical scheme has a Mach number independent diffusion.

Considering the above points, this leads to the following relaxation model

 ρt+∇⋅(ρu)=0,(ρu)t+∇⋅(ρu⊗u+π+1−M2M2ψ)=0,Et+∇⋅(u(E+M2π+(1−M2)ψ))=0,(ρπ)t+∇⋅(ρuπ+a2u)=ρε(p−π),(ρ^u)t+∇⋅(ρu⊗^u+1M2ψ)=ρε(u−^u),(ρψ)t+∇⋅(ρuψ+a2^u)=ρε(p−ψ). (2.2)

The relaxation model given in (2.2) differs from the one given in [21] in the following points:

1. In the equation for , we use instead of as proposed by the authors in [21]. This is due to the upwind discretization used in [21] which requires in order to have a Mach number independent diffusion of the numerical scheme. Here instead we use centered differences in the implicit part to ensure the Mach number independent diffusion of the numerical scheme.

2. We have simplified the model in the sense that we do not distinguish between the given Mach number and the local Mach number . This is not a restriction in application, because the choice of is given by the application as illustrated in the numerical results. Especially for the Mach number dependent shock test case in Section 7.1.2 we directly compare our results to a scheme that uses the local Mach number .

The following lemma sums up some properties of system (2.2). The proof can be found in [21].

Lemma 1

The relaxation system (2.2) is hyperbolic and is a stable diffusive approximation of (1.1) under the Mach number independent sub-characteristic condition It has the following linearly degenerate eigenvalues

 λu=u, λ±=u±aρ, λ±M=u±aMρ,

where has multiplicity 4. For , the eigenvalues have the ordering

 λ−M<λ−<λu<λ+<λ+M.

In the case of the waves given by and collapse to the waves which have then multiplicity 2 respectively.

To shorten notations we will refer to the original system (1.1) by

 wt+∇⋅f(w)=0, (2.3)

where denotes the physical variables and the flux function is given by . The relaxation model (1.1) is given by

 Wt+∇⋅F(W)=1εR(W), (2.4)

where denotes the state vector, the flux function as defined in (2.2) and the relaxation source term given by

 R(W)=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝000ρ(p−π)ρ(u−^u)ρ(p−ψ)⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠.

The relaxation equilibrium state is defined as

 Weq=(ρ,ρu,E,ρp(ρ,e),ρu,ρp(ρ,e)). (2.5)

The connection between (2.3) and (2.4) can be established, following [17], through the matrix defined by

 Q=(I2+d02+d) (2.6)

where is the dimension. For all states it is satisfied and the physical variables are then recovered by and the fluxes are connected by

3 Time semi-discrete scheme

As we have seen in Lemma 1, the largest eigenvalue of the relaxation model (2.2) tends to infinity as goes to 0. Using a time explicit scheme results in a very restrictive CFL condition. By using an IMEX approach as done in [29, 12], we can avoid the Mach number dependence of the time step.

We rewrite the relaxation system (2.2)/(2.4) in the following form:

 Wt+∇⋅F(W)+1M2∇⋅G(W)=1εR(W). (3.1)

In (3.1), we have split the flux in (2.4) into a flux function which will contain the explicit terms and a flux function which will contain the terms treated implicitly. For efficiency, we want to have as many explicit terms as possible as long as the eigenvalues of are independent of the Mach number. To avoid especially inverting a large non-linear system, we treat the non-linear advection terms explicitly. This results in the following following flux functions

 F(W)=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝ρuρu⊗u+π1+1−M2M2ψ1(E+M2π+(1−M2)ψ)uρπu+a2uρu⊗^uρψu⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠ and G(W)=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝0000ψa2M2^u⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠. (3.2)

We see, that contains the flux of the Euler equations. (1.1) whereas and only act on the relaxation variables. To obtain a time semi-discrete scheme we order the implicit and explicit steps as

 Implicit: Wt+1M2∇⋅G(W) =0, (3.3) Explicit:      Wt+∇⋅F(W) =0, (3.4) Projection:                      Wt =1εR(W). (3.5)

The relaxation source term in (3.5) is solved by projecting the variables onto the equilibrium manifold and thereby reaching the relaxation equilibrium state (2.5). The formal time semi-discrete scheme is then given by

 W(1)−Wn,eq+ΔtM2∇⋅G(W(1)) =0, (3.6) W(2)−W(1)+Δt∇⋅F(W(1)) =0, (3.7) Wn+1=W(2),eq, (3.8)

where we consider the data at time to be at relaxation equilibrium . First we solve the implicit equation (5.1) to gain , followed by the explicit step where we calculate . To get the variables at the new time level, we project onto its equilibrium state. This procedure results into a first order scheme in time.

4 Asymptotic properties

We consider as continuous limit equations the incompressible Euler equations given by

 ρ=const.ut+u⋅∇u+∇Π=0∇⋅u=0 (4.1)

with a dynamical pressure described by .

Since the solution must converge to the solution of the incompressible system (4.1) as , see for example [8, 1], its asymptotic expansion under slipping or periodic boundary conditions must satisfy

 ρ =ρ0+O(M), ρ0 =const. (4.2) u =u0+O(M), ∇⋅u0 =0 (4.3) p =p0+O(M2), p0 =const. (4.4)

We refer to (4.2), (4.3), (4.4) as well-prepared data and summarize it in the following set

 Ωwp={w∈Rd+2,∇ρ0=0,∇p0=0,∇⋅u0=0}. (4.5)

For simplicity we show the AP property for the time semi-implicit scheme. The same steps can be followed with the fully discretized scheme given in Section 5.

To show the AP property we will exploit some properties of the fast pressure obtained in the implicit step (4.6).

4.1 Asymptotic behaviour of ψ(1)

Due to the sparse structure of defined in (3.2), the implicit part reduces to solving only two coupled equations in the relaxation variables given by

 (ρ^u)t+1M2∇ψ=0,(ρψ)t+a2∇⋅^u=0, (4.6)

with the eigenvalues . As done in [10], we rewrite the coupled system (4.6) into one equation for starting from the time-semi-discrete scheme

 ρ(1)−ρnΔt =0, (4.7) (ρ^u)(1)−(ρ^u)nΔt+1M2∇ψ(1) =0, (4.8) (ρψ)(1)−(ρψ)nΔt+a2∇⋅^u(1) =0. (4.9)

To emphasize that (4.6) also depends on the density, we have included the density update (4.7) into the time-semi-discrete system. From equation (4.7) we see that . To simplify notation we define . Inserting (4.8) into (4.9) we can reduce the implicit system to only one equation with an elliptic operator for given as

 ψ(1)−Δt2a2M2τn∇⋅(τn∇ψ(1)) =ψn−Δta2τn∇⋅un. (4.10)

On the right hand side of (4.10) we have already made use of the fact that since we start from equilibrium data. We will see that the correct scaling of with respect to the Mach number is important not only for showing the AP property of the scheme, but also for the positivity of density and internal energy as well as for the Mach number independent diffusion of the fully discretized scheme.

To prevent pressure perturbations at the boundaries which would destroy the well-prepared nature of the pressure, we require boundary conditions on which preserve the scaling of the pressure in time. For a computational domain , we set

 (4.11)
Lemma 2 (Scaling of ψ(1))

Let be equipped with the boundary conditions (4.11). Then the Mach number expansion of after the first stage satisfies

 ψ(1)=p0+M2ψ2+O(M3), (4.12)

where is constant.

Proof. Let us assume that the expansion of is given by

 ψ(1)=ψ0+Mψ1+M2ψ2+O(M3). (4.13)

Since we start in , we have well-prepared data as given in (4.2)(4.3)(4.4). We insert the well-prepared data and (4.13) into the implicit update equation (4.10). Separating the terms, we find using the boundary conditions (4.11)

 {Δψ0=0 in Dψ0=p0 on ∂D.

This leads to on . Separating the terms and using that , we find

 {Δψ1=0 in Dψ1=0 on ∂D (4.14)

which leads to on . As a last step, we collect the terms and using that as well as on . It is not necessary to impose special boundary conditions for . Thus we find

 Δψ2 =0 in D. (4.15)

This shows that the fast pressure after the implicit step is still well-prepared.

4.2 Asymptotic preserving property

We show that the time discretization of (2.2) in the limit coincides with a time discretization of the incompressible Euler equations (4.1). We consider well-prepared data and the Mach number expansion of from Lemma 2. For the total energy defined in (1.2), we find the following Mach number expansion

 E=ρ0e0+M(ρ1e0+ρ0e1)+M2(|u0|2+ρ2e0+ρ1e1+ρ0e2)+O(M3) (4.16)

Inserting the Mach number expansions of and into the density, momentum and energy equation of (3.7) and considering the order terms we have

 ρn+10−ρn0Δt+∇⋅ρn0un0=0,ρn+10un+10−ρn0un0Δt+∇⋅(ρn0un0⊗un0+ψ(1)2)=0,ρn+1en+10−ρn0en0Δt+∇⋅un0(ρn0en0+ψ(1)0)=0,ρn+10en+10−ρn0en0Δt+∇⋅un0(ρn0en0+ψ(1)0)=0 (4.17)

Let us assume, the pressure at time has the Mach number expansion . Note that we have from Lemma 2 that and thus . Since we have and . Further we use . Equipped with that we can simplify (4.17) to

 ρn+10−ρn0Δt =0, (4.18) un+10−un0Δt+un0⋅∇un0+∇ψ(1)2ρn0 =0, (4.19) pn+10−pn0Δt =0. (4.20)

Especially from equations (4.18) and (4.20) we see that and are constants. Looking at the terms we have from the energy equation

 pn+11+Δt∇⋅un1=0. (4.21)

This means the density and pressure at time are well-prepared up to perturbation as in (4.2), (4.4). To be consistent with a time discretization of the incompressible Euler equations (4.1) the divergence of at time defined by has to be at least of order . To show this, we apply the divergence on the velocity update (4.19) which gives

 ∇⋅un+10=∇⋅un0+Δt∇⋅(un0⋅∇un0+Δψ(1)2ρn0). (4.22)

In the proof of Lemma 2, we have shown that on , see (4.15). Using (4.15) together with , we can simplify (4.22) to

 ∇⋅un+10=Δt∇⋅(un0⋅∇un0)=O(Δt).

In summary, we have shown the following theorem.

Theorem 3 (AP property)

Let . Then under the boundary conditions (4.11) the scheme (3.6), (3.7), (3.8) is asymptotic preserving when tends to , in the sense that if then it is and in the limit the time-semi-discrete scheme is a consistent discretization of the incompressible Euler equations (4.1).

5 Derivation of the fully discrete scheme

For simplicity, we develop the fully discretized scheme in one space dimension, but it can be straightforwardly extended to dimensions. In the implicit update (4.10), the space derivatives read

 ∇⋅(τ∇ψ)=∂x1(τ∂x1ψ)+⋯+∂xd(τ∂xdψ) and ∇⋅u=∂x1u1+⋯+∂xdud

and in the explicit part we apply dimensional splitting.

In the following we use a cartesian grid on a computational domain devided in cells of step size . We use a standard finite volume setting, where we define at time the piecewise constant functions

 w(x,tn)=wni, for x∈Ci.

Using this notation, we apply centered differences for the implicit update (4.10) and have

 ψ(1)i−Δt2Δx2a2M2τni(τni−1/2ψ(1)i−1−(τni−1/2+τni+1/2)ψ(1)i+τni+1/2ψ(1)i+1)=ψni−Δt2Δxa2τni(uni+1−uni−1), (5.1)

where .

For the explicit part, we will use a Godunov type finite volume scheme following [30] which we will describe in the following section.

5.1 Godunov type finite volume scheme

In the explicit step we consider the following equations as defined in (3.2), (3.7)

 ∂tρ+∂xρu=0,∂tρu+∂x(ρu2+π+1−M2M2ψ)=0,∂tE+∂x((E+M2π+(1−M2)ψ)u)=0,∂tρπ+∂x((ρπ+a2)u)=0,∂tρ^u+∂x(ρu^u)=0,∂tρψ+∂x(ρψu)=0. (5.2)

In the next result the structure of system (5.2) is summarized.

Lemma 4

System (5.2) admits the linear degenerate eigenvalues and , where the eigenvalue has multiplicity 4. The relaxation parameter is independent of the Mach number as well as all eigenvalues. The Riemann Invariants with respect to are

 Iu1=u, Iu2=M2π+(1−M2)ψ (5.3)

and with respect to

 I±1=u±aρ, I±2=π∓au,I±3=e−M22a2π2−1−M2a2πψ,I±4=^u,   I±5=ψ. (5.4)

Proof. We rewrite the equations (5.2) using primitive variables in non-conservative form

 Vt+A(V)Vx=0, (5.5)

where the matrix is given by

 A(V)=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝uρ00000u01ρ01−M2M20M2π+(1−M2)ψρu0000a2ρ0u000000u000000u⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠. (5.6)

It is easy to check that are eigenvalues of

. Associated to the eigenvalues, we find the respective eigenvectors

 ru1=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝0001−1M201⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠ru2=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝000010⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠ru3=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝001000⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠ru4=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝100000⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠r±=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝ρ2a2±1aM2π+(1−M2)ψa2100⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠

We see that , and . Thus all eigenvalues are linearly degenerate. A scalar function is a Riemann invariant if for all eigenvectors associated to an eigenvalue , holds. It is straightforward to check, that (5.3) and (5.4) are Riemann invariants. Since Riemann invariants are invariant under change of variables, the Riemann invariants of (5.5) are the same as for the equations in conservation form (5.2). For more details, see [26, 21].

We will follow the theory of Harten, Lax and van Leer [30] for deriving an approximate Riemann solver based on the states after the implicit step. Due to the linear-degeneracy from Lemma 4, the structure of the approximate Riemann solver, as displayed in Figure 1, is given as follows

 WRS(xt;W(1)L,W(1)R)=⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩W(1)Lxt<λ−,W∗Lλ−

To compute the intermediate states , we use the Riemann invariants as given in Lemma 4. Note that since the eigenvalues have multiplicity 1, we get the expected 5 Riemann invariants. This does not hold in general for eigenvalues with multiplicity larger than 1. Nevertheless, the invariants (5.3) and (5.4) give enough relations to determine the solution to a Riemann problem for (5.2) as shown in the following lemma.

Lemma 5

Consider an initial value problem with initial data given by

 W0(x)={WLx<0,WRx>0. (5.8)

Then the solution consists of four constant states separated by contact discontinuities with the structure given in (5.7). The solution for the states is given by

 1ρ∗L=1ρL+1a(u∗−uL),1ρ∗R=1ρR+1a(uR−u∗),u∗=u∗L=u∗R=12(uL+uR)+12a[(πL−πR)+1−M2M2(ψL−ψR)],π∗L=12(πL+πR)+a2(uL−uR)−1−M2M212(ψL−ψR),π∗R=12(πL+πR)+a2(uL−uR)+1−M2M212(ψL−ψR),e∗L=eL−12a2[(π2L−(π∗L)2)+(1−M2)(πL−π∗L)ψL],e∗R=eR−12a2[(π2R−(π∗R)2)+(1−M2)(πR−π∗R)ψR],ψ∗L,R=ψL,R,^u∗L,R=^uL,R. (5.9)

Proof. Since , we have the following order of the eigenvalues . The solution structure follows directly from the linear degeneracy of the eigenvalues given in Lemma 4 and the ordering of the eigenvalues. To derive the solution for the intermediate states one uses the Riemann invariants given in (5.3) and (5.4) and solves the resulting system of equations.

Given the solution of the Riemann problem (5.7), we define the numerical fluxes as follows

 Fi+1/2=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩F(W(1)i)λ−>0F(W∗,(1)i)λu>0>λ−F(W∗,(1)i+1)λ+>0>λuF(W(1)i+1)λ+<0 (5.10)

where the superscript emphasizes that the states after the implicit step are used. To avoid interactions between the approximate Riemann solvers at the interfaces , we have a CFL restriction on the time step of

 Δt≤12Δxmaxi|ui±a/ρi| (5.11)

which is independent of the Mach number. This leads to the following update of the explicit part

 W(2)i=W(1)i−ΔtΔx(Fi+1/2−Fi−1/2). (5.12)

At this point, we can make some remarks on the numerical scheme as it is given by (5.1),(5.12) and (3.8).

Remark 6

The flux function in (3.2) in the implicit part and the relaxation source term only act on the relaxation variables . Thus, we can immediately give the update of the physical variables as

 wn+1i=wni−ΔtΔx(f(QWRS(0;W(1)i,W(1)i+1))−f(QWRS(0;W(1)i−1,W(1)i))), (5.13)

where the matrix is the projection on the first components as defined in (2.6).

Remark 7

We see from (5.1) and (5.9), that the updates in the implicit and explicit step are independent of . Therefore it is not necessary to explicitly compute in the scheme thus reducing the computational costs.

Remark 8

For the limit case the relaxation model (2.2) reduces to a Suliciu relaxation model for the compressible Euler equations since all -terms are being cancelled and the eigenvalues and collapse. Analogously, in the scheme the implicit step becomes redundant, since in the Riemann solution (5.9) the -terms are cancelled, too. Thus the scheme reduces to a Godunov-type approximate Riemann solver based explicit scheme for the compressible Euler equations as can be found in [26].

For the order of the eigenvalues changes to . The scheme is still stable because the CFL condition (5.11) comes from the explicit eigenvalues which are now the largest ones.

5.2 Positivity of density and internal energy

For physical applications it is necessary that the density and internal energy be positive. We define the physical admissible states associated with the Euler equations (1.1) as

 Ωphy={w∈Rd+2,ρ>0,e>0}. (5.14)

The property of our scheme to preserve the domain is linked to how the fluxes are calculated. In our case it is essential that the density and internal energy of the Riemann solution (5.7),(5.9) are positive. This is shown in the following lemma.

Lemma 9

Let the initial data of the Riemann problem (5.8) be given as satisfying the boundary conditions (4.11). Then there is a relaxation parameter large enough but independent of such that .

Proof. Since the proof only concerns data after the implicit step, we will drop the superscript . We have to prove, that and . From the ordering of eigenvalues and the formula for from (5.9), we get

 1ρ∗L=1ρL+u∗−uLa≥1ρL−1ρL=0.

Analogously, we find . For the internal energy , we insert the definition of from (5.9) into to obtain a formula depending only on the left and right states

 e∗L=eL+18(uL−uR)2+12a2(−π2L+14(πL+πR+1−M2M2(ψR−ψL))2+12ψL(1−M2)(πR−πL+1−M2M2(ψR−ψL)))−14a(uL−uR)(πL+πR+1−M2M2(ψR−ψL)+(1−M2)ψL).

From Lemma 2, we know