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.
finite volume methods, Euler equations, positivity preserving, asymptotic preserving, relaxation, low Mach scheme, IMEX schemes
We consider the non-dimensional Euler equations in -space dimensions which are given by the following set of equations 
where the total energy is given by
and denotes the internal energy. The density is denoted by ,
is the velocity vector andis 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  for a study of the full Euler equations.
Standard schemes designed for compressible flows like the Roe scheme  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  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  through Jin Xin relaxation . 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 . 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  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  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 ,
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 . 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 . Following , in a first step the pressure is decomposed into a slow dynamics and a fast acoustic component
Approximating the slow and the fast pressure in the momentum equation in (1.1) by the variables and respectively, momentum equation becomes
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  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
In the equation for , we use instead of as proposed by the authors in . This is due to the upwind discretization used in  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.
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 .
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
where denotes the physical variables and the flux function is given by . The relaxation model (1.1) is given by
where denotes the state vector, the flux function as defined in (2.2) and the relaxation source term given by
The relaxation equilibrium state is defined as
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.
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
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
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
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
with a dynamical pressure described by .
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
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
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
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
Lemma 2 (Scaling of )
Let be equipped with the boundary conditions (4.11). Then the Mach number expansion of after the first stage satisfies
where is constant.
Proof. Let us assume that the expansion of is given by
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)
This leads to on . Separating the terms and using that , we find
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
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
Inserting the Mach number expansions of and into the density, momentum and energy equation of (3.7) and considering the order terms we have
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
In summary, we have shown the following theorem.
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
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
Using this notation, we apply centered differences for the implicit update (4.10) and have
For the explicit part, we will use a Godunov type finite volume scheme following  which we will describe in the following section.
5.1 Godunov type finite volume scheme
In the next result the structure of system (5.2) is summarized.
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
and with respect to
Proof. We rewrite the equations (5.2) using primitive variables in non-conservative form
where the matrix is given by
It is easy to check that are eigenvalues of
. Associated to the eigenvalues, we find the respective eigenvectors
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  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
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.
Consider an initial value problem with initial data given by
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
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
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
which is independent of the Mach number. This leads to the following update of the explicit part
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 .
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
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.
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
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
From Lemma 2, we know