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

TL;DR

This paper develops entropy-stable numerical schemes for multicomponent compressible Euler equations, addressing challenges at zero partial densities, and demonstrates their effectiveness through various interface and shock interaction simulations.

Contribution

It introduces entropy-stable schemes specifically designed for multicomponent flows, including flux formulations that handle zero partial densities and high-order reconstructions.

Findings

Exact preservation of stationary interfaces.

Ability to reproduce physical shock-bubble interactions.

Handling of pressure anomalies on moving interfaces.

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 pressure anomalies typically observed with conservative schemes [Karni, \textit{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].

Figures

Tables

Table 1. Table 1: Velocities in m / s 𝑚 𝑠 m/s of the features explained in figure 12 . The time intervals in μ ​ s 𝜇 𝑠 \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].

	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

Equations

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

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

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

$$\mathbf{u}:={\left[\begin{array}{ccccc}{\rho }_1& \dots & {\rho }_N& \rho u& \rho e^t\end{array}\right]}^T,\mathbf{f}:={\left[\begin{array}{ccccc}{\rho }_{1}u& \dots & {\rho }_{N}u& \rho u^2+p& (\rho e^t+p)u\end{array}\right]}^T.$$

$$p:=\sum _{k=1}^N{\rho }_k{r}_{k}T,r_k=\frac{R}{m_k},$$

$$\rho e=\sum _k{\rho }_k{e}_k,e_k:=e_{0k}+c_{vk}T.$$

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

$$\rho s:=\sum _{k=1}^N{\rho }_k{s}_k,s_k:=c_{vk}\ln (T)-r_k\ln ({\rho }_k),$$

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

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

$$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,$$

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

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

$$\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},$$

$$Td(\rho s)=d(\rho e)-\sum _{k=1}^N{g}_{k}d{\rho }_k,$$

$$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}.$$

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

$$\frac{\partial U}{\partial Z}=\frac{1}{T}\left[\begin{array}{ccccc}(g_1-e_1)& \dots & (g_N-e_N)& 0& -\rho c_v\end{array}\right].$$

$$\frac{\partial \mathbf{u}}{\partial Z}=\left[\begin{array}{ccccc}1& & 0& 0& 0\\ & \ddots & & ⋮& ⋮\\ 0& & 1& 0& 0\\ u& \dots & u& \rho & 0\\ e_1+\frac{1}{2}{u}^2& \dots & e_N+\frac{1}{2}{u}^2& \rho u& \rho c_v\end{array}\right].$$

$$\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}$$

$$\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}.$$

$$\mathcal{U}=\sum _{k=1}^N{r}_k{\rho }_k,\mathcal{F}=\sum _{k=1}^N{r}_k{\rho }_{k}u.$$

$$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}}$$

$$\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{)}}.$$

$$\mathbf{v}=r{\left[\begin{array}{ccc}\frac{\gamma -\bar{s}}{\gamma -1}-\frac{\rho u^2}{2p}& \frac{\rho u}{p}& -\frac{\rho }{p}\end{array}\right]}^T,\bar{s}=\ln p-\gamma \ln \rho -\ln r.$$

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

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

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

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

$$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}).$$

Full text

††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:

[TABLE]

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$:

[TABLE]

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:

[TABLE]

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:

[TABLE]

$\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:

[TABLE]

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:

[TABLE]

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:

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

The potential functions satisfy the relationships:

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

This gives:

[TABLE]

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

[TABLE]

The inverse of this matrix is given by:

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

The partial densities are inferred from the specific entropies:

[TABLE]

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. 2. 2.

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 :

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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:

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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$:

[TABLE]

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:

[TABLE]

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

[TABLE]

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:

[TABLE]

Note the identities:

[TABLE]

The jump in the potential function can be rewritten as:

[TABLE]

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

[TABLE]

The corresponding jumps then write:

[TABLE]

The remaining jumps are given by:

[TABLE]

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:

[TABLE]

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

[TABLE]

The solution under these assumptions is:

[TABLE]

Therefore we obtain:

Theorem 3.3**.**

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

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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:

[TABLE]

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:

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

At this point, the dissipation operator writes:

[TABLE]

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:

[TABLE]

$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:

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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:

[TABLE]

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

[TABLE]

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:

[TABLE]

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:

[TABLE]

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:

[TABLE]

Proof.

Since $u=0$, we have:

[TABLE]

therefore:

[TABLE]

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

[TABLE]

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

[TABLE]

This is equivalent to:

[TABLE]

Using:

[TABLE]

equation (36) can then be rewritten as:

[TABLE]

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

[TABLE]

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

[TABLE]

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):

[TABLE]

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

[TABLE]

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)

[TABLE]

implies:

[TABLE]

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:

[TABLE]

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:

[TABLE]

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

[TABLE]

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$:

[TABLE]

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:

[TABLE]

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:

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

The corresponding potential functions are unchanged because:

[TABLE]

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}}$:

[TABLE]

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

[TABLE]

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

[TABLE]

$\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:

[TABLE]

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:

[TABLE]

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:

[TABLE]

$|\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:

[TABLE]

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:

[TABLE]

$|\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:

[TABLE]

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:

[TABLE]

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:

[TABLE]

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:

[TABLE]

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:

[TABLE]

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:

[TABLE]

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.

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.