Consistent Internal Energy Based Schemes for the Compressible Euler Equations

TL;DR

This paper introduces consistent internal energy-based numerical schemes for the compressible Euler equations that ensure shock capturing, positivity, and entropy stability without Riemann solvers, applicable to various mesh types and discretizations.

Contribution

The paper develops new internal energy discretization schemes that are stable, positivity-preserving, and entropy-consistent, avoiding Riemann solvers and applicable to diverse mesh configurations.

Findings

Schemes ensure positivity of density, internal energy, and pressure.

Stability achieved without time step restrictions for pressure correction schemes.

Entropy inequality holds or is approximately satisfied, depending on the scheme.

Abstract

Numerical schemes for the solution of the Euler equations have recently been developed, which involve the discretisation of the internal energy equation, with corrective terms to ensure the correct capture of shocks, and, more generally, the consistency in the Lax-Wendroff sense. These schemes may be staggered or colocated, using either struc-tured meshes or general simplicial or tetrahedral/hexahedral meshes. The time discretization is performed by fractional-step algorithms; these may be either based on semi-implicit pressure correction techniques or segregated in such a way that only explicit steps are involved (referred to hereafter as "explicit" variants). In order to ensure the positivity of the density, the internal energy and the pressure, the discrete convection operators for the mass and internal energy balance equations are carefully designed; they use an upwind technique with respect to the material velocity only. The construction of the fluxes thus does not need any Rie-mann or approximate Riemann solver, and yields easily implementable algorithms. The stability is obtained without restriction on the time step for the pressure correction scheme and under a CFL-like condition for explicit variants: preservation of the integral of the total energy over the computational domain, and positivity of the density and the internal energy. The semi-implicit first-order upwind scheme satisfies a local discrete entropy inequality. If a MUSCL-like scheme is used in order to limit the scheme diffusion, then a weaker property holds: the entropy inequality is satisfied up to a remainder term which is shown to tend to zero with the space and time steps, if the discrete solution is controlled in L $\infty$ and BV norms. The explicit upwind variant also satisfies such a weaker property, at the price of an estimate for the velocity which could be derived from the introduction of a new stabilization term in the momentum balance. Still for the explicit scheme, with the above-mentioned MUSCL-like scheme, the same result only holds if the ratio of the time to the space step tends to zero.

Figures

Equations

$${\partial }_t\rho +\mathrm{div}(\rho \mathbf{u})=0,$$

$${\partial }_t(\rho \mathbf{u})+\mathrm{div}(\rho \mathbf{u}\otimes \mathbf{u})+\nabla p=0,$$

$${\partial }_t(\rho E)+\mathrm{div}(\rho E\mathbf{u})+\mathrm{div}(p\mathbf{u})=0,$$

$$p=(\gamma -1)\rho e,E=\frac{1}{2}|\mathbf{u}{|}^2+e,$$

$${\partial }_t\rho +\mathrm{div}(\rho \mathbf{u})=0.$$

$$\begin{array}[]{ll}\mathcal{C}(z)=\partial_{t}(\rho z)+{\rm div}(\rho z{\boldsymbol{u}})&=\rho\bigl{(}\partial_{t}z+{\boldsymbol{u}}\cdot{\boldsymbol{\nabla}}z\bigr{)}+z\bigl{(}\partial_{t}\rho+{\rm div}(\rho{\boldsymbol{u}})\bigr{)}\\[4.30554pt] &=\rho\bigl{(}\partial_{t}z+{\boldsymbol{u}}\cdot{\boldsymbol{\nabla}}z\bigr{)}.\end{array}$$

$$\varphi^{\prime}(z)\ \mathcal{C}(z)=\varphi^{\prime}(z)\,\rho\,\bigl{(}\partial_{t}z+{\boldsymbol{u}}\cdot{\boldsymbol{\nabla}}z\bigr{)}=\rho\,\bigl{(}\partial_{t}\varphi(z)+{\boldsymbol{u}}\cdot{\boldsymbol{\nabla}}\varphi(z)\bigr{)}.$$

$$\varphi^{\prime}(z)\ \mathcal{C}(z)=\partial_{t}\bigl{(}\rho\varphi(z)\bigr{)}+{\rm div}\bigl{(}\rho\varphi(z)\,{\boldsymbol{u}}\bigr{)}.$$

$$\frac{1}{\delta t}({\rho }^{n+1}-{\rho }^n)+\mathrm{div}({\rho }^{n+1}\mathbf{u})=0,$$

$$\varphi^{\prime}(z^{n+1})\ \bigl{[}\frac{1}{\delta t}\ (\rho^{n+1}z^{n+1}-\rho^{n}z^{n})+{\rm div}(\rho^{n+1}z^{n+1}{\boldsymbol{u}})\bigr{]}\\ =\frac{1}{\delta t}\ \Bigl{(}\rho^{n+1}\varphi(z^{n+1})-\rho^{n}\varphi(z^{n})\Bigr{)}+{\rm div}\bigl{(}\rho^{n+1}\varphi(z^{n+1})\,{\boldsymbol{u}}\bigr{)}+\mathcal{R}^{n},$$

$${\mathcal{R}}^n=\frac{1}{2\delta t}{\rho }^n{\phi }^{\prime \prime }(\bar{z})(z^{n+1}-z^n{)}^2,\bar{z}=\theta z^n+(1-\theta )z^{n+1},\theta \in [0,1].$$

$$\begin{array}[]{l}\displaystyle\frac{1}{\delta t}\ (\rho^{n+1}z^{n+1}-\rho^{n}z^{n})+{\rm div}(\rho^{n+1}z^{n+1}{\boldsymbol{u}})\\[8.61108pt] \displaystyle=\frac{1}{\delta t}\ \rho^{n}\ (z^{n+1}-z^{n})+\rho^{n+1}{\boldsymbol{u}}\cdot{\boldsymbol{\nabla}}z^{n+1}+z^{n+1}\Bigl{[}\frac{1}{\delta t}\ (\rho^{n+1}-\rho^{n})+{\rm div}(\rho^{n+1}{\boldsymbol{u}})\Bigr{]}\\[8.61108pt] \displaystyle=\frac{1}{\delta t}\ \rho^{n}\ (z^{n+1}-z^{n})+\rho^{n+1}{\boldsymbol{u}}\cdot{\boldsymbol{\nabla}}z^{n+1}.\end{array}$$

$$\partial_{t}(\rho E_{k})+{\rm div}\bigl{(}\rho\,E_{k}\,{\boldsymbol{u}}\bigr{)}+{\boldsymbol{\nabla}}p\cdot{\boldsymbol{u}}=0.$$

$${\partial }_t(\rho e)+\mathrm{div}(\rho e\mathbf{u})+p\mathrm{div}\mathbf{u}=0.$$

$$\frac{1}{\delta t}({\rho }^n{\tilde{\mathbf{u}}}^{n+1}-{\rho }^{n-1}{\mathbf{u}}^n)+\mathrm{div}({\rho }^n{\mathbf{u}}^n\otimes {\tilde{\mathbf{u}}}^{n+1})+{\zeta }^n\nabla p^n=0,$$

$$\frac{1}{\delta t}{\rho }^n({\mathbf{u}}^{n+1}-{\tilde{\mathbf{u}}}^{n+1})+\nabla p^{n+1}-{\zeta }^n\nabla p^n=0,$$

$$\frac{1}{\delta t}({\rho }^{n+1}-{\rho }^n)+\mathrm{div}({\rho }^{n+1}{\mathbf{u}}^{n+1})=0,$$

$$\frac{1}{\delta t}({\rho }^{n+1}{e}^{n+1}-{\rho }^n{e}^n)+\mathrm{div}({\rho }^{n+1}{e}^{n+1}{\mathbf{u}}^{n+1})+p^{n+1}\mathrm{div}{\mathbf{u}}^{n+1}=S^{n+1},$$

$$p^{n+1}=(\gamma -1){\rho }^{n+1}{e}^{n+1}.$$

$$\frac{1}{2\,\delta t}\,\bigl{(}\rho^{n}\,|\tilde{\boldsymbol{u}}^{n+1}|^{2}-\rho^{n-1}\,|{\boldsymbol{u}}^{n}|^{2}\bigr{)}+\frac{1}{2}{\rm div}\bigl{(}\rho^{n}\,|\tilde{\boldsymbol{u}}^{n+1}|^{2}{\boldsymbol{u}}^{n}\bigr{)}+\zeta^{n}{\boldsymbol{\nabla}}p^{n}\cdot\tilde{\boldsymbol{u}}^{n+1}+R^{n}_{1}=0,$$

$$R_1^n=\frac{1}{2\delta t}|{\tilde{\mathbf{u}}}^{n+1}-{\mathbf{u}}^n{|}^2.$$

$$\alpha^{n}{\boldsymbol{u}}^{n+1}+\frac{1}{\alpha^{n}}{\boldsymbol{\nabla}}p^{n+1}=\alpha^{n}\tilde{\boldsymbol{u}}^{n+1}+\frac{\zeta^{n}}{\alpha^{n}}{\boldsymbol{\nabla}}p^{n},\quad\alpha^{n}=\bigl{[}\frac{\rho^{n}}{\delta t}\bigr{]}^{1/2}$$

$$\frac{1}{2\delta t}{\rho }^n|{\mathbf{u}}^{n+1}{|}^2+\nabla p^{n+1}\cdot {\mathbf{u}}^{n+1}+R_2^n=\frac{1}{2\delta t}{\rho }^n|{\tilde{\mathbf{u}}}^{n+1}{|}^2+{\zeta }^n\nabla p^n\cdot {\tilde{\mathbf{u}}}^{n+1},$$

$$R_2^n=\frac{\delta t}{{\rho }^n}|\nabla p^{n+1}{|}^2-({\zeta }^n{)}^2\frac{\delta t}{{\rho }^n}|\nabla p^n{|}^2.$$

$$\frac{1}{2\,\delta t}\,\bigl{(}\rho^{n}\,|{\boldsymbol{u}}^{n+1}|^{2}-\rho^{n-1}\,|{\boldsymbol{u}}^{n}|^{2}\bigr{)}+\frac{1}{2}{\rm div}\bigl{(}\rho^{n}\,|\tilde{\boldsymbol{u}}^{n+1}|^{2}{\boldsymbol{u}}^{n}\bigr{)}+{\boldsymbol{\nabla}}p^{n+1}\cdot{\boldsymbol{u}}^{n+1}+R^{n}_{1}+R^{n}_{2}=0.$$

$$\zeta^{n}=\bigl{[}\frac{\rho^{n}}{\rho^{n-1}}\bigr{]}^{1/2}.$$

$$\frac{|K|}{\delta t}({\rho }_K^{n+1}-{\rho }_K^n)+\sum _{\sigma \in \mathcal{E}(K)}{F}_{K,\sigma }^{n+1}=0,$$

$${\mathcal{C}}_K(e^{n+1})=\frac{|K|}{\delta t}({\rho }_K^{n+1}{e}_K^{n+1}-{\rho }_K^n{e}_K^n)+\sum _{\sigma \in \mathcal{E}(K)}{F}_{K,\sigma }{e}_{\sigma }^{n+1},$$

$$\frac{1}{\delta t}({\rho }^{n+1}-{\rho }^n)+\mathrm{div}({\rho }^n{\mathbf{u}}^n)=0,$$

$$\frac{1}{\delta t}({\rho }^{n+1}{e}^{n+1}-{\rho }^n{e}^n)+\mathrm{div}({\rho }^n{e}^n{\mathbf{u}}^n)+p^n\mathrm{div}{\mathbf{u}}^n=S^n,$$

$$p^{n+1}=(\gamma -1){\rho }^{n+1}{e}^{n+1},$$

$$\frac{1}{\delta t}({\rho }^{n+1}{\mathbf{u}}^{n+1}-{\rho }^n{\mathbf{u}}^n)+\mathrm{div}({\rho }^n{\mathbf{u}}^n\otimes {\mathbf{u}}^n)+\nabla p^{n+1}=0.$$

$${\partial }_t\rho +\mathrm{div}(\rho \mathbf{u})=0,$$

$${\partial }_t(\rho e)+\mathrm{div}(\rho e\mathbf{u})+p\mathrm{div}(\mathbf{u})=\mathcal{R}\ge 0,$$

$$p=(\gamma -1)\rho e,$$

$$\partial_{t}\eta(\rho,e)+{\rm div}\bigl{[}\eta(\rho,e)\,{\boldsymbol{u}}\bigr{]}\leq 0.$$

Full text

1 Introduction

We address in this paper the solution of the Euler equations for an ideal gas, which read:

[TABLE]

where $t$ stands for the time, $\rho$, ${\boldsymbol{u}}$, $p$, $E$ and $e$ are the density, velocity, pressure, total energy and internal energy respectively, and $\gamma>1$ is a coefficient specific to the considered fluid. The problem is supposed to be posed over $\Omega\times(0,T)$, where $\Omega$ is an open bounded connected subset of ${\mathbb{R}}^{d}$, $1\leq d\leq 3$, and $(0,T)$ is a finite time interval. System (1) is complemented by initial conditions for $\rho$, $e$ and ${\boldsymbol{u}}$, let us say $\rho_{0}$, $e_{0}$ and ${\boldsymbol{u}}_{0}$ respectively, with $\rho_{0}>0$ and $e_{0}>0$, and by suitable boundary conditions which we suppose to be ${\boldsymbol{u}}\cdot{\boldsymbol{n}}=0$ at any time and a.e. on $\partial\Omega$, where ${\boldsymbol{n}}$ stands for the normal vector to the boundary.

Finite volume schemes for the solution of hyperbolic problems such as the system (1) generally use a collocated arrangement of the unknowns, which are associated to the cell centers, and apply a Godunov-like technique for the computation of the fluxes at the cells faces: the face is seen as a discontinuity line for the beginning-of-time-step numerical solution, supposed to be constant in the two adjacent cells; the value of the solution of the so-posed Riemann problem on the discontinuity line is computed, either exactly or approximately; the numerical solution at the end-of-time-step is computed with these values, and is a piecewise constant function (see e.g. [39, 3] for the development of such solvers). In one space dimension, this method consists, at least for exact Riemann solvers, in a projection of the exact solution. Then, thanks to the properties of the projection, this process applied to the Euler equations yields consistent schemes which preserve the non-negativity of the density and the internal energy and, for first-order variants, satisfy an entropy inequality. The price to pay is the computational cost of the evaluation of the fluxes, and the fact that this issue is intricate enough to put almost out of reach implicit-in-time formulations, which would allow to relax CFL time step constraints. In addition, preserving the accuracy for low Mach number flows is a difficult task (see e.g. [18] and references herein).

The aim here is first to review some recent schemes which follow a different route, and then prove some discrete entropy estimates and/or consistency results for these schemes. The space discretization may be colocated [25] or staggered [21, 17, 23]: in the colocated case, all unknowns are located at the center of the discretization cells, while in the staggered case, scalar variables are associated to cell centers while the velocity is associated to the faces, or, equivalently, to staggered mesh(es). The use of staggered discretization for compressible flows goes back to the MAC scheme [19], and has been the subject of a wide litterature (see [42] for a textbook and references in [21, 17, 23]). Staggered discretizations have been preferred in the open source CALIF3S [4] used for nuclear safety applications because the resulting semi-implicit schemes are asymptotically stable in the low Mach number regime [24]. Two different staggered space discretizations may be considered: either the so-called Marker-And-Cell (MAC) scheme for structured grids [20] or, for general meshes, a space discretization using degrees of freedom similar to low-order Rannacher-Turek [34] or Crouzeix-Raviart [8] finite elements (see Figure 1). With this space discretization, the use of Riemann solvers seems difficult (scalar unknowns and velocities may still be considered as piecewise constant functions, but not associated to the same partition of the computational domain). The positivity of the internal energy is thus ensured by a non-standard argument: the internal energy balance is discretized instead of the actual (total) energy balance (1c) by a positivity-preserving scheme. This strategy is known to lead to consistency problems (wrong shock speeds for instance), which are circumvented by some correction terms in the discrete internal energy correction. Until now, the use of the internal energy equation associated to a consistency correction seems to be restricted to the context of Lagrangian approaches, up to a very recent work implementing a Lagrange-remap technique on staggered meshes [9], and some recent developements extending the techniques developped here to more genral meshes [31]. Two time discretizations are proposed: a pressure correction technique and a segregated scheme involving only explicit steps. The resulting schemes offer many interesting properties: both the density and internal energy positivity are preserved, unconditionally for the pressure correction scheme and under CFL-like conditions for the segregated explicit variant, and the integral of the total energy on the computational domain is conserved (which yields a stability result); the construction of the fluxes simply relies on standard upwinding techniques of the convection operators with respect to the material velocity; finally, the space approximation, the fluxes and the choice of the internal energy balance are consistent with usual discretizations of quasi-incompressible flows, so the pressure correction scheme is asymptotic preserving by construction in the limit of vanishing low Mach number flows (see [24] for a study in the case of the barotropic Euler equations).

In addition, a discrete entropy estimate is obtained for the (upwind) pressure correction scheme, while only a conditional weak entropy estimate seems to hold for the segregated explicit variant. Note that the schemes studied here belong to a class often referred to as ”flux splitting schemes” in the literature, since they may be obtained by splitting the system by a two-step technique (usually into a ”convective” and ”acoustic” part), applying a standard scheme to each part (which, for the convection system, indeed yields, at first order, an upwinding with respect to the material velocity) and then summing both steps to obtain the final flux. Works in this direction may be found in [37, 30, 43, 29, 40], and we hope that the discussion presented on the entropy may be extended in some way to these numerical methods.

The paper is organised as follows; Section 2 is devoted to the derivation of the previously mentioned internal energy based schemes in the semi-discrete time setting. Section 3 presents some new and original results concerning some entropy estimates and/or entropy consistency which hold for both the colocated and staggered schemes, first for implicit schemes, and then for explicit schemes.

2 Derivation of the numerical schemes

2.1 A basic result on convection operators

Let $\rho$ and ${\boldsymbol{u}}$ be regular respectively scalar and vector-valued functions such that

[TABLE]

Let $z$ be a regular scalar function. Then

[TABLE]

Let $\varphi$ be a regular real function. Then:

[TABLE]

Now, reversing the computation performed in Relation (2) with $\varphi(z)$ instead of $z$ leads to

[TABLE]

The following lemma states a time semi-discrete version of this computation.

Lemma 1 (Convection operator)

Let $\rho^{n}$, $\rho^{n+1}$, $z^{n}$ and $z^{n+1}$ be regular scalar functions, let ${\boldsymbol{u}}$ be a regular vector-valued function and let $\varphi$ be a twice-differentiable real function. Let us suppose that

[TABLE]

with $\delta t$ a positive real number. Then

[TABLE]

with

[TABLE]

Proof

We first begin by deriving a discrete analogue to Identity (2):

[TABLE]

Then the result follows by multiplying this relation by $\varphi^{\prime}(z^{n+1})$, using a Taylor expansion for the first term and the same combination of partial derivative as in the continuous case for the second term, and finally, still as in the continuous cas, by performing this computation in the reverse sense with $\varphi(z^{n})$ and $\varphi(z^{n+1})$ instead of $z^{n}$ and $z^{n+1}$.

2.2 Internal energy formulation

We begin with a formal reformulation of the energy equation. Let us suppose that the solution is regular, and let $E_{k}$ be the kinetic energy, defined by $E_{k}=\frac{1}{2}\,|{\boldsymbol{u}}|^{2}$. Taking the inner product of (1b) by ${\boldsymbol{u}}$ yields, after the formal compositions of partial derivatives described in the previous section:

[TABLE]

This relation is referred to as the kinetic energy balance. Subtracting this relation to the total energy balance (1c), we obtain the so-called internal energy balance equation:

[TABLE]

Since,

as seen in the previous section, thanks to the mass balance equation, the first two terms in the left-hand side of (8) may be recast as a transport operator,

• -

and, from the equation of state, the pressure vanishes when $e=0$,

this equation implies that, if $e\geq 0$ at $t=0$ and with suitable boundary conditions, then $e$ remains non-negative at all time. The same result would hold if (8) featured a non-negative right-hand side, as for the compressible Navier-Stokes equations. Solving the internal energy balance (8) instead of the total energy balance(1c) is thus appealing, to preserve this positivity property by construction of the scheme. In addition, it avoids introducing a space discretization for the total energy which, for a staggered discretization, combines cell-centered (internal energy and density) and face-centered (velocity) variables. However, a raw discretization of a non-conservative equation derived from a conservative system (formally, i.e. supposing unrealistic regularity properties of the solution) may be non-consistent (and the numerical test presented in Section 2.6 shows that, for the problem at hand, such a scheme would be unable to capture shock solutions). To deal with this problem, we implement the following strategy:

First, we derive a discrete kinetic energy balance, by mimicking at the discrete level the computation leading to Equation (7), so as to identify the terms which are likely to lead to non-consistency: the numerical diffusion in the momentum balance equation yields dissipation terms in the kinetic energy balance which are observed to behave, when the space and time step tend to zero, as measure born by the shocks which modify the jump conditions.

• -

These terms are thus compensated in the internal energy balance.

At the fully discrete level, for staggered discretizations, the kinetic and internal energy balances are not posed on the same mesh (the dual and primal mesh respectively); however, it is possible to derive from the kinetic energy balance on the dual mesh a counterpart posed on the primal mesh and, adding to the internal energy balance yields a conservative total energy balance. The scheme then can be proven to be consistent in the Lax-Wendroff sense to the weak form of the total energy balance: for a given sequence of discrete solutions (obtained with a sequence of discretizations whith space and time steps tending to zero) controlled and converging to a limit in suitable norms (namely, uniformly bounded and converging in $L^{r}$ norms, for $r\in[1,+\infty)$), we show that the limit is a weak solution to the Euler equations [23, 22]. In the colocated case, both kinetic and internal energy balances are posed on the same mesh and a discrete local total energy balance is easily recovered [25].

2.3 The time semi-discrete pressure correction scheme

This semi-discrete pressure correction scheme takes the following general form:

[TABLE]

Solving the first equation yields a tentative velocity $\tilde{\boldsymbol{u}}^{n+1}$; this is the velocity prediction step, which is decoupled from the other equations of the system. Equations (9b)-(9e) constitute the correction step and are solved simultaneously; in the relation (9d), the term $\rho^{n+1}\,e^{n+1}$ is recast as a function of the pressure only thanks to the equation of state (1d) and the velocity ${\boldsymbol{u}}^{n+1}$ is eliminated thanks to the divergence of (9b) divided by $\rho^{n}$. The result is a nonlinear and nonconservative elliptic problem for the pressure only. This process must be performed at the fully discrete level to preserve the properties of the scheme. The coefficient $\zeta^{n}$ in Equation (9a) and the correction term $S^{n+1}$ in (9d) are chosen so as to ensure stability and consistency, as shown below. The first step of this process is to obtain a discrete kinetic energy balance. To this purpose, let us multiply (9a) by $\tilde{\boldsymbol{u}}^{n+1}$ and apply Lemma 1 component by component, with $\varphi(s)=\frac{1}{2}s^{2}$. We get:

[TABLE]

with

[TABLE]

Note that the mass balance equation (9c), which is a fundamental assumption in Lemma 1, only holds at this stage of the algorithm with the previous time step values, hence the shift of the time level of the density in (9a). Let us now recast Equation (9b) as

[TABLE]

and square this relation, to get

[TABLE]

with

[TABLE]

Summing (10) and (11) yields the kinetic energy balance that we are seeking:

[TABLE]

The coefficient $\zeta^{n}$ is then chosen in such a way that the remainder term $R^{n}_{2}$ is a difference of two consecutive time levels of the same quantity; this is the case for

[TABLE]

Supposing the control in $L^{1}(0,T,BV)$ of the pressure and in $L^{\infty}$ of the pressure and of the inverse of the density, the term $R^{n}_{2}$ may thus be seen to tend with zero with the discretization parameters in a distributional sense. The term $R^{n}_{1}$ is compensated in the internal energy balance, by choosing $S^{n+1}=R^{n}_{1}$, thus ensuring that $S^{n+1}\geq 0$. The definition of the time-discrete scheme is now complete.

2.4 The fully discrete pressure correction scheme

The fully discrete scheme is obtained from System (9) by applying the following guidelines:

The mass and internal energy balances (i.e. Equations (9c) and (9d) respectively) are discretized on the primal mesh, while the velocity prediction (9a) and correction (9b) are discretized on the dual mesh(es). The equation of state only involves cell quantities, and its expression is obtained by writing (9e) for these latter.

• -

The space arrangement of the unknowns (density discretized at the cell and velocity at the faces) yields a natural expression of the mass fluxes in the mass balance, performed by a first-order upwind scheme (with respect to the velocity). By construction, the density is thus non-negative; in fact at the discrete level, it remains positive if the initial density is positive. The discrete mass balance equation on the cell $K$ whose measure is denoted by $|K|$ takes the form:

[TABLE]

where ${\mathcal{E}}(K)$ denotes the set of edges of $K$ and $F_{K,\sigma}$ is the mass flux across $\sigma$ outward $K$.

• -

Let $\mathcal{C}_{K}(e^{n+1})$ denote the sum of the discrete time-derivative and convection operator in the internal energy balance (9d); this quantity reads:

[TABLE]

where $e^{n+1}_{\sigma}$ is the upwind approximation of $e^{n+1}$ at $\sigma$ with respect to $F^{n+1}_{K,\sigma}$ (or, equivalently, since the density is positive, with respect to the velocity). The structure of $\mathcal{C}_{K}(e^{n+1})$ (precisely speaking, the fact that $\mathcal{C}_{K}(e^{n+1})$ vanishes thanks to the mass balance if the internal energy $e^{n+1}$ is constant over $\Omega$) was shown in [27] to yield a positivity-preserving operator, and is also a necessary condition for a fully discrete version of Lemma 1 to hold; this is of course linked since both results rely on the possibility to recast $\mathcal{C}_{K}$ as a transport operator, and the positivity-preserving property of $\mathcal{C}_{K}$ may be proved by applying Lemma 1 with $\varphi(s)=\min(s,0)^{2}$. Once again, thanks to the arrangement of the unknowns, a natural discretization for ${\rm div}{\boldsymbol{u}}^{n+1}$ is available. Since $p^{n+1}$ is a function of $e^{n+1}$ (given by the equation of state) which vanishes for $e^{n+1}=0$ and since the corrective term is non-negative, we are able to show that the discrete internal energy is kept positive by the scheme.

• -

To allow to derive a discrete kinetic energy balance, the same structure is needed for the time-derivative and convection operator in the velocity prediction step (9a). This raises a difficulty since this equation is posed on the dual mesh, and thus we need an analogue of the mass balance (12) to also hold on this mesh. The way to build the face density and the mass fluxes across the faces of the dual mesh for such a relation to hold, while still ensuring the scheme consistency, is a central ingredient of the scheme; it is detailed in [13] for the MAC discretization and in [28] for unstructured discretizations.

Once the face density is defined, the discretization of the coefficient $\zeta^{n}$ is straightforward. In order to combine the discrete equivalents of ${\boldsymbol{u}}\cdot{\boldsymbol{\nabla}}p$ (kinetic energy balance) and $p\,{\rm div}{\boldsymbol{u}}$ (internal energy balance), the discrete gradient is defined as the transposed of the divergence operator with respect to the $L^{2}$ inner product (if ${\boldsymbol{u}}\cdot{\boldsymbol{\nabla}}p+p\,{\rm div}{\boldsymbol{u}}={\rm div}(p\,{\boldsymbol{u}})$, the integral of this quantity over the computational domain vanishes when the normal velocity is prescribed to zero at the boundary). Note that this definition is consistent with the usual treatment in the incompressible case, and is a key ingredient for the scheme to be asymptotic preserving in the limit of vanishing Mach number flows [24]. As in the incompressible case, it also allows to control the $L^{2}$ norm of the pressure by a weak norm of its gradient, which is central for convergence studies; with this respect, a discrete inf-sup condition is required in some sense, which is true for staggered discretizations.

2.5 A segregated variant

A variant of the proposed scheme which consists only in explicit steps (in the sense that these steps do not require the solution of any linear or non-linear algebraic system) reads, in the time semi-discrete setting:

[TABLE]

The update of the pressure before the solution of the momentum balance equation is crucial in our derivation of entropy estimates (see Section 3 below). This issue seems to be supported by numerical experiments: omitting it, we observe the appearance of non-entropic discontinuities in rarefaction waves [23].

The space discretization differs from the pressure correction scheme described in the above section in two points:

the discretization of the convection operator in the momentum balance equation (13d) is performed by the first order upwind scheme (still with respect to the material velocity ${\boldsymbol{u}}^{n}$),

• -

the corrective term $S^{n}$ is still obtained by deriving a kinetic energy balance multiplying Equation (13d) by ${\boldsymbol{u}}^{n+1}$, but its expression is quite different, due to the time-level used in the convection operator. The time-discretization is now anti-diffusive but, as usual for explicit schemes, this anti-diffusion is counterbalanced by the diffusion in the approximation of the convection (hence the upwinding) and $S^{n}$ is non-negative only under a CFL condition.

2.6 A numerical test

In this section, we reproduce a test performed in [21] to assess the behaviour of the scheme on a one dimensional Riemann problem. We choose initial conditions such that the structure of the solution consists in two shock waves, separated by the contact discontinuity, with sufficiently strong shocks to allow an easy discrimination of correct numerical solutions. These initial conditions are those proposed in [39, chapter 4], for the test referred to as Test 5. The computations are performed with the open-source software CALIF3S [4].

The density fields obtained with $h=1/2000$ (or a number of cells $n=2000$) at $t=0.035$, with and without assembling the corrective source term in the internal energy balance, together with the analytical solution, are shown on Figure 2. We observe that both schemes seem to converge, but the corrective term is necessary to obtain the right solution. Without a corrective term, one can check that the obtained solution is not a weak solution to the Euler system (Rankine-Hugoniot conditions are not verified). We also observe that the scheme is rather diffusive especially at contact discontinuities for which the beneficial compressive effect of the shocks does not apply; this may be cured in the explicit variant by implementing MUSCL-like algorithms [14].

Extensive multidimensional tests were performed in both the staggered case [17] and the colocated case [25].

3 Entropy

In the case of regular solutions to the Euler equations (1), an additional conservation law can be written for an additional quantity called entropy; however, in the presence of shock waves, the (mathematical) entropy decreases. It is now known that weak solutions of the Euler system satisfying an entropy inequality may be non unique [5]; nevertheless, entropy inequalities play an important role in providing global stability estimates [3].

When solving the Euler equations numerically, it is thus natural to design numerical schemes such that some entropy inequalities are satisfied by the approximate solutions; these inequalities should enable to prove that, as the mesh and time steps tend to 0, the limit of the approximate solutions, if it exists, satisfies an entropy inequality. A classical way of doing so is to design so-called “entropy stable schemes” [38]. Discrete entropy inequalities are known for the one dimensional case for the Godunov scheme [15] and have been derived for Roe-type schemes in the one space dimension case [26]. Entropy stability has also been proven in the multi-dimensional case for semi-discrete schemes on unstructured meshes [32, 35]. In the sequel we show that an implicit upwind scheme (at least with the upwinding with respect to the material velocity used here) is indeed entropy stable. However it is not always possible to obtain entropy stability, especially for fully discrete schemes such as the explicit schemes studied below; in this case, weaker discrete entropy inequalities or estimates are obtained which allow to fulfil our goal, namely to show that the possible limits of the approximate solutions satisfy an entropy inequality. Such a technique was used for the convergence study of a time implicit mixed finite volume–finite element scheme for the Euler-Fourier equations [11], with a special equation of state which allows to obtain a priori estimates.

Both the pressure correction scheme (9) and the segregated scheme (13) involve a discrete equivalent of the following subsystem:

[TABLE]

with the same initial and boundary conditions as for the full system (1).

The derivation of an entropy for the continuous Euler system may be deduced from the subsystem (14) in the following way. We seek an entropy function $\eta$ satisfying:

[TABLE]

To this end, we introduce the functions $\varphi_{\rho}$ and $\varphi_{e}$ defined as follows:

[TABLE]

For regular functions, the function $\eta$ defined by

[TABLE]

satisfies (15). Indeed, multiplying (14a) by $\varphi^{\prime}_{\rho}(\rho)$, a formal computation yields:

[TABLE]

Then, multiplying (14b) by $\varphi^{\prime}_{e}(e)$ yields, once again formally, since $\varphi^{\prime}_{e}(z)<0$ for $z>0$:

[TABLE]

Summing (18) and (19) and noting that $\varphi_{\rho}$ and $\varphi_{e}$ are chosen such that

[TABLE]

we obtain (15), which is an entropy balance for the Euler equations, for the specific entropy defined by (17).

In the sequel we derive some analogous discrete entropy inequalities (with a possible remainder tending to 0) for the fully discrete, time semi-implicit (and fully implicit, i.e. backward Euler, as far as System (14) only is concerned) or segregated schemes (fully explicit regarding System (14) only) presented in Section 2, with a possible upwinding limited to that of the convection terms with respect to the material velocity. Note that the entropy inequalities that we obtain here apply to both the staggered schemes [21, 17, 14] and to the colocated scheme [25] which is also based on the internal energy; indeed the entropy depends only on the mass and internal energy which are scalar unknowns located at the center of the (primal) cells in both schemes, so System (14) involves only equations posed on the primal mesh.

Depending on the time and space discretization, we obtain three types of results:

local entropy estimates, i.e. discrete analogues of (15), in which case the scheme is entropy stable,

• -

global entropy estimates, i.e. discrete analogues of:

[TABLE]

such a relation is a stability property of the scheme; this kind of relation was also proven in e.g. [7] for a higher order scheme for the 1D Euler equations;

• -

“weak local” entropy inequalities, i.e. results of the form:

[TABLE]

with $\widetilde{\mathcal{R}}$ tending to zero in a suitable sense with respect to the space and time discretization steps (or combination of both parameters), provided that the approximate solutions are controlled in reasonable norms, here, $L^{\infty}$ and BV norms. Then a ”Lax-consistency” property holds, of the form: a limit $(\bar{\rho},\bar{\boldsymbol{u}},\bar{e})$ of a convergent subsequence of approximate solutions given by the considered numerical scheme and bounded in the $L^{\infty}$ and BV norms, satisfies the following weak entropy inequality:

[TABLE]

In the sequel we address implicit schemes (Section 3.2) and segregated explicit schemes (Section 3.3). For implicit schemes, we first consider an upwind discretization for which we get a local discrete entropy inequality (Theorem 3.1), and then a MUSCL-like improvement of the discretization of the convection term in order to reduce the numerical diffusion, for which we only get a global entropy estimate and a weak local entropy inequality (Theorem 3.2). The case of explicit schemes is a little more tricky: we again consider the same two discretizations (i.e. upwind and MUSCL-like) but we first deal with the mass balance equation, then with the internal energy equation, and combine the results to address entropy inequalities.

3.1 Meshes and discrete norms

Let ${\mathcal{M}}$ be a mesh of the domain $\Omega$, supposed to be regular in the usual sense of the finite element literature (see e.g. [6]). By ${\mathcal{E}}$ and ${\mathcal{E}}(K)$ we denote the set of all $(d-1)$-faces $\sigma$ of the mesh and of the cell $K\in{\mathcal{M}}$ respectively, and we suppose that the number of the faces of a cell is bounded. The set of faces included in $\Omega$ (resp. in the boundary $\partial\Omega$) is denoted by ${\mathcal{E}}_{{\rm int}}$ (resp. ${\mathcal{E}}_{{\rm ext}}$); a face $\sigma\in{\mathcal{E}}_{{\rm int}}$ separating the cells $K$ and $L$ is denoted by $\sigma=K|L$. For $K\in{\mathcal{M}}$ and $\sigma\in{\mathcal{E}}$, we denote by $|K|$ the measure of $K$ and by $|\sigma|$ the $(d-1)$-measure of the face $\sigma$. The following quantities related to the mesh are used in the sequel:

[TABLE]

Let $(t_{n})_{0\leq n\leq N}$, with $0=t_{0}<t_{1}<\ldots<t_{N}=T$, define a partition of the time interval $(0,T)$, which we suppose uniform for the sake of simplicity, and let $\delta t=t_{n+1}-t_{n}$ for $0\leq n\leq N-1$ be the (constant) time step.

The discrete pressure, density and the internal energy unknowns are associated with the cells of the mesh ${\mathcal{M}}$; they are denoted by:

[TABLE]

In the estimates given below, we shall need some discrete norms that we now define.

Definition 1 (Discrete BV semi-norms and weak $L^{1}(0,T;(W^{1,+\infty}_{0})^{\prime})$ norm)

For a family $(z_{K}^{n})_{K\in{\mathcal{M}},0\leq n\leq N}\subset{\mathbb{R}}$, let us define the following norms of the associated piecewise constant function $z$:

[TABLE]

where $\psi_{K}^{n}$ stands for $\psi({\boldsymbol{x}}_{K},t_{n})$, with ${\boldsymbol{x}}_{K}$ the mass center of $K$. Note that this latter weak norm is the discrete equivalent of the continuous dual norm of $v\in L^{1}(\Omega\times(0,T))$, defined by

[TABLE]

Some of the proofs below are based on the following convexity result [12, Lemma 2.3]. In its formulation, and throughout the paper, $|[a,\ b]|$ stands for the interval $[\min(a,b),\ \max(a,b)]$, for any real numbers $a$ and $b$.

Lemma 2

Let $\varphi$ be a strictly convex and continuously differentiable function over an open interval $I$ of ${\mathbb{R}}$. Let $x_{K}\in I$ and $x_{L}\in I$ be two real numbers. Then the relation

[TABLE]

uniquely defines the real number $x_{KL}$ in $|[x_{K},x_{L}]|$.

Remark 1 ($x_{KL}$ for $\varphi(z)=z^{2}$)

Let us consider the specific function $\varphi(z)=z^{2}$. Then, an easy computation yields $x_{KL}=(x_{K}+x_{L})/2$ i.e. the centered approximation.

3.2 Implicit schemes

With the above notations, the space time discretization of System (14) reads:

[TABLE]

where $F_{K,\sigma}^{n+1}$ is the mass flux through the face $\sigma$, $e_{\sigma}^{n+1}$ is an approximation of the internal energy at the face $\sigma$, and $u_{K,\sigma}^{n+1}$ stands for an approximation of the normal velocity to the face $\sigma$; note that the velocity is solved in the full scheme by a space discretization of the momentum prediction and correction equations (9a)-(9b). Consistently with the boundary conditions, $u_{K,\sigma}^{n+1}$ vanishes on every external face. The mass flux $F_{K,\sigma}^{n+1}$ reads:

[TABLE]

where $\rho_{\sigma}^{n+1}$ stands for an approximation of the density on $\sigma$. Throughout the paper, we suppose that $\rho_{K}^{n}$, $e_{K}^{n}$, $\rho_{\sigma}^{n}$ and $e_{\sigma}^{n}$ are positive, for any $K\in{\mathcal{M}}$, $\sigma\in{\mathcal{E}}_{{\rm int}}$, $0\leq n\leq N$, which is verified by the solutions of the schemes presented in [21, 23, 14, 16] (of course, with positive initial conditions for $\rho$ and $e$).

The two following lemmas are straigthforward consequences of Lemmas A1 and A2 in [21] and state discrete analogues of (18) and (19) respectively which are used to obtain the entropy inequalities.

Lemma 3

Let $K\in{\mathcal{M}}$, $n$ be such that $0\leq n\leq N-1$ and let us suppose that the discrete mass balance (27a) holds. Let $\varphi$ be a twice continuously differentiable function defined over $(0,+\infty)$. Then

[TABLE]

where $\rho_{K}^{n+1/2}\in|[\rho_{K}^{n},\rho_{K}^{n+1}]|.$

Lemma 4

Let $K\in{\mathcal{M}}$ and $n$ be such that $0\leq n\leq N-1$. Let $\varphi$ be a twice continuously differentiable function defined over $(0,+\infty)$. Then:

[TABLE]

where $e_{K}^{n+1/2}\in|[e_{K}^{n},e_{K}^{n+1}]|$.

3.2.1 Upwind implicit schemes –

In this section, we suppose that the convection fluxes are approximated with a first order upwind scheme, i.e., for $\sigma\in{\mathcal{E}}_{{\rm int}}$, $\sigma=K|L$, $\rho_{\sigma}^{n+1}=\rho_{K}^{n+1}$ and $e_{\sigma}^{n+1}=e_{K}^{n+1}$ if $u_{K,\sigma}\geq 0$, $\rho_{\sigma}^{n+1}=\rho_{L}^{n+1}$ and $e_{\sigma}^{n+1}=e_{L}^{n+1}$ otherwise. In this case, the scheme (27) satisfies a local entropy estimate (i.e. a discrete analogue of Inequality (15)) which is stated in Theorem 3.1 below. Of course, this local entropy inequality also yields the global discrete inequality analogue to (21); furthermore, passing to the limit on the upwind implicit (or pressure correction) scheme applied the Euler equations, this local estimate also yields the Lax consistency, i.e. any limit $(\bar{\rho},\bar{\boldsymbol{u}},\bar{e})$ of a convergent subsequence of approximate solutions satisfies the weak entropy inequality (22).

Theorem 3.1 (Discrete entropy inequality, implicit upwind scheme)

Let $\eta$ be defined by (17), and, for $0\leq n\leq N-1$, let $\eta_{K}^{m}=\eta(\rho_{K}^{m},e_{K}^{m})$ for $m=n,\ n+1$ and $K\in{\mathcal{M}}$, and $\eta_{\sigma}^{n+1}=\eta(\rho_{\sigma}^{n+1},e_{\sigma}^{n+1})$ for $\sigma\in{\mathcal{E}}_{{\rm int}}$. Then any solution of the scheme (27) satisfies, for any $K\in{\mathcal{M}}$ and $0\leq n\leq N-1$:

[TABLE]

Proof

Let $\varphi$ be a twice continuously differentiable function. By Lemma 3, we get that (29) holds. For $\sigma\in{\mathcal{E}}_{{\rm ext}}$, thanks to the boundary conditions, the convection fluxes vanish. For $0\leq n\leq N-1$ and $K\in{\mathcal{M}}$, consider the term $(T_{m})_{K,\sigma}^{n+1}$ associated to an internal face $\sigma=K|L$ in the remainder term $(R_{m})_{K}^{n+1}$:

[TABLE]

where $\rho_{\sigma,K}^{n+1}\in|[\rho_{\sigma}^{n+1},\rho_{K}^{n+1}]|$. With the upwind choice, if $u_{K,\sigma}^{n+1}\geq 0$, $\rho_{\sigma}^{n+1}=\rho_{K}^{n+1}$ and $(T_{m})_{K,\sigma}^{n+1}$ vanishes. If $u_{K,\sigma}^{n+1}<0$ and $\varphi^{\prime\prime}$ is a non-negative function (i.e. $\varphi$ is convex), $(T_{m})_{K,\sigma}^{n+1}$ is non-negative and so is $(R_{m})_{K}^{n+1}$, for any $K\in{\mathcal{M}}$. Since $\varphi_{\rho}$ defined by (16) is indeed convex, Lemma 3 implies that any solution $\{\rho_{K}^{n},K\in{\mathcal{M}},0\leq n\leq N\}$ to Equation (27a) of the scheme satisfies, for $K\in{\mathcal{M}}$ and $0\leq n\leq N-1$:

[TABLE]

Now turning to Lemma 4, by similar arguments, the remainder term $(R_{e})_{K}^{n+1}$ in (30) is nonnegative for any regular convex function $\varphi$, for any $K\in{\mathcal{M}}$ and $0\leq n\leq N-1$. Hence, since $\varphi_{e}$ defined by Equation (16) is convex, we get that any solution to (27b) satisfies:

[TABLE]

The desired relation is then obtained by summing the inequalities (31) and (32), using (20).

3.2.2 MUSCL-like schemes –

The aim of this section is to improve the approximation of the convection fluxes in (27a) and (27b) in order to reduce the numerical diffusion, while still satisfying an entropy inequality. This leads to a condition similar to the limitation procedure which is the core of a MUSCL procedure [41]; indeed, in order to yield an entropy inequality (instead of, for a MUSCL technique, to yield a maximum principle), the approximation of the unknowns at the face must be ”sufficiently close to” the upwind approximation. The entropy inequality is then obtained only in the weak sense. The technique to reach this result consists in splitting the remainder terms appearing in Lemma 3 and 4 in two parts: the first one is non-negative under some condition for the face approximation (hence the above mentioned limitation requirement); the second one is conservative and can be bounded in a discrete negative Sobolev norm (this explains why the entropy estimate is only a weak one).

Let $\varphi_{\rho}$ and $\varphi_{e}$ be the functions defined by (16) and let $\sigma\in{\mathcal{E}}_{{\rm int}}$, $\sigma=K|L$; by Lemma 2, there exists a unique $\rho_{KL}^{n+1}\in|[\rho_{K}^{n+1},\rho_{L}^{n+1}]|$ and $e_{KL}^{n+1}\in|[e_{K}^{n+1},e_{L}^{n+1}]|$ such that

[TABLE]

Entropy estimates are obtained in Theorem 3.2 under the following conditions:

[TABLE]

where $\rho_{KL}^{n+1}$ and $e_{KL}^{n+1}$ are defined by (33). Note that these conditions are satisfied by the upwind scheme (27). They may be seen as an additional constraint to be added to the limitation of a MUSCL-like procedure (see also the conclusion of the last section of this paper).

Theorem 3.2 (Entropy inequalities, implicit MUSCL-like scheme)

Let us assume that, for $\sigma\in{\mathcal{E}}_{{\rm int}}$, $\sigma=K|L$ and for $0\leq n\leq N-1$, the approximate density $\rho_{\sigma}^{n+1}$ and internal energy $e_{\sigma}^{n+1}$ in the numerical mass fluxes (28) and in the internal energy balance (27b) satisfy the conditions (34).

Then any solution of the scheme (27) satisfies, for any $K\in{\mathcal{M}}$ and $0\leq n\leq N-1$:

[TABLE]

where the remainder term $\delta\!R_{\eta}$ satisfies $\sum_{K\in{\mathcal{M}}}|K|\ (\delta\!R_{\eta})_{K}^{n+1}=0$ so that, integrating in space (i.e. summing over the cells), the following global discrete entropy estimate holds for $0\leq n\leq N-1$:

[TABLE]

In addition, let us suppose that there exists $M>0$ such that $\rho_{K}^{n}\leq M$, $1/\rho_{K}^{n}\leq M$, $e_{K}^{n}\leq M$, $1/e_{K}^{n}\leq M$ and $|u_{K,\sigma}^{n}|\leq M$ for $K\in{\mathcal{M}}$, $\sigma\in{\mathcal{E}}(K)$ and $0\leq n\leq N$, and let us define the quantities $|\varphi_{\rho}^{\prime}|_{\infty}=\max(|\varphi^{\prime}_{\rho}(1/M)|,\ |\varphi^{\prime}_{\rho}(M)|)$ and $|\varphi_{e}^{\prime}|_{\infty}=\max(|\varphi^{\prime}_{e}(1/M)|,\ |\varphi^{\prime}_{e}(M)|)$. Then the remainder term $\delta\!R_{m}$ satisfies the following bound:

[TABLE]

Therefore, a Lax-consistency property holds; more precisely, any limit $(\bar{\rho},\bar{\boldsymbol{u}},\bar{e})$ of a converging sequence of approximate solutions bounded in the $L^{\infty}$ and BV norms satisfies (22).

Proof

Let $(\delta\varphi_{\rho})^{n+1}_{\sigma}$ be defined by:

[TABLE]

By Lemma 3, (29) holds; an easy computation shows that the term associated to the face $\sigma$ in the expression of the remainder term $(R_{m})_{K}^{n+1}$ satisfies:

[TABLE]

with \displaystyle(F^{R}_{m})_{K,\sigma}^{n+1}=|\sigma|\,\frac{1}{2}\,\Bigl{[}\varphi^{\prime}_{\rho}(\rho_{K}^{n+1})-\varphi^{\prime}_{\rho}(\rho_{L}^{n+1})\Bigr{]}\ (\rho_{\sigma}^{n+1}-\rho_{KL}^{n+1})\ u_{K,\sigma}^{n+1}. Thanks to the assumption (34a), since $\varphi^{\prime}_{\rho}$ is an increasing function, $(F^{R}_{m})_{K,\sigma}^{n+1}\geq 0$. Let us define $(\delta\!R_{m})_{K}^{n+1}$, $K\in{\mathcal{M}}$, $0\leq n\leq N-1$ by:

[TABLE]

Then, under assumption (34a), we get:

[TABLE]

Let us prove that $\delta\!R_{m}$ satisfies:

[TABLE]

Indeed, since both $\rho_{\sigma}^{n+1}$ and $\rho_{KL}^{n+1}$ lie in the interval $|[\rho_{K}^{n+1},\ \rho_{L}^{n+1}]|$, we have by convexity of $\varphi_{\rho}$:

[TABLE]

Let $\psi$ be a function of ${\mathrm{C}}^{\infty}_{c}(\Omega\times(0,T))$. We have, thanks to the conservativity of the remainder term:

[TABLE]

Therefore,

[TABLE]

which concludes the proof of (39).

Following the same line of thought for the internal energy balance, let $(\delta\varphi_{e})_{\sigma}^{n+1}$ be defined by:

[TABLE]

and $(\delta\!R_{e})_{K}^{n+1}$ the remainder term given by:

[TABLE]

Thanks to the assumption (34b) we get:

[TABLE]

In addition, $\delta\!R_{e}$ satisfies the following inequality:

[TABLE]

Combining the inequalities (38) and (42) and thanks to (20), (39) and (43) concludes the proof of the theorem.

3.3 Explicit schemes

The general form of the discrete analogue of System (14) for an explicit scheme reads:

[TABLE]

where the numerical mass flux $F_{K,\sigma}^{n}$ is still defined by (28).

Let $\rho_{KL}^{n}$ (resp. $e_{KL}^{n}$) be the real number defined by Equation (26) with $x_{K}=\rho^{n}_{K}$ (resp. $x_{K}=e^{n}_{K}$) and $x_{L}=\rho^{n}_{L}$ (resp. $x_{L}=e^{n}_{L}$) and $\varphi=\varphi_{\rho}$ (resp. $\varphi=\varphi_{e}$) and let us assume that for $\sigma\in{\mathcal{E}}_{{\rm int}}$, $\sigma=K|L$ and for $0\leq n\leq N-1$,

[TABLE]

With these two conditions, Theorem 3.3 below yields a weak discrete entropy inequality, in the sense that a remainder term exists which tends to 0 (under some conditions) with the mesh and time steps, but its sign is unknown. However, in the case of an upwind approximation of the density and the internal energy on the faces of the mesh (note that (45) and (46) are satisfied for such approximations), a local discrete entropy inequality can be obtained under the following additional conditions.

First, the normal face velocities $u_{K,\sigma}$ in the mass flux (28) are assumed to be either

computed from a discrete velocity field ${\boldsymbol{u}}$:

[TABLE]

where ${\boldsymbol{n}}_{K,\sigma}$ is the unit normal vector to $\sigma$ outward $K$ and ${\boldsymbol{u}}_{\sigma}$ is an approximation of the velocity at the face, which may be the discrete unknown itself (when the velocity degrees of freedom are those of a non conforming Crouzeix-Raviart or Rannacher-Turek approximation, see e.g. [14]) or an interpolation (for instance, for a colocated arrangement of the unknowns, as in [25]).

• -

the unknown themselves in the case of the staggered MAC scheme, since only the normal velocity is approximated in this case, see e.g. [14].

For $1\leq r$, we then define the following discrete norm:

[TABLE]

where ${\mathcal{E}}^{(i)}(K)={\mathcal{E}}(K)$ for the Crouzeix-Raviart or Rannacher-Turek case and ${\mathcal{E}}^{(i)}(K)$ is restricted to the two faces of $K$ perpendicular to the $i^{th}$ vector of the canonical basis of ${\mathbb{R}}^{d}$ in the case of the MAC scheme.

Remark 2 (Discrete $L^{r}(W^{1,r})$ norm of the velocity)

It is reasonable to suppose that, under regularity assumptions of the mesh whose precise statement depends the space approximation at hand, this norm is equivalent to the standard finite-volume discrete $L^{r}(0,T;W^{1,r})$ norm [10]; it is indeed true for usual cells (in particular, with a bounded number of faces) for staggered discretizations and for a convex interpolation of the velocity at the faces for colocated schemes. 2. 2.

Second, the following CFL conditions hold:

[TABLE]

where $\tilde{\rho}_{K}^{n+1/2}\in|[\rho_{K}^{n},\rho_{K}^{n+1}]|$, $\rho_{K,\sigma}^{n}\in|[\rho_{K}^{n},\ \rho_{L}^{n}]|$, $\tilde{e}_{K}^{n+1/2}\in|[e_{K}^{n},e_{K}^{n+1}]|$ and $e_{K,\sigma}^{n}\in|[e_{K}^{n},\ e_{L}^{n}]|$ are defined by:

[TABLE]

Theorem 3.3 (Discrete entropy inequalities, explicit schemes)

Let $\rho$ and $e$ satisfy the relations of the scheme (44). Let $M\geq 1$ and let us suppose that $\rho_{K}^{n}\leq M$, $1/\rho_{K}^{n}\leq M$, $e_{K}^{n}\leq M$, $1/e_{K}^{n}\leq M$ and $|u_{K,\sigma}|\leq M$, for $K\in{\mathcal{M}}$, $\sigma\in{\mathcal{E}}(K)$ and $0\leq n\leq N$. Assume that the discretization of the convection term in (44a) and (44b) satisfies the assumptions (45) and (46) respectively. Let $\eta$ be defined by (17). Then any solution of the scheme (44) satisfies, for any $K\in{\mathcal{M}}$ and $0\leq n\leq N-1$:

[TABLE]

where $R_{\eta}=R_{\eta,1}+R_{\eta,2}$ with:

[TABLE]

where $\underline{h}_{\mathcal{M}}$ is defined by (23), $|\varphi^{\prime}_{\rho}|_{\infty}=\max(|\varphi^{\prime}_{\rho}(1/M)|,\ |\varphi^{\prime}_{\rho}(M)|)$, $|\varphi^{\prime}_{e}|_{\infty}=\max(|\varphi^{\prime}_{e}(1/M)|,\ |\varphi^{\prime}_{e}(M)|)$, and $|\varphi^{\prime\prime}_{\rho}|_{\infty}$ and $|\varphi^{\prime\prime}_{e}|_{\infty}$ denote the maximum value taken by $\varphi^{\prime\prime}_{\rho}$ and $\varphi^{\prime\prime}_{e}$ respectively on the interval $[1/M,\ M]$.

Moreover, if the discretization of the convection term in (44a) and (44b) is upwind, and if the normal face velocities $u_{K,\sigma}$ satisfy (47), under the CFL conditions (49) and (50), we also have (with a different expression for $R_{\eta}$):

[TABLE]

where $q\geq 1$, $q^{\prime}\geq 1$ and $\dfrac{1}{q}+\dfrac{1}{q^{\prime}}=1$, $f_{\mathcal{M}}$ and $C_{\mathcal{M}}$ are defined by (24).

Proof

The results are obtained by applying the propositions 1 and 2 below with $\varphi=\varphi_{\rho}$ and $\varphi=\varphi_{e}$ respectively.

The aim of the following proposition is to derive a discrete analogue of Relation (18).

Proposition 1 (Discrete renormalized forms of the mass balance equation)

Let $\varphi$ be a twice continuously differentiable convex function from $(0,+\infty)$ to ${\mathbb{R}}$, and let $\rho$ satisfy (44a). Let $M\geq 1$ and let us suppose that $\rho_{K}^{n}\leq M$, $1/\rho_{K}^{n}\leq M$ and $|u_{K,\sigma}|\leq M$, for $K\in{\mathcal{M}}$, $\sigma\in{\mathcal{E}}(K)$ and $0\leq n\leq N$. Let $|\varphi^{\prime}|_{\infty}=\max(|\varphi^{\prime}(1/M)|,\ |\varphi^{\prime}(M)|)$ and $|\varphi^{\prime\prime}|_{\infty}$ be the maximum value taken by $\varphi^{\prime\prime}$ on the interval $[1/M,\ M]$. Assume that $\rho_{\sigma}^{n}$ satisfies (45). Then the following inequality holds:

[TABLE]

where the remainder $(R_{\rho})_{K}^{n+1}=(R_{\rho,1})_{K}^{n+1}+(R_{\rho,2})_{K}^{n+1}$ with:

[TABLE]

where $\underline{h}_{\mathcal{M}}$ is defined by (23).

Assume furthermore that the normal face velocities $u_{K,\sigma}$ satisfy (47), that the discretization of the convection term in (44a) is upwind and that the CFL condition (49) holds with $\varphi$ instead of $\varphi_{\rho}$. Then (56) still holds (with a different expression for $R_{\rho}$) and:

[TABLE]

where $q\geq 1$, $q^{\prime}\geq 1$, $\dfrac{1}{q}+\dfrac{1}{q^{\prime}}=1$, ${\lVert\cdot\rVert}_{L^{q^{\prime}}(0,T;W^{1,q^{\prime}}_{\mathcal{M}})}$ is defined by (48), $f_{\mathcal{M}}$ and $C_{\mathcal{M}}$ are defined by (24) and $C$ only depends on the maximal number of faces of the mesh cells.

Proof

Mimicking the formal computation performed at the continuous level, let us multiply (44a) by $\varphi^{\prime}(\rho_{K}^{n+1})$. We get:

[TABLE]

with

[TABLE]

By a Taylor expansion, there exists $\rho_{K}^{n+1/2}\in|[\rho_{K}^{n},\ \rho_{K}^{n+1}]|$ such that:

[TABLE]

The term $(T_{2})_{K}^{n+1}$ reads:

[TABLE]

with

[TABLE]

Thanks to assumption (45), the remainder $R_{2}$ is a sum of a non-negative part and a term tending to zero; indeed there exists $\delta\!R_{2}$ such that:

[TABLE]

This result is obtained by adapting the proof of the implicit case (indeed, up to a change of time exponents at the right-hand side from $n$ to $n+1$, the expression of $(R_{2})_{K}^{n+1}$ is the same than the second term of $(R_{m})_{K}^{n+1}$ in the expression (29), and the computation from Relation (36) up to the end of the proof of (39) may be reproduced, still with the same change of time exponents).

Let us now prove that the remainder term $R=(R_{K}^{n})_{K\in{\mathcal{M}}}^{n=0,\ldots,M}$ defined by (57) satisfies:

[TABLE]

Indeed, for $K\in{\mathcal{M}}$ and $0\leq n\leq N$, we get:

[TABLE]

where $\tilde{\rho}_{K}^{n+1/2}$ is defined by (51). Thus,

[TABLE]

which yields (60).

Let us now turn to the case where the discrete normal velocities satisfy (47) and the discretization of the density at the face $\rho_{\sigma}^{n}$ is upwind; in this case the remainder $R_{2}$ defined by (3.3) satisfies:

[TABLE]

where $\rho_{K,\sigma}^{n}$ is defined by (52). Therefore, $R_{2}$ is non-negative. Starting from Equation (61), we may now reformulate the remainder term $R_{K}^{n+1}$ as $R_{K}^{n+1}=(R_{01})_{K}^{n+1}+(R_{02})_{K}^{n+1}$ with:

[TABLE]

By Young’s inequality, the second term may be estimated as follows:

[TABLE]

Therefore, in view of the expressions (58) and (63) of $(R_{1})_{K}^{n}$ and $(R_{2})_{K}^{n}$ respectively, we get $(R_{1})_{K}^{n+1}+(R_{2})_{K}^{n+1}+(R_{02})_{K}^{n+1}\geq 0$ under the CFL condition (49) (with $\varphi$ instead of $\varphi_{\rho}$). Let us now show that

[TABLE]

with $q\geq 1$, $q^{\prime}\geq 1$ and $\dfrac{1}{q}+\dfrac{1}{q^{\prime}}=1$, where $f_{\mathcal{M}}$ and $C_{\mathcal{M}}$ is defined by (24). To this purpose, we first observe that, in the Crouzeix-Raviart or Ranncher-Turek case, since $\sum_{\sigma\in{\mathcal{E}}(K)}|\sigma|\ {\boldsymbol{n}}_{K,\sigma}=0$, we may write:

[TABLE]

where ${\boldsymbol{u}}_{K}^{n}$ stands for the mean value of the normal face velocities $({\boldsymbol{u}}_{\sigma}^{n})_{\sigma\in{\mathcal{E}}(K)}$. In the MAC case, we have $\sum_{\sigma\in{\mathcal{E}}^{(i)}(K)}|\sigma|\ u^{n}_{K,\sigma}=0$ for $1\leq i\leq d$, and thus

[TABLE]

where ${\boldsymbol{j}}^{(i)}$ stands for the $i^{th}$ vector of the canonical basis of ${\mathbb{R}}^{d}$ and $u^{n}_{K,i}$ stands for the mean value of the velocity over the two faces of ${\mathcal{E}}^{(i)}(K)$. In both cases, we obtain that:

[TABLE]

Therefore,

[TABLE]

We thus have, thanks to a Hölder estimate, for $q\geq 1$, $q^{\prime}\geq 1$ and $\dfrac{1}{q}+\dfrac{1}{q^{\prime}}=1$:

[TABLE]

Using $|\rho_{K}^{n+1}-\rho_{K}^{n}|^{q}\leq(2M)^{q-1}\ |\rho_{K}^{n+1}-\rho_{K}^{n}|$ yields (65).

The object of the following proposition is to mimick at the fully discrete level the computation of (2) and (3), applying it to $z=e$.

Proposition 2 (Inequalities derived from the internal energy balance)

Let $M\geq 1$ and let us suppose that $\rho_{K}^{n}\leq M$, $e_{K}^{n}<M$, $1/e_{K}^{n}\leq M$ and $|u_{K,\sigma}|\leq M$, for $K\in{\mathcal{M}}$, $\sigma\in{\mathcal{E}}(K)$ and $0\leq n\leq N$. Assume that the face approximation of the internal energy satisfies (45). Let $\varphi$ be a twice continuously convex function from $(0,+\infty)$ to ${\mathbb{R}}$. Then

[TABLE]

with

[TABLE]

where $|\varphi^{\prime}|_{\infty}=\max(|\varphi^{\prime}(1/M)|,\ |\varphi^{\prime}(M)|)$, $|\varphi^{\prime\prime}|_{\infty}$ stands for the maximum value taken by $\varphi^{\prime\prime}$ over the interval $[1/M,\ M]$, and $\underline{h}_{\mathcal{M}}$ is defined by (23). If, furthermore, the approximation of $e_{\sigma}^{n}$ in (66) is upwind, and the CFL condition (50) holds (with $\varphi$ instead of $\varphi_{e}$) then $(R_{e})_{K}^{n}=0$.

Proof

First, the fully discrete identity corresponding to the semi-discrete identity (5) with $z=e$ is obtained thanks to the discrete mass equation; it reads:

[TABLE]

Now let $\varphi$ be a twice continuously differentiable function from$(0,+\infty)$ to ${\mathbb{R}}$, and let us multiply the first two terms of the discrete internal energy balance (44b) by $\varphi^{\prime}(e_{K}^{n+1})$; switching from the conservative to the non conservative form, we get:

[TABLE]

with

[TABLE]

The remainder term $R_{K}^{n+1}$ is quite similar to the remainder defined by (57) in the proof of Proposition 1; following the proof of (60), we get that it satisfies:

[TABLE]

Now

[TABLE]

with:

[TABLE]

The remainder $(R_{1})_{K}^{n+1}$ may be written:

[TABLE]

where $e_{K}^{n+1/2}\in|[e_{K}^{n},\ e_{K}^{n+1}]|$. Since $\varphi$ is supposed to be convex, this term is non-negative. Let $e_{KL}^{n}$ be the real number defined by Equation (26) (and denoted in this latter relation by $x_{KL}$) with $x_{K}=e^{n}_{K}$ and $x_{L}=e^{n}_{L}$. Thanks to (45), by a computation similar to the implicit case, the remainder $R_{2}$ satisfies

[TABLE]

Switching back from the non-conservative formulation to the conservative formulation yields:

[TABLE]

which, thanks to (70), leads to (66) and (67).

Let us now suppose that the discretization of the internal energy convection term is upwind. In this case, we obtain for $(R_{2})_{K}^{n+1}$:

[TABLE]

where $e_{K,\sigma}^{n}\in|[e_{K}^{n},\ e_{L}^{n}]|$. The remainder $R_{K}^{n+1}$ yields in the upwind case:

[TABLE]

where $e_{K}^{n+1/2}\in|[e_{K}^{n},\ e_{K}^{n+1}]|$. So, thanks to the Young inequality:

[TABLE]

In view of the expressions (69) and (72) of $(R_{1})_{K}^{n+1}$ and $(R_{2})_{K}^{n+1}$ respectively, we obtain that $(R_{1})_{K}^{n+1}+(R_{2})_{K}^{n+1}+R_{K}^{n+1}\geq 0$ thanks to the CFL condition (50), which yields the result.

Theorem 3.3 deserves the following comments:

First, in the explicit case, we are able to prove neither a local nor a global discrete entropy inequality; we only obtain some weak inequalities that allow to show the consistency of the scheme, under some conditions.

• -

The convergence to zero with the space and time step of the remainders is obtained, supposing a control of discrete solutions in $L^{\infty}$ and discrete BV norms, in two cases: first when the ratio $\delta t/\underline{h}_{\mathcal{M}}$ tends to zero, second when the $L^{q}(0,T;W^{1,q}_{\mathcal{M}})$ norm of the velocity does not blow-up too quickly with the space step. To this respect, let us suppose that we implement a stabilization term in the momentum balance equation reading (in a pseudo-continuous setting, for short and to avoid the technicalities associated to the space discretization), for $1\leq i\leq d$:

[TABLE]

where $\Delta_{q}u_{i}$ is such that

[TABLE]

where $C$ is independent of $h_{\mathcal{M}}$. This kind of viscosity term may be found in turbulence models [1, 36]. Multiplying (73) by $u_{i}$ and integrating with respect to space and time yields:

[TABLE]

In this relation, the right-hand side may be controlled under $L^{\infty}$ and BV stability assumptions (remember that, at the discrete level, the BV and $W^{1,1}$ norms are the same), and we obtain an estimate on ${\lVert{\boldsymbol{u}}\rVert}_{L^{q}(0,T;W^{1,q}_{\mathcal{M}})}$ which may be used in (55). A standard first order diffusion-like stabilizing term corresponds to $q=2$ and $\alpha=1$; it yields a bound on $h_{\mathcal{M}}^{1/2}{\lVert{\boldsymbol{u}}\rVert}_{L^{2}(0,T;H^{1}_{\mathcal{M}})}$, so that (55) becomes

[TABLE]

Such a stabilization is thus not sufficient to ensure that the remainder term tends to zero. What is needed is in fact:

[TABLE]

To avoid an over-diffusion in the momentum balance, this inequality suggests to implement a non-linear stabilization with $q>2$ which, in turn, will allow $\alpha>1$. With such a trick, we will be able to obtain for first-order upwind schemes the desired ”Lax-consistency” result: the limit of a convergent sequence of solutions, bounded in $L^{\infty}$ and BV norms, and obtained with space and time steps tending to zero, satisfies a weak entropy inequality.

• -

We introduced in [33] a limitation process for a MUSCL-like algorithm for the transport equation, which consists in deriving an admissible interval for the approximation of the unknowns at the mesh faces, in convection terms, thanks to extrema preservation arguments. This limitation process has been extended to the Euler equations in [14]. The conditions (45) and (46) may easily be incorporated in this limitation: indeed, they also define an admissible interval, which is not disjoint from the MUSCL-like admissible interval of [33], since the upwind value belongs to both. A similar idea (namely restricting the choice for the face approximation in order to obtain an entropy inequality) may be found in [2].

——————————————————————————————————–