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

TL;DR

This paper introduces a second-order IMEX relaxation scheme for the Euler equations that is positivity-preserving, asymptotic-preserving, and applicable in 1D and 2D, improving numerical stability and accuracy in fluid simulations.

Contribution

The paper develops a novel all-speed second-order IMEX relaxation scheme for Euler equations that maintains positivity and asymptotic properties, with demonstrated convergence and applicability.

Findings

Scheme is positivity preserving for density and internal energy.

Scheme is asymptotic preserving towards incompressible Euler equations.

Numerical experiments confirm convergence and effectiveness in 1D and 2D.

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.

Tables

Table 1. Table 1: L 1 superscript 𝐿 1 L^{1} -error and convergence rates for the solution of the smooth Gresho vortex test at T = 0.05 𝑇 0.05 T=0.05 (non-dimensional).

$M$	$N$	$\rho$		$u_1$		$u_2$		$p$
${10}^{-1}$	20	1.810$\cdot {10}^{-3}$	—	1.729$\cdot {10}^{-2}$	—	1.729$\cdot {10}^{-2}$	—	1.921$\cdot {10}^{-3}$	—
	40	3.705$\cdot {10}^{-4}$	2.288	5.070$\cdot {10}^{-3}$	1.770	5.070$\cdot {10}^{-3}$	1.770	3.956$\cdot {10}^{-4}$	2.279
	60	1.246$\cdot {10}^{-4}$	2.688	2.403$\cdot {10}^{-3}$	1.842	2.403$\cdot {10}^{-3}$	1.842	1.343$\cdot {10}^{-3}$	2.665
	80	5.510$\cdot {10}^{-5}$	2.835	1.396$\cdot {10}^{-3}$	1.887	1.396$\cdot {10}^{-3}$	1.887	5.922$\cdot {10}^{-5}$	2.845
${10}^{-2}$	20	1.812$\cdot {10}^{-3}$	—	1.731$\cdot {10}^{-2}$	—	1.731$\cdot {10}^{-2}$	—	1.912$\cdot {10}^{-3}$	—
	40	3.582$\cdot {10}^{-4}$	2.339	5.057$\cdot {10}^{-3}$	1.775	5.057$\cdot {10}^{-3}$	1.775	3.781$\cdot {10}^{-4}$	2.337
	60	1.162$\cdot {10}^{-4}$	2.777	2.402$\cdot {10}^{-3}$	1.837	2.402$\cdot {10}^{-3}$	1.837	1.226$\cdot {10}^{-4}$	2.777
	80	4.881$\cdot {10}^{-5}$	3.014	1.386$\cdot {10}^{-3}$	1.910	1.386$\cdot {10}^{-3}$	1.910	5.151$\cdot {10}^{-5}$	3.014
${10}^{-3}$	20	1.811$\cdot {10}^{-3}$	—	1.731$\cdot {10}^{-2}$	—	1.731$\cdot {10}^{-2}$	—	1.912$\cdot {10}^{-3}$	—
	40	3.580$\cdot {10}^{-4}$	2.339	5.057$\cdot {10}^{-3}$	1.775	5.057$\cdot {10}^{-3}$	1.775	3.778$\cdot {10}^{-4}$	2.339
	60	1.162$\cdot {10}^{-4}$	2.775	2.402$\cdot {10}^{-3}$	1.836	2.402$\cdot {10}^{-3}$	1.836	1.227$\cdot {10}^{-4}$	2.775
	80	4.875$\cdot {10}^{-5}$	3.019	1.382$\cdot {10}^{-3}$	1.920	1.382$\cdot {10}^{-3}$	1.920	5.146$\cdot {10}^{-5}$	3.019

Equations

$$\begin{array}{rl}{\rho }_t+\nabla \cdot (\rho \mathbf{u})& =0\\ (\rho u{)}_t+\nabla \cdot (\rho \mathbf{u}\otimes \mathbf{u}+\frac{1}{M^2}p\mathbb{I})& =0\\ E_t+\nabla \cdot (\mathbf{u}(E+p))& =0,\end{array}$$

$$E=\rho e+\frac{1}{2}{M}^2\rho |u{|}^2$$

$$(\rho \pi {)}_t+\nabla \cdot (\rho \pi \mathbf{u})+a^2\nabla \cdot \mathbf{u}=\frac{\rho }{ϵ}(p-\pi ).$$

$$\frac{p}{M^2}=p+\frac{1-M^2}{M^2}p.$$

$$(\rho \mathbf{u}{)}_t+\nabla \cdot \left(\rho \mathbf{u}\otimes \mathbf{u}+\pi +\frac{1-M^2}{M^2}\psi \right)=0$$

$$\begin{array}{rl}{\rho }_t+\nabla \cdot (\rho u)& =0,\\ (\rho \mathbf{u}{)}_t+\nabla \cdot \left(\rho \mathbf{u}\otimes \mathbf{u}+\pi +\frac{1-M^2}{M^2}\psi \right)& =0,\\ E_t+\nabla \cdot (\mathbf{u}(E+M^2\pi +(1-M^2)\psi ))& =0,\\ (\rho \pi {)}_t+\nabla \cdot (\rho \mathbf{u}\pi +a^2\mathbf{u})& =\frac{\rho }{\epsilon }(p-\pi ),\\ (\rho \hat{\mathbf{u}}{)}_t+\nabla \cdot (\rho \mathbf{u}\otimes \hat{\mathbf{u}}+\frac{1}{M^2}\psi )& =\frac{\rho }{\epsilon }(\mathbf{u}-\hat{\mathbf{u}}),\\ (\rho \psi {)}_t+\nabla \cdot (\rho \mathbf{u}\psi +a^2\hat{\mathbf{u}})& =\frac{\rho }{\epsilon }(p-\psi ).\end{array}$$

$${\lambda }^u=\mathbf{u},\text{leavevmode}{\lambda }^{\pm }=\mathbf{u}\pm \frac{a}{\rho },\text{leavevmode}{\lambda }_M^{\pm }=\mathbf{u}\pm \frac{a}{M\rho },$$

$${\lambda }_M^{-}<{\lambda }^{-}<{\lambda }^u<{\lambda }^{+}<{\lambda }_M^{+}.$$

$$w_t+\nabla \cdot f(w)=0,$$

$$W_t+\nabla \cdot \mathcal{F}(W)=\frac{1}{\epsilon }R(W),$$

$$R(W)=\left(\begin{array}{c}0\\ 0\\ 0\\ \rho (p-\pi )\\ \rho (\mathbf{u}-\hat{\mathbf{u}})\\ \rho (p-\psi )\end{array}\right).$$

$$W^{\text{eq}}=(\rho ,\rho \mathbf{u},E,\rho p(\rho ,e),\rho \mathbf{u},\rho p(\rho ,e)).$$

$$Q=\left(\begin{array}{cc}{\mathbb{I}}_{2+d}& 0_{2+d}\end{array}\right)$$

$$W_t+\nabla \cdot F(W)+\frac{1}{M^2}\nabla \cdot G(W)=\frac{1}{\epsilon }R(W).$$

$$F(W)=\left(\begin{array}{c}\rho \mathbf{u}\\ \rho \mathbf{u}\otimes \mathbf{u}+\pi \text{mathbbm}1+\frac{1-M^2}{M^2}\psi \text{mathbbm}1\\ (E+M^2\pi +(1-M^2)\psi )\mathbf{u}\\ \rho \pi \mathbf{u}+a^2\mathbf{u}\\ \rho \mathbf{u}\otimes \hat{\mathbf{u}}\\ \rho \psi \mathbf{u}\end{array}\right)\text{ and }G(W)=\left(\begin{array}{c}0\\ 0\\ 0\\ 0\\ \psi \\ a^2{M}^2\hat{\mathbf{u}}\end{array}\right).$$

$$\text{Implicit: }{W}_t+\frac{1}{M^2}\nabla \cdot G(W)$$

$$\text{Explicit: }\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}{W}_t+\nabla \cdot F(W)$$

$$\text{Projection: }\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}{W}_t$$

$$W^{(1)}-W^{n,\text{eq}}+\frac{\Delta t}{M^2}\nabla \cdot G(W^{(1)})$$

$$W^{(2)}-W^{(1)}+\Delta t\nabla \cdot F(W^{(1)})$$

$$W^{n+1}=W^{(2),\text{eq}},$$

$$\begin{array}{rl}\rho & =\text{const.}\\ u_t+u\cdot \nabla u+\nabla \Pi & =0\\ \nabla \cdot u& =0\end{array}$$

$$\rho$$

$$u$$

$$p$$

$${\Omega }_{wp}=\left\{w\in {\mathbb{R}}^{d+2},\nabla {\rho }_0=0,\nabla p_0=0,\nabla \cdot {\mathbf{u}}_0=0\right\}.$$

$$\begin{array}{rl}(\rho \hat{\mathbf{u}}{)}_t+\frac{1}{M^2}\nabla \psi & =0,\\ (\rho \psi {)}_t+a^2\nabla \cdot \hat{\mathbf{u}}& =0,\end{array}$$

$$\frac{{\rho }^{(1)}-{\rho }^n}{\Delta t}$$

$$\frac{(\rho \hat{\mathbf{u}}{)}^{(1)}-(\rho \hat{\mathbf{u}}{)}^n}{\Delta t}+\frac{1}{M^2}\nabla {\psi }^{(1)}$$

$$\frac{(\rho \psi {)}^{(1)}-(\rho \psi {)}^n}{\Delta t}+a^2\nabla \cdot {\hat{\mathbf{u}}}^{(1)}$$

$${\psi }^{(1)}-\frac{\Delta t^2{a}^2}{M^2}{\tau }^n\nabla \cdot ({\tau }^n\nabla {\psi }^{(1)})$$

$$\begin{array}{ll}{\psi }_0^{(1)}& =p_0^n\\ {\psi }_1^{(1)}& =0\end{array}\}\text{ on }\partial D.$$

$${\psi }^{(1)}=p_0+M^2{\psi }_2+\mathcal{O}(M^3),$$

$${\psi }^{(1)}={\psi }_0+M{\psi }_1+M^2{\psi }_2+\mathcal{O}(M^3).$$

Full text

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

Andrea Thomann111Dipartimento di Scienze e Alta Tecnologia, Università degli Studi dell’Insubria, Via Valleggio 11, 22100 Como, Italy 222Correspondence to: Andrea Thomann, University of Insubria, Via Valleggio 11, 22100 Como (CO), Italy, Email: acthomann@uninsubria.it 333Marie Skłodowska-Curie fellow of the Istituto Nazionale di Alta Matematica Francesco Severi, Rome, Italy, Markus Zenk444Fakultät für Mathematik und Informatik, Universität Würzburg, Emil-Fischer-Str. 40, 97074 Würzburg, Germany, Gabriella Puppo555Dipartimento di Matematica, La Sapienza Università di Roma, Piazzale Aldo Moro 5, 00185 Roma, Italy,

Christian Klingenberg44footnotemark: 4

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 $d$-space dimensions which are given by the following set of equations [1]

[TABLE]

where the total energy is given by

[TABLE]

and $e>0$ denotes the internal energy. The density is denoted by $\rho>0$, $\mathbf{u}\in\mathbb{R}^{d}$ is the velocity vector and $M$ 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 $M$ 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 $M$ 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 $\pi$,

[TABLE]

Here $\varepsilon$ denotes the relaxation time. Equation (2.1) is derived by multiplying the density equation by $\partial_{\rho}p$. The non-linearity $\rho^{2}\partial_{\rho}p(\rho,e)$ that thereby arises is replaced by a constant parameter $a^{2}$, 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 $a>\rho\sqrt{\partial_{\rho}p(\rho,e)}>0$ [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

[TABLE]

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

[TABLE]

The evolution of the new variables $\pi$ and $\psi$ are then developed in the spirit of a Suliciu relaxation approach. The evolution for $\pi$ is given by the Suliciu relaxation equation (2.1). However, applying the standard Suliciu relaxation method also on the pressure $\psi$, 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 $\mathbf{\hat{u}}\in\mathbb{R}^{d}$ which is relaxed to the system velocity $\mathbf{u}$ and couples to the pressure $\psi$. The form of the evolution equations for $\mathbf{\hat{u}}$ and $\psi$ 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

[TABLE]

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

In the equation for $\mathbf{\hat{u}}$, we use $\frac{1}{M^{2}}$ instead of $\frac{1}{M^{4}}$ as proposed by the authors in [21]. This is due to the upwind discretization used in [21] which requires $\frac{1}{M^{4}}$ 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. 2.

We have simplified the model in the sense that we do not distinguish between the given Mach number $M$ and the local Mach number $M_{loc}$. This is not a restriction in application, because the choice of $M$ 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 $M_{loc}$.

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 $a>\rho\sqrt{\partial_{\rho}p(\rho,e)}.$ It has the following linearly degenerate eigenvalues

[TABLE]

where $\lambda^{u}$ has multiplicity 4. For $M<1$, the eigenvalues have the ordering

[TABLE]

In the case of $M=1$ the waves given by $\lambda^{\pm}$ and $\lambda_{M}^{\pm}$ collapse to the waves $\lambda^{\pm}$ which have then multiplicity 2 respectively.

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

[TABLE]

where $w=(\rho,\rho\mathbf{u},E)^{T}$ denotes the physical variables and the flux function is given by $f(w)=\begin{pmatrix}\rho\mathbf{u}\\ \rho\mathbf{u}\otimes\mathbf{u}+\frac{p}{M^{2}}\mathbb{I}\\ \mathbf{u}(E+p)\end{pmatrix}$. The relaxation model (1.1) is given by

[TABLE]

where $W=(\rho,\rho\mathbf{u},E,\rho\pi,\rho\mathbf{\hat{u}},\rho\psi)$ denotes the state vector, $\mathcal{F}$ the flux function as defined in (2.2) and $R$ the relaxation source term given by

[TABLE]

The relaxation equilibrium state is defined as

[TABLE]

The connection between (2.3) and (2.4) can be established, following [17], through the matrix $Q\in\mathbb{R}^{2+d\times 2(2+d)}$ defined by

[TABLE]

where $d$ is the dimension. For all states $W$ it is satisfied $QR(W)=0$ and the physical variables are then recovered by $w=QW$ and the fluxes are connected by $f(w)=Q\mathcal{F}(W^{\text{eq}}).$

3 Time semi-discrete scheme

As we have seen in Lemma 1, the largest eigenvalue $|\lambda_{M}^{\pm}|$ of the relaxation model (2.2) tends to infinity as $M$ 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:

[TABLE]

In (3.1), we have split the flux $\mathcal{F}$ in (2.4) into a flux function $F$ which will contain the explicit terms and a flux function $G$ 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 $F$ 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

[TABLE]

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

[TABLE]

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

[TABLE]

where we consider the data at time $t^{n}$ to be at relaxation equilibrium $W^{n,\text{eq}}$. First we solve the implicit equation (5.1) to gain $W^{(1)}$, followed by the explicit step where we calculate $W^{(2)}$. To get the variables $W^{n+1}$ at the new time level, we project $W^{(2)}$ 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

[TABLE]

with a dynamical pressure described by $\Pi$.

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

[TABLE]

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

[TABLE]

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 $\psi^{(1)}$ obtained in the implicit step (4.6).

4.1 Asymptotic behaviour of $\psi^{(1)}$

Due to the sparse structure of $G$ defined in (3.2), the implicit part reduces to solving only two coupled equations in the relaxation variables $\mathbf{\hat{u}},\psi$ given by

[TABLE]

with the eigenvalues $\lambda_{M}=\pm\frac{a}{\rho M}$. As done in [10], we rewrite the coupled system (4.6) into one equation for $\psi$ starting from the time-semi-discrete scheme

[TABLE]

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 $\rho^{(1)}=\rho^{n}$. To simplify notation we define $\tau_{i}^{n}=\frac{1}{\rho_{i}^{n}}$. Inserting (4.8) into (4.9) we can reduce the implicit system to only one equation with an elliptic operator for $\psi$ given as

[TABLE]

On the right hand side of (4.10) we have already made use of the fact that $\hat{u}^{n}=u^{n}$ since we start from equilibrium data. We will see that the correct scaling of $\psi^{(1)}$ 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 $\mathcal{O}(M)$ pressure perturbations at the boundaries which would destroy the well-prepared nature of the pressure, we require boundary conditions on $\psi$ which preserve the scaling of the pressure in time. For a computational domain $D$, we set

[TABLE]

Lemma 2** (Scaling of $\psi^{(1)}$)**

Let $w^{n}\in\Omega_{wp}$ be equipped with the boundary conditions (4.11). Then the Mach number expansion of $\psi^{(1)}$ after the first stage satisfies

[TABLE]

where $p_{0}$ is constant.

Proof. Let us assume that the expansion of $\psi^{(1)}$ is given by

[TABLE]

Since we start in $\Omega_{wp}$, we have well-prepared data as given in (4.2)(4.3)(4.4). We insert the well-prepared data $w^{n}$ and (4.13) into the implicit update equation (4.10). Separating the $\mathcal{O}(M^{-2})$ terms, we find using the boundary conditions (4.11)

[TABLE]

This leads to $\psi_{0}=p_{0}$ on $\overline{D}$. Separating the $\mathcal{O}(M^{-1})$ terms and using that $\psi_{0}=p_{0}=\text{const}$, we find

[TABLE]

which leads to $\psi_{1}=0$ on $\overline{D}$. As a last step, we collect the $\mathcal{O}(M^{-1})$ terms and using that $\psi_{0}=p_{0}$ as well as $\psi_{1}=0$ on $\overline{D}$. It is not necessary to impose special boundary conditions for $\psi_{2}$. Thus we find

[TABLE]

This shows that the fast pressure $\psi^{(1)}$ after the implicit step is still well-prepared.

4.2 Asymptotic preserving property

We show that the time discretization of (2.2) in the $M\to 0$ limit coincides with a time discretization of the incompressible Euler equations (4.1). We consider well-prepared data $w^{n}\in\Omega_{wb}$ and the Mach number expansion of $\psi^{(1)}$ from Lemma 2. For the total energy defined in (1.2), we find the following Mach number expansion

[TABLE]

Inserting the Mach number expansions of $w^{n},\psi^{(1)}$ and $E^{n}$ into the density, momentum and energy equation of (3.7) and considering the $\mathcal{O}(M^{0})$ order terms we have

[TABLE]

Let us assume, the pressure at time $t^{n+1}$ has the Mach number expansion $p^{n+1}=p_{0}^{n+1}+Mp_{1}^{n+1}+\mathcal{O}(M^{2})$. Note that we have from Lemma 2 that $\psi_{0}^{(1)}=p_{0}$ and thus $\nabla\psi_{0}^{(1)}=\nabla p_{0}=0$. Since $w^{n}\in\Omega_{wp}$ we have $\nabla\rho_{0}^{n}=0$ and $\nabla\cdot u_{0}^{n}=0$. Further we use $p_{0}=(\gamma-1)\rho_{0}e_{0}$. Equipped with that we can simplify (4.17) to

[TABLE]

Especially from equations (4.18) and (4.20) we see that $\rho_{0}^{n+1}$ and $p_{0}^{n+1}$ are constants. Looking at the $\mathcal{O}(M^{1})$ terms we have from the energy equation

[TABLE]

This means the density and pressure at time $t^{n+1}$ are well-prepared up to $\mathcal{O}(\Delta t)$ perturbation as in (4.2), (4.4). To be consistent with a time discretization of the incompressible Euler equations (4.1) the divergence of $\mathbf{u}_{0}$ at time $t^{n+1}$ defined by $\nabla\cdot\mathbf{u}_{0}^{n+1}$ has to be at least of order $\Delta t$. To show this, we apply the divergence on the velocity update (4.19) which gives

[TABLE]

In the proof of Lemma 2, we have shown that $\Delta\psi_{2}^{(1)}=0$ on $\partial D$, see (4.15). Using (4.15) together with $\nabla\cdot u^{n}_{0}=0$, we can simplify (4.22) to

[TABLE]

In summary, we have shown the following theorem.

Theorem 3** (AP property)**

Let $w^{n}\in\Omega_{wp}$. Then under the boundary conditions (4.11) the scheme (3.6), (3.7), (3.8) is asymptotic preserving when $M$ tends to [math], in the sense that if $w^{n}\in\Omega_{wp}$ then it is $w^{n+1}\in\Omega_{wp}$ and in the limit $M\to 0$ 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 $d$ dimensions. In the implicit update (4.10), the space derivatives read

[TABLE]

and in the explicit part we apply dimensional splitting.

In the following we use a cartesian grid on a computational domain $D$ devided in $N$ cells $C_{i}=(x_{i-1/2},x_{i+1/2})$ of step size $\Delta x$. We use a standard finite volume setting, where we define at time $t^{n}$ the piecewise constant functions

[TABLE]

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

[TABLE]

where $\tau_{i+1/2}=\frac{1}{2}\left(\tau_{i+1}+\tau_{i}\right)$.

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)

[TABLE]

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

Lemma 4

System (5.2) admits the linear degenerate eigenvalues $\lambda^{\pm}=u\pm\frac{a}{\rho}$ and $\lambda^{u}=u$, where the eigenvalue $\lambda^{u}$ has multiplicity 4. The relaxation parameter $a$ is independent of the Mach number as well as all eigenvalues. The Riemann Invariants with respect to $\lambda^{u}$ are

[TABLE]

and with respect to $\lambda^{\pm}$

[TABLE]

Proof. We rewrite the equations (5.2) using primitive variables $V=(\rho,u,e,\pi,\hat{u},\psi)$ in non-conservative form

[TABLE]

where the matrix $\mathcal{A}(V)$ is given by

[TABLE]

It is easy to check that $\lambda^{u},\lambda^{\pm}$ are eigenvalues of $\mathcal{A}(V)$. Associated to the eigenvalues, we find the respective eigenvectors

[TABLE]

We see that $\partial_{V}\lambda^{u}\cdot r_{i}=0$, $\forall i=1,\dots,4$ and $\partial_{V}\lambda^{\pm}\cdot r^{\pm}=0$. Thus all eigenvalues are linearly degenerate. A scalar function $I(V)$ is a Riemann invariant if for all eigenvectors $r$ associated to an eigenvalue $\lambda$, $\partial_{V}I(V)\cdot r=0$ 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 $W_{\mathcal{RS}}\left(\frac{x}{t};W_{L}^{(1)},W_{R}^{(1)}\right)$ based on the states $W^{(1)}$ 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

[TABLE]

To compute the intermediate states $W_{L,R}^{\ast}$, we use the Riemann invariants as given in Lemma 4. Note that since the eigenvalues $\lambda^{\pm}$ 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 $W=W^{(1)}$ given by

[TABLE]

Then the solution consists of four constant states separated by contact discontinuities with the structure given in (5.7). The solution for the states $W^{\ast}_{L/R}$ is given by

[TABLE]

Proof. Since $\frac{a}{\rho}>0$, we have the following order of the eigenvalues $u-\frac{a}{\rho}<u<u+\frac{a}{\rho}$. 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 $W^{\ast}_{L,R}$ 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

[TABLE]

where the superscript $(1)$ emphasizes that the states after the implicit step are used. To avoid interactions between the approximate Riemann solvers at the interfaces $x_{i+1/2}$, we have a CFL restriction on the time step of

[TABLE]

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

[TABLE]

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 $G$ in (3.2) in the implicit part and the relaxation source term $R$ only act on the relaxation variables $(\pi,\mathbf{\hat{u}},\psi)$. Thus, we can immediately give the update of the physical variables $w$ as

[TABLE]

where the matrix $Q$ is the projection on the first $2+d$ 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 $\mathbf{\hat{u}}$. Therefore it is not necessary to explicitly compute $\mathbf{\hat{u}}$ in the scheme thus reducing the computational costs.

Remark 8

For the limit case $M=1$ the relaxation model (2.2) reduces to a Suliciu relaxation model for the compressible Euler equations since all $\psi$-terms are being cancelled and the eigenvalues $\lambda^{\pm}$ and $\lambda_{M}^{\pm}$ collapse. Analogously, in the scheme the implicit step becomes redundant, since in the Riemann solution (5.9) the $\psi$-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 $M>1$ the order of the eigenvalues changes to $\lambda^{-}<\lambda_{M}^{-}<\lambda^{u}<\lambda_{M}^{+}<\lambda^{+}$. 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

[TABLE]

The property of our scheme to preserve the domain $\Omega_{phy}$ 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 $W_{L,R}^{(1)}\in\Omega_{phy}\cap\Omega_{wp}$ satisfying the boundary conditions (4.11). Then there is a relaxation parameter large enough but independent of $M$ such that $W_{\mathcal{RS}}(\frac{x}{t},W_{L}^{(1)},W_{R}^{(1)})\in\Omega_{phy}$.

Proof. Since the proof only concerns data after the implicit step, we will drop the superscript $(1)$. We have to prove, that $\rho_{L,R}^{\ast}>0$ and $e_{L,R}^{\ast}>0$. From the ordering of eigenvalues $u_{L}-\frac{a}{\rho_{L}}<u^{\ast}$ and the formula for $\rho_{L}^{\ast}$ from (5.9), we get

[TABLE]

Analogously, we find $\rho_{R}^{\ast}\geq 0$. For the internal energy $e_{L}^{\ast}$, we insert the definition of $\pi_{L}^{\ast}$ from (5.9) into $e_{L}^{\ast}$ to obtain a formula depending only on the left and right states $W_{L,R}$

[TABLE]

From Lemma 2, we know $\psi_{L,R}=p_{0}+\mathcal{O}(M^{2})$. The difference $\psi_{R}-\psi_{L}=\mathcal{O}(M^{2})$ cancels with the $1/M^{2}$ which leads to $e_{L}^{\ast}=\mathcal{O}(1)$. All possibly negative terms in $e_{L}^{\ast}$ can then be controlled by the relaxation parameter $a>0$ independent of the Mach number. The same argument holds for $e_{R}^{\ast}$.

The positivity property for the first order scheme is given in the next result.

Theorem 10** (Positivity property 1)**

Let the initial state be given as

[TABLE]

in $d$ dimensions satisfying the boundary condition described in (4.11). Then under the Mach number independent CFL condition $\frac{\Delta t}{\Delta x}\underset{i}{\max}|\lambda^{\pm}(w_{i}^{n})|<\frac{1}{2d}$, the numerical scheme defined by (5.1),(5.13) preserves the positivity of density and internal energy, that is

[TABLE]

Proof. Due to the construction of the numerical scheme, the update of the physical variables is only done through the explicit step (5.13). Therefore we can adopt the proof of Theorem 3 in [31] using the positivity of the Riemann solver from Lemma 9. The key element is, that we can write the update $w^{n+1}_{i}$ in one dimension as a convex combination of Riemann solvers which satisfy $W_{\mathcal{RS}}\in\Omega_{phy}$ according to Lemma 9:

[TABLE]

Since we are using dimensional splitting, the update in $d$ dimensions can be written as a sum of updates as given in (5.15) in each dimension and due to convexity we have $w_{i}^{n+1}\in\Omega_{phy}$. More details can be found for example in [31].

5.3 Mach number independent diffusion

Although we are using a Godunov type upwind scheme in the explicit part, our scheme does not suffer from an excessive numerical diffusion as $M$ tends to 0. As we will show in the following this is due to the well-prepared implicitly treated fast pressure $\psi^{(1)}$. In order to do so, we investigate the numerical diffusion vector $\mathcal{D}\in\mathbb{R}^{2+d}$ defined by

[TABLE]

where the Matrix $Q$ is defined in (2.6). Given well-prepared initial data $w^{n}_{i}\in\Omega_{wp}$, we have the following Mach number expansion for the physical variables as given in (4.2),(4.3) and (4.4)

[TABLE]

From Lemma 2, we have for $\psi$ after the implicit step

[TABLE]

The Mach number expansion of the states $W^{(1)}_{i},W^{(1)}_{i+1}$ used in the Riemann solver $W_{\mathcal{RS}}$ is composed of the expansions (5.21) and (5.22). Inserting them into the formulas of the intermediate states (5.9) of the Riemann solution (5.7), we have the following scaling of $W_{i}^{\ast},W_{i+1}^{\ast}$ with respect to the Mach number

[TABLE]

From (5.28) it is evident that the lowest order of $M$ in the intermediate states is $\mathcal{O}(M^{0})$. Inserting (5.28) in the interface flux gives

[TABLE]

Therefore, using (5.29) in (5.16), the diffusion vector is given by

[TABLE]

This shows that the diffusion introduced by the Riemann solver does not suffer from a $\mathcal{O}(M^{-1})$ dependent diffusion in the momentum equation.

6 Second order extension

In this section we extend the first order scheme to second order accuracy. We seek a natural extension of the first order scheme that preserves the positivity property as shown in Theorem 10.

6.1 Second order time integration scheme

For the second order extension, we use a modified two stage time integration as can be found in [25]. It is a convex combination of first order steps and is given by

[TABLE]

The relaxation equilibrium states $W^{\text{eq},n}$ and $W^{\text{eq},(2)}$ are defined as in (2.5). As in the first order case, we can rewrite the integration scheme (6.1) in terms of the update for the physical variables given in Remark 6. The first order time update (5.13) can be written as

[TABLE]

where $Q$ is defined in (2.6). Based on this we can reformulate (6.1) in terms of updates of the physical variables as

[TABLE]

where we obtain $W^{(1)}$ and $\overline{W}^{(1)}$ from the implicit step (3.6) using the update formula (5.1) to compute $\psi^{(1)}$ and $\overline{\psi}^{(1)}$. That the time integration scheme is second order accurate is numerically validated in Section 7. The time integration scheme (6.2) can be extended to variable step sizes $\Delta t_{1},\Delta t_{2}$ for each stage respectively as given in [25]. This has the advantage that the CFL criterion can be met at each stage independently. It is given by

[TABLE]

6.2 Second order space reconstruction

Following Remark 6, we focus on the update of the physical variables in the explicit step. As it is standard in the finite volume setting, we apply a reconstruction to get a higher accuracy for the interface values $w_{i+1/2}$. To get second order, we consider piecewise linear functions in the conserved variables $w=(\rho,\rho\mathbf{u},E)$ at time level $t^{n}$ and in $\psi^{(1)}$. As we are working on a Cartesian grid, we reconstruct along each dimension separately. In one dimension the linear function in $(w_{i},x_{i})$ on $(x_{i-1/2},x_{i1/2})$ is given by

[TABLE]

where $\sigma_{i}=(\sigma_{i}^{\rho},\sigma_{i}^{\rho u},\sigma_{i}^{E})$. For the multi-dimensional notation, we refer to [31]. The slopes $\sigma_{i}$ are computed from the neighbouring cells using a limiter function. To have a seconder order extension that preserves the positivity properties of the first order scheme, see Theorem 10, we choose the minmod limiter defined as

[TABLE]

and define the slopes as

[TABLE]

The same procedure is applied on $\psi_{i}^{(1)}$. Since we are using the minmod limiter on the conservative variables to determine the slope, we immediately get for $w_{i}^{n}\in\Omega_{phy}\cap\Omega_{wp}$ that for the interface values holds

[TABLE]

and the expansion of the interface values $\psi_{i-1/2}^{(1),+},\psi_{i+1/2}^{(1),-}$ with respect to the Mach number is preserved. This means by Lemma 9 that the Riemann problem applied on the interface values $w_{i+1/2}^{\mp},\psi_{i+1/2}^{\mp}$ still ensures the positivity of density and internal energy. By Theorem 10 the first order scheme has the positivity property and the second order scheme (6.2) is a convex combination of states in $\Omega_{phy}$. From this we have $w_{i}^{n+1}\in\Omega_{phy}$. We have thus proven the following result for the second order scheme:

Theorem 11** (Positivity property 2)**

Let the initial state be given as

[TABLE]

with the boundary condition described in (4.11). Then under the Mach number independent CFL condition $\frac{\Delta t}{\Delta x}\underset{i}{\max}|\lambda^{\pm}(w_{i}^{n})|<\frac{1}{2\cdot 2d}$, the numerical scheme defined by (6.1)/(6.2) preserves the positivity of density and internal energy, that is

[TABLE]

7 Numerical results

In the following section, we numerically validate the theoretical properties of the proposed scheme. Throughout the test cases, we assume an ideal gas law with the equation of state given by

[TABLE]

For solving the implicit non-symmetric linear system given by (5.1), we use the GMRES algorithm combined with a Jacobi diagonal scaling preconditioner for periodic boundary conditions used in the test cases in Section 7.2 and an preconditioner based on an incomplete LU decomposition for von Neumann boundary conditions used in the test cases from Section 7.1.

7.1 Shock test cases

To verify that our scheme captures discontinuities accurately, we perform a SOD shock tube test [32] in the regime $M=1$ and a Mach number dependent shock test case taken from [20].

7.1.1 SOD shock tube test

The computational domain for the SOD shock tube test [32] is $[0,1]$ and the initial data is given using $\gamma=1.4$ by

[TABLE]

where we place the initial discontinuity at $x=0.5$. Since the regime is compressible, we set $M=1$ and the initial data is given in dimensional form. This test also demonstrates that the collapse of the eigenvalues $\lambda^{\pm}$ and $\lambda_{M}^{\pm}$ in the case of $M=1$ as mentioned in Remark 8 is not problematic. We see in Figure 2 that the first order as well as the second order scheme captures the shock positions correctly. As expected, the second order scheme is more accurate than the first order scheme.

7.1.2 Mach number dependent shock problem

The following test case is a low Mach flow, where the contact wave travels with the Mach number $M=6.2\cdot 10^{-3}$. The initial data is taken from [20] and is given as

[TABLE]

The discontinuity is placed at $x_{0}=0.5$ on the domain $[0,1]$ with $\gamma=1.4$ and the final time is given by $T=0.25s$. To transform the initial data (7.4) into non-dimensional quantities, we define for a variable $\phi$ the relation $\phi=\phi_{r}\leavevmode\nobreak\ \hat{\phi}$, where $\phi$ denotes the dimensional variable, $\phi_{r}$ the reference value which contains the units and $\hat{\phi}$ the non-dimensional quantity. For the reference values we have the following relations

[TABLE]

We choose the scaling in space to be $x_{r}=1m$, the scaling in density to be $\rho_{r}=1\frac{kg}{m^{3}}$ and the scaling in velocity to be $u_{r}=M\frac{m}{s}$. This yields the following scaling of the sound speed $c_{r}=1\frac{m}{s}$, the time $t_{r}=Ms$ and the pressure $p_{r}=1\frac{kg}{ms^{2}}$. Then the non-dimensional initial data is given as

[TABLE]

In the simulation, we choose $M=6.2\cdot 10^{-3}$ which is the Mach number on the contact wave.

In Figure 3 we show the influence of the space and time step on the density profile computed with the first order scheme (IMEX1) and the second order scheme (IMEX2). The contact wave is always reproduced sharply whereas the acoustic waves are smoothed since the time step is chosen according to the CFL restriction (5.11) associated with the contact wave.

In Figure 4 our results computed with IMEX1 and IMEX2 are plotted against the implicit Jin Xin relaxation scheme presented in [20] and an explicit upwind scheme. Our results obtained with the IMEX schemes are in good agreement with the results of the implicit scheme. Since the time step for the explicit upwind scheme depends on the Mach number dependent acoustic waves, all waves are resolved. This is rather costly since it results in a very small time step. The implicit scheme is unconditionally stable and the time step can be chosen with respect to the desired accuracy of the numerical solution.

7.2 Gresho Vortex test

In order to demonstrate the low Mach properties, we calculate the solution to the Gresho Vortex test as given in [2]. The Gresho vortex is a steady state solution of the compressible Euler equations and the velocity field is divergence free.

The initial angular velocity in $\frac{m}{s}$ is given by

[TABLE]

where $r=\sqrt{(x-x_{0})^{2}+(y-y_{0})^{2}}$ with $x_{0}=0.5,y_{0}=0.5$ on a computational domain of $[0,1]$ and $\gamma=5/3$. The pressure distribution in $\frac{kg}{ms^{2}}$ is given by

[TABLE]

with $p_{0}=\frac{\rho_{0}u_{\Phi,\max}^{2}}{\gamma M^{2}}\frac{kg}{ms^{2}}$, where $\rho_{0}=1$. The initial density is given by $\rho=1\frac{kg}{m^{3}}$ and we transform the initial condition in non-dimensional quantities by using $x_{r}=1m$, $\rho_{r}=1\frac{kg}{m^{3}}$, $u_{r}=2\cdot 0.2\leavevmode\nobreak\ \pi\frac{m}{s}$, $p_{r}=\frac{\rho_{0}u_{0}^{2}}{\gamma M^{2}}\frac{kg}{ms^{2}}$ and $t_{r}=1\frac{m}{u_{r}}$. The computational domain is given by $[0,1]\times[0,1]$ and we use a $40\times 40$ grid with periodic boundary conditions. The results for a full turn of the vortex using the IMEX2 scheme together with the initial distribution of the Mach number are given in Figure 5. We see that even for a low resolution, the solution at $T=1$ shows only little dissipation throughout all tested Mach numbers. To further check the quality of the numerical simulation, we monitor the loss of kinetic energy. Since the Gresho vortex is a stationary solution, the kinetic energy should be preserved. In Figure 6 we show the ratio between the initial kinetic energy $E_{kin,0}$ and the kinetic energy after each time step $E_{kin,t}$ for the Mach numbers $M=10^{-2},10^{-3}$. The graphs for the different Mach numbers are indistinguishable which shows that the loss of kinetic energy does not depend on the chosen Mach number but depends on the chosen space discretization and time-step.

7.3 A smooth Gresho Vortex test

In this section we want to derive a test case to numerically validate the second order accuracy of the second order extension of Section 6. Unfortunately the Gresho vortex is not a smooth enough solution to test second order accuracy, since its velocity profile is only continuous but not continuously differentiable. Therefore we propose a smoothed velocity profile with which we calculate a pressure profile to gain a stationary vortex. A twice continuously differentiable angular velocity in $\frac{m}{s}$ with $u_{\Phi,\max}=1$ and $u_{\Phi}(0)=0$ and $u_{\Phi}(0.4)=0$ is given by

[TABLE]

with the radius $r=\sqrt{(x-0.5)^{2}+(y-0.5)^{2}}$. The profile can be easily modified to be a $C^{k}$ function, where $k\in\mathbb{N}$ denotes the degree of continuous differentiability. Under the condition that the centrifugal forces are balanced that is $\partial_{r}p=\frac{u_{\Phi}(r)^{2}}{r}$, we can calculate the pressure in $\frac{kg}{ms^{2}}$ as

[TABLE]

where

[TABLE]

As in the Gresho test case, the background pressure is scaled with the Mach number as

[TABLE]

To transform the dimensional data into non-dimensional quantities, we use the same reference values as in the Gresho vortex test case.

To show the accuracy of the IMEX2 scheme, we compute the solution of the smooth Gresho vortex test on the domain $[0,1]^{2}$ with periodic boundary conditions. In Table 1 the $L^{1}$-error between the numerical solution at $T=0.05$ and the initial configuration in non-dimensional quantities as well as the convergence rates are displayed for $M=10^{-1},10^{-2},10^{-3}$. It can be seen that we reach the expected accuracy independently of the chosen Mach number.

8 Conclusion and future developments

We have proposed an all-speed IMEX scheme for the Euler equations of gas dynamics which is based on Suliciu relaxation model. The proposed IMEX scheme is an improvement of explicit schemes since the time step is restricted by the material and not the acoustic wave speeds. It is also an improvement of implicit schemes since the implicit part consists only of one linear equation and can be solved very efficiently. The scheme has the correct numerical viscosity for all Mach numbers as shown by the Gresho vortex test case, can capture the correct shock positions as shown by the SOD shock test case as well as for Mach number dependent shock problems. In addition it is positivity preserving and shows the expected second order convergence rates.

In a future work, the scheme will be extended to the Euler equations with a gravitational source term. The aim is to design a multi-dimensional well-balanced scheme, that inherits all the properties of the IMEX scheme presented here.

Aknowledgements

G. Puppo ackwowledges the support by the GNCS-INDAM 2019 research project and A. Thomann the support of the INDAM-DP-COFUND-2015, grant number 713485. Some material of this work has been improved during the SHARK-FV conference (Sharing Higher-order Advanced Know-how on Finite Volume http://www.SHARK-FV.eu/) held in 2019.