A closed-loop green supply chain with retailers' competition and product recycling in the green environment under the cap-and-trade policy

1.   Introduction

Mathematical thermodynamics is an essential tool for studying multicomponent fluids with instabilities that arise during phase transition. Mathematical thermodynamics has notably been investigated by Guggenheim ^[1], Shapiro and Shapley ^[2], Aris ^[3], Truesdell ^[4], Krambeck ^[5], Pousin ^[6], Helluy and Mathis ^[7], and Giovangigli and Matuszewski ^[8,9]. Thermodynamics is also essential for analyzing the mathematical structure of hyperbolic-parabolic symmetrizable systems of partial differential equations modeling fluids. Hyperbolic systems have notably been studied by Friedrich and Lax ^[10], Ruggeri ^[11], Godlevski and Raviart ^[12] and Dafermos ^[13] and hyperbolic-parabolic systems by Volpert and Hudjaev ^[14], Kawashima and Shizuta ^[15], and Giovangigli et al. ^{[16,17,18,19]}. Classical multicomponent fluid thermodynamics are often built in practice from equations of states using ideal gas mixtures as the low-density limit ^[9]. Extended multicomponent thermodynamics required for diffuse interface fluid models are then obtained by adding gradient squared terms arising from molecular cohesive forces. For single species fluids, the thermodynamics of diffuse interface models has been built by Van der Walls ^[20,21], the capillary tensor derived by Korteweg ^[22] and the heat flux by Dunn and Serrin ^[23]. Diffuse interface models have been studied by de Gennes ^[24], Rowlinson ^[25], and Anderson et al. ^[26], numerical simulations with artificial interface thickening by Jamet et al. ^[27,28], and applications are presented by Gaillard et al. ^[29,30], Nayigizente et al. ^[31], and Le Calvez ^[32]. These equations have alternatively been obtained from Hamiltonian considerations by Gavrilyuk and Shugrin ^[33] and the links with the kinetic theory of dense gases addressed by Rocard ^[34], Barbante and Frezzotti ^[35], and Giovangigli ^[36]. For fluid mixtures, Cahn and Hilliard ^[37] first used a mole fraction gradient squared term in the free energy and later a density gradient in order to develop a thermodynamic formalism ^[38]. Rational thermodynamic aspects of Cahn-Hilliard fluids have been investigated by Falk ^[39], Kim and Lowengrub ^[40], Alt ^[41], Abels et al. ^[42], and Guo and Lin ^[43]. Cahn-Hilliard fluid models have further been derived from the kinetic theory of dense gas mixtures by Giovangigli ^[44]. Application of rational thermodynamics to Cahn-Hilliard type models have also presented by Lee et al. ^[45], Wang and Wise ^[46], Wang et al. ^[47], and Zhang et al. ^[48]. Mathematical aspects are discussed by Miranville ^[49], application to surface diffusion by Bretin et al. ^[50] and the physics of Cahn-Hilliard fluids with applications to water/air interfaces by Benilov ^[51]. These are strong motivations for investigating mathematically the structure of multicomponent fluid thermodynamics allowing instabilities, their construction from pressure laws, and their extended version including cohesive or capillary effects and the corresponding diffuse interface fluid models.

We initially discuss the mathematical structure and properties of classical thermodynamic functions in terms of intensive variables indispensable for fluid models. The complex situations of polymers ^[52], multitemperature flows ^[53,54], and links with statistical mechanics investigated by Fowler ^[55], Ferziger and Kaper ^[56], Keizer ^[57], Laasonen et al. ^[58], and Chen and Doolen ^[59], lay beyond the scope of this work. General physical formalisms including conservation laws and spatial transport like the thermodynamics of irreversible processes ^[60] or the GENERIC formalism of H. C. Öttinger ^[61] also lay beyond the scope of this work. A mathematical definition of thermodynamics is introduced and involves smoothness and homogeneity properties of thermodynamic functions, Gibbs' relation, and the compatibility with ideal gases at low density. We characterize thermal, mechanical and chemical thermodynamic stability. Thermal stability is generally guaranteed in practical applications and instabilities may then only be of mechanical or chemical type and determined by the spectrum of a single symmetric matrix $\Lambda$ related to the species chemical potentials or Gibbs' functions ^[9,29,30]. The natural variable is denoted by $\zeta = (T, v, \mathrm{y}_1, \ldots, \mathrm{y}_ n)^t$ where $T$ is the absolute temperature, $v$ the volume per unit mass, $\mathrm{y}_1, \ldots, \mathrm{y}_ n$ the species mass fractions supposed to be independent, and $n$ the number of species. The Gibbsian thermodynamic model is specified as $e(\zeta)$, $p(\zeta)$, and $s(\zeta)$, where $e$ denotes the energy per unit mass, $p$ the pressure, and $s$ the entropy per unit mass, over an open set $\mathcal{O}_ \zeta\subset {\mathbb R}^{2+ n}$. We also address a different structure associated with the thermodynamic variable $\mathcal{Z} = (T, \rho_1, \ldots, \rho_ n)^t$, where $\rho_1, \ldots, \rho_ n$ denotes the species masses per unit volume. Using the variable $\zeta = (T, v, \mathrm{y}_1, \ldots, \mathrm{y}_ n)^t$ leads to homogeneous thermodynamic functions—with singular entropy Hessians—and using the variable $\mathcal{Z} = (T, \rho_1, \ldots, \rho_ n)^t$ leads to nonhomogeneous thermodynamic functions and nonsingular entropy Hessian matrices when stability holds. These two frameworks are useful and have been investigated in the literature generally for ideal gas mixtures ^{[1,2,3,4,5,6,7,8,9]}. We establish the mathematical equivalence between these thermodynamic formalisms under proper transformation rules. The thermodynamic functions are not defined for all states because of nonidealities and thermodynamic instabilities—associated with phase changes—are allowed. We also discuss the thermodynamic—part of the—entropic $\mathcal{V}$ and normal $\mathcal{W}$ variables. The entropic variable $\mathcal{V}$ is closely associated with symmetrization properties for systems of conservation laws. However, the entropic variable $\mathcal{V}$ is not practical in the presence of phase changes where it is more interesting to use a thermodynamic—part of a—normal variables $\mathcal{W}$. We investigate conditions allowing the use of such a normal variable $\mathcal{W}$ in multicomponent flows.

We then study the mathematical construction of such fluid thermodynamics allowing instabilities from a pressure law, extending previous work solely devoted to the situation of supercritical fluids ^[9]. Such a procedure is often used to model multicomponent fluids using typically a cubic equation of states. Nonideal pressure laws have notably been discussed by Gillespie ^[62], Benedict et al. ^[63], Beattie ^[64], Redlich and Kwong ^[65], Soave ^[66,67], Peng and Robinson ^[68], Graboski and Daubert ^[69], Ozowkelu and Erbar ^[70], Harstadt et al. ^[71], Congiunti et al. ^[72], Colonna and Silva ^[73], and Cañas-Martin et al. ^[74,75]. We investigate mathematically natural structural assumptions on a pressure law that may then define a Gibbsian thermodynamics using ideal gas mixtures as the low-density limit. Thermodynamic stability is not required a priori and unstable states associated with phase changes are taken into account. The resulting construction appears to be more general and more natural than that obtained in previous work ^[9] only focused on supercritical states. The structural assumptions are the 0-homogeneity of the pressure $p(\zeta)$ with respect to $(v, \mathrm{y}_1, \ldots, \mathrm{y}_ n)$ and a quadratic expression of $p - p^ {\rm{id}}$ as $\rho$ goes to zero where $p^ {\rm{id}}$ denotes the ideal gas pressure. We determine the energy and entropy per unit mass and then establish that they constitute a thermodynamics according to the previous definition. We also also introduce rescaled versions of the stability matrix $\Lambda$ that behave smoothly for vanishing mass fractions.

Cubic equations of state, like the Soave-Redlich-Kwong ^[65,66] or the Peng-Robinson ^[68] pressure laws, have been found to be convenient since they may be inverted using Cardan's formula and give accurate results over the range of pressure, temperature, and mixture states of various applications. Such cubic equations of state have notably been used by Saur et al. ^[76], Meng and Yang ^[77], Oefelein ^[78], Ribert ^[79], Giovangigli et al. ^[80], and Gaillard et al. ^[29,30,81]. The Benedict-Webb-Rubin equation of state ^[63,67] or its generalizations are more accurate but uneasy to handle and contain many approximating parameters. We investigate in particular the Soave-Redlicht-Kwong equation of state and the corresponding Gibbsian thermodynamics compatible with that of ideal gas mixtures. The necessary conditions to define a thermodynamics are established and the entropy and energy per unit mass are obtained explicitly in terms of the natural variable. The matrix $\Lambda$ associated with mechanical and chemical stability is also obtained explicitly. The parameters of the Soave-Redlich-Kwong (SKR) equations of state are related to critical state for pure stable substances or to Lennard-Jones interaction potentials for unstable species.

We next investigate the mathematical structure of an extended thermodynamics associated with a simplified diffuse interface fluid model ^{[20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51]}. Diffuse interface models, like those of van der Waals or Cahn and Hilliard, are obtained in the presence of internal cohesive forces and involve density gradient squared terms in the free energies. Such diffuse interface models describe interphases—interfaces between phases—as regions with smooth variations of physical properties and have been used successfully in order to describe spinodal decomposition, droplets dynamics, three phase contact lines, and surface diffusion, as well as transcritical flames ^{[20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51]}. The simplified multicomponent flow model investigated in this paper is obtained by assuming that the capillarity coefficients are identical and derived in appendix C. In this situation, only the density gradient ${\boldsymbol\nabla}\rho$ is involved in thermodynamic functions, and this model has been used successfully in the numerical simulation of high pressure flames ^[29,30]. We further investigate the corresponding augmented formulation obtained by introducing a local representation ${\boldsymbol{w}}$ of the gradient of density following Gavrilyuk and Gouin ^[82], Benzoni et al. ^[83], Bresch et al. ^[84,85], Kotschote ^[86], and Giovangigli et al. ^[87]. The mathematical structure of the extended thermodynamics is obtained in a similar way as for classical thermodynamics. Such augmented formulations may be used in order to investigate the existence of solutions ^{[83,85,86,87]} as well as for numerical simulations of diffuse interface fluids ^[84]. We also address the mathematical structure of chemical production terms for completeness even though chemical reactions will not be considered in the numerical simulations. The chemical production rates involve, in particular, Marcelin's expression for the chemical reaction rates of progress discussed by Marcelin ^[88,89], Keizer ^[57], and Giovangigli et al. ^[8,80,90,91].

We then specialize the multicomponent diffuse interface model to the situation of self-similar strained flows discussed notably be Ribert et al. ^[79], Pons et al. ^[92], and Giovangigli and Matuszewski ^[93]. The low Mach number regime is analyzed and the simplified strained flow equations are obtained by using a similarity assumption. The evaluation of high pressure transport coefficients is briefly addressed as well as how to avoid exploding coefficients at mechanical instabilities. The numerical tools used in the simulations include application independent libraries for the evaluation of thermodynamic properties ^[94,95]. Nonlinear solvers based on Newton's method ^[96] and continuation techniques ^[97] are also used for the evaluation of equilibrium curves or critical states. Two points boundary value problems are solved using finite difference methods and adaptive grids as well as optimized chemistry and transport libraries. Newton's method are notably discussed by Deuflhard ^[96], self-adaptive grids by Smooke ^[98,99] and Oran and Boris ^[100], continuation techniques by Keller ^[97] and Giovangigli and Smooke ^[101], optimized chemistry and transport libraries by Giovangigli et al. ^[102,103] and Matuszewski ^[104], and fast and accurate evaluation of transport properties by Ern and Giovangigli ^{[105,106,107,108,109,110]} with applications to thermal diffusion ^[111].

The thermodynamic formalism is validated by comparison with experimental measurements for mixtures of ethane C$_2$H$_6$ and nitrogen N$_2$. High pressure thermodynamic data and transport coefficients as well as experimental results are obtained from the literature. We mention notably the classification of Van Konynenburg and Scott ^[112,113], the high pressure correlations of Ely and Hanley ^[114], Chung et al. ^[115], the kinetic results of Kurochkin et al. ^[116], the thermodynamic data of Kee et al. ^[117], Chase ^[118], and McBride et al. ^[119], the gradient modeling of Lin et al. ^[120], the experimental results of Youglove ^[121], Youglove and Ely ^[122], Gupta et al. ^[123], Eakin et al. ^[124], Stryjek et al. ^[125], Wisotski et al. ^[126], Japas and Frank ^[127], and the analysis of critical point of Peng and Robinson ^[128]. We first investigate specific heats with a very good agreement with the experimental measurements of Youglove and Ely ^[121,122]. High pressure and low temperature mixtures of ethane C$_2$H$_6$ and nitrogen N$_2$ may be chemically unstable. These mixtures may split between an ethane-rich gaseous-like phase and an ethane-poor liquid-like phase. The two phase chemical equilibrium curves are found to be in very good agreement with experimental results of Gupta et al. ^[123]. We also investigate the critical points of ethane-nitrogen mixtures. Following the classification of Van Konynenburg and Scott ^[112,113], the curve of critical points is known to be of Type Ⅲ ^{[124,125,126]}. The numerical simulations confirm such a type Ⅲ behavior with an overall very satisfactory agreement in comparison with experimental results of Eakin at al. ^[124], Stryjek et al. ^[125], and Wisotzki and Schneider ^[126]. We finally investigate the structure of ethane liquid-vapor interfaces and study the influence of the capillarity coefficients on the flow structure.

The thermodynamic formalism presented in this paper could naturally be used in association with other fluid models in order to investigate mathematical as well as numerical issues. We may mention, for instance, gradient flow structures that have been obtained in the absence of convection phenomena ^[129,130], numerical stability issues ^[131], or problems with boundary conditions.

Thermodynamics is investigated in Section 2 for the homogeneous case and in Section 3 for the nonhomogeneous situation. The construction of the thermodynamics is discussed in Section 4. The Soave-Redlich-Kwong equation of state is investigated in Section 5. Extended thermodynamics and diffuse interface models are presented in Section 6 and strained flow in Section 7. Numerical simulations of C$_2$H$_6$-N$_2$ mixtures are finally considered in Section 8.

2.   Homogeneous Gibbsian thermodynamics

We investigate in this section the mathematical structure of homogeneous multicomponent fluid thermodynamics compatible at low density with that of ideal gases.

2.1.   The natural variable

We investigate thermodynamics in terms of the natural variable $\zeta = (T, v, \mathrm{y}_1, \ldots, \mathrm{y}_ n)^t$ where $T$ denotes the absolute temperature, $v$ the volume per unit mass, i.e., $v = 1/\rho$ where $\rho$ is the mass density, $\mathrm{y}_1, \ldots, \mathrm{y}_ n$ the species mass fractions, $n\geqslant1$ the number of species, and ${\mathfrak S}$ the species indexing set ${\mathfrak S} = \{1, \ldots, n\}$. The thermodynamic model is specified as $e(\zeta)$, $p(\zeta)$, and $s(\zeta)$, defined on an open set $\mathcal{O}_ \zeta\subset(0, \infty)^{2+ n}$, where $e$ denotes the energy per unit mass, $p$ the pressure, and $s$ the entropy per unit mass. The open set $\mathcal{O}_ \zeta$ frequently differs from $(0, \infty)^{2+ n}$ because of real gas effects. We often commit the traditional abuse of notation of denoting by the same symbol a given quantity as function of different state variables unless confusion may arise. The superscript ${}^ {\rm{id}}$ is associated with ideal gas mixtures thermodynamics.

We denote by $\partial^{}$ the derivation operator with respect to the variable $\zeta = (T, v, \mathrm{y}_1, \ldots, \mathrm{y}_ n)^t$, by $\mathrm{d}$ the total differential operator, and for any $\lambda > 0$, we define $\zeta_\lambda^{} = (T, \lambda v, \lambda \mathrm{y}_1, \ldots, \lambda \mathrm{y}_ n)^t$. We denote by $\mathrm{y}$ the mass fraction vector $\mathrm{y} = (\mathrm{y}_1, \ldots, \mathrm{y}_ n)^t$, ${1\kern -2.7pt \mathrm{I}}$ the vector with unit components ${1\kern -2.7pt \mathrm{I}} = (1, \ldots, 1)^t$, $\langle, \rangle$ the Euclidean scalar product, and define $\textsf{ $\Sigma$ } = \{ \mathrm{y}\in(0, \infty)^ n; \ \langle \mathrm{y}, {1\kern -2.7pt \mathrm{I}}\rangle = 1\}$. Physically, the mass fraction $\mathrm{y}$ must sum up to unity $\mathrm{y} \in \textsf{ $\Sigma$ }$. However, as usual in the modeling of multicomponent flows, we assume that the mass fractions are independent and the mass constraint $\langle \mathrm{y}, {1\kern -2.7pt \mathrm{I}}\rangle = 1$ then results from the governing equations and boundary conditions ^[8,60]. The formalism developed in this paper also applies when $n = 1$, but in this special situation, the variable $\mathrm{y}_1$ is usually formally eliminated as discussed in Appendix A. We denote by $\gamma\in {\mathbb N}$, $\gamma\geqslant2$, the regularity class of thermodynamics functions.

Definition 2.1. Let $e$, $p$, and $s$, be functions of the variable $\zeta$ defined over an open set $\mathcal{O}_ \zeta \subset (0, \infty)^{2+ n}$. These functions are said to define a thermodynamics if properties $(\textsf{T}_{\textsf{0}}-\textsf{T}_\textsf{2})$ hold.

$(\textsf{T}_{\textsf{0}})$ The functions $e$, $p$, and $s$ are $C^ \gamma$ over a nonempty, connected, open set $\mathcal{O}_ \zeta \subset (0,\infty)^{2+ n}$. For any $(\lambda, \zeta) \in (0,\infty){\times} \mathcal{O}_ \zeta$, $\zeta_\lambda^{}\in \mathcal{O}_ \zeta$, $e( \zeta_\lambda^{}) = \lambda e( \zeta)$, $p( \zeta_\lambda^{}) = p( \zeta)$, and $s( \zeta_\lambda^{}) = \lambda s( \zeta)$.

$(\textsf{T}_{\textsf{1}})$ For any $\zeta\in \mathcal{O}_ \zeta$, defining $g_k = \partial^{}_{ \mathrm{y}_k}\! e - T \partial^{}_{ \mathrm{y}_k}\! s$, $k\in {\mathfrak S}$, we have Gibbs' relation:

$(\textsf{T}_{\textsf{2}})$ For any $(T, \mathrm{y}_1\ldots, \mathrm{y}_ n)^t \in(0, \infty)^{1+ n}$, there exists $v_{{{\rm{m}}}}$ so that $v > v_{{{\rm{m}}}}$ implies $(T, v, \mathrm{y}_1\ldots, \mathrm{y}_ n)^t \in \mathcal{O}_ \zeta$ and

Property $(\textsf{T}_\textsf{0})$ expresses the smoothness and the homogeneity properties of Gibbsian thermodynamics written in terms of the natural variable $\zeta$. Temperature, volume and mass fractions are assumed to be positive with $\mathcal{O}_ \zeta \subset (0,+\infty)^{2+ n}$. Property $(\textsf{T}_{\textsf{1}})$ is Gibbs' relation for fluid mixtures written with the $\zeta$ variable. With the simplified definition of the species, Gibbs functions $g_k$, $k\in {\mathfrak S}$, (2.1) are equivalent to $T \partial^{}_T s = \partial^{}_T e$ and $T \partial^{}_ v s = \partial^{}_ v e + p$. Property $\textsf{T}_{\textsf{2}})$ is the compatibility condition with ideal gases since for large $v$, we must recover the ideal gas thermodynamics. The multiplication by $v$ is required for the pressure since as $v\to\infty$, both $p$ and $p^ {\rm{id}}$ goes to zero. In comparison with previous work devoted to supercritical fluids, we do not assume that thermodynamic stability holds in the definition of the thermodynamics ^[9].

Lemma 2.2. Assuming $(\textsf{T}_{\textsf{0}}-\textsf{T}_{\textsf{1}})$ and defining the mixture Gibbs function by $g = e + p v - T s$, we have the relation $g = \sum_{k\in {\mathfrak S}} \mathrm{y}_k g_k$.

Proof. Since $e$ and $s$ are 1-homogeneous with respect to $( v, \mathrm{y}_1,\ldots, \mathrm{y}_ n)$, we obtain from Euler's relation that $e = v \partial^{}_ v e + \sum_{k\in {\mathfrak S}} \mathrm{y}_k \partial^{}_{ \mathrm{y}_k} e$ and $s = v \partial^{}_ v s + \sum_{k\in {\mathfrak S}} \mathrm{y}_k \partial^{}_{ \mathrm{y}_k} s$. On the other hand, we deduce from Gibbs' relation (2.1) that $T \partial^{}_ v s = p + \partial^{}_ v e$ and $T \partial^{}_{ \mathrm{y}_k} s = \partial^{}_{ \mathrm{y}_k} e - g_k$. Forming finally $g = e + p v - T s$, we then obtain that $g = \sum_{k\in {\mathfrak S}} \mathrm{y}_k g_k$. □

2.2.   The energy-based variable

We introduce the energy-based variable $\xi = (e, v, \mathrm{y}_1,\ldots, \mathrm{y}_ n)^t$ that may be used whenever the map $\zeta\mapsto \xi$ is at least a local diffeomorphism from a subset $\mathcal{O}_ \zeta'$ of $\mathcal{O}_ \zeta$ onto an open set $\mathcal{O}_ \xi'$, and denote by $\widetilde\partial^{}$ the derivation operator with respect to the variable $\xi = ( e, v, \mathrm{y}_1,\ldots, \mathrm{y}_ n)^t$. Both variables $\zeta$ and $\xi$ are practical and useful since $\zeta$ is naturally used in equations of state—or for compatibility conditions with ideal gases—and $\xi$ is more naturally associated with conservation laws and Gibbs' relation.

Lemma 2.3. Assume that the map $\zeta \to \xi$ is a $C^ \gamma$ local diffeomorphism from an open set $\mathcal{O}_ \zeta'$ onto an open set $\mathcal{O}'_ \xi$. Then, Gibbs' relation is locally in the form

and we have, in particular, $T \widetilde\partial^{}_ e s = 1$, $T \widetilde\partial^{}_\nu s = p$, and $T \widetilde\partial^{}_{ \mathrm{y}_y} s = - g_k$.

Proof. The map $\zeta\mapsto \xi$ is a $C^ \gamma$ local diffeomorphism, and we may use the variable $\xi$. Gibbs' relation (2.1) is then easily rewritten in the form (2.3) since $\mathrm{d} e = \partial^{}_T e \; \mathrm{d} T + \partial^{}_ v e \, \mathrm{d} v + \sum_{k\in {\mathfrak S}} \partial^{}_{ \mathrm{y}_k}\! e \, \mathrm{d} \mathrm{y}_k$. □

It is easily deduced from Gibbs relation that the equilibrium conditions between two mixture states denoted with the sharp and flat symbols are $T^\sharp = T^\flat$, $p^\sharp = p^\flat$, and $g_i^\sharp = g_i^\flat$ for all $i\in \mathfrak{S}$ ^[1,2,3,5]. We now give a sufficient conditions $(\textsf{T}_{\textsf{3}} )$, often met in practice, in such a way that the map $\zeta \mapsto \xi$ are a global diffeomorphism. The nonempty, open set $\mathcal{O}_ \zeta$ is said to be increasing with temperature if for any $(T, v, \mathrm{y}_1,\ldots, \mathrm{y}_ n)^t\in \mathcal{O}_ \zeta$ and any $T<T'$, we have $(T', v, \mathrm{y}_1,\ldots, \mathrm{y}_ n)^t\in \mathcal{O}_ \zeta$.

$(\textsf{T}_{\textsf{3}} )$ For any $\zeta\in \mathcal{O}_ \zeta$ $\partial^{}_T e( \zeta)>0$, and $\mathcal{O}_ \zeta$ is increasing with temperature.

Lemma 2.4. Assume that $(\textsf{T}_{\textsf{0}} )$ and $(\textsf{T}_{\textsf{3}} )$ hold, then the map $\zeta \mapsto \xi$ is a $C^ \gamma$ global diffeomorphism from $\mathcal{O}_ \zeta$ onto a nonempty, connected, open set $\mathcal{O}_ \xi$.

Proof. The map $\zeta\mapsto \xi$ is $C^ \gamma$ and its Jacobian matrix is invertible over $\mathcal{O}_ \zeta$ since $\partial^{}_T e>0$. From the local inversion theorem, we know that the image $\mathcal{O}_ \xi$ of $\mathcal{O}_ \zeta$ is an open set. Moreover, the map $\zeta\mapsto \xi$ is one-to-one since if $( e^\flat, v^\flat, \mathrm{y}_1^\flat,\ldots, \mathrm{y}_ n^\flat) = ( e^\sharp, v^\sharp, \mathrm{y}_1^\sharp,\ldots, \mathrm{y}_ n^\sharp)$ denoting $v = v^\flat = v^\sharp$ and $\mathrm{y}_k = \mathrm{y}_k^\flat = \mathrm{y}_k^\sharp$ for $k\in {\mathfrak S}$, then $e(T^\flat, v, \mathrm{y}_1,\ldots, \mathrm{y}_ n) = e(T^\sharp, v, \mathrm{y}_1,\ldots, \mathrm{y}_ n)$. The structure of $\mathcal{O}_ \zeta$ then guarantees that the segment $[ \zeta^\flat, \zeta^\sharp]$ given by $\bigl( \alpha T^\flat +(1-\alpha)T^\sharp, v, \mathrm{y}_1,\ldots, \mathrm{y}_ n \bigr)^t$ for $0\le \alpha\le 1$ lies in $\mathcal{O}_ \zeta$ so that

Since $\partial^{}_T e$ is positive over $[ \zeta^\flat, \zeta^\sharp]$, this relation implies that $T^\flat = T^\sharp$. The map $\zeta\mapsto \xi$ is thus a $C^ \gamma$ global diffeomorphism from $\mathcal{O}_ \zeta$ onto an open set $\mathcal{O}_ \xi$ that is nonempty and connected since $\mathcal{O}_ \zeta$ is nonempty and connected. □

Under assumptions $(\textsf{T}_{\textsf{0}})$ and $(\textsf{T}_{\textsf{3}} )$, the map $\zeta \mapsto \xi$ is a $C^ \gamma$ global diffeomorphism and it is then possible to characterize the thermodynamic properties in terms of $\xi$. We now investigate thermodynamic stability properties of fluid mixtures. From the second principle of thermodynamics, an isolated system tends to maximize its entropy. The entropy of a stable isolated system should thus be a concave function of its volume, composition variables, and internal energy. If it is not the case, the system may split between two or more phases in order to reach equilibrium. The Hessian matrix $\widetilde\partial^2_{ \xi \xi} s$ should thus be negative semi-definite with a null space reduced to ${\mathbb R} \xi$ solely arising from homogeneity. Thermodynamic stability properties naturally involve the vector of species mass fraction $\mathrm{y} = ( \mathrm{y}_1,\ldots, \mathrm{y}_ n)^t$ and the matrix $\Lambda = (\Lambda_{kl})_{k,l\in {\mathfrak S}}$ of size $n$ is defined by

as well as the matrix $\widehat\Lambda$ defined when $\langle \Lambda \mathrm{y}, \mathrm{y}\rangle\neq0$ by

We first establish a few useful properties of the stability matrix $\Lambda$.

Proposition 2.5. Assuming that $(\textsf{T}_{\textsf{0}} -\textsf{T}_{\textsf1} )$ are satisfied, for any $\zeta\in \mathcal{O}_ \zeta$, we have the relations

Proof. From Gibbs' relation, we have $T \partial^{}_T s = \partial^{}_T e$, $T \partial^{}_ v s = \partial^{}_ v e + p$, $T \partial^{}_{ \mathrm{y}_k} s = \partial^{}_{ \mathrm{y}_k} e - g_k$, and this implies the compatibility relations

The matrix $\Lambda$ is first symmetric since $\partial^{}_{ \mathrm{y}_l} g_k = \partial^2_{ \mathrm{y}_k \mathrm{y}_l}\! e - T \partial^2_{ \mathrm{y}_k \mathrm{y}_l}\! s$, and in order to establish (2.6) and (2.7), we note that $g_k(\zeta)$ is 0-homogeneous in $(v, \mathrm{y}_1, \ldots, \mathrm{y}_ n)$ since $g_k = \partial^{}_{ \mathrm{y}_k}\! e - T \partial^{}_{ \mathrm{y}_k}\! s$, $T$ is 0-homogeneous, and $e$ and $s$ are 1-homogeneous. This implies the relation $\sum_{l\in {\mathfrak S}} \mathrm{y}_l \partial^{}_{ \mathrm{y}_l} g_k = - v \partial^{}_ v g_k$ and, thus, $\sum_{k\in {\mathfrak S}} \Lambda_{kl} \mathrm{y}_l = - v \partial^{}_ v g_k/T$, and (2.6) is established by using (2.8). Upon multiplying (2.6) by $\mathrm{y}_k$ and summing over $k$, we next obtain that $T \langle\Lambda \mathrm{y}, \mathrm{y}\rangle = - v \sum_{k\in {\mathfrak S}} \mathrm{y}_k \partial^{}_{ v} g_k$. However, from the compatibility relation (2.8), we get $\partial^{}_{ v} g_k = - \partial^{}_{ \mathrm{y}_k} p$, and thanks to the 0-homogeneity of $p(\zeta)$ with respect to $(v, \mathrm{y}_1, \ldots, \mathrm{y}_ n)$, we conclude that $\sum_{k\in {\mathfrak S}} \mathrm{y}_k \partial^{}_{ v} g_k = - \sum_{k\in {\mathfrak S}} \mathrm{y}_k \partial^{}_{ \mathrm{y}_k} p = v \partial^{}_ v p$ and (2.7) is established. □

We may now characterize thermodynamic stability with the help of the matrix $\Lambda$ and the positivity of the specific heat at constant volume. Defining the free energy as $f = e - T s$, we note that $g_k = \partial^{}_{ \mathrm{y}_k} f$ as well as $\Lambda_{kl} = \partial^2_{ \mathrm{y}_k \mathrm{y}_l} f/T$ and $\Lambda = \partial^2_{ \mathrm{y} \mathrm{y}} f/T$. The energy, entropy, and free energy per unit volume are also defined by $\mathcal{E} = \rho e$, $\mathcal{S} = \rho s$, and $\mathcal{F} = \rho f$. The stability conditions obtained in the following proposition may be split between thermal stability, mechanical stability, and chemical stability conditions.

Proposition 2.6. Assuming that $(\textsf{T}_{\textsf{0}} - \textsf{T}_{\textsf{1}} )$ are satisfied and that $\zeta\mapsto \xi$ is a $C^ \gamma$ local diffeomorphism, the following statements are equivalent :

(i) $\widetilde\partial^2_{ \xi \xi} s$ is negative semi-definite with null space $N( \widetilde\partial^2_{ \xi \xi} s) = {\mathbb R} \xi$.

(ii) ${\partial_T} e >0$ and $\Lambda$ is positive definite.

(iii) ${\partial_T} e >0$, $\partial_ v p<0$, and $\widehat\Lambda$ is positive semi-definite with null space ${\mathbb R} \mathrm{y}$.

The proof, lengthy but technically not difficult, is presented in Appendix B. The property $\partial^{}_T e>0$ is usually termed the thermal stability condition, the property $\partial^{}_ v p<0$ the mechanical stability condition, and the condition that $\widehat\Lambda$ is positive semi-definite with null space $N(\widehat\Lambda) = {\mathbb R} \mathrm{y}$ is usually termed the chemical stability condition ^[1,9]. In particular, the property that $\Lambda$ is positive definite encompasses both the mechanical and chemical stability conditions.

2.3.   Other variables

We discuss in this section the pressure variable $\pi = (T, p, \mathrm{y}_1,\ldots, \mathrm{y}_ n)^t$ as well as molar properties. We denote by $\widehat\partial^{}$ the derivation operator with respect to the variable $\pi = (T, p, \mathrm{y}_1,\ldots, \mathrm{y}_ n)^t$.

Proposition 2.7. Assume that $(\textsf{T}_{\textsf{0}} - \textsf{T}_{\textsf1} )$ hold and that $\zeta\mapsto \pi$ is a local $C^ \gamma$ diffeomorphism from $\mathcal{O}_ \zeta'\subset \mathcal{O}_ \zeta$ onto an open set $\mathcal{O}_ \pi'$. Then, we have the differential relations $g_k = \widehat\partial^{}_{ \mathrm{y}_k} g$, $k\in {\mathfrak S}$, and Gibbs' relation may be written:

Proof. For any smooth function $\chi$ of $\zeta$ or $\pi$, we have the following differential relations associated with the variables $\pi$ and $\zeta$:

Letting $g = e + p v - T s$, we first obtain that $\widehat\partial^{}_T g = \widehat\partial^{}_T e + p \widehat\partial^{}_T v - s - T \widehat\partial^{}_T s$, and from (2.10), we obtain the relations $\widehat\partial^{}_T s = \partial^{}_T s + \partial^{}_ v s\; \widehat\partial^{}_T v$ and $\widehat\partial^{}_T e = \partial^{}_T e + \partial^{}_ v e\; \widehat\partial^{}_T v$. The relations (2.10) then imply that $\widehat\partial^{}_T e + p \widehat\partial^{}_T v -T \widehat\partial^{}_T s = 0$ since $\partial^{}_T e -T \partial^{}_T s = 0$ and $(\partial^{}_ v e + p -T \partial^{}_ v s) \widehat\partial^{}_T v = 0$ from $(\textsf{T}_{\textsf{1}} )$, and we have established that $\widehat\partial^{}_T g = - s$. Similarly, we have $\widehat\partial^{}_ p g = \widehat\partial^{}_ p e + v + p \widehat\partial^{}_ p v - T \widehat\partial^{}_ p s$, and from (2.11), $\widehat\partial^{}_ p e = \partial^{}_ v e\; \widehat\partial^{}_ p v$ and $\widehat\partial^{}_ p s = \partial^{}_ v s\; \widehat\partial^{}_ p v$ so that using again $(\textsf{T}_{\textsf{1}} )$, we have established that $\widehat\partial^{}_ p g = v$. Finally, we have $\widehat\partial^{}_{ \mathrm{y}_k} g = \widehat\partial^{}_{ \mathrm{y}_k} e + p \widehat\partial^{}_{ \mathrm{y}_k} v - T \widehat\partial^{}_{ \mathrm{y}_k} s$, and from (2.12), the relations $\widehat\partial^{}_{ \mathrm{y}_k} e = \partial^{}_{ \mathrm{y}_k} e + \partial^{}_ v e\; \widehat\partial^{}_{ \mathrm{y}_k} v$ and $\widehat\partial^{}_{ \mathrm{y}_k} s = \partial^{}_{ \mathrm{y}_k} s + \partial^{}_ v s\; \widehat\partial^{}_{ \mathrm{y}_k} v$, and after some algebra, we obtain so that $\widehat\partial^{}_{ \mathrm{y}_k} g = g_k$ using also $(\textsf{T}_{\textsf{1}} )$. In conclusion, we have established that $\widehat\partial^{}_T g = - s$, $\widehat\partial^{}_ p g = v$, and $\widehat\partial^{}_{ \mathrm{y}_k} g = g_k$, and (2.9) then follows. □

We may now give a new characterization of the matrix $\widehat\Lambda$ associated with chemical stability.

Proposition 2.8. Assume that $(\textsf{T}_{\textsf{0}} - \textsf{T}_{\textsf{1}} )$ hold and that $\zeta\mapsto \pi$ is a local $C^ \gamma$ diffeomorphism from $\mathcal{O}_ \zeta' \subset \mathcal{O}_ \zeta$ onto an open set $\mathcal{O}_ \pi'$. Then, for any $\zeta\in \mathcal{O}_ \zeta'$, we have the matrix identity

Proof. From Proposition 2.7, Gibbs relation may be written in the form (2.9). We thus deduce that $\widehat\partial^{}_T g = - s$, $\widehat\partial^{}_ p g = v$, $\widehat\partial^{}_{ \mathrm{y}_k} g = g_k$, and this implies the compatibility relations

Upon letting $\chi = g_l$ in (2.12), we obtain that $\widehat\partial^2_{ \mathrm{y}_k \mathrm{y}_l} g = \widehat\partial^{}_{ \mathrm{y}_k} g_l = \widehat\partial^{}_{ \mathrm{y}_l} g_k$ can be written $\widehat\partial^{}_{ \mathrm{y}_k} g_l = \partial^{}_{ \mathrm{y}_k} g_l - \widehat\partial^{}_ p g_l \partial^{}_{ \mathrm{y}_k} p$ so that from (2.8) and (2.11), $\widehat\partial^{}_{ \mathrm{y}_k} g_l = \partial^{}_{ \mathrm{y}_k} g_l + \partial^{}_ v g_k \partial^{}_ v g_l \widehat\partial^{}_ p v$. This implies the identity (2.13) thanks to $\partial^{}_ v p\; \widehat\partial^{}_ p v = 1$, (2.6), and (2.7). □

Thermodynamic stability thus involves the matrix $\widetilde\partial^2_{ \zeta \zeta} s$, mechanical and chemical stability of the matrix $\Lambda = \partial^2_{ \mathrm{y} \mathrm{y}} f$, and chemical stability of the matrix $\widehat\partial^2_{ \mathrm{y} \mathrm{y}} g$. We now investigate when the change of variable $\zeta\mapsto \pi$ is a local diffeomorphism in such a way that the pressure variable $\pi$ may be used. The nonempty, open set $\mathcal{O}_ \zeta$ is said to be increasing with volume if for any $(T, v, \mathrm{y}_1, \ldots, \mathrm{y}_ n)^t\in \mathcal{O}_ \zeta$ and any $v < v'$, we have $(T, v', \mathrm{y}_1, \ldots, \mathrm{y}_ n)^t\in \mathcal{O}_ \zeta$.

Proposition 2.9. Let $\mathcal{O}_ \zeta'$ be an open subset of $\mathcal{O}_ \zeta$ that is increasing with volume. Assume also that the mechanical stability condition $\partial_ v p < 0$ holds on $\mathcal{O}_ \zeta'$. Then, $\zeta \mapsto \pi$ is a $C^ \gamma$ diffeomorphism from $\mathcal{O}_ \zeta'$ onto an open set $\mathcal{O}_ \pi'$.

Proof. The proof is similar to that of Lemma 2.4. □

The mechanical stability condition usually does not hold when evaporation or condensation is taking place so that $\pi$ is not a good variable when there are phase changes. We essentially investigate thermodynamic properties per unit mass or per unit volume in this paper. The governing equations of fluid mixtures are actually associated with mass, momentum, and energy conservation so that mass based variables are natural. Thermodynamic properties may still be expressed equivalently per unit mole using straightforward transformation rules. The species mole fractions $\mathsf{X}_1, \ldots, \mathsf{X}_ n$ are defined by

where $m$ is the molar mass of the mixture. The factor $\sum_{i\in {\mathfrak S}} \mathrm{y}_i$ in the definition of the $m$ guarantees that $m$ is 0-homogeneous ^[8], that the mass/mole relations are invertible, and that $\sum_{i\in {\mathfrak S}} \mathrm{y}_i = \sum_{i\in {\mathfrak S}} \mathsf{X}_i$ and $\bigl(\sum_{i\in {\mathfrak S}} \mathsf{X}_i\bigr) m = \sum_{i\in {\mathfrak S}} \mathsf{X}_im_i$. As a typical example, the mixture Gibbs function per unit mole ${{\rm{G}}}$ and the species Gibbs function per unit mole ${{\rm{G}}}_i$, $i\in {\mathfrak S}$, may then be defined as

From the relation $g = \sum_{i\in {\mathfrak S}} \mathrm{y}_i g_i$, it is easily obtained that ${{\rm{G}}} = \sum_{i\in {\mathfrak S}} \mathsf{X}_i {{\rm{G}}}_i$. All thermodynamic properties per unit mass may similarly be rewritten per unit mole and the molar variables $({{\rm{E}}}, {{\rm{V}}}, \mathsf{X}_1, \ldots, \mathsf{X}_ n)^t$, $(T, {{\rm{V}}}, \mathsf{X}_1, \ldots, \mathsf{X}_ n)^t$, or $(T, p, \mathsf{X}_1, \ldots, \mathsf{X}_ n)^t$ may also be used, instead of $\xi$, $\zeta$, or $\pi$, where ${{\rm{E}}}$ denotes the energy per unit mole and ${{\rm{V}}}$ the volume per unit mole.

2.4.   Ideal gas mixtures

For completeness we describe in this section the case of ideal gas mixtures. The corresponding thermodynamic properties are denoted with the superscript ${}^ {\rm{id}}$ and involved in Property $(\textsf{T}_{\textsf{2}})$. In terms of the variable $\zeta = (T, v, \mathrm{y}_1,\ldots, \mathrm{y}_ n)^t$, the pressure of an ideal gas mixture is given by

and the internal energy by

where $e_k^ {\rm{id}}$ is the internal energy per unit mass of the $k$th species, $c_{ {{{\rm{v}}}} k}^ {\rm{id}}$ the specific heat per unit mass of the $k$ species, and $e_k^ \mathrm{st}$ the formation energy of the $k$th species per unit mass at the standard temperature $T^ \mathrm{st}$. The entropy per unit mass reads

where $s_k^ \mathrm{st}$ is the formation entropy of the $k$th species per unit mass at the standard temperature $T^ \mathrm{st}$ and standard pressure $p^ \mathrm{st}$, and $\gamma^ \mathrm{st} = p^ \mathrm{st}/ R T^ \mathrm{st}$ is the standard concentration. The enthalpy $h^ {\rm{id}} = e^ {\rm{id}} + p^ {\rm{id}}/ v$, Gibbs function $g^ {\rm{id}} = e^ {\rm{id}} + p^ {\rm{id}}/ v - T s^ {\rm{id}}$, and free energy $f^ {\rm{id}} = e^ {\rm{id}} - T s^ {\rm{id}}$ are then easily evaluated. The ideal gas specific heats, formation energies, and formation entropies satisfy the following assumptions.

$(\textsf{ID} )$ The formation energies $e^ \mathrm{st}_i$, $i\in {\mathfrak S}$, and entropies $s^ \mathrm{st}_i$, $i\in {\mathfrak S}$, are real constants. The species mass per unit mole $m_i$, $i\in {\mathfrak S}$, and the gas constant $R$ are positive constants. The species heat per unit mass $c_{ {{{\rm{v}}}} i}^ {\rm{id}}$, $i\in {\mathfrak S}$, are $C^ \gamma$ functions over $[0, \infty)$, and there exist constants $\underline{ c}_ {{{\rm{v}}}}$ and $\overline{ c}_ {{{\rm{v}}}}$ such that $0 < \underline{ c}_ {{{\rm{v}}}} \leqslant c_{ {{{\rm{v}}}} i}^ {\rm{id}} \leqslant \overline{ c}_ {{{\rm{v}}}}$ for all $T\geqslant0$ and $i\in {\mathfrak S}$.

We assume throughout the paper that Property $(\textsf{ID})$ holds. The extension up to zero temperature of specific heats, energies, and enthalpies is natural in thermodynamics. Note that the specific heats remain bounded away from zero since ideal gas mixtures are governed by Boltzmann statistics ^[8]. In the following proposition, we investigate the mathematical properties of ideal gas mixture thermodynamics and $( \textsf{T}_{\textsf{3}})$ is trivial in this situation.

Proposition 2.10. The energy per unit mass $e^ {\rm{id}}$, the ideal gas mixture pressure $p^ {\rm{id}}$, and the entropy per unit mass $s^ {\rm{id}}$ are $C^ \gamma$ functions defined on the open set $\mathcal{O}_ \zeta^ {\rm{id}} = (0,\infty)^{2+ n}$ that satisfy $\rm ($ $\textsf{T}_{\textsf{0}} - \textsf{T}_{\textsf{3}} )$. Moreover, we have $\mathcal{O}_ \pi^ {\rm{id}} = (0,\infty)^{2+ n}$ and

where $e_i^0$ denotes the energy of the $i$th species at zero temperature $e_k^0 = e_k^ \mathrm{st} - \int_0^{T^ \mathrm{st}} \!\! c_{ {{{\rm{v}}}} k}^ {\rm{id}}(\theta)\, d\theta$, $k\in {\mathfrak S}$.

Proof. The proof is straightforward ^[8]. □

3.   Nonhomogeneous Gibbsian thermodynamics

Various thermodynamic formalisms describing multicomponent fluids in terms of intensive variables may be found in the literature ^{[1,2,3,4,5,6,7,8,9]}. In nonuniform fluids, extensivity is indeed associated with volume integration and state variables can only be intensive variables like volumetric or mass densities ^[8]. Therefore, any set of $n+2$ intensive state variables—that correspond to independent extensive state variables after volume integration in spatially uniform flows—are now a priori dependent variables. Considering the species mass fractions to be independent leads to the homogeneous Gibbsian thermodynamics investigated in Section 2, and then the constraint $\langle \mathrm{y}, {1\kern -2.7pt \mathrm{I}}\rangle = 1$ is not imposed a priori but results from the conservation equations and boundary conditions ^[8].

An alternative method is to reduce the number of intensive variables by eliminating one of them. It only remains then a set of $n+1$ independent intensive state variables, and this is the formalism investigated in this section. Negative definite entropy Hessians may then be obtained when thermodynamic stability holds. The variable $\mathcal{Z} = (T, \rho_1, \ldots, \rho_ n)^t$ is used for convenience where $\rho_1, \ldots, \rho_ n$ denote the species masses per unit volume, also termed species partial densities as well as $\mathcal{X} = (\mathcal{E}, \rho_1, \ldots, \rho_ n)^t$. We also discuss the use of the thermodynamic—part of the—entropic variable $\mathcal{V}$ and of a normal variable $\mathcal{W}$. The use of a normal variable $\mathcal{W}$ may still be allowed when the matrix $\Lambda$ has at most a single negative eigenvalue, a situation generally met in practice, even though the use of the entropic variable $\mathcal{V}$ is then prohibited. We denote with calligraphic letters the properties per unit volume, and, in particular, $\mathcal{S}$ is the entropy per unit volume and $\mathcal{E}$ the energy per unit volume. In order to avoid confusion, we also denote in this section by $\mathcal{P}$ the pressure as function of $\mathcal{Z}$. We denote by ${\mathfrak d}^{}\mspace{-2mu}$ the derivation operator with respect to the variable $\mathcal{Z}$ and by $\widetilde{\mathfrak d}^{}\mspace{-2mu}$ the derivation operator with respect to the variable $\mathcal{X}$. We remind that $\gamma\in {\mathbb N}$, $\gamma\geqslant2$, is the regularity class of thermodynamic functions.

3.1.   From homogeneous to nonhomogeneous thermodynamics

We assume that a thermodynamics is available in terms of the natural variable $\zeta = (T, v, \mathrm{y}_1, \ldots, \mathrm{y}_ n)^t$ and investigate the definition and properties of thermodynamic functions in terms of the variable $\mathcal{Z} = (T, \rho_1, \ldots, \rho_ n)^t$ involving the species partial densities.

Theorem 3.1. Assume that $e$, $p$, and $s$ satisfy Properties $( \textsf{T}_{\textsf{0}}- \textsf{T}_{\textsf{2}} )$. Let $\mathcal{E}$, $\mathcal{P}$, and $\mathcal{S}$, be given by

and defined over the open set $\mathcal{O}_ \mathcal{Z} \subset (0, \infty)^{1+ n}$

Then $\mathcal{E}$, $\mathcal{P}$, and $\mathcal{S}$ satisfy Properties $( {\mathcal {T}}_{\textsf{0}} - {\mathcal {T}}_{\textsf{2}} )$.

$( {\mathcal {T}}_{\textsf{0}})$ The functions $\mathcal{E}$, $\mathcal{P}$, and $\mathcal{S}$ are $C^ \gamma$ over the nonempty, connected, open set $\mathcal{O}_ \mathcal{Z} \subset (0,\infty)^{1+ n}$.

$( {\mathcal {T}}_{\textsf{1}})$ For any $\mathcal{Z}\in \mathcal{O}_ \mathcal{Z}$, defining $g_k = {\mathfrak d}^{}\mspace{-2mu}_{ \rho_k} \mathcal{E} - T \, {\mathfrak d}^{}\mspace{-2mu}_{ \rho_k} \mathcal{S}$, $k\in {\mathfrak S}$, we have Gibbs' volumetric relation

with the constraint $\sum_{k\in {\mathfrak S}} \rho_k g_k = \mathcal{E} + \mathcal{P} - T \mathcal{S}$.

$( {\mathcal {T}}_{\textsf{2}})$ For any $(T, \mathrm{y}_1, \ldots, \mathrm{y}_ n) \in(0, \infty)^{1+ n}$, such that $\sum_{i\in {\mathfrak S}} \mathrm{y}_i = 1$, there exists a threshold density $\rho_{{{\rm{m}}}} > 0$ such that $(T, \rho \mathrm{y}_1, \ldots, \rho \mathrm{y}_ n)^t \in \mathcal{O}_ \mathcal{Z}$ for $0 < \rho < \rho_{{{\rm{m}}}}$ and

Proof. We define $\rho = \sum_{k\in {\mathfrak S}} \rho_k$, $v = 1/ \rho$, and $\mathrm{y}_k = \rho_k/ \rho$, $k\in {\mathfrak S}$, and all functions per unit mass are evaluated at $\zeta = (T, v, \mathrm{y}_1, \ldots, \mathrm{y}_ n)^t$. The fact that $\mathcal{O}_ \mathcal{Z}$ is a nonempty and connected open set is a direct consequence of the continuity of $\mathcal{Z} \to \zeta$, and the smoothness of $\mathcal{S}$, $\mathcal{E}$, and $\mathcal{P}$ is straightforward.

From $( \textsf{T}_{\textsf{0}})$ and the 1-homogeneity of energy with respect to $( v, \mathrm{y}_1, \ldots, \mathrm{y}_ n)$, we obtain, after some algebra, that

and, similarly,

This implies, in particular, that ${\mathfrak d}^{}\mspace{-2mu}_T \mathcal{E} = \rho \partial^{}_T e$, ${\mathfrak d}^{}\mspace{-2mu}_{ \rho_k} \mathcal{E} = \partial^{}_{ \mathrm{y}_k} e$, ${\mathfrak d}^{}\mspace{-2mu}_T \mathcal{S} = \rho \partial^{}_T s$, and ${\mathfrak d}^{}\mspace{-2mu}_{ \rho_k} \mathcal{S} = \partial^{}_{ \mathrm{y}_k} s$ in such a way that $g_k = {\mathfrak d}^{}\mspace{-2mu}_{ \rho_k} \mathcal{E} - T \, {\mathfrak d}^{}\mspace{-2mu}_{ \rho_k} \mathcal{S} = \partial^{}_{ \mathrm{y}_k} e - T \, \partial^{}_{ \mathrm{y}_k} s$ for $k\in {\mathfrak S}$. Then, thanks to $( \textsf{T}_{\textsf{1}} )$, since $T \partial^{}_T s = \partial^{}_T e$ and $T \partial^{}_{ \mathrm{y}_k} s = \partial^{}_{ \mathrm{y}_k} e - g_k$, we obtain $T {\mathfrak d}^{}\mspace{-2mu}_T \mathcal{S} = {\mathfrak d}^{}\mspace{-2mu}_T \mathcal{E}$ and $T {\mathfrak d}^{}\mspace{-2mu}_{ \rho_k} \mathcal{S} = {\mathfrak d}^{}\mspace{-2mu}_{ \rho_k} \mathcal{E} - g_k$ for $k\in {\mathfrak S}$, and this now implies Gibbs' relation for volumetric densities (3.4). In addition, the relation $\sum_{k\in {\mathfrak S}} \rho_k g_k = \mathcal{E} + \mathcal{P} - T \mathcal{S}$ is a consequence of Lemma 2.2 and $( {\mathcal {T}}_{\textsf{1}} )$ is established.

Furthermore, $( {\mathcal {T}}_{\textsf{2}} )$ is a simple rewriting of $(\textsf{T}_{\textsf2} )$ with $e$, $p$, and $s$ evaluated at $T, v, \mathrm{y}_1,\ldots, \mathrm{y}_ n$, and the proof is complete. □

The functions $\mathcal{E}$, $\mathcal{P}$, and $\mathcal{S}$ are said to define a thermodynamics if Properties $( {\mathcal {T}}_{\textsf{0}}- {\mathcal {T}}_{\textsf{2}} )$ hold. Note that there are not anymore homogeneity properties in terms of the variable $\mathcal{Z}$ in $( {\mathcal {T}}_{\textsf{0}})$. There is also a constraint associated with Gibbs' relation (3.4) in terms of $\mathcal{S}$, $\mathcal{E}$, and $\rho_1,\ldots, \rho_ n$ since the variables are per unit volume. This constraint essentially replaces the differential relation $T \partial^{}_T s = \partial^{}_T e + p$ associated with volume differentiation. The free energy per unit volume $\mathcal{F}$ is the given by $\mathcal{F} = \mathcal{E} - T \mathcal{S}$ and the Gibbs function per unit volume by $\mathcal{G} = \mathcal{F} + p = \mathcal{E} + p - T \mathcal{S}$.

Noting that ${\mathfrak d}^{}\mspace{-2mu}_T \mathcal{E} = \rho \partial^{}_T e$, we now investigate a sufficient conditions similar to Property $( \textsf{T}_{\textsf{3}})$ and often met in practice in such way that $\mathcal{Z} \mapsto \mathcal{X}$ is a global diffeomorphism. In this situation, it is possible to make the change of variable $\mathcal{Z}\mapsto \mathcal{X}$ and to rewrite Gibbs relation with the $\mathcal{X}$ variable.

$( {\mathcal {T}}_{\textsf{3}})$ For any $\mathcal{Z}\in \mathcal{O}_ \mathcal{Z}$ ${\mathfrak d}^{}\mspace{-2mu}_T \mathcal{E}( \mathcal{Z})>0$, $\mathcal{O}_ \mathcal{Z}$ is increasing with temperature.

Proposition 3.2. Assuming that Properties $( \textsf{T}_{\textsf{0}})$ and $( \textsf{T}_{\textsf{3}})$ hold, $( {\mathcal {T}}_{\textsf{0}})$ and $( {\mathcal {T}}_{\textsf{3}} )$ also hold. In this situation, the map $\mathcal{Z} \mapsto \mathcal{X}$ is a $C^ \gamma$ diffeomorphism from $\mathcal{O}_ \mathcal{Z}$ onto a nonempty, connected, open set $\mathcal{O}_ \mathcal{X}$ and the volumetric Gibbs relation may be written:

Proof. The proof is similar to that of the homogeneous case. □

We now investigate thermodynamic stability properties in the framework of the new thermodynamic formalism and note that the matrix $\Lambda$ may be rewritten as

since the Gibbs functions $g_l(\zeta)$ are 0-homogeneous. We have the following equivalence result concerning stability.

Proposition 3.3. Assume that $( {\mathcal {T}}_{\textsf{0}}-{\mathcal {T}}_{\textsf{2}})$ are satisfied and that the map $\mathcal{Z} \mapsto \mathcal{X}$ is a $C^ \gamma$ local diffeomorphism from an open set $\mathcal{O}_ \mathcal{Z}'$ onto an open set $\mathcal{O}_ \mathcal{X}'$. Then, for any $\mathcal{Z}\in \mathcal{O}_ \mathcal{Z}'$, the following statements are equivalent:

$(i)$ $\widetilde{\mathfrak d}^2_{ \mathcal{X} \mathcal{X}} \mathcal{S}$ is negative definite.

$(ii)$ ${\mathfrak d}^\mspace{-2mu}_T \mathcal{E} >0$ and $\Lambda$ is positive definite.

Proof. From Gibbs relation (3.8) in the $\mathcal{X}$ variable, we obtain that $\widetilde{\mathfrak d}^{}\mspace{-2mu}_ \mathcal{E} \mathcal{S} = \frac{1}{T}$, $\widetilde{\mathfrak d}^{}\mspace{-2mu}_{ \rho_k} \mathcal{S} = -\frac{ g_k}{T} = -\frac{ g_k}{T}$, and this implies the compatibility relations $\widetilde{\mathfrak d}^{}\mspace{-2mu}_ \mathcal{E} \bigl(\frac{- g_k}{T}\bigr) = \widetilde{\mathfrak d}^{}\mspace{-2mu}_{ \rho_k} \bigl(\frac{1}{T}\bigr)$. Moreover, for any function $\chi$ of $\mathcal{Z}$ and $\mathcal{X}$, we have the differential relations

We now evaluate the Hessian matrix of $\mathcal{S}$. We first have $\widetilde{\mathfrak d}^2_{ \mathcal{E} \mathcal{E}} \mathcal{S} = \widetilde{\mathfrak d}^{}\mspace{-2mu}_ \mathcal{E}(\frac{1}{T})$ so that $\widetilde{\mathfrak d}^2_{ \mathcal{E} \mathcal{E}} \mathcal{S} = - \widetilde{\mathfrak d}^{}\mspace{-2mu}_ \mathcal{E} T/T^2$. Similarly, we have $\widetilde{\mathfrak d}^2_{ \mathcal{E} \rho_k} \mathcal{S} = \widetilde{\mathfrak d}^{}\mspace{-2mu}_{ \rho_k}(\frac{1}{T})$ so that $\widetilde{\mathfrak d}^2_{ \mathcal{E} \rho_k} \mathcal{S} = - \widetilde{\mathfrak d}^{}\mspace{-2mu}_{ \rho_k} T/T^2$. Letting $\chi = T$ in (3.11), we get $\widetilde{\mathfrak d}^{}\mspace{-2mu}_{ \rho_k} T = - \widetilde{\mathfrak d}^{}\mspace{-2mu}_ \mathcal{E}^{}T \, {\mathfrak d}^{}\mspace{-2mu}_{ \rho_k} \mathcal{E}$, and since ${\mathfrak d}^{}\mspace{-2mu}_{ \rho_k} \mathcal{E} = \partial^{}_{ \mathrm{y}_k} e$, we have established that $\widetilde{\mathfrak d}^2_{ \mathcal{E} \rho_k} \mathcal{S} = \widetilde{\mathfrak d}^{}\mspace{-2mu}_ \mathcal{E} T\, \partial^{}_{ \mathrm{y}_k} e/T^2$. Moreover, with $\widetilde{\mathfrak d}^2_{ \rho_k \rho_l} \mathcal{S} = - \widetilde{\mathfrak d}^{}\mspace{-2mu}_{ \rho_k} (\frac{ g_l}{T}) = - \widetilde{\mathfrak d}^{}\mspace{-2mu}_{ \rho_k} (\frac{ g_l}{T})$ and (3.11), we deduce that $\widetilde{\mathfrak d}^2_{ \rho_k \rho_l} \mathcal{S} = - {\mathfrak d}^{}\mspace{-2mu}_{ \rho_k} (\frac{ g_l}{T}) - {\mathfrak d}^{}\mspace{-2mu}_T^{} (\frac{ g_l}{T}) \; \widetilde{\mathfrak d}^{}\mspace{-2mu}_{ \rho_k} T$ so that $\widetilde{\mathfrak d}^2_{ \rho_k \rho_l} \mathcal{S} = - \frac{ \partial^{}_{ \mathrm{y}_k} g_l}{ \rho\, T} - \frac{ \widetilde{\mathfrak d}^{}\mspace{-2mu}_ \mathcal{E} T\, \partial^{}_{ \mathrm{y}_k} \! e \; \partial^{}_{ \mathrm{y}_l} \! e}{T^2}$, $k, l\in {\mathfrak S}$. We have thus established that

Letting $\mathfrak{f}^ \mathcal{E} = (1, - \partial^{}_{ \mathrm{y}_1} e, \ldots, - \partial^{}_{ \mathrm{y}_1} e)^t$, we obtain that for any $\mathsf{X} = (\mathsf{X}_ \mathcal{E}, \mathsf{X}_1, \ldots, \mathsf{X}_ n)^t\in {\mathbb R}^{1+ n}$

where $\Lambda_{kl} = \partial^{}_{ \mathrm{y}_k} g_l/T$, $k, l\in {\mathfrak S}$ and $\mathsf{X}_ \mathrm{y} = (\mathsf{X}_1, \ldots, \mathsf{X}_ n)^t$.

The equivalence between $(i)$ et $(ii)$ is then a consequence of the quadratic relation (3.14). □

Remark 3.4. The functions $\mathcal{E}$, $\mathcal{P}$, and $\mathcal{S}$ may also be defined by $\mathcal{E}(T, \rho_1, \ldots, \rho_ n) = e \bigl(T, 1, \rho_1, \ldots, \rho_ n \bigr)$, $\mathcal{P}(T, \rho_1, \ldots, \rho_ n) = p \bigl(T, 1, \rho_1, \ldots, \rho_ n \bigr)$, and $\mathcal{S}(T, \rho_1, \ldots, \rho_ n) = s \bigl(T, 1, \rho_1, \ldots, \rho_ n \bigr)$, and the open set $\mathcal{O}_ \mathcal{Z}$ by $\mathcal{O}_ \mathcal{Z} = \{\; (T, \rho_1, \ldots, \rho_ n) \in (0, \infty)^{1+ n}; \ \bigl(T, 1, \rho_1, \ldots, \rho_ n \bigr)\in \mathcal{O}_ \zeta \; \}$ thanks to the homogeneity properties of $e$, $s$, $p$, and $\mathcal{O}_ \zeta$.

3.2.   From nonhomogeneous to homogeneous thermodynamics

We assume that a thermodynamics is available in terms of the variable $\mathcal{Z} = (T, \rho_1, \ldots, \rho_ n)^t$ and we investigate the definition and properties of thermodynamic functions in terms of the natural variable $\zeta = (T, v, \mathrm{y}_1, \ldots, \mathrm{y}_ n)^t$.

Theorem 3.5. Assume that $\mathcal{E}$, $\mathcal{P}$, and $\mathcal{S}$ satisfy Properties $( {\mathcal {T}}_{\textsf{0}}-{\mathcal {T}}_{\textsf{2}})$. Let $e$, $p$, and $s$ be given by

and defined over the open set

Then $e$, $p$, and $s$ are $C^ \gamma$ functions of the variable $\zeta = (T, v, \mathrm{y}_1, \ldots, \mathrm{y}_ n)^t$ over $\mathcal{O}_ \zeta\subset (0, \infty)^{2+ n}$ and satisfy Properties $(\textsf{T}_{\textsf{0}}-\textsf{T}_{\textsf{2}})$.

Proof. The fact that $\mathcal{O}_ \zeta$ is open, nonempty, and connected and that $e$, $p$, and $s$ are smooth is straightforward to establish. The homogeneity properties of $\mathcal{O}_ \zeta$, $e$, $p$, and $s$ are also a direct consequence of their definition (3.15)–(3.18).

In order to establish Gibbs relation (2.1), we note that

However, from Gibbs' relation in $( {\mathcal {T}}_{\textsf{2}})$, we obtain that $T {\mathfrak d}^{}\mspace{-2mu}_T \mathcal{S} = {\mathfrak d}^{}\mspace{-2mu}_T \mathcal{E}$. From the definition of $g_k = {\mathfrak d}^{}\mspace{-2mu}_{ \rho_k}\! \mathcal{E} - T {\mathfrak d}^{}\mspace{-2mu}_{ \rho_k}\! \mathcal{S}$, $k\in {\mathfrak S}$ and from the natural constraint $\sum_{k\in {\mathfrak S}} \rho_k g_k = \mathcal{E} + \mathcal{P} - T \mathcal{S}$, we obtain that

so in forming $T \mathrm{d} s - \mathrm{d} e$, we directly obtain $(\textsf{T}_{\textsf{1}})$.

Finally, $(\textsf{T}_{\textsf{2}})$ is a simple transformation of $( {\mathcal {T}}_{\textsf{2}})$ and the proof is complete. □

When thermal stability holds with respect to the thermodynamic in the variable $\mathcal{Z}$, it then holds with respect to the variable $\zeta$.

Proposition 3.6. Assuming Property $( {\mathcal {T}}_{\textsf{3}})$, Property $(\textsf{T}_{\textsf{3}})$ holds.

Proof. The proof is straightforward since ${\mathfrak d}^{}\mspace{-2mu}_T \mathcal{E} = \rho \partial^{}_T e$. □

3.3.   Entropic and normal variables

Entropic variables are naturally involved in symmetrization methods for systems of conservation laws ^{[10,11,12,13,14,15,16,17,18,19]}. They naturally appear in the expression of entropy production rates as well as in the definition of chemical equilibrium ^[5,8,9]. Normal variables, that split between hyperbolic and parabolic components are also involved in a priori estimates and existence theorems for hyperbolic-parabolic systems of conservation laws ^{[11,14,15,16,17,18,19]}.

The thermodynamic—part of the—entropic variable is defined by $\mathcal{V} = - ( {\mathfrak d}^{}\mspace{-2mu}_ \mathcal{X} \mathcal{S})^t$ and easily evaluated to be

The mathematical entropy $\sigma$ is defined by $\sigma = - \mathcal{S}$ in such a way that it is a convex function around stable states and this yields the minus sign in the definition of $\mathcal{V}$ ^{[10,11,12,13,14,15,16]}. The entropic variable $\mathcal{V} = ({\mathfrak d}^{}\mspace{-2mu}_ \mathcal{X} \sigma)^t$ is associated with symmetrization methods and

is notably symmetric positive definite around stable points. The mathematical entropy $\sigma$ is generally compatible with convective, diffusive, and chemical sources for fluid systems ^[8,18]. Entropic variables lead to symmetrized system of conservation laws that may be used to derive a priori estimates or to establish local or global existence results ^{[10,11,12,13,14,15,16,17,18,19]}.

A major difficulty associated with the entropic variable, however, is that the map $\mathcal{X}\mapsto \mathcal{V}$ cannot generally be one-to-one—and, thus, a global diffeomorphism—in the presence of phase changes. The equilibrium conditions ^[1] between two different phases $\mathcal{X}^\flat \ne \mathcal{X}^\sharp$ indeed implies that $\mathcal{V}^\flat = \mathcal{V}^\sharp$ and the map $\mathcal{X}\mapsto \mathcal{V}$ thus cannot be one-to-one. In other words, using the entropic variables $\mathcal{V}$ is not a good idea for multiphase fluids. In addition, entropic variables are used in the study of hyperbolic-parabolic systems modeling fluid mix hyperbolic and parabolic components. In order to solve these difficulties, it is mandatory to use normal variables ^[8,15].

We thus introduce the thermodynamic—part of a—normal variable in the form

and the following property is required in order to use this normal variable $\mathcal{W}$.

$( {\mathcal {T}}_{\textsf{4}})$ The map $\mathcal{Z} \mapsto \mathcal{W}$ is a $C^{ \gamma-1}$ diffeomorphism from $\mathcal{O}_ \mathcal{Z}$ onto an open set $\mathcal{O}_ \mathcal{W}$.

We first deduce important consequences of $({\mathcal {T}}_{\textsf{4}})$ and next investigate sufficient conditions that imply the property $({\mathcal {T}}_{\textsf{4}})$.

Proposition 3.7. Assume that $({\mathcal {T}}_{\textsf{0}}- {\mathcal {T}}_{\textsf{2}})$ and $({\mathcal {T}}_{\textsf{4}})$ hold. Then, the matrix $\Lambda$ is positive definite on the hyperplane ${1\kern -2.7pt \mathrm{I}}^\perp$, and the following direct sum holds:

Proof. Since the map $\mathcal{Z} \mapsto \mathcal{W}$ is a $C^{ \gamma-1}$ diffeomorphism, the Jacobian matrix ${\mathfrak d}^{}\mspace{-2mu}_ \mathcal{Z} \mathcal{W}$ is invertible on the open set $\mathcal{O}_ \mathcal{Z}$. Keeping in mind that $\mathcal{Z} = (T, \rho_1, \ldots, \rho_ n)^t$, the Jacobian matrix ${\mathfrak d}^{}\mspace{-2mu}_ \mathcal{Z} \mathcal{W}$ is of size $n + 1$ and in the form

where $0_{i, j}$ denotes a zero block with $i$ lines and $j$ columns, ${\ast}_{i, j}$ a nonzero block with $i$ lines and $j$ columns, and where the ${\mathcal J}$ block is a square matrix of dimension $n$.

The lines $2 \le j \le n$ of ${\mathcal J}$ have their coefficients ${\mathcal J}_{i, j}$ given by ${\mathcal J}_{i, j} = \partial_{\rho_j} (\frac{ g_i- g_1}{T})$ in such a way that $\rho {\mathcal J}_{i, j} = \Lambda_{i, j}-\Lambda_{1, j}$ for $2 \le i \le n$ and $1 \le j \le n$. Subtracting now the first line of ${\mathcal J}$ multiplied by $\partial_{\rho_1} (\frac{ g_i- g_1}{T})$ from the $i$th line of ${\mathcal J}$, and denoting by $\mathsf{e}_i$, $i\in \mathfrak{S}$, the canonical base vectors of ${\mathbb R}^ n$, we have

where the ${\mathcal J}'$ block is a square matrix of size $n-1$ and, by construction, $\det({\mathcal J}) = \det({\mathcal J}')$ in such a way that ${\mathcal J}'$ remains invertible. The submatrix ${\mathcal J}'$ of size $n-1$ has its coefficients given by $\rho {\mathcal J}'_{i, j} = \Lambda_{i, j} - \Lambda_{1, j} - \Lambda_{i, 1} + \Lambda_{1, 1}$, $2 \le i, j\le n$, and remains invertible over $\mathcal{O}_ \mathcal{Z}$. Letting then $\mathcal{v}_i = \mathsf{e}_i - \mathsf{e}_1$, for $2 \le i \le n$, we note that the vectors $\mathcal{v}_i$, for $2 \le i \le n$, are spanning the hyperplane ${1\kern -2.7pt \mathrm{I}}^\perp$ and that

since $(\Lambda \mathcal{v}_k)_i = \Lambda_{ik} - \Lambda_{i1} = \Lambda_{ki} - \Lambda_{1i} = {\mathcal J}_{ki}/\rho$. Therefore, the matrix ${\mathcal J}'$ of size $n-1$ may be written

and is thus a Gramm matrix up to the $\rho$ scalar factor that remains invertible over $\mathcal{O}_ \mathcal{Z}$. This implies that the quadratic form associated with $\Lambda$ remains definite on ${1\kern -2.7pt \mathrm{I}}^\perp$ over the open set $\mathcal{O}_ \mathcal{Z}$. Since this quadratic form is positive definite for small $\rho$, and since $\mathcal{O}_ \mathcal{Z}$ is connected, we deduce that $\Lambda$ remains positive definite on the hyperplane ${1\kern -2.7pt \mathrm{I}}^\perp$ over $\mathcal{O}_ \mathcal{Z}$. Since ${\mathcal J}$ remains nonsingular and since $(\Lambda \mathcal{v}_k)_i = \Lambda_{ik} - \Lambda_{i1} = \Lambda_{ki} - \Lambda_{1i} = \rho {\mathcal J}_{ki}$, we also directly obtain that the direct sum (3.21) holds. □

An important consequence of $({\mathcal {T}}_{\textsf{4}})$ is that the matrix $\Lambda$ has at most one isolated negative eigenvalue. That is, when $({\mathcal {T}}_{\textsf{4}})$ holds, however unstable the thermodynamic state of a mixture is, there is at most one negative eigenvalue, and we also have the following corollary.

Corollary 3.8. Assume that Properties $({\mathcal {T}}_{\textsf{0}}- {\mathcal {T}}_{\textsf{4}})$ hold. Then, thermodynamic stability holds if, and only if, $\det\Lambda>0$.

The property that there is at most one negative eigenvalue for the matrix $\Lambda$ has been verified numerically for various mixtures associated with hydrogen combustion as well as for H$_2$-N$_2$ mixtures, O$_2$-H$_2$O mixtures, and C$_2$H$_6$-N$_2$ mixtures. For such binary mixtures, the property $\Lambda_{2,2} + \Lambda_{1,1} - \Lambda_{1,2} - \Lambda_{2,1} >0$ has also been verified numerically. We are now interested in the converse property of $({\mathcal {T}}_{\textsf{4}})$, assuming that the slices of $\mathcal{O}_ \mathcal{Z}$ along ${1\kern -2.7pt \mathrm{I}}^\perp$ are convex.

Proposition 3.9. Assume that $({\mathcal {T}}_{\textsf{0}}- {\mathcal {T}}_{\textsf{2}})$ hold, that the matrix $\Lambda$ is positive definite on the hyperplane ${1\kern -2.7pt \mathrm{I}}^\perp$, that for any $\mathcal{Z}\in \mathcal{O}_ \mathcal{Z}$ the set $\mathcal{O}_ \mathcal{Z}\cap( \mathcal{Z} + {1\kern -2.7pt \mathrm{I}}^\perp)$ is convex, and that the direct sum (3.21) holds. Then, Property $({\mathcal {T}}_{\textsf{4}})$ holds.

Proof. Since $\Lambda$ is positive definite on the hyperplane ${1\kern -2.7pt \mathrm{I}}^\perp$, and since we have the direct sum (3.21), we conclude that the Jacobian matrix ${\mathfrak d}^{}\mspace{-2mu}_ \mathcal{Z} \mathcal{W}$ is invertible. The map $\mathcal{Z} \mapsto \mathcal{W}$ is thus a local $C^{ \gamma-1}$ diffeomorphism by the local inversion theorem, and we only have to establish that this map is one-to-one.

Assuming that $\mathcal{Z}^\flat$ and $\mathcal{Z}^\sharp$ are such that $\mathcal{W}^\flat = \mathcal{W}^\sharp$, we first deduce from (3.20) that $T^\flat = T^\sharp$ and $\rho^\flat = \rho^\sharp$ so that $\mathcal{Z}^\sharp - \mathcal{Z}^\flat \in {1\kern -2.7pt \mathrm{I}}^\perp$. Since $\mathcal{O}_ \mathcal{Z}\cap( \mathcal{Z}^\flat + {1\kern -2.7pt \mathrm{I}}^\perp)$ is convex, we further obtain that the whole segment $\mathcal{Z}^\flat + \alpha( \mathcal{Z}^\sharp- \mathcal{Z}^\flat)$, $0 \le \alpha \le 1$, lies in $\mathcal{O}_ \mathcal{Z}$ since $\mathcal{Z}^\flat\in \mathcal{O}_ \mathcal{Z}$ and $\mathcal{Z}^\sharp- \mathcal{Z}^\flat \in {1\kern -2.7pt \mathrm{I}}^\perp$. We may thus write that

and multiplying scalarly by $\mathcal{Z}^\sharp- \mathcal{Z}^\flat$, we obtain that

and this yields

Letting $\varrho = (\rho_1, \ldots, \rho_ n)^t$, and using that $T^\flat = T^\sharp$ and $\rho^\flat = \rho^\sharp$, we then deduce that

and this implies that $\varrho^\sharp- \varrho^\flat = 0$ since $\varrho^\sharp- \varrho^\flat \in {1\kern -2.7pt \mathrm{I}}^\perp$ and $\Lambda$ is positive definite over ${1\kern -2.7pt \mathrm{I}}^\perp$, and the proof is complete. □

3.4.   Ideal gas mixtures

We discuss for completeness the case of ideal gas mixtures in terms of $\mathcal{Z} = (T, \rho_1, \ldots, \rho_ n)^t$ ^[8]. The pressure law is given by

and the energy per unit volume reads

where $\mathcal{E}_k^ {\rm{id}}$ denotes the internal energy per unit mass of the $k$th species, $e_k^ \mathrm{st}$ the formation energy of the $k$th species per unit mass at the standard temperature $T^ \mathrm{st}$, and $c_{ {{{\rm{v}}}} k}^ {\rm{id}}$ the specific heat per unit mass of the $k$ species. Similarly, the entropy per unit volume is given by

where $s_k^ \mathrm{st}$ is the formation entropy of the $k$th species per unit mass at the standard temperature $T^ \mathrm{st}$ and standard pressure $p^ \mathrm{st}$, and $\gamma^ \mathrm{st} = p^ \mathrm{st}/ R T^ \mathrm{st}$ is the standard concentration. The enthalpy $\mathcal{H}^ {\rm{id}} = \mathcal{E}^ {\rm{id}} + \mathcal{P}^ {\rm{id}}$, Gibbs function $\mathcal{G}^ {\rm{id}} = \mathcal{E}^ {\rm{id}} + \mathcal{P}^ {\rm{id}} - T \mathcal{S}^ {\rm{id}}$, and free energy $\mathcal{F}^ {\rm{id}} = \mathcal{E}^ {\rm{id}} -T \mathcal{S}^ {\rm{id}}$ are then easily evaluated. The assumptions required for ideal gas mixtures in this volumetric framework are yet Property $(\textsf{ID})$.

4.   Thermodynamics from equations of state

Thermodynamics should preferably be deduced from molecular theories like the kinetic theory of dense gases or statistical mechanics ^[55,56]. Various relevant expressions have been obtained for mixtures of dense gases but are not of the most practical form. Molecular dynamics ^[58], Lattice Boltzmann methods ^[59], and experimental measurements may also be used to study thermodynamic properties of high pressure fluid mixtures. In order to bridge between experimental measurements or molecular simulations and macroscopic fluid thermodynamics, it has been found appropriate to introduce accurate approximate pressure law and to build the corresponding thermodynamics ^{[62,63,64,65,66,67,68,69,70,71,72,73,74,75]}.

In this section we study when a pressure law may be used for constructing a Gibbsian thermodynamics compatible with that of ideal gas mixtures at low density. To this purpose, natural conditions are introduced on the pressure law $p( \zeta)$ and on the open set $\mathcal{O}_ \zeta$ where it is defined. These conditions do not involve anymore stability-related properties associated with supercritical fluids as in previous work ^[9]. The structural conditions are essentially 0-homogeneity of pressure as well as a second order behavior $p - p^ {\rm{id}} = O( \rho^2)$ at low densities $\rho$, where $p^ {\rm{id}}$ is the ideal gas pressure. We first determine the energy and entropy per unit mass and then show that they constitute a thermodynamics. The matching with ideal gases is due to Van der Waals and Gillespie ^[62] according to Beattie ^[64]. We also study more closely rescaled versions of the stability matrix $\Lambda$ that determine thermodynamically unstable states.

4.1.   Identifying the thermodynamic functions

We discuss in this section the energy, entropy, and species Gibbs functions corresponding to a given pressure law and compatibility at small density with that of ideal gas mixtures. We consider a pressure law $p( \zeta)$ as function of $\zeta = (T, v, \mathrm{y}_1,\ldots, \mathrm{y}_ n)^t$ defined on a connected open set $\mathcal{O}_ \zeta\subset(0,\infty)^{ n+2}$ and define the pressure corrector $\phi( \zeta) = p( \zeta) - p^ {\rm{id}}( \zeta)$ so that

The ideal gas pressure $p^ {\rm{id}}$ has been introduced in Section 2 and discussed in Sections 2.4 and 3.4. The open set $\mathcal{O}_ \zeta$ is said to be increasing with volume if for any $(T, v, \mathrm{y}_1, \ldots, \mathrm{y}_ n)^t \in \mathcal{O}_ \zeta$ whenever $v' > v$, we have $(T, v', \mathrm{y}_1, \ldots, \mathrm{y}_ n)^t \in \mathcal{O}_ \zeta$. We introduce the following properties concerning the pressure law (4.1).

$( {\textsf {P}}_{\textsf{0}})$ The pressure $p$ is a $C^{ \gamma+1}$ function of $\zeta$ over the nonempty, connected, open set $\mathcal{O}_ \zeta\subset(0, \infty)^{ n+2}$. For any $(\lambda, \zeta)\in(0, \infty)\times \mathcal{O}_ \zeta$, we have $\zeta_\lambda^{}\in \mathcal{O}_ \zeta$ and $p(\zeta_\lambda) = p(\zeta)$.

$( {\textsf {P}}_{\textsf{1}})$ The open set $\mathcal{O}_ \zeta\subset(0, \infty)^{ n+2}$ is increasing with volume.

$( {\textsf {P}}_{\textsf{2}})$ There exists $\underline{ v} > 0$ such that for any $(T, \mathrm{y}_1, \ldots, \mathrm{y}_ n)^t \in(0, \infty)^{1+ n}$ with $\sum_{i\in {\mathfrak S}} \mathrm{y}_i = 1$, and any $v\geqslant \underline{ v}$, we have $\zeta = (T, v, \mathrm{y}_1, \ldots, \mathrm{y}_ n)^t\in \mathcal{O}_ \zeta$. The corrector $\phi$ and its partial derivatives $\partial^ \beta_ \zeta\phi$ of order $\vert \beta\vert \leqslant 1+ \gamma$ have a smooth extension to the set $T > 0$, $v\geqslant \underline{ v}$, $\mathrm{y}_i\geqslant0$, $\langle \mathrm{y}, {1\kern -2.7pt \mathrm{I}} \rangle = 1$. There exist positive continuous functions of temperature $\textsf{c}(\gamma, T)$ depending on $\gamma$, such that for any $T > 0$, $v\geqslant \underline{ v}$, and $\mathrm{y}_i\geqslant0$, $i\in {\mathfrak S}$, $\langle \mathrm{y}, {1\kern -2.7pt \mathrm{I}} \rangle = 1$, we have

We also introduce a sufficient condition that will ensure that the map $\zeta\mapsto \xi$ is one-to-one.

$( {\textsf {P}}_{\textsf{3}})$ For any $\zeta\in \mathcal{O}_ \zeta$, $\partial^2_T\phi(\zeta) \le 0$, $\mathcal{O}_ \zeta$ is increasing with temperature.

The 0-homogeneity property of pressure $(\textsf{P}_{\textsf{0}})$ is natural since we seek to obtain a thermodynamics such that $(\textsf{T}_{\textsf{0}})$ holds. The structure of the open set $\mathcal{O}_ \zeta$ with $(\textsf{P}_{\textsf{1}})$ is also required in order to allow integrations with respect to $v$ and to write the compatibility with ideal gases. In $(\textsf{P}_{\textsf{2}})$, we have used the traditional notation $\partial^ \beta = \partial^{ \beta_T}_T \partial^{ \beta_ v}_ v \partial^{ \beta_1}_{ \mathrm{y}_1} \cdots \partial^{ \beta_ n}_{ \mathrm{y}_ n}$ where $\beta = ( \beta_T, \beta_ v, \beta_1, \ldots, \beta_ n) \in {\mathbb N}^{2+ n}$ is a multiindex. The condition $(\textsf{P}_{\textsf{2}})$ guarantees that the corrector $\phi$ and its partial derivatives of order lower or equal to $1+ \gamma$ decrease at least quadratically with respect to $v$ as $v \to \infty$, insuring convergence of various integrals. The second order estimates (4.2) are naturally obtained by expanding the pressure law $p( \zeta)$ in a series of density $\rho = 1/ v$ in the neighborhood of the ideal gas limit $\rho\to0$ since the linear term is associated with the ideal gas contribution $p^ {\rm{id}}$. Property $(\textsf{P}_{\textsf{3}})$ is often met in practice and will guarantee that the thermal thermodynamic stability condition holds. We first obtain the energy and entropy associated with a pressure law satisfying $(\textsf{P}_{\textsf{0}}-\textsf{P}_{\textsf{2}})$ and next establish that they constitute a thermodynamics as in Definition 2.1.

Theorem 4.1. Let $p$ be a pressure law in the form (4.1) with $\phi$ satisfying $(\textsf{P}_{\textsf{0}}-\textsf{P}_{\textsf{2}})$ and assume that there exists a corresponding Gibbsian thermodynamics as in Definition 2.1. Then, the energy $e$, entropy $s$, and species Gibbs functions $g_k$, $k\in {\mathfrak S}$, are given by

The triplet function $e$, $p$, and $s$ defined on $\mathcal{O}_ \zeta$ is then a thermodynamics in the sense of Definition 2.1, and Properties $(\textsf{T}_{\textsf{0}}-\textsf{T}_{\textsf{2}})$ hold. Moreover, if $(\textsf{P}_{\textsf{3}})$ holds, then $(\textsf{T}_{\textsf{3}})$ also holds.

Proof. We first obtain from the compatibility relations (2.8) that $\partial^{}_ v e = T^2 \partial^{}_T\bigl(\frac{p}{T}\bigr)$ so that $\partial^{}_ v e = T^2 \partial^{}_T\bigl(\frac{\phi}{T}\bigr)$ since $\partial^{}_T\bigl({ p^ {\rm{id}}}/{T}\bigr) = 0$ from (2.17). Integrating with respect to $v'$ over $[ v, v'']$, keeping in mind that $\mathcal{O}_ \zeta$ in increasing with respect to volume, we obtain

Using the integrability properties (4.2) and $\lim_{ v''\to\infty} e(T, v'', \mathrm{y}_1, \ldots, \mathrm{y}_ n) = e^ {\rm{id}}(T, \mathrm{y}_1, \ldots, \mathrm{y}_ n)$ from $(\textsf{T}_{\textsf{2}})$, we may pass to the limit $v''\to\infty$ and we obtain (4.3).

Similarly, from $\partial^{}_ v(T s- e) = p$ and $\partial^{}_ v(T s^ {\rm{id}}- e^ {\rm{id}}) = p^ {\rm{id}}$, we deduce $\partial^{}_ v(T s- e) - \partial^{}_ v(T s^ {\rm{id}}- e^ {\rm{id}}) = \phi$. Integrating with respect to $v'$ over $[ v, v'']$ yields

Passing to the limit $v''\to\infty$, using (4.2) and $(\textsf{T}_{\textsf{2}})$, and noting that $\bigl( T s - T s^ {\rm{id}} - e + e^ {\rm{id}} \bigr) (T, v'', \mathrm{y}_1, \ldots, \mathrm{y}_ n)$ goes to zero as $v''\to\infty$ yields (4.4).

We now establish that $(\textsf{T}_{\textsf{0}}-\textsf{T}_{\textsf{2}})$ hold using (4.3) and (4.4). Using the regularity of $\phi$ and the estimates (4.2), we may differentiate (4.3) and (4.4) with respect to $T$, $v$, and $\mathrm{y}_k$, $k\in {\mathfrak S}$, and deduce that $e$ and $s$ are $C^ \gamma$ functions over $\mathcal{O}_ \zeta$. Deriving with respect to the mass fraction $\mathrm{y}_k$, we may calculate $g_k = \partial^{}_{ \mathrm{y}_k} e - T \partial^{}_{ \mathrm{y}_k} s$ and obtain (4.5), and it is further established that $g_k$, $k\in {\mathfrak S}$, are $C^{ \gamma-1}$ functions over $\mathcal{O}_ \zeta$.

We also note that if $\zeta\in \mathcal{O}_ \zeta$ and $\lambda>0$, then $\zeta_\lambda\in \mathcal{O}_ \zeta$ from $(\textsf{P}_{\textsf{0}})$. Taking into account the homogeneity properties of $e^ {\rm{id}}$, $s^ {\rm{id}}$, and $g^ {\rm{id}}_k$, $k\in {\mathfrak S}$, and from $(\textsf{P}_{\textsf{0}})$, we easily deduce that $e$ and $s$ are 1-homogeneous over $\mathcal{O}_ \zeta$ and that $g_k$, $k\in {\mathfrak S}$ are 0-homogeneous over $\mathcal{O}_ \zeta$. The corrector $\phi$ is indeed 0-homogeneous, so that $\partial^{}_T\phi$ and $T \partial^{}_T\phi - \phi$ are 0-homogeneous, and the corresponding integrals $\int_ v^\infty \partial^{}_T\phi\, d v'$ and $\int_ v^\infty (T \partial^{}\phi - \phi)\,d v'$ are then 1-homogeneous. Similarly, $\partial^{}_{ \mathrm{y}_k}\phi$ is (-1)-homogeneous so that $\int_ v^\infty \partial^{}_{ \mathrm{y}_k}\phi\,d v'$ is 0-homogeneous and property $(\textsf{T}_{\textsf{0}})$ holds.

In order to establish the Gibbs relations (2.1), we obtain from (4.3), (4.4), and the estimates (4.2) that

and since $\partial^{}_T e^ {\rm{id}} = T \partial^{}_T s^ {\rm{id}}$, we obtain that $\partial^{}_T e = T \partial^{}_T s$. Similarly, from (4.3), (4.4), and the estimates (4.2), we have

Using now $\partial^{}_ v e^ {\rm{id}} = T \partial^{}_ v s^ {\rm{id}} - p^ {\rm{id}}$, we deduce that $\partial^{}_ v e = T \partial^{}_ v s^ {\rm{id}} - p^ {\rm{id}} + T \partial^{}_T\phi - \phi$, and thus $\partial^{}_ v e = T \partial^{}_ v s - p$ which implies $(\textsf{T}_{\textsf{1}})$. Keeping in mind that $g_k = \partial^{}_{ \mathrm{y}_k} e - T \partial^{}_{ \mathrm{y}_k} s$ by definition, and we have established properties $(\textsf{T}_{\textsf{1}})$.

The property $(\textsf{T}_{\textsf{2}})$ is then guaranteed by construction of the thermodynamics (4.1), (4.3), and (4.4) and from the estimates (4.2) since $(T, v, \mathrm{y}_1,\ldots, \mathrm{y}_ n)^t$ belongs to $\mathcal{O}_ \zeta$ for $v$ large enough. The fact that $(\textsf{P}_{\textsf{3}})$ implies $(\textsf{T}_{\textsf{3}})$ is finally a consequence of (4.3) since

$\partial_T e^ {\rm{id}} = \sum_{k\in {\mathfrak S}} \mathrm{y}_k c_{ {{{\rm{v}}}} k}^ {\rm{id}} > 0$ and $\partial^2_T\phi \le 0$ from $(\textsf{P}_{\textsf{3}})$, so that $\partial_T e>0$ and $(\textsf{T}_{\textsf{3}})$ is established. □

This construction of the thermodynamics appears to be more general and natural than that devoted to supercritical states where it was mixed with stability properties ^[9]. The main tools are the fluid equation of state, Gibbs' relation that yields the partial derivatives $\partial^{}_ v e = T^2 \partial^{}_T(p/T)$ and $\partial^{}_ v(T s- e) = p$, and the ideal mixture low-density limit that then plays the role of a boundary condition as $v\to\infty$. Note that both differential relations $\partial^{}_ v e = T^2 \partial^{}_T(p/T)$ and $\partial^{}_ v(T s- e) = p$ are consequences of $T \partial^{}_T s = \partial^{}_T e$ and $T \partial^{}_ v s = \partial^{}_ v e + p$ that are directly obtained from Gibbs relation.

4.2.   Rescaled stability matrices

From Theorem 4.1 and Proposition 2.6, the thermodynamic stability domain where the mixture is locally stable is the domain where $\partial^{}_T e>0$ and $\Lambda$ is positive definite. The condition $\partial^{}_T e>0$ is the thermal stability condition whereas the positive definiteness of the matrix $\Lambda$ includes both the mechanical and chemical stability conditions. Using (2.4) and (4.5), the matrix $\Lambda = \partial^2_{ \mathrm{y} \mathrm{y}} f/T$ is easily shown to be

In the low density limit $v\to\infty$, the specific heat $\partial^{}_T e = c_{ {{{\rm{v}}}}}$ tends to $c_{ {{{\rm{v}}}}}^ {\rm{id}}$, which is positive from $(\textsf{ID})$, and the matrix $\Lambda$ goes to $\lim_{v\to\infty} \Lambda = R \ {{{\rm{diag}}}}\bigl( \frac{1}{m_1 \mathrm{y}_1},\ldots, \frac{1}{m_{ n} \mathrm{y}_{ n}} \bigr)$, which is positive definite.

We now establish some results about the stability matrix $\Lambda$ that are convenient in order to locate the stability domain, especially as some mass fractions goes to zero. The spectrum of the matrix $\Lambda$ is related to that of matrices $\overline\Lambda$, $\widehat\Lambda_ \mathsf{r}$, and $\widehat\Lambda_ \mathsf{l}$, that behave more regularly for zero mass fractions. For any diagonalizable matrix $A$, we denote by $\mathsf{d}_{}^+(A)$, $\mathsf{d}_{}^0(A)$, and $\mathsf{d}_{}^-(A)$ the number of positive, zero and negative eigenvalues of $A$, respectively.

Proposition 4.2. Assume that Properties $(\textsf{P}_{\textsf{0}}-\textsf{P}_{\textsf{2}})$ hold. For any $\zeta \in \mathcal{O}_ \zeta$, let $\Lambda$ be defined as in (2.4) and $\varPi$ be the diagonal matrix

The matrices $\overline\Lambda$, $\widehat\Lambda_ \mathsf{r}$, and $\widehat\Lambda_ \mathsf{l}$, defined by

are then $C^{ \gamma-1}$ functions of $\zeta\in \mathcal{O}_ \zeta$. Furthermore, if we assume that $v > \underline{ v}$, then $\overline\Lambda$ admits a continuous extension to nonnegative mass fractions whereas $\widehat\Lambda_ \mathsf{r}$ and $\widehat\Lambda_ \mathsf{l}$ admit $C^{ \gamma-1}$ extensions to nonnegative mass fractions. Finally, for any $\zeta\in \mathcal{O}_ \zeta$, we have $\mathsf{d}_{}^+(\Lambda) = \mathsf{d}_{}^+(\overline\Lambda)$, $\mathsf{d}_{}^0(\Lambda) = \mathsf{d}_{}^0(\overline\Lambda)$, $\mathsf{d}_{}^-(\Lambda) = \mathsf{d}_{}^-(\overline\Lambda)$, and the matrices $\overline\Lambda$, $\widehat\Lambda^ \mathsf{r}$, and $\widehat\Lambda_ \mathsf{l}$ have the same spectrum.

Proof. The proof is presented in ^[9]. □

We deduce from Proposition 4.2 that $\mathsf{d}_{}^+(\Lambda) = \mathsf{d}_{}^+(\overline\Lambda) = \mathsf{d}_{}^+(\widehat\Lambda_ \mathsf{r}) = \mathsf{d}_{}^+(\widehat\Lambda_ \mathsf{l})$, $\mathsf{d}_{}^0(\Lambda) = \mathsf{d}_{}^0(\overline\Lambda) = \mathsf{d}_{}^0(\widehat\Lambda_ \mathsf{r}) = \mathsf{d}_{}^0(\widehat\Lambda_ \mathsf{l})$, and $\mathsf{d}_{}^-(\Lambda) = \mathsf{d}_{}^-(\overline\Lambda) = \mathsf{d}_{}^-(\widehat\Lambda_ \mathsf{r}) = \mathsf{d}_{}^-(\widehat\Lambda_ \mathsf{l})$, in such a way that any of the matrices $\Lambda$, $\overline\Lambda$, $\widehat\Lambda_ \mathsf{r}$, or $\widehat\Lambda_ \mathsf{l}$, may be used over $\mathcal{O}_ \zeta$ to determine the stability domain.

It is interesting to discuss the limiting situation where some mass fractions vanish. Assume for the sake of simplicity that the positive mass fractions are the first $n^+$ components of $\mathrm{y}$ with $1 \leqslant n^+ < n$. The general case is easily reduced to this situation by introducing permutation matrices ^[8]. We then have $\mathrm{y}_k > 0$ for $1\leqslant k \leqslant n^+$ and $\mathrm{y}_k = 0$ for $n^++1\leqslant k \leqslant n$. Let's denote by $n^0 = n- n^+$ the number of zero mass fractions, $\mathrm{y}^+$ the vector of positive mass fractions, and $\mathrm{y}^0$ the vector of zero mass fractions. The indexing set of positive mass fractions is also denoted by ${\mathfrak S}^+$ and that of zero mass fractions by ${\mathfrak S}^0$. The decomposition ${\mathbb R}^ n_{} = {\mathbb R}^{ n^+}{\times}\, \, {\mathbb R}^{ n^0}$ induces a partitionning of vectors and any $\mathsf{X}\in {\mathbb R}^ n$ may be written in the form $\mathsf{X} = (\mathsf{X}^+, \mathsf{X}^0)$, and we have, for instance, ${1\kern -2.7pt \mathrm{I}}^+\in {\mathbb R}^{ n^+}$ and ${1\kern -2.7pt \mathrm{I}}^+ = (1, \ldots, 1)^t$. This partitionning of vectors induces a partitionning of matrices, and for any matrix $A\in {\mathbb R}^{ n, n}$, we denote by $A^{++}$, $A^{+0}$, $A^{0+}$, $A^{00}$, the corresponding blocks so that $(A\, \mathsf{X})^+ = A^{++} \mathsf{X}^+ + A^{+0} \mathsf{X}^0$ and $(A\, \mathsf{X})^0 = A^{0+} \mathsf{X}^+ + A^{00} \mathsf{X}^0$. From (4.6), only the blocks $\Lambda^{++}$, $\Lambda^{+0}$, and $\Lambda^{0+}$, of the matrix $\Lambda$ are defined for nonnegative mass fractions whereas the matrices $\overline\Lambda$, $\widehat\Lambda_ \mathsf{r}$, and $\widehat\Lambda_ \mathsf{l}$, are well-defined for nonnegative mass fractions. Naturally, the $++$ blocks are associated with the ${\mathfrak S}^+$ sub-mixture, that is, would be obtained by solely considering the species subset indexed by ${\mathfrak S}^+$.

Proposition 4.3. Assume that Properties $( \textsf{P}_{\textsf{0}} - \textsf{P}_{\textsf{1}})$ hold. Assume that $T>0$, $v > \underline{ v}$, and that the mass fraction vector is in the form $\mathrm{y} = ( \mathrm{y}^+, \mathrm{y}^0)$, where $\mathrm{y}^+\in {\mathbb R}^{ n^+}$, $\mathrm{y}^0\in {\mathbb R}^{ n^0}$, $\mathrm{y}^+ = ( \mathrm{y}_1,\ldots, \mathrm{y}_{ n^+})^t$, $\mathrm{y}^0 = (0,\ldots,0)^t$, $y_i>0$, $1\leqslant i \leqslant n^+$, and $\langle \mathrm{y}, {1\kern -2.7pt \mathrm{I}}\rangle = \langle \mathrm{y}^+, {1\kern -2.7pt \mathrm{I}}^+\rangle = 1$. Then, we have the block decompositions

where ${\mathbb I}$ denotes the identity matrix of size $n^0$ and the symbol $0$ denotes various rectangular matrices with zero entries so that $\overline\Lambda^{+0}=0$, $\overline\Lambda^{0+}=0$, $\widehat\Lambda_ \mathsf{r}^{+0}=0$, and $\widehat\Lambda_ \mathsf{l}^{0+}=0$. Moreover, denoting by

the upper left block of the matrix $\varPi$ we have the relations $\overline\Lambda^{++} = \varPi^{++} \Lambda^{++} \varPi^{++}$, $\widehat\Lambda_ \mathsf{r}^{++} = \Lambda^{++} (\varPi^{++})^2$, $\widehat\Lambda_ \mathsf{r}^{0+} = \Lambda^{0+} (\varPi^{++})^2$, $\widehat\Lambda_ \mathsf{l}^{++} = (\varPi^{++})^2 \Lambda^{++}$, $\widehat\Lambda_ \mathsf{l}^{+0} = (\varPi^{++})^2 \Lambda^{+0}$. In particular, we have $\mathsf{d}_{}^+(\Lambda^{++}) = \mathsf{d}_{}^+(\overline\Lambda^{++})$, $\mathsf{d}_{}^0(\Lambda^{++}) = \mathsf{d}_{}^0(\overline\Lambda^{++})$, $\mathsf{d}_{}^-(\Lambda^{++}) = \mathsf{d}_{}^-(\overline\Lambda^{++})$, and the matrices $\overline\Lambda^{++}$, $\widehat\Lambda_ \mathsf{r}^{++}$, and $\widehat\Lambda_ \mathsf{l}^{++}$ have the same spectrum.

Proof. These properties are easily obtained from (4.6) and the definition of the rescaled matrices using block manipulations. □

A practical consequence of Proposition 4.3 is that the stability of multicomponent mixtures may conveniently be investigated with the help of the matrices $\widehat\Lambda_ \mathsf{r}$ or $\widehat\Lambda_ \mathsf{l}$ which, behave smoothly for nonnegative mass fractions. The whole stability domain may then be obtained over the set $\mathrm{y}_i\geqslant0$, $i\in {\mathfrak S}$, $\langle \mathrm{y}, {1\kern -2.7pt \mathrm{I}}\rangle = 1$, including automatically all mixture states with zero mass fractions. As a special case, we may consider the limiting situation of pure species states. Denoting by $\mathsf{e}^1, \ldots, \mathsf{e}^ n$ the canonical base vectors of ${\mathbb R}^ n$ and by $\zeta^i$ the state $\zeta^i = (T, v, \mathsf{e}^i)^t$ for $\mathrm{y} = \mathsf{e}^i$, the following proposition discusses the thermodynamic stability in the neighborhood of $\mathsf{e}^i$.

Proposition 4.4. Assume that for $\zeta = \zeta^i$ we have $\partial^{}_T e(\zeta^i) = c_{ {{{\rm{v}}}} i} > 0$ and $\partial^{}_ v p(\zeta^i) < 0$. Then, thermodynamic stability holds in the neighborhood of $\zeta^i$.

Proof. We refer the reader to ^[9]. □

The modified matrices $\widehat\Lambda_ \mathsf{r}$ and $\widehat\Lambda_ \mathsf{l}$ may thus be conveniently used to build the whole stability diagram depending on the mixture state $\zeta$.

5.   Thermodynamics from SRK

We first present in Section 5.1 the Soave-Redlich-Kwong equation of state and obtain in Section 5.2 the corresponding thermodynamics.

5.1.   SRK pressure law

The Soave-Redlich-Kwong pressure law is in the form ^[65,66]

where $p$ denotes the pressure, $R$ the perfect gas constant, $v$ the volume per unit mass, and $a$ and $b$ the attractive and repulsive parameters per unit mass. These parameters $a (T, \mathrm{y}_1, \ldots, \mathrm{y}_{ n})$ and $b (\mathrm{y}_1, \ldots, \mathrm{y}_{ n})$ may be evaluated by using the Van der Waals mixing rules:

where the species attractive parameters $\alpha_i = \alpha_i(T)$, $i\in \mathfrak{S}$, are functions of temperature and the species repulsive parameters $b_i$, $i\in \mathfrak{S}$, are constants. The validity of this equation of state (5.1) and of the mixing rules (5.2) have been thoroughly investigated by comparison with NIST data by Congiunti et al. ^[72] and with results of Monte Carlo simulations by Colonna and Silva ^[73] and Cañas-Marín et al. ^[74,75]. This pressure law has been used in high pressure combustion models by Meng and Yang ^[77], Ribert et al. ^[79], and Giovangigli et al. ^[80,81].

The following mathematical properties are assumed on the coefficients $\alpha_i(T) = \sqrt{ a_i(T)}$ and $b_i$, $i\in {\mathfrak S}$, where $\gamma\in {\mathbb N}$ and $\gamma\geqslant2$.

$(\textsf{SRK})$ For any $i\in {\mathfrak S}$, $\alpha_i\in C^0[0, \infty)\cap C^{1+ \gamma}(0, \infty)$, $\alpha_i(0) > 0$, $\alpha_i\geqslant0$, $\lim_{+\infty} \alpha_i = 0$, $\partial^{}_T \alpha_i \leqslant0$, and $\partial^2_{TT} \alpha_i \geqslant0$ over $(0, \infty)$, and the parameters $b_i$, $i\in {\mathfrak S}$, are positive constants.

From physical considerations, the coefficient $\alpha_i = \sqrt{ a_i}$ should be positive for small temperature, be a decreasing function of temperature ^[69,70], and go to zero as $T\to\infty$ since we must recover the ideal gas limit, and all these properties are in $(\textsf{SRK})$. The following proposition is easily deduced from $(\textsf{SRK})$.

Proposition 5.1. Assume that $(\textsf{SRK})$ holds. Then, $a = a (T, \mathrm{y}) \in C^0\bigl([0,\infty)^{1+ n}\bigr) \cap C^{1+ \gamma}\bigl((0,\infty)^{1+ n}\bigr)$, $a(0, \mathrm{y})>0$, $a\geqslant0$, $\lim_{T\to+\infty} a(T, \mathrm{y}) = 0$, $\partial^{}_T a \leqslant0$, and $\partial_{TT}^2 a \geqslant0$ over $(0,\infty){\times}(0,\infty)^ n$.

Proof. This is a consequence of (5.2) and $(\textsf{SRK})$. □

5.2.   SRK thermodynamic properties

The thermodynamics associated with the SRK pressure law from Theorem 4.1 is discussed in this section.

Proposition 5.2. The corrector $\phi$ corresponding to the pressure law (5.1) may be written as

Under assumptions $(\textsf{SRK})$, the corrector $\phi$ is such that $(\textsf{P}_{\textsf{0}} -\textsf{P}_{\textsf{2}})$ hold on $\mathcal{O}_ \zeta$ given by

with notably $\underline{ v} = 2 \max_{i\in {\mathfrak S}} b_i$.

Proof. The proof is presented in ^[9]. □

From Properties $(\textsf{P}_{\textsf{0}} -\textsf{P}_{\textsf{2}})$ and Proposition 4.1, we may now identify the thermodynamic functions. It is further possible to calculate the thermodynamic functions explicitly, and from (5.3), the mixture internal energy $e$ is found in the form

where $e_i^{ {\rm{id}}} = e_i^{ {\rm{id}}} (T)$ denotes the ideal gas specific energy of the $i$th species. The mixture entropy $s$ is similarly given by

where $s_i^{ {\rm{id}}\star} = s_i^{ {\rm{id}}\star} (T)$ denotes the ideal gas specific entropy of the $i$th species at pressure $p^ \mathrm{st}$. Finally, the matrix $\Lambda$ has its coefficients given by

We now establish that property $(\textsf{P}_{\textsf{3}})$ holds under $(\textsf{SRK})$ so that thermal stability automatically holds.

Proposition 5.3. Assuming $(\textsf{SRK})$, the corrector $\phi = p - p^ {\rm{id}}$ is such that $\partial^2_{TT}\phi\le0$ and $\mathcal{O}_ \zeta$ is increasing with respect to temperature $T$ so that $(\textsf{P}_{\textsf{3}})$ holds.

Proof. By deriving (5.1) with respect to temperature, we obtain $\partial^2_{TT}\phi = - { \partial^2_{TT} a}/{ v ( v+ b)}$. On the other hand, using (5.2) we have $\partial^2_{TT} a = \sum_{i,j\in {\mathfrak S}}$$\mathrm{y}_i \mathrm{y}_j$$( \alpha_i \partial^2_{TT} \alpha_j + 2 \partial^{}_T \alpha_i \partial^{}_T \alpha_j + \alpha_i \partial^2_{TT} \alpha_j)$, and since $\alpha_i \partial^2_{TT} \alpha_j\ge0$ and $\partial^{}_T \alpha_i \partial^{}_T \alpha_j \ge 0$, we obtain that $\partial^2_{TT} a \ge 0$ and $\partial^2_{TT}\phi \le0$ from $(\textsf{SRK})$. Moreover, the open set $\mathcal{O}_ \zeta$ defined by (5.4) is independent of $T$ and $(\textsf{P}_{\textsf{3}})$ is established. □

The attractive parameters $\alpha_i$, $i \in \mathfrak{S}$ involved in the SRK pressure law are taken in practice in a special form ^[80] discussed in the following proposition.

Proposition 5.4. Assume that $\alpha_i$ is in the form

where $\alpha_{ {{{\rm{c}}}}, i}$, ${\rm{s}}_{i}$, and $T_{ {{{\rm{c}}}}, i}$ are positive constants and $A\in C^{1+ \gamma}({\mathbb R})$. Further assume that $A'\geqslant0$, $A''\geqslant0$, $A(0) = 0$, $-1 \leqslant A$, and $\lim_{-\infty} A = -1$. Then, Property $(\textsf{SRK})$ holds.

Proof. Since $A$ is $C^{1+ \gamma}( {\mathbb R})$, it is easily checked that, for any $i\in {\mathfrak S}$, we have $\alpha_i \in C^0[0,\infty)\cap C^{1+ \gamma}(0,\infty)$ and for $T>0$ we have

From the assumption on $A$, we deduce that $\partial^{}_T \alpha_i(T) \leqslant0$ and $\partial_{TT}^2 \alpha_i(T)\geqslant0$ so that $\alpha_i$ is nonincreasing and convex. Moreover, we have $\alpha_i(0) = \alpha_{ {{{\rm{c}}}}, i} \bigl(1+ A({\rm{s}}_{i}) \bigr)$ and $A({\rm{s}}_{i})\geqslant0$ since $A$ is nondecreasing and $A(0) = 0$ so that $\alpha_i(0) > 0$. Finally, we also have $\lim_{T\to+\infty} \alpha_i(T) = 0$ since $\lim_{x\to-\infty} A(x) = -1$, and this completes the proof. □

The parameter $\alpha_{ {{{\rm{c}}}}, i}$ is the value of $\alpha_i$ at the critical temperature $T_{ {{{\rm{c}}}}, i}$ and ${\rm{s}}_{i}$ is associated with the geometry of the molecule ^[113]. As an example of function $A$, one may use $A(x) = x$ when $x\geqslant0$, and $A(x) = x/(1+x^4)^{1/4}$ when $x\leqslant0$. It is then checked that $A\in C^4({\mathbb R})$, $A'\geqslant0$, $A''\geqslant0$, $A(0) = 0$, $-1\leqslant A$, and $\lim_{-\infty} A = -1$, and the assumptions of Proposition 5.4 hold with $\gamma = 3$. The Soave coefficients correspond to $A(x) = \vert 1+x\vert - 1$ and neither satisfy property $(\textsf{SRK})$ nor the assumptions of Proposition 5.4. It is possible to truncate those coefficients by using $A(x) = \vert 1+x\vert^+ - 1$, where $y^+ = \max(y,0)$, but neither the Soave coefficients nor their truncated version are smooth enough for mixtures of gases since they introduce small jumps in the temperature derivative of the attractive factor $\alpha_i = \sqrt{ a_i}$ and thereby of the attractive parameter $a = \sum_{i,j\in {\mathfrak S}} \mathrm{y}_i \mathrm{y}_j \alpha_i \alpha_j$ at the crossing temperatures $T_{ {{{\rm{c}}}},i} ( 1+ {1}/{ {\rm{s}}_{i}})^2$, $i\in {\mathfrak S}$. Numerical spurious behavior has been observed because of these small discontinuities of temperature derivative of the attractive factor ^[80].

5.3.   SRK mechanical stability

We address the mechanical stability condition for the SRK pressure law. For the sake of conciseness, we introduce $r = R \sum_{i\in {\mathfrak S}} \mathrm{y}_i/m_i$ so that the equation of state is written in the form

Proposition 5.5. For any $\mathrm{y} \in(0, \infty)^ n$ and $v > b$, there exists a unique temperature $T^\star(v, \mathrm{y}) > 0$ such that $\partial^{}_ v p(T^\star, v, \mathrm{y}) = 0$. This temperature satisfies the nonlinear equation

where $a^\star = a(T^\star, \mathrm{y})$, $T^\star$ is a $C^ \gamma$ function of $v, \mathrm{y}_1, \ldots, \mathrm{y}_ n$, and we have the relations

Finally, for any fixed $\mathrm{y} \ (0, \infty)^ n$, we have $\lim_{ v\to b} T^\star = 0$ and $\lim_{ v\to\infty} T^\star = 0$.

Proof. The proof is presented in ^[9]. □

For a frozen mixture $\mathrm{y}$, the curve $v \to T^\star(v, \mathrm{y})$ is the mechanical stability curve. Inside the domain defined by $T^\star$, mechanical stability does not hold as established in the following corollary.

Corollary 5.6. We have the equivalence $\partial^{}_ v p(T, v, \mathrm{y}) > 0 \Longleftrightarrow T < T^\star(v, \mathrm{y})$, and since the derivative $\partial^{}_T p$ at fixed $v, \mathrm{y}_1, \ldots, \mathrm{y}_ n$ is always positive, defining $p^\star (v, \mathrm{y}) = p \bigl(T^\star(v, \mathrm{y}), v, \mathrm{y} \bigr)$, we also have the equivalence $\partial^{}_ v p(T, v, \mathrm{y}) > 0 \Longleftrightarrow p(T, v, \mathrm{y}) < p^\star(v, \mathrm{y})$.

The mechanical stability condition $\partial^{}_ v p(T, v, \mathrm{y}) < 0$ may equivalently be written $T > T^\star(v, \mathrm{y})$ or $p(T, v, \mathrm{y}) > p^\star(v, \mathrm{y})$. We now discuss the maximum value of $T^\star$ when $v$ is varying, and the corresponding point is a critical point with frozen mass fractions ^[1].

Proposition 5.7. For any $\mathrm{y}\in(0, \infty)^ n$, there exists a unique maximum positive value of $T^\star$ for $v\in(b, \infty)$ reached for $v = \theta_ {{{\rm{c}}}} b$, where $\theta_ {{{\rm{c}}}} = 1/(\sqrt[3]{2} -1)$. For any fixed $\mathrm{y}\in(0, \infty)^ n$, define the functions

keeping in mind that $b = \sum_{i\in {\mathfrak S}} \mathrm{y}_i b_i$. Then, $T^{\star\star}(\mathrm{y})$ and $p^{\star\star}(\mathrm{y})$ correspond to the critical temperature and pressure of the frozen mixture with mass fractions $\mathrm{y}$. Moreover, we have $p^\star (v, \mathrm{y}) \leqslant p^{\star\star}(\mathrm{y})$ so that $\partial^{}_ v p > 0$ implies $p(T, v, \mathrm{y}) < p^{\star\star}(\mathrm{y})$. Finally, since the set $\overline { \Sigma } = \{ \mathrm{y}\in[0, \infty)^ n; \ \langle \mathrm{y}, {1\kern -2.7pt \mathrm{I}}\rangle = 1\}$ is compact and $p \bigl(T^\star(\theta_ {{{\rm{c}}}} v, \mathrm{y}), \theta_ {{{\rm{c}}}} v, \mathrm{y} \bigr)$ is a smooth function of $\mathrm{y}$, there exists a maximum pressure $\max_{ \mathrm{y}\in\overline { \Sigma}} p^{\star\star}(\mathrm{y})$ above which the fluid is mechanically stable for all mixtures.

Proof. We refer the reader to ^[9]. □

Critical points of mixtures have a more complex behavior than that of frozen mixtures. The complex fluid phase diagrams have been classified by Van Konynenburg and Scott ^[112,113].

6.   Extended thermodynamics and cohesive fluids

We present the Van der Waals free energy, the extended thermodynamics, and the diffuse interface fluid equations ^{[20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51]}. These equations model fluids with cohesive internal forces and include a Korteweg type pressure tensor and a Dunn and Serrin type heat flux. The gradient squared terms in the extended free energy are associated with long range molecular interactions, and the simplified diffuse interface model is obtained by assuming that the capillary coefficients of species pairs are identical, as detailed in Appendix C. Cahn-Hilliard models with different capillarity coefficients ^[44], gradient flows obtained in the absence of convection ^[129,130], and numerical issues ^[131], as well as boundary condition issues lay out of the scope of the present work.

6.1.   Van der Waals gradient energy

The free energy per unit volume $\mathcal{F}^ {\rm{ex}}$ is of Van der Waals type:

where $\mathcal{F}^ {}$ denotes the classical or bulk free energy. The superscript ${}^ {\rm{ex}}$ is used to denote extended thermodynamic properties whereas classical bulk phase thermodynamic properties—that do not involve gradients—are denoted without superscripts. The free energy $\mathcal{F}^ {}$ only depends on the absolute temperature $T$ and on the species densities $\rho_1, \ldots, \rho_n$, whereas the gradient term ${\tfrac{1}{2}} \varkappa\vert \boldsymbol{\nabla}\mspace{-3mu}\rho\vert^2$ in $\mathcal{F}^ {\rm{ex}}$ represents an interfacial excess free energy. This gradient term may be interpreted from attractive long range molecular interactions and the capillary coefficient $\varkappa$ related to the species interaction potentials ^[25,44].

From the definition (6.1) and the classical relation $\mathrm{d} \mathcal{F}^ {} = - \mathcal{S}^ {} \mathrm{d} T + \sum_{i\in \mathfrak{S}} g_i^ {} \mathrm{d}\rho_i$, it is obtained that $\mathrm{d} \mathcal{F}^ {\rm{ex}} = - (\mathcal{S}^ {} - {\tfrac{1}{2}}\partial_T \varkappa\vert {\boldsymbol\nabla}\rho\vert^2) \mathrm{d} T + \sum_{i\in \mathfrak{S}} g_i^ {} \mathrm{d}\rho_i + \varkappa \boldsymbol{\nabla}\mspace{-3mu}\rho {\boldsymbol{\cdot}} \mathrm{d} \boldsymbol{\nabla}\mspace{-3mu}\rho$ where $\mathcal{S}^ {}$ denotes the bulk entropy per unit volume and $g_i^ {}$ the bulk Gibbs function of the $i$th species per unit mass. From the thermodynamic relations $\partial_T \mathcal{F}^ {\rm{ex}} = - \mathcal{S}^ {\rm{ex}}$ and $\partial_{\rho_i} \mathcal{F}^ {\rm{ex}} = g_i^ {\rm{ex}}$, $i\in \mathfrak{S}$, where $\mathcal{S}^ {\rm{ex}}$ denotes the extended entropy per unit volume and $g_i^ {\rm{ex}}$ the Gibbs function per unit mass of the $i$th species, it is then obtained that $\mathcal{S}^ {\rm{ex}} = \mathcal{S}^ {} - {\tfrac{1}{2}}\partial_T \varkappa\vert {\boldsymbol\nabla}\rho\vert^2$ and $g_i^ {\rm{ex}} = g_i^ {}$.

From the expression (6.1) of the free energy $\mathcal{F}^ {\rm{ex}}$ and the identity $\mathcal{S}^ {\rm{ex}} = \mathcal{S}^ {} - {\tfrac{1}{2}}\partial_T \varkappa\vert {\boldsymbol\nabla}\rho\vert^2$, the Gibbs function $\mathcal{G}^ {\rm{ex}}$ is found to be given by its standard value $\mathcal{G}^ {\rm{ex}} = \mathcal{G}^ {}$ whereas the energy per unit volume $\mathcal{E}^ {\rm{ex}}$ and the pressure $p^ {\rm{ex}}$ are given by

Letting $e^ {\rm{ex}}$, $s^ {\rm{ex}}$, $f^ {\rm{ex}}$, and $g^ {\rm{ex}}$ denote, respectively, the energy, entropy, free energy, and Gibbs function per unit mass, by $\rho = \sum_{i\in \mathfrak{S}}\rho_i$ the mass density, and $\mathrm{y}_i$ the mass fraction of the $i$th species, we have $\mathcal{E}^ {\rm{ex}} = \rho e^ {\rm{ex}}$, $\mathcal{S}^ {\rm{ex}} = \rho s^ {\rm{ex}}$, $\mathcal{F}^ {\rm{ex}} = \rho f^ {\rm{ex}}$, $\mathcal{G}^ {\rm{ex}} = \rho g^ {\rm{ex}}$, $\rho_i = \rho \mathrm{y}_i$, $i\in \mathfrak{S}$. The volumetric Gibbs relation for the extended entropy $\mathcal{S}^ {\rm{ex}}$ finally reads

where $\mathrm{d}$ is the differentiation operator.

6.2.   Cohesive fluids

The species mass, momentum, and energy conservation equations are in the form

where $\partial_t$ denotes the time derivative operator, $\boldsymbol{\nabla}$ the spatial differential operator, $\rho_i$ the partial density of the $i$th species, $\omega_i$ the molar production rate of the $i$th species, ${\boldsymbol{u}}$ the fluid velocity, $\boldsymbol{\mathcal{J}}_i$ the mass flux of the $i$th species, $m_i$ the molar mass of the $i$th species, $\omega_i$ the molar production rate of the $i$th species, $\boldsymbol{\mathcal{P}}$ the total pressure tensor, $\mathcal{E}^ {\rm{ex}}$ the internal energy per unit volume, and $\boldsymbol{\mathcal{Q}}$ the total heat flux.

The diffusive fluxes $\boldsymbol{\mathcal{J}}_i$, $i\in \mathfrak{S}$, and the chemical source terms $\omega_i$, $i\in \mathfrak{S}$, satisfy the mass conservation constraints $\sum_{i\in \mathfrak{S}} \boldsymbol{\mathcal{J}}_i = 0$ and $\sum_{i\in \mathfrak{S}} m_i \omega_i = 0$ so that by summing the species Eq (6.4), we recover the total mass conservation equation

A conservation equation for the internal energy $\mathcal{E}^ {\rm{ex}}$ may also be derived by subtracting the kinetic energy balance equation as in standard fluids.

The transport fluxes obtained with identical capillarity coefficients are in the form

with the dissipative fluxes given by

where $\mathfrak{v}$ denotes the volume viscosity, $\eta$ the shear viscosity, and $L_{ij}$, $i, j\in \mathfrak{S}\cup\{ \mathrm{e}\}$, the mass and heat transport coefficients as derived in Appendix C. The matrix of mass and heat transport coefficients $L$ defined by $L = (L_{ij})_{i, j\in \mathfrak{S}\cup\{ \mathrm{e}\}}$ is symmetric positive semi-definite with null space spanned by the vector $(1, \ldots, 1, 0)^t$ as for ordinary fluids ^[44,60].

The dissipative fluxes $\boldsymbol{\mathcal{J}}_i$, $i\in \mathfrak{S}$, and $\boldsymbol{q}$ are in their thermodynamic form in (6.11) and (6.12) and expressed directly in terms of the gradients of the chemical potentials $g_j^ {}/T$, $j\in \mathfrak{S}$, keeping in mind that $g_j^ {\rm{ex}} = g_j^ {}$, and of the usual thermodynamic thermal variable $-1/T$. In the pressure tensor $\boldsymbol{\mathcal{P}}$, the new contributions from the Korteweg tensor ^[22] are $- \rho \boldsymbol{\nabla} {\boldsymbol{\cdot}}(\varkappa \boldsymbol{\nabla}\mspace{-3mu}\rho) \boldsymbol{I}$ and $\varkappa \boldsymbol{\nabla}\mspace{-3mu}\rho {\otimes} \boldsymbol{\nabla}\mspace{-3mu}\rho$, and the pressure $p^ {\rm{ex}}$ also differs from $p^ {}$. The extra capillary contribution in the total heat flux is the Dunn and Serrin heat flux $\varkappa \rho \boldsymbol{\nabla} {\boldsymbol{\cdot}} {\boldsymbol{u}} \, \boldsymbol{\nabla}\mspace{-3mu}\rho$ ^[23]. These extra capillary terms do not produce entropy at variance with diffusion, chemistry, and thermal conductivity.

6.3.   Production rates

In this section, we describe for completeness the chemical production rates even though no chemical reactions will be considered in the numerical simulations. We consider an arbitrary complex chemical reaction mechanism in the form

where $\nu_{ij}^ {{{\rm{f}}}}$ and $\nu_{ij}^ {{{\rm{b}}}}$ are the forward and backward stoichiometric coefficients of the $i$th species in the $j$th reaction, $\mathfrak{M}_i$ the symbol of the $i$th molecule, ${\mathfrak R} = \{1, \ldots, n^{{{\rm{r}}}}\}$ the indexing set of chemical reactions, and $n^{{{\rm{r}}}}$ the number of chemical reactions. The production rates are given by

where $\nu_{ij} = \nu_{ij}^ {{{\rm{b}}}} - \nu_{ij}^ {{{\rm{f}}}}$ is the net stoichiometric coefficient of the $i$th species in the $j$th reaction and $\tau_j$ is the rate of progress of the $j$th reaction. We introduce the following notation for convenience

where $\mu_i = m_i g_i /RT$ denotes the reduced bulk chemical potential of the $i$th species. The Marcelin's reaction rates of progress given by statistical thermodynamics ^{[8,57,80,88,89,90,91]} then read

where $\mathcal{K}_j$ is the constant of the $j$th chemical reaction.

Using Gibbs relation and the conservation equations, the balance equation governing $\mathcal{S}^ {\rm{ex}}$ is then obtained in the form

Using the expression of the transport fluxes and the production rates the balance of entropy is finally

where we have let $g_ \mathrm{e} = -1$ for convenience, $\boldsymbol{\mathsf S}$ denotes the deviatoric part of the strain rate tensor $\boldsymbol{\mathsf S} = \boldsymbol{\nabla} {\boldsymbol{u}} + \boldsymbol{\nabla} {\boldsymbol{u}}^t - {\tfrac{2}{3}} \boldsymbol{\nabla} {\boldsymbol{\cdot}} {\boldsymbol{u}} \boldsymbol{I}$, and $\vert A\vert^2$ the Frobenius norm $\vert A\vert^2 = A:A = \sum_{ij} a_{ij}^2$ of a tensor $A = (a_{ij})$, so that entropy production appears as a sum of nonnegative terms.

6.4.   Augmented formulation

The equations governing cohesive fluids (6.4)–(6.6) may also be reformulated by using an extra unknown vector ${\boldsymbol{w}}$ representing locally the gradient of density ^{[82,83,84,85,86,87]}:

We may also consider that the spatial dimension $d$ has eventually been reduced with $1 \le d \le 3$. The governing equation for ${\boldsymbol{w}}$ is obtained by applying the differential operator ${\boldsymbol\nabla}$ to (6.7) and reads

where $\mathcal{D} = \{1, \ldots, d\}$ denotes the indexing set of spatial directions, $u_i$ the velocity in the $i$th spatial direction, and $\partial_i^{}$ the derivative in the $i$th spatial direction with ${\boldsymbol{u}} = (u_1, \ldots, u_ d)^t$ and ${\boldsymbol\nabla} = (\partial_1, \ldots, \partial_ d)^t$ in $d$ dimensions.

The augmented conservative $\mathsf{u}$ and natural $\mathsf{z}$ variables for multicomponent mixtures are then

and are vectors of ${\mathbb R}^ \mathsf{n}$ with $\mathsf{n} = n + 2 d + 1$. The vectors of ${\mathbb R}^ \mathsf{n}$ are decomposed as ${\mathbb R} ^ \mathsf{n} = {\mathbb R}^ n\times {\mathbb R}^ d\times {\mathbb R}^ d \times {\mathbb R}$ and the transposition operation is made by block for such column vectors for the sake of notational simplicity. The equations governing the multicomponent augmented system are given by

where ${\mathsf F}_{i}$ denotes the augmented convective flux in the $i$th direction, ${\mathsf F}_{i}^{ {\rm{c}}}$ the augmented cohesive flux in the $i$th direction, ${\mathsf F}_{i}^{ {\rm{d}}}$ the augmented dissipative flux in the $i$th direction, and $\Omega$ the chemistry source term. These fluxes are identified by using the Legendre transform of entropy as in the single species situation ^[87] and found to be

where $\mathcal{P}^{{\rm{d}}}_{ij}$, $i, j\in \mathcal{D}$ denote the components of the viscous tensor $\boldsymbol{\mathcal{P}}^{{\rm{d}}}\!\!\! = (\mathcal{P}^{{\rm{d}}}_{ij})_{i, j\in \mathcal{D}}$, $\boldsymbol{\mathcal{P}}^{{\rm{d}}}\!\!\!_i$, $i\in \mathcal{D}$, the vectors $\boldsymbol{\mathcal{P}}^{{\rm{d}}}\!\!\!_i = (\mathcal{P}^{{\rm{d}}}_{i1}, \ldots, \mathcal{P}^{{\rm{d}}}_{i d})^t$, $\boldsymbol{\mathsf{e}}_i$, $i\in \mathcal{D}$, the base vectors in ${\mathbb R}^ d$, ${\boldsymbol{w}} = (w_1, \ldots, w_ d)^t$, the spatial components of ${\boldsymbol{w}}$, $\boldsymbol{q} = (q_1, \ldots, q_ d)^t$, the spatial components of the dissipative heat flux $\boldsymbol{q}$, $\boldsymbol{\mathcal{J}}_i = (\mathcal{J}_{1i}, \ldots \mathcal{J}_{ d i})^t$, the spatial components of the mass flux $\boldsymbol{\mathcal{J}}_i$ of the $i$th species and $\Omega$ the source term. The resulting augmented system is easily found to be equivalent to the original system by using calculus identities in a similar way as for single species fluids and the cohesive flux ${\mathsf F}_{i}^{ {\rm{c}}}$ is similar to that of single species fluids ^[87]. It is interesting to note that the relavant pressure is then $p^ {\rm{ex}} + \varkappa\vert {\boldsymbol{w}}\vert^2 = p + {\tfrac{1}{2}} \varkappa\vert {\boldsymbol{w}}\vert^2$.

As for single species fluids, the cohesive ${\mathsf F}_{i}^{ {\rm{c}}}$ and dissipative ${\mathsf F}_{i}^{ {\rm{d}}}$ fluxes are linear in the gradients, and may thus be written:

where ${\mathsf B}^{ {\rm{d}}}_{ij}$ and ${\mathsf B}^{ {\rm{c}}}_{ij}$ are uniquely defined matrices. We may further introduce the Jacobian ${\mathsf A}_{i} = \partial_ \mathsf{u}{\mathsf F}_{i}$ of the convective flux ${\mathsf F}_{i}$ and the governing Eq (6.19) may finally be written as an augmented quasilinear second order system of partial differential equations

The coefficient matrices ${\mathsf A}_{i} = \partial_ \mathsf{u}{\mathsf F}_{i}$, ${\mathsf B}^{ {\rm{d}}}_{ij}$, and ${\mathsf B}^{ {\rm{c}}}_{ij}$, for $i, j\in \mathcal{D}$, have at least regularity $C^{ \gamma-1}$ over the open set $\mathcal{O}_ \mathsf{u}$. As for single species fluids, the original systems (6.4)–(6.6), that involves third order derivatives of density in the momentum equation, has been rewritten in the form of a quasilinear second order system. This hyperbolic-parabolic second order structure has been found convenient in order to derive a priori estimates, and to establish existence results as well as for numerical simulations ^{[82,83,84,85,86,87]}. The corresponding symmetrized equations have symmetric dissipative terms issued from ${\mathsf B}^{ {\rm{d}}}_{ij}$, $i, j\in \mathcal{D}$, and antisymmetric cohesive terms issued from ${\mathsf B}^{ {\rm{c}}}_{ij}$, $i, j\in \mathcal{D}$ ^[87].

6.5.   Extended thermodynamics

Using the local representation ${\boldsymbol{w}}$ of the density gradient, the thermodynamic functions are now local functions in the form

The matrix $\Lambda$ closely associated with the thermodynamic stability of multicomponent mixtures ^[9] may also be written as

The extended thermodynamic functions depend on $\mathsf{z} = \bigl(\rho_1, \ldots, \rho_ n, {\boldsymbol{w}}, {\boldsymbol{u}}, T \bigr)^t$ over an open set $\mathcal{O}_ \mathsf{z}\subset {\mathbb R}^ \mathsf{n}$. The open set $\mathcal{O}_ \mathsf{z}$ is said to be increasing with temperature if $(\rho_1, \ldots, \rho_ n, {\boldsymbol{w}}, {\boldsymbol{u}}, T)^t\in \mathcal{O}_ \mathsf{z}$ implies that $(\rho_1, \ldots, \rho_ n, {\boldsymbol{w}}, {\boldsymbol{u}}, T')^t\in \mathcal{O}_ \mathsf{z}$ for any $T'\ge T$. The thermodynamics functions $\mathcal{E}^ {\rm{ex}}$, $p^ {\rm{ex}}$, and $\mathcal{S}^ {\rm{ex}}$ defined over $\mathcal{O}_ \mathsf{z}$ are said to define an extended thermodynamics if the following properties hold.

$(\textsf{H}_{\textsf{0}} )$ $\mathcal{E}^ {\rm{ex}}$, $\mathcal{P}^ {\rm{ex}}$, $\mathcal{S}^ {\rm{ex}}$ are $C^ \gamma$ functions of $\mathsf{z}$ defined over $\mathcal{O}_ \mathsf{z} \subset (0, \infty)^ n\times {\mathbb R}^ d\times {\mathbb R}^ d\times(0, \infty)$. The open set $\mathcal{O}_ \mathsf{z}$ is nonempty and connected and $\varkappa = \varkappa(T)$ is a $C^{ \gamma+1}$ function of temperature $T$. If $(T, \rho_1, \ldots, \rho_ n)^t\in \mathcal{O}_ \mathcal{Z}$, then $(\rho_1, \ldots, \rho_ n, 0, 0, T)^t\in \mathcal{O}_ \mathsf{z}$, and if $(\rho_1, \ldots, \rho_ n, {\boldsymbol{w}}, {\boldsymbol{u}}, T)^t\in \mathcal{O}_ \mathsf{z}$, then $(T, \rho_1, \ldots, \rho_ n)^t\in \mathcal{O}_ \mathcal{Z}$.

$(\textsf{H}_{\textsf{1}} )$ Letting $g_k^ {\rm{ex}} = g = \partial_{\rho_k} \mathcal{E}^ {\rm{ex}} - T \partial_{\rho_k} \mathcal{S}^ {\rm{ex}}$, we have the extended volumetric Gibbs' relation $T \mathrm{d} \mathcal{S}^ {\rm{ex}} = \mathrm{d} \mathcal{E}^ {\rm{ex}} - \sum_{i\in \mathfrak{S}} g_i^ {\rm{ex}} \mathrm{d}\rho_i - \varkappa {\boldsymbol{w}} {\boldsymbol{\cdot}} \mathrm{d} {\boldsymbol{w}}$ and the constraint $\mathcal{E}^ {\rm{ex}} + p^ {\rm{ex}} - T \mathcal{S}^ {\rm{ex}} = \sum_{i\in \mathfrak{S}} \rho_i g_i^ {\rm{ex}}$.

$(\textsf{H}_{\textsf{2}} )$ For any $(T, \mathrm{y}_1, \ldots, \mathrm{y}_ n) \in(0, \infty)^{1+ n}$, such that $\sum_{i\in {\mathfrak S}} \mathrm{y}_i = 1$, and any ${\boldsymbol{u}}\in {\mathbb R}^ d$, there exists a threshold density $\rho_{{{\rm{m}}}} > 0$ such that $(\rho \mathrm{y}_1, \ldots, \rho \mathrm{y}_ n, 0, {\boldsymbol{u}}, T)^t \in \mathcal{O}_ \mathsf{z}$ for $0 < \rho < \rho_{{{\rm{m}}}}$ and

Assumption $(\textsf{H}_{\textsf{0}})$ describes the regularity class of extended thermodynamic functions, $(\textsf{H}_{\textsf{1}})$ is associated with the extended Gibbs relation for volumetric quantities, $(\textsf{H}_{\textsf{2}})$ and with the compatibility with ideal gases. These properties are natural extensions of those concerning the classical thermodynamics $( {\mathcal {T}}_{\textsf{0}}- {\mathcal {T}}_{\textsf{2}})$. We now introduce that Property $(\textsf{H}_{\textsf{3}})$ associated with the thermal stability condition.

$(\textsf{H}_{\textsf{3}})$ The open set $\mathcal{O}_ \mathsf{z}$ is increasing with temperature $T$, the specific heat is positive $\partial_T \mathcal{E}^ {\rm{ex}} >0$, and the capillarity coefficient is positive $\varkappa(T)>$ over $\mathcal{O}_ \mathsf{z}$.

The thermodynamic stability of the augmented model is now analyzed in the following proposition.

Proposition 6.1. Assuming that $(\textsf{H}_{\textsf{0}}-\textsf{H}_{\textsf{3}})$ hold, the Hessian matrix $\partial^2_{ \mathsf{u}, \mathsf{u}} \mathcal{S}^ {\rm{ex}}$ is negative definite if, and only if, $\Lambda$ is positive definite, $\partial_T \mathcal{E}^ {\rm{ex}}>0$, and $\varkappa>0$.

Proof. This is a consequence of the expression

where $\varrho = (\rho_1, \ldots, \rho_ n)^t$, ${1\kern -2.7pt \mathrm{I}} = (1, \ldots, 1)^t$, and $\mathsf{X} = (\mathsf{X}_\varrho, \mathsf{X}_ {\boldsymbol{w}}, \mathsf{X}_ {\boldsymbol{u}}, \mathsf{X}_ \mathcal{E})^t$ denotes an arbitrary vector of ${\mathbb R}^ \mathsf{n}$ with $\mathsf{X}_\varrho = (\mathsf{X}_1, \ldots, \mathsf{X}_ n)^t$. This expression then directly yields that stability is equivalent to $\partial_T \mathcal{E} > 0$, $\varkappa > 0$, and $\Lambda$ is positive definite, keeping in mind that $T$ and $\rho$ are positive. □

An important consequence of Proposition 6.1 is that the classical and extended thermodynamics have the same unstable states solely determined by the spectrum of $\Lambda = \Lambda^ {\rm{ex}}$ keeping in mind that thermal stability holds with $( {\mathcal {T}}_{\textsf{3}})$ and $(\textsf{H}_{\textsf{3}})$. Extra assumptions may also be introduced concerning the normal variable in the following form:

$(\textsf{H}_{\textsf{4}})$ The map $(\rho_1,\ldots,\rho_ n,T) \mapsto \Bigl(\rho,\frac{ g_2 - g_1}{T},\ldots, \frac{ g_ n - g_1}{T},T\Bigr)^t$ is a $C^{ \gamma-1}$ diffeomorphism from $\mathcal{O}_ \mathcal{Z}$ onto an open set $\mathcal{O}_ \mathcal{W}$.

Assuming $(\textsf{H}_{\textsf{0}}-$ $\textsf{H}_{\textsf{4}})$, we then deduce that the matrix $\Lambda$ is positive definite if, and only if, $\det\Lambda>0$. Therefore, in this situation, keeping in mind that thermal stability also holds, $\partial^2_{ \mathsf{u} \mathsf{u}} \mathcal{S}$ is negative definite if, and only if, $\det\Lambda>0$ and thermodynamic stability has been reduced to a scalar condition as for single species fluids where it is reduced to the mechanical stability condition ^[87].

The species of the mixture are assumed to be constituted by atoms, and we denote by ${\mathfrak a}_{il}$ the number of $l$th atom in the $i$th species, ${\mathfrak A} = \{1,\ldots, n^{{{\rm{a}}}}\}$ the set of atom indices, and $n^{{{\rm{a}}}}\geqslant1$ the number of atoms—or elements—in the mixture. The atomic vectors ${\mathfrak a}_l$, $l\in {\mathfrak A}$, are defined by ${\mathfrak a}_l = ( {\mathfrak a}_{1l},\ldots, {\mathfrak a}_{ n l})^t$. The assumptions for chemistry terms $(\textsf{CH})$ and transport coefficients $(\textsf{TRP}_)$ may also be written in the following form.

$(\textsf{CH})$ The atom conservation relations $\langle \nu_j, \mathfrak{a}_l \rangle = 0$ hold for $j\in {\mathfrak R},l\in {\mathfrak A}$ and $\mathcal{m} = \sum_{l\in {\mathfrak A}} {\widetilde m}_l \, \mathfrak{a}_{l}$. The rate constants $\mathcal{K}_j$, $j\in {\mathfrak R}$ are $C^ \gamma$ positive functions of $T>0$.

$(\textsf{TRP})$ The coefficients $\mathfrak{v}$, $\eta$, and the matrix $L$ are $C^ \gamma$ functions over $\mathcal{O}_ \mathsf{z}$. We have $\eta>0$, $\mathfrak{v}\ge0$, $\mathfrak{v} + \eta (1-\tfrac{2}{ d}) >0$, $L$ is symmetric positive semi-definite, and $N(L) = {\mathbb R}(1,\ldots,1,0)^t$.

Under the assumptions $(\textsf{H}_{\textsf{0}}-\textsf{H}_{\textsf{4}})$ and $(\textsf{TRP})$, it is possible to introduce a normal form for the system of equations with symmetric dissipative terms, antisymmetric capillary terms, and to obtain local existence theorems ^[87].

7.   Strained flows with diffuse interfaces

We investigate the low Mach number approximation of diffuse interface fluids and next present the self-similar strained flow equations that will be solved in the next section.

7.1.   Small Mach number setting

We assume that the Mach number is small and consider the situation where the interface is planar. We also assume that the capillarity coefficient $\varkappa$ is independent of temperature so that it is a constant. Proceeding as for classical fluids ^[8,29], we may expand the bulk pressure $p^ {}$ in powers of the square of the Mach number $\epsilon$:

with $\overline p^ {}$ denoting the zeroth-order thermodynamic pressure and $\widetilde p^ {}$ the fluid perturbation. The classic thermodynamic properties may then be evaluated at the zeroth-order pressure $\overline p^ {}$. The zeroth order pressure tensor is then in the form $\overline {\boldsymbol{\mathcal{P}}} = \overline p^ {} \boldsymbol{I}$-${\tfrac{1}{2}} \varkappa \vert \boldsymbol{\nabla}\mspace{-3mu}\rho\vert^ 2 \boldsymbol{I}$+$\varkappa \boldsymbol{\nabla}\mspace{-3mu}\rho {\otimes} \boldsymbol{\nabla}\mspace{-3mu}\rho$-$\varkappa \rho \Delta\rho \boldsymbol{I}$, and substituting this expansion in the momentum equation, it is found at zeroth-order that

This low Mach number relation is also the slow flow limit obtained by neglecting both convective and diffusive terms in the momentum equation. All other unknowns other than pressure are also expanded in terms of the square of the Mach number and—in order to simplify notation—are denoted as their zeroth-order value. The perturbed pressure $\widetilde p$ plays a role in the momentum equation, and from $p^ {\rm{ex}} = p^ {} - {\tfrac{1}{2}} \varkappa \vert \boldsymbol{\nabla}\mspace{-3mu}\rho\vert^2$ it is obtained that $p^ {\rm{ex}} = \overline p^ {} - {\tfrac{1}{2}} \varkappa \vert \boldsymbol{\nabla}\mspace{-3mu}\rho\vert^2 + \epsilon^ 2 \widetilde p$. For standard fluids where $\varkappa = 0$, we then recover that the ambiant pressure $\overline p^ {}$ is spatially uniform ^[8] since $\boldsymbol{\nabla} {\boldsymbol{\cdot}}(\overline p^ {} \boldsymbol{I}) = \boldsymbol{\nabla} \overline p^ {} = 0$. When capillary effects are included, the small Mach number relation (7.2) is not practical without further simplifying assumptions.

We assume that the interface is planar so that all quantities only depend on a normal coordinate $y$ and on time with $\rho_i = \rho_i(t, y)$, $i\in \mathfrak{S}$, $\rho = \rho(t, y)$, $\overline p^ {} = \overline p^ {}(t, y)$, and $T = T(t, y)$. In this framework, the zeroth order pressure tensor $\overline {\boldsymbol{\mathcal{P}}}$ is in the form $\overline {\boldsymbol{\mathcal{P}}} = \bigl(\overline p^ {} - {\tfrac{1}{2}} \varkappa \rho^{\prime \, 2} - \varkappa\rho \rho'' \bigr) \boldsymbol{I} + \varkappa \rho^{\prime \, 2} \mathsf{e}_ y{\otimes} \mathsf{e}_ y$, where $\mathsf{e}_ y$ denotes the unit vector normal to the interface, $\rho' = \partial_ y\rho(t, y)$, and $\rho'' = \partial_{ y y}^2\rho = \Delta\rho(t, y)$. From the normal component of the zeroth order momentum, we deduce that $\overline p^ {} + {\tfrac{1}{2}} \varkappa \rho^{\prime \, 2} - \varkappa\rho\rho''$ is a constant denoted by $p^\infty = \overline p^ {} + {\tfrac{1}{2}} \varkappa \rho^{\prime \, 2} - \varkappa\rho\rho''$. This pressure $p^\infty$ is simply the constant pressure of the bulk phases far from the interface. The resulting zeroth-order pressure tensor is then given by

so that normal to the interface, the effective pressure is constant given by $p^\infty$ whereas there are attractive forces tangential to the interface. The corresponding energy $\int_{-\infty}^{\infty} \varkappa \rho^{\prime \, 2} \, d y$ may then be interpreted as surface tension ^[25,27,28].

For planar interfaces, only the perturbed pressure $\widetilde p^ {}$ is playing a role in the tangential momentum conservation equation and not the perturbed density $\widetilde\rho$. The zeroth-order energy conservation equation in terms of the bulk enthalpy $h^ {\rm{ex}} = h^ {}$ also takes the simplified form

with a term ${\boldsymbol{u}} {\boldsymbol{\cdot}} \boldsymbol{\nabla} \overline p^ {}$ that is not anymore a priori negligible as for standard fluids. Note incidentally that there is a typographic error in reference ^[29] with $\boldsymbol{\mathcal{Q}}$ instead of $\boldsymbol{q}$ in the bulk enthalpy equation as noted in ^[30].

7.2.   Strained flows

The geometry is considered to be either two-dimensional or axisymmetric and the spatial coordinates are denoted by $(x, y)$ where $x$ is parallel to the interface and $y$ is normal to it. The corresponding components of the velocity vector ${\boldsymbol{u}}$ are denoted by ${\boldsymbol{u}} = ({\rm{U}}, {\rm{V}})^t$. The solution is assumed to have the following self-similar structure ^{[8,29,30,76,92,93]}

where $J$ is a reduced tangential pressure curvature. It is deduced then from (6.11) and (6.12) that the heat and mass fluxes only have a normal component $\boldsymbol{\mathcal{J}}_i = \mathcal{J}_i \mathsf{e}_ y$, $i\in \mathfrak{S}$, and $\boldsymbol{q} = q \mathsf{e}_ y$ and $\mathcal{J}_i(t, y)$, $i\in \mathfrak{S}$, and $q(t, y)$ are the normal flux components. Using such expressions in the equations in the low number approximation, the normal momentum equation is found to uncouple and the tangential momentum equation has its terms proportional to $x$. The resulting system forms a two-point boundary value problem given by

where $\delta = 0$ in a two-dimensional geometry and $\delta = 1$ for axisymmetric flows. Incidentally there is also a typographic error in the strained flow enthalpy equation of reference ^[29], but the correct enthalpy Eq (7.7) has properly been used in the simulations ^[29,30]. The dissipative normal mass and heat fluxes are given by (6.11) and (6.12) in their thermodynamic form

and it has been established in Section 7.1 that

where $p^ \infty$ denotes the constant ambient pressure far from the interface. The boundary conditions are in the form

where the superscript $+$ is associated with the fluid coming from the positive side on the right, the superscript $-$ to the fluid coming from the negative side on the left, and $\alpha$ is the strain rate related to the pressure curvature with the relation

The stagnation plane is finally located at the origin with

In the situation where the flow is steady, the governing equations are then simplified in the form

and form a steady two-point boundary value problem. Note that with the small Mach number relation (7.10) we do not need the augmented formulation in order to obtain a set of equations only involving only second order derivatives.

7.3.   Dissipative transport fluxes

In order to evaluate the multicomponent fluxes, the reduced potential $\mu_i = m_i g_i/RT$ may be decomposed in the form

so as to isolate the singular component $\log \mathsf{X}_i$ when the mole fraction $\mathsf{X}_i$ goes to zero whereas the component $\mu_i^ \mathrm{sm}$ remains smooth. We may thus evaluate $\boldsymbol{\nabla}{\mu_{{i}}}$ in the form $\boldsymbol{\nabla} \mathsf{X}_i/ \mathsf{X}_i + \boldsymbol{\nabla}\mu_i^ \mathrm{sm}$ and the smooth part $\mu_i^ \mathrm{sm}$ is given by

where $g_i^{ {\rm{id}}\star}$ denotes the ideal gas Gibbs function of the $i$th species at the standard pressure $p^ \mathrm{st}$.

The mass and heat transport coefficient matrix $L$ of size $n+1$ may be written ^[18,80]:

where $\mathcal{D}$ is a matrix of size $n$, $\mathcal{h}$ a vector of size $n$, $\lambda$ the thermal conductivity, and $\langle \; , \; \rangle$ the Euclidean scalar product. When thermodynamic stability holds, the matrix coefficients $\mathcal{D}_{ij}$, $i, j\in \mathfrak{S}$, and the vector $\mathcal{h}_i$, $i\in \mathfrak{S}$, are given by

where $D_{ij}$, $i, j\in \mathfrak{S}$, denote the multicomponent diffusion coefficients, $h_i$, $i\in \mathfrak{S}$, the species specific enthalpies defined by $h_i = \widehat\partial^{}_{ \mathrm{y}_i} h$, $i\in \mathfrak{S}$, where the derivative is thus taken at constant pressure, and $\widetilde\chi_i$, $i\in \mathfrak{S}$, the reduced thermal diffusion ratios ^[18,80]. The multicomponent diffusion coefficients $D_{ij}$ are obtained with an approximate inversion ^{[8,105,106,107,108,109,110,111]} using

where $\mathrm{y}$ denotes the mass fraction vector, ${1\kern -2.7pt \mathrm{I}} = (1, \ldots, 1)^t$ the vector with unit components, $\Delta$ the Stefan-Maxwell matrix

$\mathsf{X}_k$ the mole fraction of the $k$th species, and ${\mathcal D}^ {{{\rm{bin}}}}_{kl}$, the binary coefficient for the species pair $(k, l)$, $k\neq l$. The binary diffusion coefficients ${\mathcal D}^ {{{\rm{bin}}}}_{kl}$, $k, l\in \mathfrak{S}$, $k\neq l$, are calculated by using the kinetic theory of dense hard sphere mixtures ^[80,116] and read ${\mathcal D}_{kl}^ {{{\rm{bin}}}} = {\mathcal D}_{kl}^{ {{{\rm{bin}}}}\, {\rm{id}}}/\Upsilon_{kl}$, $k, l\in \mathfrak{S}$, $k\neq l$, where ${\mathcal D}_{kl}^{ {{{\rm{bin}}}}\, {\rm{id}}}$ denotes the ideal gas binary diffusion coefficient and $\Upsilon_{ij}$ the steric factor

We have denoted here by $\sigma_{ij}$ the collision diameter of the species pair $(i, j)$, ${\mathfrak n}_k = \rho_k {\mathcal N}/m_k$ the particle number density of the $k$th species, and by ${\mathcal N}$ the Avogadro number. The viscosity and thermal conductivity may be calculated from the correlations of Ely and Hanley ^[114] or Chung et al. ^[80,115], and the thermal diffusion ratios $\widetilde\chi_i$, $i\in \mathfrak{S}$ are evaluated as for ideal gases ^[80].

The properties of the matrix $L$ are derived from the properties of the matrix $D = (D_{ij})_{i, j\in \mathfrak{S}}$ which is symmetric positive semi-definite with its null space $N(D) = {\mathbb R} \mathrm{y}$ ^{[8,107,108,109,110,111]}. When thermodynamic stability does not hold, the expression (7.22) is still used with the same binary multicomponent diffusion coefficients, but some care must be taken in order to evaluate the species' specific enthalpies $h_i = \widehat\partial^{}_{ \mathrm{y}_i} h$, $i\in \mathfrak{S}$. We indeed have the identity

so that at a mechanical thermodynamic unstable point where $\partial^{}_\rho \overline p^ {} = 0$, the species' specific enthalpies $h_k$, $k\in \mathfrak{S}$, are exploding and, thus, also some components of the matrix $L$. This is an expected behavior since it is known that singularities in the transport coefficients may arise at critical points. It is therefore necessary to control the size of the vector $\mathcal{h}$. To this aim, the derivative $\partial^{}_{ \mathrm{y}_k} \overline p^ {}$ is rewritten in the form

where the residual $\mathcal{R}_k$ is given by

The interest of this formulation (7.25) is that $\mathcal{R}_k$ vanishes for the pure species state $\mathrm{y}_k = 1$ and $\mathrm{y}_l = 0$, $l\in \mathfrak{S}$, $l\neq k$. The species enthalpies are then written in the modified form $\widetilde h_k$, $k\in \mathfrak{S}$, with

where $\mathfrak{f}$ is a smooth approximation of the function $\mathfrak{f}(x) = \max({\rm m}, x)$ for a positive constant ${\rm m} > 0$ which depends on the mixture upon consideration. Such revised enthalpies are such that $\widetilde h_k = h_k$ away from the vaporizing zone, and when the $k$th species is the only species present, it remain always bounded. The points where the matrix $L$ is effectively changed only concern zones of the fluid where temperature is essentially constant so that the stabilization of the matrix $L$ does not modify the physics involved. For similar reasons, we did not include in the model the complex behavior of transport coefficients near critical points. Note finally that the traditional formulation of transport fluxes is not recommended and the thermodynamical formulations (7.8) and (7.9) are to be preferred since they eliminate the many exploding quantities at mechanical instability points.

7.4.   Nonlinear solvers

Thermodynamic functions and chemical production rates may conveniently be evaluated by using application independent libraries like CHEMKIN ^[94] or CANTERA ^[95] since this allows us, in particular, to develop application independent software suites.

Nonlinear equations typically associated with multiphase equilibrium or critical points may be solved by using Newton's method or any of its variants ^[96]. Once a first solution is obtained, continuation techniques may then be used to generate solution branches depending on a parameter ^[97,101]. Continuation techniques reparameterize the arc of solutions using pseudo-arclength coordinates ^[97,101]. These methods have notably been used in order to investigate the location of nontrivial zero eigenvalues of entropy Hessian matrices or for the determination of thermodynamic stability domains for binary mixtures.

7.5.   Two-point boundary value problems

The two-point boundary value problem (7.4)–(7.15) associated with strained diffuse interfaces has been discretized by using finite differences on staggered grids. The normal velocity is discretized at midpoints whereas the other quantities are discretized at node points. The set of discretized equations may be solved by using Newton's method and self-adaptive grids ^[98,99]. The grid points are located by equidistributing a weight function associated with gradients or curvatures ^[97,98]. The Jacobian matrices have a block-tridiagonal structure and are solved by direct methods. Pseudo unsteady iterations may also be used to bring the initial guess into the domain of convergence of steady Newton's iterations ^[98,99]. Continuation techniques may also be used to generate families of diffuse interface strained flows ^[101].

Evaluating aero-thermochemistry quantities is computationally expensive since they involve multiple sums and products. Highly optimized thermochemistry and transport routines have been used and further extended to the high pressure domain ^{[102,103,104,105,106,107,108,109,110,111]}.

8.   Validation of the formalism

In order to validate the thermodynamic formalism developed in the previous sections, numerical simulations are compared with experimental measurements. We investigate, more specifically, thermodynamic properties of ethane and nitrogen and thermodynamic stability of C$_2$H$_6$-N$_2$ mixtures. We also simulate numerically strained flows of ethane with diffuse interfaces in order to illustrate the use of the extended formalism.

8.1.   Parameters for the SRK equation and capillarity

The pure species attractive and repulsive parameters per unit mass $\alpha_i = \sqrt{ a_i}$ and $b_i$, $i\in {\mathfrak S}$, associated with the SRK equation of state (5.1) and (5.2), can be evaluated from the critical data ^[80]. Assuming that the $i$th species is chemically stable, the attractive parameter $\alpha_{ {{{\rm{c}}}}, i}$ at the critical temperature $T_{ {{{\rm{c}}}}, i}$ of Proposition 5.4 and the repulsive parameter $b_i$ are evaluated in the form

where $T_{ {{{\rm{c}}}}, i}$ and $p_{ {{{\rm{c}}}}, i}$ are the critical temperature and pressure of the $i$th species. The attractive and repulsive parameters of chemically stable species like C$_2$H$_6$, O$_2$, N$_2$, or H$_2$O, or metastable species like H$_2$O$_2$, may thus be determined from such critical state conditions. The quantity ${\rm{s}}_{i}$ from Proposition 5.4 may also be expressed in the form

where $\varpi_{{i}}$ is the $i$th species acentric factor ^[113]. The attractive and repulsive parameters and the acentric factors that have been used for ethane C$_2$H$_6$ and nitrogen N$_2$ are presented in Table 1.

This procedure, however, cannot be generalized to chemically unstable species like H, O, or OH radicals since critical states do not exist for such species ^[80]. Assuming that the $i$th species is a Lennard-Jones gas, however, it is possible to estimate ^[74,75] a pseudo-critical state and to deduce that

where ${\rm{N}}$ is the Avogadro number, and $\sigma_{{i}}$ and $\epsilon_{{i}}$ are the molecular diameter and Lennard-Jones potential well depth of the $i$th species. Note that this procedure is less accurate for stable species ^[80]. Finally, the thermodynamic properties of ideal gases may be evaluated from the SANDIA database, the NIST/JANAF Thermochemical Tables, or from the NASA coefficients ^{[113,117,118,119]}. It is important to note that the required parameters for the SRK equation of state or the ideal gas mixtures are obtained from independent physical measurements or quantum simulations of fundamental fluid properties in such a way that there are no adjustable parameters in the thermodynamic model.

The capillarity coefficient may be obtained from experimental measurements of surface tension but also increased in order to artificially thicken the diffuse interfaces as done in this paper. Thickening an interface may generally modify surface tension forces so that there is a trade-off between grid limitation and the resulting artificially increased surface tension.

8.2.   Thermodynamic properties of nitrogen

We investigate the specific heat at constant pressure $c_p(T, p)$ of pure nitrogen as a function of temperature $T$ at $p = 10$ MPa. For a given pressure $p_ {{{\rm{ref}}}}$ and temperature $T_ {{{\rm{ref}}}}$, for a given pure species with $n = 1$, $\mathfrak{S} = \{1\}$, $\mathrm{y}_1 = 1$, the corresponding specific volume(s) $v$ may be obtained by solving the cubic equation $p(T_ {{{\rm{ref}}}}, v, \mathrm{y}_1) = p_ {{{\rm{ref}}}}$. Denoting by $v_ {{{\rm{inf}}}}$ and $v_ {{{\rm{sup}}}}$ the lowest and highest solutions, the fluid is liquid with volume $v_ {{{\rm{inf}}}}$ if $g(T_ {{{\rm{ref}}}}, v_ {{{\rm{inf}}}}, \mathrm{y}_1) < g(T_ {{{\rm{ref}}}}, v_ {{{\rm{sup}}}}, \mathrm{y}_1)$, gaseous with volume $v_ {{{\rm{sup}}}}$ if $g(T_ {{{\rm{ref}}}}, v_ {{{\rm{sup}}}}, \mathrm{y}_1) < g(T_ {{{\rm{ref}}}}, v_ {{{\rm{inf}}}}, \mathrm{y}_1)$, and both states coexist when $g(T_ {{{\rm{ref}}}}, v_ {{{\rm{inf}}}}, \mathrm{y}_1) = g(T_ {{{\rm{ref}}}}, v_ {{{\rm{sup}}}}, \mathrm{y}_1)$. The constant pressure heat capacity is then taken to be the partial derivative with respect to $T_ {{{\rm{ref}}}}$ of $h\bigl(T_ {{{\rm{ref}}}}, v_ {{{\rm{inf}}}}(T_ {{{\rm{ref}}}}, p_ {{{\rm{ref}}}}), \mathrm{y}_1\bigr)$ or $h\bigl(T_ {{{\rm{ref}}}}, v_ {{{\rm{sup}}}}(T_ {{{\rm{ref}}}}, p_ {{{\rm{ref}}}}), \mathrm{y}_1\bigr)$ depending on the gaseous or liquid state of the fluid with $h = e + p/\rho$.

The constant pressure heat capacity $c_p(T, p)$ of nitrogen N$_2$ is presented in Figure 1. The constant pressure heat capacity of ideal gases is also presented and only agrees with the nonideal heat capacity and sufficiently large temperature. There is also a very good agreement with the experimental data from Youglove ^[121] and Youglove and Ely ^[122].

The saturation pressure $p_ {{{\rm{sat}}}}$ of nitrogen as a function of temperature $T$ is obtained by solving the nonlinear equation $g(T_ {{{\rm{ref}}}}, v_ {{{\rm{inf}}}}, \mathrm{y}_1) = g(T_ {{{\rm{ref}}}}, v_ {{{\rm{sup}}}}, \mathrm{y}_1)$, where $v_ {{{\rm{inf}}}}$ and $v_ {{{\rm{sup}}}}$ are the lowest and highest solutions of the cubic equation $p(T_ {{{\rm{ref}}}}, v, \mathrm{y}_1) = p_ {{{\rm{ref}}}}$, and we then have $p_ {{{\rm{ref}}}} = p_ {{{\rm{sat}}}}(T_ {{{\rm{ref}}}})$. The saturation pressure of nitrogen as a function of temperature is presented in Figure 2, and very good agreement is obtained with the experimental data from Youglove ^[121]. The mechanical stability limits are also presented on this plot and the saturation curve logically stays within the unstable zone. A similar curve has been obtained for ethane with very good agreement with the experimental data from Youglove and Ely ^[122], and these simulations validate the nonideal thermodynamics.

8.3.   Thermodynamic of ethane and nitrogen mixtures

We now consider a mixture of ethane and nitrogen and the species indexing set may be taken in both forms $\mathfrak{S} = \{1, 2\} = \{ {\rm C}_2{\rm H}_6, {\rm N}_2 \}$. Thermal stability holds for the SRK thermodynamics and it has been checked that the quantity $\Lambda_{2, 2} + \Lambda_{1, 1} - \Lambda_{1, 2} - \Lambda_{2, 1}$ remains positive for mixtures of ethane and nitrogen in such a way that thermodynamic stability is equivalent to $\det\Lambda > 0$ as in Corollary 3.8. The thermodynamic stability domain of the mixture may thus be determined by solving the equation $\det\Lambda(T, v, \mathrm{y}) = 0$ and the mechanical stability domain—of the mixture—by solving the equation $\partial_ v p(T, v, \mathrm{y}) = 0$. We have used a continuation in the variable $v$ for solving the system

that is, used nonlinear solvers and continuation methods ^[96,98,101] in order to generate the whole stability domain. We also used in practice the rescaled matrix $\widehat\Lambda^ \mathsf{r}$ of Proposition 4.2 that is more convenient than $\Lambda$. The thermodynamic stability domain and the mechanical stability domain are presented at temperature $T = 200$ K in the plane $(\mathrm{y}_{{\rm C}_2{\rm H}_6}, v)$ on Figure 3. The unstable zones are inside the curves and the mechanical unstable domain is naturally interior to the thermodynamic stability domain.

Multiphase equilibrium mixtures of ethane and nitrogen have also been investigated. For a given $T_ {{{\rm{ref}}}}$, we have used a continuation in terms of the $v^\sharp$ variable for solving the system of equations:

The equilibrium curve between ethane and nitrogen at $T = 220$ K is presented in Figure 4 in the $(\mathrm{y}_{{\rm C}_2{\rm H}_6}, p)$ plane. The experimental measurements are taken from Youglove ^[121] and Youglove and Ely ^[122] with an overall very good agreement with the simulations. The top of the numerical equilibrium curve is a critical point of the binary mixture, i.e., a maximum of pressure with respect to the relative concentration.

We have also investigated the critical points of ethane-nitrogen mixtures. The critical points of pure mixtures are solutions of $\partial_ v p = 0$ and $\partial_ v^2 p = 0$, but the situation of mixtures is more complex. Critical points may be obtained by using a continuation in the $v_ {{{\rm{c}}}}$ of the system of equations ^[128]

Only the relevant branch starting from the critical point of pure ethane given by $T_ {{{\rm{c}}}} = 305.3$ K and the $p_ {{{\rm{c}}}} = 4.9$ bar is investigated ^[112]. The critical diagram in the plane $(T, p)$ is presented in Figure 5 and found to be of type Ⅲ according to the classification of Konynenburg and Scott ^[112,113]. Binary O$_2$-H$_2$O mixtures also have such a type Ⅲ critical point diagram ^[29,127]. A very good agreement is found with the experimental measurements of Eakin at al. ^[124], as well as those of Stryjek et al. ^[125]. In the lower temperature part of the curve, the tendencies are well reproduced with a temperature shift of 30 K in comparison with the measurements of Wisotzki and Schneider ^[126], and all these simulations validate the nonideal mixture thermodynamics. The critical diagram in the plane $(\mathrm{y}_{{\rm C}_2{\rm H}_6}, p)$ is also presented in Figure 6, and we notably observe a turning point with respect to the ethane mass fraction $\mathrm{y}_{{\rm C}_2{\rm H}_6}$ that is easily taken into account by using continuation methods.

8.4.   Strained flows of ethane

We investigate a steady strained flow of ethane governed by the system of Eqs (7.16)–(7.19) with the boundary conditions (7.11)–(7.15). The flow parameters are taken to be

The capillarity coefficient is varied to be

For these pseudo-vaporizing interfaces, C$_2$H$_6$ is the only chemical species and no species equations are solved. The mass density $\rho(y)$ normal to the interface is presented as function of $y$ in Figure 7 for the various capillarity coefficients. The liquid-like phase is on the righthand side and the gaseous phase is on the left side.

The capillary parameter $\varkappa$ has an important impact on the interface structure. The lower the value of $\varkappa$, the sharper the interface, and when $\varkappa\to0$, the infinitely thin interface—discontinuous—solution is recovered. The corresponding local pressure $\overline p$ curves, at subcritical ambient pressure $p^\infty$, then essentially follows an isotherm curve of the SRK pressure law at the interface as analyzed in ^[29,30]. This would not be the case at supercritical pressure where the local pressure variation in the interfacial region goes to zero with the capillary parameter $\varkappa$ ^[29,30].

9.   Conclusions

We have investigated the mathematical structure of multicomponent fluid thermodynamic allowing instabilities. We have also investigated the construction of such thermodynamics from equations of states with a special emphasis on the SRK state law. Numerical simulation of C$_2$H$_6$-N$_2$ chemical instabilities as well as strained fronts have shown the applicability of the resulting thermodynamics. Overall, a very good agreement has been found with experimental data of thermodynamic properties without any adjustable parameters.

The thermodynamic formalism presented in the paper could potentially be used in association with any compressible mixture fluid model. It may be used, for instance, to investigate gradient flow structures—valid in the absence of convection phenomena and numerical stability issues, as well as Cahn-Hilliard models. The direct numerical simulation of multidimensional flows with diffuse interfaces using an augmented formulation of the system of partial differential equations would also be of high scientific interest.

Author contributions

V. Giovangigli: Conceptualization, investigation, software, writing; Y. Le Calvez: Investigation, software, numerical simulations, data curation; G. Ribert: Investigation, resources. All authors have read and approved the final version of the manuscript for publication.

Acknowledgments

This work was supported by the ANR INSIDE project, grant ANR-19-CE05-0037-02 of the French Agence Nationale de la Recherche.

Conflict of interest

On the behalf of the authors there are no conflict of interest with the authors.

Appendix

A.   The single species case

In this section, we discuss for completeness the situation where there is only one species with $n = 1$ and $\mathfrak{S} = \{1\}$. It is first possible to apply the general formalism to this special situation, but further simplifications are usually introduced to abridge the model by eliminating formally $\mathrm{y}_1$. Since the mass fractions are assumed to be independent, the mass fraction of the single species $\mathrm{y}_1$ is indeed not a priori unity and must be taken into account in the general Gibbsian homogeneous formalism. The natural variable then reads $\zeta = (T, v, \mathrm{y}_1)^t$, the thermodynamic is defined with $e(\zeta)$, $p(\zeta)$, and $s(\zeta)$, and Gibbs relation is in the form

We may, however, define for convenience the specific volume $v_1$ by setting

From the 1-homogeneity of the thermodynamic functions $e(\zeta)$ and $s(\zeta)$ and the 0-homogeneity of $p(\zeta)$, we then obtain that $e(T, v, \mathrm{y}_1) = \mathrm{y}_1 e(T, v_1, 1)$, $p(T, v, \mathrm{y}_1) = p(T, v_1, 1)$, and $s(T, v, \mathrm{y}_1) = \mathrm{y}_1 s(T, v_1, 1)$. We may thus define the densities $e_1(T, v_1) = e(T, v_1, 1)$, $p_1(T, v_1) = p(T, v_1, 1)$, and $s_1(T, v_1) = s(T, v_1, 1)$, in such a way that

From the definition $g_1 = \partial^{}_{ \mathrm{y}_1}\! e - T \partial^{}_{ \mathrm{y}_1}\! s$, we also deduce that

so that

We have $g_1 = e_1 + v_1 p_1 - T s_1$, and since $g_1$ is 0-homogeneous, we also have $g_1(\zeta) = g(T, v_1, 1)$.

From Gibbs relation $T \mathrm{d} s = \partial^{}_T e \; \mathrm{d} T + (p + \partial^{}_ v e) \, \mathrm{d} v + (\partial^{}_{ \mathrm{y}_1}\! e - g_1) \, \mathrm{d} \mathrm{y}_1$, we now deduce, after some algebra, that

and this yields the usual relation

that is, the classical Gibbs relation without the homogeneous framework. In other words, the thermodynamics functions generally used when $n = 1$ are $e_1$, $p$, $s_1$ functions of $(T, v_1)$ that differ formally from $e$, $p$, $s$ that are functions of functions of $(T, v, \mathrm{y}_1)$, but the link is easily made.

B.   Proof of Propostion 2.6

We present in this appendix the proof of Proposition 2.6 and denote by $\mathsf{f}^ e$, $\mathsf{f}^ v$, $\mathsf{f}^1, \ldots, \mathsf{f}^ n$ the natural basis vector of ${\mathbb R}^{2+ n}$.

Proof. We first calculate the coefficients of the Hessian matrix $\widetilde\partial^2_{ \xi \xi} s$ when $\zeta\mapsto \xi$ is a $C^ \gamma$ local diffeomorphism. From $\widetilde\partial^{}_ e s = 1/T$, we get that $T^2 \widetilde\partial^2_{ e e}s = - \widetilde\partial^{}_ e T$ and $T^2 \widetilde\partial^2_{ e v}s = - \widetilde\partial^{}_ v T$. Moreover, for any function $\chi$ of $\zeta$ or $\xi$, we have the differential relations

and letting $\chi = T$ in (B.2), we deduce that $\widetilde\partial^{}_ v T = - (\widetilde\partial^{}_ e T) \partial^{}_ v e$ and, finally, $T^2 \widetilde\partial^2_{ e v}s = (\widetilde\partial^{}_ e T) \partial^{}_ v e$. Similarly, we have $T^2 \widetilde\partial^2_{ e \mathrm{y}_k}s = - \widetilde\partial^{}_{ \mathrm{y}_k} T$, and using $\chi = T$ in (B.3), we deduce $\widetilde\partial^{}_{ \mathrm{y}_k} T = -(\widetilde\partial^{}_ e T) \partial^{}_{ \mathrm{y}_k} e$ so that $T^2 \widetilde\partial^2_{ e \mathrm{y}_k}s = (\widetilde\partial^{}_ e T) \partial^{}_{ \mathrm{y}_k} e$. Furthermore, we have $\widetilde\partial^2_{ v v}s = \widetilde\partial^{}_ v \bigl(\frac{p}{T}\bigr)$, and using (B.2), we deduce that $\widetilde\partial^{}_ v\bigl(\frac{p}{T}\bigr) = \partial^{}_ v\bigl(\frac{p}{T}\bigr) + \partial^{}_T\bigl(\frac{p}{T}\bigr) \widetilde\partial^{}_ v T$. From the compatibility conditions (2.8) and $\widetilde\partial^{}_ v T = - (\widetilde\partial^{}_ e T) \partial^{}_ v e$, we thus obtain $\widetilde\partial^2_{ v v}s = \partial^{}_ v\bigl(\frac{p}{T}\bigr) - (\widetilde\partial^{}_ e T) \frac{(\partial^{}_ v e)^2}{T^2}$. Similarly, $\widetilde\partial^2_{ v \mathrm{y}_k}s = \widetilde\partial^{}_{ \mathrm{y}_k} \bigl(\frac{p}{T}\bigr)$, and thanks to (B.3), we have $\widetilde\partial^{}_{ \mathrm{y}_k}\bigl(\frac{p}{T}\bigr) = \partial^{}_{ \mathrm{y}_k}\bigl(\frac{p}{T}\bigr) + \partial^{}_T\bigl(\frac{p}{T}\bigr) \widetilde\partial^{}_{ \mathrm{y}_k}T$. From the compatibility conditions (2.8) and $\widetilde\partial^{}_{ \mathrm{y}_k} T = - (\widetilde\partial^{}_ e T) \partial^{}_{ \mathrm{y}_k} e$, we get

Since $\widetilde\partial^2_{ \mathrm{y}_k \mathrm{y}_l}s = \widetilde\partial^{}_{ \mathrm{y}_l} \bigl(\frac{- g_k}{T}\bigr)$, and using $\chi = \frac{- g_k}{T}$ in (B.3), we also obtain that

Defining $\mathsf{\widehat f}^e = \bigl(1, - \partial^{}_ v e, - \partial^{}_{ \mathrm{y}_1} e, \ldots, - \partial^{}_{ \mathrm{y}_ n} e \bigr)^t$, we have established that for any $\mathsf{X} = (\mathsf{X}_ e, \mathsf{X}_ v, \mathsf{X}_1, \ldots, \mathsf{X}_ n)^t\in {\mathbb R}^{2+ n}$

We now prove that $(i)$ implies $(ii)$. Since $\xi = e \mathsf{f}^ e + v \mathsf{f}^ v + \sum_{i\in {\mathfrak S}} \mathrm{y}_i \mathsf{f}^{ \mathrm{y}_i}$ and $\mathsf{f}^ e$ are not proportional we deduce that $\bigl\langle (\widetilde\partial^2_{ \xi \xi}s) \mathsf{f}^ e, \mathsf{f}^ e \bigr\rangle < 0$ so that $\widetilde\partial^{}_ e T = - T^2 \bigl\langle (\widetilde\partial^2_{ \xi \xi}s) \mathsf{f}^ e, \mathsf{f}^ e \bigr\rangle > 0$, and this yields $\partial^{}_T e > 0$ since $\widetilde\partial^{}_ e T\; \partial^{}_T e = 1$. Letting $\mathsf{f} = \mathsf{f}^ v + (\partial^{}_ v e) \mathsf{f}^e$, we have $\langle \, \mathsf{\widehat f}^e, \mathsf{f} \rangle = 0$ and $\mathsf{f} = \mathsf{f}^ v + (\partial^{}_ v e) \mathsf{f}^e$, and since the mass fractions are positive, $\xi$ cannot be proportional. Thus, $\partial^{}_ v p/T = \bigl\langle (\widetilde\partial^2_{ \xi \xi}s) \mathsf{f}, \mathsf{f} \bigr\rangle < 0$, and we have proved that $\partial^{}_ v p < 0$. In the same vein, assume that $\mathsf{X}_1, \ldots, \mathsf{X}_ n$ are arbitrary with $(\mathsf{X}_1, \ldots, \mathsf{X}_ n)\neq (0, \ldots, 0)$, consider $\mathsf{f} = \sum_{i\in {\mathfrak S}} \mathsf{X}_i \mathsf{f}^i + \mathsf{X}_e \mathsf{f}^e$, and select $\mathsf{X}_e$ such that $\langle \, \mathsf{\widehat f}^e, \mathsf{f} \rangle = 0$, that is, $\mathsf{X}_e = \sum_{i\in {\mathfrak S}} \mathsf{X}_i \partial^{}_{ \mathrm{y}_i} e$. Then, $\mathsf{f}$ and $\xi$ cannot be proportional since the volume per unit mass is positive, so that $\bigl\langle (\widetilde\partial^2_{ \xi \xi}s) \mathsf{f}, \mathsf{f}\bigr\rangle < 0$. This yields that $\sum_{k, l\in {\mathfrak S}} \partial^{}_{ \mathrm{y}_l} \bigl(\frac{ g_k}{T} \bigr)\, \mathsf{X}_k \mathsf{X}_l = - \bigl\langle (\widetilde\partial^2_{ \xi \xi}s) \mathsf{f}, \mathsf{f}\bigr\rangle > 0$ and the matrix $\Lambda$ is positive definite.

We now prove that $(ii)$ implies $(i)$. From the identity (2.7), we deduce that $\partial^{}_ v p < 0$. Defining $\mathsf{\widehat f}^ v = \bigl(0, 1, - \frac{ \partial^{}_ v g_1}{ \partial^{}_ v p}, \ldots, - \frac{ \partial^{}_ v g_ n}{ \partial^{}_ v p} \bigr)^t$ and using (B.4), we easily deduce that for $\mathsf{X} = (\mathsf{X}_ e, \mathsf{X}_ v, \mathsf{X}_1, \ldots, \mathsf{X}_ n)^t$

Thanks to (2.6) and (2.7), this identity is rewritten in the form

where $\mathsf{X}_ \mathrm{y} = (\mathsf{X}_1, \ldots, \mathsf{X}_ n)^t$. Since $\Lambda$ is positive definite, we first note that $\Lambda - \frac{\Lambda \mathrm{y}{\otimes}\Lambda \mathrm{y}}{ \langle\Lambda \mathrm{y}, \mathrm{y}\rangle}$ is positive semi-definite with null space ${\mathbb R} \mathrm{y}$ since for any $\mathsf{X}_ \mathrm{y} = (\mathsf{X}_1, \ldots, \mathsf{X}_ n)^t$, we have

Since $\partial^{}_T e > 0$, $\partial^{}_ v p < 0$, and $\Lambda - \frac{\Lambda \mathrm{y}{\otimes}\Lambda \mathrm{y}}{ \langle \Lambda \mathrm{y}, \mathrm{y}\rangle}$ is positive semi-definite with null space ${\mathbb R} \mathrm{y}$, it is obtained from (B.6) that $\widetilde\partial_{ \xi \xi}^2 s$ is negative semi-definite. Moreover, any vector $\mathsf{X}$ in the null space of $\widetilde\partial_{ \xi \xi}^2 s$ is such that $\langle \, \mathsf{\widehat f}^ e, \mathsf{X} \rangle = 0$, $\langle \, \mathsf{\widehat f}^ v, \mathsf{X} \rangle = 0$, $(\Lambda - \frac{\Lambda \mathrm{y}{\otimes}\Lambda \mathrm{y}}{ \langle \Lambda \mathrm{y}, \mathrm{y}\rangle} \bigr) \mathsf{X}_ \mathrm{y} = 0$. Since $\mathsf{X}_ \mathrm{y} = (\mathsf{X}_1, \ldots, \mathsf{X}_ n)^t$ is in the null space of $\Lambda - \frac{\Lambda \mathrm{y}{\otimes}\Lambda \mathrm{y}}{ \langle \Lambda \mathrm{y}, \mathrm{y}\rangle}$ we first deduce that $\mathsf{X}_ \mathrm{y} = \lambda \mathrm{y}$ for some $\lambda\in {\mathbb R}$. From the expression of $\mathsf{\widehat f}^ v$, we first obtain that $\mathsf{X}_ v = \sum_{i\in {\mathfrak S}} \mathsf{X}_i \partial^{}_ v g_i/ \partial^{}_ v p$ so that $\mathsf{X}_ v = \lambda \sum_{i\in {\mathfrak S}} \mathrm{y}_i \partial^{}_ v g_i/ \partial^{}_ v p$. Using then the compatibility relations (2.8), we have $- \partial^{}_ v g_i = \partial^{}_{ \mathrm{y}_i} p$ so that $\mathsf{X}_ v = - \lambda \sum_{i\in {\mathfrak S}} \mathrm{y}_i \partial^{}_{ \mathrm{y}_i} p/ \partial^{}_ v p$ and, finally, $\mathsf{X}_ v = \lambda v$ because of the 0-homogeneity of $p$. In addition, from the expression of $\mathsf{\widehat f}^ e$, we get that $\mathsf{X}_ e = \mathsf{X}_ v \partial^{}_ v e + \sum_{i\in {\mathfrak S}} \mathsf{X}_i \partial^{}_{ \mathrm{y}_i} e_i$ so that $\mathsf{X}_ e = \lambda v \partial^{}_ v e + \lambda \sum_{i\in {\mathfrak S}} \mathrm{y}_i \partial^{}_{ \mathrm{y}_i} e_i$, and $\mathsf{X}_ e = \lambda e$ because of the 1-homogeneity of $e$. We thus conclude that $\mathsf{X}$ is proportional to $\xi$ so that $N(\widetilde\partial_{ \xi \xi}^2 s) = {\mathbb R} \xi$.

The proof that $(ii)$ and $(iii)$ are equivalent is a consequence of (2.7) and of the identity

where $\mathsf{X}_ \mathrm{y} = (\mathsf{X}_1, \ldots, \mathsf{X}_ n)^t$. This identity may be obtained from (B.7) or by decomposing $\mathsf{X}_ \mathrm{y}$ in the form $\mathsf{X}_ \mathrm{y} = \alpha \mathrm{y} + z$ with $\langle \Lambda \mathrm{y}, z\rangle = 0$. □

C.   The simplified multicomponent model

We derive in this appendix the multicomponent flow model investigated in previous sections from a kinetic theory Cahn-Hilliard model ^[44]. The mixture fluid conservation Eqs (6.4)–(6.6) presented in Section 6.2 are also valid for Cahn-Hilliard type fluids so that only the transport fluxes need to be investigated. With the Cahn-Hilliard fluid model derived from the kinetic theory—indeed from the BBGKY hierarchy—the pressure tensor $\boldsymbol{\mathcal{P}}$ and the heat flux $\boldsymbol{\mathcal{Q}}$ are found in the form ^[44]

where $p^ {\rm{ex}}$ denotes the pressure, $\boldsymbol{I}$ the unit tensor, $\varkappa_{ij} = \varkappa_{ij}(T)$, $i, j\in \mathfrak{S}$, the species pair capillary coefficients that only depend on temperature, $\boldsymbol{\mathcal{P}}^{{\rm{d}}}$ the viscous tensor, and $\boldsymbol{q}$ the dissipative heat flux. The species pair capillarity coefficients $\varkappa_{ij}$ may be related to the species pair interaction potentials ^[44]. The dissipative fluxes $\boldsymbol{\mathcal{P}}^{{\rm{d}}}$, $\boldsymbol{\mathcal{J}}_i$, and $\boldsymbol{q}$ are given by

where $\mathfrak{v}$ denotes the bulk viscosity, $\eta$ the shear viscosity, and $L = (L_{ij})_{i, j\in \mathfrak{S}\cup\{ \mathrm{e}\}}$ the mass and heat diffusion matrix. These coefficients $\mathfrak{v}$, $\eta$, and $L$ are defined in the framework of the kinetic theory ^[44] and are not phenomenological coefficients as in the theory of irreversible processes.

The free energy $\mathcal{F}^ {\rm{ex}} = \mathcal{F}^ {}(\rho_1, \ldots, \rho_ n, T) + {\tfrac{1}{2}} \sum_{i, j\in \mathfrak{S}} \varkappa_{ij} \boldsymbol{\nabla}\mspace{-3mu}\rho_i {\boldsymbol{\cdot}} \boldsymbol{\nabla}\mspace{-3mu}\rho_j$ is of Van der Waals type, and the corresponding thermodynamic properties are found in the form

where $\mathcal{E}^ {\rm{ex}}$ is the energy per unit volume, $p^ {\rm{ex}}$ the pressure, $\mathcal{S}^ {\rm{ex}}$ the entropy per unit volume, and $g_i^ {\rm{ex}}$ the Gibbs function per unit mass of the $i$th species. The extended properties are denoted with the ${}^ {\rm{ex}}$ superscript whereas the corresponding classical properties that only depend on $(\rho_1, \ldots, \rho_ n, T)$ do not have any superscript. The extended Gibbs relation may further be written:

In order to simplify the full Cahn-Hilliard model derived in the framework of the kinetic theory, we assume that all the species pair capillarity coefficients are equal:

The governing equations may then be simplified with the relations $\sum_{i, j\in \mathfrak{S}} \varkappa_{ij} \boldsymbol{\nabla}\mspace{-3mu}\rho_i {\otimes} \boldsymbol{\nabla}\mspace{-3mu}\rho_j = \varkappa \boldsymbol{\nabla}\mspace{-3mu}\rho {\otimes} \boldsymbol{\nabla}\mspace{-3mu}\rho$, $\sum_{i, j\in \mathfrak{S}} \rho_i \boldsymbol{\nabla} {\boldsymbol{\cdot}}(\varkappa_{ij} \boldsymbol{\nabla}\mspace{-3mu}\rho_j) = \rho \boldsymbol{\nabla} {\boldsymbol{\cdot}}(\varkappa \boldsymbol{\nabla}\mspace{-3mu}\rho)$, and $\sum_{i, j\in \mathfrak{S}} \varkappa_{ij} \rho_i \boldsymbol{\nabla}\mspace{-3mu}\rho_j \boldsymbol{\nabla} {\boldsymbol{\cdot}} {\boldsymbol{u}} = \rho \boldsymbol{\nabla}\mspace{-3mu}\rho \boldsymbol{\nabla} {\boldsymbol{\cdot}} {\boldsymbol{u}}$, as well as the thermodynamic properties with $\sum_{i, j\in \mathfrak{S}} \varkappa_{ij} \boldsymbol{\nabla}\mspace{-3mu}\rho_i {\boldsymbol{\cdot}} \boldsymbol{\nabla}\mspace{-3mu}\rho_j = \varkappa \vert \boldsymbol{\nabla}\mspace{-3mu}\rho\vert^2$. Various terms may further be eliminated from the heat flux with $\sum_{i, j\in \mathfrak{S}} \varkappa_{ij} \boldsymbol{\nabla}\mspace{-3mu}\rho_j \, \boldsymbol{\nabla} {\boldsymbol{\cdot}} \boldsymbol{\mathcal{J}}_i = 0$ and $\sum_{i, j\in \mathfrak{S}} \boldsymbol{\nabla} {\boldsymbol{\cdot}}(\varkappa_{ij} \boldsymbol{\nabla}\mspace{-3mu}\rho_j) \boldsymbol{\mathcal{J}}_i = 0$ since $\sum_{i\in \mathfrak{S}} \boldsymbol{\mathcal{J}}_i = 0$. The Cahn-Hilliard type driving forces are also species independent, since $\sum_{l\in \mathfrak{S}} \varkappa_{il} \boldsymbol{\nabla}\mspace{-3mu}\rho_l = \varkappa \boldsymbol{\nabla}\mspace{-3mu}\rho$, and are thus in the null space of the mass and heat diffusion matrix $L$. The Cahn-Hilliard driving forces therefore disappear from the multicomponent fluxes, and it is finally obtained that

The thermodynamic properties also read

and (C.1)–(C.5) is the multicomponent model investigated in the previous sections.