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 hyperbolicelliptic incompressible Euler system in the limit of zero Mach number can be found.
It is widely known in the literature that standard explicit Gudunovtype 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 abovementioned 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 MACtype 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 socalled soundproof 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], semiimplicit timestepping techniques provide a systematic approach to derive AP schemes; see, e.g. [5, 7, 9, 24, 30], for some semiimplicit 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 semiimplicit time discretisations. The property of being AP of a numerical approximation is primarily that of the particular time discretisation used. Implicitexplicit RungeKutta (IMEXRK) schemes [14]
offer a precise and robust approach to define high order semiimplicit 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 timedependent partial differential equations. The decisive step in the implementation of an IMEXRK method to a system of conservation laws, such as the Euler or shallow water equations, is a splitting of the fluxes into the socalled stiff and nonstiff 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 Godunovtype 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 Gudunovtype 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 abovementioned estimate in the linear case is the invariance of a socalled ‘wellprepared space’ of constant densities and divergencefree velocities. In the light of this invariance property, Dellacherie established that the origin of inaccuracies of standard Godunovtype schemes in two and threedimensional geometries using the notion of firstorder 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 IMEXRK timestepping procedure. We prove that both the time semidiscrete and spacetime 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 wellprepared space as in [10]. It is shown theoretically as well as numerically that both the time semidiscrete and spacetime 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 IMEXRK time semidiscrete schemes, and in Section 5 we show their AP property and asymptotic accuracy. Spacetime fullydiscrete 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)  
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) 
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 pluggingin the ansatz
(2.3) 
for all dependent variables. Doing a scale analysis, and using appropriate boundary conditions, cf. e.g. [22], yields the multiscale representation
(2.4) 
for the density and the incompressible Euler system:
(2.5)  
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 nondimensionalised, isentropic Euler system (2.1) can be recast in a nonconservative, evolution form as
(2.6) 
where
(2.7) 
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) 
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)  
(2.10) 
Note that is the subspace of spatially constant densities and divergencefree velocities. In other words, it is the ‘incompressible’ subspace and hence, hereafter, is referred to as the ‘wellprepared’ space. Having defined and , the following HelmholtzLeray decomposition can be immediately written.
(2.11) 
Thus, for any , there exists a unique and , such that . Let us also define the HelmholtzLeray projection of any onto the space as
(2.12) 
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)  
(2.14) 
and let be a solution of the projected subsystem:
(2.15)  
(2.16) 
where . Then, for , there holds the estimate:
(2.17) 
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) 
where
(2.19) 
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) 
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) 
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)  
(2.24) 
and let be a solution of the IVP:
(2.25)  
(2.26) 
where . Let be the HelmholtzLeray decomposition of . Then the following holds.

,

is the solution of (2.23) with initial condition .
Moreover, there holds the energy conservation:
(2.27) 
and as a consequence, the following estimate holds:
(2.28) 
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)  
(2.30) 
for the unknowns . Since the boundary conditions are periodic, and (2.29) is linear, the divergencefree condition (2.30) on now forces to be a constant. Hence, , and (2.29) then reduces to a linear transport equation for , and the divergencefree 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)  
(2.32) 
which is assumed to be wellposed in , with a linear spatial differential operator. Then the following conclusions hold.
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 IMEXRK 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

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

the stability constraints on are independent of , called the asymptotic stability.
In other words, the following diagram commutes:
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 wellprepared. 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 semidiscrete as well as spacetime fullydiscrete schemes. In the following sections we establish that the IMEXRK schemes considered in this paper possess both the AP and AA properties.
4. Time Semidiscrete Scheme
4.1. ImplicitExplicit (IMEX) RungeKutta (RK) Time Discretisation
The IMEXRK schemes are designed for the numerical integration of stiff systems of ordinary differential equations (ODEs) of the form
(4.1) 
where , and is usually known as the stiffness parameter. The functions and are called, respectively, the nonstiff 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 IMEXRK scheme, cf., e.g. [1, 25], updates to through intermediate stages:
(4.2)  
(4.3) 
Let us denote , , and . The above IMEXRK 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 socalled diagonally implicit RungeKutta (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 IMEXRK scheme as follows:
(4.4)  
(4.5)  
(4.6) 
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 typeA and typeCK schemes which are defined below; see [17, 26] for details.
Definition 4.1.
An IMEXRK scheme with the Butcher tableau given in Figure 1 is said to be of

typeA, if the matrix is invertible;

typeCK, if the matrix , can be written as
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 IMEXRK schemes we refer the reader to [14, 17, 25, 26] and the references therein.
0  0  0 

1  1  0 
1  0 
0  0  0 

1  0  1 
0  1 
0  0  0  0 

0  0  
0  
0  0  
0  
4.2. IMEXRK 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 nonstiff part, via
(4.7) 
where is the momentum. Applying IMEXRK time discretisation, i.e. treating explicitly and implicitly, results in the following semidiscrete scheme.
Definition 4.2.
The stage of an stage IMEXRK scheme for the Euler system (2.1) is defined as
(4.8)  
(4.9) 
The numerical solution at time is given by
(4.10)  
(4.11) 
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) 
for . Here, we have denoted by
(4.13)  
(4.14) 
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 Semidiscrete Scheme
The goal of this section is to establish the AP property and asymptotic accuracy of the time semidiscrete IMEXRK 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 wellprepared, i.e. and have the form:
(5.1)  
(5.2) 
where and . Then at the intermediate stages, and defined by (4.8) and (4.9) satisfy and , i.e. they are wellprepared as well, and as a consequence of this the numerical solutions and are wellprepared. In other words, the semidiscrete 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) 
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) 
We Integrate the above equation (5.4) over a spatial domain and use the Gauss’ divergence theorem to get
(5.5) 
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) 
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)  
(5.8) 
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) 
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) 
As , and for all , we have
(5.11) 
Considering the terms in the final mass update (4.10) gives
(5.12) 
Since and , we immediately get from (5.12) that
(5.13) 
and similarly for . To get the divergencefree 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) 
Since for the intermediate stages the zeroth order densities are constants, and velocities are divergencefree, using (5.11) in (5.14) yields the divergencefree condition:
(5.15) 
Remark 5.2.
Note that the above proof is valid for any stage IMEXRK scheme. However, it uses the fact that for , which is the property of an RK scheme of typeA. The result also holds for a typeCK scheme in which the first step is trivial, and then on the proof follows similar lines as that of a typeA scheme.
5.1.2. Linearised stability and 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 Godunovtype 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 Godunovtype 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 IMEXRK 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 timesteps are independent of .
In order to establish the asymptotic stability of the IMEXRK 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 IMEXRK time discrete scheme for the linear wave equation system (2.18) is defined as follows.
Definition 5.3.
The stage of an stage IMEXRK scheme for the wave equation system (2.18) is given by
(5.17)  
(5.18) 
The numerical solutions and at time are defined as
(5.19)  
(5.20) 
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 energydissipative, and the timesteps 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 wellprepared space invariant.
Theorem 5.4.
Proof.
We consider a general second order accurate IMEXRK 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 IMEXRK 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 secondorder 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) 
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) 
Now using the CauchySchwarz inequality in (5.23) and (5.24), and using the inequalities thus obtained in (5.22) we finally obtain
(5.25) 
Comments
There are no comments yet.