InstituteComputationalMathematics Comparison of some Entropy Conservative Numerical Fluxes for the Euler Equations Hendrik Ranocha 7 1 0 January 10, 2017 2 n a J Entropy stability is a well-known design principle for numerical methods in gas 9 dynamics. Entropy conservative numerical fluxes can be used as ingredients in two ] kinds of schemes: Firstly, as building blocks in the subcell flux differencing form of A Fisher and Carpenter (2013) and secondly (enhanced by dissipation) as numerical N surface fluxes in finite volume like schemes. In this article, the flux differencing . h theory is extended, guaranteeing high-order for symmetric and consistent numerical t fluxes. Additionally, an extension to simplices using a genuinely multi-dimensional a m framework of summation-by-parts operators is investigated for the first time in this [ framework. Moreover, several new entropy conservative numerical fluxes are devel- oped and compared extensively with existing ones, both in theory and in numerical 1 v tests, using the flux differencing and the finite volume framework. 4 6 2 1 Introduction 2 0 . Entropystabilityhaslongbeenusedasdesignprinciplefornumericalschemesforgasdynamics. 1 0 As an ingredient, entropy conservative numerical fluxes can be used in two kinds of application: 7 They can be used as volume fluxes in the flux differencing framework of Fisher and Carpenter 1 (2013a) and – enhanced with additional dissipation operators – as numerical fluxes in a finite : v volume framework. i X In this article, the theory of the flux differencing form by Fisher and Carpenter (2013a) is r extended. Thus, high order can be guaranteed not only for the entropy conservative flux of a Tadmor (1987), but also for every consistent and symmetric numerical flux, as observed in all numerical experiments known to the author. Secondly, for the first time, a genuinely multi- dimensional formulation of generalised summation-by-parts operators is used to investigate the extension of the framework to simplices. Afterwards, the construction of affordable entropy conservative fluxes is briefly reviewed and severalnewentropyconservativenumericalfluxesareconstructed–withandwithoutthekinetic energy preserving property of Jameson (2008). Finally,thenumericalfluxesareenhancedwithseveraldissipationoperatorsandarecompared both theoretically and in diverse test cases. Included in this comparison are the classical local Lax-Friedrichs flux and the Suliciu relaxation solver described by Bouchut (2004). Thisarticleisorganisedasfollows. Atfirst,somewell-knownpropertiesoftheEulerequations aresummedupinsection2inordertofixthenotationandforfurtherreference. Afterwards,the extensionofthefluxdifferencingtheoryofFisherandCarpenter(2013a)ispresentedinsection3. 1 Thereafter, several entropy conservative numerical fluxes are constructed in sections 4 and 5. To get numerical surface fluxes usable in finite volume methods, the addition of dissipation is discussed in section 6, especially with regard to positivity preservation. After that, the methods are tested and compared using several problems in section 7. Finally, the results are summed up in section 8, conclusions are drawn and remaining open problems are formulated. 2 Euler equations In this section, some well known properties of the Euler equations in two space dimensions are given in order to fix the notation and refer to them later. The Euler equations are       (cid:37) (cid:37)v (cid:37)v x y (cid:37)v   (cid:37)v2 +p   (cid:37)v v  ∂  x+∂  x +∂  x y  = 0, t(cid:37)vy x (cid:37)vxvy  y (cid:37)vy2+p  (1) (cid:37)e ((cid:37)e+p)v ((cid:37)e+p)v x y (cid:124) (cid:123)(cid:122) (cid:125) (cid:124) (cid:123)(cid:122) (cid:125) (cid:124) (cid:123)(cid:122) (cid:125) =u =fx(u) =fy(u) where (cid:37) is the density of the gas, v = (v ,v ) its speed, (cid:37)v the momentum, e the specific total x y energy, and p the pressure. The total energy (cid:37)e can be decomposed into the internal energy (cid:37)ε and the kinetic energy 1(cid:37)v2, i.e. 2 1 (cid:37)e = (cid:37)ε+ (cid:37)v2. (2) 2 For a perfect gas, (cid:18) (cid:19) 1 p = (cid:37)RT = (γ 1)(cid:37)ε = (γ 1) (cid:37)e (cid:37)v2 , (3) − − − 2 where R is the gas constant, T the (absolute) temperature, and γ the ratio of specific heats. For air, γ = 1.4 will be used, unless stated otherwise. The (mathematical) entropy (scaled by a constant for convenience, as chosen inter alia by Ismail and Roe (2009) and Chandrashekar (2013)) used is (cid:37)s U = , (4) −γ 1 − where the (physical) specific entropy is given by p s = log = logp γlog(cid:37). (5) (cid:37)γ − With the associated entropy flux (cid:37)s F = Uv = v, (6) −γ 1 − smooth solutions fulfil ∂ U +∂ F +∂ F = 0, and the entropy inequality t x x y y ∂ U +∂ F +∂ F 0 (7) t x x y y ≤ will be used as an additional admissibility criterion for weak solutions. For (cid:37),p > 0, the entropy U(u) is strictly convex, and the entropy variables (cid:32) (cid:33)T γ s (cid:37)v2 (cid:37)v (cid:37)v (cid:37) x y w = U (u) = , , , (8) (cid:48) γ 1 − γ 1 − 2p p p −p − − can be used interchangeably with the conservative variables u. The flux potentials ψ = (cid:37)v , ψ = (cid:37)v (9) x x y y (cid:0) (cid:1) fulfil ψ (w) = f u(w) and F = w f ψ . x(cid:48)/y x/y x/y · x/y − x/y 2 3 Flux differencing form Entropy conservative two-point numerical fluxes can be used to create a high-order method by the approach of Fisher and Carpenter (2013a,b). Using diagonal-norm nodal SBP bases including boundary nodes, they constructed high-order numerical fluxes from given symmetric and consistent two-point fluxes fvol. Using the derivative matrix D with entries D , the corresponding discretisation of the i,k divergence ∂ f at x in one dimension reads x i p (cid:88) 2D fvol, (10) i,k i,k k=0 wherefvol = fvol(u ,u ). Ifthetwo-pointfluxfvol usedistheentropyconservativeoneproposed i,k i k by Tadmor (1987), (cid:90) 1 (cid:16) (cid:17) fvol(u ,u ) = f u(cid:0)w(u )+t(w(u ) w(u ))(cid:1) dt, (11) i k i k i − 0 Fisher and Carpenter (2013a, Theorem 3.1) showed that the resulting discretisation (10) of ∂ f(x ) is of the same order as the SBP derivative operator D. x i Additionally, Fisher and Carpenter (2013a, Theorem 3.2) showed that by using a consistent and symmetric two-point flux fvol that is also entropy conservative in the sense of Tadmor (1987), i.e. [[w]] fvol [[ψ]] = 0, the discretisation (10) is also entropy conservative. Here and in · − the following, [[a]] = a a denotes the jump of a quantity. + − − In the following, extensions to generalised SBP operators in multiple dimensions as well as generalisations / variations of the Theorems 3.1 and 3.2 of Fisher and Carpenter (2013a) will be developed. 3.1 Multiple dimensions In order to apply the methods to multi-dimensional problems, it is common to use rectangular gridsandtensorproductbasesusingtheresultsoftheone-dimensionalcase. However,genuinely multi-dimensional SBP operators can be constructed as well and the extension of the flux differencingformwithitshigh-orderaswellasentropyconservingpropertiesforone-dimensional nodal diagonal-norm SBP operators by Fisher and Carpenter (2013a, Theorems 3.1 and 3.2) can be described as follows. Multi-dimensional SBP operators can be described in a numerical setting (more adopted to finite difference methods) or an analytical setting (more common in finite element / discontin- uous Galerkin methods). Here, the numerical setting of Hicken, Fern´andez, and Zingg (2016) will be described using the notation of the analytical setting of Ranocha (2016a, Section 6). An SBP operator on a d dimensional element Ω consists of the following components. Derivative operators D , i 1,...,d , approximating the partial derivative in the i-th i • ∈ { } coordinate direction. These are required to be exact for polynomials of degree p. ≤ A mass matrix M , approximating the L scalar product on Ω via 2 • (cid:90) uTM v = u, v u, v = uv, (12) (cid:104) (cid:105)M ≈ (cid:104) (cid:105)L2(Ω) Ω where u,v are functions on Ω and u,v their approximations in the SBP basis. A restriction operator R performing interpolation to the boundary ∂Ω of Ω. • A boundary mass matrix B approximating the L scalar product of ∂Ω via 2 • (cid:68) (cid:69) (cid:90) u TBv = u , v u , v = u v , (13) B B B B B ≈ (cid:104) B B(cid:105)L2(∂Ω) ∂Ω B B where u ,v are functions on ∂Ω and u ,v their approximations in the SBP basis. B B B B 3 Multiplication operators N , i 1,...,d , performing multiplication of functions on the i • ∈ { } boundary ∂Ω with the i-th component of the outer unit normal. Together, the restriction and boundary operators approximate • (cid:90) uTRTBN Rv uvn , (14) i i ≈ ∂Ω where n is the i-th component of the outer unit normal n, and this approximation has to i be exact for polynomials of degree p. ≤ Finally, the SBP property • M D +D TM = RTBN R (15) i i i has to be fulfilled, mimicking integration by parts (i.e. the divergence theorem) on a discrete level (cid:90) (cid:90) (cid:90) u(∂ v)+ (∂ u)v uTM D v+uTD TM v = uTRTBN Rv uvn . (16) i i i i i i ≈ ≈ Ω Ω ∂Ω Inthefollowing,theSBPbasesusedtorepresentfunctionsonthevolumeΩandtheboundary ∂Ωarenodalbases. Nonlinearoperationswillbeperformedpointwiseandoftenthemassmatrix M will be assumed to be diagonal, corresponding to a quadrature rule on Ω, as described inter alia by Hicken and Zingg (2013) and Hicken, Fern´andez, and Zingg (2016, Theorem 3.2). In one dimension, an entropy conservative numerical flux fvol in the sense of Tadmor (1987, 2003) has to fulfil (w w ) fvol(u ,u ) (ψ ψ ) = 0, where w are the entropy variables i k i k i k − · − − (cid:0) (cid:1) w = U (u) and ψ is the flux potential obeying ∂ ψ(w) = f u(w) . A natural extension to (cid:48) w (cid:80) multipledimensionsistousefluxesf inthej-thcoordinatedirectionsuchthatdivf = ∂ f j j j j (cid:0) (cid:1) and corresponding flux potentials ψ obeying ∂ ψ (w) = f u(w) . Then, the j-th entropy flux j w j j F = w f ψ fulfils j j j · − F (u) =w (u) f +w f (u) ψ (u) j(cid:48) (cid:48) · j · j(cid:48) − j(cid:48) =w (u) f +w f (u) ∂ ψ (u) w (u) (17) (cid:48) · j · j(cid:48) − w j · (cid:48) =U (u) f (u), (cid:48) · j(cid:48) since w = U (u) and ∂ ψ = f . Thus, smooth solutions u of (cid:48) w j j (cid:88) ∂ u+ ∂ f (u) = 0 (18) t j j j fulfil (cid:88) (cid:88) ∂ U =U (u) ∂ u = U (u) ∂ f (u) = U (u) f (u) ∂ u t (cid:48) · t − (cid:48) · j j − (cid:48) · j(cid:48) · j j j (19) (cid:88) (cid:88) = F (u) ∂ u = ∂ F (u), − j(cid:48) · j − j j j j and the entropy inequality (cid:88) ∂ U + ∂ F (u) 0 (20) t j j ≤ j can be used as additional admissibility criterion. The natural extension of the one-dimensional flux differencing form is then given as follows. Choose consistent and symmetric numerical volume fluxes fvol, j 1,...,d . Then, the j ∈ { } approximation of ∂ f at x is given by j j i (cid:20) (cid:21) (cid:88) 2 D fvol(u ,u ). (21) j j i k k i,k Thatis, droppingtheindexj, ithasthesameformasinonespacedimension, usingappropriate derivative matrices D and corresponding volume fluxes fvol consistent with the flux f in the j j j j-th coordinate direction. 4 3.2 Order of approximation ThebasicargumentofFisherandCarpenter(2013a,Theorem3.1)canbegeneralisedasfollows. Assume that for fixed w , the j-th component of the volume flux can be written as i 1 (cid:88) fvol(w ,w ) = f (w )+ f (w ) (w w )+ c (w w )α, (22) j i k j i 2 j(cid:48) i · k − i α k − i α 2 | |≥ wheremulti-indexnotationhasbeenusedandwdenotesanyvariable,i.e. conservativevariables, primitive variables, or entropy variables. Note that a general Taylor expansion of fvol(w , ) j i · around w would read i (cid:88) fvol(w ,w ) = fvol(w ,w )+∂ fvol(w ,w ) (w w )+ c (w w )α. (23) j i k j i i 2 j i i · k − i α k − i α 2 | |≥ Using this, the volume discretisation at x (21) can be rewritten (dropping the indices of the i derivative operator and the flux) as p p p p (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) 2D fvol(w ,w ) = 2D f(w )+ D f (w ) (w w )+ D c (w w )α. i,k i k i,k i i,k (cid:48) i k i i,k α k i · − − k=0 k=0 k=0 k=0 α 2 | |≥ (24) Since the derivative is exact for constants, i.e. D1 = 0, the first sum on the right hand side of (24) vanishes. By the same reason, the second sum can be rewritten as p p (cid:88) (cid:88) D f (w ) (w w ) = f (w ) D w (25) i,k (cid:48) i k i (cid:48) i i,k k · − · k=0 k=0 and is therefore of the desired order. Finally, the third sum in (24) is a higher order correction to the product rule. Writing (cid:18) (cid:19) (cid:88) (cid:88) (cid:88) α c (w w )α = c wβ( w )α β (26) α k − i α β k − i − α 2 α 2 β α | |≥ | |≥ ≤ for multi-indices α,β, it becomes (cid:18) (cid:19) p (cid:88) (cid:88) α (cid:88) c ( w )α βD wβ. (27) α β − i − i,k k α 2 β α k=0 | |≥ ≤ By the product rule, a smooth function w of x satisfies d d ∂xwβ = ∂x(cid:16)w1β1...wnβn(cid:17) = (cid:88)βjw1β1...wjβj−11wjβj−1wjβ+j+11...wnβn∂xwj = (cid:88)βjwβ−ej∂xwj, − j=1 j=1 (28) where e is the j-th unit vector, (e ) = δ . Thus, the third sum in (24) is an approximation of j j l jl the same order as the derivative matrix D to (cid:18) (cid:19) p (cid:88) (cid:88) α (cid:88) c ( w )α βD wβ α β − i − i,k k α 2 β α k=0 | |≥ ≤ (cid:18) (cid:19) d p ≈ (cid:88) cα (cid:88) αβ (−wi)α−β(cid:88)βjwiβ−ej (cid:88)Di,kwk,j (29) α 2 β α j=1 k=0 | |≥ ≤ d (cid:18) (cid:19) p = (cid:88) cα(cid:88)(cid:88) αβ βj(−1)α−βwiα−ej (cid:88)Di,kwk,j, α 2 j=1β α k=0 | |≥ (cid:124)≤ (cid:123)(cid:122) (cid:125) 5 where 1 is the vector with components 1 of the same size as w and w is the j-th component i k,j of the vector w approximating w at x = x . The sum depending on β vanishes, since k k (cid:18) (cid:19) (cid:18) (cid:19) (cid:88) α (cid:88) α ∂ ( 1+w)α = ∂ ( 1)α βwβ = ( 1)α ββ wβ ej wj − wj β − − β − − j − β α β α ≤ ≤ (30) (cid:18) (cid:19) w==10 = αj( 1+1)α−ej = (cid:88) α ( 1)α−ββj. ⇒ − β − β α ≤ Thus, (21) is an approximation of the same order as D to ∂ f (w) at x , if (22) is fulfilled. j xj j i Finally, the assumption (22) is true, if the numerical flux fvol is consistent and symmetric, i.e. w: fvol(w,w) = f(w) w ,w : fvol(w ,w ) = fvol(w ,w ). (31) 1 2 1 2 2 1 ∀ ∧ ∀ Denoting the partial derivative with respect to the l-th component of the second argument of fvol as ∂ fvol(w,w), m 2,l m fvol(w,w+δe ) fvol(w,w) ∂ fvol(w,w) =lim m l − m 2,l m δ 0 δ → (32) fvol(w+δe ,w) fvol(w,w) =lim m l − m = ∂ fvol(w,w), δ 0 δ 1,l m → (cid:16) (cid:17) and the gradient has the form fvol(w,w) = ∂ fvol(w,w),∂ fvol(w,w) . Thus, the ∇(·,l) m 2,l m 2,l m directional derivative in direction 1 (1,1)T is given by √2 fvol(w+ δ e ,w+ δ e ) fvol(w,w) 2 ∂ fvol(w,w) = 1 (1,1) fvol(w,w) = lim m √2 l √2 l − m √2 2,l m √2 ·∇(·,l) m δ 0 δ → f (w+ δ e ) f (w) = lim m √2 l − m = 1 lim fm(w+δel)−fm(w) = 1 ∂ f (w). l m δ 0 δ √2 δ 0 δ √2 → → (33) Therefore, ∂ fvol(w,w) = 1∂ f (w), as required. This proves the following generalisation of 2,l m 2 l m Theorem 3.1 of Fisher and Carpenter (2013a) Theorem 1. If the numerical flux fvol is smooth, consistent with f , and symmetric, the flux j j differencing form (21) is an approximation to ∂ f (w) of the same order as the SBP derivative xj j matrix D . j 3.3 Entropy conservation TheentropyconservationprovedbyFisherandCarpenter(2013a,Theorem3.2)canbeextended similarly — at least in a setting based on ’small’ elements and looking at entropy conservation across elements. An extension of the subcell entropy conservation property as proved by Fisher and Carpenter (2013a, Theorem 3.2) is not aspired here. Dropping the index j in (21), the following assumptions will be used. The volume flux fvol is symmetric, consistent with the flux f, and entropy conservative in • the sense (w w ) fvol = ψ ψ , where fvol = fvol(u ,u ), w are the entropy variables, i− k · i,k i− k i,k i k and ψ is the flux potential. A nodal SBP basis with diagonal mass matrix M is used. • TheboundaryoperatorRTBN R isdiagonal,i.e. thereareenoughnodesontheboundary • to get the required exactness of integration. These assumptions on the SBP bases are fulfilled for tensor product Lobatto-Legendre nodes. 6 (cid:82) Usingtheseassumptionsinasemidiscretesetting, therateofchangeofthetotalentropy U Ω is discretely approximated as (cid:90) d U wTM ∂ u, (34) t dt ≈ Ω and (dropping again the index j of the derivative operator) the flux difference form (21) yields due to the SBP property (15) (cid:88) (cid:104) (cid:105) (cid:88) (cid:104) (cid:105) 2w M D fvol = w M D +RTBN R DTM fvol. i· i,k i,k i· − i,k i,k (35) i,k i,k Since the mass matrix M is diagonal, the volume term can be written as (cid:88) (cid:104) (cid:105) (cid:88) w M D DTM fvol = (M D M D )w fvol i· − i,k i,k ii ik − kk ki i· i,k i,k i,k (36) (cid:88) = M D (w w ) fvol, ii ik i− k · i,k i,k where the indices i,j have been exchanged in the second part of the sum, using the symmetry of fvol. Then, by entropy conservation (w w ) fvol = ψ ψ , i− k · i,k i− k (cid:88) (cid:88) (cid:88) M D (w w ) fvol = M D (ψ ψ ) = M D ψ ii ik i− k · i,k ii ik i− k − ii ik k i,k i,k i,k (cid:88)(cid:104) (cid:105) (cid:88)(cid:104) (cid:105) = M D ψ = RTBN R DTM ψ − ik k − − ik k (37) i,k i,k (cid:88)(cid:104) (cid:105) = RTBN R ψ , k − ik i,k since the derivative D is exact for constants, i.e. D1 = 0. The boundary term can be written, using that RTBN R is diagonal, as (cid:88) (cid:104) (cid:105) (cid:88)(cid:104) (cid:105) w RTBN R fvol = RTBN R w fvol i· i,k i,k k,k k · k,k i,k k (38) (cid:88)(cid:104) (cid:105) = RTBN R w f , k k k,k · k since the volume flux fvol is consistent with the flux f. Therefore, the total expression becomes (cid:88) (cid:104) (cid:105) (cid:88)(cid:104) (cid:105) 2w M D fvol = RTBN R (w f ψ ) i· i,k i,k k,k(cid:124) k · (cid:123)k(cid:122)− k(cid:125) i,k k =Fk =(cid:88)(cid:104)RTBN R(cid:105) F (39) k i,k i,k =1TRTBN RF, since the entropy flux F is given by F = w f ψ. This results in a consistent discretisation of · − (cid:90) (cid:90) w ∂ f = F n (40) j j j j · Ω ∂Ω for each j and entropy conservation follows. This proves the following generalisation / variation of Theorem 3.2 of Fisher and Carpenter (2013a) Theorem 2. If the numerical fluxes fvol are consistent with f , symmetric, and entropy con- j j servative, the nodal mass matrix M is diagonal, and the boundary operators RTBN R are j diagonal, too, the flux differencing form (21) is entropy conservative. Remark 3. To the author’s knowledge, there are no known SBP operators on simplices in general with diagonal RTBN R. In the framework of Hicken, Fern´andez, and Zingg (2016), j this operator is called E and they mention (Remark 4 in section 4.2) that they have not been j able to get diagonal operators that are sufficiently accurate. However, using tensor products of Lobatto-Legendre nodes in cubes, these operators are diagonal. 7 4 Entropy conservative fluxes In the semidiscrete setting of Tadmor (1987, 2003), an entropy conservative numerical flux fnum has to obey [[w]] fnum,j [[ψ ]] = 0, (41) j · − where w are the entropy variables (8), fnum,j is the numerical flux in direction j, ψ is the flux j potential in direction j, and [[a]] = a a (42) + − − denotes the jump of a quantity. Since the flux f is the gradient of the potential ψ , i.e. f = j j j ∂ ψ , the condition (41) for an entropy conservative flux determines fnum,j as an appropriate w j mean value of f . Indeed, the entropy conservative flux proposed by Tadmor (1987, Equation j (4.6a)) has the form of an integral mean (cid:90) 1 (cid:16) (cid:17) fnum,j(w ,w ) = f u(cid:0)w +s(w w )(cid:1) ds. (43) + j + − s=0 − − − However, this integral mean value is difficult to compute in general. Tadmor (2003, Theorem 6.1) proposed another integral mean based on a piecewise linear path in phase space to compute an integral mean similar to (43). Nevertheless, another approach will be used here. Following the well-known proverb “Differentiation is mechanics, integration is art.”, the in- tegral mean can be exchanged by some kind of differential mean. Sadly, there is no differential mean value theorem giving some kind of numerical flux fulfilling (41) directly in general. How- ever, the mean value theorem can be used for scalar variables. In this way, the affordable, entropy conservative numerical fluxes of Ismail and Roe (2009) and Chandrashekar (2013) can be constructed using the same general approach. This general procedure can be described as Express the flux potentials ψ (9) and the entropy variables w (8) using the chosen set of • variables. Express the jumps of ψ,w as products of some mean values and jumps of the chosen • variables using some kind of product / chain rule as in the mean value theorem. There are several mean values that can be used for this task. The simplest one is the arithmetic mean a +a + a = − , (44) {{ }} 2 with corresponding product and chain rule [[ab]] = a [[b]]+ b [[a]], [[a2]] = 2 a [[a]]. (45) {{ }} {{ }} {{ }} This is enough to get some entropy conservative fluxes for the shallow water equations, since the entropy variables w and the flux potential ψ can both be expressed as polynomials in both the primitive variables and the entropy variables as described by Ranocha (2016b). However, this is not true for the Euler equations. Therefore, other means have to be used. Roe (2006) proposed the logarithmic mean a a + a = − − , (46) {{ }}log loga loga + − − described by Ismail and Roe (2009), including a numerically stable implementation. The corre- sponding chain rule reads as 1 [[loga]] = [[a]]. (47) a {{ }}log Another criterion for a numerical flux is the kinetic energy preservation. The kinetic energy 1(cid:37)v2 obeys (for smooth solutions) 2 (cid:18) (cid:19) (cid:18) (cid:19) (cid:18) (cid:19) 1 1 1 ∂ (cid:37)v2 +∂ (cid:37)v2v +∂ (cid:37)v2v +v ∂ p+v ∂ p = 0. (48) t x x y y x x y y 2 2 2 8 In order to mimic this behaviour discretely in one space dimension, Jameson (2008) formulated the condition fnum = v fnum+pnum, (49) (cid:37)v {{ }} (cid:37) wherepnum isaconsistentnumericalfluxapproximatingthepressure. However, everyconsistent numerical flux fnum can be written in this form if some differences are accepted, i.e. if pnum = (cid:37)v fnum v fnum is accepted as numerical approximation of the pressure. (cid:37)v −{{ }} (cid:37) In the following, some entropy conservative numerical fluxes are presented. The fluxes of Roe (2006) and Chandrashekar (2013) in section 4.1 and 4.2 are well-known in the literature while the other ones are new. (cid:113) (cid:113) 4.1 (cid:37), (cid:37)v,√(cid:37)p as variables p p The entropy conservative flux of Ismail and Roe (2009); Roe (2006) can be derived using the variables (cid:114) (cid:114) (cid:114) (cid:37) (cid:37) (cid:37) z := , z := v , z := v , z := √(cid:37)p. (50) 1 2 x 3 y 5 p p p In these variables, the flux potential (9) and the entropy variables (8) are given by ψ = z z , ψ = z z , (51) x 2 5 y 3 5 (cid:18) (cid:19)T γ s 1 1 w = z2 z2, z z , z z , z2 , s = (γ +1)logz (γ 1)logz . γ 1 − γ 1 − 2 2 − 2 3 1 2 1 3 − 1 − 1− − 5 − − (52) Thus, the jumps can be expressed using (45) and (47) as 1 1 1 γ +1 1 1 [[w ]] = [[s]] [[z2]] [[z2]] = [[logz ]]+[[logz ]] [[z2]] [[z2]] 1 − γ 1 − 2 2 − 2 3 γ 1 1 5 − 2 2 − 2 3 − − γ +1 1 1 = [[z ]]+ [[z ]] z [[z ]] z [[z ]], 1 5 2 2 3 3 γ 1 z z −{{ }} −{{ }} − {{ 1}}log {{ 5}}log [[w ]] =[[z z ]] = z [[z ]]+ z [[z ]], 2 1 2 {{ 1}} 2 {{ 2}} 1 (53) [[w ]] =[[z z ]] = z [[z ]]+ z [[z ]], 3 1 3 1 3 3 1 {{ }} {{ }} [[w ]] = [[z2]] = 2 z [[z ]], 4 − 1 − {{ 1}} 1 [[ψ ]] =[[z z ]] = z [[z ]]+ z [[z ]], x 2 5 2 5 5 2 {{ }} {{ }} [[ψ ]] =[[z z ]] = z [[z ]]+ z [[z ]], y 3 5 3 5 5 3 {{ }} {{ }} and the entropy conservation conditions [[w]] fnum,x/y [[ψ ]] = 0 (41) become x/y · − (cid:32) (cid:33) γ +1 1 0 = fnum,x+ z fnum,x+ z fnum,x 2 z fnum,x [[z ]] γ 1 z (cid:37) {{ 2}} (cid:37)vx {{ 3}} (cid:37)vy − {{ 1}} (cid:37)e 1 − {{ 1}}log (cid:16) (cid:17) (cid:16) (cid:17) + z fnum,x+ z fnum,x z [[z ]]+ z fnum,x+ z fnum,x [[z ]] −{{ 2}} (cid:37) {{ 1}} (cid:37)vx −{{ 5}} 2 −{{ 3}} (cid:37) {{ 1}} (cid:37)vy 3 (cid:32) (cid:33) 1 + fnum,x z [[z ]], z (cid:37) −{{ 2}} 5 {{ 5}}log (54) (cid:32) (cid:33) γ +1 1 0 = fnum,y + z fnum,y + z fnum,y 2 z fnum,y [[z ]] γ 1 z (cid:37) {{ 2}} (cid:37)vx {{ 3}} (cid:37)vy − {{ 1}} (cid:37)e 1 − {{ 1}}log (cid:16) (cid:17) (cid:16) (cid:17) + z fnum,y + z fnum,y [[z ]]+ z fnum,y + z fnum,y z [[z ]] −{{ 2}} (cid:37) {{ 1}} (cid:37)vx 2 −{{ 3}} (cid:37) {{ 1}} (cid:37)vy −{{ 5}} 3 (cid:32) (cid:33) 1 + fnum,y z [[z ]]. z (cid:37) −{{ 3}} 5 {{ 5}}log 9 Thus, the fluxes  f(cid:37)num,x = {{z2}}{{z5}}log,   z z f(cid:37)nvuxm,x = {{{{z21}}}}f(cid:37)num,x+ {{{{z51}}}}, fnum,x z f(cid:37)nvuym,x = {{{{z31}}}}f(cid:37)num,x,   1γ +1 1 1 z 1 z f(cid:37)neum,x = 2γ −1{{z1}}{{z1}}logf(cid:37)num,x+ 2{{{{z21}}}}f(cid:37)nvuxm,x+ 2{{{{z31}}}}f(cid:37)nvuym,x, (55)  f(cid:37)num,y = {{z3}}{{z5}}log,   z f(cid:37)nvuxm,y = {{{{z21}}}}f(cid:37)num,y, fnum,y z z f(cid:37)nvuym,y = {{{{z31}}}}f(cid:37)num,y + {{{{z51}}}},   1γ +1 1 1 z 1 z f(cid:37)neum,y = 2γ −1{{z1}}{{z1}}logf(cid:37)num,y + 2{{{{z21}}}}f(cid:37)nvuxm,y + 2{{{{z31}}}}f(cid:37)nvuym,y, proposed (in one space dimension) by Roe (2006) and Ismail and Roe (2009) can be seen to be entropy conservative and consistent. However, by this choice of variables z, the pressure influences the numerical density flux. As explained by Derigs, Winters, Gassner, and Walch (2017), this can lead to problems if there are discontinuities in the pressure, see also the numerical tests in section 7. 4.2 (cid:37),v,β as variables Using the inverse of the temperature 1 (cid:37) β = = , (56) 2RT 2p Chandrashekar (2013) derived some entropy conservative fluxes. The flux potential and the entropy variables are ψ = (cid:37)v , ψ = (cid:37)v , (57) x x y y (cid:18) (cid:19)T γ s p w = βv2, 2βv , 2βv , 2β , s = log = logβ (γ 1)log(cid:37) log2. γ 1 − γ 1 − x y − (cid:37)γ − − − − − − (58) 4.2.1 Variant 1 Writing the jumps using the chain rules (45) and (47) as 1 1 [[w ]] = [[s]] [[βv2]] = [[log(cid:37)]]+ [[logβ]] [[βv2]] 1 − γ 1 − γ 1 − − − 1 1 1 = [[(cid:37)]]+ [[β]] v2 [[β]] v2 [[β]] 2 β v [[v ]] 2 β v [[v ]], (cid:37) γ 1 β −{{ x}} −{{ y}} − {{ }}{{ x}} x − {{ }}{{ y}} y {{ }}log − {{ }}log [[w ]] =2[[βv ]] = 2 β [[v ]]+2 v [[β]], 2 x x x {{ }} {{ }} [[w ]] =2[[βv ]] = 2 β [[v ]]+2 v [[β]], 3 y y y {{ }} {{ }} [[w ]] = 2[[β]], 4 − [[ψ ]] =[[(cid:37)v ]] = (cid:37) [[v ]]+ v [[(cid:37)]], x x x x {{ }} {{ }} [[ψ ]] =[[(cid:37)v ]] = (cid:37) [[v ]]+ v [[(cid:37)]], y y y y {{ }} {{ }} (59) 10