# 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 accurate IMEX finite volume schemes for the compressible Euler equations in the zero Mach number limit is presented. It is well known that in the zero Mach number limit fast acoustic waves are filtered out, which results in the incompressible Euler equations containing only slow advection waves. In order to account for the fast and slow waves, the nonlinear fluxes in the Euler equations are split into stiff and non-stiff components, respectively. The time discretisation is performed by an IMEX Runge-Kutta method, therein the stiff terms are treated implicitly and the non-stiff terms explicitly. In the space discretisation, a Rusanov-type central flux is used for the non-stiff part, and simple central differencing for the stiff part. Both the time semi-discrete and space-time fully-discrete schemes are shown to be asymptotic preserving, i.e. both these schemes are consistent with the incompressible system in the zero Mach number limit, and their stability characteristics are independent of the Mach number. The numerical experiments confirm that the schemes achieve uniform second order convergence with respect to the Mach number. A notion of accuracy at low Mach numbers, termed as the asymptotic accuracy, is introduced in terms of the invariance of a well-prepared space of constant densities and divergence-free velocities. It is shown theoretically as well as numerically that the proposed schemes are asymptotically accurate.

## Authors

• 1 publication
• 2 publications
• ### Analysis of an Asymptotic Preserving Low Mach Number Accurate IMEX-RK Scheme for the Wave Equation System

In this paper the analysis of an asymptotic preserving (AP) IMEX-RK fini...
09/28/2019 ∙ by Arun K. R., et al. ∙ 0

• ### Asymptotic properties of a class of linearly implicit schemes for weakly compressible Euler equations

In this paper we derive and analyse a class of linearly implicit schemes...
11/27/2020 ∙ by Václav Kučera, et al. ∙ 0

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

We present an implicit-explicit finite volume scheme for the Euler equat...
07/19/2019 ∙ by Andrea Thomann, et al. ∙ 0

• ### High Order Semi-implicit WENO Schemes for All Mach Full Euler System of Gas Dynamics

In this paper, we propose high order semi-implicit schemes for the all M...
06/04/2021 ∙ by Sebastiano Boscarino, et al. ∙ 0

• ### Asymptotic structure of cosmological Burgers flows in one and two space dimensions: a numerical study

We study the cosmological Burgers model, as we call it, which is a nonli...
11/17/2019 ∙ by Yangyang Cao, et al. ∙ 0

• ### A novel structure preserving semi-implicit finite volume method for viscous and resistive magnetohydrodynamics

In this work we introduce a novel semi-implicit structure-preserving fin...
12/21/2020 ∙ by Francesco Fambri, et al. ∙ 0

• ### Numerical investigation of Mach number consistent Roe solvers for the Euler equations of gas dynamics

While traditional approaches to prevent the carbuncle phenomenon in gas ...
02/12/2021 ∙ by Friedemann Kemm, 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.

## 1. Introduction

Many physical phenomena in hydrodynamics and magnetohydrodynamics are governed by the compressible Euler equations which represent the fundamental conservation principles of mass, momentum and energy. In many problems from meteorology, geophysics, combustion etc., one often encounters a situation where a characteristic fluid velocity is much lesser than a corresponding sound velocity in the medium. In such cases, the ratio of these two velocities typically gives rise to a singular perturbation parameter, such as the Mach number or the Froude number. When the motion of the fluid is humble in comparison to the motion of sound waves in the medium, i.e. when the Mach number tends to zero, the fluid can be treated as almost incompressible. From a mathematical point of view, one can say that solutions of the compressible Euler equations are close to those of their incompressible counter parts as the Mach number approaches zero. In many seminal works, e.g. [18, 19, 28], a rigorous convergence analysis of the solutions of the purely hyperbolic compressible Euler system to those of the mixed hyperbolic-elliptic incompressible Euler system in the limit of zero Mach number can be found.

It is widely known in the literature that standard explicit Gudunov-type finite volume compressible flow solvers suffer from a lot pathologies at low Mach numbers. The stiffness arising from stringent CFL stability restrictions, loss of accuracy due to the creation of spurious waves, lack of stability due to the dependence of numerical viscosity on the Mach number, and inability to respect the transitional behaviour of the continuous system of governing equations are some of the commonly observed ailments, to name but a few. The stiffness due to CFL restrictions severely imposes the timesteps to be as small as the order of the Mach number, resulting in a drastic slowdown of the solver in a low Mach number regime. The inaccuracy of the numerical solution, and inability of the numerical method to respect the transitional behaviour of the equations are mainly due to the loss of information in the passage from continuous to discrete level. We refer the reader to [10, 13, 20], and the references therein for a detailed account of the above-mentioned anomalies.

Design of stable and accurate numerical schemes for weakly compressible or low Mach number flows is a challenging task due to the aforementioned difficulties and more. In the literature, several approaches to develop numerical schemes which work for small as well as large Mach numbers, addressing one or more of the above issues, can be found. In [4], Bijl and Wesseling used the standard MAC-type finite difference method for the incompressible Euler equations to derive a scheme which can simulate flows at a wide range of Mach numbers. Another approach is due to Munz et al. [23]

, in which the wellknown SIMPLE method for incompressible flows is extended to weakly compressible flows using the insight gained from an asymptotic analysis. Recently, Feistauer and Kučera

[12] developed an all Mach number flow solver in the context of high order discontinuous Galerkin methods. We refer the interested reader also to, e.g. [3, 21, 29], for a different approach using the so-called sound-proof models which eliminate sound waves completely and provide good approximations for low speed atmospheric flows.

The convergence of solutions of the compressible Euler equations to those of the incompressible Euler equations in the limit of zero Mach number, and the associated challenges in numerical approximation is an active area of research. The consistency of the limit of a compressible flow solver, if it exists, with the incompressible limit system, and the stiffness arising from stringent stability requirements are usually addressed in the framework of the famous asymptotic preserving (AP) methodology. The notion of AP schemes was initially introduced by Jin [15] for kinetic transport equations; see also [8, 16] for a comprehensive review of this subject. The AP framework provides a systematic and robust tool for developing numerical discretisation techniques for singular perturbation problems. The procedure takes into account the transitional behaviour of the governing equations of the problem, and its stability requirements are independent of the singular perturbation parameter. Hence, an AP scheme for the compressible Euler equations automatically transforms to an incompressible solver when the Mach number goes to zero. As discussed in [8], semi-implicit time-stepping techniques provide a systematic approach to derive AP schemes; see, e.g. [5, 7, 9, 24, 30], for some semi-implicit AP schemes for the Euler or shallow water equations.

The presence of singular perturbation parameters gives rise to multiple scales in time as well as space. As mentioned above, a widely used strategy to resolve multiple timescales, and to derive an AP scheme is the use of semi-implicit time discretisations. The property of being AP of a numerical approximation is primarily that of the particular time discretisation used. Implicit-explicit Runge-Kutta (IMEX-RK) schemes [14]

offer a precise and robust approach to define high order semi-implicit AP schemes for singularly perturbed ordinary differential equations. These schemes have been extensively used and analysed, e.g. in

[1, 25, 26]

for stiff systems of ordinary differential equations, and also time-dependent partial differential equations. The decisive step in the implementation of an IMEX-RK method to a system of conservation laws, such as the Euler or shallow water equations, is a splitting of the fluxes into the so-called stiff and non-stiff components in order that the resulting scheme is AP.

Recently, in [10, 11], the authors have presented the results of a detailed study on the accuracy of Godunov-type schemes at low Mach numbers. It has been shown in these papers that the origin of inaccuracies can be attributed to several factors: the number of space dimensions, particular discretisation strategies and numerical diffusion in the scheme used, cell geometry chosen, and so on. Based on the work of Schochet [28], Dellacherie in [10]

has presented the results of an extensive study on the accuracy of Gudunov-type schemes at low Mach numbers on Cartesian meshes. At the core of this analysis is an estimate from

[28] which states that when the Mach number is small, a solution to the compressible Euler equations is close to that of the incompressible equations whenever the corresponding initial data are close, cf. also [18, 19]. Further, an analysis of the nature of solutions of the linear wave equation system at low Mach numbers using Hodge decompositions showed that a sufficient condition to ensure the above-mentioned estimate in the linear case is the invariance of a so-called ‘well-prepared space’ of constant densities and divergence-free velocities. In the light of this invariance property, Dellacherie established that the origin of inaccuracies of standard Godunov-type schemes in two and three-dimensional geometries using the notion of first-order modified equations. Finally, a sufficient condition to avoid the generation of spurious waves, and to ensure accuracy at low Mach numbers for nonlinear schemes for the Euler equations is proposed by modifying the numerical viscosity, and changing the discretisation of the pressure gradient term.

In this paper, we study a class of second order accurate finite volume schemes for the isentropic Euler equations, employing the IMEX-RK time-stepping procedure. We prove that both the time semi-discrete and space-time fully discrete schemes are asymptotically consistent with the limiting incompressible system in the zero Mach number limit, and that their stability constraints are independent of the Mach number which implies the AP property. We also define a notion of accuracy at low Mach numbers, hereafter designated as the asymptotic accuracy, in terms of the invariance of the well-prepared space as in [10]. It is shown theoretically as well as numerically that both the time semi-discrete and space-time fully discrete schemes are asymptotically accurate. The rest of this paper is organised as follows. In Section 2 we mention the zero Mach number limits of the isentropic Euler and the linear wave equation systems, and recall the main convergence results from [10]. Section 3 is dedicated to the defintion of an AP scheme which is also accurate low Mach numbers. In Section 4 we introduce the IMEX-RK time semi-discrete schemes, and in Section 5 we show their AP property and asymptotic accuracy. Space-time fully-discrete schemes, and their analysis are presented in Section 6. In Section 7 numerical evidences are provided to support the theoretical claims made. Finally, we close the paper with some concluding remarks in Section 8.

## 2. Zero Mach Number Limits of the Euler and the Wave Equation Systems

### 2.1. Zero Mach Number Limit of the Euler System

The scaled, compressible Euler equations are given by

 (2.1) ∂tρ+∇⋅(ρu) =0, ∂t(ρu)+∇⋅(ρu⊗u)+∇pε2 =0,

where the independent variables are time and space , and the dependent variables are the density and velocity of the fluid . In order to close the system (2.1), an isentropic equation of state , where is a positive constant, is assumed. Here, the parameter is the reference Mach number defined by

 (2.2) ε:=urefcref,

with and being a reference fluid speed, and a reference sound speed, respectively.

In [18], Klainerman and Majda rigorously investigated the limit of solutions of (2.1) as . Formally, the asymptotic nature of solutions to (2.1) can be studied by plugging-in the ansatz

 (2.3) f(t,x)=f(0)(t,x)+εf(1)(t,x)+ε2f(2)(t,x)

for all dependent variables. Doing a scale analysis, and using appropriate boundary conditions, cf. e.g. [22], yields the multiscale representation

 (2.4) ρ(t,x)=ρ(0)(t)+ε2ρ(2)(t,x),

for the density and the incompressible Euler system:

 (2.5) =0, ∇⋅u(0) =0

for the velocity . Here, the second order pressure survives as the incompressible pressure.

Loosely speaking, the results of [18, 19] shows that when the initial data are almost incompressible, the solutions of the compressible Euler equations (2.1) are good approximations of those of the incompressible Euler equations (2.5). However, in order to make a precise statement, we recall the following convergence result due to Schochet [28]. Note that the non-dimensionalised, isentropic Euler system (2.1) can be recast in a non-conservative, evolution form as

 (2.6) ∂tU+H(U)+1ε2L(U)=0,

where

 (2.7) U:=(ρu), H(U):=(u⋅∇ρ(u⋅∇)u), L(U):=(ε2ρ∇⋅up′(ρ)∇ρ).

Here, is the convective operator with a time scale of order and is the acoustic operator with a time scale of the order of . The system (2.6) is supplied with periodic boundary conditions

 (2.8) U(0,x)=U0(x), x∈Td,

where denotes the -dimensional torus.

In [10, 28], a solution of (2.1) is sought in the space . Let us define the subspaces

 (2.9) E :={U∈L2(Td)1+d:∇ρ=0, and ∇⋅u=0}, (2.10) ~E :={U∈L2(Td)1+d:∫Tdρdx=0, and ∇×u=0}.

Note that is the subspace of spatially constant densities and divergence-free velocities. In other words, it is the ‘incompressible’ subspace and hence, hereafter, is referred to as the ‘well-prepared’ space. Having defined and , the following Helmholtz-Leray decomposition can be immediately written.

 (2.11) E⊕~E=L2(Td)1+d, and  E⊥~E.

Thus, for any , there exists a unique and , such that . Let us also define the Helmholtz-Leray projection of any onto the space as

 (2.12) PU:=^U.

In [28], Schochet has proved the following crucial estimate for the solutions of (2.1) in the limit , and we restate this important theorem as done in [10].

###### Theorem 2.1 ([10, 28]).

Let be a solution of the initial value problem for the system (2.6), i.e.

 (2.13) ∂tU+H(U)+1ε2L(U) =0, t>0, x∈Td, (2.14) U(0,x) =U0(x), x∈Td,

and let be a solution of the projected subsystem:

 (2.15) ∂t¯U+PH(¯U) =0, t>0, x∈Td, (2.16) ¯U(0,x) =^U0(x), x∈Td,

where . Then, for , there holds the estimate:

 (2.17) ∥ρ0−^ρ0∥=O(ε2), ∥u0−^u0∥=O(ε)⟹∥ρ(t)−¯ρ(t)∥=O(ε2), ∥u(t)−¯u(t)∥=O(ε) for all t>0.
###### Remark 2.2.

It is easy to see that the projected subsystem (2.15)-(2.16) is nothing but the incompressible Euler system (2.5). The essence of the above theorem is that when the initial data are well prepared, i.e. they are taken in , solutions of the compressible Euler system (2.1) approximate those of the incompressible Euler system (2.5) when . The theorem, clearly, is in agreement with the results obtained in [18, 19].

### 2.2. Zero Mach Number Limit of the Linear Wave Equation System

The rest of this section is devoted to the analysis of the linear wave equation system as a prototype model of the compressible Euler equations (2.6). Linearising the system of equations (2.6) about a constant state , we get the following wave equation system with advection:

 (2.18) ∂tU+H(U)+1ε2L(U)=0,

where

 (2.19) U:=(ρu), H(U):=(u––⋅∇ρ(u––⋅∇)u), L(U):=⎛⎜⎝ε2ρ––∇⋅ua––2ρ––∇ρ⎞⎟⎠.

In order to present an analysis of the solutions of (2.18), we first introduce a few definitions. For any and in , we consider the innerproduct

 (2.20) (U1,U2):=a––2ρ––ε2⟨ρ1,ρ2⟩+ρ––⟨u1,u2⟩,

where the innerproducts appearing on the right hand side are the usual innerproducts in . The norm generated by the above inner product is given by

 (2.21)

Clearly, the norm is equivalent to the usual norm in . We also consider the energies

 (2.22) E:=|||U|||2, Ein:=∣∣∣∣∣∣¯U∣∣∣∣∣∣2, Eac:=∣∣∣∣∣∣~U∣∣∣∣∣∣2,

which are, respectively, the total energy, the incompressible energy, and the acoustic or compressible energy. Since , we have . In the following, we state the corresponding linearised version of Schochet’s result, i.e. Theorem 2.1.

###### Theorem 2.3 ([10]).

Let be a solution of the IVP for the wave equation system (2.18), i.e.

 (2.23) ∂tU+H(U)+1ε2L(U) =0, t>0, x∈Td, (2.24) U(0,x) =U0(x), x∈Td,

and let be a solution of the IVP:

 (2.25) ∂t¯U+H(¯U) =0, t>0, x∈Td, (2.26) ¯U(0,x) =^U0(x), x∈Td,

where . Let be the Helmholtz-Leray decomposition of . Then the following holds.

1. ,

2. is the solution of (2.23) with initial condition .

Moreover, there holds the energy conservation:

 (2.27) Ein(t)=Ein(0) and Eac(t)=Eac(0), for all t>0,

and as a consequence, the following estimate holds:

 (2.28) ∥ρ0−^ρ0∥=O(ε2), ∥u0−^u0∥=O(ε)⟹∥ρ(t)−¯ρ(t)∥=O(ε2), ∥u(t)−¯u(t)∥=O(ε) for all t>0.
###### Remark 2.4.

Carrying out an analogous asymptotic analysis using the ansatz (2.3) for the wave equation yields the same multiscale form (2.4) for the density and the system

 (2.29) ∂tu(0)+(u––⋅∇)u(0)+a––2ρ––∇ρ(2) =0, (2.30) ∇⋅u(0) =0

for the unknowns . Since the boundary conditions are periodic, and (2.29) is linear, the divergence-free condition (2.30) on now forces to be a constant. Hence, , and (2.29) then reduces to a linear transport equation for , and the divergence-free condition (2.30) is now required only at time . Thus, the projected subsystem (2.25)-(2.26) and the system (2.29)-(2.30) are equivalent, and are incompressible.

The above theorem lies at the centre of the analysis of numerical schemes presented in [10]. However, numerical discretisation procedures introduce numerical diffusion, dispersion or higher order derivative terms, depending on the order of the scheme under consideration. Hence, the Theorem 2.3 has to be generalised to accomodate more general spatial differential operators arising from the discretisation. Such a generalisation is presented in [10] which also gives a sufficient condition to ensure the estimate (2.28).

###### Theorem 2.5 ([10]).

Let be a solution of the IVP:

 (2.31) ∂tU+FxU =0, t>0, x∈Td, (2.32) U(0,x) =U0(x), x∈Td,

which is assumed to be well-posed in , with a linear spatial differential operator. Then the following conclusions hold.

1. The solution satisfies the estimate

 (2.33) ∥U0−PU0∥=O(ε)⟹∥U(t)−¯U(t)∥=O(ε), for all t>0,

where is a solution of (2.31) with the initial condition . However, we don’t have the apriori estimate for all .

2. When the operator leaves invariant, i.e. whenever implies for all , then satisfies the estimate (2.33), and in addition we have

 (2.34) ∥U0−PU0∥=O(ε)⟹∥U(t)−PU(t)∥=O(ε) for all t>0.
###### Remark 2.6.

Note that in Theorem 2.1 and 2.3 the estimate on the difference is , whereas in [10] a scaled variable is used instead of , and hence the corresponding estimate is .

## 3. Asymptotic Preserving and Asymptotically Accurate Schemes

As mentioned in the introduction, the AP property must be necessary for any numerical scheme to survive in the passage of the limit . However, the accuracy of the scheme also deteriorates in the low Mach number limit, and hence it is essential for the scheme to maintain its accuracy. The two primary aims of the present work are to design and test an IMEX-RK finite volume scheme which is AP, and to maintain the accuracy in the limit . In what follows, we formally state the definitions of the AP property and asymptotic accuracy.

### 3.1. AP Property

Numerical solution of singular perturbation problems is, in general, a difficult task mainly because of the existence of multiple scales, the dependence of the stability characteristics on the perturbation parameter, and the transitional nature of the governing equations. A numerical approximation scheme may not resolve the scales consistently, its stability requirements might deteriorate, and in the singular limit the scheme might approximate a completely different set of equations than the actual limiting system. AP schemes provide a robust approach as a remedy against these problems. AP schemes were originally developed in the context of kinetic transport equations; see [8, 15, 16], and the references therein for more details.

###### Definition 3.1.

Let denote a singularly perturbed problem with the perturbation parameter . Let denote the limiting system of for . A discretisation of , with being the discretisation parameter, is called AP, if

1. is a consistent discretisation of the problem , called the asymptotic consistency, and

2. the stability constraints on are independent of , called the asymptotic stability.

In other words, the following diagram commutes:

 Pεhh→0−−−−→Pε⏐⏐↓ε→0⏐⏐↓ε→0P0hh→0−−−−→P0

### 3.2. Asymptotic Accuracy

It has to be noted that the Theorem 2.5, and its conclusions are valid for a linear scheme for a linear system of equations. Nonetheless, motivated by this result, and analogous considerations from [10], we propose the asymptotic accuracy of a scheme to be its ability to preserve the estimate (i) in Theorem 2.5. In other words, a numerical scheme for the compressible Euler system (2.1) should possess the ability to maintain solutions close to the incompressible solutions in for all times, whenever the initial data are well-prepared. In [10], such a notion was designated to be a scheme being accurate at low Mach numbers. Since the invariance of the well prepared space is a sufficient condition to maintain the estimate in the linear case, we make the following definition.

###### Definition 3.2.

A numerical approximation for the compressible Euler system (2.1) is said to be asymptotically accurate (AA), if it leaves the incompressible subspace invariant.

###### Remark 3.3.

It has to be noted both the AP and AA properties are to be satisfied by the time semi-discrete as well as space-time fully-discrete schemes. In the following sections we establish that the IMEX-RK schemes considered in this paper possess both the AP and AA properties.

## 4. Time Semi-discrete Scheme

### 4.1. Implicit-Explicit (IMEX) Runge-Kutta (RK) Time Discretisation

The IMEX-RK schemes are designed for the numerical integration of stiff systems of ordinary differential equations (ODEs) of the form

 (4.1) y′=f(t,y)+1εg(t,y),

where , and is usually known as the stiffness parameter. The functions and are called, respectively, the non-stiff part and the stiff part of the system (4.1). Stiff systems of ODEs of the type (4.1) typically arise in the modelling of semiconductor devices, study of kinetic equations, theory of hyperbolic systems with relaxation etc.; see, e.g. [14], for a comprehensive treatment of such systems.

Let be a numerical solution of (4.1) at time and let denote a fixed timestep. An -stage IMEX-RK scheme, cf., e.g. [1, 25], updates to through intermediate stages:

 (4.2) Yi =yn+Δti−1∑j=1~ai,jf(tn+~cjΔt,Yj)+Δts∑j=1ai,j1εg(tn+cjΔt,Yj), 1≤i≤s, (4.3) yn+1 =yn+Δts∑i=1~ωif(tn+~ciΔt,Yi)+Δts∑i=1ωi1εg(tn+ciΔt,Yi).

Let us denote , , and . The above IMEX-RK scheme (4.2)-(4.3) can be symbolically represented by the double Butcher tableau:

The matrices and are matrices such that resulting scheme is explicit in and implicit in . However, in order to reduce the number of implicit evaluations in the intermediate stages (4.2), we consider only the so-called diagonally implicit Runge-Kutta (DIRK) schemes in which for , and for . The coefficients and and the weights and are fixed by the order conditions; see [25, 26] for details. For the sake of completion, and to draw reference for the analysis carried out later, we write down the order conditions for a first order and second order IMEX-RK scheme as follows:

 (4.4) ~ci=i−1∑j=1~ai,j, ci=i∑j=1ai,j, (4.5) s∑i=1~ωi=1, s∑i=1ωi=1, (4.6) s∑i=1~ωi~ci=12, s∑i=1ωici=12, s∑i=1~ωici=12, s∑i=1ωi~ci=12.

The conditions (4.4), (4.5) and (4.6) are, respectively, the consistency, the first order, and the second order order conditions. To further simplify the analysis of the schemes presented in this paper, we restrict ourselves only to two types of DIRK schemes, namely the type-A and type-CK schemes which are defined below; see [17, 26] for details.

###### Definition 4.1.

An IMEX-RK scheme with the Butcher tableau given in Figure 1 is said to be of

• type-A, if the matrix is invertible;

• type-CK, if the matrix , can be written as

 A=(00αAs−1),

where and is invertible.

In the results presented in the later sections, we have used both the first order Euler(1,1,1), and second order ARS(2,2,2) schemes for time discretisations; see Figure 2 for their Butcher tableaus. Here, in the triplet , is the number of stages of the implicit part, the number gives the number of stages for the explicit part and gives the overall order of the scheme. For a detailed study of IMEX-RK schemes we refer the reader to [14, 17, 25, 26] and the references therein.

### 4.2. IMEX-RK Time Discretisation of the Euler System

Based on the asymptotic analysis and the convergence results presented in Section 2, we can split the flux functions in the Euler system (2.1) into a stiff and a non-stiff part, via

 (4.7) G(U):=(qp(ρ)ε2), F(U):=(0q⊗qρ),

where is the momentum. Applying IMEX-RK time discretisation, i.e. treating explicitly and implicitly, results in the following semi-discrete scheme.

###### Definition 4.2.

The stage of an -stage IMEX-RK scheme for the Euler system (2.1) is defined as

 (4.8) ρk =ρn−Δtak,l∇⋅ql, (4.9) qk =qn−Δt~ak,ℓ∇⋅(qℓ⊗qℓρℓ)−Δtak,l∇p(ρl)ε2, k=1,…,s.

The numerical solution at time is given by

 (4.10) ρn+1 =ρn−Δtωk∇⋅qk, (4.11) qn+1 =qn−Δt~ωk∇⋅(qk⊗qkρk)−Δtωk∇p(ρk)ε2.

In the above, and throughout the rest of this paper, we follow the convention that a repeated index always denotes the summation with respect to that index. The indices and are used to denote, respectively, the summation in the explicit and implicit terms, i.e. they assume values in the sets and .

Though the intermediate stages (4.8)-(4.9) consist of two implicit steps, the resolution of the scheme (4.8)-(4.9) is quite easy. As in [9], we eliminate between (4.8) and (4.9), and obtain the nonlinear elliptic equation:

 (4.12) −Δt2a2k,kε2Δp(ρk)+ρk=^ρk−Δtak,k∇⋅^qk

for . Here, we have denoted by

 (4.13) ^ρk :=ρn−Δtak,ℓ∇⋅qℓ, (4.14) ^qk :=qn−Δt~ak,ℓ∇⋅(qℓ⊗qℓρℓ)−Δtak,ℓε2∇p(ρℓ)

those terms that can be explicitly evaluated. Once is known after solving the elliptic equation (4.12), can be explicitly evaluated from (4.9). Finally, the updates and are calculated explicitly using (4.10)-(4.11) with the values obtained from the intermediate stages.

###### Remark 4.3.

In order to further reduce the nonlinear nature of the elliptic equation (4.12), and to make the numerical implementation simpler, in our computations we transform the nonlinear elliptic equation for into a semilinear elliptic equation for the pressure .

## 5. Analysis of Time Semi-discrete Scheme

The goal of this section is to establish the AP property and asymptotic accuracy of the time semi-discrete IMEX-RK scheme (4.8)-(4.11) in sense of Definitions 3.1 and 3.2.

### 5.1. AP Property

Establishing the AP property consists of showing the consistency of the scheme with that of the incompressible system as with its stability requirements independent of . In order to show the former, we perform an asymptotic analysis of the time discrete scheme (4.8)-(4.11), and to show the latter, we use a linear stability analysis.

#### 5.1.1. Asymptotic Consistency

###### Theorem 5.1.

Suppose that the data at time are well-prepared, i.e. and have the form:

 (5.1) ρn =ρn(0)+ε2ρn(2), (5.2) un =un(0)+εun(1),

where and . Then at the intermediate stages, and defined by (4.8) and (4.9) satisfy and , i.e. they are well-prepared as well, and as a consequence of this the numerical solutions and are well-prepared. In other words, the semi-discrete scheme (4.8)-(4.11) is asymptotically consistent with the incompressible limit system in the sense of Definition 3.1.

###### Proof.

We use induction on the number of stages to prove the theorem. To this end, we prove the consistency for the first stage, i.e. , which corresponds to a fully implicit step. We plugin the ansatz (2.3) for each of the dependent variables in (4.8)-(4.11). Equating to zero the terms in system (4.8)-(4.9) yields

 (5.3) ∇p(ρ1(0))=0.

Hence, the zeroth order density is spatially constant. In an analogous way we can show that the first order density is also spatially constant.

Next, we consider the terms in the mass update (4.8) to obtain

 (5.4) ρ1(0)−ρn(0)Δt=−a11ρ1(0)∇⋅u1(0).

We Integrate the above equation (5.4) over a spatial domain and use the Gauss’ divergence theorem to get

 (5.5) |Ω|ρ1(0)−ρn(0)Δt=−a11ρ1(0)∫Ω∇⋅u1(0)dx=−a11ρ1(0)∫∂Ων⋅u1(0)dσ.

If the boundary conditions of the problem are periodic or wall, then the right most integral in the above equation (5.5) vanishes; yielding . As the zeroth order density is a constant, as well as are constants. Using this in equation (5.4), we obtain the divergence condition:

 (5.6) ∇⋅u1(0)=0

for the leading order velocity . This completes the proof for .

To prove the result for onwards, we rewrite the intermediate stages (4.8)-(4.9) in the form

 (5.7) ρk =^ρk−Δtak,k∇⋅qk, (5.8) qk =^qk−Δtak,kε2∇p(ρk).

Now proceeding as in the case of , we can obtain . Hence, for the stage the zeroth order density is same as that of , and in turn we also obtain the divergence condition:

 (5.9) ∇⋅uk(0)=0

for .

Finally, to prove that the numerical solution at time is also well prepared, we consider the terms in the elliptic equation (4.12) and get

 (5.10) −Δt2a2k,kΔpk(2)+ρk(0)=^ρk(0)−Δtak,k∇⋅^qk(0).

As , and for all , we have

 (5.11) Δpk(2)=0.

Considering the terms in the final mass update (4.10) gives

 (5.12) ρn+1(0)=ρn(0)−Δtωk∇⋅qk(0).

Since and , we immediately get from (5.12) that

 (5.13) ρn+1(0)=ρn(0)=const.,

and similarly for . To get the divergence-free condition for the zeroth order velocity of the numerical solution, we equate to zero the terms in the final velocity update (4.11), and take its divergence to get

 (5.14) ∇⋅qn+1(0)=∇⋅qn(0)−Δt~ωk∇2:⎛⎝qk(0)⊗qk(0)ρk(0)⎞⎠−ΔtωkΔpk(2).

Since for the intermediate stages the zeroth order densities are constants, and velocities are divergence-free, using (5.11) in (5.14) yields the divergence-free condition:

 (5.15) ∇⋅un+1(0)=∇⋅un(0)=0.

Summarising the above steps, the asymptotic limit of scheme (4.8)-(4.11) is given by

 (5.16) ρn+1(0) =const., un+1(0) =un(0)−Δt~ωk∇⋅(uk(0)⊗uk(0))−Δtωk∇pk(2), ∇⋅un+1(0) =0.

Clearly, (5.16) is a consistent discretisation of the incompressible limit system (2.5). Hence, the semi-discrete scheme (4.8)-(4.11) is asymptotically consistent. ∎

###### Remark 5.2.

Note that the above proof is valid for any -stage IMEX-RK scheme. However, it uses the fact that for , which is the property of an RK scheme of type-A. The result also holds for a type-CK scheme in which the first step is trivial, and then on the proof follows similar lines as that of a type-A scheme.

#### 5.1.2. Linearised L2-stability and E-invariance

In this section, we analyse the correction terms arising from the time discretisations, and their effect on the asymptotic accuracy and stability of the resulting scheme. In order to analyse these correction terms, we follow the standard modified partial differential equation approach; see also [13] for a related analysis on stability. It has been shown in [10] that explicit Godunov-type schemes in a low Mach number regime are accurate only in one dimension, and that they are inaccurate in two and three dimensions. Specifically, the analysis in [10] shows that there is an acoustic time scale so that a Godunov-type scheme suffers from inaccuracies in a time interval of . A cure proposed in [10] is to remove the numerical diffusion terms from the momentum updates, and to use central differencing for the pressure gradient term. However, diffusion terms are required due to stability constraints and deleting them may not be a feasible option always; see also [2] for a related discussion. In the following, we show that using IMEX-RK time discretisation has the following advantages: no inaccuracy arises in two and three dimensions, it is not necessary to delete the diffusion coefficients to maintain stability, and the time-steps are independent of .

In order to establish the asymptotic stability of the IMEX-RK time discretisation for the isentropic Euler system (2.1), we carry out a thorough linear stability analysis for the linear wave equation system (2.18) as a linearised model. We believe that the results of linear stability analysis holds good also for the nonlinear Euler system (2.1). Analogous to Section 4, an IMEX-RK time discrete scheme for the linear wave equation system (2.18) is defined as follows.

###### Definition 5.3.

The stage of an -stage IMEX-RK scheme for the wave equation system (2.18) is given by

 (5.17) ρk =ρn−Δt~ak,ℓ(u––⋅∇)ρℓ−Δtak,lρ––∇⋅ul, (5.18) uk =un−Δt~ak,ℓ(u––⋅∇)uℓ−Δtak,la––2ρ––ε2∇ρl.

The numerical solutions and at time are defined as

 (5.19) ρn+1 =ρn−Δt~ωk(u––⋅∇)ρk−Δtωkρ––∇⋅uk, (5.20) un+1 =un−Δt~ωk(u––⋅∇)uk−Δtωka––2ρ––ε2∇ρk.

In the following, we derive the modified partial differential equation (MPDE) corresponding to a general second order accurate time discrete scheme of the form (5.17)-(5.20). We show that the solution of the MPDE is energy-dissipative, and the time-steps are independent of . Hence, the time discrete scheme (4.8)-(4.11) for the Euler system is linearly asymptotically -stable. In addition, we also prove that the MPDE keeps the well-prepared space invariant.

###### Theorem 5.4.

The time-discrete IMEX-RK scheme (4.8)-(4.11) is linearly -stable as long as the timesteps are bounded.

###### Proof.

We consider a general second order accurate IMEX-RK time discretisation in (5.17)-(5.20). Expanding the unknown functions in a Taylor series, and making use of the conditions (4.4)-(4.6) for the IMEX-RK coefficients, results in the following MPDE:

 (5.21)

Here, the operators and contain the third and fourth derivatives of the unknown functions and , respectively. Since these expressions are quite long, we present them only in the Appendix. After using the second-order order conditions (4.4)-(4.6), the MPDE (5.21) is free from first and second order derivative terms.

The rate of change of the energy , defined in (2.22), is given by

 (5.22) dEdt=2(U,∂tU)=2a––2ρ––ε2⟨ρ,∂tρ⟩+2ρ––⟨u,∂tu⟩.

We use the MPDE (5.21) in (5.22), and integrate the resulting terms by parts. It is easy to see that the third order derivatives do not contribute as they all vanish in view of periodic boundary conditions, and hence we get

 (5.23)
 (5.24) ⟨u,∂tu⟩=Δt3{b(4)1a––2ε2⟨u,(u––⋅∇)2∇(∇⋅u)⟩+b(4)2a––2ρ––ε2⟨u,(u––⋅∇)3∇ρ⟩+ b(4)3a––4ρ––ε4⟨u,(u––⋅∇)∇Δρ⟩+b(4)4a––4ε4⟨u,∇Δ(∇⋅u)⟩−124⟨u,(u––⋅∇)4u⟩}.

Now using the Cauchy-Schwarz inequality in (5.23) and (5.24), and using the inequalities thus obtained in (5.22) we finally obtain

 (5.25) dEdt≤2Δt3[b(4)1a––4ρ––ε4∥(u––⋅∇)∇ρ∥2+b(4)2a––2ε2(∥(u––⋅∇)2ρ∥2+∥(u––⋅∇)∇⋅u∥2)+b(4)3a––4ε4(∥(u––⋅∇)∇⋅u∥2+∥Δρ∥2)+b(4)4a––6ρ––ε6∥Δρ∥2−a––224ρ––ε2∥(u––⋅∇)2ρ∥2+ b(4)1ρ––a––2ε2∥(u––⋅∇)∇⋅u∥2+b(4)4ρ––