A skew-symmetric energy and entropy stable formulation of the compressible Euler equations

We show that a speciﬁc skew-symmetric form of nonlinear hyperbolic problems leads to energy and entropy bounds. Next, we exemplify by considering the compressible Euler equations in primitive variables, transform them to skew-symmetric form and show how to obtain energy and entropy estimates. Finally we show that the skew-symmetric formulation lead to energy and entropy stable discrete approximations if the scheme is formulated on summation-by-parts form


Introduction
The energy method applied to linear initial boundary value problems (IBVPs) lead to well posed boundary conditions and energy estimates [1,2,3,4,5,6,7,8]. The energy analysis uses integration-by-parts (IBP) as the main tool and leads to an estimate in an L 2 equivalent norm. Symmetric matrices and a symmetrising matrix are normally required which can be hard to find in the nonlinear case (although exeptions exist [9]). The most common procedure to obtain estimates for nonlinear problems is to use the entropy stability theory [10,11,12,13,14,15,16,17,18]. The entropy analysis needs both IBP and the chain rule and aims for conservation of mathematical entropy (a convex function more or less related to the physical entropy). In combination with certain sign requirements (for compressible flow the density and temperature must be positive) it leads to L 2 estimates, otherwise not. It can with relative ease be applied to many nonlinear equations without specific symmetry requirements on the involved matrices.
In this paper we will apply and extend the general stability theory developed in [19]. 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 L 2 estimates. The only requirement for the energy bound is that a certain skew-symmetric form of the governing equations exist. It was shown in [19] that this form exists for the velocity-divergence form of the incompressible Euler equations and could be derived for the shallow water equations (SWEs). One 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 and use the compressible Euler equations as an example. Once the skew-symmetric formulation is obtained, an energy bound follows by applying IBP, without using the chain rule and without sign requirements. Although we focus on energy estimates, we show that the new formulation also allows for a mathematical (or generalised) entropy conservation and bound. By discretising the equations in space using summation-by-parts (SBP) operators [20,21], nonlinear stability follows by discretely mimicking the IBP procedure.
Skew-symmetric formulations for parts or the whole set of governing flow equations have drawn interest previously. In [22], a similar set of variables, were considered. L 2 continuous estimates were obtained for a subset of the variables. The discrete aspect was discussed but not analysed. Aiming for simulation of turbulence, also [23] considered similar variables. Focus was on various conservation properties which were deduced from an a' priori assumption that an energy bound existed. The discrete conservation and stability aspects were discussed, and supporting calculations were provided, but no proofs were given. A related ambition using conventional primitive variables were provided in [24]. Focus was again on conservation properties, in particular on preservation of kinetic energy in aeroacoustic calculations and supporting calculations were provided. Also in [25] primitive variables were used and skew-symmetry was targeted in order to preserved moments in plasma physics calculations. Neither in [24] nor [25] were proofs provided.
The papers [22,23,24,25] (and references therein) contain fragments of the general theory for nonlinear hyperbolic problems presented in [19] and in this paper. This theory include the compressible Euler equations which we use in this paper to exemplify the whole chain of actions leading to nonlinear stability. As stated above, The only requirement for the continuous and discrete bounds is that a certain skew-symmetric form of the governing equations exist. We show in detail how to arrive at this formulation. The remaining part of paper is organised as follows: In Section 2 we shortly reiterate the main theoretical findings in [19] and outline the general procedure for obtaining energy and entropy bounds. After that, we leave the general theory and we move to the compressible Euler equations which we use as the prime example. 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 and entropy bound. Section 4 illustrate the relation between the new continuous skew-symmetric formulation and stability of the numerical SBP based semi-discrete scheme. A summary and conclusions are provided in Section 5.

The main theory
The general theory extending the energy method to the nonlinear case in [19] is shortly summarised here. Consider the following general hyperbolic IBVP augmented with homogeneous boundary conditions L p U = 0 at the boundary ∂Ω. In (2.1), the Einstein summation convention is used and P is a symmetric positive definite (or semi-definite) time-independent matrix that defines an energy norm (or semi-norm) U 2 P = Ω U T P U dΩ. We assume that U and V are smooth. The n × n matrices A i , B i , C are smooth functions (each matrix element is smooth) of the n component vector V , but otherwise arbitrary. Note that (2.1) encapsulates both linear (V = U ) and nonlinear (V = U ) problems.

Energy analysis
The following two concepts are essential for a proper treatment of (2.1). Definition 2.1. The problem (2.1) is energy conserving if U 2 P = Ω U T P U dΩ only changes due to boundary effects. It is energy bounded if U 2 P < ∞ as t → ∞.
holds. It is energy bounded if it is energy conserving and the boundary conditions L p U = 0 are such that

3)
Proof. The energy method applied to (2.1) yields where (n 1 , .., n k ) T is the outward pointing unit normal. The terms on the right-hand side of (2.4) are cancelled by ( Remark 2.3. For linear problems, a minimal number of boundary conditions that lead to a bound is a necessary and sufficient condition for well-posedness. For nonlinear problems this is not the case. A minimal number of boundary conditions that lead to a bound is a necessary but not a sufficient condition [6,8,9].
We will in Section 3 show a detailed derivation of steps 1-3 for the compressible Euler equations. Steps 4-5 will be shortly reviewed. Before that, we will show that also a specific mathematical entropy is conserved.

Entropy analysis
For non-smooth nonlinear solutions U , the formulation (2.1) interpreted in a weak sense also allows for an entropy conservation law. Proposition 2.4. The IBVP (2.1) together with the conditions (2.2) leads to the entropy conservation law where S = U T P U/2 is the mathematical (or generalised) entropy and F i = U T A i U are the entropy fluxes.
Proof. Multiplication of (2.1) from the left with U T yields The right-hand side of (2.7) is cancelled by (2.2) leading to the entropy conservation relation (2.5).
Remark 2.5. The entropy conservation law (2.5) holds for smooth solutions. For discontinuous solutions it holds in a distributional sense. The non-standard compatibility conditions in this case reads The entropy S is convex (S UU = P ) and identical to the energy. A similar identity between the energy and entropy was found in [9] for the SWEs. In the following we will use energy to denote both quantities, but sometimes remind the reader by writing out both notations explicitly.

The new skew-symmetric form of the compressible Euler equations
We will go through the list in Remark 2.2 and start with the choice of norm and new dependent variables.

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 p/(γ −1), kinetic energy ρ(u 2 + v 2 )/2 and potential energy ρgh. We have E = p/(γ − 1) + ρ(u 2 + v 2 )/2 + ρgh, where p is the pressure, γ is the ratio of specific heats, ρ is the density, u, v are the velocities in the (x, y) direction respectively and gh is the gravity times the height. Note that gh has the dimension velocity 2 . Based on this observation we choose the dependent variables and the diagonal norm matrix (inspired by E) to be: where α 2 , β 2 , θ 2 are positive and non-dimensional but yet unknown. We have also introduced the arbitrary and constant velocity w. The total energy E connects to the quadratic form Φ T P Φ on which we base the norm Φ 2 P by observing that all terms involved have dimension density × velocity 2 . In the following, we set w=1 without restriction. This completes the first task in Remark 2.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 This completes the second task in Remark 2.2.

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 P in (3.1), we need to solve (3.6) In (3.6) we have two systems of ordinary differential equations for the 32 entries a ij , b ij , in the matricesÃ,B and the 3 unknowns on the diagonal in P . 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Ã,B and the norm P are coupled and solved for together. To exemplify the procedure we consider the x-direction and later directly provide the results for the y-direction.
We start with row 1 in (3.6) for the x-direction. The following equation holds There are no entries involving φ 3 , φ 4 on the righthand side (RHS) of the equation. Hence we make the ansatz a 13 =ã 13 φ 3 , a 31 =ã 31 φ 3 , a 14 =ã 14 φ 4 , a 41 =ã 41 φ 4 . This ansatz cancels all terms related to φ 3 , φ 4 if a 31 = −2ã 13 andã 41 = −2ã 14 withã 13 ,ã 14 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 η = −1, a 11 = α 2 u leaving no terms on the RHS in (3.8). Hence we get the same type of solution as for φ 3 , φ 4 , i.e.ã 21 = −2ã 12 . This provide the matrixÃ with a determined first row and column as (3.9) Next we consider row 2 in (3.6). The following equation holds where a 12 and a 21 are already determined. Since there are no entries involving φ 3 on the RHS of the equation, we find (as for row 1) that a 23 =ã 23 φ 3 and a 32 = −2ã 32 φ 3 is a solution. By inspecting the terms related to φ 4 we see that a 24 = 0 and a 42 = 4β 2 φ4 φ1 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 a 12 and a 21 as well as making the ansatz a 22 = (ã 22 φ 1 + ψ 2 u) based on the RHS in (3.11) withã 22 , ψ 2 as free parameters we find the relation with the solution η = 1/3, ψ 2 = β 2 ,ã 22 =ã 12 . We now have determined also the second row and columñ 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 a 33 = (ã 33 φ 1 + ψ 3 u) (for the same reason as in step 2) and finds thatã 13 =ã 23 =ã 33 = 0 and ψ 3 = θ 2 must hold. The resulting matrix after the third step is (3.14) In the fourth step most of the matrix is already determined, and hence we directly state the final equation which after the ansatz a 44 = (ã 44 φ 1 + ψ 4 u) and realising thatã 34 = 0 must hold becomes The solution is given by β 2 = (γ − 1)/2, ψ 4 = γ − 2,ã 44 =ã 14 . Note that this gives the first information about the norm matrix P via the requirement for β 2 . The final matrixÃ and related norm P becomes To be precise, the third element in P is obtained in the derivation of the matrixB given below This completes the third task in Remark 2.2.

Task 4: Applying the energy method such that only boundary terms remain
We multiply (3.4) with 2Φ T P from the left, use (3.6) and integrate over the domain Ω. By using Greens formula and Proposition 2.1 we find where (n x , n y ) T is the outward pointing unit normal from the boundary ∂Ω. The relation (3.18) shows that energy (and entropy) is conserved in the sense that it only changes due to boundary effects. The matricesÃ,B and the norm matrix P contain 5 undetermined parametersã 12 ,ã 14 ,b 13 ,b 14 and α 2 . Except for the parameter α 2 , 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ã 12 ,ã 14 ,b 13 ,b 14 and find showing that the free parameters do not influence the energy rate. The remaining boundary contraction is where we introduced the normal velocity u n = n x u + n y v. This completes the fourth task in Remark 2.2.

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 [9], [19] 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 [9], and [19].
We start by rotating the velocities to be normal (u n = n x u + n y v) or aligned (u τ = −n y u + n x v) with the boundary. By inserting these transformation into (3.20) we obtain the rotated boundary contraction where the rotated solution is Φ r = (φ 1 , φ 1 u n , φ 1 u τ , φ 4 ) T . 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 , c is the speed of sound and M n = u n /c, the normal Mach number. The relations (3.22) indicate that for outflow (u n > 0) we have two situations. When M n > b there are 4 positive eigenvalues and no boundary condition is required. For M n < b, 3 eigenvalues are positive, 1 is negative and 1 boundary condition seem to be required. For inflow (u n < 0) we have the reversed situation with 4 negative eigenvalues and four required boundary conditions for M n > b, which goes down to 3 negative eigenvalues and 3 required boundary conditions for M n < b. Remark 3.1. The shift from subsonic to supersonic flow at M n = 1 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 b ≈ 0.976 if γ = 1.4, Maybe even more interesting is that b ≡ 1 for γ = √ 2 ≈ 1.414. However, considering eigenvalues is not sufficient for nonlinear problems [9], [19]. 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, see [26,27,28,29,30] for some previous results on this matter. This completes the fifth and final task in Remark 2.2.

The final form of the governing equations
The derivations above focused on stability and resulted in matricesÃ,B that were functions of the constant norm matrix P as seen in (3.6). The final governing equations can be transformed to where This removes the dependence of the norm P in the matrices which take the form (without free parameters) (3.27)

Some open questions
It is interesting to consider the energy (and entropy) rate. By combining (3.18) and (3.23) we find d dt This means the rate of change in the domain Ω of the energy (Φ T P Φ) is increasing or decreasing due to the transport of energy in or out of the domain plus an additional amount due to pressure work (γ − 1)u n p.
As we have seen theã 12 ,ã 14 ,b 13 ,b 14 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 α 2 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 p (and hence temperature). The first interpretation considers the problem in a physical way which means that ρ and p 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 φ 1 and φ 4 both ρ and p will always be positive. The bounds on the new variables directly lead to bounds on the original variables, by squaring them.

A stable energy and entropy conserving numerical approximation
To exemplify the straightforward construction of stable schemes based on the new formulation, we consider a summation-by-parts (SBP) approximation of (3.24),(3.25) as given in where Φ = ( Φ T 1 , Φ T 2 , ..., Φ T n ) T include approximations of Φ = (φ 1 , φ 2 , ..., φ n ) T in each node. The matrix elements of A 1 , A 2 , B 1 , B 2 are matrices with node values of the matrix elements in A 1 , A 2 , B 1 , B 2 injected on the diagonal as exemplified below Moreover D x = I n ⊗ D x ⊗ I y and D y = I n ⊗ I x ⊗ D y where D x,y = P −1 x,y Q x,y are 1D SBP difference operators, P x,y are positive definite diagonal quadrature matrices, Q x,y satisfies the SBP constraint Q x,y + Q T x,y = B x,y = diag[−1, 0, ..., 0, 1], ⊗ denotes the Kronecker product and I 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 U T B x V and U T B y V contain numerical integration along rectangular domain boundaries. In (4.3) we have usedP = I n ⊗ P x ⊗ P y , B x = (I n ⊗ B x ⊗ P y ) and B y = (I n ⊗ P x ⊗ B y ). The discrete energy method (multiply (4.1) from the left with 2 Φ T PP) where P = P ⊗ I x ⊗ I y yields where we have used that P commutes with D x , D y . Next, the discrete relations corresponding to (3.25), The semi-discrete energy rate in (4.7) mimics the continuous result in (3.18) and hence the scheme is energy and entropy 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 are all that is needed for stability.

Summary, conclusions and outlook
We have shown that a specific skew-symmetric form of linear and nonlinear problem leads to energy and entropy bounds for the compressible Euler equations. The skew-symmetric formulation automatically produced energy and entropy 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 2(γ − 1)/(γ(2 − γ)) for γ = 1.4. It also shoved that this ratio is identically one for γ = √ 2 ≈ 1.414. This paper together with [19] 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, easy to apply nonlinear schemes follows if summation-by parts operators are used for the discretisation. No additional requirements, such as chain rules or sign requirements are needed.
In future work we will continue the study of nonlinear boundary conditions, and especially it's relation to linear boundary procedures. We will also investigate presently unknown numerical advantages and disadvantages, such as stiffness, robustness, coarse mesh effects etc. In addition we will include dissipative effects, stemming from viscous terms in the Euler case, and bottom effects for the SWEs.