DeepAI

# Robust second-order approximation of the compressible Euler equations with an arbitrary equation of state

This paper is concerned with the approximation of the compressible Euler equations supplemented with an arbitrary or tabulated equation of state. The proposed approximation technique is robust, formally second-order accurate in space, invariant-domain preserving, and works for every equation of state, tabulated or analytic, provided the pressure is nonnegative. An entropy surrogate functional that grows across shocks is proposed. The numerical method is verified with novel analytical solutions and then validated with several computational benchmarks seen in the literature.

• 1 publication
• 5 publications
• 8 publications
• 4 publications
• 1 publication
11/03/2022

### On the invariant region for compressible Euler equations with a general equation of state

The state space for solutions of the compressible Euler equations with a...
07/19/2019

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

We present an implicit-explicit finite volume scheme for the Euler equat...
09/13/2020

### Second-order invariant domain preserving approximation of the compressible Navier–Stokes equations

We present a fully discrete approximation technique for the compressible...
08/27/2020

### On Fake Accuracy Verification

In this paper, we reveal a mechanism behind a fake accuracy verification...
03/12/2018

### A second-order well-balanced finite volume scheme for Euler equations with gravity

We present a well-balanced, second order, Godunov-type finite volume sch...
01/20/2022

### Analytic Adjoint Solutions for the 2D Incompressible Euler Equations Using the Green's Function Approach

The Green's function approach of Giles and Pierce is used to build the l...
06/30/2020

### Efficient parallel 3D computation of the compressible Euler equations with an invariant-domain preserving second-order finite-element scheme

We discuss the efficient implementation of a high-performance second-ord...

## 1 Introduction

The objective of this paper is to continue the work started in Clayton et al. [6] where we introduced an explicit and invariant-domain preserving method to approximate the Euler equations equipped with an equation of state that is either tabulated or is given by an expression that is so involved that elementary Riemann problems cannot be solved exactly or efficiently. The method proposed in [6]

is invariant-domain preserving in the sense that it guarantees that the density and the internal energy of the approximation are positive. A key ingredient of the method is a local approximation of the equation of state using a covolume ansatz (Noble-Abel equation of state) from which upper bounds on the maximum wave speed are derived for every elementary Riemann problem. The downside is that the method is first-order accurate in space. In the present paper we propose a technique that increases the accuracy in space to second-order and preserves the invariant-domain properties. We also introduce a functional that can be used as an entropy surrogate. The functional in question is shown to be increasing across shocks in the Riemann problem involved in the construction of the local maximum wave speed. We show how to use this functional to limit the internal energy from below. This feature is useful for equations of state that are tabulated or interpolated from experimental data since in this case no natural notion of entropy is available.

It is sometimes possible, though possibly expensive, to solve Riemann problems when the equation of state is analytic. For instance this is done in Colella and Glaz [7, §1], Ivings et al. [23], Quartapelle et al. [32]. This cannot be done with tabulated equations of state because the information on the pressure is incomplete. Several attempts to develop methods working with an arbitrary or tabulated equation of state have been reported in the literature. One way to do so consists of using approximate Riemann solvers like those found in Dukowicz [8], [7, §2], Roe and Pike [34], Pike [31], and Lee et al. [25]. One can also simplify the Riemann problem by using flux splitting techniques as those reported in Toro et al. [39]. We also refer to Saurel et al. [35], Banks [1], Dumbser and Casulli [9], Dumbser et al. [10] where approximation techniques are developed using approximate Riemann solvers for various equations of state. Some of these techniques guarantee positivity of the density, but little else is guaranteed in general. The method introduced in Clayton et al. [6] is based on a graph viscosity technique using upper bounds on the maximum wave speed in the Riemann problem. Instead of using the two-shock approximation of the Riemann solution, as done in most methods based on approximate solvers, the method proposed in [6]

approximates the pressure in the Riemann fan by the covolume equation of state and subsequently estimates guaranteed upper bounds on the maximum wave speed. This in turn ensures that the internal energy in the proposed algorithm is positive in addition to the density being positive. If it happens that the equation of state is of covolume type, then the method also preserves the minimum principle on the specific entropy. A technique based on similar principles is reported in

Wang and Li [40] where the authors use a stiffened gas equation of state to approximate the pressure in the Riemann fan. To the best of our knowledge, the method proposed in the present paper is among the very first ones that are provably invariant-domain preserving for complex or tabulated equations of state and high-order accurate in space.

The paper is organized as follows. In Section 2, we introduce the mathematical model of interest and the corresponding notation. We also briefly discuss the assumptions we make on the equation of state and recall results from [6] that are needed for this work (we recall the first-order approximation of the Euler equations in §2.3). In Section 3

, we construct a provisional update that is higher-order accurate in space. This update is based on a high-order graph viscosity using an entropy commutator and an activation function. Then, in Section

4, we apply a novel convex limiting technique that corrects the invariant-domain violations of the provisional higher-order method. The final result is an approximation technique that is robust, formally second-order accurate in space, provably invariant-domain preserving, and works for every equation of state (tabulated or analytic) that satisfies the mild assumptions stated in §2.2. Finally in Section 5, the method is verified with analytical solutions and published benchmarks and is validated with experiments. A short conclusion is given in §6.

## 2 Preliminaries

We formulate the problem and introduce notation in this section. We also recall essential results from [6] for completeness.

### 2.1 The Euler equations

Let be a bounded polyhedron in . Given some initial data and initial time , we look for solving the compressible Euler equations in the weak sense: equationparentequation

 (1a) ∂tρ+\DIV(\bvρ)=0 \aet>t0, \bx∈\Dom, (1b) ∂t\bbm+\DIV(\bv⊗\bbm+p(\bu)\polId)=\bzero \aet>t0, \bx∈\Dom, (1c) ∂tE+\DIV(\bv(E+p(\bu)))=0 \aet>t0, \bx∈\Dom.

The components of the dependent variable

(considered to be a column vector) are the density,

, the momentum, , and the total mechanical energy, . We also introduce the velocity and the specific internal energy, . To simplify the notation later on we introduce the flux , where is the identity matrix.

### 2.2 Equation of state

In (1) the function denotes the pressure which we assume to be given by an oracle. This means we assume to have no a priori knowledge of the function itself apart from some mild structural assumptions stated below. We only assume that for a given state we are able to retrieve a pressure in a suitable way, for example by evaluating an arbitrary analytic expression, or by deriving a value from tabulated experimental data. More precisely, in the numerical illustrations reported in §5, we consider oracles given by analytic functions for the ideal gas, van der Waals, Jones–Wilkins–Lee and Mie-Gruneisen equations of state. We also use the SESAME database developed at Los Alamos National Laboratory [27] to test the method with experimental tabulated data.

Throughout the paper we assume that the domain of definition of the oracle is the set given by

 (2) \calB(b)\eqq{\bu\eqq(ρ,\bbm,E)∈Rd+2\st0<ρ, 0<1−bρ,0

We henceforth refer to as the admissible set. One of the objectives of the paper is to guarantee that the approximation is high-order accurate and leaves invariant. The inequality found in the definition of is the so called maximum compressibility condition. The constant can be set to zero if the user has no a priori knowledge about the maximum compressibility of the fluid under consideration. We recall, however, that a large class of analytic equations of state, such as the Noble-Abel, the Mie-Gruneisen (with the Hugoniot locus as the reference curve), and the Noble-Abel-Stiffened-Gas equations of state all involve a maximum compressibility constant. Finally, we assume that the oracle returns a non-negative pressure,

 (3) p:\calB(b)→\Real≥0.

The assumption can be weakened, but for the sake of simplicity, we refrain from doing so in this paper. We leave this extension for future works.

### 2.3 First-order time and space approximation

The time and space discretization proposed in Clayton et al. [6] is based on [16]. This method is in some sense a discretization-agnostic generalization of an algorithm introduced by Lax [24, p. 163]. We denote by the current time, , and we denote by the current time step size; that is . Without going into the details of the space approximation, we assume that the current approximation of is a collection of states , where the index set

is used to enumerate all the degrees of freedom of the approximation, and

is in for all . The update at is obtained as follows:

 (4) mi\dt(\bsfUi\upLnp−\bsfUni)+∑j∈\calI(i)\polf(\bsfUnj)\bcij−∑j∈\calI(i)∖{i}dij\upLn(\bsfUnj−\bsfUni)=\bzero.

The quantity is called the lumped mass and we assume that for all . The vector encodes the space discretization. The index set is called the local stencil. This set collects only the degrees of freedom in that interact with (). We assume that the -valued coefficients are such that approximates at some grid point , and the consistency error in space scales optimally with respect to the mesh size for the considered approximation setting. We require that the method be conservative; more precisely, we assume that

 (5) \bcij=−\bcjiand∑j∈\calI(i)\bcij=\bzero.

Concrete expressions for and are given in [19, §4] for continuous and discontinuous finite elements as well as for finite volumes. The computations reported at the end of this paper are done with piecewise linear continuous finite elements. But to stay general, we continue with the abstract discretization-agnostic notation introduced above.

For completeness we now recall how the graph viscosity is defined in [6]. Given and , we set . For , we set , , , and . The non-negativity assumption on the pressure (3) and the assumptions on the density () implies that

 (6) γZ\eqqΓZϱZ≥1.

Then we consider the following Riemann problem:

 (7) ∂t⎛⎜ ⎜ ⎜⎝ρmEΓ⎞⎟ ⎟ ⎟⎠+∂x⎛⎜ ⎜ ⎜ ⎜ ⎜⎝m1ρm2+\calpcovmρ(E+\calpcov)mρΓ⎞⎟ ⎟ ⎟ ⎟ ⎟⎠=0,with\calpcov(ρ,m,\calE,Γ)\eqqΓρ−11−bρ(E−m22ρ),

with left data and right data . This problem is well-posed because . Its complete solution is given in Clayton et al. [6, §4]. We denote by the maximum wave speed in this problem. Let be a nontrivial convex subset of . We say that is an invariant set for (7) if for every pair of Riemann data in the solution of (7) takes values in . We then have:

###### Theorem 1 ([6, Thm. 4.6])

Let . Let be a convex invariant set for (7). Assume that for all . For all , let be any positive number larger than or equal to . Let

 (8) dij\upLn\eqqmax(\wlambda(\bnij,\bsfUni,\bsfUnj)∥\bcij∥ℓ2,\wlambda(\bnji,\bsfUnj,\bsfUni)∥\bcji∥ℓ2).

Assume that is small enough so that . Let be the update defined in (4). Then .

###### Remark 1 (Invariant domains)

Theorem 1 asserts that every convex invariant set in is invariant by the update procedure (4) with the artificial viscosity defined in (8). This means that if for all , then for all . We say that the method is invariant-domain preserving for . Notice in particular that the method is invariant-domain preserving for the admissible set . The admissible set may not be the smallest invariant domain. For instance, if the oracle admits a mathematical entropy , then the approximation defined above is also invariant-domain preserving for the set ; see [16, Cor. 4.2]. By slight abuse of terminology and provided the context is unambiguous we will simply call a method invariant-domain preserving without quantifying the precise convex set.

A source code to compute is available online [5]. A notable drawback of the graph viscosity (8) is that it reduces the space accuracy of the method to first order. The remainder of the paper is concerned with constructing a higher-order approximation and a respective novel convex limiting technique that is invariant-domain preserving.

###### Remark 2 (Pressure approximation)

Notice that the oracle is only invoked to compute the left and right values of the pressure in (7). The pressure in the Riemann fan is approximated by two covolume equations of state. There is one covolume equation of state for each side of the contact discontinuity since on the left of the contact wave and on the right.

## 3 Provisional higher-order method

In this section, we introduce a provisional higher-order method that extends the update (4) to second order when using linear finite elements for the discretization. A convex limiting technique for this provisional update is introduced in §4. The provisional method is based on the work reported in [18, 28].

### 3.1 The method

Higher-order accuracy in space requires using the consistent mass matrix instead of the lumped mass matrix for the discretization of the time derivative. By reducing dispersive errors, the consistent mass matrix is known to yield superconvergence at the grid points; see Christon et al. [4], [14], [28, Sec. 3.4], Thompson [38]. Let be the mass matrix with entries , then can be approximated by where , is the Kronecker symbol, and . This approximation has been shown in [4, 14] to be superconvergent for piecewise-linear continuous finite elements. Let with , then using this approximation we have . Notice that , because

. The skew-symmetry of the summand in

is used in §4 for limiting purposes.

Let denote the high-order update (here the superindex reminds us that the update is higher order accurate). Then, for every , the provisional high-order update is given by: equationparentequation

 (9a) miτ(\bsfUH,n+1i−\bsfUni)=\bRni+∑j∈\calI(i)(bij\bRnj−bji\bRni), (9b) with\bRni\eqq∑j∈\calI(i)(−\polf(\bsfUnj)\SCAL\bcij+dij\upHn(\bsfUnj−\bsfUni)).

The high-order graph viscosity coefficient , defined in Section 3.2, shares the same properties as its low-order counterpart which are necessary for conservation:

 (10) dij\upHn=dji\upHn,0≤dij\upHn≤dij\upLn,for% every j∈\calI(i),i∈\calV.

### 3.2 Entropy commutator

Using the entropy-viscosity methodology introduced in Guermond et al. [17], we now discuss the construction of the high-order graph viscosity coefficient that is used in the high-order update (9b). A complicating factor in this construction is that the pressure is given by an oracle. We thus do not have a priori knowledge on the equation of state and entropies of the system.

Recall that is an entropy pair for the Euler equations if

where is the flux of the system (1). Since one may not have access to entropies for tabulated equations of state, we are going to use at every and every time one entropy pair associated with the following flux

 (12) \polfi,n(\bu)\eqq⎛⎜ ⎜⎝\bbm\bv⊗\bbm+\calpcovi,n(\bu)\polId\bv(E+\calpcovi,n(\bu))⎞⎟ ⎟⎠,\calpcovi,n(\bu)\eqq(γmin,ni−1)ρe(\bu)1−bρ.

Here, with (see (6)). We use the following shifted Harten entropy pair in the numerical tests reported in §5:

 (13) ηi,n(\bu)\eqq⎛⎝ρ2e(\bu)(1−bρ)1−γmin,ni⎞⎠1γmin,ni+1−ρϱniηi,nref,\bFi,n(\bu)\eqq\bbmρηi,n(\bu),

where . Then, we estimate “entropy production” by inserting the approximate solution into a discrete counterpart of (11) which we write as follows:

where are the shape functions of the finite element approximation. Notice that the above identity holds true for every smooth function and for every . Substituting into (14), we estimate the local entropy residual as follows:

The residual (15) can be thought of as a measure of how well the discrete solution satisfies (14) in each local stencil . We then define the normalized entropy residual as follows: equationparentequation

 (16a) Rni \eqq\absNniDi+ϵDmax,Dmax=maxi∈\calVDi,ϵ=10−2, (16b) Dni \eqq∣∣∑j∈\calI(i)\bFi,n(\bsfUnj)\SCAL\bcij∣∣+∣∣∑j∈\calI(i)(\GRAD\buηi,n(\bsfUni))\tr\polfi,n(\bsfUnj)\SCAL\bcij∣∣.

By definition, the normalized residual has values in and, for piecewise linear finite elements, behaves like where is the typical meshsize. Finally we set

 (17) dij\upHn\eqqdij\upLnmax(ψ(Rni),ψ(Rnj)),

where the activation function is defined as follows:

 (18) ψ(x)\eqq4x30−(x+x0)(x−2x0){(x−2x0)−ReLU(x−2x0)}4x30,

and is the rectified linear activation function. Notice that , and for all for a chosen . As a result, one recovers if the entropy residual is larger than . Note that for (the identity also holds for ). When we have . The numerical tests reported in §5 are done with (). An activation function with the same purpose is used in Persson and Peraire [30, Eq. (8)]. Up to a translation, the activation function therein behaves like in the interval .

###### Remark 3 (Entropy shift)

The entropy shift considered in (13) is motivated by the observation that the numerator is unchanged by the change of variable for all and for every entropy associated with the flux (12). The constant is chosen in (13) so that . This entropy shift was first is introduced in [18, §3.4].

## 4 Convex Limiting

The high-order update (9) is not guaranteed to be invariant-domain preserving; in particular, it is not guaranteed to stay in the admissible set . In order to correct this defect we now discuss a new convex limiting technique that re-establishes invariant-domain preservation for the final high-order update, .

### 4.1 Key observation

A key observation is that one can rewrite (4) as follows: equationparentequation

 (19a) \bsfUi\upLnp=(1−∑j∈\calI(i)∖{i}2\dtdij\upLnmi)\bsfUni+∑j∈\calI(i)∖{i}2\dtdij\upLnmi¯¯¯¯¯¯¯¯¯¯¯¯¯¯\bsfUij\upn, (19b) with\oline\bsfUnij\eqq12(\bsfUni+\bsfUnj)−12dij\upLn(\polf(\bsfUnj)−\polf(\bsfUni))\SCAL\bcij.

That is, the low-order update is a convex combination (under the appropriate CFL condition) of the local state and the auxiliary states ,

 (20) \bsfUi\upLnp∈\conv{¯¯¯¯¯¯¯¯¯¯¯¯¯¯\bsfUnij\stj∈\calI(i)}.

The main result established in [6, Thm. 4.6] (and summarized in Theorem 1) is that under the CFL condition stated in Theorem 1 and the definition (8) for , the states are in provided this is already the case of the states . This is done by proving that the states are space averages of the solution to a Riemann problem. We refer the reader to [16, Sec. 3.3], [18, Sec. 3.2] and [19, Sec. 3.2] where this is discussed in detail. We are going to use the states to define local bounds in space and time to perform the limiting of the high-order states .

### 4.2 Entropy surrogate

Our goal is to use the methodology introduced in [17] to perform the convex limiting of the update . In this context the use of an oracle with little a priori knowledge on the equation of state poses a significant challenge as it makes it impossible to properly define an entropy. We resolve the impasse by introducing an artificial surrogate entropy that has the right mathematical properties for the convex limiting methodology to be applied; see Theorem 2.

For any admissible state and , we define

 (21) S(\bu,γ):=(ρe)(\bu)ργ(1−bρ)γ−1,

where and are the density and specific internal energy of the state . Furthermore, for every index , we set

 (22) γni\eqq1+\sfpni(1−bϱni)ϱni\sfeni,Γni\eqqϱniγni,

where , and . The following result is the key motivation for the definition of an entropy surrogate.

###### Lemma 1

For all , assume that . For all , all , all , all left data and right data in the extended Riemann problem (7), and all , we set

 (23) Ψij(\bu):=ρe(\bu)−Sminijργij(1−bρ)1−γij,whereSminij\eqqmin(S(\bsfUni,γij),S(\bsfUnj,γij)).

Then, increases across shocks in the solution of the extended Riemann problem (7) (if a shock wave exists).

###### Proof

We omit the superscript in the proof to simplify the notation. The solution to the extended Riemann problem (7), , is given in Clayton et al. [6]. In particular, we have that if and if . Here, , and is the speed of the contact wave. Let , with the convention that the index is for the left state and is for the right state. Assume that the -wave is a shock wave. Let be the state before the shock and let be an arbitrary state connected to through a shock curve. With the notation for the specific volume (not to be confused with the time step), and since along the wave curve, the Rankine-Hugoniot condition implies that (see Godlewski and Raviart [12, Eq. (4.8), p. 144]), , where by slight abuse of notation we renamed the pressure in (7) by setting . We then infer that only depends on along the shock curve; more precisely, we have

 e(\bu)=\sfeZ1−(γZ−1)(τ−τZ)2(τZ−b)1+(γZ−1)(τ−τZ)2(τ−b)\qqer(τ).

Notice that this function is well defined only on the interval with

 (24) τ∞Z\eqq(γZ−1)τZ+2bγZ+1,

and it is nonegative on the interval with . Notice that since we assumed that . We now show that the function is nonnegative and increasing on the shock curve for all and . This will then prove the assertion for the choices and , because and

 0

. Setting , we have ; hence, the pressure, , is a monotone increasing function of along shock curves. By definition of shock curves, starting from the state the pressure increases, we conclude that also increases along shock curves, or . Actually, the pressure is finite only in the range ; hence, showing that is nonnegative increasing on shock curves, is equivalent to showing that

 (τ∞Z,τZ]∋τ↦g(τ)\eqqτ−1r(τ)−cτ−1(τ−b)1−γ

is nonnegative decreasing on shock curves. Recall that the specific entropy for the covolume equation of state is given by . A fundamental property of the specific entropy is that it is an increasing function along shocks. That is, is a decreasing function over the interval . This also means that is a decreasing function. We now have

 g(τ)=τ−1r(τ)−cτ−1(τ−b)1−γ=ς(τ)τ−1(τ−b)1−γZ−cτ−1(τ−b)1−γ.

Let be the defined by . Then we have

 g(τ)=(τ−b)1−γτ(˜ς(τ)−c).

Computing the derivative of , we see,

 g′(τ)=(τ−b)1−γτ˜ς′(τ)−(γ−1)τ+(τ−b)τ2(τ−b)γ(˜ς(τ)−c).

Note that

 ˜ς′(τ)=(τ−b)γ−γZς′(τ)+(γ−γZ)(τ−b)γ−γZ−1ς(τ).

Then as for , , and , and as is a decreasing function. Thus, by the choice of the constants and , we have that and . Hence, is nonnegative increasing on shock curves. This completes the proof.

We now define the surrogate entropy. For all , we set equationparentequation

 (25a) γmin,ni \eqqminj∈\calI(i)γnj,Smin,ni\eqqmin(minj∈\calI(i)S(\bsfUnj;γmin,ni),minj∈\calI(i)S(¯¯¯¯¯¯¯¯¯¯¯¯¯¯\bsfUnij;γ