Formulation of Entropy-Stable schemes for the multicomponent compressible Euler equations

Ayoub Gouasmi Karthik Duraisamy University of Michigan, Department of Aerospace Engineering, Ann Arbor, MI, USA Scott M. Murman NASA Ames Research Center, NASA Advanced Supercomputing Division, Moffett field, CA, USA

Abstract

In this work, Entropy-Stable (ES) schemes are formulated for the multicomponent compressible Euler equations. Entropy-conservative (EC) and ES fluxes are derived. Particular attention is paid to the limit case of zero partial densities where the structure required by ES schemes breaks down (the entropy variables are no longer defined). It is shown that while an EC flux is well-defined in this limit, a well-defined upwind ES flux requires appropriately averaged partial densities in the dissipation matrix. A similar result holds for the high-order TecNO reconstruction. However, this does not prevent the numerical solution from developing negative partial densities or internal energy. Numerical experiments were performed on one-dimensional and two-dimensional interface and shock-interface problems. The present scheme exactly preserves stationary interfaces. On moving interfaces, it produces spurious pressure oscillations typically observed with conservative schemes [Karni, J. Comput. Phys., 112 (1994) 1]. We find that these anomalies, which are not present in the single-component case, violate neither entropy stability nor a minimum principle of the specific entropy. Finally, we show that the scheme is able to reproduce the physical mechanisms of the two-dimensional shock-bubble interaction problem [Haas & Sturtevant, J. Fluid Mech. 181 (1987) 41, Quirk & Karni, J. Fluid Mech. 318 (1996) 129].

keywords:

Entropy-Stable \sepMulticomponent Euler equations \sepFinite volume method \sepNonlinear hyperbolic conservation law

^†^†journal:

1 Introduction

Entropy-Stable (ES) schemes Tadmor have been gaining interest over the past decade, especially in the context of under-resolved simulations of compressible turbulent flows using high-order accurate numerical methods HO ; Diosady ; Pazner ; Fernandez_ES . ES schemes are attractive because they provide stability in an integral sense: the total entropy of the discrete system can be made non-decreasing (2nd principle) or conserved in which case the scheme is termed Entropy-Conservative (EC). ES schemes were initially motivated by the fact that some hyperbolic systems of conservation laws:

\frac{\partial\mathbf{u}}{\partial t}+\frac{\partial\mathbf{f}}{\partial x}=0,

(1)

where $\mathbf{u}$ and $\mathbf{f}$ are the state and flux vectors, admit a convex extension Friedrichs ; Harten in the sense that they imply an additional conservation equation for an entropy $U$ :

\frac{\partial U}{\partial t}+\frac{\partial F}{\partial x}=0,

(2)

where $(U,\ F)=(U(\mathbf{u}),F(\mathbf{u}))\in\mathbb{R}^{2}$ is an entropy pair satisfying $U_{,\mathbf{u}}^{T}\mathbf{f}_{,\mathbf{u}}=F_{,\mathbf{u}}^{T}$ (for equation (1) to imply equation (2)) and $U$ is convex. Entropy functions play an important role in the stability analysis of systems of conservation laws (see Tadmor_acta and references therein). The entropy inequality:

\frac{\partial U}{\partial t}+\frac{\partial F}{\partial x}<0,

(3)

considered in the sense of distributions, is typically used to eliminate non-physical weak solutions. For the compressible Euler equations, this leads to the well known entropy conditions which must be satisfied across a shock Lax . Note however that the entropy inequality (3) alone is not enough to uniquely determine the correct weak solution NES1 ; NES2 ; NES3 ; Gouasmi1 .
Tadmor Tadmor introduced finite-volume discretizations of equation (1) which are locally consistent with either equation (2) or inequality (3) for a given entropy pair. This is achieved in two fundamental steps. First, an EC flux, which does not introduce any entropy production or loss, is derived by solving a scalar entropy conservation condition Tadmor . Second, a dissipation operator is added to the EC flux in order to produce entropy across discontinuities. Despite some developments Barth ; LeFloch ; Tadmor_acta , an issue that hindered the use of EC/ES schemes is the non-closed form of the first EC flux Tadmor (path integral in phase space) in the general nonlinear systems case. Roe Roe1 addressed it by showing, for the compressible Euler equations, how to solve the EC condition algebraically to obtain a closed-form EC flux. Additionally, he proposed an upwind-like ES dissipation operator Roe1 that has become a default choice when constructing ES schemes. From there, a number a developments followed, mostly focused on high-order discretizations Fisher ; Fried ; Fjordholm ; Pazner ; Fernandez_ES . To the best of the authors’ knowledge, the most concrete improvement ES schemes brought about in Computational Fluid Dynamics (CFD) is in improving the robustness of under-resolved turbulent flow calculations using high-order methods Diosady ; Pazner ; Fernandez_ES . Efforts are being spent on the development of code infrastructures which take advantage of the benefits of ES schemes. An example is the eddy solver eddy0 ; Diosady developed at NASA Ames Research Center for the simulation of turbulent separated flows. The present work is part of an effort to expand the field of ES schemes towards multi-physics applications eddy and identify relevant research directions.
In this work, we consider the multicomponent ( $N$ species) compressible Euler equations consisting of the conservation equations for partial densities, momentum and total energy. This system can be viewed as the compressible Euler equations (conservation of total mass, momentum and total energy) complemented with $N-1$ species conservation equations. This observation motivated early multicomponent schemes such as the one by Habbal et al. Habbal , where the Roe scheme Roe is applied to the Euler part of the equations and the $N-1$ remaining equations are treated separately. In Larrouturou Larrouturou , such approaches are termed uncoupled as opposed to fully coupled approaches, which treat the multicomponent system as a whole. An example of a fully coupled approach is the extension of the Roe scheme by Fernandez & Larroutouru Fernandez , and Abgrall Abgrall0 . It might be then tempting to use an existing entropy-stable scheme for the Euler part and evolve the remaining $N-1$ species equations with another scheme. A case in point can be found in Derigs et al. Derigs (section 3.8). While this approach has the benefit of simplicity (minimal programming and computational effort), it is lacking from a theoretical viewpoint. This approach implicitly assumes that the entropy of the single component system is an admissible entropy for the multicomponent system, meaning that it is a conserved quantity in the smooth regime and a convex function of the conserved variables. The necessity of a fully coupled approach to develop entropy-stable schemes, in the sense of Tadmor Tadmor , for multicomponent flows is motivated by mathematical and physical arguments.
At this juncture, it is important to note that the schemes we are interested in this work achieve entropy-stability in a specific way. That is to say that there is more than one way that a scheme can be made stable in an entropy sense, and hence be called ’entropy-stable’. Osher’s family of E-schemes Osher and Barth’s space-time discontinuous galerkin schemes Barth are conservative schemes which satisfy an entropy inequality, but their construction is different. There are also non-conservative schemes which can be designed to conserve or produce entropy Castro . Entropy stability can be understood in a different way as well. The scheme introduced by Ma et al. Ihme is termed entropy-stable because it enforces a minimum principle of the specific entropy, proved by Tadmor for the Euler equations Tadmor . In their work, stability is sought in a point-wise sense (a scheme which preserves the positivity of density and satisfies the minimum principle cannot crash in principle), not in an integral sense. Integral stability and point-wise stability are both important concepts, and in principle they do not imply each other. In either case, it is important to emphasize that the correct formulation of these schemes depends on the structure of the equations they are applied to. Chalot et al. Chalot and Giovangigli Giovangigli demonstrated that the multicomponent compressible Euler equations do possess the structure that ES schemes require. Additionally, a minimum principle of the mixture’s specific entropy was recently shown in Gouasmi et al. Gouasmi4 . It is not clear to the authors whether these results extend to the systems considered in Ihme .
Integral entropy stability is not the only trait of ES schemes. An interesting feature of Tadmor’s framework is that the amount of entropy produced by the scheme can be precisely quantified and controlled to some extent. Ismail & Roe’s work on the carbuncle problem Ismail ; IsmailThesis relies on this feature to delve into the question of appropriate entropy production across shocks. There is a sense of how entropy is being managed locally, and neither the choice of EC flux nor the choice of dissipation operator is unique. However, this freedom in the design of the scheme leaves the user the outstanding challenge of figuring out what is the right way to manage entropy in a discrete flow field. This is an extremely complex question, the scope of which is by far not limited to shocks, but also includes rarefactions Gouasmi1 , contact discontinuities (of relevance to this work) and the under-resolution of turbulence MurmanCTR . The construction of our ES scheme follows a classic procedure Tadmor ; Roe : from the definition of the hyperbolic system, the entropy pair and the corresponding entropy variables, we derive the potential flux functions, a baseline EC flux Tadmor ; Roe1 and a scaling matrix Merriam ; Barth needed to define the entropy-stable dissipation operator of Roe Roe1 . Hence, our work does not explore this question in depth, but it does provides a lucid view of how the current state of the art of ES schemes fares in multicomponent flows.
We have come across a few obstacles during the course of this work. The first one is theoretical: the structure required by ES schemes collapses when one of the partial densities is zero. The entropy $U$ , which we take as the opposite of the mixture’s thermodynamic entropy Chalot ; Giovangigli , is no longer convex and the entropy variables, which are key in constructing ES schemes, are no longer defined. Upon closer examination, we observed that the EC flux we derived is well-defined in this limit and still satisfies the Entropy Conservation condition. We also found that the dissipation operator remains defined provided that the averaged partial densities involved in the dissipation matrix are evaluated in a certain way. The second issue is that while the overall scheme is always defined, there is no guarantee that it will not produce negative densities or pressure, even at first order, at the next time step. Last but not least, numerical experiments showed that while the ES scheme can handle shocks and stationary contact discontinuities correctly, it fails to preserve pressure equilibrium and constant velocity when a moving interface is simulated. This is in fact a longtime and well-known Abgrall ; Karni failure of conservative finite-volume schemes on multicomponent flows (there are no such anomalies in the single component case). We found that these anomalies neither violate entropy-stability nor a minimum principle on the specific entropy of the mixture Gouasmi4 .
The paper is organized as follows. In section 2, we present the modeling assumptions, the governing equations and their symmetrization using the entropy variables Chalot ; Giovangigli . Section 3 is dedicated to the construction of ES schemes for multicomponent flows. We formulate an EC flux and an ES interface flux based on upwinding Roe1 ; Ismail and discuss their definition in the limit of vanishing partial densities. In section 4, we discuss how this limit impacts standard high-order ES discretizations. In section 5, the numerical scheme is tested on one-dimensional and two-dimensional interface and shock-interface problems.

2 Governing equations, entropy variables and symmetrization

The governing equations describe the conservation of species mass, momentum and total energy. In 1D, that is system (1) with the vector conserved variables $\mathbf{u}$ and the vector of fluxes $\mathbf{f}$ defined by:

\mathbf{u}:=\begin{bmatrix}\rho_{1}&\ldots&\rho_{N}&\rho u&\rho e^{t}\end{bmatrix}^{T},\ \mathbf{f}:=\begin{bmatrix}\rho_{1}u&\ldots&\rho_{N}u&\rho u^{2}+p&(\rho e^{t}+p)u\end{bmatrix}^{T}.

$\rho_{k}$ denotes the partial density of species k, $\rho:=\sum_{k=1}^{N}\rho_{k}$ denotes the total density, $u$ denotes the velocity and $e^{t}$ denotes the specific total energy. The pressure $p$ is given by the ideal gas law:

p:=\sum_{k=1}^{N}\rho_{k}r_{k}T,\ r_{k}=\frac{R}{m_{k}},

where $m_{k}$ is the molar mass of species k and $R$ is the gas constant. The temperature $T$ is determined by the internal energy $\rho e:=\rho e^{t}-(\rho u)^{2}/(2\rho)$ which in this work is modeled following a calorically perfect gas assumption:

\rho e=\sum_{k}\rho_{k}e_{k},\ e_{k}:=e_{0k}+c_{vk}T.

(4)

For species k, $e_{0k}$ is a constant and $c_{vk}$ is the constant volume specific heat. $T$ is computed by solving equation (4). An additional conservation equation Chalot ; Giovangigli can be derived from the governing equations:

\frac{\partial\rho s}{\partial t}+\frac{\partial\rho us}{\partial x}=0.

(5)

$\rho s$ is the thermodynamic entropy of the mixture given by:

\displaystyle\rho s:=\sum_{k=1}^{N}\rho_{k}s_{k},\ s_{k}:=c_{vk}\ln({T})-r_{k}\ln(\rho_{k}),

where $s$ denotes the specific entropy of the mixture and $s_{k}$ denotes the specific entropy of species $k$ . Equation (5) can be rewritten in the form of equation (2) with $(U,\ F)=(-\rho s,\ -\rho us)$ . For $\rho_{k}>0$ and $T>0$ , $U$ is a convex function of the conserved variables and $(U,\ F)$ is a valid entropy pair for the multicomponent system Chalot ; Giovangigli .
At this point, we draw the attention of the reader to the difference between thermodynamic entropy $\rho s$ and mathematical entropy $U$ , which is a general concept proper to hyperbolic PDEs. For the Burgers equation for instance, $U=\frac{1}{2}u^{2}$ is an entropy. In our context, the mathematical entropy is the opposite of the thermodynamic entropy, and the statement of integral entropy stability can be interpreted either as dissipation of (mathematical) entropy or as production of (thermodynamic) entropy. We adopt the latter throughout this paper.
For systems of conservation laws with an entropy pair, Mock Mock proved that the mapping $\mathbf{u}\rightarrow\mathbf{v}$ where $\mathbf{v}$ denotes the entropy variables defined as:

\mathbf{v}:=\bigg{(}\frac{\partial U}{\partial\mathbf{u}}\bigg{)}^{T},

symmetrizes the original system. This means that, starting from the quasi-linear form of (1):

\frac{\partial\mathbf{u}}{\partial t}+A\frac{\partial\mathbf{u}}{\partial x}=0,\ A:=\frac{\partial\mathbf{f}}{\partial\mathbf{u}},

(6)

where $A$ is the flux Jacobian, the change of variables $\mathbf{u}\rightarrow\mathbf{v}$ turns (6) into a system of the form:

H\frac{\partial\mathbf{v}}{\partial t}+B\frac{\partial\mathbf{u}}{\partial x}=0,\ H:=\frac{\partial\mathbf{u}}{\partial\mathbf{v}},\ B:=\frac{\partial\mathbf{f}}{\partial\mathbf{v}}=AH,

(7)

where $B$ and $H$ are symmetric, and $H$ is symmetric positive definite (such systems are called symmetric hyperbolic and are well appreciated in the analysis of PDEs Friedrichs ; Giovangigli ). The matrices $A$ and $H$ will be derived in section 3.3. For completeness, we mention that Godunov Godunov showed that the converse of Mock’s result holds as well, namely that if there exists a mapping which symmetrizes the system, then there exists a valid entropy pair.
The fact that the entropy variables symmetrize the system (or equivalently, the fact that the matrix $H$ symmetrizes the flux Jacobian $A$ from the right) will be a key component in constructing the ES dissipation operator in section 3.3. In addition to the entropy variables, the construction of ES schemes also involve potential functions $(\mathcal{U},\mathcal{F})$ defined by:

\mathcal{U}:=\mathbf{v}\cdot\mathbf{u}-U,\ \mathcal{F}:=\mathbf{v}\cdot\mathbf{f}-F.

The potential functions satisfy the relationships:

\mathbf{u}=\bigg{(}\frac{\partial\mathcal{U}}{\partial\mathbf{v}}\bigg{)}^{T},\ \mathbf{f}=\bigg{(}\frac{\partial\mathcal{F}}{\partial\mathbf{v}}\bigg{)}^{T}.

In order to derive them, we must first derive the entropy variables. Following Chalot ; Giovangigli , we use a chain rule:

\frac{\partial U}{\partial\mathbf{u}}=\frac{\partial U}{\partial Z}\bigg{(}\frac{\partial\mathbf{u}}{\partial Z}\bigg{)}^{-1},\ Z:=\begin{bmatrix}\rho_{1}&\ldots&\rho_{N}&u&T\end{bmatrix}^{T},

where $Z$ denotes the vector of primitive variables. The Gibbs identity can be written as:

Td(\rho s)=d(\rho e)-\sum_{k=1}^{N}g_{k}d\rho_{k},

(8)

where $g_{k}:=h_{k}-Ts_{k}$ is the Gibbs function of species k and $h_{k}:=e_{k}+r_{k}T$ is the specific enthalpy of species $k$ . We have:

d(\rho e)=\sum_{k=1}^{N}e_{k}d\rho_{k}+\rho c_{v}dT,\ \rho c_{v}:=\sum_{k=1}^{N}\rho_{k}c_{vk}.

(9)

Combining equations (9) and (8), one obtains:

d(\rho s)=\frac{1}{T}\bigg{(}\sum_{k=1}^{N}(e_{k}-g_{k})d\rho_{k}+\rho c_{v}dT\bigg{)}.

This gives:

\frac{\partial U}{\partial Z}=\frac{1}{T}\begin{bmatrix}(g_{1}-e_{1})&\ldots&(g_{N}-e_{N})&0&-\rho c_{v}\end{bmatrix}.

(10)

The Jacobian of the mapping $Z\rightarrow\mathbf{u}$ is given by:

\frac{\partial\mathbf{u}}{\partial Z}=\begin{bmatrix}1&&0&0&0\\ &\ddots&&\vdots&\vdots\\ 0&&1&0&0\\ u&\ldots&u&\rho&0\\ e_{1}+\frac{1}{2}u^{2}&\ldots&e_{N}+\frac{1}{2}u^{2}&\rho u&\rho c_{v}\end{bmatrix}.

(11)

The inverse of this matrix is given by:

\bigg{(}\frac{\partial\mathbf{u}}{\partial Z}\bigg{)}^{-1}=\begin{bmatrix}1&&0&0&0\\ &\ddots&&\vdots&\vdots\\ 0&&1&0&0\\ -u\rho^{-1}&\ldots&-u\rho^{-1}&\rho^{-1}&0\\ (\frac{1}{2}u^{2}-e_{1})(\rho c_{v})^{-1}&\ldots&(\frac{1}{2}u^{2}-e_{N})(\rho c_{v})^{-1}&-u(\rho c_{v})^{-1}&(\rho c_{v})^{-1}\end{bmatrix}

(12)

Combining equations (12) and (10) yields the entropy variables Chalot ; Giovangigli :

\mathbf{v}=\bigg{(}\frac{\partial U}{\partial\mathbf{u}}\bigg{)}^{T}=\frac{1}{T}\begin{bmatrix}g_{1}-\frac{1}{2}u^{2}&\ldots&g_{N}-\frac{1}{2}u^{2}&u&-1\end{bmatrix}^{T}.

(13)

From this expression, the potential functions $(\mathcal{U},\ \mathcal{F})$ can be derived:

\mathcal{U}=\sum_{k=1}^{N}r_{k}\rho_{k},\ \mathcal{F}=\sum_{k=1}^{N}r_{k}\rho_{k}u.

(14)

We conclude this section with two remarks:

Denote $\mathbf{v}=[v_{1,k}\ \dots\ v_{1,N}\ v_{2}\ v_{3}]^{T}$ . In order to derive the mapping from entropy variables to primitive variables, one can first compute the temperature, velocity, gibbs functions and specific entropies as:

T:=-\frac{1}{v_{3}},\ u:=-\frac{v_{2}}{v_{3}},\ g_{k}:=-\frac{v_{1,k}}{v_{3}}+\frac{1}{2}\bigg{(}\frac{v_{2}}{v_{3}}\bigg{)}^{2},\ s_{k}(T,\rho_{k})=c_{pk}-v_{1,k}+\frac{1}{2}\frac{v_{2}^{2}}{v_{3}}

The partial densities are inferred from the specific entropies:

\rho_{k}:=\exp{\bigg{(}\frac{m_{k}}{R}(c_{vk}ln(T)-s_{k})\bigg{)}}=\exp{\Bigg{(}\frac{m_{k}}{R}\bigg{(}-c_{vk}ln(v_{3})-c_{pk}+v_{1,k}-\frac{1}{2}\frac{v_{2}^{2}}{v_{3}}\bigg{)}\Bigg{)}}.

The requirement that $\rho_{k}>0$ and $T>0$ manifests in the definition of the entropy variables, which require the evaluation of $\ln(\rho_{k})$ and $\ln(T)$ . On the other hand, it is interesting to note that if one works with the entropy variables instead of the conservative variables, the requirement that $\rho_{k}>0$ and $T>0$ boils down to the single requirement that $v_{3}<0$ . The remaining entropy variables can be of any sign in principle. The authors are not aware of any physical consideration which would impose the sign of $g_{k}-\frac{1}{2}u^{2}$ , namely the difference between gibbs energy and kinetic energy.

In the compressible Euler case with $e_{0k}=0$ , it is easy to show, using the ideal gas law $p=\rho rT$ and the Mayer relation $c_{p}-c_{v}=r$ that the vector of entropy variables (13) simplifies to :

\mathbf{v}=r\begin{bmatrix}\frac{\gamma-\bar{s}}{\gamma-1}-\frac{\rho u^{2}}{2p}&\frac{\rho u}{p}&-\frac{\rho}{p}\end{bmatrix}^{T},\ \bar{s}=\ln p-\gamma\ln\rho-\ln r.

This differs by a constant factor $r$ from the expression that is usually used when designing entropy stable schemes for the compressible Euler equations (see Barth ; Roe1 for instance). The trivial difference comes from a different choice of entropy pair $(U,\ F)$ .

3 Formulation

The first subsection covers the definitions and fundamental results relevant to the construction of ES schemes. The following subsections cover the details relevant to the multicomponent system.

3.1 Background

Here, the subscript $j$ denotes the cell index, and the subscript $j+\frac{1}{2}$ refers to the interface between cell $j$ and $j+1$ . We also use the notation $[.]_{j+\frac{1}{2}}$ to refer to the jump between cell $j+1$ and $j$ .
EC and ES schemes are defined as follows:

Definition 3.1 (Tadmor ).

The finite volume scheme:

\frac{d}{dt}\mathbf{u}_{j}(t)+\frac{1}{\Delta x}(\mathbf{f}_{j+\frac{1}{2}}-\mathbf{f}_{j-\frac{1}{2}})=0,

(15)

where $\mathbf{f}_{j+\frac{1}{2}}$ is a consistent interface flux, is called entropy conservative if it satisfies the equation:

\frac{d}{dt}U(\mathbf{u}_{j}(t))+\frac{1}{\Delta x}[F_{j+\frac{1}{2}}-F_{j-\frac{1}{2}}]=0,

where $F_{j+\frac{1}{2}}$ is a consistent entropy interface flux, and entropy stable if it satisfies the inequality:

\frac{d}{dt}U(\mathbf{u}_{j}(t))+\frac{1}{\Delta x}[F_{j+\frac{1}{2}}-F_{j-\frac{1}{2}}]<0.

(16)

The first step in the construction of an ES scheme is to construct an EC scheme, following:

Theorem 3.1 (Tadmor ).

The finite volume scheme (15) is entropy conservative if and only if the interface flux $\mathbf{f}_{j+\frac{1}{2}}$ satisfies the condition:

[\mathbf{v}]_{j+\frac{1}{2}}\cdot\mathbf{f}_{j+\frac{1}{2}}=[\mathcal{F}]_{j+\frac{1}{2}},

(17)

where $\mathcal{F}_{j}=\mathcal{F}(\mathbf{u}_{j})$ , is the potential function evaluated at cell $j$ . In this case, $\mathbf{f}_{j+\frac{1}{2}}$ is called an entropy-conservative flux and the corresponding entropy flux $F_{j+\frac{1}{2}}$ is given by the formula:

F_{j+\frac{1}{2}}=\frac{1}{2}(\mathbf{v}_{j}+\mathbf{v}_{j+1})\cdot\mathbf{f}_{j+\frac{1}{2}}-\frac{1}{2}(\mathcal{F}_{j}+\mathcal{F}_{j+1}).

For scalar PDEs, $\mathbf{v}$ is a scalar and the entropy conservation condition (17) has only one solution $\mathbf{f}_{j+\frac{1}{2}}=[\mathcal{F}]_{j+\frac{1}{2}}/[\mathbf{v}]_{j+\frac{1}{2}}$ . For systems, equation (17) does not uniquely determine the interface flux. The first two EC fluxes proposed by Tadmor Tadmor ; Tadmor_acta have the inconvenience of either not having a closed-form Tadmor expression or being relatively expensive to compute Tadmor_acta . Using algebraic manipulations analogous to that of Roe , Roe Roe1 proposed a simple, closed-form EC flux for the Euler equations that is more popular. The EC flux we derive for the multicomponent system in section 3.2 uses the same method.
An ES scheme is built by complementing an EC scheme with appropriately defined dissipation operators as follows:

Theorem 3.2 (Tadmor ).

The finite-volume scheme (15) with the interface flux $\mathbf{f}_{j+\frac{1}{2}}$ defined as:

\mathbf{f}_{j+\frac{1}{2}}=\mathbf{f}^{*}_{j+\frac{1}{2}}-D_{j+\frac{1}{2}}[\mathbf{v}]_{j+\frac{1}{2}}.

where $f^{*}_{j+\frac{1}{2}}$ satisfies the entropy conservation condition (17), and $D_{j+\frac{1}{2}}$ is a positive definite matrix, satisfies

\frac{d}{dt}U(\mathbf{u}_{j}(t))+\frac{1}{\Delta x}[F_{j+\frac{1}{2}}-F_{j-\frac{1}{2}}]=-\mathcal{E}_{j}.

(18)

where $\mathcal{E}_{j}$ is given by:

\mathcal{E}_{j}=\frac{1}{2\Delta x}([\mathbf{v}]_{j+\frac{1}{2}}^{T}D_{j+\frac{1}{2}}[\mathbf{v}]_{j+\frac{1}{2}}+[\mathbf{v}]_{j-\frac{1}{2}}^{T}D_{j-\frac{1}{2}}[\mathbf{v}]_{j-\frac{1}{2}})>0,

(19)

and is therefore entropy stable. In this case, the interface entropy flux $F_{j+\frac{1}{2}}$ is given by:

F_{j+\frac{1}{2}}=F^{*}_{j+\frac{1}{2}}-\frac{1}{2}(\mathbf{v}_{j}+\mathbf{v}_{j+1})D_{j+\frac{1}{2}}[\mathbf{v}]_{j+\frac{1}{2}},

(20)

where $F^{*}_{j+\frac{1}{2}}$ is the entropy flux associated with $f^{*}_{j+\frac{1}{2}}$ .

Summing over all cells and assuming periodic boundary conditions leads to the integral entropy stability statement:

\frac{d}{dt}\sum_{j}U(\mathbf{u}_{j})=-\sum_{j}\mathcal{E}_{j}<0.

(21)

This integral stability statement is obtained as a consequence of the local statement (18), which in itself is not a stability statement. It does not necessarily imply for instance that in every cell $j$ :

\frac{d}{dt}U(\mathbf{u}_{j})<0.

Looking at equation (18), we see that the local variation of $U$ in time depends on both $\mathcal{E}_{j}$ , which is determined by the dissipation operator, and the interface entropy flux $F_{j+\frac{1}{2}}$ , which according to equation (20) is determined by both the EC flux and the dissipation operator.
The second principle of thermodynamics is often invoked to defend the consistency of ES schemes with physics. First, this principle only applies to closed systems, hence it can hardly be invoked when the boundary conditions are non-trivial. Second, the second principle is not a local statement. This means that it does not explicitly support the local inequality (16).

3.2 Entropy-conservative flux

In compact notation, the entropy conservation condition (17) writes:

[\mathbf{v}]\cdot\mathbf{f}^{*}=[\mathcal{F}],

(22)

where $\mathbf{f}^{*}=[f_{1,1}\ \dots\ f_{1,N}\ f_{2}\ f_{3}]$ denotes the interface flux. Define the set of algebraic variables:

\mathbf{z}=\begin{bmatrix}\rho_{1}&\ldots&\rho_{N}&u&\frac{1}{T}\end{bmatrix}=\begin{bmatrix}z_{1,1}&\ldots&z_{1,N}&z_{2}&z_{3}\end{bmatrix}.

The first step of Roe’s technique Roe1 consists in rewriting the jump terms in equation (22) as linear combinations of the jumps in the algebraic variables. For a given quantity $a$ , let $\overline{a}$ and $a^{ln}$ be the arithmetic and logarithmic averages, respectively:

\overline{a}=\frac{1}{2}(a_{L}+a_{R}),\ a^{ln}=\begin{cases}a_{L},\ \ \mbox{if $a_{L}=a_{R}$}\\ 0,\ \ \mbox{if $a_{L}=0$ or $a_{R}=0$}\\ \frac{a_{R}-a_{L}}{\ln(a_{R})-\ln(a_{L})},\ \ \mbox{else}\end{cases}

Note the identities:

[ab]=\overline{a}[b]+\overline{b}[a],\ \ \mbox{and}\ \ [\ln(a)]=[a]/a^{ln}.

The jump in the potential function can be rewritten as:

[\mathcal{F}]=\sum_{k=1}^{N}r_{k}[\rho_{k}u]=\bigg{(}\sum_{k=1}^{N}r_{k}\overline{z_{1,k}}\bigg{)}[z_{2}]+\sum_{k=1}^{N}r_{k}\overline{z_{2}}[z_{1,k}].

(23)

For the jump in entropy variables, we need to examine the first N components. For $1\leq k\leq N$ :

\displaystyle\frac{1}{T}(g_{k}-\frac{1}{2}u^{2})=\frac{e_{0k}}{T}+c_{vk}+r_{k}-c_{vk}\ln T+r_{k}\ln(\rho_{k})-\frac{1}{2}\frac{u^{2}}{T}=e_{0k}z_{3}+c_{vk}+r_{k}+c_{vk}\ln(z_{3})+r_{k}\ln(z_{1,k})-\frac{1}{2}z_{2}^{2}z_{3}.

The corresponding jumps then write:

\displaystyle\bigg{[}\frac{1}{T}(g_{k}-\frac{1}{2}u^{2})\bigg{]}=

\displaystyle\ [z_{1,k}]\frac{r_{k}}{z_{1,k}^{ln}}-[z_{2}]\overline{z_{2}}\ \overline{z_{3}}+[z_{3}](e_{0k}+\frac{c_{vk}}{z_{3}^{ln}}-\frac{1}{2}\overline{z_{2}^{2}})

(24)

The remaining jumps are given by:

\bigg{[}\frac{u}{T}\bigg{]}=\overline{z_{3}}[z_{2}]+\overline{z_{2}}[z_{3}],-\bigg{[}\frac{1}{T}\bigg{]}=-[z_{3}].

(25)

Using equations (23), (24) and (25), the entropy conservation condition (22) can be rewritten as the requirement that a linear combination of the jumps in the algebraic variables equals zero:

\sum_{k=1}^{N}[z_{1,k}]\bigg{(}\frac{r_{k}}{z_{1,k}^{ln}}f_{1,k}\bigg{)}+[z_{2}]\bigg{(}(-\overline{z_{3}}\ \overline{z_{2}})\sum_{k=1}^{N}f_{1,k}+\overline{z_{3}}f_{2}\bigg{)}+[z_{3}]\bigg{(}\sum_{k=1}^{N}(e_{0k}+c_{vk}\frac{1}{z_{3}^{ln}}-\frac{1}{2}\overline{z_{2}^{2}})f_{1,k}+\overline{z_{2}}f_{2}-f_{3}\bigg{)}=\\ [z_{2}]\bigg{(}\sum_{k=1}^{N}r_{k}\overline{z_{1,k}}\bigg{)}+\sum_{k=1}^{N}[z_{1,k}]r_{k}\overline{z_{2}}.

The second step turns this scalar condition into a system of $N+3$ equations by invoking the independence of these jumps:

\displaystyle r_{k}\frac{1}{z_{1,k}^{ln}}=r_{k}\overline{z_{2}}f_{1,k},\ 1\leq k\leq N,\ \ (-\overline{z_{3}}\ \overline{z_{2}})\sum_{k=1}^{N}f_{1,k}+\overline{z_{3}}f_{2}=\bigg{(}\sum_{k=1}^{N}r_{k}\overline{z_{1,k}}\bigg{)},\ \ \sum_{k=1}^{N}(e_{0k}+c_{vk}\frac{1}{z_{3}^{ln}}-\frac{1}{2}\overline{z_{2}^{2}})f_{1,k}+\overline{z_{2}}f_{2}-f_{3}=0.

The solution under these assumptions is:

	$\displaystyle f_{1,k}=$	$\displaystyle\ \overline{z_{2}}z_{1,k}^{ln},$
	$\displaystyle f_{2}=$	$\displaystyle\ \frac{1}{\overline{z_{3}}}\bigg{(}\sum_{k=1}^{N}r_{k}\overline{z_{1,k}}\bigg{)}+\overline{z_{2}}\sum_{k=1}^{N}f_{1,k},$
	$\displaystyle f_{3}=$	$\displaystyle\ \sum_{k=1}^{N}(e_{0k}+c_{vk}\frac{1}{z_{3}^{ln}}-\frac{1}{2}\overline{z_{2}^{2}})f_{1,k}+\overline{z_{2}}f_{2}.$

Therefore we obtain:

Theorem 3.3.

The interface flux $f^{*}=[f_{1,1}\ \dots\ f_{1,N}\ f_{2}\ f_{3}]$ defined by:

$\displaystyle f_{1,k}=$	$\displaystyle\ \rho_{k}^{ln}\overline{u},\ 1\leq k\leq N$
$\displaystyle f_{2}=$	$\displaystyle\ \frac{1}{\overline{1/T}}\bigg{(}\sum_{k=1}^{N}r_{k}\overline{\rho_{k}}\bigg{)}+\overline{u}\sum_{k=1}^{N}f_{1,k},$	(26)
$\displaystyle f_{3}=$	$\displaystyle\ \sum_{k=1}^{N}(e_{0k}+c_{vk}\frac{1}{(1/T)^{ln}}-\frac{1}{2}\overline{u^{2}})f_{1,k}+\overline{u}f_{2}.$

is entropy conservative for the multicomponent system defined in section 2.

This EC flux is the multicomponent version of Chandrasekhar’s EC flux Chandrasekhar for the compressible Euler equations. Chandrasekhar’s flux was designed with the additional property of being Kinetic-Energy Preserving (KEP) in the sense of Jameson Jameson , meaning that the kinetic energy equation is satisfied by the finite volume scheme (15) in a semi-discrete sense (in the same spirit as EC schemes). This property can be useful in turbulent flow simulations Subbareddy . For the compressible Euler equations, Jameson Jameson showed that this is achieved if the momentum flux $f^{\rho u}$ has the form $f^{\rho u}=\tilde{p}+\overline{u}f^{\rho}$ , where $f^{\rho}$ is the mass flux and $\tilde{p}$ is any consistent pressure average. The extension of Jameson’s analysis Jameson to the multicomponent case is straightforward and it can be showed that if the momentum flux has the same form as in the single component case (with $f^{\rho}$ denoting the total mass flux), the KEP property is achieved. The EC flux we derived qualifies, with $\tilde{p}$ given by:

\tilde{p}=\frac{1}{\overline{(1/T)}}\sum_{k=1}^{N}r_{k}\overline{\rho_{k}}.

Note in passing that if $[p]=0$ across the discontinuity, then $\tilde{p}$ defined above is exactly the pressure on each side.
The flux (3.3) is well-defined in the limit case $\rho_{k}=0$ , where the entropy variable $v_{1,k}$ is undefined, but is the entropy conservation condition still met? Equation (17) writes

\sum_{k=1}^{N}[v_{1,k}]f_{1,k}+[v_{2}]f_{2}+[v_{3}]f_{3}=[\mathcal{F}].

The jumps in $v_{1,k}$ are undefined because they involve jumps in $\ln(\rho_{k})$ . However, we note that the first term:

	$\displaystyle\sum_{k=1}^{N}[v_{1,k}]f_{1,k}=$	$\displaystyle\ \sum_{k=1}^{N}\bigg{[}\frac{h_{k}}{T}-\frac{1}{2}\frac{u^{2}}{T}\bigg{]}(\rho_{k}^{ln}\overline{u})\ -\ \sum_{k}[s_{k}](\rho_{k}^{ln}\overline{u})$
	$\displaystyle=$	$\displaystyle\ \sum_{k=1}^{N}\bigg{[}\frac{h_{k}}{T}-\frac{1}{2}\frac{u^{2}}{T}\bigg{]}(\rho_{k}^{ln}\overline{u})\ -\ \sum_{k}c_{vk}[\ln T](\rho_{k}^{ln}\overline{u})+\ \sum_{k=1}^{N}r_{k}[\ln\rho_{k}](\rho_{k}^{ln}\overline{u})$
	$\displaystyle=$	$\displaystyle\ \sum_{k=1}^{N}\bigg{[}\frac{h_{k}}{T}-\frac{1}{2}\frac{u^{2}}{T}\bigg{]}(\rho_{k}^{ln}\overline{u})\ -\ \sum_{k}c_{vk}[\ln T](\rho_{k}^{ln}\overline{u})+\ \sum_{k=1}^{N}r_{k}[\rho_{k}]\overline{u}$

is well-defined at $\rho_{k}=0$ , because the logarithmic averages $\rho_{k}^{ln}$ compensated for the logarithmic jumps $[\ln\rho_{k}]$ . We thus find that the entropy conservation condition is satisfied even in the limit $\rho_{k}=0$ .
Note that the EC flux does not transfer mass across an interface separating two different species.

3.3 Entropy-stable flux

3.3.1 Upwind-type dissipation operator

Let $A=R\Lambda R^{-1}$ be an eigendecomposition of $A$ . A popular choice for the dissipation operator consists of recasting the upwind operator of Roe Roe $\rightarrow\frac{1}{2}R|\Lambda|R^{-1}[\mathbf{u}]$ in terms of the entropy variables. With the differential relation $d\mathbf{u}=Hd\mathbf{v}$ , this leads to:

D[\mathbf{v}]=\frac{1}{2}R|\Lambda|R^{-1}H[\mathbf{v}].

For the compressible Euler equations, Merriam (Merriam , section 7.3) pointed out that there exists a scaling of the columns of $R$ such that $H=RR^{T}$ , which ultimately leads to a dissipation matrix $D=R|\Lambda|R^{T}$ which has the desirable property of being positive definite (Theorem 3.2). Barth Barth generalized this result to hyperbolic systems with a convex extension. The generalization takes the form of an eigenscaling theorem (Theorem 4 in Barth ) which states that for any diagonalizable matrix $A$ symmetrized on the right by a symmetric positive definite matrix $H$ , there exists a symmetric block diagonal matrix $T$ that block scales the eigenvectors $R$ of $A$ in such a manner that:

A=(RT)\Lambda(RT)^{T}\ \ \mbox{and}\ \ H=(RT)(RT)^{T}.

(27)

The dimensions of the blocks of $T$ correspond to the multiplicities of the eigenvalues of A. The second identity in equation (27) provides an explicit expression for the squared scaling matrix $T^{2}=R^{-1}HR^{-T}$ .
We now proceed to derive the scaling matrix for the multicomponent Euler system. First, the flux Jacobian is given by:

A=\frac{\partial\mathbf{f}}{\partial\mathbf{u}}=\begin{bmatrix}u(1-Y_{1})&\ldots&-uY_{1}&Y_{1}&0\\ &\ddots&&\vdots&0\\ -uY_{N}&\ldots&u(1-Y_{N})&Y_{N}&0\\ \frac{(\gamma-3)}{2}u^{2}+d_{1}&\ldots&\frac{(\gamma-3)}{2}u^{2}+d_{N}&u(3-\gamma)&\gamma-1\\ u(d_{1}-h^{t}+\frac{u^{2}}{2}(\gamma-1))&\ldots&u(d_{N}-h^{t}+\frac{u^{2}}{2}(\gamma-1))&h^{t}-u^{2}(\gamma-1)&u\gamma\end{bmatrix},

where $d_{k}=h_{k}-\gamma e_{k}$ and $h^{t}=h+\frac{1}{2}u^{2}$ . We have $A=R\Lambda R^{-1}$ with:

\displaystyle R=\begin{bmatrix}1&&0&Y_{1}&Y_{1}\\ &\ddots&&\vdots&\vdots\\ 0&&1&Y_{N}&Y_{N}\\ u&\ldots&u&u+a&u-a\\ k-\frac{d_{1}}{\gamma-1}&\ldots&k-\frac{d_{N}}{\gamma-1}&h^{t}+ua&h^{t}-ua\end{bmatrix},\ \Lambda=diag(u,\ \ldots,\ u,\ u+a,\ u-a),\ a^{2}=\gamma rT,

where $k=\frac{u^{2}}{2}$ . The Jacobian $H$ of the mapping $\mathbf{v}\rightarrow\mathbf{u}$ is given by Giovangigli :

H=\frac{\partial\mathbf{u}}{\partial\mathbf{v}}=\begin{bmatrix}\rho_{1}/r_{1}&&0&u\rho_{1}/r_{1}&\rho_{1}e_{1}^{t}/r_{1}\\ &\ddots&&\vdots&\vdots\\ 0&&\rho_{N}/r_{N}&u\rho_{N}/r_{N}&\rho_{N}e_{N}^{t}/r_{1}\\ u\rho_{1}/r_{1}&\ldots&u\rho_{N}/r_{N}&\rho T+u^{2}\sum_{k=1}^{N}\rho_{k}/r_{k}&u(\rho T+\sum_{k=1}^{N}\rho_{k}e_{k}^{t}/r_{k})\\ \rho_{1}e_{1}^{t}/r_{1}&\ldots&\rho_{N}e_{N}^{t}/r_{N}&u\rho T+u\sum_{k=1}^{N}\rho_{k}e_{k}^{t}/r_{k}&\rho T(u^{2}+c_{v}T)+\sum_{k=1}^{N}\rho_{k}(e_{k}^{t})^{2}/r_{k}\end{bmatrix},

(28)

where $e_{k}^{t}=e_{k}+k$ . The squared scaling matrix is given by:

T^{2}=R^{-1}HR^{-T}=\frac{\rho}{\gamma r}diag(T^{2Y},\ 1/2,\ 1/2),

(29)

where $T^{2Y}\in\mathbb{R}^{N\times N}$ is given by:

	$\displaystyle T^{2Y}_{ii}=$	$\displaystyle\ (\gamma-1)Y_{i}^{2}+\sum_{k\neq i}(\gamma r_{k}/r_{i})Y_{k}Y_{i},\ 1\leq i\leq N,$
	$\displaystyle T^{2Y}_{ij}=$	$\displaystyle\ -Y_{i}Y_{j},\ 1\leq i\neq j\leq N.$

At this point, the dissipation operator writes:

D[\mathbf{v}]=\frac{1}{2}R|\Lambda|T^{2}R^{T}[\mathbf{v}],

(30)

and qualifies for the ES scheme because the matrix $R|\Lambda|T^{2}R^{T}$ is positive definite ( $T^{2}$ and $|\Lambda|$ commute). However, as will be seen in section 4.2, a scaled form (27) of the dissipation operator is necessary.
For $N=2$ , we have:

T^{2Y}=\begin{bmatrix}(\gamma-1)Y_{1}^{2}+(\gamma r_{2}/r_{1})Y_{1}Y_{2}&-Y_{1}Y_{2}\\ -Y_{1}Y_{2}&(\gamma-1)Y_{2}^{2}+(\gamma r_{1}/r_{2})Y_{1}Y_{2}\end{bmatrix}.

$T^{2Y}$ is symmetric, real-valued with non-negative eigenvalues therefore there exists $T^{Y}$ with the same properties such that $T^{2Y}=(T^{Y})^{2}=T^{Y}(T^{Y})^{T}$ ( $T^{Y}$ is the square root of $T^{2Y}$ ). This matrix can be derived using an eigenvalue decomposition. However the expression of $T^{Y}$ is quite complicated. The square root of $T^{2Y}$ is not necessary to proceed. For $N=2$ , $T^{2Y}$ can be rewritten as:

T^{2Y}=\mathcal{T}^{Y}(\mathcal{T}^{Y})^{T},\ \mathcal{T}^{Y}=\begin{bmatrix}-\sqrt{Y_{1}Y_{2}}\sqrt{\gamma r_{2}/r_{1}}&Y_{1}\sqrt{\gamma-1}\\ \sqrt{Y_{1}Y_{2}}\sqrt{\gamma r_{1}/r_{2}}&Y_{2}\sqrt{\gamma-1}\end{bmatrix}.

(31)

$\mathcal{T}^{Y}$ is not the square root of $T^{2Y}$ , however it is enough to obtain a scaled formulation (27) because:

T^{2}=\mathcal{T}\mathcal{T}^{T},\ \mathcal{T}=\sqrt{\frac{\rho}{\gamma r}}diag(\mathcal{T}^{Y},\ 1/\sqrt{2},\ 1/\sqrt{2}),

and $\mathcal{T}$ commutes with $|\Lambda|$ therefore the dissipation operator can be rewritten as:

D[\mathbf{v}]=\frac{1}{2}(R\mathcal{T})|\Lambda|(R\mathcal{T})^{T}[\mathbf{v}].

(32)

For $N>2$ , the expression for $T^{Y}$ becomes even more complicated. For $N=3$ , the alternative (31) to $T^{Y}$ we proposed for $N=2$ becomes:

T^{2Y}=\mathcal{T}^{Y}(\mathcal{T}^{Y})^{T},\ \mathcal{T}^{Y}=\begin{bmatrix}-\sqrt{Y_{1}Y_{2}}\sqrt{\gamma r_{2}/r_{1}}&-\sqrt{Y_{1}Y_{3}}\sqrt{\gamma r_{3}/r_{1}}&0&-Y_{1}\sqrt{\gamma-1}\\ \sqrt{Y_{1}Y_{2}}\sqrt{\gamma r_{1}/r_{2}}&0&-\sqrt{Y_{2}Y_{3}}\sqrt{\gamma r_{3}/r_{2}}&-Y_{2}\sqrt{\gamma-1}\\ 0&\sqrt{Y_{1}Y_{3}}\sqrt{\gamma r_{1}/r_{3}}&\sqrt{Y_{2}Y_{3}}\sqrt{\gamma r_{2}/r_{3}}&-Y_{3}\sqrt{\gamma-1}\end{bmatrix}.

There is one more column compared to the $N=2$ case. The form (32) still holds except that $R\mathcal{T}\in\mathbb{R}^{3\times 4}$ instead of $\mathbb{R}^{3\times 3}$ and $|\Lambda|\in\mathbb{R}^{4\times 4}$ diagonal with an extra $|u|$ term. For $N$ species, the ”pseudo” scaling matrix $\mathcal{T}^{Y}$ we described will be in $\mathbb{R}^{N\times(N(N-1)/2+1)}$ .

3.3.2 Average state

To complete the definition of the dissipation operator, an average state (referred to with the superscript $*$ ) needs to be specified.
We showed in section 3.2 that the EC flux (3.3) is well defined in the limit $\rho_{k}=0$ . What about the dissipation operator $R|\Lambda|T^{2}R^{T}[\mathbf{v}]$ ? At first glance, the presence of the jump term $[\mathbf{v}]$ is problematic, because $[v_{1,k}]$ is undefined. For $N=2$ , the squared scaling matrix can be rewritten as:

T^{2}=\overline{T}^{2}\mathcal{R},\ \overline{T}^{2}=\frac{1}{\gamma r}\begin{bmatrix}(\gamma-1)Y_{1}+\gamma r_{2}/r_{1}Y_{2}&-Y_{2}&0&0\\ -Y_{1}&(\gamma-1)Y_{2}+\gamma r_{1}/r_{2}Y_{1}&0&0\\ 0&0&\frac{1}{2}&0\\ 0&0&0&\frac{1}{2}\end{bmatrix},\ \mathcal{R}=\begin{bmatrix}\rho_{1}&0&0&0\\ 0&\rho_{2}&0&0\\ 0&0&\rho&0\\ 0&0&0&\rho\end{bmatrix}.

(33)

Denoting $D_{k}=k-\frac{d_{k}}{\gamma-1}$ , we have:

\mathcal{R}R^{T}[\mathbf{v}]=\begin{bmatrix}\rho_{1}^{*}&0&0&0\\ 0&\rho_{2}^{*}&0&0\\ 0&0&\rho^{*}&0\\ 0&0&0&\rho^{*}\end{bmatrix}\begin{bmatrix}1&0&u^{*}&D_{1}^{*}\\ 0&1&u^{*}&D_{2}^{*}\\ Y_{1}^{*}&Y_{2}^{*}&u^{*}+a^{*}&(h^{t})^{*}+u^{*}a^{*}\\ Y_{1}^{*}&Y_{2}^{*}&u^{*}-a^{*}&(h^{t})^{*}-u^{*}a^{*}\end{bmatrix}\begin{bmatrix}\big{[}v_{1,1}\big{]}\\ \big{[}v_{1,2}\big{]}\\ \big{[}v_{2}\big{]}\\ \big{[}v_{3}\big{]}\end{bmatrix}=\\ \begin{bmatrix}\rho_{1}^{*}\big{[}v_{1,1}\big{]}+\rho_{1}^{*}u\big{[}v_{2}\big{]}+\rho_{1}^{*}D_{1}^{*}\big{[}v_{3}\big{]}\\ \rho_{2}^{*}\big{[}v_{1,2}\big{]}+\rho_{2}^{*}u\big{[}v_{2}\big{]}+\rho_{2}^{*}D_{2}^{*}\big{[}v_{3}\big{]}\\ \rho_{1}^{*}\big{[}v_{1,1}\big{]}+\rho_{2}^{*}\big{[}v_{1,2}\big{]}+\rho^{*}(u^{*}+a^{*})\big{[}v_{2}\big{]}+\rho^{*}((h^{t})^{*}+u^{*}a^{*})\big{[}v_{3}\big{]}\\ \rho_{1}^{*}\big{[}v_{1,1}\big{]}+\rho_{2}^{*}\big{[}v_{1,2}\big{]}+\rho^{*}(u^{*}-a^{*})\big{[}v_{2}\big{]}+\rho^{*}((h^{t})^{*}-u^{*}a^{*})\big{[}v_{3}\big{]}\end{bmatrix}

(34)

We can see that $\mathcal{R}R^{T}[\mathbf{v}]$ is well-defined if $\rho_{k}^{*}[v_{1,k}]$ is well-defined as well. Given that:

\rho_{k}^{*}[v_{1,k}]=\bigg{(}-c_{vk}[ln(T)]-\bigg{[}\frac{u^{2}}{2T}\bigg{]}\bigg{)}\rho_{k}^{*}+\frac{R}{m_{k}}[\ln\rho_{k}]\rho_{k}^{*},

we see that with $\rho_{k}^{*}=\rho_{k}^{ln}$ the dissipation operator is well-defined. For total density, one might be tempted to take $\rho^{*}=\sum_{k=1}^{N}\rho_{k}^{*}=\sum_{k=1}^{N}\rho_{k}^{ln}$ . This definition makes $\rho^{*}=0$ possible, which is undesirable given that $Y_{k}^{*}=\rho_{k}^{*}/\rho^{*}$ . $\rho^{*}=\rho^{\ln}$ is a safer choice.
The exact resolution of stationary contact discontinuities is a desirable property in the calculation of boundary and shear layers (even though it might produce carbuncles on blunt-body calculations, see Quirk_carbuncle paragraph 2.4.). In this case, $[p]=0,\ [u]=\overline{u}=0$ and the EC flux we derived reduces to:

f_{1,k}=0,\ f_{2}=\frac{1}{\overline{1/T}}\bigg{(}\sum_{k=1}^{N}\frac{R}{m_{k}}\overline{\rho_{k}}\bigg{)},\ f_{3}=0.

From the ideal gas law we can state that $f_{2}$ is exactly the pressure on both sides of the contact. The dissipation term must therefore cancel out if we want the ES scheme to exactly preserve stationary contact discontinuities.

Lemma 3.4.

The dissipation operator given by equation (30) vanishes at stationary contact discontinuities if the averaged state satisfies the relationship:

\rho^{*}h^{*}=\sum_{k=1}^{N}\rho_{k}^{*}\bigg{(}e_{0k}+c_{vk}\frac{1}{(1/T)^{ln}}\bigg{)}+\overline{p},\ \rho_{k}^{*}=\rho_{k}^{ln}.

(35)

Proof.

Since $u=0$ , we have:

\displaystyle R^{T}=\begin{bmatrix}1&&0&0&d_{1}^{*}/(\gamma-1)\\ &\ddots&&\vdots&\vdots\\ 0&&1&0&d_{N}^{*}/(\gamma-1)\\ Y_{1}^{*}&\ldots&Y_{N}^{*}&a^{*}&h^{*}\\ Y_{1}^{*}&\ldots&Y_{N}^{*}&-a^{*}&h^{*}\end{bmatrix},\ [\mathbf{v}]=\begin{bmatrix}\big{[}\frac{g_{1}}{T}\big{]}&\ldots&\big{[}\frac{g_{N}}{T}\big{]}&0&-\big{[}\frac{1}{T}\big{]}\end{bmatrix}^{T},

therefore:

R^{T}[\mathbf{v}]=\begin{bmatrix}\big{[}\frac{g_{1}}{T}\big{]}-d_{1}^{*}/(\gamma-1)\big{[}\frac{1}{T}\big{]}\\ \vdots\\ \big{[}\frac{g_{N}}{T}\big{]}-d_{N}^{*}/(\gamma-1)\big{[}\frac{1}{T}\big{]}\\ \sum_{k=1}^{N}Y_{k}^{*}\big{[}\frac{g_{k}}{T}\big{]}-\big{[}\frac{1}{T}\big{]}h^{*}\\ \sum_{k=1}^{N}Y_{k}^{*}\big{[}\frac{g_{k}}{T}\big{]}-\big{[}\frac{1}{T}\big{]}h^{*}\end{bmatrix}.

$|u^{*}|=0$ so the product of the eigenvalue matrix $|\Lambda|$ and the squared scaling matrix $T^{2}$ simplifies to:

|\Lambda|T^{2}=\frac{\rho^{*}a^{*}}{2\gamma^{*}r^{*}}\begin{bmatrix}0&0&0&0\\ 0&0&0&0\\ 0&0&1&0\\ 0&0&0&1\end{bmatrix}.

Therefore, the dissipation term $R|\Lambda|T^{2}R^{T}[\mathbf{v}]$ cancels out if:

\sum_{k=1}^{N}Y_{k}^{*}\bigg{[}\frac{g_{k}}{T}\bigg{]}-\bigg{[}\frac{1}{T}\bigg{]}h^{*}=0.

This is equivalent to:

\sum_{k=1}^{N}\rho_{k}^{*}\bigg{(}e_{0k}\bigg{[}\frac{1}{T}\bigg{]}-\big{[}s_{k}\big{]}\bigg{)}-\bigg{[}\frac{1}{T}\bigg{]}\rho^{*}h^{*}=0

(36)

Using:

-[s_{k}]=c_{vk}\bigg{[}\frac{1}{T}\bigg{]}\frac{1}{(1/T)^{ln}}+r_{k}\frac{[\rho_{k}]}{\rho^{ln}},

equation (36) can then be rewritten as:

\sum_{k=1}^{N}\bigg{(}\frac{\rho_{k}^{*}}{\rho_{k}^{ln}}\bigg{)}r_{k}[\rho_{k}]-\bigg{[}\frac{1}{T}\bigg{]}\bigg{(}\rho^{*}h^{*}-\sum_{k=1}^{N}\rho_{k}^{*}\bigg{(}e_{0k}+c_{vk}\frac{1}{(1/T)^{ln}}\bigg{)}\bigg{)}=0

(37)

The ideal gas law along with the assumption of constant pressure allows us to relate the jumps in partial densities and temperature:

\sum_{k=1}^{N}r_{k}[\rho_{k}]=\overline{p}\bigg{[}\frac{1}{T}\bigg{]}.

(38)

If $\rho_{k}^{*}=\rho_{k}^{ln}$ , then using equation (38), equation (37) simplifies to:

\bigg{[}\frac{1}{T}\bigg{]}\bigg{(}\bar{p}+\sum_{k=1}^{N}\rho_{k}^{*}\bigg{(}e_{0k}+c_{vk}\frac{1}{(1/T)^{ln}}\bigg{)}-\rho^{*}h^{*}\bigg{)}=0.

This leads to the condition (35) ∎

The remaining averages are taken as $u^{*}=\overline{u},\ T^{*}=1/(1/T)^{ln},\ r^{*}=\bar{r},\ \gamma^{*}=\overline{\gamma}$ and $a^{*}=\sqrt{\gamma^{*}r^{*}T^{*}}$ .

3.4 Additional considerations

3.4.1 Time integration

The first-order finite volume scheme we derived is ES at the semi-discrete level only. Entropy stability or entropy conservation at the fully discrete level can be obtained using a variety of techniques Barth ; LeFloch ; Tadmor_acta ; DeepThesis ; Diosady ; Gouasmi1 ; Gouasmi2 ; Fried which can be applied to the multicomponent compressible Euler system. However, this typically requires implicit time-integration schemes. For simplicity, we use explicit Runge-Kutta schemes in time, which do not guarantee entropy stability at the fully discrete level.
For the sake of completeness, is it worthwhile to recall the two fundamental results by Tadmor Tadmor_acta at the fully discrete level. Denote $\lambda=\Delta t/\Delta x$ , and consider the Forward Euler (FE) scheme in time applied to (15):

\mathbf{u}_{j}^{n+1}-\mathbf{u}_{j}^{n}+\lambda(\mathbf{f}_{j+\frac{1}{2}}^{n}-\mathbf{f}_{j-\frac{1}{2}}^{n})=0,

(39)

where the interface flux can be either EC or ES. Tadmor showed that (39) implies the following equation for entropy:

U(\mathbf{u}_{j}^{n+1})-U(\mathbf{u}_{j}^{n})+\lambda(F_{j+\frac{1}{2}}^{n}-F_{j-\frac{1}{2}}^{n})=-\mathcal{E}_{j}^{n}+\mathcal{E}_{j}^{FE},\ \mathcal{E}_{j}^{FE}>0.

(40)

If the interface flux is EC, $\mathcal{E}_{j}^{n}=0$ and the fully discrete scheme violates the entropy inequality in every cell. The numerical solution will develop non-physical oscillations growing in time. If the interface flux is ES, then the fully discrete entropy stability depends on the relative magnitudes of $\mathcal{E}^{n}_{j}$ (entropy production in space) and $\mathcal{E}^{FE}$ (entropy loss in time). In the scalar case, time step conditions (similar to CFL conditions) such that the fully discrete scheme is entropy-stable can be derived Zakerzadeh . The systems case is more complicated. In practice, a small enough time step is enough to ensure that the entropy losses in time do not overtake the entropy produced in space.
Tadmor also showed that the Backward Euler (BE) scheme in time applied to (15)

\mathbf{u}_{j}^{n+1}-\mathbf{u}_{j}^{n}+\lambda(\mathbf{f}_{j+\frac{1}{2}}^{n+1}-\mathbf{f}_{j-\frac{1}{2}}^{n+1})=0,

(41)

implies:

U(\mathbf{u}_{j}^{n+1})-U(\mathbf{u}_{j}^{n})+\lambda(F_{j+\frac{1}{2}}^{n+1}-F_{j-\frac{1}{2}}^{n+1})=-\mathcal{E}_{j}^{n+1}-\mathcal{E}_{j}^{BE},\ \mathcal{E}_{j}^{BE}>0.

(42)

In this case, the entropy inequality is always satisfied, even if an EC flux is used in space. The expressions for $\mathcal{E}_{j}^{BE}$ and $\mathcal{E}_{j}^{FE}$ can be found in Tadmor_acta (section 7).

3.4.2 Positivity

We showed in sections 3.2 and 3.3 that despite the fact that the entropy variables are undefined in the limit $\rho_{k}\rightarrow 0$ , the EC flux is well-defined and that the ES flux remains well-defined if the averaged partial densities are properly chosen (we show a similar result for high-order TecNO schemes in section 4.2). This does not guarantee that the resulting scheme will not produce negative densities and/or negative pressures. In the one-dimensional test problems, the first order semi-discrete ES scheme fortunately did not produce negative densities or pressure but the high-order TecNO scheme (section 4.2) systematically did. In the original setup of the two-dimensional shock-bubble interaction problem SBI_Marquina ; Quirk , the first-order scheme had the same issue.
There are schemes which are conservative, entropy stable and can preserve the correct sign of density and pressure. The Godunov scheme as well as the Lax-Friedrichs scheme qualify if an appropriate CFL condition is met Tadmor2 . Note also the recent work of Guermond et al. Guermond_IDP ; Guermond_IDP2 . A common trait of these schemes is that they take root in the notion of a Riemann problem and the assumption that there exists a solution satisfying all entropy inequalities. These schemes also require an algorithm capable of computing the maximum speed of propagation (or an upper bound) for the Riemann problem at each interface. One such algorithm is discussed in Guermond & Popov Guermond_speed . We could in principle adopt a hybrid approach where these schemes are used in areas where ours fails to maintain positive densities and pressure. Whether this can effectively be accomplished is left for future work.

3.4.3 Construction for thermally perfect gases

Here we briefly discuss the construction of an ES scheme in the case where the specific heats are not constant but functions of temperature. In this configuration, the specific internal energies and entropies are defined as:

e_{k}:=e_{0k}+\int_{0}^{T}c_{vk}(\tau)d\tau,\ s_{k}:=\int_{0}^{T}\frac{c_{vk}(\tau)}{\tau}d\tau-r_{k}\ln\rho_{k}.

The structure that ES schemes build on Giovangigli ; Chalot is still present. The multicomponent system still admits an additional conservation equation for the thermodynamic entropy of the mixture $\rho s$ and $(U,F)=(-\rho s,\ -\rho us)$ is a valid entropy pair for $\rho_{k}>0,\ T>0$ . In practice, the specific heats are represented using polynomial interpolation:

c_{vk}:=c_{k0}+\sum_{j=1}^{P}c_{kj}T^{j},

where $c_{kj},\ 0\leq j\leq P$ are constants. This gives:

e_{k}=e_{0k}+c_{k0}T+\sum_{j=1}^{p}\frac{c_{kj}}{j+1}T^{j+1},\ s_{k}=c_{k0}\ln T+\sum_{j=1}^{P}\frac{c_{kj}}{j}T^{j}-r_{k}\ln\rho_{k}.

The expressions of the entropy variables and potential functions remain unchanged. Regarding the construction of the EC flux, we have for $1\leq k\leq N$ :

[v_{1,k}]=\bigg{[}\frac{1}{T}\bigg{(}g_{k}-\frac{u^{2}}{2}\bigg{)}\bigg{]}=\bigg{[}\frac{e_{k}}{T}-s_{k}\bigg{]}-\bigg{[}\frac{u^{2}}{2T}\bigg{]}=e_{0k}\bigg{[}\frac{1}{T}\bigg{]}-c_{k0}[\ln T]-\sum_{j=1}^{P}\frac{c_{kj}}{j(j+1)}[T^{j}]+r_{k}[\ln\rho_{k}]-\frac{\overline{u^{2}}}{2}\bigg{[}\frac{1}{T}\bigg{]}-\overline{\frac{1}{T}}\overline{u}[u]

(43)

For a given quantity $a$ , let’s define the product operator $a^{\times}=a_{L}a_{R}$ . For $a\neq 0$ , we have the following identity:

[a]=-a^{\times}\ \bigg{[}\frac{1}{a}\bigg{]}

It can be easily showed, by induction for instance, that for $j\geq 1$ , there exists an averaging operator $f_{j}(a)$ , consistent with $a^{j-1}$ , such that $[a^{j}]=jf_{j}(a)[a]$ ( $f_{1}(a)=1$ , $f_{2}(a)=\overline{a}$ , $f_{3}(a)=(2/3)\overline{a}\ \overline{a}+(1/3)\overline{a^{2}}$ and so on using product rules). Equation (43) can therefore be rewritten as:

	$\displaystyle[v_{1,k}]=$	$\displaystyle\ e_{0k}\bigg{[}\frac{1}{T}\bigg{]}-c_{k0}[\ln T]+\sum_{j=1}^{P}\frac{c_{kj}}{(j+1)}(f_{j}(T)T^{\times})\bigg{[}\frac{1}{T}\bigg{]}+r_{k}[\ln\rho_{k}]-\frac{\overline{u^{2}}}{2}\bigg{[}\frac{1}{T}\bigg{]}-\overline{\frac{1}{T}}\overline{u}[u]$
	$\displaystyle=$	$\displaystyle\ \bigg{[}\frac{1}{T}\bigg{]}\bigg{(}e_{0k}+c_{k0}\frac{1}{(1/T)^{ln}}+\sum_{j=1}^{P}\frac{c_{kj}}{(j+1)}(f_{j}(T)T^{\times})-\frac{\overline{u^{2}}}{2}\bigg{)}+[\rho_{k}]\frac{r_{k}}{\rho_{k}^{ln}}-[u]\overline{\bigg{(}\frac{1}{T}\bigg{)}}\overline{u}$

Repeating the procedure outlined in section 3.2, we obtain an EC flux:

$\displaystyle f_{1,k}=$	$\displaystyle\ \rho_{k}^{ln}\overline{u},$
$\displaystyle f_{2}=$	$\displaystyle\ \frac{1}{\overline{1/T}}\bigg{(}\sum_{k=1}^{N}r_{k}\overline{\rho_{k}}\bigg{)}+\overline{u}\sum_{k=1}^{N}f_{1,k},$	(44)
$\displaystyle f_{3}=$	$\displaystyle\ \sum_{k=1}^{N}\bigg{(}e_{0k}+c_{k0}\frac{1}{(1/T)^{ln}}+\sum_{j=1}^{P}\frac{c_{kj}}{(j+1)}(f_{j}(T)T^{\times})-\frac{1}{2}\overline{u^{2}}\bigg{)}f_{1,k}+\overline{u}f_{2}.$

which is consistent ( $f_{j}(T)T^{\times}$ is consistent with $T^{j+1}$ ) and differs from the expression we obtained in the calorically perfect gas case (equation (3.3)) in the total energy component only.
We will not go into the details of the upwind dissipation operator. Barth’s eigenscaling theorem applies because $H$ still symmetrizes $A$ from the right. Therefore the scaling matrix exists ( $T^{2}=R^{-1}HR^{-T}$ ) and the dissipation operator constructed in section 3.3 can be constructed for thermally perfect gases as well.

3.4.4 Total mass form

Instead of solving the conservation equations for the N partial densities $\rho_{k}$ , one might want to solve for the conservation of the total density $\rho$ and N-1 partial densities. The state and flux vectors are then:

\mathbf{\tilde{u}}:=\begin{bmatrix}\rho&\rho_{2}&\ldots&\rho_{N}&\rho u&\rho e^{t}\end{bmatrix}^{T},\ \mathbf{\tilde{f}}:=\begin{bmatrix}\rho u&\rho_{2}u&\ldots&\rho_{N}u&\rho u^{2}+p&(\rho e^{t}+p)u\end{bmatrix}^{T}.

The entropy variables in this configuration can be easily obtained using a chain rule:

\mathbf{\tilde{v}}^{T}=-\frac{\partial\rho s}{\partial\mathbf{\tilde{u}}}=-\frac{\partial\rho s}{\partial\mathbf{u}}\bigg{(}\frac{\partial\mathbf{\tilde{u}}}{\partial\mathbf{u}}\bigg{)}^{-1}=\mathbf{v}^{T}\bigg{(}\frac{\partial\mathbf{\tilde{u}}}{\partial\mathbf{u}}\bigg{)}^{-1}=\mathbf{v}^{T}\bigg{(}\frac{\partial\mathbf{u}}{\partial\mathbf{\tilde{u}}}\bigg{)}

$\mathbf{\tilde{u}}$ and $\mathbf{u}$ only differ in the first component, $\rho_{1}=\rho-\sum_{k=2}^{N}\rho_{k}$ therefore:

\frac{\partial\mathbf{u}}{\partial\mathbf{\tilde{u}}}=\begin{bmatrix}1&-1&\ldots&-1&0&0\\ &1&&0&0&0\\ &&\ddots&0&0&0\\ 0&0&\ldots&1&0&0\\ 0&0&\ldots&0&1&0\\ 0&0&\ldots&0&0&1\end{bmatrix},\ \mathbf{\tilde{v}}=\frac{1}{T}\begin{bmatrix}(g_{1}-\frac{1}{2}u^{2})&(g_{2}-g_{1})&\ldots&(g_{N}-g_{1})&u&-1\end{bmatrix}^{T}.

The corresponding potential functions are unchanged because:

\mathbf{\tilde{v}}\cdot\mathbf{\tilde{f}}=\mathbf{v}^{T}\bigg{(}\frac{\partial\mathbf{u}}{\partial\mathbf{\tilde{u}}}\bigg{)}\bigg{(}\frac{\partial\mathbf{\tilde{u}}}{\partial\mathbf{u}}\bigg{)}\mathbf{f}=\mathbf{v}\cdot\mathbf{f},\ \mathbf{\tilde{v}}\cdot\mathbf{\tilde{u}}=\mathbf{v}^{T}\bigg{(}\frac{\partial\mathbf{u}}{\partial\mathbf{\tilde{u}}}\bigg{)}\bigg{(}\frac{\partial\mathbf{\tilde{u}}}{\partial\mathbf{u}}\bigg{)}\mathbf{u}=\mathbf{v}\cdot\mathbf{u}

Accordingly, an EC flux $\mathbf{\tilde{f}_{EC}}$ for the total mass form can simply be obtained by mapping an EC flux in the first form $\mathbf{f_{EC}}$ :

\mathbf{\tilde{f}_{EC}}=\bigg{(}\frac{\partial\mathbf{\tilde{u}}}{\partial\mathbf{u}}\bigg{)}\mathbf{f}_{EC}.

Likewise, applying the same mapping to an ES flux $\mathbf{f^{*}}$ in the first form given by:

\mathbf{f^{*}}=\mathbf{f_{EC}}-D[\mathbf{v}]

where $D$ is positive definite, results in a flux $\mathbf{\tilde{f}^{*}}$ given by:

\mathbf{\tilde{f}^{*}}=\bigg{(}\frac{\partial\mathbf{\tilde{u}}}{\partial\mathbf{u}}\bigg{)}\mathbf{f^{*}}=\mathbf{\tilde{f}_{EC}}-\bigg{(}\frac{\partial\mathbf{\tilde{u}}}{\partial\mathbf{u}}\bigg{)}D[\mathbf{v}]=\mathbf{\tilde{f}_{EC}}-\bigg{(}\frac{\partial\mathbf{\tilde{u}}}{\partial\mathbf{u}}\bigg{)}D\bigg{(}\frac{\partial\mathbf{\tilde{u}}}{\partial\mathbf{u}}\bigg{)}^{T}[\mathbf{\tilde{v}}]=\mathbf{\tilde{f}_{EC}}-\tilde{D}[\mathbf{\tilde{v}}].

$\tilde{D}$ is positive definite by congruence therefore $\mathbf{\tilde{f}^{*}}$ is entropy stable.

4 High-order discretizations

In this section, we are essentially interested in how the fundamental issues highlighted in the previous sections manifest in a high-order setting. We examine two high-order ES formulations: TecNO schemes Fjordholm and Discontinuous Galerkin Reed ; Cockburn (DG) schemes discretizing the entropy variables Barth ; Hughes . High-order ES schemes are not limited to these two options (formulations based on Summation-By-Parts operators Fisher ; Fried for instance are actively being developed), but the issues their formulation raises are no different.

4.1 Discontinuous Galerkin

In Hughes , Hughes et al. showed, under the assumption of exact numerical quadrature, that continuous finite element solutions to the compressible Navier-Stokes equations become consistent with the entropy equation when the entropy variables are discretized instead of conservative variables. Furthermore, they showed that with suitably defined dissipation operators, the discrete solution satisfies the Clausius-Duhem inequality.
For Discontinuous Galerkin discretizations of the Compressible Euler Equations, the same result follows if ES fluxes are used. Yet, the formulation of ES DG schemes by Barth Barth does not proceed along these lines. The entropy variables are discretized but entropy-stable fluxes are defined in a different way. This is discussed in C.
The fact that the entropy variables are undefined in the limit $\rho_{k}\rightarrow 0$ poses a daunting problem, unless the flow configuration is such that one can expect $\rho_{k}>0$ at all times. A way around this issue (other than not discretizing the entropy variables) has not been found by the authors. There might be other entropy functions for which the corresponding entropy variables are well-defined in this limit. However, in the context of the multicomponent Navier-Stokes equations (including viscous stresses, heat conduction, multicomponent diffusion), the a priori entropy stability resulting from discretizing the entropy variables might be lost Hughes ; Giovangigli_Mat .

4.2 TecNO schemes

TecNO schemes (Fjordholm et al. Fjordholm ) are high-order entropy-stable finite volume schemes that combine the high-order EC flux formulation of LeFloch et al. LeFloch , the stencil selection procedure of ENO schemes HartenENO ; ShuOsher and entropy-stable dissipation operators Tadmor ; Roe1 .
The dissipation term of a first order ES flux typically takes the form:

D[\mathbf{v}]=(RT)|\Lambda|(RT)^{T}[\mathbf{v}],

where $R$ is the matrix of right eigenvectors of $A$ and $T$ is a scaling matrix. The reconstruction used by TecNO schemes is motivated by the sign property of the ENO reconstruction. It was shown Fjordholm2 that for any vector $\mathbf{w}\in\mathbb{R}^{N}$ , the ENO reconstruction is such that:

\left<\mathbf{w}\right>=B[\mathbf{w}],\ B=diag([b_{0},\ \dots,\ b_{N-1}]),\ b_{i}\geq 0,

where $\left<\mathbf{w}\right>$ and $[\mathbf{w}]$ are the reconstructed and initial jumps, respectively. Let $\mathbf{w}$ be such that $[\mathbf{w}]=(RT)^{T}[\mathbf{v}]$ , then the dissipation operator $\tilde{D}[\mathbf{w}]=(RT)|\Lambda|\left<\mathbf{w}\right>$ is ES because:

(RT)|\Lambda|\left<\mathbf{w}\right>=(RT)|\Lambda|B[\mathbf{w}]=(RT)(|\Lambda|B)(RT)^{T}[\mathbf{v}]

$|\Lambda|B$ is a positive diagonal matrix, therefore the high-order dissipation operator $\tilde{D}[\mathbf{w}]$ is ES. The TecNO approach does require the knowledge of the scaled eigenvectors $(RT)$ . If the dissipation operator is only known in the form:

D[\mathbf{v}]=R|\Lambda|T^{2}R^{T}[\mathbf{v}],

then choosing $\mathbf{w}$ such that $[\mathbf{w}]=(T^{2}R^{T})[\mathbf{v}]$ for the ENO reconstruction will result in a high-order dissipation operator:

\tilde{D}[\mathbf{w}]=R|\Lambda|\left<\mathbf{w}\right>=R(|\Lambda|BT^{2})R^{T}[\mathbf{v}].

$|\Lambda|$ , $B$ and $T^{2}$ are all symmetric and at least positive semi-definite, however their product is no necessarily positive semi-definite. If the squared scaling matrix is diagonal then entropy stability is preserved. In the general case of a block-diagonal scaling matrix, entropy stability is no longer ensured. The same problem would arise if $\mathbf{w}$ was chosen such that $[\mathbf{w}]=R^{T}[\mathbf{v}]$ (the question would be whether $(|\Lambda|T^{2}B)$ is positive definite).
If the dissipation operator is expressed according to equation (32) with the pseudo-scaling matrix $\mathcal{T}$ given by equation (31), the vector of reconstructed variables $\mathbf{w}$ such that $=[\mathbf{w}]=(R\mathcal{T})^{T}[\mathbf{v}]$ will have as many components as the number of rows of $\mathcal{T}$ . For $N>2$ , the number of rows of $\mathcal{T}$ grows as $\mathcal{O}(N^{2})$ , so for a large enough $N$ , one might prefer working with $T$ instead of $\mathcal{T}$ and avoid having to reconstruct too many variables.
In section 3.3, we showed that the dissipation term expressed as $R|\Lambda|T^{2}R^{T}$ is well defined in the limit $\rho_{k}\rightarrow 0$ provided that $\rho_{k}^{*}=\rho_{k}^{ln}$ . This was possible because one could extract from $T^{2}$ a diagonal matrix $\mathcal{R}$ (see equation (33)) of partial densities. The TecNO algorithm requires the isolated evaluation of the scaled entropy variables defined by the jump relation $[\mathbf{w}]=(R\mathcal{T})^{T}[\mathbf{v}]$ or $[\mathbf{w}]=(RT)^{T}[\mathbf{v}]$ . The matrix $\mathcal{R}$ that was extracted from the squared scaling matrix, might not be extracted from $T$ or $\mathcal{T}$ without leaving $1/\rho_{k}^{*}$ terms behind. However, $\mathcal{R}^{1/2}$ can be extracted from the pseudo-scaling matrix $\mathcal{T}$ and since:

[\ln(\rho_{k})]=[2\ \ln(\sqrt{\rho_{k}})]=2\ \frac{[\sqrt{\rho_{k}}]}{(\sqrt{\rho_{k}})^{ln}},

it can be shown that $[\mathbf{w}]=(R\mathcal{T})^{T}[\mathbf{v}]$ is well-defined in the limit $\rho_{k}\rightarrow 0$ , provided that $\rho_{k}^{*}=((\sqrt{\rho_{k}})^{ln})^{2}$ . The decomposition:

T^{2}=\mathcal{R}^{1/2}\widetilde{T}^{2}(\mathcal{R}^{1/2})^{T},\ \widetilde{T}^{2}=\frac{1}{\gamma r}\begin{bmatrix}(\gamma-1)Y_{1}+\gamma r_{2}/r_{1}Y_{2}&-\sqrt{Y_{1}Y_{2}}&0&0\\ -\sqrt{Y_{1}Y_{2}}&(\gamma-1)Y_{2}+\gamma r_{1}/r_{2}Y_{1}&0&0\\ 0&0&\frac{1}{2}&0\\ 0&0&0&\frac{1}{2}\end{bmatrix},

(45)

suggests that a TecNO reconstruction based on the scaling matrix $T$ would still be defined in the limit with the same averaging. Note that this average is not compatible with the stationary contact preservation condition (35) we derived in section 3.3. That is because equation (35) requires $\rho_{k}^{*}=\rho_{k}^{ln}$ .

5 Numerical experiments

In this section, we present and discuss numerical results on 1D and 2D test problems involving interfaces and shocks. A 3D formulation of the scheme (entropy variables, EC flux, ES flux) is provided in B.
In the 1D problems, the first-order finite volume scheme with the ES flux in space and forward Euler in time with a CFL of 0.3 on three grids (100, 300, and 1000 cells, respectively) was applied. All the figures in the next three subsections use the same legend as figure 1-(a). A fourth-order TecNO scheme is used in the 2D problem (section 5.4).

5.1 Moving interface

The first test problem is the advection of a contact discontinuity (constant velocity and constant pressure) separating two different species. The initial conditions are given by:

\displaystyle\begin{cases}(\rho_{1},\ \rho_{2},\ u,\ p)=&(0.1,\ 0.0,\ 1.0,\ 1.0),\ 0\leq x\leq 0.5,\\ (\rho_{1},\ \rho_{2},\ u,\ p)=&(0.0,\ 1.0,\ 1.0,\ 1.0),\ 0.5<x\leq 1.\end{cases}

with $\gamma_{1}=1.4,\ \gamma_{2}=1.6$ , $c_{v1}=c_{v2}=1$ . The velocity and pressure profiles at $t=0.022s$ and $t=0.1s$ are shown in figures 1 and 2. These profiles show overshoots and undershoots that are typically observed with conservative schemes (see figure 4 in Abgrall & Karni Abgrall for instance). These anomalies are often termed oscillations in the literature Abgrall ; Karni ; Abgrall1 ; SBI_Marquina ; SBI_Johnsen . Upon closer examination, we can see three wave packets propagating at different speeds. The one moving to the left has the fastest propagation speed, roughly $-3$ . The remaining two are moving to the right with propagation speeds of roughly $1$ (the speed of the contact) and $2$ . From figure 3, which shows the acoustic eigenvalues $u\pm a$ , it clearly appears that that the first and third wave packets are acoustic.
Figure 4 shows the total entropy $\rho s=\rho_{1}s_{1}+\rho_{2}s_{2}$ profiles at $t=0.02s$ and $t=0.1s$ . The wave structure of the pressure and velocity anomalies is more apparent, and we can see that each wave is carrying an spurious increase in entropy. Figure 6 shows the evolution of the total entropy over time (the contributions from the boundaries were removed). This shows that the anomalies observed do not violate entropy stability. On the contrary, it appears that inappropriate production of entropy is the issue. Additionally, we find (A.2) that these anomalies are still present when the upwind dissipation operator is discarded and the Backward Euler scheme is used to ensure entropy stability (recall the discussion in section 3.4.1 and equation (42)). In the same vein, we find that these anomalies do not violate a minimum principle of the specific entropy Gouasmi4 (see figure 5(b)).
That these anomalies are not linked to entropy stability or the minimum entropy principle does not come as a surprise considering that they were already observed with the Godunov scheme Abgrall ; Karni , which is both entropy stable and satisfies a minimum entropy principle Tadmor2 .
It is important to understand that these anomalies are intrinsic to multicomponent flows. There are no such anomalies on moving contacts in the compressible Euler equations (A.1). Furthermore, these anomalies are produced by first-order schemes. They are therefore different than the oscillations typically observed with high-order schemes when discontinuous solutions are sought. Intrigued readers are strongly encouraged to read early studies of this problem Karni ; Abgrall ; Abgrall .
For this problem and the following, high-order TecNO schemes were found to produce non-physical values of density and pressure.

5.2 Two-species shock-tube problem

We simulate a shock-tube problem with two species. The initial conditions are given by:

\displaystyle\begin{cases}(\rho_{1},\ \rho_{2},\ u,\ p)=&(1,\ 0,\ 0,\ 1),\ 0\leq x\leq 0.5,\\ (\rho_{1},\ \rho_{2},\ u,\ p)=&(0,\ 0.125,\ 0,\ 0.1),\ 0.5<x\leq 1.0,\end{cases}

with $\gamma_{1}=1.4,\ \gamma_{2}=1.6$ and $c_{v1}=c_{v2}=1$ . The exact solution to the multicomponent shock-tube problem is almost the same as the exact solution in the single-component case. The composition of the gas changes across the contact discontinuity but does not change across rarefaction waves and shock waves Fezoui .
Figures 7(a)-(d) show the velocity, pressure, total density and specific heat ratio profiles at $t=0.2s$ . There is a good agreement with the exact solution, and we do not observe pressure and velocity oscillations around the moving contact which separates the two species.

5.3 Shock-interface interaction

We simulate a test problem from Quirk & Karni Karni which consists of a shock tube filled with air, where a shock wave moves to the right and eventually meets a stationary bubble of helium at pressure equilibrium. The initial conditions are given by:

\displaystyle\begin{cases}(\rho_{1},\ \rho_{2},\ u,\ p)=&\ (1.3765,\ 0,\ 0.3948,\ 1.57),\ 0\leq x\leq 0.25,\ \mbox{Post-shock, air},\\ (\rho_{1},\ \rho_{2},\ u,\ p)=&\ (1.,\ 0,\ 0.,\ 1.),\ 0.25\leq x\leq 0.4,\ \mbox{Pre-shock, air},\\ (\rho_{1},\ \rho_{2},\ u,\ p)=&\ (0.,\ 0.139,\ 0.,\ 1.),\ 0.4\leq x\leq 0.6,\ \mbox{Pre-shock, helium bubble},\\ (\rho_{1},\ \rho_{2},\ u,\ p)=&\ (1.,\ 0,\ 0.,\ 1.),\ 0.6\leq x\leq 1,\ \mbox{Post-shock, air}.\end{cases}

For air $c_{v1}=0.72,\ \gamma_{1}=1.4$ . For helium $c_{v2}=2.42,\ \gamma_{2}=1.67$ . In Quirk , this is problem is used to highlight the better behavior of Karni’s non-conservative scheme Karni over a conservative scheme using the Roe flux (see figure 2 in Quirk ). In a similar spirit, we compared our semi-discrete ES scheme with the Roe scheme. Figure 8 shows the pressure profile at $t=0.35s$ obtained with each scheme. As expected, the solution with Roe’s scheme is rife with oscillations unlike the solution with the present ES scheme which is free of these oscillations. The cause of this improvement is not entropy stability, but the property of preserving stationary contact discontinuities. Figure 9 shows the pressure profile before the right-moving shock a couple of instants before it meets the helium bubble. Roe’s scheme does not preserve stationary contacts and therefore produces pressure anomalies which eventually pollute the solution at $t=0.35$ (figure 8). This problem is a good illustration of the importance of treating interfaces properly in the simulation of multicomponent compressible flows. The results also suggest that the current ES scheme will not produce oscillations on shock-interface interaction problems if the interface is stationary. This encouraged us to consider the next problem.

5.4 Shock-bubble interaction

A test case that is commonly used in the development of numerical schemes for compressible multicomponent flows SBI_Kawai ; SBI_Marquina ; SBI_Johnsen ; Quirk is the interaction of a shock wave with a cylindrical gas inhomogeneity. This problem is a two-dimensional analog of the three-dimensional shock-induced mixing concept proposed by Marble et al. Marble in the context of supersonic scramjet design. This problem is also used in experimental and computational investigations of the Richtmyer-Meshkov instability Richtmyer ; Meshkov .
Validating the present ES scheme against experimental data is beyond the scope of the present work. In this section, we are essentially interested in the ability of the scheme to simulate the physics relevant to this classic problem. For this purpose, we tried to reproduce the results of Marquina & Mulet SBI_Marquina . The computational domain (ABCD) is shown in figure 10. A Mach $M_{S}=1.22$ shock wave, positioned at $x=275\ mm$ , moves to the left through quiescent air (species 1, $\gamma_{1}=1.4$ and $r_{1}=0.287\ 10^{3}\ J.kg.^{-1}.K^{-1}$ ) and eventually meets a cylindrical bubble, centered at $(x,y)=[225,\ 0]\ mm$ , filled with helium contaminated with 28% of air ( $\gamma_{2}=1.647,r_{2}=1.578\ 10^{3}\ J.kg^{-1}.K^{-1}$ ). The flow is assumed to be symmetric about the shock-tube axis (BC), therefore only the upper half of the physical domain is considered. Reflecting boundary conditions are applied on the top (AD) and bottom (BC) boundaries. The boundary conditions upstream (AB) and downstream (CD) the shock are not crucial in this problem Karni so we simply extrapolate the flow, as in SBI_Marquina .
Since the current scheme is unable to guarantee positive partial densities and pressure, the initial conditions from SBI_Marquina had to be modified. We set:

	$\displaystyle\mbox{Region I:}\ \ (\rho_{1},\ \rho_{2},\ u,\ v,\ p)=$	$\displaystyle\ (\delta\rho,\ 1.225(r_{1}/r_{2})-\delta\rho,\ 0.,\ 0.,\ 101325),$
	$\displaystyle\mbox{Region II:}\ \ (\rho_{1},\ \rho_{2},\ u,\ v,\ p)=$	$\displaystyle\ (1.225-\delta\rho,\ \delta\rho,\ 0.,\ 0.,\ 101325),$
	$\displaystyle\mbox{Region III:}\ \ (\rho_{1},\ \rho_{2},\ u,\ v,\ p)=$	$\displaystyle\ (1.6861-\delta\rho,\ \delta\rho,\ -113.5243,\ 0.,\ 159060).$

with $\delta\rho=0.03$ (units for density, velocity and pressure are $kg/m^{3}$ , $m/s$ and Pa, respectively). This setup differs from the original in two aspects. First, the composition of the gas in regions I, II and III will not be the same. Second, regions I and II are in pressure and temperature equilibrium in the original setup whereas in ours, temperature equilibrium is lost. These differences make quantitative comparisons, notably with the experimental data of Haas & Sturtevant SBI_Haas , difficult to carry out. However, we expect the physics to remain similar qualitatively (see Picone & Boris Picone who studied this problem using a single gas flow model).
For this simulation, we used a 4-th order TecNO scheme in space with a 4-th order explicit Runge-Kutta scheme in time. We used a $4000\times 400$ grid and set the CFL number to $0.3$ . Figure 11 shows pseudo-schlieren images of the density gradients at different times after the shock reached the bubble. These are in good qualitative agreement with those produced by Marquina & Mulet, figure 7 in SBI_Marquina (see also Quirk & Karni Quirk , figure 9). We refer to these two references for a detailed discussion of the physical mechanisms at work.
In figure 12-(a), we show an x–t diagram of the position of the key features of the shock-bubble interaction. These features are explained in figure 12-(b). The positions of these features are obtained by looking at inflection points of horizontal sections of the shading function $\phi$ used in figure 11. The upstream bubble interface is tracked on a section at a height $20$ mm from the axis. The incident shock is tracked on a section at $5$ mm from the top wall. The remaining features are tracked on a section along the symmetry line. The x-t diagram from Marquina & Mulet, figure 5 in SBI_Marquina , shows similar trends. The mean velocities of these features are calculated from their visually straight trajectories using linear regression, and displayed in Table 1.

	VS	VR	VT	Vui1	Vui2	Vdi	Vj
Gouasmi et al.	422	954	377	185	105	138	228
Marquina and Mulet SBI_Marquina	414	943	373	176	111	153	229

Table 1: Velocities in

m/s

of the features explained in figure 12. The time intervals in

\mu s

for computing each velocity are: VS [3.66, 62.64], VR [3.66, 52.81], VT [52.81, 141.26], Vui1 [3.66, 141.26], Vui2 [146.18, 254.29], Vdi [141.26, 254.29], Vj [146.18, 254.29].

6 Conclusions

We formulated entropy-stable schemes for the multicomponent compressible Euler equations. This effort built on the theoretical ground laid out by Chalot et al. Chalot and Giovangigli Giovangigli and followed a procedure pioneered by Tadmor Tadmor : we first derived a baseline EC flux and complemented it with an entropy-producing upwind dissipation operator Roe1 . We showed that in the limit $\rho_{k}\rightarrow 0$ , the EC flux is well-defined. This also holds for the upwind dissipation operator and the TecNO reconstruction provided that the averaged partial densities are well-chosen. Unfortunately, this does not prevent the scheme from producing negative values of partial densities and pressure. We also derived a condition on the averaged state of the dissipation matrix so that the ES scheme can exactly preserve stationary contact discontinuities.
It is a well-known issue that conservative schemes are subject to pressure oscillations in moving interface configurations. Numerical experiments showed that the ES scheme we constructed is no exception. We stress that these anomalies, which are not present in the single component case, violate neither entropy stability nor a minimum principle of the specific entropy Gouasmi4 . The remedies to the pressure oscillations problem typically consist in giving up on conservation of total energy Abgrall1 ; Karni ; Billet ; Abgrall , which can impair the ability of the scheme to properly capture shocks. A compromise between ensuring entropy stability and the proper treatment of moving interfaces could perhaps be achieved with the EC/ES schemes for non-conservative hyperbolic systems developed by Castro et al. Castro , even though non-conservative schemes have their own lot of issues Hou ; Abgrall_NC .
We remain careful to not make any peremptory statement about ES schemes. This work highlighted issues associated with an ES formulation that is more of a default choice than a well-grounded one. The choices of EC flux and dissipation operator are not unique and have to be explored in more depth. Ultimately, we hope that this effort will encourage further work on the proper way to manage entropy locally in a discrete field.

Acknowledgments

Ayoub Gouasmi and Karthik Duraisamy were funded by the AFOSR through grant number FA9550-16-1-030 (Tech. monitor: Fariba Fahroo). Ayoub Gouasmi would like to thank Laslo Diosady and Philip Roe for helpful conversations, and Eitan Tadmor for discussions on the minimum entropy principle Tadmor2 .
Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center.

References

(1) Friedrichs, K. O., and Lax, P. D. : Systems of Conservation Equations with a Convex Extension, P. Natl. Acad. Sci. U.S.A, 68 (8) pp. 1686-1688, 1971.
(2) Lax, P. : Shock Waves and Entropy. Contributions to Nonlinear Functional Analysis, pp. 603–634, 1971.
(3) Harten A. : On the symmetric form of systems of conservation laws with entropy, J. Comput. Phys., 49 (1) pp. 151-164 , 1983.
(4) Mock, M. S. : Systems of conservation laws of mixed type, J. Hyperbol. Differ. Eq., 70 (1) pp. 70-88, 1980.
(5) Chiodaroli, E., A counterexample to well-posedness of entropy solutions to the compressible euler system, J. Hyperbol. Differ. Eq., 11, pp. 493-519, 2014.
(6) de Lellis, D., and Szekelyhidi, L. : On admissibility criteria for weak solutions of the euler equations, Arch. Ration. Mech. An., 195, pp. 225-260, 2010.
(7) Elling, V., The carbuncle phenomenon is incurable, Acta Math. Sci., 29 (6), pp. 1647-1656, 2009.
(8) Hughes, T.J. R., Franca, L. P., and Mallet, M. : A new finite element formulation for computational fluid dynamics: I. Symmetric forms of the compressible Euler and Navier-Stokes equations and the second law of thermodynamics, Comput. Method. Appl. M., 54 (2) pp. 223-234, 1986.
(9) Reed, W. H., and Hill, T. R. : Triangular mesh methods for the neutron transport equation, Technical Report LA-UR-73-479, Los Alamos Scientific Laboratory, 1973.
(10) Cockburn B., Karniadakis G.E., and Shu C.W. : The Development of Discontinuous Galerkin Methods. In: Discontinuous Galerkin Methods. Lecture Notes in Computational Science and Engineering, vol 11. Springer, Berlin, Heidelberg, 2000.
(11) Shu, C.W., and Osher, S. : Efficient implementation of essentially non-oscillatory shock capturing schemes, J. Comput. Phys., 77 (2) pp. 439-471, 1988.
(12) Giovangigli, V. : Multicomponent flow modeling, Chapters 1 and 8, Birkhauser, Boston, 1999.
(13) Giovangigli, V., and Matuszewski, L. : Structure of Entropies in Dissipative Multicomponent Fluids, Kin. Rel. Models, 6 pp. 373-406, 2013.
(14) Chalot, F., Hughes, T.J.R. and Shakib, F. : Symmetrization of conservation laws with entropy for high-temperature hypersonic computations, Comput. Syst. Engrg., 1 (2-4) pp. 495-521, 1990.
(15) Tadmor, E. : The numerical viscosity of entropy stable schemes for systems of conservation laws. I, Math. Comput., 49 pp. 91-103, 1987.
(16) Tadmor, E. : The entropy dissipation by numerical viscosity in nonlinear conservative difference schemes, in: Nonlinear Hyperbolic Problems. Lecture Notes in Mathematics, 1270. Springer, Berlin, Heidelberg, 1987.
(17) Tadmor, E. : A minimum entropy principle in the gas dynamics equations, Appl. Numer. Math., 2 (3-5) pp. 211-219, 1986.
(18) Tadmor, E. : Entropy stability theory for difference approximations of nonlinear conservation laws and related time-dependent problems, Acta Numer., 12, pp. 451–512, 2003
(19) Zhong, W., and Tadmor, E. : Entropy Stable Approximations of Navier-Stokes equations with no artificial numerical viscosity, J. Hyperbol. Differ. Eq., 3 (3) pp. 529-559, 2006.
(20) Guermond, J.L., Popov, B. : Invariant Domain and First-Order Continuous Finite Element Approximation for hyperbolic systems, SIAM J. Numer. Anal., 54 (4), pp. 2466-2489, 2016.
(21) Guermond, J.L., Nazarov, M., Popov, B., Tomas, I. : Second-order invariant domain preserving approximation of the Euler equations using convex limiting, SIAM J. Sci. Comput., 40 (5), pp. 3211–3239, 2018.
(22) Guermond, J.L., and Popov, B. : Fast estimation from above of the maximum wave speed in the Riemann problem for the Euler equations, J. Comput. Phys., 328 pp. 908-926, 2016.
(23) Roe, P. L. : Approximate Riemann Solvers, Parameter vectors, and difference schemes, J. Comput. Phys., 43 (2) pp. 357-372, 1981.
(24) Roe, P.L. : Affordable, entropy consistent flux functions. In: Eleventh International Conference on Hyperbolic Problems: Theory, Numerics and Applications, 2006.
(25) Merriam, M. : An Entropy-Based Approach to Nonlinear Stability, NASA Technical Memorandum 101086, 1989.
(26) Godunov, S.K. : An interesting class of quasilinear systems, Dokl. Akad. Nauk. SSSR, 139 pp. 521-523, 1961.
(27) Barth, T.J. : Numerical Methods for Gasdynamic Systems on Unstructured Meshes. In: Kröner D., Ohlberger M., Rohde C. (eds) An Introduction to Recent Developments in Theory and Numerics for Conservation Laws. Lecture Notes in Computational Science and Engineering, vol 5. Springer, Berlin, Heidelberg, 1999.
(28) Barth, T.J. : Simplified Discontinuous Galerkin Methods for Systems of Conservation Laws with Convex Extension. In: Cockburn B., Karniadakis G.E., Shu CW. (eds) Discontinuous Galerkin Methods. Lecture Notes in Computational Science and Engineering, vol 11. Springer, Berlin, Heidelberg, 2000.
(29) Barth, T.J. : On the Role of Involutions in the Discontinuous Galerkin Discretization of Maxwell and Magnetohydrodynamic Systems, in: Arnold D.N., Bochev P.B., Lehoucq R.B., Nicolaides R.A., Shashkov M. (eds) Compatible Spatial Discretizations. The IMA Volumes in Mathematics and its Applications, vol 142. Springer, New York, 2006.
(30) Ismail, F., and Roe, P.L. : Affordable, entropy-consistent Euler flux functions II: Entropy production at shocks, J. Comput. Phys., 228 (15) pp. 5410-5436, 2009.
(31) Ismail, F. : Toward a reliable prediction of shocks in hypersonic flow: resolving carbuncles with entropy and vorticity control, PhD thesis, University of Michigan, 2006.
(32) Fjordholm, U.S., Mishra, S., and Tadmor, E. : Arbitrarily High-order Accurate Entropy Stable Essentially Nonoscillatory Schemes for Systems of Conservation Laws, SIAM J. Numer. Anal., 50 (2) pp. 544-573, 2012.
(33) Fjordholm, U.S., Mishra, S., and Tadmor, E. : ENO Reconstruction and ENO Interpolation Are Stable, Found. Comput. Math., 13 (2) pp. 139-159, 2012.
(34) LeFloch, P.G., Mercier, J. M., and Rhode, C. : Fully Discrete, Entropy Conservative Schemes of Arbitrary Order, SIAM J. Numer. Anal., 40 (5) pp. 1968-1992, 2002.
(35) Harten, A., Engquist, B., Osher, S., and Chakravarthy, S. : Uniformly high order essentially non-oscillatory schemes III, J. Comput. Phys., 71 (2) pp. 231-303, 1987.
(36) Osher, S. : Riemann Solvers, the Entropy Condition, and Difference, SIAM J. Numer. Anal., 21 (2) pp. 217–235, 1984.
(37) Abgrall, R. : Généralisation du solveur de Roe pour le calcul d’écoulements de mélanges de gaz parfaits à concentrations variables. La Recherche Aérospatiale, 6, pp 31–43, 1988.
(38) Billet, G., and Abgrall, R. : An adaptive shock-capturing algorithm for solving unsteady reactive flows, Comput. Fluids, 32 (10) 1473-1495, 2003.
(39) Fernandez, G., and Larrouturou, B. : Hyperbolic Schemes for Multi-Component Euler Equations. In: Ballmann J., Jeltsch R. (eds) Nonlinear Hyperbolic Equations — Theory, Computation Methods, and Applications. Notes on Numerical Fluid Mechanics, vol 24. Vieweg-Teubner Verlag, 1989.
(40) Larrouturou, B. : How to preserve the mass fractions positivity when computing compressible multi-component flows, J. Comput. Phys., 95 (1) 59-84, 2001.
(41) Larrouturou, B., and Fezoui, L. : On the equations of multi-component perfect of real gas inviscid flow, Nonlinear Hyperbolic Problems, pp.69-98, 1989.
(42) Habbal, A., Dervieux, A., Guillard, H., and Larrouturou, B. : Explicit calculation of reactive flows with an upwind finite element hydro-dynamical code. Technical Report 690, INRIA, 1987.
(43) Abgrall, R., and Karni, S. : Computations of Compressible Multifluids, J. Comput. Phys., 169 (2) pp. 594-623, 2000.
(44) Karni, S. : Multicomponent flow calculation by a consistent primitive algorithm. J. Comput. Phys., 112 (1) pp. 31-43, 1994.
(45) Abgrall, R. : How to prevent pressure oscillations in multicomponent flow calculations: a quasi-conservative approach. J. Comput. Phys., 125 (1) pp. 150-160, 1996.
(46) Richtmyer, R.D. : Taylor instability in shock acceleration of compressible fluids, Comm. Pure Appl. Math., 13 (2) pp. 297-319, 1960.
(47) Meshkov, E. E. : Instability of the interface of two gases accelerated by a shock wave, Fluid Dynam., 4 (5) pp. 101-104, 1969.
(48) Picone, J. M., and Boris, J. P. : Vorticity generation by shock propagation through bubbles in a gas, J. Fluid Mech., 189 pp. 23-51, 1988.
(49) Marble, F. E, Hendricks, G.J., and Zukoski, E. : Progress towards shock enhancement of supersonic combustion processes, 23rd Joint Propulsion Conference, 1987.
(50) Karni, S., and Quirk, J.J. : On the dynamics of a shock-bubble interaction, J. Fluid Mech., 318 pp. 129-163, 1996.
(51) Quirk, J.J. : A contribution to the great Riemann solver debate, Int. J. Numer. Fl., 18 pp. 555-574, 1994.
(52) Kawai, S., and Terashima, H. : A high-resolution scheme for compressible multicomponent flows with shock waves, Int. J. Numer. Fl., 66 (10) pp. 1207-1225, 2011.
(53) Haas, J.F., and Sturtevant, B. : Interactions of weak shock waves with cylindrical and spherical gas inhomogeneities, J. Fluid Mech., 181 pp. 41-76, 1987.
(54) Marquina, A., and Mulet, P. : A flux-split algorithm applied to conservative models for multicomponent compressible flows, J. Comput. Phys., 185 (1) pp. 120-138, 2003.
(55) Johnsen, E., and Colonius, T. : Implementation of WENO schemes in compressible multicomponent flow problems, J. Comput. Phys., 219 (2) pp. 715-732, 2006.
(56) Jameson, A. : Formulation of kinetic energy preserving conservative schemes for gas dynamics and direct numerical simulation of one-dimensional viscous compressible flow in a shock tube using entropy and kinetic energy preserving schemes, J. Sci. Comput., 34(2) pp. 188–208, 2008.
(57) Chandrasekhar, P. : Kinetic energy preserving and entropy stable finite volume schemes for compressible Euler and Navier–Stokes equations, Commun. Comput. Phys., 14 pp. 1252–1286, 2013.
(58) Subbareddy, P., and Candler, G.V. : A fully discrete, kinetic energy consistent finite-volume scheme for compressible flows, J. Comput. Phys., 228 (5) pp. 1347-1364, 2009.
(59) Derigs, D., Winters, A. R., Gassner, G. J., and Walch, S. : A novel high-order, entropy stable, 3D AMR MHD solver with guaranteed positive pressure, J. Comput. Phys., 317 pp. 223-256, 2016.
(60) Ma, P.C., Lv, Y., and Ihme, M. : An entropy-stable hybrid scheme for simulations of transcritical real-fluid flows, J. Comput. Phys., 340 pp. 330-357, 2017.
(61) Castro, M.J., Fjordholm, U.S., Mishra, S. and Parés, C. : Entropy conservative and entropy stable schemes for nonconservative hyperbolic systems, SIAM J. Numer. Anal., 51(3) pp. 1371-1391, 2012.
(62) Gouasmi, A., Murman, S.M., and Duraisamy, K. : Entropy conservative schemes and the receding flow problem, J. Sci. Comput., 78 (2) pp. 971-994, 2018.
(63) Gouasmi, A., Duraisamy, K. and Murman, S.M. : On entropy stable temporal fluxes, arXiv:1807.03483, 2018.
(64) Gouasmi, A., Duraisamy, K. and Murman, S.M., and Tadmor, E. : A minimum entropy principle in the compressible multicomponent Euler equations, ESAIM-Math. Model. Num., accepted, 2019.
(65) Murman, S.M., and Frontin, C. : Analysis of numerical dissipation in entropy-stable schemes for turbulent flows, Center for Turbulence Research, Proceedings of the Summer Program 2018.
(66) Murman, S.M., Diosady, L.T., Garai, A., and Ceze, M. : A Space-Time Discontinuous-Galerkin Approach for Separated Flows, 54th AIAA Aerospace Sciences Meeting, 2016.
(67) Diosady, L. T., and Murman, S. M. : Higher-Order Methods for Compressible Turbulent Flows Using Entropy Variables, 53rd AIAA Aerospace Sciences Meeting, 2015.
(68) Carton de Wiart, C., Diosady, L.T., Garai, A., Burgess, N., Blonigan, P., and Murman, S.M. : Design of a modular monolithic implicit solver for multi-physics applications, AIAA SciTech Forum, 2018.
(69) Pazner, W., and Persson, P-O : Analysis and Entropy Stability of the Line Based Discontinuous Galerkin Method, J. Sci. Comput., 80 (1) pp. 376-402, 2019.
(70) Fernandez, P., Nguyen, N-C, and Peraire, J. : Entropy-stable hybridized discontinuous Galerkin methods for the compressible Euler and Navier-Stokes equations, arXiv:1808.05066, 2018.
(71) Fisher, T.C. and Carpenter, M.H. : High-order entropy stable finite difference schemes for nonlinear conservation laws: Finite domains, J. Comput. Phys., 252(1) pp 518-557, 2013.
(72) Friedrich, L., Schnücke, G., Winters, A.R., Del Rey Fernández, D.C., Gassner, G. J., and Carpenter, M.H. : Entropy Stable Space–Time Discontinuous Galerkin Schemes with Summation-by-Parts Property for Hyperbolic Conservation Laws, J. Sci. Comput., 2019.
(73) Ray, D. : Entropy-stable finite difference and finite-volume schemes for compressible flows, PhD thesis, TIFR, 2017.
(74) Hou, T.Y., and Le Floch : Why nonconservative schemes converge to wrong solutions: error analysis, Math. Comput., 62 (206), pp. 497-530, 1994.
(75) Abgrall, R., and Karni, S. : A comment on the computation of non-conservative products, J. Comput. Phys., 229(8), pp. 2759-2763, 2010.
(76) Zakerzadeh, H., and, Fjordholm, U.S. : High-order accurate, fully discrete entropy stable schemes for scalar conservation laws, IMA J. Numer. Anal., 36 (2), pp. 633–654, 2016.
(77) Slotnick J, Khodadoust A, Alonso J, Darmofal D, Gropp W, et al. : CFD Vision 2030 study: a path to revolutionary computational aerosciences. NASA Tech. Rep. CR-2014-218178, Langley Res. Cent., Hampton, VA.

Appendix A Additional results - Moving interface problem

A.1 Single component case

Here we present numerical results for a moving contact discontinuity in the compressible Euler equations. The initial conditions are given by:

\displaystyle\begin{cases}(\rho,\ u,\ p)=&(0.1,\ 1.,\ 1.),\ -1.0\leq x\leq 0.0,\\ (\rho,\ u,\ p)=&(1.,\ 1.,\ 1.),\ 0.0<x\leq 1.0,\end{cases}

with $\gamma=1.4$ and $c_{v}=1$ . The anomalies observed in the multicomponent case are not present. The velocity and pressure remain constant at all times.

A.2 Conserving entropy in space, producing entropy in time

The implicit scheme did not converge in the original setup (most likely due to the generation of negative densities/pressure). We therefore considered a different setup, with non-zero partial densities, for which the pressure oscillations problem is still present.

\displaystyle\begin{cases}(\rho_{1},\ \rho_{2},\ u,\ p)=&(0.3,\ 0.15,\ 1.,\ 1.),\ 0\leq x\leq 0.5,\\ (\rho_{1},\ \rho_{2},\ u,\ p)=&(0.15,\ 1.,\ 1.,\ 1.),\ 0.5<x\leq 1.0,\end{cases}

with $\gamma_{1}=1.4,\ \gamma_{2}=1.6$ and $c_{v1}=c_{v2}=1$ . A grid of 200 cells is used. Figure 14 shows the pressure profiles obtained with an EC flux and an ES flux in space, respectively, for two different CFL numbers. In the first case, the entropy production of the scheme only comes from the stabilization of the Backward Euler time scheme, which grows with $\Delta t$ . The high frequency oscillations observed are typically observed when too little dissipation is added to EC schemes Gouasmi1 ; Zhong . In the event that the scheme is EC at the fully discrete level (see Gouasmi1 for an example), these oscillations will be present but will not increase in magnitude over time. This is not desirable.

The main conclusion we draw is that the pressure anomalies remain present when the upwind dissipation typically used in space is replaced by the dissipation of Backward Euler. The magnitude of the pressure anomalies is higher when both the interface flux and time scheme are ES.

Appendix B Three-dimensional version

The three-dimensional multicomponent Euler equations are given by:

\frac{\partial\mathbf{u}}{\partial t}+\frac{\partial\mathbf{f_{1}}}{\partial x_{1}}+\frac{\partial\mathbf{f_{2}}}{\partial x_{2}}+\frac{\partial\mathbf{f_{3}}}{\partial x_{3}}=0.

(46)

The state vector $\mathbf{u}$ and flux vectors $\mathbf{f_{1}},\ \mathbf{f_{2}}$ and $\mathbf{f_{3}}$ are defined by:

	$\displaystyle\mathbf{u}=$	$\displaystyle\ \begin{bmatrix}\rho_{1}&\ldots&\rho_{N}&\rho u&\rho v&\rho w&\rho e^{t}\end{bmatrix}^{T},$
	$\displaystyle\mathbf{f_{1}}=$	$\displaystyle\ \begin{bmatrix}\rho_{1}u&\ldots&\rho_{N}u&\rho u^{2}+p&\rho uv&\rho uw&(\rho e^{t}+p)u\end{bmatrix}^{T},$
	$\displaystyle\mathbf{f_{2}}=$	$\displaystyle\ \begin{bmatrix}\rho_{1}v&\ldots&\rho_{N}v&\rho uv&\rho v^{2}+p&\rho vw&(\rho e^{t}+p)v\end{bmatrix}^{T},$
	$\displaystyle\mathbf{f_{3}}=$	$\displaystyle\ \begin{bmatrix}\rho_{1}w&\ldots&\rho_{N}v&\rho uw&\rho vw&\rho w^{2}+p&(\rho e^{t}+p)w\end{bmatrix}^{T}.$

The conservation equation for entropy writes:

\frac{\partial U}{\partial t}+\frac{\partial F_{1}}{\partial x_{1}}+\frac{\partial F_{2}}{\partial x_{2}}+\frac{\partial F_{3}}{\partial x_{3}}=0,\ U=-\rho s,\ F_{1}=-u\rho s,\ F_{2}=-v\rho s,F_{3}=-w\rho s.

(47)

The vector of entropy variables is:

\mathbf{v}=\frac{1}{T}\begin{bmatrix}g_{1}-\frac{1}{2}(u^{2}+v^{2}+w^{2})&\ldots&g_{N}-\frac{1}{2}(u^{2}+v^{2}+w^{2})&u&v&w&-1\end{bmatrix}.

(48)

The potential flux in space $\mathcal{U}$ is unchanged. There are now three spatial potential functions to work with. They are given by:

\mathcal{F}_{1}=\mathbf{v}\cdot\mathbf{f_{1}}-F_{1}=\sum_{k=1}^{N}\frac{R}{m_{k}}\rho_{k}u,\ \mathcal{F}_{2}=\mathbf{v}\cdot\mathbf{f_{2}}-F_{2}=\sum_{k=1}^{N}\frac{R}{m_{k}}\rho_{k}v,\ \mathcal{F}_{3}=\mathbf{v}\cdot\mathbf{f_{3}}-F_{3}=\sum_{k=1}^{N}\frac{R}{m_{k}}\rho_{k}w.

(49)

It can be easily shown that an entropy-conservative flux across an interface of normal $\mathbf{n}=(n_{1},n_{2},n_{3})$ , denoted $\mathbf{f^{*}}$ must satisfy:

[\mathbf{v}]\cdot\mathbf{f^{*}}=[n_{1}\mathcal{F}_{1}+n_{2}\mathcal{F}_{2}+n_{3}\mathcal{F}_{3}].

Using the same method as in the 1D case, one obtains $\mathbf{f^{*}}=[f_{1,1}\ \dots\ f_{1,N}\ f_{2}\ f_{3}\ f_{4}\ f_{5}]$ with:

$\displaystyle f_{1,k}=$	$\displaystyle\ \rho_{k}^{ln}\overline{u_{n}},$
$\displaystyle f_{2}=$	$\displaystyle\ \frac{n_{1}}{\overline{1/T}}\bigg{(}\sum_{k=1}^{N}r_{k}\overline{\rho_{k}}\bigg{)}+\overline{u}\sum_{k=1}^{N}f_{1,k},$
$\displaystyle f_{3}=$	$\displaystyle\ \frac{n_{2}}{\overline{1/T}}\bigg{(}\sum_{k=1}^{N}r_{k}\overline{\rho_{k}}\bigg{)}+\overline{v}\sum_{k=1}^{N}f_{1,k},$	(50)
$\displaystyle f_{4}=$	$\displaystyle\ \frac{n_{3}}{\overline{1/T}}\bigg{(}\sum_{k=1}^{N}r_{k}\overline{\rho_{k}}\bigg{)}+\overline{w}\sum_{k=1}^{N}f_{1,k},$
$\displaystyle f_{5}=$	$\displaystyle\ \sum_{k=1}^{N}(e_{0k}+c_{vk}\frac{1}{(1/T)^{ln}}-\frac{1}{2}\overline{u^{2}+v^{2}+w^{2}})f_{1,k}+\overline{u}f_{2}+\overline{v}f_{3}+\overline{w}f_{4}.$

$u_{n}=n_{1}u+n_{2}v+n_{3}w$ is the velocity normal to the interface. The temporal Jacobian Giovangigli is given by:

	$\displaystyle H=\begin{bmatrix}\frac{\rho_{1}}{r_{1}}&&0&\frac{\rho_{1}}{r_{1}}u&\frac{\rho_{1}}{r_{1}}v&\frac{\rho_{1}}{r_{1}}w&\frac{\rho_{1}}{r_{1}}(e_{1}^{t})\\ &\ddots&&\vdots&\vdots&\vdots&\vdots\\ &&\frac{\rho_{N}}{r_{N}}&\frac{\rho_{N}}{r_{N}}u&\frac{\rho_{N}}{r_{N}}v&\frac{\rho_{N}}{r_{N}}w&\frac{\rho_{N}}{r_{N}}(e_{N}^{t})\\ &&&\rho T+u^{2}S_{1}&uvS_{1}&uwS_{1}&u(\rho T+S_{2})\\ &&&&\rho T+v^{2}S_{1}&vwS_{1}&v(\rho T+S_{2})\\ &&&&&\rho T+w^{2}S_{1}&w(\rho T+S_{2})\\ sym&&&&&&\rho T(2k+c_{v}T)+S_{3}.\end{bmatrix},$
	$\displaystyle S_{1}=\sum_{k}\frac{\rho_{k}}{r_{k}},\ S_{2}=\sum_{k}\frac{\rho_{k}}{r_{k}}(e_{k}^{t}),\ S_{3}=\sum_{k}\frac{\rho_{k}}{r_{k}}(e_{k}^{t})^{2}$

Next is the scaling matrix. Let $A_{n}$ be the flux Jacobian in the normal direction:

A_{n}=n_{1}\frac{\partial\mathbf{f_{1}}}{\partial\mathbf{u}}+n_{2}\frac{\partial\mathbf{f_{2}}}{\partial\mathbf{u}}+n_{3}\frac{\partial\mathbf{f_{3}}}{\partial\mathbf{u}}.

A general expression for the eigenvector matrix $R$ such that $A_{n}=R\Lambda R^{-1}$ is the following:

	$\displaystyle R=\begin{bmatrix}1&&&0&0&Y_{1}&Y_{1}\\ &\ddots&&\vdots&\vdots&\vdots&\vdots\\ &&1&0&0&Y_{N}&Y_{N}\\ u&\ldots&u&&&u+an_{1}&u-an_{1}\\ v&\ldots&v&\mathbf{r_{I}}&\mathbf{r_{II}}&v+an_{2}&v-an_{2}\\ w&\ldots&w&&&w+an_{3}&w-an_{3}\\ k-\frac{d_{1}}{\gamma-1}&\ldots&k-\frac{d_{N}}{\gamma-1}&&&h^{t}+u_{n}a&h^{t}-u_{n}a\end{bmatrix},$
	$\displaystyle\Lambda=diag(\begin{bmatrix}u_{n}&\ldots&u_{n}&u_{n}&u_{n}+a&u_{n}-a\end{bmatrix}),$

where $\mathbf{r_{I}},\mathbf{r_{II}}\in\mathbb{R}^{4\times 1}$ are such that:

\mathbf{r_{I}},\mathbf{r_{II}}\in span\bigg{\{}\mathbf{e_{1}}=\begin{bmatrix}0\\ -an_{3}\\ an_{2}\\ -a(n_{3}v-n_{2}w)\end{bmatrix},\ \mathbf{e_{2}}=\begin{bmatrix}-an_{3}\\ 0\\ an_{1}\\ -a(n_{3}u-n_{1}w)\end{bmatrix},\ \mathbf{e_{3}}=\begin{bmatrix}-an_{2}\\ an_{1}\\ 0\\ -a(n_{2}u-n_{1}v)\end{bmatrix}\bigg{\}},

and $rank(\mathbf{r_{I}},\ \mathbf{r_{II}})=2$ . If $n_{1}\neq 0$ , take $(\mathbf{r_{I}},\mathbf{r_{II}})=(\mathbf{e_{3}},\mathbf{e_{2}})$ . The squared scaling matrix is given by:

T^{2}=R^{-1}HR^{-T}=\frac{\rho}{\gamma r}diag(\begin{bmatrix}T^{2Y}&T^{2N}&1/2&1/2\end{bmatrix}),

(51)

where $T^{Y}=\mathcal{T}^{Y}(\mathcal{T}^{Y})^{T}\in\mathbb{R}^{N\times N}$ is the same as in the 1D setting, and $T^{2N}\in\mathbb{R}^{2}$ is given by:

T^{2N}=\frac{1}{n_{1}^{2}}\begin{bmatrix}n_{1}^{2}+n_{3}^{2}&-n_{2}n_{3}\\ -n_{2}n_{3}&n_{1}^{2}+n_{2}^{2}\end{bmatrix}.

If $n_{2}^{2}+n_{3}^{2}=0$ , then $T^{2N}=I_{2\times 2}$ , otherwise $T^{2N}=T^{N}(T^{N})^{T}$ with:

T^{N}=\frac{1}{n_{1}\sqrt{n_{2}^{2}+n_{3}^{2}}}\begin{bmatrix}-n_{3}&n_{1}n_{2}\\ n_{2}&n_{1}n_{3}\end{bmatrix}.

If $n_{1}=0$ , then if $n_{2}\neq 0$ , take $(\mathbf{r_{I}},\mathbf{r_{II}})=(\mathbf{e_{3}},-\mathbf{e_{1}})$ . equation (51) holds with:

T^{2N}=\frac{1}{n_{2}^{2}}\begin{bmatrix}n_{2}^{2}+n_{3}^{2}&-n_{1}n_{3}\\ -n_{1}n_{3}&n_{2}^{2}+n_{1}^{2}\end{bmatrix}.

Since $n_{1}=0$ , this simplifies to $T^{2N}=diag([1+(n_{3}/n_{2})^{2}\ 1])$ and $T^{N}=diag([\sqrt{1+(n_{3}/n_{2})^{2}}\ 1])$ . If $n_{1}=0$ and $n_{2}=0$ , then $n_{3}\neq 0$ . In this final case, take $(\mathbf{r_{I}},\mathbf{r_{II}})=(\mathbf{e_{2}},\mathbf{e_{1}})$ . This leads to $T^{2N}=n_{3}^{2}diag([1\ 1])$ and $T^{N}=n_{3}diag([1\ 1])$ .
In each case, we can write $T^{2}=\mathcal{T}\mathcal{T}^{T}$ with:

\mathcal{T}=\sqrt{\frac{\rho}{\gamma r}}diag(\begin{bmatrix}\mathcal{T}^{Y}&T^{N}&1/\sqrt{2}&1/\sqrt{2}\end{bmatrix}).

The entropy-stable flux $\mathbf{f^{*}}$ writes:

\mathbf{f^{*}}=\mathbf{f_{EC}}-\frac{1}{2}(R\mathcal{T})|\Lambda|(R\mathcal{T})^{T}[\mathbf{v}].

Appendix C Some perspective on the fluxes of Tadmor and Barth

In Barth Barth , the entropy variables are discretized but the entropy stability of an inviscid interface flux $\mathbf{f^{*}}$ is defined by the condition that $\mathbf{f^{*}}$ produces more entropy than the so-called Symmetric Mean-Value (SMV) flux he proposed. This flux actually has two distinct expressions Barth ; Barth2 . In chronological order, they are given by:

	$\displaystyle\mathbf{f_{SMV,I}}=$	$\displaystyle\ \frac{1}{2}(\mathbf{f}(\mathbf{v}_{L})+\mathbf{f}(\mathbf{v}_{R}))-\frac{1}{2}\bigg{(}\int_{0}^{1}(1-\theta)(\|A(\mathbf{\overline{v}}(\theta))\|_{H}\ +\ \|A(\mathbf{\overline{\overline{v}}}(\theta))\|_{H})d\theta\bigg{)}[\mathbf{v}],$
	$\displaystyle\mathbf{f_{SMV,II}}=$	$\displaystyle\ \frac{1}{2}(\mathbf{f}(\mathbf{v}_{L})+\mathbf{f}(\mathbf{v}_{R}))-\frac{1}{2}\bigg{(}\int_{0}^{1}\|A(\mathbf{\overline{v}}(\theta))\|_{H}d\theta\bigg{)}[\mathbf{v}],$
	$\displaystyle\mathbf{\overline{v}}(\theta)=$	$\displaystyle\ \mathbf{v}_{L}+\theta[\mathbf{v}],\ \mathbf{\overline{\overline{v}}}(\theta)=\mathbf{v}_{R}-\theta[\mathbf{v}],\ \|A\|_{H}=\|A\|H=R\|\Lambda\|R^{T}.$

Barth Barth ; Barth2 defines an entropy-stable flux by the condition that it produces more entropy than the SMV flux:

[\mathbf{v}]\cdot\mathbf{f^{*}}\leq[\mathbf{v}]\cdot\mathbf{f_{SMV}}.

The SMV flux is entropy stable, therefore it does not define a genuine threshold for entropy stability. Tadmor Tadmor defines an entropy-stable flux by comparison with an entropy-conservative flux $\mathbf{f_{EC}}$ he derived. This condition writes:

[\mathbf{v}]\cdot\mathbf{f^{*}}\leq[\mathbf{v}]\cdot\mathbf{f_{EC}}=[\mathcal{F}],\ \mathbf{f_{EC}}=\int_{0}^{1}\mathbf{f}(\mathbf{\overline{v}}(\theta))d\theta.

Using integration by parts, this flux can be cast in viscosity form Tadmor :

	$\displaystyle\mathbf{f_{EC}}=$	$\displaystyle\ (\theta-\frac{1}{2})\mathbf{f}(\mathbf{\overline{v}}(\theta))\bigg{\|}_{\theta=0}^{\theta=1}-\int_{0}^{1}(\theta-\frac{1}{2})\frac{d}{d\theta}\mathbf{f}(\mathbf{\overline{v}}(\theta))d\theta$
	$\displaystyle=$	$\displaystyle\ \frac{1}{2}(\mathbf{f}(\mathbf{v}_{L})+\mathbf{f}(\mathbf{v}_{R}))-\frac{1}{2}\bigg{(}\int_{0}^{1}(2\theta-1)A(\mathbf{\overline{v}}(\theta))H(\mathbf{\overline{v}}(\theta))d\theta\bigg{)}[\mathbf{v}].$

Let $Q_{SMV}$ and $Q_{EC}$ be the matrices such that:

\mathbf{f_{EC}}=\frac{1}{2}(\mathbf{f}(\mathbf{v}_{L})+\mathbf{f}(\mathbf{v}_{R}))-\frac{1}{2}Q_{EC}[\mathbf{v}],\ \mathbf{f_{SMV}}=\frac{1}{2}(\mathbf{f}(\mathbf{v}_{L})+\mathbf{f}(\mathbf{v}_{R}))-\frac{1}{2}Q_{SMV}[\mathbf{v}].

Using:

[\mathbf{v}]\cdot AH[\mathbf{v}]=[\mathbf{v}]\cdot R\Lambda R^{T}[\mathbf{v}]\leq[\mathbf{v}]\cdot R|\Lambda|R^{T}[\mathbf{v}]=[\mathbf{v}]\cdot|A|_{H}[\mathbf{v}],

it can be easily shown that:

[\mathbf{v}]\cdot Q_{SMV,II}[\mathbf{v}]\geq[\mathbf{v}]\cdot Q_{EC}[\mathbf{v}]\ \rightarrow[\mathbf{v}]\cdot\mathbf{f_{SMV,II}}\leq[\mathbf{v}]\cdot\mathbf{f_{EC}}.

Note that a flux similar to $\mathbf{f_{SMV,II}}$ was hinted at by Tadmor (Example 3.2 in Tadmor1 ). Noting that $\mathbf{\overline{v}}(\theta)=\mathbf{\overline{\overline{v}}}(1-\theta)$ , and applying the variable change $\theta\leftarrow 1-\theta$ to half of $Q_{EC}$ leads to another expression for $Q_{EC}$ :

\displaystyle Q_{EC}=\int_{0}^{1}(\theta-\frac{1}{2})A(\mathbf{\overline{v}}(\theta))H(\mathbf{\overline{v}}(\theta))d\theta+\int_{0}^{1}(\frac{1}{2}-\theta)A(\mathbf{\overline{\overline{v}}}(\theta))H(\mathbf{\overline{\overline{v}}}(\theta))d\theta,

which resembles $Q_{SMV,I}$ . From here, it is not clear whether

[\mathbf{v}]\cdot Q_{EC}[\mathbf{v}]\leq[\mathbf{v}]\cdot Q_{SMV,I}[\mathbf{v}],

which would imply:

[\mathbf{v}]\cdot\mathbf{f_{SMV,I}}\leq[\mathbf{v}]\cdot\mathbf{f_{EC}}.

In a later publication BarthMHD , Barth defines an EC/ES condition:

[\mathbf{v}]\cdot(\mathbf{f^{*}}-\mathbf{f}(\mathbf{v}(\theta)))\leq 0,\ \forall\theta\in[0\ 1],

(52)

which is slightly different from Tadmor’s:

[\mathbf{v}]\cdot\mathbf{f^{*}}-[\mathbf{\mathcal{F}}]\leq 0.

(53)

However using the definition of the potential function $\mathcal{F}$ :

\frac{\partial\mathcal{F}}{\partial\mathbf{v}}^{T}=\mathbf{f}

and a corollary of the fundamental theorem of calculus:

[\mathbf{\mathcal{F}}]=[\mathbf{v}]\cdot\bigg{(}\int_{0}^{1}\frac{\partial\mathcal{F}}{\partial\mathbf{v}}^{T}(\mathbf{\overline{v}}(\theta))d\theta\bigg{)}.

Tadmor’s condition (53) can be rewritten as:

\int_{0}^{1}[\mathbf{v}]\cdot(\mathbf{f}^{*}-\mathbf{f}(\mathbf{\overline{v}}(\theta)))d\theta\leq 0,

(54)

which is the requirement that some path-integral is negative. Barth’s condition (52) is the more stringent requirement that the integrand in (54) is negative.