Two Decades Old Entropy Stable Method for the Euler Equations Revisited

The objective of this paper is to prove for the first time that the entropy split scheme of Olsson and Oliger (Energy and maximum norm estimates for nonlinear conservation laws. RIACS Technical Report 94.01, 1994; Gerritsen and Olsson, J Comput Phys 129:245–262, 1996; Yee et al., J Comput Phys 162:33–81, 2000) is the recent definition of an entropy stable method for central differencing with SBP operators for both periodic and non-periodic boundary conditions for nonlinear Euler equations. The proof is to replace the spatial derivatives by summation-by-parts (SBP) difference operators in the entropy split form of the equations using the physical entropy of the Euler equations. The numerical boundary closure follows directly from the SBP operator. No additional numerical boundary procedure is required. In contrast, Tadmor-type entropy conserving schemes Tadmor (Acta Numer 12:451–512, 2003) using mathematical entropies do not naturally come with a numerical boundary closure and a generalized SBP operator has to be developed Roanocha (J Comput Phys 362:20–48, 2018).

norm estimate readily without an added numerical dissipation term for smooth flows. For flows containing discontinuities the Yee et al. nonlinear filter approach [10-12, 14, 15, 22-25] is employed at isolated computed locations. After each full time step of the entropy split method to suppress spurious oscillations while maintaining accuracy on the remaining flow field. Since the nonlinear filter step is executed as an Euler time discretization at isolated location after the completion of a full time step of the entropy stable central scheme, entropy conservation/stability is valid almost everywhere. The efficiency and performance of the entropy stable split schemes using the physical entropies are compared with Tadmor-type entropy conservative method [18] using mathematical entropies for long time integration of a 2D smooth flows and a 3D direct numerical simulation (DNS) of turbulence with shocklets. It is found that Tadmor-type entropy conservative methods required twice the CPU time than the entropy stable split schemes using the same order of the central scheme. Comparisons among the three skew-symmetric splittings (entropy splitting [19,20,25], Ducros et al. splitting [1] and the Kennedy and Grubber splitting [5]) on their nonlinear stability and accuracy performance without added numerical dissipations for smooth flows is included. See [16] for additional details and comparison.
Remarks It is noted that the Hughes et al. formulation [4] using the Harten's idea [3] but solving the flow equations in nonconservative form in terms of the entropy variables is completely different from the entropy split schemes. The entropy split scheme solve the entropy splitting form of the Euler flux derivatives consisting of a one parameter family of conservative and a non-conservative portions in terms of the entropy variables. If the parameter satisfies the energy estimate, entropy stability is immediate. The entropy split scheme has been generalized from a perfect gas to a thermally perfect gas and gas flows consisting of linear combination of perfect gases [21,25]. In addition, these high order schemes have been formulated in time varying deforming curvilinear grids with free-stream preservation [17,21].

Entropy Splitting of the Euler Flux Derivatives
We consider the 3D equations of inviscid compressible gas dynamics with conserved variables q = (ρ ρu ρv ρw e) T and fluxes in an arbitrary direction k = (k 1 k 2 k 3 ) with |k| 2 = 1, and whereû = k 1 u+k 2 v +k 3 w. The total energy is related to the pressure p by the ideal gas law, e = p γ −1 + 1 2 ρ|u| 2 , where γ > 1 is a given constant, and |u| 2 = u 2 +v 2 +w 2 .
An entropy is a convex function, E(q), of the conserved variables that allows an additional conservation law, when the solution is smooth. The entropy fluxes in the x-, y-, and z-directions are denoted by F , G, and H , respectively. The entropy variables are defined by v = ∇ q E (the notation E q for the gradient will sometimes be used). The convexity of E ensures that these are well-defined. The Entropy conservation law (2) follows if the relation v T ∂f ∂q = ∇ q F for the x-direction fluxes, and similarly for the yand zdirections, holds. Moreover, the entropy variables symmetrize the equations; ∂f/∂v is a symmetric matrix.
Harten [3] considered the class of entropies where α is a parameter. To ensure that E is convex, i.e., that the matrix E q,q is positive definite, α is required to satisfy α > 0 or α < −γ . The full range for α was given in [25], while [3] only considered α > 0, and [2] used only the special case α = 1 − 2γ from α < −γ . The corresponding entropy flux in the direction where s denotes pρ −γ . The conserved variables are homogeneous functions of the entropy variables (4), See [3,16] for the proof. The range of α, where E q,q is positive definite, translates to β satisfying β < − γ γ −1 or β > 0. Entropy splitting of the Euler flux derivative in the x-direction with the yand z-directions suppressed [2,25] is written as a weighted sum of a conservative part, f x , and a non-conservative part, f v v x , as Replacing f x by this split flux derivative gives The entropy splitting weights the non-conservative portion of the flux derivative by 1 1+β . This means that the range β > 0 corresponds to a weight that is less than 1, whereas negative β leads, unphysically, to a weight that is greater than 1. The global entropy conservation can be rewritten as an L 2 -like estimate. The entropy time derivative can be rewritten as by using the homogeneity (5). Due to a page limit, see [16] for further discussion. Note that it is necessary to bound the eigenvalues of E −1 q,q in order to make the L 2like norm a valid estimate.

Semi-Discrete Entropy Split Discretization of the Euler Equations
Consider the 1D compressible gas dynamic equations discretized on a domain a < x < b by a uniform grid x j = (j − 1)Δx + a, j = 1, . . . , N, and grid spacing . Define the semi-discrete entropy split approximation where D is a SBP difference operator. With entropy split scheme, we will always mean the entropy split form of Eqs. (8) discretized in space by a summation-byparts finite difference operator. The flux Jacobian matrix with respect to the entropy variables, f v , is symmetric. The SBP scalar product is denoted by where ω j > 0 are weights that are different from 1 only at a few points near the boundaries. The operator D satisfies the SBP property but is otherwise arbitrary. In the most common case D is a standard SBP centered difference operator, but other operators are possible. A zero velocity, u 1 = 0, u N = 0, boundary condition is enforced, corresponding to wall boundaries. Thanks to the SBP property of the difference approximation the derivation of entropy conservation for the continuous problem can be carried over to the discretization.
The scheme (9) can be written where the projection P sets u 1 = 0 and u N = 0. Because P 2 = P , applying P to both sides of (11) gives that i.e., that P q = q if the initial data satisfy the boundary conditions. For the entropy where we use that P v = v. This is due to the second component of v is zero when the x-velocity, u, is zero, and the orthogonality (P v, (I − P )r) h = 0. The entropy equation is now of the same form as for the continuous problem, but replacing with integration-by-parts by summation-by-parts gives Entropy conservation follows by observing that F = uE, so that the boundary conditions imply that F 1 = F N = 0.
If the boundary conditions are periodic, no SBP modification of the difference operator is needed. Entropy conservation is proved with periodic boundary conditions by direct application of the same technique as above. It can be shown that the result carries over directly to the semi-discrete approximation, since only time derivatives are used in the proof. Hence, the L 2 -like estimate is obtained for the approximation (9). It can be shown that Tadmor-type entropy conservative discretization using the Harten entropy and high order central spatial differencings are also entropy conservative methods. See Sjögreen and Yee [16] for the proof.

Numerical Experiments
More extensive numerical experiments are reported in the extended version of this paper [16]. Previous studies using SBP boundary closures for non-periodic boundary conditions can be found in [25]. Here selected summary results are presented.

Test Case 1: 2D Compressible Euler Simulation of Smooth Flow: Isentropic Vortex Convection
The compressible Euler equations in two space dimensions are solved with initial data where r 2 = x 2 + y 2 , β = 5, γ = 1.4, u ∞ = 1, and v ∞ = 0. The exact solution is the initial data translated, u( The computational domain is 0 ≤ x ≤ 18, 0 ≤ y ≤ 18 with periodic boundary conditions. The center of the vortex is chosen to be (x 0 , y 0 ) = (9,9). The problem is solved in time with the classical fourth-order accurate explicit Runge-Kutta method to time t = 72, which corresponds to four revolutions of the vortex across the domain.
Comparisons of high order classical central split schemes with high order DRP schemes with grid refinements are reported in [13]. Due a space limitation only one grid with maximum and L 2 error norm compared with the exact solution is shown in Fig. 1. Here C08-DS represents eighth-order central differencing applied to the Ducros et al. splitting form of the Euler flux derivatives. The corresponding eighth-order entropy splitting, entropy conservative method and Kennedy Grubber splitting are indicated by "C08-ES", "C08-EC" and "C08-KGS". If the computed solutions by "C08-DS", "C08-ES", "C08-EC" and "C08-KGS" are nonlinearly filtered by a dissipative portion of WENO7 (seventh-order weighted essentially nonoscillatory spatial method) with an adaptive flow sensor, they are indicated by C08-DS+WENO7FI, C08-ES+WENO7FI, C08-EC+WENO7FI, and C08-KGS+WENO7FI [14,15,[22][23][24][25]. For the smooth flow without any turbulent structure, β = 1 for the entropy split scheme. The β parameter studies are reported in [9,16,25]. In general, for compressible shock-free turbulence and turbulence with shocklets, β lies somewhere in the range 1.5 < β < 2.5. In general, the optimal β is problem dependent. A general conclusion is that β should not be very large or below 1.
Other high resolution dissipative shock-capturing methods are also candidates for the nonlinear filter approach as well as other optimal WENO or ENO methods. However, with a good control of the numerical dissipation away from discontinuities, there is no need to use the more complicated and more CPU intensive shock-capturing methods. The non-split C08 without any added numerical dissipation diverges shortly after time evolution. Results by WENO5 or WENO7 are very diffusive with large maximum or L 2 errors. For this smooth long time integration flow, entropy splitting is the most accurate method.

Test Case 2: 3D Isotropic Turbulence with Eddy Shocklets
The second numerical test problem computes decaying compressible isotropic turbulence with eddy shocklets. For high enough turbulent Mach numbers weak shocks (shocklets) develop from the turbulent motion. Here the initial turbulent Mach number is 0.6. The Navier-Stokes equations are solved using γ = 1.4. The computational domain is a cube with side length 2π and periodic boundary conditions in all three directions. The initial datum is a random divergence free velocity field, u i,0 , i = 1, 2, 3, that satisfies The computations were made with u rms,0 = 1 and k 0 = 4. The angular brackets denote averaging over the entire computational domain. The density and pressure fields are initially constant. The Taylor-scale Reynolds number, Re λ,0 , is 100. See [6] for definitions of the quantities and more details about the set up of the problem. The simulation is run to the final time 4. Figure 2 shows the comparison of two splitting methods (DS and KGS), ES (entropy splitting and entropy stable) and EC (entropy conservative) using the same nonlinear filter. The time evolution of the domain averaged kinetic energy (upper left), enstrophy (upper right), temperature variance (lower left), and dilatation (lower right) are compared. All four forms of the nonlinear filter method provide similar resolution. All four schemes without the nonlinear filter are stable but not as accurate as the nonlinear filter versions. Over all, DS splitting is slightly less CPU intensive than ES. KGS skew-symmetric splitting is more CPU intensive than DS and ES. The EC method is around two times more expensive than DS. In addition, as the order of these methods increases, the gain in efficiency (CPU) by entropy split schemes increases. Although entropy split methods are not in conservation form but entropy conservative, Sect. 4 showed that they perform well on problems with shocklets. Over all, Extension of the entropy split scheme to other equations of state (nonperfect gas) and the MHD can be found in the original 2000 Yee et al. [25] paper. The entropies (3) can be used to construct entropy conserving schemes in conservative form. See [16] for the derivation.
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.