Entropy Conserving and Kinetic Energy Preserving Numerical Methods for the Euler Equations Using Summation-by-Parts Operators

Summation-by-parts operators can be used in the context of finite difference and discontinuous Galerkin methods to create discretisations mimicking properties given at the continuous level such as entropy conservation. Recently, there has been some interest in schemes for the Euler equations that additionally preserve the kinetic energy. However, some these methods resulted in undesired and unexpected changes of the kinetic energy in numerical experiments of Gassner et al. (J Comput Phys 327:39–66, 2016). Here, analytical insights into kinetic energy preservation are given and new entropy conservative and kinetic energy preserving numerical fluxes are proposed.


Introduction
Considering the solution of hyperbolic conservation laws, high order methods can be very efficient, providing accurate numerical solutions with relatively low computational effort [21]. In order to make use of this accuracy, stability has to be established. Mimicking estimates obtained on the continuous level via integrationby-parts, summation-by-parts (SBP) operators [22,37] can be used. In short, SBP operators are discrete derivative operators equipped with a compatible quadrature providing a discrete analogue of the L 2 norm. The compatibility of discrete integration and differentiation mimics integration-by-parts on a discrete level. Combined with the weak enforcement of boundary conditions via simultaneous approximation terms (SATs) [1], highly efficient and stable semidiscretisations can be obtained at least for linear problems, see e.g. [6,14,39] and references cited therein.
In recent years, there has been an enduring and increasing interest in the basic ideas of SBP operators and their application in various frameworks including finite volume (FV) [25,26], discontinuous Galerkin (DG) [2,4,10,11,13,20,27,28,30], and the recent flux reconstruction/correction procedure via reconstruction framework [15,16,42] as described in [31,32]. While there is only a limited amount of well-posedness theory for nonlinear conservation laws, mimicking properties such as entropy stability semidiscretely has received much interest. Building on the seminal work of Tadmor [40,41], entropy stability of second order schemes using symmetric numerical fluxes has been investigated, resulting in well-defined properties that numerical fluxes have to satisfy in order to result in entropy conservative H. Ranocha ( ) TU Braunschweig, Institute Computational Mathematics, Braunschweig, Germany e-mail: h.ranocha@tu-bs.de schemes. Decomposing general semidiscretisations into a non-dissipative central part and an additional dissipative part, suitable artificial dissipation or filtering can be added afterwards, cf. [7,9,38]. Second order methods based on symmetric numerical fluxes can be extended to high order in a conservative way, cf. [4,7,28] and [8,23,[34][35][36].
Another property of numerical methods for the Euler equations that has received much interest in the literature concerns the kinetic energy. A structural property of numerical fluxes described by Jameson [18] has been used to construct socalled kinetic energy preserving (KEP) numerical fluxes inter alia by Chandrashekar [3]. However, schemes using these fluxes do not preserve the kinetic energy as expected in numerical experiments by Gassner et al. [12]. They had to change the discretisation of the pressure to reduce undesired changes of the kinetic energy. However, this resulted in a loss of entropy conservation. Motivated by these results, some analytical insights into this behaviour have been developed in [29,Section 7.4] and will be presented here.
This chapter is structured as follows. At first, some basic results about SBP operators and corresponding semidiscretisations of hyperbolic conservation laws are reviewed in Sect. 2. Afterwards, the Euler equations are considered in Sect. 3. After demonstrating that the property that has been used to characterise numerical fluxes as KEP is not well-defined, the new concept of KEP numerical methods is introduced. Moreover, a numerical flux that is both entropy conservative and kinetic energy preserving in the new sense is developed. Thereafter, results of a numerical experiment comparing entropy conservative numerical fluxes are described in Sect. 4. Finally, a brief summary is given in Sect. 5.

Discretisations Using Summation-by-Parts Operators
Consider the Euler equations in two space dimensions where ρ is the density, v the velocity, e the specific total energy, and p the pressure. For a perfect gas, p = (γ − 1) ρe − 1 2 ρv 2 . The usual entropy is U = − ρs γ −1 , where s = log p − γ log ρ is the specific (physical) entropy.
With the entropy fluxes F j fulfilling ∂ u U · ∂ u f j = ∂ u F j , smooth solutions of the Euler equations in d space dimensions satisfy ∂ t U(u) + d j =1 ∂ j F j (u) = 0 and the entropy inequality is used as additional admissibility criterion for weak solutions, cf. [5].
In order to discretise (1), the domain Ω is divided into several non-overlapping sub-domains Ω l ⊆ Ω and SBP operators will be used on each element. SBP operators consist of discrete derivative operators D j , approximating the partial derivative in direction j , and a symmetric and positive definite mass/norm matrix M, Moreover, an interpolation operator R approximates the restriction of functions on Ω l to the boundary ∂Ω l and a symmetric and positive definite boundary mass matrix B approximate the L 2 (∂Ω l ) scalar product. Representing the multiplication by the j -th component of the outer unit normal ν at ∂Ω l by the diagonal matrix n j , the SBP property has to be satisfied in order to mimic integration-by-parts discretely via Semidiscretisation of (1) will be constructed as follows. Each sub-domain Ω l ⊆ Ω is mapped onto a reference element and all computations are performed there. On each element, the resulting semidiscretisation is of the form where the volume terms VOL discretise the flux divergence in the interior of Ω l and the surface terms SURF couple elements or impose boundary conditions. Here, u is the vector of the nodal values of the numerical solution at specified nodes ξ i in Ω l and a collocation approach is used. Thus, nonlinear operations are performed pointwise and the discrete fluxes f j are given by their nodal values f . As in (nodal) discontinuous Galerkin methods, the surface terms will be built using numerical fluxes f num,j in the j -th coordinate direction as Finally, the volume terms are constructed using symmetric (two-point) numerical fluxes f vol,j (volume fluxes) that are consistent with f j as where VOL i is the volume term at ξ i [7]. If f vol,j are smooth fluxes, the discretisation (7) is of the same order of accuracy as the derivative matrices D j [4,28]. Moreover, if the mass matrix M is diagonal, this approximation can be written in a conservative form [7]. Finally, if the boundary operators R T Bn j R are also diagonal and f vol,j are entropy conservative in the sense of Tadmor [40, 41], the semidiscretisation (5) is entropy conservative/stable across elements if the numerical surface fluxes f num,j are entropy conservative/stable. Moreover, some results on the kinetic energy can be transferred as well [12]. In the following, the focus will lie on the fluxes f vol,j .

Euler Equations and Kinetic Energy
The kinetic energy E kin = 1 2 ρv 2 fulfils (for sufficiently smooth solutions) Jameson [18] investigated the kinetic energy in a one-dimensional semidiscrete setting using finite volume methods. To simplify the notation, this setup will be used in the following; its extension to multiple dimensions is straightforward. Jameson proposed to mimic (8) semidiscretely by using numerical momentum fluxes of the is the numerical density flux, and p num is a consistent numerical approximation of the pressure. Later, this has been used as a kind of "definition" of kinetic energy preserving (KEP) numerical fluxes, e.g. in [3,12]. However, this is not a well-defined concept, cf. [28,29]. Indeed, every numerical momentum flux can be written as Since the numerical fluxes are consistent, p num := f num ρv − {{v}}f num ρ is a consistent approximation of the pressure. The insufficiency of the condition f num ρv = {{v}}f num ρ + p num is in accordance with observations of Gassner et al. [12]. They investigated a Taylor-Green vortex problem and compared several numerical fluxes for the Euler equations. There, numerical fluxes of the form f num ρv = {{v}}f num ρ + p num with p num = {{p}} resulted in a clear loss of kinetic energy compared to other KEP fluxes using the arithmetic average p num = {{p}} as approximation of the pressure. They observed that "the discretisation of the pressure plays a crucial role for the kinetic energy" and that the choice of the arithmetic average p num = {{p}} "seems to be important for the kinetic energy equation" [12, Section 4.2]. However, they had no (theoretical) explanations for this observation.

New Approach to Kinetic Energy Preservation
By a heuristic argument, the balance law (8) may not be suitable in the incompressible limit: Indeed, for smooth solutions, (8) can be rewritten as which becomes a conservation law for smooth solutions of the incompressible Euler equations due to div(v) = 0 or an energy inequality similar to the entropy inequality (2). Since the kinetic energy is plays a crucial role in the incompressible limit [24], the second form (10) might be considered the "better" one. Thus, a semidiscretisation mimicking this equation might be desirable near the incompressible limit.

Theorem 1 (Corollary 7.5 of [29]) If a kinetic energy preserving numerical flux is used in a semidiscrete FV method, the resulting semidiscrete kinetic energy equation mimics both the conservative and the non-conservative terms of Eq. (10).
Proof (Sketch) Using the chain rule in a one dimensional finite volume setting, the time derivative of the kinetic energy in cell i becomes Using the momentum flux f num ρv = {{v}}f num ρ + {{p}} in the volume terms (7) in one dimension, the arithmetic average of the pressure yields the volume term Dp, i.e. a straightforward discretisation of ∂ x p. Analogous results hold in multiple space dimensions, cf. Sect. 2.

Entropy Conservative and KEP Numerical Fluxes
Since entropy stability has received much interest and the entropy conservative numerical fluxes of [3,17] are not KEP in the sense of Definition 1, it is interesting whether both concepts can be fulfilled simultaneously. The logarithmic mean value {{ρ}} log = [[ρ]]/[[log ρ]] has been proposed by Roe [33] in the context of entropy conservative numerical fluxes and is described in [17]. Many useful entropy conservative numerical density fluxes are of the form f num ρ = {{ρ}} log {{v}}, e.g. the one presented in [3]. This form seems to be preferable, since positivity preservation of the density can be achieved using local Lax-Friedrichs/Rusanov dissipation operators [28, Section 6.2]. Using this ansatz for f num ρ and Definition 1, the following entropy conservative and kinetic energy preserving numerical flux (f num,y analogously) has been constructed in [29,Section 7.4]

Numerical Results
Since the kinetic energy is an important quantity for the incompressible Euler equations, a Taylor-Green vortex given by v y (t, x, y) = − cos(x) sin(y), p(t, x, y) = 100 γ + cos(2x) + cos(2y) 4 , for (x, y) ∈ [0, 2π ] 2 with periodic boundary conditions is considered, which is a stationary solution of the incompressible Euler equations. Using tensor product Lobatto bases for polynomials of degree p = 5 on N = 16 elements per coordinate direction, the numerical solutions have been computed in the time interval t ∈ [0, 30] with the fourth order, ten-stage, strong stability preserving Runge-Kutta method of [19]. The time step Δt has been chosen as Δt = cfl min Δx/(2p + 1)λ , where λ is the greatest absolute value of the eigenvalues of f and the minimum is taken over all cells and nodes. As in [12], the given numerical fluxes have been used for both the volume terms (7) and as surface fluxes in (6), without additional dissipation. The evolution of the entropy U and the kinetic energy E kin using a CFL number cfl = 0.9 for the entropy conservative fluxes of Ismail and Roe [17], Chandrashekar [3], and the new flux (11) are visualised in Fig. 1. As can be seen there, the entropy remains approximately constant and the kinetic energy oscillates uniformly until t ≈ 20. Afterwards, the kinetic energy drops for the fluxes of [3,17] and there is a relative change of the entropy of order 10 −5 . Contrary, there is no visible change for the new flux (11).
The entropy loss for the fluxes of Ismail and Roe [17] and Chandrashekar [3] is caused by the time integration scheme, as can be seen in Fig. 2, where the time step is reduced by an order of magnitude (cfl = 0.09). However, the behaviour of the kinetic energy is nearly unchanged.

Summary and Discussion
Using summation-by-parts operators, high order numerical schemes with specific properties can be constructed using symmetric (two-point) numerical fluxes. While several "kinetic energy preserving" methods have been proposed, they have been characterised by a property of the numerical fluxes that is not well-defined. Such numerical fluxes resulted in schemes that did not preserve the kinetic energy as expected [12]. Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence and indicate if changes were made. The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.