An Asymptotic Preserving Time Integrator for Low Mach Number Limits of the Euler Equations with Gravity

TL;DR

This paper introduces an asymptotic preserving time integrator for Euler equations under gravity, effectively capturing low Mach number limits like incompressible and Boussinesq regimes with stability and validated numerical results.

Contribution

It develops a semi-implicit AP scheme that handles stiff and non-stiff terms, ensuring consistency with low Mach limit systems and stability in gravitational flows.

Findings

Scheme is consistent with incompressible and Boussinesq limits.

Linear stability analysis confirms $L^2$-stability.

Numerical experiments validate theoretical predictions.

Abstract

We consider two distinguished asymptotic limits of the Euler equations in a gravitational field, namely the incompressible and Boussinesq limits. Both these limits can be obtained as singular limits of the Euler equations under appropriate scaling of the Mach and Froude numbers. We propose and analyse an asymptotic preserving (AP) time discretisation for the numerical approximation of the Euler system in these asymptotic regimes. A key step in the construction of the AP scheme is a semi-implicit discretisation of the fluxes and the source term. The non-stiff convective terms are treated explicitly whereas the stiff pressure-gradient and source term are implicit. The implicit terms are combined to get a nonlinear elliptic equation. We show that the overall scheme is consistent with the respective limit system when the Mach number goes to zero. A linearised stability analysis confirms the $L^2$-stability of the proposed scheme. The results of numerical experiments validate the theoretical findings.

Equations

$${\partial }_t\rho +\nabla \cdot (\rho u)$$

$${\partial }_t(\rho u)+\nabla \cdot (\rho u\otimes u)+\frac{\nabla p}{{\mathrm{Ma}}^2}$$

$$f(t,x)=f_{(0)}(t,x)+\epsilon f_{(1)}(t,x)+{\epsilon }^2{f}_{(2)}(t,x).$$

$${\partial }_t{u}_{(0)}+\nabla \cdot \left(u_{(0)}\otimes u_{(0)}\right)+\nabla p_{(2)}$$

$$\nabla \cdot u_{(0)}$$

$${\partial }_t{u}_{(0)}+\nabla \cdot \left(u_{(0)}\otimes u_{(0)}\right)+\nabla p_{(2)}$$

$$\nabla \cdot u_{(0)}$$

$${\rho }_{(1)}=1-\frac{x_3}{\gamma }.$$

$$\frac{{\rho }^{n+1}-{\rho }^n}{\Delta t}+\nabla \cdot q^{n+1}$$

$$\frac{q^{n+1}-q^n}{\Delta t}+\nabla \cdot {\left(\frac{q\otimes q}{\rho }\right)}^n+\frac{\nabla p({\rho }^{n+1})}{{\epsilon }^2}$$

$$-\frac{\Delta t^2}{{\epsilon }^2}\Delta P({\rho }^{n+1})-\frac{\Delta t^2}{{\epsilon }^{\alpha }}\nabla \cdot \left({\rho }^{n+1}{e}_3\right)+{\rho }^{n+1}={\rho }^n-\Phi ({\rho }^n,u^n),$$

$$\Phi ({\rho }^n,u^n):=\Delta t\nabla \cdot q^n+\Delta t^2{\nabla }^2:{\left(\frac{q\otimes q}{\rho }\right)}^n$$

$$-\nabla \cdot u_{(0)}^{n+1}=\frac{{\rho }_{(0)}^{n+1}-{\rho }_{(0)}^n}{{\rho }_{(0)}^{n+1}\Delta t}.$$

$$-\frac{1}{|\Omega |}{\int }_{\partial \Omega }{u}_{(0)}^{n+1}\cdot \nu d\sigma =\frac{{\rho }_{(0)}^{n+1}-{\rho }_{(0)}^n}{{\rho }_{(0)}^{n+1}\Delta t}.$$

$$\nabla \cdot u_{(0)}^{n+1}=0.$$

$$\frac{u_{(0)}^{n+1}-u_{(0)}^n}{\Delta t}+\nabla \cdot (u_{(0)}^n\otimes u_{(0)}^n)+p_{(2)}^{n+1}$$

$$\nabla \cdot u_{(0)}^{n+1}$$

$${\partial }_t\rho +(\bar{u}\cdot \nabla )\rho +\bar{\rho }\nabla \cdot u$$

$${\partial }_{t}u+(\bar{u}\cdot \nabla )u+\frac{{\bar{a}}^2}{\bar{\rho }{\epsilon }^2}\nabla \rho$$

$$\frac{{\rho }^{n+1}-{\rho }^n}{\Delta t}+(\bar{u}\cdot \nabla ){\rho }^n+\bar{\rho }\nabla \cdot u^{n+1}$$

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

$$G(\Delta t,\xi )=G(0,\xi )+\mathcal{O}(\Delta t)\text{ as }\Delta t\to 0.$$

$${\hat{U}}^{n+1}=G(\Delta t,\xi ){\hat{U}}^n,$$

$$G(\Delta t,\xi )=\gamma \left(\begin{array}{ccc}1& -i\Delta t\bar{\rho }{\xi }_1& -i\Delta t\bar{\rho }{\xi }_2\\ -i\Delta t\lambda {\xi }_1& 1+\Delta t^2\bar{\rho }\lambda {\xi }_2^2& -\Delta t^2\bar{\rho }\lambda {\xi }_1{\xi }_2\\ -i\Delta t\lambda {\xi }_2& -\Delta t^2\bar{\rho }\lambda {\xi }_1{\xi }_2& 1+\Delta t^2\bar{\rho }\lambda {\xi }_1^2\end{array}\right),\lambda =\frac{{\bar{a}}^2}{\bar{\rho }{\epsilon }^2}\text{and}\gamma =\frac{1-i\Delta t(\bar{u}\cdot \xi )}{1+\Delta t^2\bar{\rho }\frac{{\bar{a}}^2}{{\epsilon }^2}|\xi {|}^2}.$$

$$G(\Delta t,\xi )-G(0,\xi )=\Delta t\left(\begin{array}{ccc}-\frac{\Delta t\bar{\rho }\lambda |\xi {|}^2+i(\bar{u}\cdot \xi )}{1+\Delta t^2\bar{\rho }\lambda |\xi {|}^2}& -i\bar{\rho }{\xi }_1\gamma & -i\bar{\rho }{\xi }_2\gamma \\ -i\lambda {\xi }_1\gamma & -\frac{\Delta t\bar{\rho }\lambda {\xi }_1^2+i(\bar{u}\cdot \xi )(1+\Delta t^2\bar{\rho }\lambda {\xi }_2^2)}{1+\Delta t^2\bar{\rho }\lambda |\xi {|}^2}& -\Delta t^2\bar{\rho }\lambda {\xi }_1{\xi }_2\gamma \\ -i\lambda {\xi }_2\gamma & -\Delta t^2\bar{\rho }\lambda {\xi }_1{\xi }_2\gamma & -\frac{\Delta t\bar{\rho }\lambda {\xi }_2^2+i(\bar{u}\cdot \xi )(1+\Delta t^2\bar{\rho }\lambda {\xi }_1^2)}{1+\Delta t^2\bar{\rho }\lambda |\xi {|}^2}\end{array}\right).$$

$$\rho (0,x_1,x_2)$$

$$q_1(0,x_1,x_2)$$

$$q_2(0,x_1,x_2)$$

Full text

∗ Corresponding author: S. Samantaray

School of Mathematics

Indian Institute of Science Education and Research Thiruvananthapuram

Thiruvananthapuram - 695551, India

(Communicated by the associate editor name)

1. Introduction

The presence of sound/acoustic waves poses a major challenge in atmospheric and meteorological flow computations due to their fast characteristic time scales. Hence, in most of the practical computations, one relies on the so-called ‘sound-proof’ models in which the sound waves are eliminated. The incompressible equations, Boussinesq equations, pseudo-incompressible equations, anelastic equations etc. are sound-proof models frequently used in the literature, to name but a few. The derivation and analysis of sound-proof models, study of their regimes of validity etc. are topics of active research even today; see, e.g., [2] and the references cited therein for more details.

A powerful and systematic method to derive a sound-proof model is an asymptotic analysis of the Euler equations in which one or more of the non-dimensional quantities, such as the Mach, Froude or Rossby numbers, assume the role of limiting parameters [4]. However, from a mathematical point of view, a sound-proof model is often recognised as a singular limit of the Euler equations under appropriate scalings. In addition, sound-proof equation systems are typically of hyperbolic-elliptic in nature, as opposed to the purely hyperbolic compressible Euler equations. On the other hand, from a numerical point of view, approximation of singular limits poses several challenges: stiffness arising from stringent stability requirements, reduction of order of accuracy due to the presence of limiting parameters and so on.

The goal of the present work is to obtain the incompressible and Boussinesq equations as two distinguished singular limits of the Euler equations in a gravitational field under appropriate scalings of the Mach and Froude numbers. We present their numerical resolution via the so-called asymptotic preserving (AP) methodology. An AP discretisation for a singularly perturbed problem in general is a one which reduces to a consistent discretisation of the limit model when the limits of perturbation parameters are taken. In addition, the stability requirements of the discretisation should remain independent of the perturbation parameters; see [6]. A key step in the construction of our AP scheme is a semi-implicit time discretisation based on a splitting of the flux and source terms into stiff and non-stiff terms. We show the asymptotic consistency of the scheme with the incompressible and Boussinesq limits as the Mach number approaches zero. As a first step towards the stability of the scheme in the asymptotic regime, we perform an $L^{2}$-stability analysis of the proposed scheme on a linearised model, namely the wave equation system. The results of our numerical experiments presented here clearly validate the AP nature of the proposed scheme.

2. Isentropic Euler System with Gravity and Its Asymptotic

Limits

We consider the scaled, isentropic compressible Euler equations with gravity:

[TABLE]

where $\rho>0$ is the density and $u\in\mathbb{R}^{3}$ is the velocity vector. Here, $\nabla$, $\nabla\cdot$ and $\otimes$ are respectively the gradient, divergence and tensor product operators and $e_{3}$ is the unit vector in the $x_{3}$-direction. We assume a simplified equation of state of an isentropic process, therein the pressure is related to density via $p=P(\rho)=\rho^{\gamma}$, where $\gamma$ is a constant. In (2.1)-(2.2), the non-dimensional parameters $\mathrm{Ma}$ and $\mathrm{Fr}$ are respectively, the reference Mach and Froude numbers.

The goal of the present work is the numerical approximation of some distinguished asymptotic limits of the Euler system (2.1)-(2.2) which models slow convection in a highly stratified medium; see, e.g. [2, 4] for more details. In order to describe these asymptotic regimes, in the following, we consider two important scalings of $\mathrm{Ma}$ and $\mathrm{Fr}$ in terms of an infinitesimal parameter $\varepsilon$.

• •

$\mathrm{Ma}=\varepsilon$ and $\mathrm{Fr}=1$. In this case, the pressure gradient term dominates the gravity term and we obtain the low Mach number limit.

• •

$\mathrm{Ma}=\varepsilon$ and $\mathrm{Fr}=\sqrt{\varepsilon}$. In this case, the gravitational term is also significant we derive the Boussinesq limit.

As a first step towards the derivation of the low Mach and Boussinesq limits, we expand all the dependent variables using the following three-term ansatz:

[TABLE]

We do not intent to provide the details of the derivation, but refer the interested reader to [4] for more details.

2.1. Zero Mach Number Limit

We set $\mathrm{Ma}=\varepsilon$ and $\mathrm{Fr}=1$ in (2.1)-(2.2) and let $\varepsilon\to 0$ to obtain the zero Mach number limit model:

[TABLE]

The above system (2.4)-(2.5) is the standard incompressible Euler system for the unknowns $u_{(0)}$ and $p_{(2)}$.

Remark 2.1**.**

Throughout our analysis and the numerical experiments presented in this paper, we assume either periodic or wall boundary conditions. As a consequence, the leading order density $\rho_{(0)}$ is a constant and the leading order velocity $u_{(0)}$ is divergence-free. Therefore, both the zero Mach and Boussinesq limits fall in the category of ‘sound-proof’ models.

2.2. Boussinesq Limit

Now we set $\mathrm{Ma}=\varepsilon$ and $\mathrm{Fr}=\sqrt{\varepsilon}$ in (2.1)-(2.2). Letting $\varepsilon\to 0$ yields the Boussinesq model:

[TABLE]

Since the first order density $\rho_{(1)}$ appears in (2.6)-(2.7), we need a closure relation. Using the multiscale ansatz (2.3) in the equation of state $p=\rho^{\gamma}$ and using the hydrostatic balance $\nabla p_{(1)}=-\rho_{(0)}e_{3}$ gives

[TABLE]

Remark 2.2**.**

It has be noted that both zero Mach and the Boussinesq limit systems are hyperbolic-elliptic in nature.

3. Semi-implicit Time Discretisation

In this section we present the time discretisation of the Euler system (2.1)-(2.2) based on implicit-explicit (IMEX) Runge Kutta (RK) schemes. These schemes were originally designed for stiff ordinary differential equations; see .e.g. [5] and the references therein.

Let $0=t^{0}<t^{1}<\cdots<t^{n}<t^{n+1}<\cdots$ be an increasing sequence of times and let $\Delta t$ be the uniform time-step. Let us denote by $f^{n}(x)$, the approximation to the value of any function $f$ at time $t^{n}$, i.e. $f^{n}(x)\sim f(t^{n},x)$.

A first order accurate semi-discrete scheme for the Euler equations (2.1)-(2.2) is defined as

[TABLE]

Here, $q=\rho u$ denotes the momentum and $\alpha\in\{0,1\}$ is a parameter so that $\alpha=0$ corresponds to the low Mach limit and $\alpha=1$ corresponds to the Boussinesq limit. Though the scheme (3.1)-(3.2) consists of a fully implicit step (3.1) and a semi-implicit step (3.2), its numerical resolution is fairly simple. Eliminating $q^{n+1}$ between (3.1) and (3.2) yields the nonlinear elliptic equation:

[TABLE]

where the known expression $\Phi$ is given by

[TABLE]

with $:$ denoting the contracted product. Solving the elliptic equation (3.3) yields the updated density $\rho^{n+1}$. The velocity $u^{n+1}$ can then be updated using (3.2), which is now an explicit evaluation. Hence, the scheme (3.1)-(3.2) consists of solving the elliptic equation (3.3), followed by an explicit evaluation of (3.2).

4. Asymptotic Preserving Property

A numerical scheme for a singular perturbation problem, such as the Euler system (2.1)-(2.2), may not resolve the existing multiple scales in space and time. In addition, when the perturbation parameter goes to zero, the scheme may approximate a completely different set of equations than the actual limiting systems. An asymptotic preserving (AP) scheme is the one which is consistent with the limiting set of equations in the singular limit; see [6] for a review of AP schemes.

Theorem 4.1**.**

The time semi-discrete scheme (3.1)-(3.2) for $\alpha=0$ is asymptotically consistent with the low Mach number model as $\varepsilon\to 0$.

Proof.

First, we apply the same ansatz (2.3) for all the dependent variables at times $t^{n}$ and $t^{n+1}$ in the semi-discrete scheme (3.1)-(3.2) and balance the like-powers of $\varepsilon$. The lowest order terms gives $\nabla P(\rho^{n+1}_{(0)})=0$ and the equation of state $P(\rho)=\rho^{\gamma}$ then yields that $\rho^{n+1}_{(0)}$ is constant. Therefore, from the mass update (3.1) we get

[TABLE]

We integrate the above equation (4.1) over a domain $\Omega$ and use Gauss’ divergence theorem to obtain:

[TABLE]

Hence, the leading order density $\rho_{(0)}$ rises or falls only due to compressions or expansions at the boundary. The temporal variations in $\rho_{(0)}$ can produce nonzero divergences in the leading order velocity $u_{(0)}$. It can be proved that the integral on the left hand side of (4.2) vanishes under most of the physically relevant boundary conditions. In this case, we obtain $\rho^{n+1}_{(0)}=\rho^{n}_{(0)}$ and this in turn enforces the divergence constraint at $t^{n+1}$ as

[TABLE]

Combining (4.3) and the $\mathcal{O}(1)$ terms in (3.2), we have the following limiting system:

[TABLE]

The above system (4.4)-(4.5) is clearly a consistent discretisation of the low Mach number limit system (2.4)-(2.5). ∎

Theorem 4.2**.**

The time semi-discrete scheme (3.1)-(3.2) for $\alpha=1$ is asymptotically consistent with the Boussinesq model.

Proof.

The proof is similar to that of Theorem 4.1 and hence omitted. ∎

5. $L^{2}$ Stability Analysis of the Semi-discrete Scheme

The aim of this section is to present the results of an $L^{2}$-stability analysis of the semi-discrete scheme (3.1)-(3.2). To this end, we consider the homogeneous linear wave equation system:

[TABLE]

as a simplified model of the Euler system (2.1)-(2.2). Here, $(\bar{\rho},\bar{u})$ is a linearisation state and $\bar{a}$ is a linearisation state for the sound velocity. Applying the AP methodology introduced in (3.1)-(3.2) to (5.1)-(5.2) yields the semi-discrete scheme:

[TABLE]

In the following, we use a stability result due to Richtmyer; see e.g. [7, 8] for details. Note that any difference scheme of the form $B_{1}U^{n+1}=B_{2}U^{n}$, where $B_{1},B_{2}$ are $p\times p$ matrices, independent of $t$ and $x$, and $U^{n}\in\mathbb{R}^{p}$ is the approximation to the original solution at time $t^{n}$, can be reduced to $\hat{U}^{n+1}=G(\Delta t,\xi)\hat{U}^{n}$ in the Fourier variable $\xi$. Here, $G(\Delta t,\xi)$ is the Fourier transform of the matrix $(B_{1})^{-1}B_{2}$ and is called the amplification matrix. The stability result due to Richtmyer states that

Theorem 5.1**.**

A difference scheme given by $B_{1}U^{n+1}=B_{2}U^{n}$ is stable if

1. (i)

the elements of $G(0,\xi)$ are bounded for all $\xi\in\mathbb{L}$, where $\mathbb{L}$ is a lattice where $\xi$ varies, 2. (ii)

$\lVert G(0,\xi)\rVert\leq 1$* and* 3. (iii)

$G(\Delta t,\xi)$* is Lipschitz continuous at $\Delta t=0$ in the sense that*

[TABLE]

Using the above theorem, we have the following stability result.

Theorem 5.2**.**

The semi-discrete scheme (5.3)-(5.4) is $L^{2}$-stable.

Proof.

Taking the Fourier transform of (5.3)-(5.4) and re-arranging the terms gives

[TABLE]

where

[TABLE]

Now, $G(0,\xi)$ reduces to the $3\times 3$ identity matrix and hence conditions (i) and (ii) of Theorem 5.1 are automatically satisfied. Further,

[TABLE]

Note that the matrix on the right hand side in (5.7) is bounded for every bounded lattice $\mathbb{L}$. Hence, by Theorem 5.1, the semi-discrete scheme (5.3)-(5.4) is $L^{2}$-stable. ∎

6. Numerical Experiments

We do not intend to discuss the space discretisation in detail as we use employ standard techniques. We use a finite volume approach to approximate the semi-discrete scheme (3.1)-(3.2). The explicit flux terms are approximated by a Rusanov-type flux whereas the implicit terms by simple central differences. The nonlinear elliptic equation (3.3) is solved iteratively after discretisation of the derivatives by central differences.

In the following, we consider a test problem in two dimensions to demonstrate the AP property of the scheme. We take the well-prepared initial data given in [1] which reads

[TABLE]

The computational domain $[0,1]\times[0,1]$ is divided into $50\times 50$ mesh points and we apply periodic boundary conditions on all four sides. The CFL number is set to 0.45 and we perform the computations up to a final time $T=1.0$. The parameter $\varepsilon$ is set to 0.1. Note that our CFL condition is independent of $\varepsilon$.

In Figures 1 and 2 we plot the density, $x_{1}$-velocity and the divergence of the velocity at times $t=0$ and $t=1$, for the low Mach and Boussinesq cases, respectively. It can be noted from the figures that in both the cases the density converges to the constant value 1 and the divergence approach 0. This is in conformity with the AP nature of the scheme in both the cases.

7. Conclusion

An AP semi-implicit time discretisation is proposed for the numerical approximation of the isentropic Euler equations with gravity in the low Mach number and Boussinesq limits. The schemes are theoretically shown to be asymptotically consistent as well as linearly stable. The results of numerical experiments provide a justification to AP nature of the scheme.

Acknowledgement

The authors thank Arnab Das Gupta for several useful discussions on the topic.