A new energy stable formulation of the compressible Euler equations

by   Jan Nordström, et al.
Linköping University

We show that a specific skew-symmetric form of hyperbolic problems leads to energy conservation and an energy bound. Next, the compressible Euler equations is transformed to this skew-symmetric form and it is explained how to obtain an energy estimate. Finally we show that the new formulation lead to energy stable and energy conserving discrete approximations if the scheme is formulated on summation-by-parts form.



There are no comments yet.


page 1

page 2

page 3

page 4


Nonlinear and Linearised Primal and Dual Initial Boundary Value Problems: When are they Bounded? How are they Connected?

Linearisation is often used as a first step in the analysis of nonlinear...

Entropy stable numerical approximations for the isothermal and polytropic Euler equations

In this work we analyze the entropic properties of the Euler equations w...

Analysis of a Backward Euler-type Scheme for Maxwell's Equations in a Havriliak-Negami Dispersive Medium

For the Maxwell's equations in a Havriliak-Negami (H-N) dispersive mediu...

A Class of A Stable Summation by Parts Time Integration Schemes

Since integration by parts is an important tool when deriving energy or ...

Entropy Stable Numerical Fluxes for Compressible Euler Equations which are Suitable for All Mach Numbers

We propose two novel two-state approximate Riemann solvers for the compr...

Energy of Computing on Multicore CPUs: Predictive Models and Energy Conservation Law

Energy is now a first-class design constraint along with performance in ...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

The energy method applied to linear IBVPs lead to well posed boundary conditions and energy estimates Kreiss (1970); Kreiss and Lorenz (1989); Gustafsson and Sundstrom (1978); Gustafsson et al. (1995); Oliger and Sundström (1978); Nordström and Hagstrom (2020); Nordström (2017); Nordström and Svärd (2005). The most common procedure to obtain estimates for nonlinear problems is to use the entropy stability theory Godunov (1961); Vol’pert (1967); Kružkov (1970); Dafermos (1973); Lax (1973); Harten (1983); Tadmor (1984, 1987); Tadmor,Eitan (2003). The main character of these two approaches are (roughly speaking): The energy analysis uses integration-by-parts (IBP) as the main tool and leads to an estimate in an equivalent norm. Normally, symmetric matrices and a symmetrising matrix is required which can be hard to find in the nonlinear case (although exeptions exist Nordström and Winters (2021)

). The entropy analysis needs both IBP and the chain rule and aims for conservation of the mathematical entropy, a convex functions related to the physical entropy. In combination with certain sign requirements (for compressible flow the density and temperature must be positive) it leads to

estimates, otherwise not. It can with relative ease be applied to nonlinear equations without specific symmetry requirements on the involved matrices.

In this note we will apply the general stability theory developed in Nordström (2022). This theory is valid for both linear and nonlinear problems and extend the use of the energy method to nonlinear problems. It is very direct, easy to understand and leads to estimates. The only requirement for the energy bound is that a certain skew-symmetric form of the governing equations exist. It was shown in Nordström (2022) that this form exists for the velocity-divergence form of the incompressible Euler equations and could be derived for the shallow water equations. The drawback with this procedure is that the required skew-symmetric form can be quite complicated to derive. In this paper we will explain this procedure in detail for the compressible Euler equations. Once the skew-symmetric formulation is obtained, an energy bound follows by applying IBP, without using the chain rule and without sign requirements. By discretising the equations using summation-by-parts (SBP) operators Svärd and Nordström (2014); Fernández et al. (2014), nonlinear stability follows by discretely mimicking the IBP procedure.

The remaining part of paper is organised as follows: In Section 2 we shortly reiterate the main theoretical findings in Nordström (2022) and outline the general procedure. We proceed in Section 3 to choose an appropriate solution norm on which we base our choice of new dependent variables. With a suitable form of the norm, we derive new governing equations in the new variables. Next we rewrite the governing equations in skew-symmetric form and show how to get an energy bound. Section 4 illustrate the relation between the new continuous formulation and stability of the numerical scheme. A summary and conclusions are provided in Section 5.

2 The main theory

The general theory extending the energy method to the nonlinear case in Nordström (2022) is shortly summarised here. Consider the following general hyperbolic IBVP


augmented with homogeneous boundary conditions at the boundary . In (2.1), the Einstein summation convention is used and is a symmetric positive definite (or semi-definite) time-independent matrix that defines an energy norm (or semi-norm) . We assume that and are smooth. The matrices are smooth functions (each matrix element is smooth) of the

component vector

, but otherwise arbitrary. Note that (2.1) encapsulates both linear () and nonlinear () problems. The following two concepts are essential for a proper treatment of (2.1).

Definition 2.1.

The problem (2.1) is energy conserving if only changes due to boundary effects. It is energy bounded if as .

Proposition 2.1.

The IBVP (2.1) for linear () and nonlinear () is energy conserving if


holds. It is energy bounded if it is energy conserving and the boundary conditions are such that


The energy method applied to (2.1) yields


where is the outward pointing unit normal. The terms on the right-hand side of (2.4) are cancelled by (2.2) leading to energy conservation. If in addition (2.3) holds, an energy bound is obtained. ∎

Remark 2.2.

Proposition 2.1 shows that whatever form the original IBVP has, energy boundedness and energy conservation can be proved if it can be rewritten in the form given by (2.1)-(2.2). The procedure to arrive at an energy estimate involve the following steps.

  1. Find an appropriate energy norm from which one can choose new dependent variables.

  2. Derive a new set of governing equations in the new variables from the standard governing equations.

  3. Transform the new set of equations into a skew-symmetric formulation as in (2.1).

  4. Apply the energy method such that only boundary terms remain as in (2.4).

  5. Find a minimal number of boundary conditions that limits the resulting boundary terms as in (2.3).

3 The new skew-symmetric form of the compressible Euler equations

We start by considering the choice of norm and choose new dependent variables.

3.1 Task 1: Find an appropriate norm and new dependent variables

The total energy in a two-dimensional compressible ideal gas is a combination of internal energy , kinetic energy and potential energy . We have , where is the pressure, is the ratio of specific heats, is the density, are the velocities in the direction respectively and is the gravity times the height. Note that has the dimension . Based on this observation we choose the dependent variables and the diagonal norm matrix (inspired by E) to be:


where are positive but yet unknown. The total energy connects to the quadratic form on which we base the norm by observing that all terms involved have dimension if have dimension and are non-dimensional. This completes the first task in Remark 2.2.

3.2 Task 2: Rewriting the compressible Euler equations in new dependent variables

The compressible Euler equations in primitive form using the pressure and the ideal gas law is


Repeated use of the chain rule and the introduction of the new variables in (3.1) leads to the new equations


For convenience we have used the relations . A matrix-vector form of (3.3) is




This completes the second task in Remark 2.2.

3.3 Task 3: Transforming the new set equations into a skew-symmetric form

Proposition 2.1 implies that we can proceed each direction separately. By taking into account that we aim for an energy estimate using the (yet unknown) matrix in (3.1), we need to solve


In (3.6

t( we have wo systems of ordinary differential equations for the 32 entries

in the matrices and the 3 unknowns on the diagonal in . The number of equations is 8 and we have 4 variables which leads to 32 relations. A simple counting argument implies that the problem is likely solvable. Another important observation is that both the matrices and the norm are coupled and solved for together. To exemplify the procedure we consider the -direction and later directly provide the results for the -direction.

We start with row 1 in (3.6) for the -direction. The following equation holds


There are no entries involving on the righthand side (RHS) of the equation. Hence we make the ansatz . This ansatz cancels all terms related to if and with as arbitrary constants. The remaining terms on the RHS are


with as a free parameter. Equating the terms in (3.7) using (3.8) yields leaves no terms on the RHS in (3.8). Hence we get the same type of solution as for , i.e. . This provide the matrix with a determined first row and column as


Next we consider row 2 in (3.6). The following equation holds


where and are already determined. Since there are no entries involving on the RHS of the equation, we find (as for row 1) that and is a solution. By inspecting the terms related to we see that and is the only possible solution. The remaining terms on the RHS of (3.10) are rewritten as


with as a free parameter. By inserting the known values of and as well as making the ansatz based on the RHS in (3.11) with as free parameters we find the relation


with the solution . We now have determined also the second row and column


The procedure is now clear. One proceeds row by row with a decreasing number of new entries to determine.

In the third step one makes the ansatz (for the same reason as in step 2) and finds that and must hold. The resulting matrix after the third step is


In the fourth step most of the matrix is already determined, and hence we directly state the final equation which after the ansatz and realising that must hold becomes


The solution is given by . Note that this gives the first information about the norm matrix via the requirement for . The final matrix and related norm becomes


To be precise, the third element in is obtained in the derivation of the matrix given below


This completes the third task in Remark 2.2.

3.4 Task 4: Applying the energy method such that only boundary terms remain

We multiply (3.4) with from the left, use (3.6) and integrate over the domain . By using Greens formula and Proposition 2.1 we find


where is the outward pointing unit normal from the boundary . The relation (3.18) shows that energy is conserved in the sense that it only changes due to boundary effects.

The matrices and the norm matrix contain 5 undetermined parameters and . Except for the parameter , which is part of the norm, the energy rate cannot depend on these parameters since they are not present in (3.4), (3.5) and (3.6). Hence as a sanity check we compute the boundary contraction involving only the terms multiplied by and find


showing that the free parameters do not influence the energy rate. The remaining boundary contraction is


where we introduced the normal velocity . This completes the fourth task in Remark 2.2.

3.5 Task 5: The choice of nonlinear boundary conditions

A nonlinear and linear analysis may lead to a different number and type of boundary conditions required for an energy bound. This was discussed extensively in Nordström and Winters (2021), Nordström (2022) where the boundary matrix was found to be different in the linear and nonlinear case, and also to have a different meaning. For completeness we will repeat part of that discussion here. For more details we refer the reader to Nordström and Winters (2021), Nordström (2022).

We start by rotating the velocities to be normal () or aligned () with the boundary. By inserting these transformation into (3.20) we obtain the rotated boundary contraction


where the rotated solution is

. By considering the eigenvalues of the matrix we see that they indicate a similar but not identical sign pattern as in the linear case. We find the eigenvalues


where , is the speed of sound and , the normal Mach number.

The relations (3.22) indicate that for outflow () we have two situations. When there are 4 positive eigenvalues and no boundary condition is required. For , 3 eigenvalues are positive, 1 is negative and 1 boundary condition seem to be required. For inflow () we have the reversed situation with 4 negative eigenvalues and four required boundary conditions for , which goes down to 3 negative eigenvalues and 3 required boundary conditions for .

Remark 3.1.

The shift from subsonic to supersonic flow at is generally assumed to be crucial and to modify the number of boundary conditions in a linear analysis. Interestingly, here in the nonlinear analysis the shift occur when if , Maybe even more interesting is that for .

However, considering eigenvalues is not sufficient for nonlinear problems Nordström and Winters (2021),Nordström (2022). Expanding (3.21) give


which proves that no boundary conditions are necessary in the outflow case. It also proves that specifying the normal velocity to zero is correct at a solid wall. This completes the fifth and final task in Remark 2.2.

3.6 The final form of the governing equations

The derivations above focused on stability and resulted in matrices that were functions of the constant norm matrix as seen in (3.6). The final governing equations can be transformed to




This removes the dependence of the norm in the matrices which take the form (without free parameters)


3.7 Some open questions

It is interesting to consider the energy rate. By combining and expanding (3.18) and (3.23) we find


This means the rate of change in the volume of the energy () is increasing or decreasing due to the transport of energy in or out of the domain plus an additional amount due to pressure work .

As we have seen the does not contribute to the rate of change in the energy but could be part of the scheme by populating the matrices in (3.26) and (3.27). It is an open question whether they will modify the spectrum and hence time-integration procedure. The parameter in the norm could be any positive number that defines a reasonable norm. It has no influence on the scheme.

There are two possible interpretations of the sign requirements for the density and the pressure (and hence temperature). The first interpretation considers the problem in a physical way which means that and must be positive, otherwise the new variables involving square roots do not exist. In the second opposite interpretation, the new variables are considered as the ones defining the original variables. With this point of view, the sign problem vanishes since by squaring and both and will always be positive.

4 A stable and energy conserving numerical approximation

To exemplify the straightforward construction of stable and energy conserving schemes based on the new formulation, we consider a summation-by-parts (SBP) approximation of (3.24),(3.25) as given in


where include approximations of in each node. The matrix elements of are matrices with node values of the matrix elements in injected on the diagonal as exemplified below


Moreover and where are 1D SBP difference operators, are positive definite diagonal quadrature matrices, satisfies the SBP constraint , denotes the Kronecker product and with subscripts denote identity matrices. All matrices have appropriate sizes such that the matrix-matrix and matrix-vector operations are defined. Based on the 1D SBP operators, the 2D SBP relations mimicking integration by parts are given by


where and contain numerical integration along rectangular domain boundaries. In (4.3) we have used , and .

The discrete energy method (multiply (4.1) from the left with ) where yields


where we have used that commutes with . Next, the discrete relations corresponding to (3.25),


and the notation transforms (4.4) to


Finally, the SBP relations (4.3), and (the matrices consist of diagonal blocks) yield


The semi-discrete energy rate in (4.7) mimics the continuous result in (3.18) and hence the scheme is energy conserving. Stability can be obtained by adding a proper dissipative boundary treatment.

Remark 4.1.

It is irrelevant whether the problem is linear or nonlinear. The skew-symmetric formulation, an SBP discretisation and a proper boundary treatment lead to stability and energy conservation.

5 Summary and conclusions

We have shown that a specific skew-symmetric form of linear and nonlinear problem leads to an energy bound and energy conservation for the compressible Euler equations. The skew-symmetric formulation automatically produced energy conserving and energy stable numerical schemes for the compressible Euler equations if these are formulated on summation-by-parts form.

The derivation shoved that the skew-symmetric formulation required a coupled derivation of the matrices and the norm matrix. The derivations also indicated that the shift in number of boundary conditions might not occur precisely at Mach number = 1, but at a slightly lower number given by for . It also shoved that this ratio is identically one for .

This paper together with Nordström (2022) have shown that the incompressible Euler equations, the shallow water equations and the compressible Euler equations can all be transformed to skew-symmetric form. Once in that form stable, energy conserving, easy to apply schemes follows if summation–by parts operators are used for the discretisation. The estimates are obtained by using integration-by-parts or summation-by-parts. No additional requirements, such as chain rules or sign requirements are needed.


Many thanks to my colleagues Fredrik Laurén and Andrew R. Winters for helpful comments on the manuscript. Jan Nordström was supported by Vetenskapsrådet, Sweden [award no. 2018-05084 VR] and the Swedish e-Science Research Center (SeRC).


  • C. M. Dafermos (1973) The entropy rate admissibility criterion for solutions of hyperbolic conservation laws. J. Differ. Equations 14 (2), pp. 202–212. Cited by: §1.
  • D. C. D. R. Fernández, J. E. Hicken, and D. W. Zingg (2014)

    Review of summation-by-parts operators with simultaneous approximation terms for the numerical solution of partial differential equations

    Computers & Fluids 95, pp. 171–196. Cited by: §1.
  • S. K. Godunov (1961) An interesting class of quasilinear systems. In Dokl. Acad. Nauk SSSR, Vol. 11, pp. 521–523. Cited by: §1.
  • B. Gustafsson and A. Sundstrom (1978) Incompletely parabolic problems in fluid dynamics. SIAM J. Appl. Math. 35 (2), pp. 343–357. Cited by: §1.
  • B. Gustafsson, H. Kreiss, and J. Oliger (1995) Time dependent problems and difference methods. Vol. 24, JWS. Cited by: §1.
  • A. Harten (1983)