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

TL;DR

This paper develops second order IMEX finite volume schemes for the Euler equations that are asymptotic preserving and accurate at low Mach numbers, effectively bridging compressible and incompressible flow regimes.

Contribution

It introduces a novel IMEX scheme that is both asymptotic preserving and asymptotically accurate for low Mach number flows, with rigorous analysis and validation.

Findings

Schemes achieve uniform second order convergence across Mach numbers.

Numerical experiments confirm asymptotic accuracy at low Mach regimes.

Theoretical analysis supports the asymptotic preserving property.

Abstract

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. 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. 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. The asymptotic accuracy is concerned with the closeness of the compressible solution with that of its incompressible counterpart in a low Mach number regime. It is shown theoretically as well as numerically that the proposed schemes are asymptotically accurate.

Tables

Table 1. Table 1. L 1 superscript 𝐿 1 L^{1} and L 2 superscript 𝐿 2 L^{2} errors in u 1 subscript 𝑢 1 u_{1} and u 2 subscript 𝑢 2 u_{2} and the EOC computed for Problem  7.1 . Top to bottom: ε = 10 − 6 𝜀 superscript 10 6 \varepsilon=10^{-6} to ε = 10 − 1 𝜀 superscript 10 1 \varepsilon=10^{-1} .

$N$	$L^1$ error in $u_1$	EOC	$L^2$ error in $u_1$	EOC	$L^1$ error in $u_2$	EOC	$L^2$ error in $u_2$	EOC
10	4.9120e-03		8.8088e-03		7.9419e-03		1.7131e-02
20	1.3454e-03	1.8682	2.5016e-03	1.8160	2.6670e-03	1.5742	5.7609e-03	1.5722
40	3.1818e-04	2.0800	6.5165e-04	1.9406	6.5744e-04	2.0203	1.4305e-03	2.0097
80	8.0467e-05	1.9834	1.8591e-04	1.8094	1.6067e-04	2.0327	3.6866e-04	1.9561
10	4.9125e-03		8.8099e-03		7.9430e-03		1.7131e-02
20	1.3445e-03	1.8693	2.5017e-03	1.8162	2.6657e-03	1.5751	5.7604e-03	1.5724
40	3.1752e-04	2.0822	6.5147e-04	1.9411	6.5705e-04	2.0205	1.4302e-03	2.0099
80	7.7912e-05	2.0270	1.8559e-04	1.8116	1.5926e-04	2.0446	3.6858e-04	1.9562
10	4.9125e-03		8.8099e-03		7.9430e-03		1.7131e-02
20	1.3445e-03	1.8693	2.5017e-03	1.8162	2.6657e-03	1.5751	5.7604e-03	1.5724
40	3.1752e-04	2.0822	6.5147e-04	1.9411	6.5705e-04	2.0205	1.4302e-03	2.0099
80	7.7911e-05	2.0270	1.8558e-04	1.8116	1.5926e-04	2.0446	3.6858e-04	1.9562
10	4.9125e-03		8.8099e-03		7.9430e-03		1.7131e-02
20	1.3445e-03	1.8693	2.5016e-03	1.8162	2.6657e-03	1.5751	5.7604e-03	1.5724
40	3.1740e-04	2.0827	6.5149e-04	1.9411	6.5702e-04	2.0205	1.4302e-03	2.0099
80	7.7824e-05	2.0280	1.8558e-04	1.8117	1.5923e-04	2.0448	3.6854e-04	1.9563
10	4.9117e-03		8.8096e-03		7.9411e-03		1.7130e-02
20	1.3433e-03	1.8704	2.5006e-03	1.8167	2.6647e-03	1.5753	5.7578e-03	1.5729
40	3.1573e-04	2.0890	6.5096e-04	1.9416	6.5528e-04	2.0237	1.4272e-03	2.0123
80	8.1065e-05	1.9615	1.8943e-04	1.7808	1.5865e-04	2.0462	3.6843e-04	1.9537
10	4.9009e-03		8.7985e-03		7.9211e-03		1.7111e-02
20	1.3395e-03	1.8713	2.5124e-03	1.8081	2.6297e-03	1.5907	5.7013e-03	1.5856
40	3.4569e-04	1.9541	7.1801e-04	1.8070	6.4199e-04	2.0342	1.4288e-03	1.9964
80	1.1985e-04	1.5282	2.6469e-04	1.4396	1.7083e-04	1.9099	3.9443e-04	1.8570

Table 2. Table 2. L 1 superscript 𝐿 1 L^{1} and L 2 superscript 𝐿 2 L^{2} errors in u 1 subscript 𝑢 1 u_{1} and u 2 subscript 𝑢 2 u_{2} , and the AOC for Problem  7.2 . Top to bottom: ε = 10 − 6 𝜀 superscript 10 6 \varepsilon=10^{-6} to ε = 10 − 4 𝜀 superscript 10 4 \varepsilon=10^{-4} .

$N$	$L^1$ error in $u_1$	AOC	$L^2$ error in $u_1$	AOC	$L^1$ error in $u_2$	AOC	$L^2$ error in $u_2$	AOC
10	7.5900e-01		9.3643e-01		7.5900e-01		9.3644e-01
20	2.7311e-01	1.4746	3.3303e-01	1.4915	2.7311e-01	1.4746	3.3303e-01	1.4915
40	6.4266e-02	2.0874	7.5131e-02	2.1482	6.4266e-02	2.0874	7.5131e-02	2.1482
80	1.5283e-02	2.0721	1.7248e-02	2.1230	1.5283e-02	2.0721	1.7248e-02	2.1230
10	7.5900e-01		9.3643e-01		7.5900e-01		9.3643e-01
20	2.7308e-01	1.4748	3.3299e-01	1.4917	2.7308e-01	1.4748	3.3299e-01	1.4917
40	6.4183e-02	2.0890	7.5043e-02	2.1497	6.4183e-02	2.0890	7.5043e-02	2.1497
80	1.5198e-02	2.0783	1.7154e-02	2.1292	1.5198e-02	2.0783	1.7154e-02	2.1292
10	7.5900e-01		9.3643e-01		7.5900e-01		9.3643e-01
20	2.7302e-01	1.4751	3.3292e-01	1.4920	2.7302e-01	1.4751	3.3292e-01	1.4920
40	6.3787e-02	2.0977	7.4589e-02	2.1582	6.3787e-02	2.0977	7.4589e-02	2.1582
80	1.3905e-02	2.1976	1.5690e-02	2.2490	1.3905e-02	2.1976	1.5690e-02	2.2490

Equations

$$\displaystyle\partial_{t}\rho+\nabla\cdot(\rho\mbox{\boldmath$u$})$$

$$\displaystyle\partial_{t}(\rho\mbox{\boldmath$u$})+\nabla\cdot(\rho\mbox{\boldmath$u$}\otimes\mbox{\boldmath$u$})+\frac{\nabla p}{\varepsilon^{2}}$$

$$\epsilon :=\frac{u_{\mathrm{ref}}}{c_{\mathrm{ref}}},$$

$$f(t,\mbox{\boldmath$x$})=f_{(0)}(t,\mbox{\boldmath$x$})+\varepsilon f_{(1)}(t,\mbox{\boldmath$x$})+\varepsilon^{2}f_{(2)}(t,\mbox{\boldmath$x$})$$

$$p(t,\mbox{\boldmath$x$})=p_{(0)}(t)+\varepsilon^{2}p_{(2)}(t,\mbox{\boldmath$x$}),$$

$$\displaystyle\partial_{t}\mbox{\boldmath$u$}_{(0)}+\nabla\cdot\left(\mbox{\boldmath$u$}_{(0)}\otimes\mbox{\boldmath$u$}_{(0)}\right)+\nabla p_{(2)}$$

$$\displaystyle\nabla\cdot\mbox{\boldmath$u$}_{(0)}$$

$${\partial }_{t}U+\mathcal{H}(U)+\frac{1}{{\epsilon }^2}\mathcal{L}(U)=0,$$

$$U:=\begin{pmatrix}\rho\\ \mbox{\boldmath$u$}\end{pmatrix},\ \mathcal{H}(U):=\begin{pmatrix}\mbox{\boldmath$u$}\cdot\nabla\rho\\ (\mbox{\boldmath$u$}\cdot\nabla)\mbox{\boldmath$u$}\end{pmatrix},\ \mathcal{L}(U):=\begin{pmatrix}\varepsilon^{2}\rho\nabla\cdot\mbox{\boldmath$u$}\\ \frac{p^{\prime}(\rho)}{\rho}\nabla\rho\end{pmatrix}.$$

$$U(0,\mbox{\boldmath$x$})=U_{0}(\mbox{\boldmath$x$}),\ \mbox{\boldmath$x$}\in\mathbb{T}^{d},$$

$$\mathcal{E}$$

$$\tilde{\mathcal{E}}$$

$$\mathcal{E}\oplus \tilde{\mathcal{E}}=L^2({\mathbb{T}}^d{)}^{1+d},\text{mbox}and\mathcal{E}\perp \tilde{\mathcal{E}}.$$

$$\mathbb{P}U:=\hat{U}.$$

$${\partial }_{t}U+\mathcal{H}(U)+\frac{1}{{\epsilon }^2}\mathcal{L}(U)$$

$$\displaystyle U(0,\mbox{\boldmath$x$})$$

$${\partial }_t\bar{U}+\mathbb{P}\mathcal{H}(\bar{U})$$

$$\displaystyle\bar{U}(0,\mbox{\boldmath$x$})$$

$$\lVert\rho_{0}-\hat{\rho}_{0}\rVert=\mathcal{O}(\varepsilon^{2}),\ \lVert\mbox{\boldmath$u$}_{0}-\hat{\mbox{\boldmath$u$}}_{0}\rVert=\mathcal{O}(\varepsilon)\implies\lVert\rho(t)-\bar{\rho}(t)\rVert=\mathcal{O}(\varepsilon^{2}),\ \lVert\mbox{\boldmath$u$}(t)-\bar{\mbox{\boldmath$u$}}(t)\rVert=\mathcal{O}(\varepsilon)\ \mbox{for all}\ t>0.$$

$${\partial }_{t}U+H(U)+\frac{1}{{\epsilon }^2}L(U)=0,$$

$$\displaystyle U:=\begin{pmatrix}\rho\\ \mbox{\boldmath$u$}\end{pmatrix},\ H(U):=\begin{pmatrix}\underline{\mbox{\boldmath$u$}}\cdot\nabla\rho\\ (\underline{\mbox{\boldmath$u$}}\cdot\nabla)\mbox{\boldmath$u$}\end{pmatrix},\ L(U):=\begin{pmatrix}\varepsilon^{2}\underline{\rho}\nabla\cdot\mbox{\boldmath$u$}\\ \frac{\underline{a}^{2}}{\underline{\rho}}\nabla\rho\end{pmatrix}.$$

$$(U_{1},U_{2}):=\frac{\underline{a}^{2}}{\underline{\rho}\varepsilon^{2}}\langle\rho_{1},\rho_{2}\rangle+\underline{\rho}\langle\mbox{\boldmath$u$}_{1},\mbox{\boldmath$u$}_{2}\rangle,$$

$${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|U\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}:=\left[\frac{\underline{a}^{2}}{\underline{\rho}\varepsilon^{2}}\lVert\rho\rVert^{2}+\underline{\rho}\lVert\mbox{\boldmath$u$}\rVert^{2}\right]^{\frac{1}{2}}.$$

$$E:={\left|\left|\left|U\right|\right|\right|}^2,E_{in}:={\left|\left|\left|\bar{U}\right|\right|\right|}^2,E_{ac}:={\left|\left|\left|\tilde{U}\right|\right|\right|}^2,$$

$${\partial }_{t}U+H(U)+\frac{1}{{\epsilon }^2}L(U)$$

$$\displaystyle U(0,\mbox{\boldmath$x$})$$

$${\partial }_t\bar{U}+H(\bar{U})$$

$$\displaystyle\bar{U}(0,\mbox{\boldmath$x$})$$

$$E_{in}(t)=E_{in}(0)\text{and}{E}_{ac}(t)=E_{ac}(0),\text{for all}t>0,$$

$$\lVert\rho_{0}-\hat{\rho}_{0}\rVert=\mathcal{O}(\varepsilon^{2}),\ \lVert\mbox{\boldmath$u$}_{0}-\hat{\mbox{\boldmath$u$}}_{0}\rVert=\mathcal{O}(\varepsilon)\implies\lVert\rho(t)-\bar{\rho}(t)\rVert=\mathcal{O}(\varepsilon^{2}),\ \lVert\mbox{\boldmath$u$}(t)-\bar{\mbox{\boldmath$u$}}(t)\rVert=\mathcal{O}(\varepsilon)\ \mbox{for all}\ t>0.$$

$$\displaystyle\partial_{t}\mbox{\boldmath$u$}_{(0)}+(\underline{\mbox{\boldmath$u$}}\cdot\nabla)\mbox{\boldmath$u$}_{(0)}+\frac{\underline{a}^{2}}{\underline{\rho}}\nabla\rho_{(2)}$$

$$\displaystyle\nabla\cdot\mbox{\boldmath$u$}_{(0)}$$

$$\displaystyle\partial_{t}U+\mathcal{F}_{\mbox{\boldmath$x$}}U$$

$$\displaystyle U(0,\mbox{\boldmath$x$})$$

$$\Vert U_0-\mathbb{P}{U}_0\Vert =\mathcal{O}(\epsilon )⟹\Vert U(t)-\bar{U}(t)\Vert =\mathcal{O}(\epsilon ),\text{for all}t>0,$$

Full text

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 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 counterparts as the Mach number approaches zero. In many seminal works, e.g. [25, 26, 35], 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 [14, 19, 27] 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 [5], 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. [30], 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 [18] 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. [4, 28, 36], 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 [22] for kinetic transport equations; see also [12, 23] 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 [12], semi-implicit time-stepping techniques provide a systematic approach to derive AP schemes; see, e.g. [6, 11, 13, 31, 38], 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 offer a precise and robust approach to define high order semi-implicit AP schemes for singularly perturbed ordinary differential equations (ODEs). In the literature, one can find several references where these schemes have been extensively used and analysed for stiff systems of ODEs, and also time-dependent partial differential equations. We refer the interested reader to [2, 32, 33] for the application of IMEX schemes to stiff ODEs, to [8] for an error analysis of IMEX schemes, and to [10, 11, 13, 16, 17, 20, 37] for an application of the IMEX strategy to low Mach number Euler equations. However, 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 [14, 15], 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 is predominantly due to the formation of spurious acoustic waves in the discrete solution. There are several factors which can attributed to this, such as 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 [35], Dellacherie in [14] 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 [35] which states that when the Mach number is small, a solution to the compressible Euler equations is “$\varepsilon$-close” to that of the incompressible equations whenever the corresponding initial data are close, where $\varepsilon$ is a reference Mach number, cf. also [25, 26]. Further, an analysis of the nature of solutions of the linear wave equation system at low Mach numbers using Hodge decompositions shows 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 identified 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. The flux-splitting required to design these IMEX schemes is analogous to that used in [10, 16, 37]. 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 [14]. In fact, we show that the invariance is the key property for the numerical solution to mimic the $\varepsilon$-closeness as in the case of the exact solution. As result of the $\varepsilon$-closeness, the error remains bounded, and we obtain the accuracy in the stiff limit $\varepsilon=0$, justifying the name asymptotic accuracy. 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 [14]. Section 3 is dedicated to the defintion of an AP and AA scheme. In Section 4, we introduce the IMEX-RK time semi-discrete schemes, and in Section 5 we show their AP property by giving a sufficient condition for $L^{2}$-stability, and the asymptotic accuracy. Space-time fully-discrete schemes, and their analysis containing a similar stability condition, 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

[TABLE]

where the independent variables are time $t>0$ and space \mbox{\boldmathx}\in\mathbb{R}^{d},d=1,2,3, and the dependent variables are the density \rho=\rho(t,\mbox{\boldmathx})>0 and velocity of the fluid \mbox{\boldmathu}=\mbox{\boldmathu}(t,\mbox{\boldmathx})\in\mathbb{R}^{d}. In order to close the system (2.1), an isentropic equation of state $p(\rho)=\rho^{\gamma}$, where $\gamma$ is a positive constant, is assumed. Here, the parameter $\varepsilon$ is a reference Mach number defined by

[TABLE]

with $u_{\mathrm{ref}}$ and $c_{\mathrm{ref}}$ being a reference fluid speed, and a reference sound speed, respectively.

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

[TABLE]

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

[TABLE]

for the pressure and the incompressible Euler system:

[TABLE]

for the velocity \mbox{\boldmathu}_{(0)}. Here, the second order pressure $p_{(2)}$ survives as the incompressible pressure.

Loosely speaking, the results of [25, 26] shows that when the initial data are almost incompressible, the solutions of the compressible Euler equations (2.1) are good approximations to those of the incompressible Euler equations (2.5). However, in order to make a more precise statement, and to carry out the analysis presented later, we recall the following convergence result due to Schochet [35], and some related results. Note that the non-dimensionalised, isentropic Euler system (2.1) can be recast into a non-conservative, evolution form as

[TABLE]

where

[TABLE]

Here, $\mathcal{H}$ is the convective operator with a time scale of order $1$ and $\mathcal{L}/\varepsilon^{2}$ is the acoustic operator with a time scale of the order of $\varepsilon$. The system (2.6) is supplied with periodic boundary conditions

[TABLE]

where $\mathbb{T}^{d}$ denotes the $d$-dimensional torus.

Following [14, 35], in order to study the solutions U=(\rho,\mbox{\boldmathu})^{T} of (2.1) in the space $L^{2}(\mathbb{T}^{d})^{1+d}$, we consider the subspaces

[TABLE]

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

[TABLE]

Thus, for any $U\in L^{2}(\mathbb{T}^{d})^{1+d}$, there exists a unique $\hat{U}\in\mathcal{E}$ and $\tilde{U}\in\tilde{\mathcal{E}}$, such that $U=\hat{U}+\tilde{U}$. Let us also define the Helmholtz-Leray projection of any $U\in L^{2}(\mathbb{T}^{d})^{1+d}$ onto the space $\mathcal{E}$ as

[TABLE]

In [35], Schochet proved the following crucial estimate for the solutions of (2.1) in the limit $\varepsilon\to 0$, and we restate this important theorem as done in [14].

Theorem 2.1** ([14, 35]).**

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

[TABLE]

and let $\bar{U}$ be a solution of the projected subsystem:

[TABLE]

where $\hat{U}_{0}:=\mathbb{P}U_{0}$. Then, for $\varepsilon\ll 1$, there holds the estimate:

[TABLE]

Remark 2.2*.*

It is easy to see that the projected subsystem (2.15)-(2.16) is equivalent to 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 $\mathcal{E}$, solutions of the compressible Euler system (2.1) approximate those of the incompressible Euler system (2.5) when $\varepsilon\to 0$. The theorem, clearly, is in agreement with the results obtained in [25, 26].

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

In the rest of this section, we consider the linear wave equation system as a prototype model of the compressible Euler equations (2.6), and present an analogous analysis of the zero Mach number limit. Linearising the system of equations (2.6) about a constant state (\underline{\rho},\underline{\mbox{\boldmathu}}), we get the following wave equation system with advection:

[TABLE]

where

[TABLE]

Before carrying out an analysis of the solutions of (2.18), we introduce a few definitions. For any U_{1}=(\rho_{1},\mbox{\boldmathu}_{1}) and U_{2}=(\rho_{2},\mbox{\boldmathu}_{2}) in $L^{2}(\mathbb{T}^{d})^{1+d}$, we define the inner-product

[TABLE]

where the inner-products appearing on the right hand side are the usual $L^{2}$ inner-products. The norm generated by the above inner-product is given by

[TABLE]

Clearly, the norm ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$ is equivalent to the usual norm in $L^{2}(\mathbb{T}^{d})^{1+d}$. We also consider the energies

[TABLE]

which are, respectively, the total energy, the incompressible energy, and the acoustic or compressible energy. Since $\mathcal{E}\perp\tilde{\mathcal{E}}$, we have $E=E_{in}+E_{ac}$. In the following, we state the corresponding linearised version of Schochet’s result, i.e. Theorem 2.1.

Theorem 2.3** ([14]).**

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

[TABLE]

and let $\bar{U}$ be a solution of the IVP:

[TABLE]

where $\hat{U}_{0}:=\mathbb{P}U_{0}$. Let $U=\hat{U}+\tilde{U}$ be the Helmholtz-Leray decomposition of $U$. Then the following holds.

1. (i)

$\hat{U}=\bar{U}$, 2. (ii)

$\tilde{U}$* is the solution of (2.23) with initial condition $\tilde{U}_{0}:=(\mathbb{I}-\mathbb{P})U_{0}$.*

Moreover, there holds the energy conservation:

[TABLE]

and as a consequence, the following estimate holds:

[TABLE]

Remark 2.4*.*

Carrying out an analogous asymptotic analysis using the ansatz (2.3) for the wave equation yields the system:

[TABLE]

for the unknowns (\rho_{(2)},\mbox{\boldmathu}_{(0)}), analogous to (2.5). Since the boundary conditions are periodic, and (2.29) is linear, the divergence-free condition (2.30) on \mbox{\boldmathu}_{(0)} now forces $\rho_{(2)}$ to be a constant. Hence, $\rho\equiv\mathrm{const.}$, and (2.29) then reduces to a linear transport equation for \mbox{\boldmathu}_{(0)}, and the divergence-free condition (2.30) is now required only at time $t=0$. Thus, the projected subsystem (2.25)-(2.26) and the system (2.29)-(2.30) are equivalent, and are incompressible.

Theorem 2.3 lies at the centre of the analysis of numerical schemes presented in [14]. 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 has to be generalised to accomodate more general spatial differential operators arising from the discretisation. Such a generalisation is presented in [14] which also gives a sufficient condition to ensure the estimate (2.28).

Theorem 2.5** ([14]).**

Let $U$ be a solution of the IVP:

[TABLE]

which is assumed to be well-posed in $L^{\infty}\left([0,\infty);L^{2}(\mathbb{T}^{d})^{1+d}\right)$, with \mathcal{F}_{\mbox{\boldmathx}} a linear spatial differential operator. Then the following conclusions hold.

1. (i)

The solution $U$ satisfies the estimate

[TABLE]

where $\bar{U}$ is a solution of (2.31) with the initial condition $\bar{U}(0)=\mathbb{P}U_{0}$. However, we don’t have the apriori estimate $\lVert U(t)-\mathbb{P}U(t)\rVert=\mathcal{O}(\varepsilon)$ for all $t>0$. 2. (ii)

When the operator \mathcal{F}_{\mbox{\boldmathx}} leaves $\mathcal{E}$ invariant, i.e. whenever $U_{0}\in\mathcal{E}$ implies $U(t)\in\mathcal{E}$ for all $t>0$, then $U$ satisfies the estimate (2.33), and in addition we have

[TABLE]

Remark 2.6*.*

Note that in Theorem 2.1 and 2.3 the estimate on the difference $\lVert\rho-\bar{\rho}\rVert$ is $\mathcal{O}(\varepsilon^{2})$, whereas in [14] a scaled variable is used instead of $\rho$, and hence the corresponding estimate is $\mathcal{O}(\varepsilon)$.

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 $\varepsilon\to 0$. However, it is wellknown from the literature, e.g. [14, 15, 19], that the accuracy of a scheme can also deteriorate in the low Mach number limit. Hence, it is essential for an AP scheme to be not only consistent with the stiff limit but also to maintain its accuracy. In addition, an AP scheme is expected to be accurate in the transient regimes as well. The two primary aims of the present work are to design and analyse, and test an IMEX-RK finite volume scheme which is AP while maintaining the accuracy uniformly in $\varepsilon$. In what follows, we formally state the definitions of the AP property, and the asymptotic accuracy (AA). The key principle behind the AA property is an estimate which states that it is essential for a numerical solution to be $\varepsilon$-close to the incompressible solution as in the continuous asymptotics. Following Jin [23], we show how the AA property of a scheme suffices to obtain the error estimates expected from an AP scheme.

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 parameters, and the transitional nature of the governing equations. A classical numerical approximation scheme may not resolve the scales, 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 [12, 22, 23], and the references therein for more details.

Definition 3.1**.**

Let $\mathcal{P}^{\varepsilon}$ denote a singularly perturbed problem with the perturbation parameter $\varepsilon$. Let $\mathcal{P}^{0}$ denote the limiting system of $\mathcal{P}^{\varepsilon}$ when $\varepsilon\to 0$. A discretisation $\mathcal{P}_{h}^{\varepsilon}$ of $\mathcal{P^{\varepsilon}}$, with $h$ being the discretisation parameter, is called AP if

1. (i)

$\mathcal{P}^{0}_{h}$ is a consistent discretisation of the problem $\mathcal{P}^{0}$, called the asymptotic consistency, and 2. (ii)

the stability constraints on $h$ are independent of $\varepsilon$, called the asymptotic stability.

In other words, the following diagram commutes:

[TABLE]

Jin in [23] has presented a systematic discussion regarding the error estimates expected from an AP scheme. The solutions of the continuous systems $\mathcal{P}^{\varepsilon}$ and $\mathcal{P}^{0}$ satisfy the following error estimate:

[TABLE]

Suppose $\mathcal{P}_{h}^{\varepsilon}$ is a $p^{th}$-order accurate approximation to $\mathcal{P}^{\varepsilon}$ for a fixed $\varepsilon$ with discretisation parameter $h$. Caused by the presence of the singular perturbation parameter $\varepsilon$, a classical numerical approximation is typically expected to give the following error estimate:

[TABLE]

If the scheme $\mathcal{P}_{h}^{\varepsilon}$ is assumed to be AP, it may give the following estimates

[TABLE]

Combining (3.1), (3.3) and (3.4), by using the triangle inequality, for an AP scheme we have the error estimate:

[TABLE]

Comparing the estimates (3.5) and (3.2), it can be concluded that unlike a classical scheme, the error $e_{ap}$ of an AP scheme remains bounded as $\varepsilon\to 0$. The estimate (3.3), which plays a crucial role to obtain the error bound (3.5), makes sure that the numerical solution of the singular perturbation problem $\mathcal{P}^{\varepsilon}$ remains $\varepsilon$-close to the numerical solution of the limit problem $\mathcal{P}^{0}$, which is a property exhibited by the exact solutions of $\mathcal{P}^{\varepsilon}$ and $\mathcal{P}^{0}$.

In a different context, Dellacherie in [14] has investigated the origin of inaccuracies exhibited by explicit Godunov-type schemes, arising from the creation of spurious acoustic waves. The results of the above study reveal that the inaccuracies originate due to the inability of the schemes to maintain the estimate (3.3) even if the initial datum is well-prepared. It has been shown in [14], cf. also Theorem 2.5, that the $\mathcal{E}$-invariance is a sufficient condition for linear schemes to satisfy the estimate (3.3), provided the initial datum is well-prepared.

3.2. Asymptotic Accuracy

A numerical scheme for the compressible Euler system should possess the ability to maintain the solutions close to the incompressible solutions in $\mathcal{E}$ for all times, whenever the initial data are well-prepared. It has to be noted that Theorem 2.5, and its conclusions are valid only for a linear scheme applied to the wave equation system. A scheme maintaining the estimate (2.34) from Theorem 2.5 was designated to be a low Mach accurate scheme in [14], and the $\mathcal{E}$-invariance is a sufficient condition at least in the linear case. It should be taken into account that the estimates (2.34) and (3.3) are the same in the case of the wave equation. Therefore, motivated by this result, the discussion on AP schemes, and analogous considerations from [14], we propose the asymptotic accuracy of a scheme to be its ability to preserve the well-prepared space $\mathcal{E}$.

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 $\mathcal{E}$ invariant.

Remark 3.3*.*

Even though the AA property is a sufficient condition to obtain the estimate (3.3), it is very useful as it gives an easy and verifiable criterion for the wave equation system, leading towards the AP property.

The estimate (3.4) for an AP scheme gives $p^{th}$-order accuracy in the stiff limit $\varepsilon\to 0$. Thus, if the solutions of the modified partial differential equation (MPDE) resulting from a $p^{th}$-order linear scheme applied to the wave equation leaves $\mathcal{E}$ invariant, then such a scheme will avoid inaccuracies. In addition, the scheme will reduce to a $p^{th}$-order accurate discretisation of the limit system in the stiff limit.

Remark 3.4*.*

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

[TABLE]

where $y\colon\mathbb{R}\to\mathbb{R}^{n},\ f,g\colon\mathbb{R}\times\mathbb{R}^{n}\to\mathbb{R}^{n}$, and $0<\varepsilon\ll 1$ is usually known as the stiffness parameter. The functions $f$ and $g$ 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. [21], for a comprehensive treatment of such systems.

Let $y^{n}$ be a numerical solution of (4.1) at time $t^{n}$, and let $\Delta t$ denote a fixed timestep. An $s$-stage IMEX-RK scheme, cf., e.g. [2, 32], updates $y^{n}$ to $y^{n+1}$ through $s$ intermediate stages:

[TABLE]

Let us denote $\tilde{A}=(\tilde{a}_{i,j})$, $A=(a_{i,j})$, $\tilde{c}=(\tilde{c}_{1},\tilde{c}_{2},\ldots,\tilde{c}_{s}),\ c=(c_{1},c_{2},\ldots,c_{s}),\ \tilde{\omega}=(\tilde{\omega}_{1},\tilde{\omega}_{2},\ldots,\tilde{\omega}_{s})$ and $\omega=(\omega_{1},\omega_{2},\ldots,\omega_{s})$. The above IMEX-RK scheme (4.2)-(4.3) can be symbolically represented by the double Butcher tableau:

The matrices $A$ and $\tilde{A}$ are $s\times s$ matrices such that resulting scheme is explicit in $f$ and implicit in $g$. 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 $\tilde{a}_{i,j}=0$ for $j\geq i$, and $a_{i,j}=0$ for $j>i$. The coefficients $\tilde{c}_{i}$ and $c_{i}$, and the weights $\tilde{\omega}_{i}$ and $\omega_{i}$ are fixed by the order conditions; see, e.g. [21, 32, 33], 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:

[TABLE]

The conditions (4.4), (4.5) and (4.6) are, respectively, the consistency, the first-order, and the second-order order conditions.

Definition 4.1**.**

An IMEX-RK scheme with the Butcher tableau given in Figure 1 is said to be globally stiffly accurate (GSA), if

[TABLE]

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 [8, 24] for more details.

Definition 4.2**.**

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

• •

type-A, if the matrix $A$ is invertible;

• •

type-CK, if the matrix $A\in\mathbb{R}^{s\times s},\ s\geq 2$, can be written as

[TABLE]

where $\alpha\in\mathbb{R}^{s-1}$ and $A_{s-1}\in\mathbb{R}^{s-1\times s-1}$ is invertible.

In the results presented in the later sections, we have used predominantly the first order Euler(1,1,1), and second order ARS(2,2,2) schemes for time discretisations; see the Appendix for the Butcher tableaux of the several variants of the IMEX-RK schemes we have considered in our numerical case studies. Here, in the triplet $(s,\sigma,p)$, $s$ is the number of stages of the implicit part, the number $\sigma$ gives the number of stages for the explicit part and $p$ gives the overall order of the scheme. For a detailed account of IMEX-RK schemes we refer the reader to [2, 8, 24, 32, 33] 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. Denoting by W=(\rho,\mbox{\boldmathq})^{T}, the vector of conservative variables, the fluxes are split via

[TABLE]

where \mbox{\boldmathq}=\rho\mbox{\boldmathu} is the momentum. Applying the IMEX-RK time discretisation, i.e. treating $F$ explicitly and $G$ implicitly, results in the following semi-discrete scheme.

Definition 4.3**.**

The $k^{th}$ stage of an $s$-stage IMEX-RK scheme for the Euler system (2.1) is defined as

[TABLE]

The numerical solution at time $t^{n+1}$ is given by

[TABLE]

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 $\ell$ and $l$ are used to denote, respectively, the summation in the explicit and implicit terms, i.e. they assume values in the sets $\{1,2,\ldots,k-1\}$ and $\{1,2,\ldots,k\}$.

Though the intermediate stages (4.9)-(4.10) consist of two implicit steps, the resolution of the scheme (4.9)-(4.10) is quite simple. As in [13], we eliminate $q^{k}$ between (4.9) and (4.10), and obtain the nonlinear elliptic equation:

[TABLE]

for $\rho^{k}$. Here, we have denoted by

[TABLE]

those terms that can be explicitly evaluated. Once $\rho^{k}$ is known after solving the elliptic equation (4.13), \mbox{\boldmathq}^{k} can be explicitly evaluated from (4.10). Finally, the updates $\rho^{n+1}$ and \mbox{\boldmathq}^{n+1} are calculated explicitly using (4.11)-(4.12) with the values obtained from the intermediate stages; see also [9, 13, 16, 17, 37] for related approaches.

Remark 4.4*.*

In order to further reduce the nonlinear nature of the elliptic equation (4.13), and to make the numerical implementation simpler, in our computations we transform the nonlinear elliptic equation for $\rho^{k}$ into a semi-linear elliptic equation for the pressure $p^{k}:=p(\rho^{k})$.

5. Analysis of the 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.9)-(4.12) 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 $\varepsilon\to 0$, with its stability requirements independent of $\varepsilon$. In order to show the former, we perform an asymptotic analysis of the time semi-discrete scheme (4.9)-(4.12), 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 $t^{n}$ are well-prepared, i.e. $\rho^{n}$ and \mbox{\boldmathu}^{n} admit the decomposition:

[TABLE]

where $\nabla\rho^{n}_{(0)}=0$ and \nabla\cdot\mbox{\boldmathu}^{n}_{(0)}=0. Then for a GSA scheme, if the intermediate solutions $\rho^{k}$ and \mbox{\boldmathu}^{k} defined by (4.9) and (4.10) admit the decomposition (2.3), then they must satisfy $\nabla\rho^{k}_{(0)}=0$ and \nabla\cdot\mbox{\boldmathu}^{k}_{(0)}=0, i.e. $\rho^{k}$ and \mbox{\boldmathu}^{k} are well-prepared as well. As a consequence, if the numerical solutions $\rho^{n+1}$ and \mbox{\boldmathu}^{n+1} admit the decomposition (2.3), then they are also well-prepared. In other words, the semi-discrete scheme (4.9)-(4.12) 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 begin with, we prove the consistency for the first stage, i.e. $k=1$, which corresponds to a fully implicit step. We plugin the ansatz (2.3) for each of the dependent variables in (4.9)-(4.12). Equating to zero the $\mathcal{O}{(\varepsilon^{-2})}$ terms in system (4.9)-(4.10) yields

[TABLE]

Hence, the zeroth order density $\rho^{1}_{(0)}$ is spatially constant. In an analogous way, we can show that the first order density $\rho^{1}_{(1)}$ is also spatially constant.

Next, we consider the $\mathcal{O}(1)$ terms in the mass update (4.9) to obtain

[TABLE]

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

[TABLE]

If the boundary conditions of the problem are periodic or wall, then the right most integral in the above equation vanishes, yielding $\rho^{1}_{(0)}=\rho^{n}_{(0)}$. As the zeroth order density $\rho^{n}_{(0)}$ is a constant, $\rho^{1}_{(0)}$ as well as $\rho^{1}_{(1)}$ are constants. Using this in equation (5.4), we obtain the divergence condition:

[TABLE]

for the leading order velocity \mbox{\boldmathu}^{1}_{(0)}. This completes the proof for $k=1$.

To prove the result for $k=2$ onwards, we rewrite the intermediate stages (4.9)-(4.10) in the form

[TABLE]

Proceeding as in the case of $k=1$, we can obtain $\rho^{k}_{(0)}=\rho^{n}_{(0)}$. Hence, for the $k^{th}$ stage, the zeroth order density $\rho^{k}_{(0)}$ is same as that of $\rho^{n}_{(0)}$, and in turn we also obtain the divergence condition:

[TABLE]

for \mbox{\boldmathu}^{k}_{(0)}.

Since the scheme under consideration is GSA, the numerical solution at time $t^{n+1}$ is same as the solution at the $s^{th}$ stage. Summarising the above steps, the asymptotic limit of scheme (4.9)-(4.12) is given by

[TABLE]

Clearly, (5.10) is a consistent discretisation of the incompressible limit system (2.5). Hence, the semi-discrete scheme (4.9)-(4.12) is asymptotically consistent. ∎

Remark 5.2*.*

Note that the above proof is valid for any $s$-stage GSA IMEX-RK scheme. However, it uses the fact that $a_{k,k}\neq 0$ for $k=1,2,\ldots,s$, which is the property of an RK scheme of type-A. The result also holds for a type-CK GSA 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 $L^{2}$-stability Analysis

In this subsection, we analyse the correction terms arising from the time discretisation, and their effect on the asymptotic stability of the resulting scheme. In order to analyse these correction terms, we follow the standard MPDE approach; see also [19] for a related analysis on stability.

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 $L^{2}$-stability analysis for the wave equation system (2.18) as a linearised model. As a result, we obtain a sufficient condition for $L^{2}$-stability; see [1] for the analysis of a first order IMEX scheme. We believe that the results of linear $L^{2}$-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 $k^{th}$ stage of an $s$-stage IMEX-RK scheme for the wave equation system (2.18) is given by

[TABLE]

The numerical solutions $\rho^{n+1}$ and \mbox{\boldmathu}^{n+1} at time $t^{n+1}$ are defined as

[TABLE]

In the following, we derive the MPDE corresponding to a general second order accurate time discrete scheme of the form (5.11)-(5.14). We show that the solution of the MPDE is energy-dissipative under a sufficient condition involving only the RK coefficients, and hence, the timesteps are independent of $\varepsilon$; see also [1]. Thus, the time-discrete scheme (4.9)-(4.12) for the Euler system is linearly asymptotically $L^{2}$-stable.

Theorem 5.4**.**

Consider a second-order time-discrete IMEX-RK scheme of the form (4.9)-(4.12), and let the coefficients $b^{(4)}_{1},b^{(4)}_{2},b^{(4)}_{3}$ and $b^{(4)}_{4}$ be defined by (8.1) in the Appendix. Then, the scheme is linearly $L^{2}$-stable as long as the timesteps are bounded, and the coefficients $b^{(4)}_{1},b^{(4)}_{2},b^{(4)}_{3}$ and $b^{(4)}_{4}$ are negative.

Proof.

We consider a general second order accurate IMEX-RK time discretisation in (5.11)-(5.14). 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:

[TABLE]

Here, the operators $B^{(3)}$ and $B^{(4)}$ contain the third and fourth derivatives of the unknown functions $\rho$ and $u$, 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.15) is free from first and second order derivatives.

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

[TABLE]

We use the MPDE (5.15) in (5.16), 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 the periodic boundary conditions, and hence we get

[TABLE]

[TABLE]

Applying the Cauchy-Schwarz inequality in (5.17) and (5.18), and using the inequalities thus obtained in (5.16), finally yields

[TABLE]

It is easy to see that the right hand side of the above inequality is negative under the hypotheses of the theorem. ∎

Remark 5.5*.*

The coefficients $b^{(4)}$ given in the Appendix are for a general second order IMEX-RK scheme. It can be seen that all the second order variants we have considered, namely ARS(2,2,2), PR(2,2,2), Jin(2,2,2) and RK2CN(2,2,2) are linearly $L^{2}$-stable by verifying the hypotheses of the above theorem.

5.2. $\mathcal{E}$-invariance and Asymptotic Accuracy

In this subsection, we prove the $\mathcal{E}$-invariance of the IMEX-RK time discrete scheme (4.9)-(4.12) and its asymptotic accuracy. It has been shown in [14] 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 [14] shows that there is an acoustic time scale $\tau_{ac}$ so that a Godunov-type scheme suffers from inaccuracies beyond a time interval of $\mathcal{O}{(\tau_{ac})}$ due to the creation of spurious acoustic waves. A cure proposed in [14] is to remove the numerical diffusion terms from the MPDE corresponding to 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 [3] for a related discussion.

In the following, we prove the $\mathcal{E}$-invariance of the time semi-discrete scheme. However, analysing the numerical solutions or the solutions of an MPDE resulting from a general second order IMEX time discretisation is difficult, and hence we perform the corresponding analysis on the linear wave equation system. The results of linear analysis shows that the MPDE leaves the well-prepared space $\mathcal{E}$ invariant without requiring any change in the diffusion matrices. As a result, it follows from Theorem 2.5 that the solution of the MPDE are $\varepsilon$-close to the incompressible solution, satisfying the estimates (2.34) or (3.3) which is crucial to maintain the asymptotic accuracy.

Theorem 5.6**.**

Suppose that at time $t^{n}$ the numerical solution (\rho^{n},\mbox{\boldmathu}^{n}) to the compressible Euler system (2.1) is in $\mathcal{E}$, i.e. $\nabla\rho^{n}(x)=0$ and \nabla\cdot\mbox{\boldmathu}^{n}(x)=0, for all $x$. Then at time $t^{n+1}$, the numerical approximation (\rho^{n+1},\mbox{\boldmathu}^{n+1}) obtained from the scheme (4.9)-(4.12) satisfy

[TABLE]

In other words, the semi-discrete scheme (4.9)-(4.12) keeps the well-prepared space $\mathcal{E}$ invariant.

Proof.

We use induction on the number of stages to prove the theorem. As a first step, we prove the well-preparedness for the first-stage, i.e. $k=1$, which corresponds to a fully implicit step. Multiplying (4.10) by $\varepsilon^{2}$ and letting $\varepsilon\to 0$, we get

[TABLE]

Hence, the density $\rho^{1}$ at the first-stage is spatially constant. Next, we consider the mass update (4.9) and obtain

[TABLE]

We Integrate the above equation (5.22) over a domain $\Omega$, and use the Gauss’ divergence theorem to get

[TABLE]

Using periodic or wall boundary conditions, we can see that the right most integral in the above equation (5.5) vanishes. As a result, $\rho^{1}=\rho^{n}$. Since the density $\rho^{n}$ at time $t^{n}$ is a constant, we immediately obtain that $\rho^{1}$ is also a constant. Using this in equation (5.22) yields the divergence condition

[TABLE]

for the velocity \mbox{\boldmathu}^{1} at the first-stage. This completes the proof for $k=1$.

For $k=2$ and then on, we consider the intermediate stages (4.9)-(4.10) in the form (5.7)-(5.8). Now proceeding as in the case of $k=1$, we can obtain $\rho^{k}=\rho^{n}$. Hence, for the $k^{th}$-stage, the density$\rho^{k}$ is same as that of $\rho^{n}$, and in turn we also obtain the divergence condition:

[TABLE]

As the scheme is assumed to be GSA, the numerical solution (\rho^{n+1},\mbox{\boldmathu}^{n+1}) obtained from the scheme (4.9)-(4.12) lives in $\mathcal{E}$, whenever the initial data (\rho^{n},\mbox{\boldmathu}^{n}) lives in $\mathcal{E}$. In other words, the IMEX-RK scheme (4.9)-(4.12) leaves the well-prepared space $\mathcal{E}$ invariant. ∎

Remark 5.7*.*

It has to be noted that the analysis carried out in the above theorem doesn’t take into account the correction terms introduced by the time discretisation. In order to substantiate the claims of the formal analysis thus performed in Theorem 5.6, and to obtain the estimate (2.34) for the solution of the MPDE (5.15), we prove its $\mathcal{E}$-invariance in the following proposition.

Proposition 5.8**.**

The modified equation system (5.15) for the IMEX-RK scheme (4.9)-(4.12) applied to the wave equation system is $\mathcal{E}$-invariant, i.e. if the initial data (\rho(0,\mbox{\boldmathx}),\mbox{\boldmathu}(0,\mbox{\boldmathx})) is in $\mathcal{E}$, then the solution (\rho(t,\mbox{\boldmathx}),\mbox{\boldmathu}(t,{x})) of the MPDE (5.15) is in $\mathcal{E}$ for all times $t>0$.

Proof.

Since the problem (2.6) is given on the torus $\mathbb{T}^{d}$, we can assume that it is given in the whole of $\mathbb{R}^{d}$ by periodic extension. It is convenient to work in the Fourier variables as it simplifies the calculations. Taking the Fourier transform of (5.15) in space we obtain

[TABLE]

where the matrices $\hat{B}^{(3)}$ and $\hat{B}^{(4)}$ are, respectively, the Fourier transforms of the matrices $B^{(3)}$ and $B^{(4)}$ introduced in Theorem 5.4.

Note that (\rho,\mbox{\boldmathu})\in\mathcal{E} if, and only if, (\hat{\rho},\hat{\mbox{\boldmathu}})\in\hat{\mathcal{E}}, where $\hat{\mathcal{E}}$ is given by

[TABLE]

Hence, the quantities of interest are $\hat{\rho}$ and \mbox{\boldmath\xi}\cdot\hat{\mbox{\boldmathu}}. Since the first equation of (5.26) is already an equation for $\hat{\rho}$, taking the dot product of the second equation for \hat{\mbox{\boldmathu}} in (5.26) with $\xi$ and combining with the first equation of (5.26) for $\hat{\rho}$ yields a linear system:

[TABLE]

where M(\mbox{\boldmath\xi}) is a matrix consisting of polynomials in $\xi$ obtained from (5.26). Let us suppose that (\rho(0,\cdot),\mbox{\boldmathu}(0,\cdot))\in\mathcal{E}, i.e. \hat{\rho}(0,\mbox{\boldmath\xi})=\mbox{\boldmath\xi}\cdot\hat{\mbox{\boldmathu}}(0,\mbox{\boldmath\xi})=0. The solution of (5.28) then clearly shows that \hat{\rho}(t,\mbox{\boldmath\xi})=\mbox{\boldmath\xi}\cdot\hat{\mbox{\boldmathu}}(t,\mbox{\boldmath\xi})=0. Hence, solution of (5.15) lives in $\mathcal{E}$ for all times $t>0$, if it was in $\mathcal{E}$ at time $t=0$. ∎

Remark 5.9*.*

Combining the above results, we can conclude that the time semi-discrete IMEX-RK scheme (4.9)-(4.12) is AP and AA in the sense of Definitions 3.1 and 3.2.

6. Analysis of the Space-time Fully-discrete Scheme

This section is devoted to the analysis of the fully-discrete IMEX-RK scheme resulting after a space discretisation by the finite volume method. For the sake of simplicity, we consider only two-dimensional (2-D) problems, and we assume that the given Cartesian spatial domain $\Omega$ is discretised into rectangular mesh cells of lengths $\Delta x_{1}$ and $\Delta x_{2}$ in $x_{1}$- and $x_{2}$-directions, respectively. For notational conveniences, let us also introduce the following finite difference operators $\mu$ and $\delta$, defined via

[TABLE]

in the $x_{1}$-direction, with analogous definitions in the $x_{2}$-direction.

Given an approximation $W_{i,j}^{n}$ to the piecewise constant cell averages of the conserved variable $W$ at time $t^{n}$, we reconstruct all the conservative variables using a standard MUSCL-type piecewise linear interpolant

[TABLE]

in every computational cell. Here, $W_{i,j}^{\prime}$ and $W_{i,j}^{\backprime}$ are, respectively, the discrete slopes in $x_{1}$- and $x_{2}$-directions. It has to be noted that the analysis presented in Section 2 can be carried out to the fully-discrete scheme, only if we restrict ourselves to smooth solutions. In such problems, the discrete slopes $W_{i,j}^{\prime}$ and $W_{i,j}^{\backprime}$ are approximated using central differences, i.e.

[TABLE]

However, if the solution under consideration is discontinuous, we use slope limitting techniques to obtain $W_{i,j}^{\prime}$ and $W_{i,j}^{\backprime}$, e.g.

[TABLE]

in the $x_{1}$-direction. Here, the $\mathrm{CWENO}$ function is defined by

[TABLE]

Let $F_{1}$ and $F_{2}$ be, respectively, the non-stiff fluxes in $x_{1}$- and $x_{2}$-directions, and let $G_{1}$ and $G_{2}$ be, respectively, the stiff fluxes in $x_{1}$- and $x_{2}$-directions, cf. (4.8). The eigenvalues of the Jacobians of $F_{1}(W)$ and $F_{2}(W)$ with respect to $W$ can be obtained as $\lambda_{1,1}=0,\lambda_{1,2}=u_{1},\lambda_{1,3}=2u_{1}$, and $\lambda_{2,1}=0,\lambda_{2,2}=u_{2},\lambda_{2,3}=2u_{2}$. The timestep $\Delta t$ at time $t^{n}$ is computed by the CFL condition:

[TABLE]

where $\nu<1$ is the given CFL number. Note that the above condition is almost like the advective CFL condition, and is independent of $\varepsilon$.

6.1. Space-time Fully-discrete Scheme

Discretising the semi-discrete scheme (4.9)-(4.12) using the finite volume method, and replacing the physical fluxes by suitable numerical fluxes yields the fully-discrete scheme. In order to introduce the numerical fluxes, we first rewrite the IMEX time semi-discrete scheme (4.9)-(4.12) in terms of the stiff and non-stiff flux functions, via

[TABLE]

Definition 6.1**.**

The $k^{th}$ stage of an $s$-stage fully-discrete IMEX-RK scheme for the Euler system (2.1) is defined as

[TABLE]

with the numerical solution at time $t^{n+1}$ given by

[TABLE]

where the repeated index $m$ takes values in $\{1,2\}$, $\nu_{m}:=\frac{\Delta t}{\Delta x_{m}}$ denote the mesh ratios, and $\mathcal{F}_{m}$ and $\mathcal{G}_{m}$ are, respectively, the numerical fluxes used to approximate the physical fluxes $F_{m}$ and $G_{m}$. In this paper, we have used a simple Rusanov-type flux for $\mathcal{F}_{m}$ and central flux for $\mathcal{G}_{m}$, defined, e.g. along $x_{1}$- and $x_{2}$-directions as

[TABLE]

Here, $W^{k\pm}_{i+\frac{1}{2},j}$ and $W^{k\pm}_{i,j+\frac{1}{2}}$ are the interpolated states obtained using the piecewise linear reconstructions, and the wave-speeds are computed as, e.g. in the $x_{1}$-direction

[TABLE]

The rest of this section is devoted to establishing the AP and AA properties of the fully-discrete scheme (6.8)-(6.9) as done in Section 5 for the case of the time semi-discrete scheme.

6.2. AP Property

6.2.1. Asymptotic Consistency

Theorem 6.2**.**

Suppose that the data at time $t^{n}$ are well-prepared, i.e. $\rho^{n}$ and \mbox{\boldmathu}^{n} admit the decomposition:

[TABLE]

where $\hat{\nabla}\rho^{n}_{(0),i,j}=0$ and \hat{\nabla}\cdot\mbox{\boldmathu}^{n}_{(0),i,j}=0. Here, $\hat{\nabla}$ is the discrete gradient introduced by the implicit terms, i.e. by replacing the derivatives by central differences. Then for a GSA scheme, if the intermediate solutions $\rho^{k}_{i,j}$ and \mbox{\boldmathu}^{k}_{i,j} defined by (6.8) admit the same decomposition (6.12)-(6.13), they satisfy $\hat{\nabla}\rho^{k}_{(0),i,j}=0$ and \hat{\nabla}\cdot\mbox{\boldmathu}^{k}_{(0),i,j}=0, i.e. they are well-prepared as well. As a consequence, if the numerical solutions $\rho^{n+1}_{i,j}$ and \mbox{\boldmathu}^{n+1}_{i,j} admit the decomposition (6.12)-(6.13), then they are also well-prepared. In other words, the fully discrete scheme (6.8)-(6.9) is asymptotically consistent in the sense of Definition 3.1.

Proof.

We follow similar lines as in the proof of Theorem 5.1 to establish the asymptotic consistency of the fully-discrete scheme. We plugin the ansatz (2.3) in a discrete form for all the dependent variables in (6.8)-(6.9). Equating to zero the $O(\varepsilon^{-2})$ terms in the first stage, i.e. $k=1$, which is a fully implicit step, gives

[TABLE]

Hence, the zeroth order density $\rho^{1}_{(0),i,j}$ is constant for all $i,j$, if we assume periodic or wall boundary conditions. In an analogous way, we can show that the first order density $\rho^{1}_{(1),i,j}$ is constant for all $i,j$.

Next, we consider the $\mathcal{O}(1)$ terms in the mass update of (6.8) to obtain

[TABLE]

If the boundary conditions of the problem are periodic or wall, it can easily be seen that the right hand side of equation (6.15), when summed over the entire domain, vanishes to yield

[TABLE]

As a result, $\rho^{1}_{(0),i,j}=\rho^{n}_{(0),i,j}$, for all $i,j$. Since the zeroth order density $\rho^{n}_{(0)}$ is a constant, $\rho^{1}_{(0),i,j}$ as well as $\rho^{1}_{(1),i,j}$ are constants. Using this in equation (6.15), we obtain the discrete divergence condition:

[TABLE]

for the leading order velocity \mbox{\boldmathu}^{1}_{(0)}, i.e. the discrete divergence vanishes for all $i,j$. This completes the proof for $k=1$.

For $k=2$ and then on, we can mimic the steps carried out in the proof of Theorem 5.1 for each of the intermediate stages, and conclude that the discrete gradient of the density $\rho^{k}_{(0)}$, and the discrete divergence of the velocity \mbox{\boldmathu}^{k}_{(0)} vanish. Finally, we can prove that the numerical solution at time $t^{n+1}$ is also well-prepared along similar lines as in the proof of Theorem 5.1.

Summarising the above steps, we can clearly see that the limt of the scheme (6.9) is a consistent discretisation of the incompressible limit system (2.5). Hence, the fully-discrete scheme (6.8)-(6.9) is asymptotically consistent. ∎

6.2.2. Linearised $L^{2}$-stability Analysis

We now proceed to establish the linear $L^{2}$-stability of the fully-discrete IMEX-RK scheme (6.8)-(6.9) by applying it to the linear wave equation system. Analogous to Subsection 5.1.2, we obtain a sufficient condition for stability which involves only the RK coefficients. The derivation of MPDE corresponding to a space-time fully-discrete general second order IMEX scheme turns out to be extremely complicated and hence, we restrict ourselves only to a first order IMEX discretisation; see [32] for an example of a first order IMEX scheme other than Euler(1,1,1).

Theorem 6.3**.**

A first order space-time fully-discrete IMEX-RK scheme given by (6.8)-(6.9) is linearly $L^{2}$-stable under a CFL-like condition $\Delta t=C\min(\Delta x_{1},\Delta x_{2})$, where $C$ is a constant independent of $\varepsilon$, if

[TABLE]

where the coefficients $b_{1}^{(2)}$ and $b_{3}^{(2)}$ are defined by (8.1) in Appendix.

Proof.

The MPDE for the IMEX-RK scheme (6.8)-(6.9) can be obtained as

[TABLE]

where the operator $\mathcal{B}^{(2)}$ is defined in Appendix. As done in the case of the time semi-discrete scheme in Theorem 5.4, we use the MPDE (6.19) in the equation (5.16) to calculate the rate of change of energy. The right hand side of (5.16) then gives

[TABLE]

[TABLE]

Regrouping the terms in the above equations after applying the Cauchy-Schwarz and AM-GM inequalities, and then using the inequalities thus obtained in (5.16), finally yields

[TABLE]

If the coefficients $-b^{(2)}_{1}$ and $b^{(2)}_{3}$ are positive, it can easily be seen that the overall expression on the right hand side of the above inequality is negative under the following CFL-like condition:

[TABLE]

∎

Remark 6.4*.*

Note that $-b^{(2)}_{1}$ and $b^{(2)}_{3}$ are positive for Euler(1,1,1). An analogous $L^{2}$-stability analysis is carried out in [17] and a condition similar to (6.23) is obtained. Further, the analysis in [17] also show that more stringent conditions are to be enforced to obtain $L^{\infty}$-stability.

6.3. $\mathcal{E}$-invariance and Asymptotic Accuracy

As done in Section 5 for the time semi-discrete scheme, we first formally prove that the fully-discrete scheme (6.8)-(6.9) leaves $\mathcal{E}$ invariant. To obtain the estimate (2.34), we also carry out an analogous MPDE analysis of the linearised scheme. Since the proofs of the following results follow similar lines as that of their semi-discrete counterparts, we do not intent to present the details.

Theorem 6.5**.**

Suppose that at time $t^{n}$ the numerical solution (\rho^{n}_{i,j},\mbox{\boldmathu}^{n}_{i,j}) to the compressible Euler system (2.1) is in $\mathcal{E}$, i.e. $\rho^{n}_{i,j}=\mathrm{const.}$ and \hat{\nabla}\cdot\mbox{\boldmathu}^{n}_{i,j}=0 for all $i,j$. Then, at time $t^{n+1}$, the numerical approximation (\rho^{n+1}_{i,j},\mbox{\boldmathu}^{n+1}_{i,j}) obtained from the scheme (6.8)-(6.9) satisfy

[TABLE]

In other words, the fully-discrete scheme (6.8)-(6.9) keeps the well-prepared space $\mathcal{E}$ invariant.

Proposition 6.6**.**

The modified equation system (6.19) for a first order IMEX-RK scheme (6.8)-(6.9), applied to the linear wave equation system (2.18), is $\mathcal{E}$-invariant, i.e. if the initial data (\rho(0,\mbox{\boldmathx}),\mbox{\boldmathu}(0,\mbox{\boldmathx})) is in $\mathcal{E}$, then the solution (\rho(t,\mbox{\boldmathx}),\mbox{\boldmathu}(t,\mbox{\boldmathx})) of the MPDE (6.19) is in $\mathcal{E}$ for all times $t>0$.

7. Numerical Case Studies

In this section, we present the results of our numerical experiments in order to substantiate the claims made in the previous sections. In the first test problem, we consider an advecting vortex for which a smooth exact solution is available. We make use of the exact solution to test the experimental order of convergence (EOC) for different $\varepsilon$. The results obtained clearly demonstrate the uniform second order convergence with respect to $\varepsilon$. In addition, the results also show that the dissipation of the scheme is independent of $\varepsilon$. In the second test problem, we start with an exact analytical solution of the incompressible equations, and use this solution to measure the convergence of the numerical solution to the incompressible solution. The convergence study clearly yields second order convergence for very small values of $\varepsilon$, or the stiff accuracy in the incompressible regime. In the last problem, we consider a discontinuous solution corresponding to a cylindrical explosion therein we set $\varepsilon=1$. The scheme captures the shock wave and gives a good performance in the compressible regime as well.

7.1. Experimental Order of Convergence

Drawing inspiration from the traveling vortex problem studied in [7], we appropriate the initial conditions for the isentropic Euler system as follows.

[TABLE]

where r=\lVert\mbox{\boldmathx}-(0.5,0.5)\rVert, $\Gamma=1.5$, $\omega=4\pi$, and $k(r)=2\cos r+2r\sin r+\frac{1}{8}\cos(2r)+\frac{1}{4}r\sin(2r)+\frac{3}{4}r^{2}$. Here, $\Gamma$ is a parameter known as the vortex intensity, $r$ denotes the distance from the core of the vortex, and $\omega$ is an angular wave frequency specifying the width of the vortex. The Mach number $\varepsilon$ is controlled by adjusting the value of $\eta$ via the relation $\varepsilon={0.6\eta}/{\sqrt{110}}$.

The above problem admits an exact solution $(\rho,u_{1},u_{2})(t,x_{1},x_{2})=(\rho,u_{1},u_{2})(0,x_{1}-0.6t,x_{2})$. The computational domain $\Omega=[0,1]\times[0,1]$ is successively divided to $10\times 10,20\times 20$ up to $80\times 80$ square mesh cells. All the four boundaries are set to be periodic. The computations are performed up to a time $T=0.1$ using the second order ARS(2,2,2) scheme. The CFL number is set as $0.45$. The EOC computed in $L^{1}$ and $L^{2}$ norms using the above exact solution, for $\varepsilon$ ranging in $\{10^{-6},10^{-5},\ldots,10^{-1}\}$, are given in Table 1. The table clearly shows uniform second order convergence of the scheme with respect to $\varepsilon$.

Next, we give in Figure 2, the pcolor plot of the initial Mach number distribution at time $T=0$. In Figure 3, the pcolor plots of physical Mach numbers at time $T=1.67$, which is the taken by the vortex to complete one cycle, for the above range of values of $\varepsilon$ is given. Comparing each of the Mach profiles in Figure 3 with the initial one in Figure 2, it is very evident that the vortex does not deform or degrade much as it advects. Visually, one cannot observe any differences among the plots in Figure 3, which confirms that the dissipation of the vortex is independent of $\varepsilon$. This fact is confirmed also by the behaviour of the relative kinetic energy, and the vorticity for different $\varepsilon$. In Figure 4, the distribution of the relative kinetic energy as a function of time from $T=0$ to $T=1.67$, and the cross section of the vorticity at $x_{2}=0$ and $T=1.67$ as a function of $x_{1}$ is presented. The plots show that the relative kinetic energy is almost equal to one, independent of $\varepsilon$. It is quite remarkable that even for $\varepsilon=10^{-6}$, the energy dissipation is only $0.03\%$. Similarly, the vorticity plot also confirms that the dissipation is almost independent of $\varepsilon$.

7.2. Asymptotic Order of Convergence

The aim of this case study is to demonstrate the asymptotic convergence of the scheme, i.e. its ability to converge to the incompressible solution as $\varepsilon\to 0$, and to show that the convergence rate is two.

We consider the following exact solution:

[TABLE]

of the incompressible system (2.5) as given in [34], with $\rho_{(0)}(t,x_{1},x_{2})=1$. The computational domain $\Omega=[0,1]\times[0,1]$ is divided into $10\times 10,20\times 20,40\times 40,80\times 80$ mesh cells and the CFL number is set to $0.45$. The data are initialised using (7.1) at $t=0$ and the computations are done up to a final time $T=3$. The boundaries are periodic everywhere. The EOC computed using the exact solution (7.1) as the reference solution is termed as the asymptotic order of convergence (AOC). Our simulation results show that the density remains exactly as $\rho\equiv 1$, and we use the velocity components to measure the AOC. In Table 2, we present the AOC computed using three very small values of $\varepsilon$, namely $10^{-6},10^{-5}$ and $10^{-4}$. It is very clear from the data that the scheme indeed converges to the incompressible solution with second order accuracy as $\varepsilon\to 0$. It is also worth remarking that the present results are in conformity with the fact that the variant ARS(2,2,2) used here is GSA [2]. We have also tested non-GSA variants, such as PR(2,2,2), and the results show only marginal differences.

7.3. Cylindrical Explosion Problem

We consider a 2-D cylindrical explosion problem motivated by [9, 17] for the Euler equations (2.1) with a linear relation between the pressure and density $P(\rho)=\rho$.

The computations are carried out on the square $\Omega=[-1,1]^{2}$. The initial density profile reads

[TABLE]

In (7.2), $r=\sqrt{x_{1}^{2}+x_{2}^{2}}$ is the distance of any point $(x_{1},x_{2})$ from the origin $(0,0)$. At time $T=0$, the velocity of the fluid is given by

[TABLE]

where the coefficient $\alpha(x_{1},x_{2})$ is given by $\alpha:=\text{max}(0,1-r)(1-e^{-16r^{2}})$ and $(u_{1},u_{2})$ is set to $(0,0)$, if $r<10^{-15}$. All the boundaries are assumed to be periodic. The domain is divided uniformly by a $100\times 100$ mesh. We have used the JIN(2,2,2) variant to do the time discretisation.

In order to simulate a compressible regime, we first set $\varepsilon=1$. The surface plots of the density, and quiver plots of the velocity field at times $T=0.1,0.24$ and $0.5$ are given in Figure 5. Clearly, the circular shockwave moving outwards can be observed in the both the figures, confirming the good performance of the scheme in the fully compressible case.

In the same problem, we set $\varepsilon=10^{-3}$ to go to the weakly compressible or nearly incompressible regime. The surface plots of the density and divergence of the velocity is given Figure 6. It can be seen that the density is equal to one, and the velocity divergence is very small, indicating the convergence of the solution to the incompressible one.

8. Concluding Remarks

We have presented a class of second order accurate IMEX-RK finite volume schemes for the compressible Euler equations in the zero Mach number limit. Guided by the results of asymptotic analysis, the nonlinear fluxes in the Euler equations are split into stiff and non-stiff parts in which stiff terms correspond to fast acoustic waves and the non-stiff terms to slow advection waves. In the time discretisation, the stiff terms are treated implicitly, and the non-stiff terms explicitly. The resulting time semi-discrete scheme is shown to be AP and AA. In order to derive a space-time fully discrete scheme, a Rusanov-type central flux used for the non-stiff terms and central differences for stiff terms. The fully-discrete scheme is also shown to be AP and AA. A sufficient for $L^{2}$-stability, involving only the RK coefficients, overcoming the stiffness due to CFL restrictions is proposed for the time semi-discrete as well as space-time fully-discrete schemes. The numerical experiments confirm that the scheme achieves uniform second order convergence with respect to the Mach number, its dissipation is independent of the Mach number, and it converges to the incompressible solution with second order convergence as the Mach number goes to zero.

Acknowledgements

The authors thank the anonymous referees whose comments have lead to a significant improvement of the manuscript.

Appendix

First, let us introduce the following shorthands.

[TABLE]

With the above notations, the entries in the matrices $\mathcal{B}^{(2)},B^{(3)}$ and $B^{(4)}$ read

[TABLE]

[TABLE]

[TABLE]