# Multicomponent thermodynamics with instabilities and diffuse interfaces fluids

Vincent Giovangigli, Yoann Le Calvez, Guillaume Ribert

► To cite this version:

Vincent Giovangigli, Yoann Le Calvez, Guillaume Ribert. Multicomponent thermodynamics with instabilities and diffuse interfaces fluids. AIMS Mathematics, 2024, 9, pp.25979 - 26034. 10.3934/math.20241270. hal-04967684

HAL Id: hal-04967684

https://hal.science/hal-04967684v1

Submitted on 26 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Research article

# Multicomponent thermodynamics with instabilities and diffuse interfaces fluids

Vincent Giovangigli¹,*, Yoann Le Calvez¹ and Guillaume Ribert²

¹ CMAP, CNRS, École Polytechnique, 91128 Palaiseau, France

² CORIA, CNRS, Av. de l'Université, 76801 Saint-Etienne-du-Rouvray, France

* Correspondence: Email: vincent.giovangigli@polytechnique.edu.

**Abstract:** We investigated the mathematical structure of Gibbsian multicomponent thermodynamics with instabilities. We analyzed the construction of such thermodynamics from a pressure law using ideal gases as the low density limit. The fluid mixtures were allowed to have mechanically and chemically unstable states that were excluded in previous work on supercritical fluids, and the Soave-Redlich-Kwong cubic equation of state was specifically considered. We also investigated the mathematical structure of extended thermodynamics in the presence of cohesive forces—capillary effects—for a simplified diffuse interface fluid model. The thermodynamic formalism was validated by comparison with experimental data for mixtures of ethane and nitrogen. Very good agreement with experimental data was obtained for specific heats, multiphase equilibrium, and critical points, and we also analyzed the structure of strained jets of ethane.

**Keywords:** multicomponent thermodynamics; mathematical structure; extended thermodynamics; diffuse interface; strained flows

**Mathematics Subject Classification:** 80A17, 76B45, 35Q35

## 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–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, y_1, \dots, y_n)^t$ where $T$ is the absolute temperature, $v$ the volume per unit mass, $y_1, \dots, 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 $O_\zeta \subset \mathbb{R}^{2+n}$. We also address a different structure associated with the thermodynamic variable $Z = (T, \rho_1, \dots, \rho_n)^t$, where $\rho_1, \dots, \rho_n$ denotes the species masses per unit volume. Using the variable $\zeta = (T, v, y_1, \dots, y_n)^t$ leads to homogeneous thermodynamic functions—with singular entropy Hessians—and using the variable $Z = (T, \rho_1, \dots, \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–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 $V$ and normal $W$ variables. The entropic variable $V$ is closely associated with symmetrization properties for systems of conservation laws. However, the entropic variable $V$ is not practical in the presence of phase changes where it is more interesting to use a thermodynamic—part of a—normal variables $W$. We investigate conditions allowing the use of such a normal variable $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, y_1, \dots, y_n)$ and a quadratic expression of $p - p^{\text{id}}$ as $\rho$ goes to zero where $p^{\text{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 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-Redlich-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–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—
inter 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–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 $\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 $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–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–110] with applications to thermal diffusion [111].

The thermodynamic formalism is validated by comparison with experimental measurements for mixtures of ethane $C_2H_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_2H_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 III [124–126]. The numerical simulations confirm such a type III 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₂H₆-N₂ 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, y_1, \dots, 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, $y_1, \dots, y_n$ the species mass fractions, $n \ge 1$ the number of species, and $\mathfrak{S}$ the species indexing set $\mathfrak{S} = \{1, \dots, n\}$. The thermodynamic model is specified as $e(\zeta)$, $p(\zeta)$, and $s(\zeta)$, defined on an open set $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 $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 $id$ is associated with ideal gas mixtures thermodynamics.

We denote by $\partial$ the derivation operator with respect to the variable $\zeta = (T, v, y_1, \dots, y_n)^t$, by $d$ the total differential operator, and for any $\lambda > 0$, we define $\zeta_\lambda = (T, \lambda v, \lambda y_1, \dots, \lambda y_n)^t$. We denote by $y$ the mass fraction vector $y = (y_1, \dots, y_n)^t$, $\mathbb{I}$ the vector with unit components $\mathbb{I} = (1, \dots, 1)^t$, $\langle \cdot, \cdot \rangle$ the Euclidean scalar product, and define $\Sigma = \{y \in (0, \infty)^n; \langle y, \mathbb{I} \rangle = 1\}$. Physically, the mass fraction $y$ must sum up to unity $y \in \Sigma$. However, as usual in the modeling of multicomponent flows, we assume that the mass fractions are independent and the mass constraint $\langle y, \mathbb{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 $y_1$ is usually formally eliminated as discussed in Appendix A. We denote by $\gamma \in \mathbb{N}$, $\gamma \ge 2$, 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 $O_\zeta \subset (0, \infty)^{2+n}$. These functions are said to define a thermodynamics if properties $(T_0-T_2)$ hold.
(T_0) The functions e, p, and s are C^γ over a nonempty, connected, open set O_ζ ⊂ (0, ∞)^{2+n}. For any (λ, ζ) ∈ (0, ∞)×O_ζ, ζ_λ ∈ O_ζ, e(ζ_λ) = λe(ζ), p(ζ_λ) = p(ζ), and s(ζ_λ) = λs(ζ).

(T_1) For any $\zeta \in O_\zeta$, defining $g_k = \partial_{y_k} e - T \partial_{y_k} s$, $k \in \mathfrak{S}$, we have Gibbs' relation:

$$ T ds = \partial_T e \, dT + (p + \partial_v e) \, dv + \sum_{k \in \mathfrak{S}} (\partial_{y_k} e - g_k) \, dy_k. \quad (2.1) $$

(T_2) For any $(T, y_1, \dots, y_n)^t \in (0, \infty)^{1+n}$, there exists $v_m$ so that $v > v_m$ implies $(T, v, y_1, \dots, y_n)^t \in O_\zeta$ and

$$ \lim_{v \to \infty} (e - e^{id}) = 0, \quad \lim_{v \to \infty} v(p - p^{id}) = 0, \quad \lim_{v \to \infty} (s - s^{id}) = 0. \quad (2.2) $$

Property (T_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 $O_\zeta \subset (0, +\infty)^{2+n}$. Property (T_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 (T_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^{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 (T₀-T₁) and defining the mixture Gibbs function by $g = e + pv - Ts$, we have the relation $g = \sum_{k \in \mathfrak{S}} y_k g_k$.*

*Proof.* Since $e$ and $s$ are 1-homogeneous with respect to $(v, y_1, \dots, y_n)$, we obtain from Euler's relation that $e = v \partial_v e + \sum_{k \in \mathfrak{S}} y_k \partial_{y_k} e$ and $s = v \partial_v s + \sum_{k \in \mathfrak{S}} y_k \partial_{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_{y_k} s = \partial_{y_k} e - g_k$. Forming finally $g = e + pv - Ts$, we then obtain that $g = \sum_{k \in \mathfrak{S}} y_k g_k$. $\square$

## 2.2. The energy-based variable

We introduce the energy-based variable $\xi = (e, v, y_1, \dots, y_n)^t$ that may be used whenever the map $\zeta \mapsto \xi$ is at least a local diffeomorphism from a subset $O'_\zeta$ of $O_\zeta$ onto an open set $O'_\xi$, and denote by $\tilde{\partial}$ the derivation operator with respect to the variable $\xi = (e, v, y_1, \dots, 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^γ local diffeomorphism from an open set $O'_\zeta$ onto an open set $O'_\xi$. Then, Gibbs' relation is locally in the form*

$$ T ds = de + pdv - \sum_{k \in \mathfrak{S}} g_k dy_k, \quad (2.3) $$

and we have, in particular, $T \tilde{\partial}_e s = 1$, $T \tilde{\partial}_v s = p$, and $T \tilde{\partial}_{y_\nu} 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 $de = \partial_T e \,dT + \partial_v e \,dv + \sum_{k \in \mathcal{S}} \partial_{y_k} e \,dy_k$. $\square$

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^b$, $p^\sharp = p^b$, and $g_i^\sharp = g_i^b$ for all $i \in \mathcal{S}$ [1–3, 5]. We now give a sufficient conditions ($\mathcal{T}_3$), often met in practice, in such a way that the map $\zeta \mapsto \xi$ are a global diffeomorphism. The nonempty, open set $O_\zeta$ is said to be increasing with temperature if for any $(T, v, y_1, \dots, y_n)^t \in O_\zeta$ and any $T < T'$, we have $(T', v, y_1, \dots, y_n)^t \in O_\zeta$.

($\mathcal{T}_3$) For any $\zeta \in O_\zeta$, $\partial_T e(\zeta) > 0$, and $O_\zeta$ is increasing with temperature.

**Lemma 2.4.** *Assume that ($\mathcal{T}_0$) and ($\mathcal{T}_3$) hold, then the map $\zeta \mapsto \xi$ is a $C^\gamma$ global diffeomorphism from $O_\zeta$ onto a nonempty, connected, open set $O_\xi$.*

*Proof.* The map $\zeta \mapsto \xi$ is $C^\gamma$ and its Jacobian matrix is invertible over $O_\zeta$ since $\partial_T e > 0$. From the local inversion theorem, we know that the image $O_\xi$ of $O_\zeta$ is an open set. Moreover, the map $\zeta \mapsto \xi$ is one-to-one since if $(e^b, v^b, y_1^b, \dots, y_n^b) = (e^\sharp, v^\sharp, y_1^\sharp, \dots, y_n^\sharp)$ denoting $v = v^b = v^\sharp$ and $y_k = y_k^b = y_k^\sharp$ for $k \in \mathcal{S}$, then $e(T^b, v, y_1, \dots, y_n) = e(T^\sharp, v, y_1, \dots, y_n)$. The structure of $O_\zeta$ then guarantees that the segment $[\zeta^b, \zeta^\sharp]$ given by $(\alpha T^b + (1-\alpha)T^\sharp, v, y_1, \dots, y_n)^t$ for $0 \le \alpha \le 1$ lies in $O_\zeta$ so that

$$e^\sharp - e^b = \int_{T^b}^{T^\sharp} \partial_T e(\theta, v, y_1, \dots, y_n) \, d\theta = 0.$$

Since $\partial_T e$ is positive over $[\zeta^b, \zeta^\sharp]$, this relation implies that $T^b = T^\sharp$. The map $\zeta \mapsto \xi$ is thus a $C^\gamma$ global diffeomorphism from $O_\zeta$ onto an open set $O_\xi$ that is nonempty and connected since $O_\zeta$ is nonempty and connected. $\square$

Under assumptions ($\mathcal{T}_0$) and ($\mathcal{T}_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 $\tilde{\partial}_{\xi\xi}^2 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 $y = (y_1, \dots, y_n)^t$ and the matrix $\Lambda = (\Lambda_{kl})_{k,l \in \mathcal{S}}$ of size $n$ is defined by

$$\Lambda_{kl} = \frac{\partial_{y_k} g_l}{T} = \frac{\partial_{y_l} g_k}{T}, \quad k, l \in \mathcal{S}, \qquad (2.4)$$

as well as the matrix $\tilde{\Lambda}$ defined when $\langle \Lambda y, y \rangle \neq 0$ by

$$\tilde{\Lambda} = \Lambda - \frac{\Lambda y \otimes \Lambda y}{\langle \Lambda y, y \rangle}. \qquad (2.5)$$

We first establish a few useful properties of the stability matrix $\Lambda$.
**Proposition 2.5.** *Assuming that (T₀-T₁) are satisfied, for any ζ ∈ O_ζ, we have the relations*

$$ \Lambda y = \frac{v}{T}(\partial_{y_1} p, \dots, \partial_{y_n} p)^t, \qquad (2.6) $$

$$ \langle \Lambda y, y \rangle = -\frac{v^2}{T} \partial_v p. \qquad (2.7) $$

*Proof.* From Gibbs' relation, we have $T\partial_T s = \partial_T e$, $T\partial_v s = \partial_v e + p$, $T\partial_{y_k} s = \partial_{y_k} e - g_k$, and this implies the compatibility relations

$$ \partial_T \left( \frac{p}{T} \right) = \frac{\partial_v e}{T^2}, \quad \partial_T \left( \frac{-g_k}{T} \right) = \frac{\partial_{y_k} e}{T^2}, \quad -\partial_v g_k = \partial_{y_k} p. \qquad (2.8) $$

The matrix $\Lambda$ is first symmetric since $\partial_{y_l} g_k = \partial_{y_k}^2 e - T \partial_{y_k}^2 s$, and in order to establish (2.6) and (2.7), we note that $g_k(\zeta)$ is 0-homogeneous in $(v, y_1, \dots, y_n)$ since $g_k = \partial_{y_k} e - T \partial_{y_k} s$, $T$ is 0-homogeneous, and $e$ and $s$ are 1-homogeneous. This implies the relation $\sum_{l \in \mathfrak{Z}} y_l \partial_{y_l} g_k = -v \partial_v g_k$ and, thus, $\sum_{k \in \mathfrak{Z}} \Lambda_{kl} y_l = -v \partial_v g_k / T$, and (2.6) is established by using (2.8). Upon multiplying (2.6) by $y_k$ and summing over $k$, we next obtain that $T \langle \Lambda y, y \rangle = -v \sum_{k \in \mathfrak{Z}} y_k \partial_v g_k$. However, from the compatibility relation (2.8), we get $\partial_v g_k = -\partial_{y_k} p$, and thanks to the 0-homogeneity of $p(\zeta)$ with respect to $(v, y_1, \dots, y_n)$, we conclude that $\sum_{k \in \mathfrak{Z}} y_k \partial_v g_k = - \sum_{k \in \mathfrak{Z}} y_k \partial_{y_k} p = v \partial_v p$ and (2.7) is established. $\square$

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-Ts$, we note that $g_k = \partial_{y_k} f$ as well as $\Lambda_{kl} = \partial_{y_k y_l}^2 f / T$ and $\Lambda = \partial_{yy}^2 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 (T₀-T₁) are satisfied and that $\zeta \mapsto \xi$ is a C^γ local diffeomorphism, the following statements are equivalent :*

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

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

(iii) $\partial_T e > 0$, $\partial_v p < 0$, and $\tilde{\Lambda}$ is positive semi-definite with null space $\mathbb{R}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 $\tilde{\Lambda}$ is positive semi-definite with null space $N(\tilde{\Lambda}) = \mathbb{R}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, y_1, \dots, y_n)^t$ as well as molar properties. We denote by $\tilde{\partial}$ the derivation operator with respect to the variable $\pi = (T, p, y_1, \dots, y_n)^t$.

**Proposition 2.7.** *Assume that (T₀-T₁) hold and that $\zeta \mapsto \pi$ is a local C^γ diffeomorphism from $O'_\zeta \subset O_\zeta$ onto an open set $O_\pi$. Then, we have the differential relations $g_k = \tilde{\partial}_{y_k} g$, $k \in \mathfrak{Z}$, and Gibbs' relation may be written:*

$$ dg = -sdT + vdp + \sum_{k \in \mathfrak{Z}} g_k dy_k. \qquad (2.9) $$
*Proof.* For any smooth function $\chi$ of $\zeta$ or $\pi$, we have the following differential relations associated with the variables $\pi$ and $\zeta$:

$$ \widehat{\partial}_{T}\chi = \partial_{T}\chi + \partial_{v}\chi \widehat{\partial}_{T}v, \qquad \partial_{T}\chi = \widehat{\partial}_{T}\chi + \widehat{\partial}_{p}\chi \partial_{T}p, \qquad (2.10) $$

$$ \widehat{\partial}_{p}\chi = \partial_{v}\chi \widehat{\partial}_{p}v, \qquad \partial_{v}\chi = \widehat{\partial}_{p}\chi \partial_{v}p, \qquad (2.11) $$

$$ \widehat{\partial}_{y_k}\chi = \partial_{y_k}\chi + \partial_{v_k}\chi \widehat{\partial}_{y_k}v, \qquad \partial_{y_k}\chi = \widehat{\partial}_{y_k}\chi + \widehat{\partial}_{p_k}\chi \partial_{y_k}p. \qquad (2.12) $$

Letting $g = e + pv - Ts$, 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 $(\Gamma_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 $(\Gamma_1)$, we have established that $\widehat{\partial}_{p}g = v$. Finally, we have $\widehat{\partial}_{y_k}g = \widehat{\partial}_{y_k}e + p\widehat{\partial}_{y_k}v - T\widehat{\partial}_{y_k}s$, and from (2.12), the relations $\widehat{\partial}_{y_k}e = \partial_{y_k}e + \partial_{v_k}e \widehat{\partial}_{y_k}v$ and $\widehat{\partial}_{y_k}s = \partial_{y_k}s + \partial_{v_k}s \widehat{\partial}_{y_k}v$, and after some algebra, we obtain so that $\widehat{\partial}_{y_k}g = g_k$ using also $(\Gamma_1)$. In conclusion, we have established that $\widehat{\partial}_{T}g = -s$, $\widehat{\partial}_{p}g = v$, and $\widehat{\partial}_{y_k}g = g_k$, and (2.9) then follows. $\square$

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

**Proposition 2.8.** Assume that $(\Gamma_0-\Gamma_1)$ hold and that $\zeta \mapsto \pi$ is a local $C^\gamma$ diffeomorphism from $O'_\zeta \subset O_\zeta$ onto an open set $O'_\pi$. Then, for any $\zeta \in O'_\zeta$, we have the matrix identity

$$ \widehat{\partial}_{yy}^2 g = T \left( \Lambda - \frac{\Lambda y \otimes \Lambda y}{\langle \Lambda y, y \rangle} \right). \qquad (2.13) $$

*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}_{y_k}g = g_k$, and this implies the compatibility relations

$$ \widehat{\partial}_{T}v = -\widehat{\partial}_{p}s, \qquad \widehat{\partial}_{T}g_k = -\widehat{\partial}_{y_k}s, \qquad \widehat{\partial}_{p}g_k = \widehat{\partial}_{y_k}v. \qquad (2.14) $$

Upon letting $\chi = g_l$ in (2.12), we obtain that $\widehat{\partial}_{y_k y_l}^2 g = \widehat{\partial}_{y_k}g_l = \widehat{\partial}_{y_l}g_k$ can be written $\widehat{\partial}_{y_k}g_l = \partial_{y_k}g_l - \widehat{\partial}_{p}g_l \partial_{y_k}p$ so that from (2.8) and (2.11), $\widehat{\partial}_{y_k}g_l = \partial_{y_k}g_l + \partial_{v_k}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). $\square$

Thermodynamic stability thus involves the matrix $\widehat{\partial}_{\zeta\zeta}^2 s$, mechanical and chemical stability of the matrix $\Lambda = \partial_{yy}^2 f$, and chemical stability of the matrix $\widehat{\partial}_{yy}^2 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 $O_\zeta$ is said to be increasing with volume if for any $(T, v, y_1, \dots, y_n)^t \in O_\zeta$ and any $v < v'$, we have $(T, v', y_1, \dots, y_n)^t \in O_\zeta$.

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

*Proof.* The proof is similar to that of Lemma 2.4. $\square$
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 $x_1, \dots, x_n$ are defined by

$$x_i = \frac{y_i m}{m_i}, \quad i \in \mathcal{S}, \quad \frac{\sum_{i \in \mathcal{S}} y_i}{m} = \sum_{i \in \mathcal{S}} \frac{y_i}{m_i}, \qquad (2.15)$$

where $m$ is the molar mass of the mixture. The factor $\sum_{i \in \mathcal{S}} 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 \mathcal{S}} y_i = \sum_{i \in \mathcal{S}} x_i$ and $(\sum_{i \in \mathcal{S}} x_i)m = \sum_{i \in \mathcal{S}} x_i m_i$. As a typical example, the mixture Gibbs function per unit mole $G$ and the species Gibbs function per unit mole $G_i$, $i \in \mathcal{S}$, may then be defined as

$$G = mg, \quad G_i = m_i g_i, \quad i \in \mathcal{S}. \qquad (2.16)$$

From the relation $g = \sum_{i \in \mathcal{S}} y_i g_i$, it is easily obtained that $G = \sum_{i \in \mathcal{S}} x_i G_i$. All thermodynamic properties per unit mass may similarly be rewritten per unit mole and the molar variables $(E, V, x_1, \dots, x_n)^t$, $(T, V, x_1, \dots, x_n)^t$, or $(T, p, x_1, \dots, x_n)^t$ may also be used, instead of $\xi, \zeta$, or $\pi$, where $E$ denotes the energy per unit mole and $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 $id$ and involved in Property (T$_2$). In terms of the variable $\zeta = (T, V, y_1, \dots, y_n)^t$, the pressure of an ideal gas mixture is given by

$$p^{id} = \frac{RT}{V} \sum_{k \in \mathcal{S}} \frac{y_k}{m_k}, \qquad (2.17)$$

and the internal energy by

$$e^{id} = \sum_{k \in \mathcal{S}} y_k e_k^{id}, \quad e_k^{id} = e_k^{st} + \int_{T}^{T} c_{vk}^{id}(\theta) d\theta, \qquad (2.18)$$

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

$$s^{id} = \sum_{k \in \mathcal{S}} y_k s_k^{id}, \quad s_k^{id} = s_k^{st} + \int_{T}^{T} \frac{c_{vk}^{id}(\theta)}{\theta} d\theta - \frac{R}{m_k} \log \frac{y_k}{vm_k \gamma^{st}}, \qquad (2.19)$$

where $s_k^{st}$ is the formation entropy of the $k$th species per unit mass at the standard temperature $T^{st}$ and standard pressure $p^{st}$, and $\gamma^{st} = p^{st}/RT^{st}$ is the standard concentration. The enthalpy $h^{id} = e^{id} + p^{id}/V$, Gibbs function $g^{id} = e^{id} + p^{id}/V - T s^{id}$, and free energy $f^{id} = e^{id} - T s^{id}$ are then easily evaluated. The ideal gas specific heats, formation energies, and formation entropies satisfy the following assumptions.
(ID) The formation energies $e_i^{st}$, $i \in \mathfrak{S}$, and entropies $s_i^{st}$, $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_{vi}^{id}$, $i \in \mathfrak{S}$, are $C^\gamma$ functions over $[0, \infty)$, and there exist constants $c_v$ and $\bar{c}_v$ such that $0 < c_v \le c_{vi}^{id} \le \bar{c}_v$ for all $T \ge 0$ and $i \in \mathfrak{S}$.

We assume throughout the paper that Property (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 $(T_3)$ is trivial in this situation.

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

$$O_{\xi}^{id} = \{\xi = (\xi_e, \xi_v, \xi_1, \dots, \xi_n)^t \text{ with } \xi_v > 0, \xi_1 > 0, \dots, \xi_n > 0, \xi_e > \sum_{i \in \mathfrak{S}} \xi_i e_i^0\},$$

where $e_i^0$ denotes the energy of the $i$th species at zero temperature $e_k^0 = e_k^{st} - \int_0^{T^{st}} c_{vk}^{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–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 y, \mathbf{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 $Z = (T, \rho_1, \dots, \rho_n)^t$ is used for convenience where $\rho_1, \dots, \rho_n$ denote the species masses per unit volume, also termed species partial densities as well as $X = (\mathcal{E}, \rho_1, \dots, \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 $Z$. We denote by $\mathfrak{D}$ the derivation operator with respect to the variable $Z$ and by $\mathfrak{D}$ the derivation operator with respect to the variable $X$. We remind that $\gamma \in \mathbb{N}$, $\gamma \ge 2$, 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, y_1, \dots, y_n)^t$ and investigate the definition and properties of thermodynamic functions in terms of the variable $Z = (T, \rho_1, \dots, \rho_n)^t$ involving the species partial densities.

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

$$ \mathcal{E}(T, \rho_1, \dots, \rho_n) = \left( \sum_{i \in \mathfrak{Z}} \rho_i \right) e \left( T, \frac{1}{\sum_{i \in \mathfrak{Z}} \rho_i}, \frac{\rho_1}{\sum_{i \in \mathfrak{Z}} \rho_i}, \dots, \frac{\rho_n}{\sum_{i \in \mathfrak{Z}} \rho_i} \right), \quad (3.1) $$

$$ \mathcal{P}(T, \rho_1, \dots, \rho_n) = p \left( T, \frac{1}{\sum_{i \in \mathfrak{Z}} \rho_i}, \frac{\rho_1}{\sum_{i \in \mathfrak{Z}} \rho_i}, \dots, \frac{\rho_n}{\sum_{i \in \mathfrak{Z}} \rho_i} \right), \quad (3.2) $$

$$ \mathcal{S}(T, \rho_1, \dots, \rho_n) = \left( \sum_{i \in \mathfrak{Z}} \rho_i \right) s \left( T, \frac{1}{\sum_{i \in \mathfrak{Z}} \rho_i}, \frac{\rho_1}{\sum_{i \in \mathfrak{Z}} \rho_i}, \dots, \frac{\rho_n}{\sum_{i \in \mathfrak{Z}} \rho_i} \right), \quad (3.3) $$

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

$$ O_Z = \{ (T, \rho_1, \dots, \rho_n) \in (0, \infty)^{1+n} ; (T, \frac{1}{\sum_{i \in \mathfrak{Z}} \rho_i}, \frac{\rho_1}{\sum_{i \in \mathfrak{Z}} \rho_i}, \dots, \frac{\rho_n}{\sum_{i \in \mathfrak{Z}} \rho_i}) \in O_{\zeta} \}. $$

Then $\mathcal{E}, \mathcal{P}$, and $\mathcal{S}$ satisfy Properties $(\mathcal{T}_0-\mathcal{T}_2)$.

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

(T1) For any $Z \in O_Z$, defining $g_k = \mathfrak{D}_{\rho_k} \mathcal{E} - T \mathfrak{D}_{\rho_k} \mathcal{S}$, $k \in \mathfrak{Z}$, we have Gibbs' volumetric relation

$$ T d\mathcal{S} = \mathfrak{D}_T \mathcal{E} dT + \sum_{k \in \mathfrak{Z}} (\mathfrak{D}_{\rho_k} \mathcal{E} - g_k) d\rho_k, \quad (3.4) $$

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

(T2) For any $(T, y_1, \dots, y_n) \in (0, \infty)^{1+n}$, such that $\sum_{i \in \mathfrak{Z}} y_i = 1$, there exists a threshold density $\rho_m > 0$ such that $(T, \rho y_1, \dots, \rho y_n)^t \in O_Z$ for $0 < \rho < \rho_m$ and

$$ \lim_{\rho \to 0} \frac{1}{\rho}(\mathcal{E} - \mathcal{E}^{\text{id}}) = 0, \quad \lim_{\rho \to 0} \frac{1}{\rho}(\mathcal{P} - \mathcal{P}^{\text{id}}) = 0, \quad \lim_{\rho \to 0} \frac{1}{\rho}(\mathcal{S} - \mathcal{S}^{\text{id}}) = 0. \quad (3.5) $$

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

From $(\mathcal{T}_0)$ and the 1-homogeneity of energy with respect to $(v, y_1, \dots, y_n)$, we obtain, after some algebra, that

$$ d\mathcal{E} = \rho \partial_T e dT + \sum_{k \in \mathfrak{Z}} \partial_{y_k} e d\rho_k, \quad (3.6) $$

and, similarly,

$$ d\mathcal{S} = \rho \partial_T s dT + \sum_{k \in \mathfrak{Z}} \partial_{y_k} s d\rho_k. \quad (3.7) $$
This implies, in particular, that $\mathfrak{d}_T\mathcal{E} = \rho\partial_T e$, $\mathfrak{d}_{\rho_k}\mathcal{E} = \partial_{y_k}e$, $\mathfrak{d}_T S = \rho\partial_T s$, and $\mathfrak{d}_{\rho_k}S = \partial_{y_k}s$ in such a way that $g_k = \mathfrak{d}_{\rho_k}\mathcal{E} - T\mathfrak{d}_{\rho_k}S = \partial_{y_k}e - T\partial_{y_k}s$ for $k \in \mathfrak{Z}$. Then, thanks to $(\mathcal{T}_1)$, since $T\partial_T s = \partial_T e$ and $T\partial_{y_k}s = \partial_{y_k}e - g_k$, we obtain $T\mathfrak{d}_T S = \mathfrak{d}_T\mathcal{E}$ and $T\mathfrak{d}_{\rho_k}S = \mathfrak{d}_{\rho_k}\mathcal{E} - g_k$ for $k \in \mathfrak{Z}$, and this now implies Gibbs' relation for volumetric densities (3.4). In addition, the relation $\sum_{k \in \mathfrak{Z}} \rho_k g_k = \mathcal{E} + \mathcal{P} - TS$ is a consequence of Lemma 2.2 and $(\mathcal{T}_1)$ is established.

Furthermore, $(\mathcal{T}_2)$ is a simple rewriting of $(\mathcal{T}_2)$ with $e, p$, and $s$ evaluated at $T, v, y_1, \dots, y_n$, and the proof is complete. $\square$

The functions $\mathcal{E}, \mathcal{P}$, and $S$ are said to define a thermodynamics if Properties $(\mathcal{T}_0-\mathcal{T}_2)$ hold. Note that there are not anymore homogeneity properties in terms of the variable $Z$ in $(\mathcal{T}_0)$. There is also a constraint associated with Gibbs' relation (3.4) in terms of $S, \mathcal{E}$, and $\rho_1, \dots, \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 given by $\mathcal{F} = \mathcal{E} - TS$ and the Gibbs function per unit volume by $\mathcal{G} = \mathcal{F} + p = \mathcal{E} + p - TS$.

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

($\mathcal{T}_3$) For any $Z \in O_Z$, $\mathfrak{d}_T\mathcal{E}(Z) > 0$, $O_Z$ is increasing with temperature.

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

$$ T dS = d\mathcal{E} - \sum_{k \in \mathfrak{Z}} g_k d\rho_k. \qquad (3.8) $$

*Proof.* The proof is similar to that of the homogeneous case. $\square$

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

$$ \Lambda_{kl} = \frac{\rho \mathfrak{d}_{\rho_k} g_l}{T} = \frac{\rho \mathfrak{d}_{\rho_l} g_k}{T}, \qquad (3.9) $$

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}_0-\mathcal{T}_2)$ are satisfied and that the map $Z \mapsto X$ is a $C^\gamma$ local diffeomorphism from an open set $O'_Z$ onto an open set $O'_X$. Then, for any $Z \in O'_Z$, the following statements are equivalent:

(i) $\tilde{\mathfrak{d}}_{XX'}^2 S$ is negative definite.

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

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

$$ \tilde{\mathfrak{d}}_{\mathcal{E}}\chi = \mathfrak{d}_T\chi \tilde{\mathfrak{d}}_{\mathcal{E}}T, \qquad \mathfrak{d}_T\chi = \tilde{\mathfrak{d}}_{\mathcal{E}}\chi \mathfrak{d}_T\mathcal{E}, \qquad (3.10) $$
$$ \tilde{\mathfrak{D}}_{\rho_k} \chi = \mathfrak{D}_{\rho_k} \chi + \mathfrak{D}_T \chi \tilde{\mathfrak{D}}_{\rho_k} T, \quad \mathfrak{D}_{\rho_k} \chi = \tilde{\mathfrak{D}}_{\rho_k} \chi + \tilde{\mathfrak{D}}_{\mathcal{E}} \chi \mathfrak{D}_{\rho_k} \mathcal{E}. \qquad (3.11) $$

We now evaluate the Hessian matrix of $S$. We first have $\tilde{\mathfrak{D}}_{\mathcal{E}\mathcal{E}}^2 S = \tilde{\mathfrak{D}}_{\mathcal{E}}(\frac{1}{T})$ so that $\tilde{\mathfrak{D}}_{\mathcal{E}\mathcal{E}}^2 S = -\tilde{\mathfrak{D}}_{\mathcal{E}} T/T^2$. Similarly, we have $\tilde{\mathfrak{D}}_{\mathcal{E}\rho_k}^2 S = \tilde{\mathfrak{D}}_{\rho_k}(\frac{1}{T})$ so that $\tilde{\mathfrak{D}}_{\mathcal{E}\rho_k}^2 S = -\tilde{\mathfrak{D}}_{\rho_k} T/T^2$. Letting $\chi = T$ in (3.11), we get $\tilde{\mathfrak{D}}_{\rho_k} T = -\tilde{\mathfrak{D}}_{\mathcal{E}} T \mathfrak{D}_{\rho_k} \mathcal{E}$, and since $\mathfrak{D}_{\rho_k} \mathcal{E} = \partial_{y_k} e$, we have established that $\tilde{\mathfrak{D}}_{\mathcal{E}\rho_k}^2 S = \tilde{\mathfrak{D}}_{\mathcal{E}} T \partial_{y_k} e/T^2$. Moreover, with $\tilde{\mathfrak{D}}_{\rho_k \rho_l}^2 S = -\tilde{\mathfrak{D}}_{\rho_k}(\frac{g_l}{T}) = -\tilde{\mathfrak{D}}_{\rho_k}(\frac{g_l}{T})$ and (3.11), we deduce that $\tilde{\mathfrak{D}}_{\rho_k \rho_l}^2 S = -\tilde{\mathfrak{D}}_{\rho_k}(\frac{g_l}{T}) - \tilde{\mathfrak{D}}_T(\frac{g_l}{T}) \tilde{\mathfrak{D}}_{\rho_k} T$ so that $\tilde{\mathfrak{D}}_{\rho_k \rho_l}^2 S = -\frac{\partial_{y_k} g_l}{\rho T} - \frac{\tilde{\mathfrak{D}}_{\mathcal{E}} T \partial_{y_k} e \partial_{y_l} e}{T^2}$, $k, l \in \mathfrak{S}$. We have thus established that

$$ \tilde{\mathfrak{D}}_{\mathcal{E}\mathcal{E}}^2 S = -\frac{\tilde{\mathfrak{D}}_{\mathcal{E}} T}{T^2}, \quad \tilde{\mathfrak{D}}_{\mathcal{E}\rho_k}^2 S = \frac{\tilde{\mathfrak{D}}_{\mathcal{E}} T \partial_{y_k} e}{T^2}, \quad k \in \mathfrak{S}, \qquad (3.12) $$

$$ \tilde{\mathfrak{D}}_{\rho_k \rho_l}^2 S = -\frac{\partial_{y_k} g_l}{\rho T} - \frac{\tilde{\mathfrak{D}}_{\mathcal{E}} T \partial_{y_k} e \partial_{y_l} e}{T^2}, \quad k, l \in \mathfrak{S}. \qquad (3.13) $$

Letting $\mathfrak{f}^\mathcal{E} = (1, -\partial_{y_1} e, \dots, -\partial_{y_1} e)^t$, we obtain that for any $x = (x_\mathcal{E}, x_1, \dots, x_n)^t \in \mathbb{R}^{1+n}$

$$ \langle (\tilde{\mathfrak{D}}_{x_k}^2 S) x, x \rangle = -\tilde{\mathfrak{D}}_{\mathcal{E}} T \frac{\langle \mathfrak{f}^\mathcal{E}, x \rangle^2}{T^2} - \frac{1}{\rho} \langle \Lambda x_y, x_y \rangle, \qquad (3.14) $$

where $\Lambda_{kl} = \partial_{y_k} g_l/T$, $k, l \in \mathfrak{S}$ and $x_y = (x_1, \dots, x_n)^t$.

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

**Remark 3.4.** The functions $\mathcal{E}$, $\mathcal{P}$, and $S$ may also be defined by $\mathcal{E}(T, \rho_1, \dots, \rho_n) = e(T, 1, \rho_1, \dots, \rho_n)$, $\mathcal{P}(T, \rho_1, \dots, \rho_n) = p(T, 1, \rho_1, \dots, \rho_n)$, and $S(T, \rho_1, \dots, \rho_n) = s(T, 1, \rho_1, \dots, \rho_n)$, and the open set $O_z$ by $O_z = \{ (T, \rho_1, \dots, \rho_n) \in (0, \infty)^{1+n}; (T, 1, \rho_1, \dots, \rho_n) \in O_\zeta \}$ thanks to the homogeneity properties of $e$, $s$, $p$, and $O_\zeta$.

## 3.2. From nonhomogeneous to homogeneous thermodynamics

We assume that a thermodynamics is available in terms of the variable $Z = (T, \rho_1, \dots, \rho_n)^t$ and we investigate the definition and properties of thermodynamic functions in terms of the natural variable $\zeta = (T, v, y_1, \dots, y_n)^t$.

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

$$ e(T, v, y_1, \dots, y_n) = v \mathcal{E} \left( T, \frac{y_1}{v}, \dots, \frac{y_n}{v} \right), \qquad (3.15) $$

$$ p(T, v, y_1, \dots, y_n) = \mathcal{P} \left( T, \frac{y_1}{v}, \dots, \frac{y_n}{v} \right), \qquad (3.16) $$

$$ s(T, v, y_1, \dots, y_n) = v S \left( T, \frac{y_1}{v}, \dots, \frac{y_n}{v} \right), \qquad (3.17) $$

and defined over the open set

$$ O_\zeta = \{ (T, v, y_1, \dots, y_n) \in (0, \infty)^{2+n}; (T, \frac{y_1}{v}, \dots, \frac{y_n}{v}) \in O_z \}. \qquad (3.18) $$

Then $e$, $p$, and $s$ are $C^\gamma$ functions of the variable $\zeta = (T, v, y_1, \dots, y_n)^t$ over $O_\zeta \subset (0, \infty)^{2+n}$ and satisfy Properties $(\mathcal{T}_0-\mathcal{T}_2)$.
*Proof.* The fact that $O_ζ$ is open, nonempty, and connected and that $e, p$, and $s$ are smooth is straightforward to establish. The homogeneity properties of $O_ζ$, $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

$$ds = v \partial_T S dT + (S - \sum_{k \in \mathcal{S}} \rho_k \partial_{\rho_k} S) dv + \sum_{k \in \mathcal{S}} \partial_{\rho_k} S dy_k,$$

$$de = v \partial_T \mathcal{E} dT + (\mathcal{E} - \sum_{k \in \mathcal{S}} \rho_k \partial_{\rho_k} \mathcal{E}) dv + \sum_{k \in \mathcal{S}} \partial_{\rho_k} \mathcal{E} dy_k.$$

However, from Gibbs’ relation in ($\mathcal{T}_2$), we obtain that $T \partial_T S = \partial_T \mathcal{E}$. From the definition of $g_k = \partial_{\rho_k} \mathcal{E} - T \partial_{\rho_k} S$, $k \in \mathcal{S}$ and from the natural constraint $\sum_{k \in \mathcal{S}} \rho_k g_k = \mathcal{E} + \mathcal{P} - T S$, we obtain that

$$T(S - \sum_{k \in \mathcal{S}} \rho_k \partial_{\rho_k} S_k) - (\mathcal{E} - \sum_{k \in \mathcal{S}} \rho_k \partial_{\rho_k} \mathcal{E}_k) = \mathcal{P}, \quad T \partial_{\rho_k} S - \partial_{\rho_k} \mathcal{E} = -g_k,$$

so in forming $Tds - de$, we directly obtain $(\mathcal{T}_1)$.

Finally, $(\mathcal{T}_2)$ is a simple transformation of $(\mathcal{T}_2)$ and the proof is complete. $\square$

When thermal stability holds with respect to the thermodynamic in the variable $Z$, it then holds with respect to the variable $ζ$.

**Proposition 3.6.** *Assuming Property (T$_3$), Property (T$_3$) holds.*

*Proof.* The proof is straightforward since $\partial_T \mathcal{E} = \rho \partial_T e$. $\square$

### 3.3. Entropic and normal variables

Entropic variables are naturally involved in symmetrization methods for systems of conservation laws [10–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–19].

The thermodynamic—part of the—entropic variable is defined by $\mathcal{V} = -(\partial_\chi S)^t$ and easily evaluated to be

$$\mathcal{V} = -(\tilde{\partial}_\chi S)^t = \frac{1}{T} (-1, g_1, \dots, g_n)^t. \qquad (3.19)$$

The mathematical entropy $\sigma$ is defined by $\sigma = -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–16]. The entropic variable $\mathcal{V} = (\partial_\chi \sigma)^t$ is associated with symmetrization methods and

$$\tilde{A}_0 = \tilde{\partial}_\chi \mathcal{V} = \tilde{\partial}_{\chi\chi}^2 \sigma,$$

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–19].
A major difficulty associated with the entropic variable, however, is that the map $\chi \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 $\chi^b \neq \chi^\sharp$ indeed implies that $\mathcal{V}^b = \mathcal{V}^\sharp$ and the map $\chi \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

$$ \mathcal{W} = (T, \rho, \frac{g_2 - g_1}{T}, \dots, \frac{g_n - g_1}{T})^t, \qquad (3.20) $$

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

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

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

**Proposition 3.7.** Assume that ($\mathcal{T}_0-\mathcal{T}_2$) and ($\mathcal{T}_4$) hold. Then, the matrix $\Lambda$ is positive definite on the hyperplane $\mathbb{I}^\perp$, and the following direct sum holds:

$$ \mathbb{R}\mathbb{I} \oplus \Lambda(\mathbb{I}^\perp) = \mathbb{R}^n. \qquad (3.21) $$

*Proof.* Since the map $\mathcal{Z} \mapsto \mathcal{W}$ is a $C^{\gamma-1}$ diffeomorphism, the Jacobian matrix $\mathcal{D}_\mathcal{Z}\mathcal{W}$ is invertible on the open set $O_\mathcal{Z}$. Keeping in mind that $\mathcal{Z} = (T, \rho_1, \dots, \rho_n)^t$, the Jacobian matrix $\mathcal{D}_\mathcal{Z}\mathcal{W}$ is of size $n+1$ and in the form

$$ \mathcal{D}_\mathcal{Z}\mathcal{W} = \begin{pmatrix} 1 & 0_{1,n} \\ *_{n,1} & \mathcal{J} \end{pmatrix}, \qquad \mathcal{J} = \begin{pmatrix} 1 & 1 & \cdots & 1 \\ \partial_{\rho_1}\left(\frac{g_2-g_1}{T}\right) & \partial_{\rho_2}\left(\frac{g_2-g_1}{T}\right) & \cdots & \partial_{\rho_n}\left(\frac{g_2-g_1}{T}\right) \\ \vdots & \vdots & & \vdots \\ \partial_{\rho_1}\left(\frac{g_n-g_1}{T}\right) & \partial_{\rho_2}\left(\frac{g_n-g_1}{T}\right) & \cdots & \partial_{\rho_n}\left(\frac{g_n-g_1}{T}\right) \end{pmatrix}, $$

where $0_{i,j}$ denotes a zero block with $i$ lines and $j$ columns, $*_{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_1-g_1}{T})$ from the $i$th line of $\mathcal{J}$, and denoting by $e_i$, $i \in \mathfrak{S}$, the canonical base vectors of $\mathbb{R}^n$, we have

$$ \mathcal{J} - \sum_{2 \le i \le n} \partial_{\rho_1} \left( \frac{g_i - g_1}{T} \right) e_i \otimes \mathbb{I} = \begin{pmatrix} 1 & 1 & \cdots & 1 \\ 0_{n-1,1} & \mathcal{J}' & & \end{pmatrix}, $$

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 $O_\mathcal{Z}$. Letting then $v_i = e_i - e_1$, for $2 \le i \le n$, we note that the vectors $v_i$, for $2 \le i \le n$, are spanning the hyperplane $\mathbb{I}^\perp$ and that

$$ \langle \Lambda v_i, v_j \rangle = \Lambda_{i,j} - \Lambda_{1,j} - \Lambda_{i,1} + \Lambda_{1,1}, \quad 2 \le i, j \le n, $$
since $(\Lambda 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

$$ \mathcal{J}' = \rho \begin{pmatrix} \langle \Lambda v_2, v_2 \rangle & \cdots & \langle \Lambda v_2, v_n \rangle \\ \vdots & \ddots & \vdots \\ \langle \Lambda v_n, v_2 \rangle & \cdots & \langle \Lambda v_n, v_n \rangle \end{pmatrix}, $$

and is thus a Gramm matrix up to the $\rho$ scalar factor that remains invertible over $O_Z$. This implies that the quadratic form associated with $\Lambda$ remains definite on $\mathbb{I}^{\perp}$ over the open set $O_Z$. Since this quadratic form is positive definite for small $\rho$, and since $O_Z$ is connected, we deduce that $\Lambda$ remains positive definite on the hyperplane $\mathbb{I}^{\perp}$ over $O_Z$. Since $\mathcal{J}$ remains nonsingular and since $(\Lambda 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. $\square$

An important consequence of $(\mathcal{T}_4)$ is that the matrix $\Lambda$ has at most one isolated negative eigenvalue. That is, when $(\mathcal{T}_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}_0-\mathcal{T}_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₂-N₂ mixtures, O₂-H₂O mixtures, and C₂H₆-N₂ 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}_4)$, assuming that the slices of $O_Z$ along $\mathbb{I}^{\perp}$ are convex.

**Proposition 3.9.** Assume that $(\mathcal{T}_0-\mathcal{T}_2)$ hold, that the matrix $\Lambda$ is positive definite on the hyperplane $\mathbb{I}^{\perp}$, that for any $Z \in O_Z$ the set $O_Z \cap (Z + \mathbb{I}^{\perp})$ is convex, and that the direct sum (3.21) holds. Then, Property $(\mathcal{T}_4)$ holds.

*Proof.* Since $\Lambda$ is positive definite on the hyperplane $\mathbb{I}^{\perp}$, and since we have the direct sum (3.21), we conclude that the Jacobian matrix $\delta_Z W$ is invertible. The map $Z \mapsto 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 $Z^b$ and $Z^\sharp$ are such that $W^b = W^\sharp$, we first deduce from (3.20) that $T^b = T^\sharp$ and $\rho^b = \rho^\sharp$ so that $Z^b - Z^\sharp \in \mathbb{I}^{\perp}$. Since $O_Z \cap (Z^b + \mathbb{I}^{\perp})$ is convex, we further obtain that the whole segment $Z^b + \alpha(Z^\sharp - Z^b)$, $0 \le \alpha \le 1$, lies in $O_Z$ since $Z^b \in O_Z$ and $Z^\sharp - Z^b \in \mathbb{I}^{\perp}$. We may thus write that

$$ W^{\sharp} - W^{b} = \int_{0}^{1} \delta_{Z} W(Z^{b} + \alpha(Z^{\sharp} - Z^{b})) (Z^{\sharp} - Z^{b}) d\alpha = 0, $$

and multiplying scalarly by $Z^\sharp - Z^b$, we obtain that

$$ \int_{0}^{1} \langle \delta_{Z} W(Z^{b} + \alpha(Z^{\sharp} - Z^{b})) (Z^{\sharp} - Z^{b}) d\alpha, Z^{\sharp} - Z^{b} \rangle = 0, $$

and this yields

$$ \sum_{i,j \in \Xi} \int_{0}^{1} \frac{\partial_{\rho_i} g_j}{T} (Z^b + \alpha(Z^\sharp - Z^b)) d\alpha (\rho_i^\sharp - \rho_i^b)(\rho_j^\sharp - \rho_j^b) = 0. \quad (3.22) $$
Letting $\varrho = (\rho_1, \dots, \rho_n)^t$, and using that $T^b = T^\sharp$ and $\rho^b = \rho^\sharp$, we then deduce that

$$ \int_0^1 \langle \Lambda(Z^b + \alpha(Z^\sharp - Z^b)) d\alpha \ (\varrho^\sharp - \varrho^b), \varrho^\sharp - \varrho^b \rangle = 0, $$

and this implies that $\varrho^\sharp - \varrho^b = 0$ since $\varrho^\sharp - \varrho^b \in \mathbb{I}^\perp$ and $\Lambda$ is positive definite over $\mathbb{I}^\perp$, and the proof is complete. $\square$

### 3.4. Ideal gas mixtures

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

$$ \mathcal{P}^{\text{id}} = RT \sum_{k \in \mathcal{Z}} \frac{\rho_k}{m_k}, \qquad (3.23) $$

and the energy per unit volume reads

$$ \mathcal{E}^{\text{id}} = \sum_{k \in \mathcal{Z}} \rho_k e_k^{\text{id}}, \quad e_k^{\text{id}} = e_k^{\text{st}} + \int_{T^{\text{st}}}^{T} c_{vk}^{\text{id}}(\theta) d\theta, \qquad (3.24) $$

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

$$ S^{\text{id}} = \sum_{k \in \mathcal{Z}} \rho_k s_k^{\text{id}}, \quad s_k^{\text{id}} = s_k^{\text{st}} + \int_{T^{\text{st}}}^{T} \frac{c_{vk}^{\text{id}}(\theta)}{\theta} d\theta - \frac{RT}{m_k} \log \frac{\rho_k}{m_k \gamma^{\text{st}}}, \qquad (3.25) $$

where $s_k^{\text{st}}$ is the formation entropy of the $k$th species per unit mass at the standard temperature $T^{\text{st}}$ and standard pressure $p^{\text{st}}$, and $\gamma^{\text{st}} = p^{\text{st}}/RT^{\text{st}}$ is the standard concentration. The enthalpy $H^{\text{id}} = \mathcal{E}^{\text{id}} + \mathcal{P}^{\text{id}}$, Gibbs function $G^{\text{id}} = \mathcal{E}^{\text{id}} + \mathcal{P}^{\text{id}} - TS^{\text{id}}$, and free energy $F^{\text{id}} = \mathcal{E}^{\text{id}} - TS^{\text{id}}$ are then easily evaluated. The assumptions required for ideal gas mixtures in this volumetric framework are yet Property (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–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 $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^{\text{id}} = O(\rho^2)$ at low densities $\rho$, where $p^{\text{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, y_1, \dots, y_n)^t$ defined on a connected open set $O_\zeta \subset (0, \infty)^{n+2}$ and define the pressure corrector $\phi(\zeta) = p(\zeta) - p^{\text{id}}(\zeta)$ so that

$$p = p^{\text{id}} + \phi. \qquad (4.1)$$

The ideal gas pressure $p^{\text{id}}$ has been introduced in Section 2 and discussed in Sections 2.4 and 3.4. The open set $O_\zeta$ is said to be increasing with volume if for any $(T, v, y_1, \dots, y_n)^t \in O_\zeta$ whenever $v' > v$, we have $(T, v', y_1, \dots, y_n)^t \in O_\zeta$. We introduce the following properties concerning the pressure law (4.1).

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

(P$_1$) The open set $O_\zeta \subset (0, \infty)^{n+2}$ is increasing with volume.

(P$_2$) There exists $\underline{v} > 0$ such that for any $(T, y_1, \dots, y_n)^t \in (0, \infty)^{1+n}$ with $\sum_{i \in \mathcal{I}} y_i = 1$, and any $v \ge \underline{v}$, we have $\zeta = (T, v, y_1, \dots, y_n)^t \in O_\zeta$. The corrector $\phi$ and its partial derivatives $\partial_\zeta^\beta \phi$ of order $|\beta| \le 1 + \gamma$ have a smooth extension to the set $T > 0$, $v \ge \underline{v}$, $y_i \ge 0$, $\langle y, \mathbb{I} \rangle = 1$. There exist positive continuous functions of temperature $c(\gamma, T)$ depending on $\gamma$, such that for any $T > 0$, $v \ge \underline{v}$, and $y_i \ge 0$, $i \in \mathcal{I}$, $\langle y, \mathbb{I} \rangle = 1$, we have

$$|\partial_{\zeta}^{\beta} \phi| \le \frac{c(\gamma, T)}{v^{\beta v+2}}, \quad 0 \le |\beta| \le 1 + \gamma. \qquad (4.2)$$

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

(P$_3$) For any $\zeta \in O_\zeta$, $\partial_T^2 \phi(\zeta) \le 0$, $O_\zeta$ is increasing with temperature.

The 0-homogeneity property of pressure (P$_0$) is natural since we seek to obtain a thermodynamics such that (T$_0$) holds. The structure of the open set $O_\zeta$ with (P$_1$) is also required in order to allow integrations with respect to $v$ and to write the compatibility with ideal gases. In (P$_2$), we have used the traditional notation $\partial^\beta = \partial_T^{\beta_T} \partial_v^{\beta_v} \partial_{y_1}^{\beta_1} \dots \partial_{y_n}^{\beta_n}$ where $\beta = (\beta_T, \beta_v, \beta_1, \dots, \beta_n) \in \mathbb{N}^{2+n}$ is a multiindex. The condition (P$_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 \to 0$ since the linear term is associated with the ideal gas contribution $p^{\text{id}}$. Property (P$_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 (P$_0$-P$_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 $(P_0-P_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

$$e = e^{\text{id}} - \int_{v}^{\infty} T^2 \partial_{T} \left( \frac{\phi}{T} \right) (T, v', y_1, \dots, y_n) dv', \quad (4.3)$$

$$s = s^{\text{id}} - \int_{v}^{\infty} \partial_{T} \phi(T, v', y_1, \dots, y_n) dv', \quad (4.4)$$

$$g_k = g_k^{\text{id}} + \int_{v}^{\infty} \partial_{y_k} \phi(T, v', y_1, \dots, y_n) dv', \quad k \in \mathfrak{S}. \quad (4.5)$$

The triplet function $e$, $p$, and $s$ defined on $O_{\zeta}$ is then a thermodynamics in the sense of Definition 2.1, and Properties $(T_0-T_2)$ hold. Moreover, if $(P_3)$ holds, then $(T_3)$ also holds.

*Proof.* We first obtain from the compatibility relations (2.8) that $\partial_v e = T^2 \partial_T(\frac{p}{T})$ so that $\partial_v e = T^2 \partial_T(\frac{\phi}{T})$ since $\partial_T(p^{\text{id}}/T) = 0$ from (2.17). Integrating with respect to $v'$ over $[v, v'']$, keeping in mind that $O_{\zeta}$ in increasing with respect to volume, we obtain

$$e(T, v, y_1, \dots, y_n) - e(T, v'', y_1, \dots, y_n) = - \int_{v}^{v''} T^2 \partial_T \left( \frac{\phi}{T} \right) (T, v', y_1, \dots, y_n) dv'$$

Using the integrability properties (4.2) and $\lim_{v'' \to \infty} e(T, v'', y_1, \dots, y_n) = e^{\text{id}}(T, y_1, \dots, y_n)$ from $(T_2)$, we may pass to the limit $v'' \to \infty$ and we obtain (4.3).

Similarly, from $\partial_v(Ts - e) = p$ and $\partial_v(Ts^{\text{id}} - e^{\text{id}}) = p^{\text{id}}$, we deduce $\partial_v(Ts - e) - \partial_v(Ts^{\text{id}} - e^{\text{id}}) = \phi$. Integrating with respect to $v'$ over $[v, v'']$ yields

$$(Ts - Ts^{\text{id}} - e + e^{\text{id}})(T, v, y_1, \dots, y_n) - (Ts - Ts^{\text{id}} - e + e^{\text{id}})(T, v'', y_1, \dots, y_n) = - \int_{v}^{v''} \phi(T, v', y_1, \dots, y_n) dv'$$

Passing to the limit $v'' \to \infty$, using (4.2) and $(T_2)$, and noting that $(Ts - Ts^{\text{id}} - e + e^{\text{id}})(T, v'', y_1, \dots, y_n)$ goes to zero as $v'' \to \infty$ yields (4.4).

We now establish that $(T_0-T_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 $y_k$, $k \in \mathfrak{S}$, and deduce that $e$ and $s$ are $C^{\gamma}$ functions over $O_{\zeta}$. Deriving with respect to the mass fraction $y_k$, we may calculate $g_k = \partial_{y_k} e - T \partial_{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 $O_{\zeta}$.

We also note that if $\zeta \in O_{\zeta}$ and $\lambda > 0$, then $\zeta_{\lambda} \in O_{\zeta}$ from $(P_0)$. Taking into account the homogeneity properties of $e^{\text{id}}$, $s^{\text{id}}$, and $g_k^{\text{id}}$, $k \in \mathfrak{S}$, and from $(P_0)$, we easily deduce that $e$ and $s$ are 1-homogeneous over $O_{\zeta}$ and that $g_k, k \in \mathfrak{S}$ are 0-homogeneous over $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 dv'$ and $\int_{v}^{\infty} (T \partial_T \phi - \phi) dv'$ are then 1-homogeneous. Similarly, $\partial_{y_k} \phi$ is (-1)-homogeneous so that $\int_{v}^{\infty} \partial_{y_k} \phi dv'$ is 0-homogeneous and property $(T_0)$ holds.

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

$$\partial_T e = \partial_T e^{\text{id}} - T \int_{v}^{\infty} \partial_T^2 \phi \, dv', \quad \partial_T s = \partial_T s^{\text{id}} - \int_{v}^{\infty} \partial_T^2 \phi \, dv'$$
and since $\partial_T e^{id} = T \partial_T s^{id}$, we obtain that $\partial_T e = T \partial_T s$. Similarly, from (4.3), (4.4), and the estimates (4.2), we have

$$ \partial_v e = \partial_v e^{id} + T \partial_T \phi - \phi, \quad \partial_v s = \partial_v s^{id} + \partial_T \phi. $$

Using now $\partial_v e^{id} = T \partial_v s^{id} - p^{id}$, we deduce that $\partial_v e = T \partial_v s^{id} - p^{id} + T \partial_T \phi - \phi$, and thus $\partial_v e = T \partial_v s - p$ which implies $(T_1)$. Keeping in mind that $g_k = \partial_{y_k} e - T \partial_{y_k} s$ by definition, and we have established properties $(T_1)$.

The property $(T_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, y_1, \dots, y_n)'$ belongs to $O_\zeta$ for $v$ large enough. The fact that $(P_3)$ implies $(T_3)$ is finally a consequence of (4.3) since

$$ \partial_T e = \partial_T e^{id} - T \int_v^\infty \partial_T^2 \phi \, dv', $$

$\partial_T e^{id} = \sum_{k \in \mathcal{Z}} y_k c_{v_k}^{id} > 0$ and $\partial_T^2 \phi \le 0$ from $(P_3)$, so that $\partial_T e > 0$ and $(T_3)$ is established. $\square$

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 (Ts - 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 (Ts - 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_{yy}^2 f/T$ is easily shown to be

$$ \Lambda_{kl} = \frac{R}{y_k m_k} \delta_{kl} + \frac{1}{T} \int_v^\infty \partial_{y_k y_l}^2 \phi \, dv', \quad k, l \in \mathcal{Z}. \qquad (4.6) $$

In the low density limit $v \to \infty$, the specific heat $\partial_T e = c_v$ tends to $c_v^{id}$, which is positive from (ID), and the matrix $\Lambda$ goes to $\lim_{v \to \infty} \Lambda = R \operatorname{diag}(\frac{1}{m_1 y_1}, \dots, \frac{1}{m_n y_n})$, 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 $\tilde{\Lambda}$, $\tilde{\Lambda}_r$, and $\tilde{\Lambda}_l$, that behave more regularly for zero mass fractions. For any diagonalizable matrix $A$, we denote by $d^+(A)$, $d^0(A)$, and $d^-(A)$ the number of positive, zero and negative eigenvalues of $A$, respectively.

**Proposition 4.2.** Assume that Properties $(P_0-P_2)$ hold. For any $\zeta \in O_\zeta$, let $\Lambda$ be defined as in (2.4) and $\Pi$ be the diagonal matrix

$$ \Pi = \operatorname{diag}(\sqrt{m_1 y_1}, \dots, \sqrt{m_1 y_n}) / \sqrt{R}. $$

The matrices $\tilde{\Lambda}$, $\tilde{\Lambda}_r$, and $\tilde{\Lambda}_l$, defined by

$$ \tilde{\Lambda} = \Pi \wedge \Pi, \quad \tilde{\Lambda}_r = \Lambda \Pi^2, \quad \tilde{\Lambda}_l = \Pi^2 \Lambda, $$
are then $C^{\gamma-1}$ functions of $\zeta \in O_\zeta$. Furthermore, if we assume that $v > v_*$, then $\bar{\Lambda}$ admits a continuous extension to nonnegative mass fractions whereas $\hat{\Lambda}_r$ and $\hat{\Lambda}_l$ admit $C^{\gamma-1}$ extensions to nonnegative mass fractions. Finally, for any $\zeta \in O_\zeta$, we have $d^+(\Lambda) = d^+(\bar{\Lambda})$, $d^0(\Lambda) = d^0(\bar{\Lambda})$, $d^-(\Lambda) = d^-(\bar{\Lambda})$, and the matrices $\bar{\Lambda}$, $\hat{\Lambda}^r$, and $\hat{\Lambda}_l$ have the same spectrum.

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

We deduce from Proposition 4.2 that $d^+(\Lambda) = d^+(\bar{\Lambda}) = d^+(\hat{\Lambda}_r) = d^+(\hat{\Lambda}_l)$, $d^0(\Lambda) = d^0(\bar{\Lambda}) = d^0(\hat{\Lambda}_r) = d^0(\hat{\Lambda}_l)$, and $d^-(\Lambda) = d^-(\bar{\Lambda}) = d^-(\hat{\Lambda}_r) = d^-(\hat{\Lambda}_l)$, in such a way that any of the matrices $\Lambda$, $\bar{\Lambda}$, $\hat{\Lambda}_r$, or $\hat{\Lambda}_l$, may be used over $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 $y$ with $1 \le n^+ < n$. The general case is easily reduced to this situation by introducing permutation matrices [8]. We then have $y_k > 0$ for $1 \le k \le n^+$ and $y_k = 0$ for $n^+ + 1 \le k \le n$. Let's denote by $n^0 = n - n^+$ the number of zero mass fractions, $y^+$ the vector of positive mass fractions, and $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 $x \in \mathbb{R}^n$ may be written in the form $x = (x^+, x^0)$, and we have, for instance, $\mathbb{I}^+ \in \mathbb{R}^{n^+}$ and $\mathbb{I}^+ = (1, \dots, 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 $(Ax)^+ = A^{++}x^+ + A^{+0}x^0$ and $(Ax)^0 = A^{0+}x^+ + A^{00}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 $\bar{\Lambda}$, $\hat{\Lambda}_r$, and $\hat{\Lambda}_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 (P₀-P₁) hold. Assume that $T > 0$, $v > v_*$, and that the mass fraction vector is in the form $y = (y^+, y^0)$, where $y^+ \in \mathbb{R}^{n^+}$, $y^0 \in \mathbb{R}^{n^0}$, $y^+ = (y_1, \dots, y_{n^+})^t$, $y^0 = (0, \dots, 0)^t$, $y_i > 0$, $1 \le i \le n^+$, and $\langle y, \mathbb{I} \rangle = \langle y^+, \mathbb{I}^+ \rangle = 1$. Then, we have the block decompositions

$$ \bar{\Lambda} = \begin{pmatrix} \bar{\Lambda}^{++} & 0 \\ 0 & \mathbb{I} \end{pmatrix}, \quad \hat{\Lambda}_r = \begin{pmatrix} \hat{\Lambda}_r^{++} & 0 \\ \hat{\Lambda}_r^{0+} & \mathbb{I} \end{pmatrix}, \quad \hat{\Lambda}_l = \begin{pmatrix} \hat{\Lambda}_l^{++} & \hat{\Lambda}_l^{+0} \\ 0 & \mathbb{I} \end{pmatrix}, $$

where $\mathbb{I}$ denotes the identity matrix of size $n^0$ and the symbol $0$ denotes various rectangular matrices with zero entries so that $\bar{\Lambda}^{+0} = 0$, $\hat{\Lambda}^{+0} = 0$, $\hat{\Lambda}_r^{+0} = 0$, and $\hat{\Lambda}_l^{+0} = 0$. Moreover, denoting by

$$ \Pi^{++} = \mathrm{diag}(\sqrt{m_1 y_1}, \dots, \sqrt{m_{n^+} y_{n^+}}) / \sqrt{R}, $$

the upper left block of the matrix $\Pi$ we have the relations $\bar{\Lambda}^{++} = \Pi^{++}\Lambda^{++}\Pi^{++}$, $\hat{\Lambda}_r^{++} = \Lambda^{++}(\Pi^{++})^2$, $\hat{\Lambda}_l^{0+} = \Lambda^{0+}(\Pi^{++})^2$, $\hat{\Lambda}_l^{++} = (\Pi^{++})^2\Lambda^{++}$, $\hat{\Lambda}_l^{+0} = (\Pi^{++})^2\Lambda^{+0}$. In particular, we have $d^+(\Lambda^{++}) = d^+(\bar{\Lambda}^{++})$, $d^0(\Lambda^{++}) = d^0(\bar{\Lambda}^{++})$, $d^-(\Lambda^{++}) = d^-(\bar{\Lambda}^{++})$, and the matrices $\bar{\Lambda}^{++}$, $\hat{\Lambda}_r^{++}$, and $\hat{\Lambda}_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 $\hat{\Lambda}_r$ or $\hat{\Lambda}_l$ which, behave smoothly for
nonnegative mass fractions. The whole stability domain may then be obtained over the set $y_i \ge 0$, $i \in \mathfrak{S}$, $\langle y, \mathbf{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 $\mathbf{e}^1, \dots, \mathbf{e}^n$ the canonical base vectors of $\mathbb{R}^n$ and by $\zeta^i$ the state $\zeta = (T, v, \mathbf{e}^i)^t$ for $y = \mathbf{e}^i$, the following proposition discusses the thermodynamic stability in the neighborhood of $\mathbf{e}^i$.

**Proposition 4.4.** Assume that for $\zeta = \zeta^i$ we have $\partial_T e(\zeta^i) = c_{vi} > 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 $\tilde{\Lambda}_r$ and $\tilde{\Lambda}_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]

$$p = \sum_{i \in \mathfrak{S}} \frac{y_i RT}{m_i v - b} - \frac{a}{v(v + b)}, \qquad (5.1)$$

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, y_1, \dots, y_n)$ and $b(y_1, \dots, y_n)$ may be evaluated by using the Van der Waals mixing rules:

$$a = \sum_{i,j \in \mathfrak{S}} y_i y_j \alpha_i \alpha_j \qquad b = \sum_{i \in \mathfrak{S}} y_i b_i, \qquad (5.2)$$

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 \ge 2$.

(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 \ge 0$, $\lim_{+\infty} \alpha_i = 0$, $\partial_T \alpha_i \le 0$, and $\partial_{TT}^2 \alpha_i \ge 0$ 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 (SRK). The following proposition is easily deduced from (SRK).
**Proposition 5.1.** Assume that (SRK) holds. Then, $a = a(T, y) \in C^0([0, \infty)^{1+n}) \cap C^{1+\gamma}((0, \infty)^{1+n})$, $a(0, y) > 0$, $a \ge 0$, $\lim_{T \to +\infty} a(T, y) = 0$, $\partial_T a \le 0$, and $\partial_{TT}^2 a \ge 0$ over $(0, \infty) \times (0, \infty)^n$.

*Proof.* This is a consequence of (5.2) and (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

$$ \phi = \left( \sum_{i \in \mathcal{I}} \frac{RTy_i}{m_i} \right) \frac{b}{v(v-b)} - \frac{a}{v(v+b)}. \qquad (5.3) $$

Under assumptions (SRK), the corrector $\phi$ is such that $(P_0-P_2)$ hold on $O_\zeta$ given by

$$ O_{\zeta} = \{ \zeta \in (0, \infty)^{2+n}; v > b \}, \qquad (5.4) $$

with notably $\underline{v} = 2 \max_{i \in \mathcal{I}} b_i$.

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

From Properties $(P_0-P_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

$$ e = \sum_{i \in \mathcal{I}} y_i e_i^{id} + (T \partial_T a - a) \frac{1}{b} \log(1 + \frac{b}{v}), \qquad (5.5) $$

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

$$ s = \sum_{i \in \mathcal{I}} y_i s_i^{id*} - \sum_{i \in \mathcal{I}} \frac{y_i R}{m_i} \log\left(\frac{y_i R T}{m_i (v-b) p^{st}}\right) + \frac{\partial_T a}{b} \log\left(1 + \frac{b}{v}\right), \qquad (5.6) $$

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

$$ \Lambda_{ij} = \frac{R \delta_{ij}}{m_i y_i} + \frac{R}{v-b} \left( \frac{b_i}{m_j} + \frac{b_j}{m_i} \right) + \sum_{k \in \mathcal{I}} \frac{y_k}{m_k} \frac{R}{(v-b)^2} b_i b_j - \frac{2 \alpha_i \alpha_j}{T b} \log\left(1 + \frac{b}{v}\right) \\ + \frac{2}{T} \sum_{k \in \mathcal{I}} y_k (\alpha_i \alpha_k b_j + \alpha_j \alpha_k b_i) \left( \frac{1}{b^2} \log\left(1 + \frac{b}{v}\right) - \frac{1}{b(v+b)} \right) \\ + \frac{1}{T} a b_i b_j \left( - \frac{2}{b^3} \log\left(1 + \frac{b}{v}\right) + \frac{2}{b^2(v+b)} + \frac{1}{b(v+b)^2} \right), \quad i, j \in \mathcal{I}. \qquad (5.7) $$

We now establish that property $(P_3)$ holds under (SRK) so that thermal stability automatically holds.
**Proposition 5.3.** *Assuming (SRK), the corrector $\phi = p - p^{id}$ is such that $\partial_{TT}^2 \phi \le 0$ and $O_\zeta$ is increasing with respect to temperature T so that (P$_3$) holds.*

*Proof.* By deriving (5.1) with respect to temperature, we obtain $\partial_{TT}^2 \phi = -\partial_{TT}^2 a/v(v + b)$. On the other hand, using (5.2) we have $\partial_{TT}^2 a = \sum_{i,j \in \mathcal{Z}} y_i y_j (\alpha_i \partial_{TT}^2 \alpha_j + 2\partial_T \alpha_i \partial_T \alpha_j + \alpha_i \partial_{TT}^2 \alpha_j)$, and since $\alpha_i \partial_{TT}^2 \alpha_j \ge 0$ and $\partial_T \alpha_i \partial_T \alpha_j \ge 0$, we obtain that $\partial_{TT}^2 a \ge 0$ and $\partial_{TT}^2 \phi \le 0$ from (SRK). Moreover, the open set $O_\zeta$ defined by (5.4) is independent of T and (P$_3$) is established. $\square$

The attractive parameters $\alpha_i$, $i \in \mathcal{Z}$ 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*

$$ \frac{\alpha_i}{\alpha_{c,i}} = 1 + A(s_i(1 - \sqrt{T}/\sqrt{T_{c,i}})), \qquad (5.8) $$

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

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

$$ \frac{\partial_T \alpha_i}{\alpha_{c,i}} = - \frac{s_i}{2(TT_{c,i})^{1/2}} A'(s_i(1 - \sqrt{T}/\sqrt{T_{c,i}})), $$

$$ \frac{\partial_{TT}^2 \alpha_i}{\alpha_{c,i}} = \frac{s_i}{4T^{3/2}T_{c,i}^{1/2}} A'(s_i(1 - \sqrt{T}/\sqrt{T_{c,i}})) + \frac{s_i^2}{4TT_{c,i}} A''(s_i(1 - \sqrt{T}/\sqrt{T_{c,i}})). $$

From the assumption on A, we deduce that $\partial_T \alpha_i(T) \le 0$ and $\partial_{TT}^2 \alpha_i(T) \ge 0$ so that $\alpha_i$ is nonincreasing and convex. Moreover, we have $\alpha_i(0) = \alpha_{c,i}(1 + A(s_i))$ and $A(s_i) \ge 0$ 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. $\square$

The parameter $\alpha_{c,i}$ is the value of $\alpha_i$ at the critical temperature $T_{c,i}$ and $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 \ge 0$, and $A(x) = x/(1 + x^4)^{1/4}$ when $x \le 0$. It is then checked that $A \in C^4(\mathbb{R})$, $A' \ge 0$, $A'' \ge 0$, $A(0) = 0$, $-1 \le A$, and $\lim_{-\infty} A = -1$, and the assumptions of Proposition 5.4 hold with $\gamma = 3$. The Soave coefficients correspond to $A(x) = |1 + x| - 1$ and neither satisfy property (SRK) nor the assumptions of Proposition 5.4. It is possible to truncate those coefficients by using $A(x) = |1 + x|^{\gamma} - 1$, where $\gamma = \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 \mathcal{Z}} y_i y_j \alpha_i \alpha_j$ at the crossing temperatures $T_{c,i}(1 + 1/s_i)^2$, $i \in \mathcal{Z}$. 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 \mathcal{Z}} y_i/m_i$ so that the equation of state is written in the form

$$ p = \frac{rT}{v-b} - \frac{a}{v(v+b)}. \qquad (5.9) $$
**Proposition 5.5.** For any $y \in (0, \infty)^n$ and $v > b$, there exists a unique temperature $T^*(v, y) > 0$ such that $\partial_v p(T^*, v, y) = 0$. This temperature satisfies the nonlinear equation

$$ T^* = \frac{a^*(v-b)^2(2v+b)}{r v^2(v+b)^2}, \qquad (5.10) $$

where $a^* = a(T^*, y)$, $T^*$ is a $C^\gamma$ function of $v, y_1, \dots, y_n$, and we have the relations

$$ \partial_{T_v}^2 p \partial_v T^* + \partial_{vv}^2 p = 0, \quad \partial_{T_v}^2 p \partial_{y_k} T^* + \partial_{v y_k}^2 p = 0, \quad k \in \mathcal{Z}. \qquad (5.11) $$

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

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

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

**Corollary 5.6.** We have the equivalence $\partial_v p(T, v, y) > 0 \iff T < T^*(v, y)$, and since the derivative $\partial_T p$ at fixed $v, y_1, \dots, y_n$ is always positive, defining $p^*(v, y) = p(T^*(v, y), v, y)$, we also have the equivalence $\partial_v p(T, v, y) > 0 \iff p(T, v, y) < p^*(v, y)$.

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

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

$$ T^{**}(y) = T^{*}(\theta_{c}b, y), \quad p^{**}(y) = p^{*}(\theta_{c}b, y) = p(T^{*}(\theta_{c}b, y), \theta_{c}b, y), \qquad (5.12) $$

keeping in mind that $b = \sum_{i \in \mathcal{Z}} y_i b_i$. Then, $T^{**}(y)$ and $p^{**}(y)$ correspond to the critical temperature and pressure of the frozen mixture with mass fractions $y$. Moreover, we have $p^*(v, y) \le p^{**}(y)$ so that $\partial_v p > 0$ implies $p(T, v, y) < p^{**}(y)$. Finally, since the set $\bar{\Sigma} = \{y \in [0, \infty)^n; \langle y, \mathbb{I} \rangle = 1\}$ is compact and $p(T^*(\theta_c v, y), \theta_c v, y)$ is a smooth function of $y$, there exists a maximum pressure $\max_{y \in \bar{\Sigma}} p^{**}(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–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}^{\text{ex}}$ is of Van der Waals type:

$$ \mathcal{F}^{\text{ex}} = \mathcal{F} + \frac{1}{2} \kappa |\nabla \rho|^2, \qquad (6.1) $$

where $\mathcal{F}$ denotes the classical or bulk free energy. The superscript 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, \dots, \rho_n$, whereas the gradient term $\frac{1}{2}\kappa|\nabla\rho|^2$ in $\mathcal{F}^{\text{ex}}$ represents an interfacial excess free energy. This gradient term may be interpreted from attractive long range molecular interactions and the capillary coefficient $\kappa$ related to the species interaction potentials [25, 44].

From the definition (6.1) and the classical relation $d\mathcal{F} = -SdT + \sum_{i \in \mathcal{I}} g_i d\rho_i$, it is obtained that $d\mathcal{F}^{\text{ex}} = -(S - \frac{1}{2}\partial_T \kappa |\nabla \rho|^2)dT + \sum_{i \in \mathcal{I}} g_i d\rho_i + \kappa \nabla \rho \cdot d\nabla \rho$ where $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}^{\text{ex}} = -S^{\text{ex}}$ and $\partial_{\rho_i} \mathcal{F}^{\text{ex}} = g_i^{\text{ex}}$, $i \in \mathcal{I}$, where $S^{\text{ex}}$ denotes the extended entropy per unit volume and $g_i^{\text{ex}}$ the Gibbs function per unit mass of the $i$th species, it is then obtained that $S^{\text{ex}} = S - \frac{1}{2}\partial_T \kappa |\nabla \rho|^2$ and $g_i^{\text{ex}} = g_i$.

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

$$ \mathcal{E}^{\text{ex}} = \mathcal{E} + \frac{1}{2}(\kappa - T\partial_T \kappa)|\nabla \rho|^2, \qquad p^{\text{ex}} = p - \frac{1}{2}\kappa|\nabla \rho|^2. \qquad (6.2) $$

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

$$ T d\mathcal{S}^{\text{ex}} = d\mathcal{E}^{\text{ex}} - \sum_{i \in \mathcal{I}} g_i d\rho_i - \kappa \nabla \rho \cdot d\nabla \rho, \qquad (6.3) $$

where $d$ is the differentiation operator.

## 6.2. Cohesive fluids

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

$$ \partial_t \rho_i + \nabla \cdot (\rho_i \mathbf{u}) + \nabla \cdot \mathcal{J}_i = m_i \omega_i, \qquad i \in \mathcal{I} = \{1, \dots, n\}, \qquad (6.4) $$

$$ \partial_t (\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) + \nabla \cdot \mathcal{P} = 0, \qquad (6.5) $$

$$ \partial_t (\mathcal{E}^{\text{ex}} + \frac{1}{2}\rho|\mathbf{u}|^2) + \nabla \cdot (\mathbf{u}(\mathcal{E}^{\text{ex}} + \frac{1}{2}\rho|\mathbf{u}|^2)) + \nabla \cdot (\mathbf{Q} + \mathcal{P} \cdot \mathbf{u}) = 0, \qquad (6.6) $$

where $\partial_t$ denotes the time derivative operator, $\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, $\mathbf{u}$ the fluid velocity, $\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, $\mathcal{P}$ the total pressure tensor, $\mathcal{E}^{\text{ex}}$ the internal energy per unit volume, and $\mathbf{Q}$ the total heat flux.
The diffusive fluxes $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}} 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

$$ \partial_t \rho + \nabla \cdot (\rho u) = 0. \qquad (6.7) $$

A conservation equation for the internal energy $\mathcal{E}^{\text{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

$$ \mathcal{P} = p^{\text{ex}} \mathbf{I} + \kappa \nabla \rho \otimes \nabla \rho - \rho \nabla \cdot (\kappa \nabla \rho) \mathbf{I} + \mathcal{P}^{\text{d}}, \qquad (6.8) $$

$$ \mathcal{Q} = \kappa \rho \nabla \cdot \mathbf{u} \nabla \rho + \mathbf{q}, \qquad (6.9) $$

with the dissipative fluxes given by

$$ \mathcal{P}^{\text{d}} = - \nu \nabla \cdot \mathbf{u} \mathbf{I} - \eta (\nabla \mathbf{u} + \nabla \mathbf{u}^t - \frac{2}{3} \nabla \cdot \mathbf{u} \mathbf{I}), \qquad (6.10) $$

$$ \mathcal{J}_i = - \sum_{j \in \mathfrak{S}} L_{ij} \nabla \left( \frac{g_j}{T} \right) - L_{ie} \nabla \left( \frac{-1}{T} \right), \quad i \in \mathfrak{S}, \qquad (6.11) $$

$$ \mathbf{q} = - \sum_{i \in \mathfrak{S}} L_{ei} \nabla \left( \frac{g_i}{T} \right) - L_{ee} \nabla \left( \frac{-1}{T} \right), \qquad (6.12) $$

where $\nu$ denotes the volume viscosity, $\eta$ the shear viscosity, and $L_{ij}$, $i, j \in \mathfrak{S} \cup \{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 \{e\}}$ is symmetric positive semi-definite with null space spanned by the vector $(1, \dots, 1, 0)^t$ as for ordinary fluids [44, 60].

The dissipative fluxes $J_i$, $i \in \mathfrak{S}$, and $\mathbf{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^{\text{ex}} = g_j$, and of the usual thermodynamic thermal variable $-1/T$. In the pressure tensor $\mathcal{P}$, the new contributions from the Korteweg tensor [22] are $-\rho \nabla \cdot (\kappa \nabla \rho) \mathbf{I}$ and $\kappa \nabla \rho \otimes \nabla \rho$, and the pressure $p^{\text{ex}}$ also differs from $p$. The extra capillary contribution in the total heat flux is the Dunn and Serrin heat flux $\kappa \rho \nabla \cdot \mathbf{u} \nabla \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

$$ \sum_{i \in \mathfrak{S}} v_{ij}^{\text{f}} \mathfrak{M}_i \rightleftharpoons \sum_{i \in \mathfrak{S}} v_{ij}^{\text{b}} \mathfrak{M}_i, \quad j \in \mathbb{R}, \qquad (6.13) $$

where $v_{ij}^{\text{f}}$ and $v_{ij}^{\text{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, $\mathbb{R} = \{1, \dots, n^r\}$ the indexing set of chemical reactions, and $n^r$ the number of chemical reactions. The production rates are given by

$$ \omega_i = \sum_{j \in \mathbb{R}} v_{ij} \tau_j, $$
where $v_{ij} = v_{ij}^b - v_{ij}^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

$$v_j = \begin{pmatrix} v_{1j} \\ \vdots \\ v_{nj} \end{pmatrix}, \quad v_j^f = \begin{pmatrix} v_{1j}^f \\ \vdots \\ v_{nj}^f \end{pmatrix}, \quad v_j^b = \begin{pmatrix} v_{1j}^b \\ \vdots \\ v_{nj}^b \end{pmatrix}, \quad \mu = \begin{pmatrix} \mu_1 \\ \vdots \\ \mu_n \end{pmatrix},$$

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–91] then read

$$\tau_j = \mathcal{K}_j(\exp\langle v_j^f, \mu \rangle - \exp\langle v_j^b, \mu \rangle), \quad j \in \mathbb{R}, \qquad (6.14)$$

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

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

$$\partial_t S^{\text{ex}} + \nabla \cdot (u S^{\text{ex}}) + \nabla \cdot \left(\frac{Q}{T} - \frac{\kappa \rho \nabla \rho \nabla \cdot u}{T} - \sum_{i \in \Xi} \frac{g_i}{T} \mathcal{J}_i\right) = -\frac{1}{T} (\mathcal{P} - p \mathcal{I} - (\kappa \nabla \rho \otimes \nabla \rho - \rho \nabla \cdot (\kappa \nabla \rho) \mathcal{I})) : \nabla u \\ - (\mathcal{Q} - \kappa \rho \nabla \rho \nabla \cdot u) \cdot \nabla \left(\frac{-1}{T}\right) - \sum_{i \in \Xi} \mathcal{J}_i \cdot \nabla \left(\frac{g_i}{T}\right) - \sum_{i \in \Xi} \frac{g_i m_i \omega_i}{T}. \qquad (6.15)$$

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

$$\partial_t S^{\text{ex}} + \nabla \cdot (u S^{\text{ex}}) + \nabla \cdot \left(\frac{q}{T} - \sum_{i \in \Xi} \frac{g_i}{T} \mathcal{J}_i\right) = \frac{v}{T} (\nabla \cdot u)^2 + \frac{\eta}{2T} |S|^2 \\ + \sum_{i,j \in \Xi \cup \{e\}} L_{ij} \nabla \left(\frac{g_i}{T}\right) \cdot \nabla \left(\frac{g_j}{T}\right) + \sum_{j \in \mathcal{R}} R \mathcal{K}_j \left(\langle v_j^f, \mu \rangle - \langle v_j^b, \mu \rangle\right) \left(\exp\langle v_j^f, \mu \rangle - \exp\langle v_j^b, \mu \rangle\right),$$

where we have let $g_e = -1$ for convenience, $S$ denotes the deviatoric part of the strain rate tensor $S = \nabla u + \nabla u' - \frac{2}{3} \nabla \cdot u \mathcal{I}$, and $|A|^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 $w$ representing locally the gradient of density [82–87]:

$$w = \nabla \rho. \qquad (6.16)$$

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

$$\partial_t w + \sum_{i \in \mathcal{D}} \partial_i (w u_i + \rho \nabla u_i) = 0, \qquad (6.17)$$
where $\mathcal{D} = \{1, \dots, 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 $\mathbf{u} = (u_1, \dots, u_d)^t$ and $\nabla = (\partial_1, \dots, \partial_d)^t$ in $d$ dimensions.

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

$$ \mathbf{u} = (\rho_1, \dots, \rho_n, \mathbf{w}, \rho\mathbf{u}, \mathcal{E} + \frac{1}{2}\rho|\mathbf{u}|^2)^t, \quad \mathbf{z} = (\rho_1, \dots, \rho_n, \mathbf{w}, \mathbf{u}, T)^t, \qquad (6.18) $$

and are vectors of $\mathbb{R}^n$ with $n = n + 2d + 1$. The vectors of $\mathbb{R}^n$ are decomposed as $\mathbb{R}^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

$$ \partial_t \mathbf{u} + \sum_{i \in \mathcal{D}} \partial_i (\mathbf{F}_i + \mathbf{F}_i^c + \mathbf{F}_i^d) = \mathbf{\Omega}, \qquad (6.19) $$

where $\mathbf{F}_i$ denotes the augmented convective flux in the $i$th direction, $\mathbf{F}_i^c$ the augmented cohesive flux in the $i$th direction, $\mathbf{F}_i^d$ the augmented dissipative flux in the $i$th direction, and $\mathbf{\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

$$ \mathbf{F}_i = (\rho_1 u_i, \dots, \rho_n u_i, \mathbf{w} u_i, \rho \mathbf{u} u_i + (p^{\text{ex}} + \kappa |\mathbf{w}|^2) \mathbf{e}_i, (\mathcal{E}^{\text{ex}} + \frac{1}{2} \rho |\mathbf{u}|^2 + p^{\text{ex}} + \kappa |\mathbf{w}|^2) u_i)^t, \quad i \in \mathcal{D} \qquad (6.20) $$

$$ \mathbf{F}_i^c = (0, \dots, 0, \rho \nabla u_i, -\rho \nabla(\kappa w_i), \rho \kappa \mathbf{w} \cdot \nabla u_i - \rho \mathbf{u} \cdot \nabla(\kappa w_i))^t, \qquad (6.21) $$

$$ \mathbf{F}_i^d = (\mathcal{J}_{li}, \dots, \mathcal{J}_{ni}, 0_{d,1}, \mathcal{P}_i^d, q_i + \sum_{j \in \mathcal{D}} \mathcal{P}_{ij}^d u_j)^t, \qquad (6.22) $$

$$ \mathbf{\Omega} = (m_1 \omega_1, \dots, m_n \omega_n, 0, 0, 0)^t, \qquad (6.23) $$

where $\mathcal{P}_{ij}^d, i, j \in \mathcal{D}$ denote the components of the viscous tensor $\mathcal{P}^d = (\mathcal{P}_{ij}^d)_{i,j \in \mathcal{D}}, \mathcal{P}_i^d, i \in \mathcal{D}$, the vectors $\mathcal{P}_i^d = (\mathcal{P}_{i1}^d, \dots, \mathcal{P}_{id}^d)^t, \mathbf{e}_i, i \in \mathcal{D}$, the base vectors in $\mathbb{R}^d, \mathbf{w} = (w_1, \dots, w_d)^t$, the spatial components of $\mathbf{w}$, $\mathbf{q} = (q_1, \dots, q_d)^t$, the spatial components of the dissipative heat flux $\mathbf{q}$, $\mathcal{J}_i = (\mathcal{J}_{li}, \dots, \mathcal{J}_{di})^t$, the spatial components of the mass flux $\mathcal{J}_i$ of the $i$th species and $\mathbf{\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 $\mathbf{F}_i^c$ is similar to that of single species fluids [87]. It is interesting to note that the relevant pressure is then $p^{\text{ex}} + \kappa |\mathbf{w}|^2 = p + \frac{1}{2}\kappa |\mathbf{w}|^2$.

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

$$ \mathbf{F}_i^d = - \sum_{j \in \mathcal{D}} \mathbf{B}_{ij}^d \partial_j \mathbf{u}, \quad \mathbf{F}_i^c = - \sum_{j \in \mathcal{D}} \mathbf{B}_{ij}^c \partial_j \mathbf{u}, \quad i \in \mathcal{D}, \qquad (6.24) $$

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

$$ \partial_t \mathbf{u} + \sum_{i \in \mathcal{D}} \mathbf{A}_i(\mathbf{u}) \partial_i \mathbf{u} - \sum_{i,j \in \mathcal{D}} \partial_i (\mathbf{B}_{ij}^d(\mathbf{u}) \partial_j \mathbf{u}) - \sum_{i,j \in \mathcal{D}} \partial_i (\mathbf{B}_{ij}^c(\mathbf{u}) \partial_j \mathbf{u}) = \mathbf{\Omega}(\mathbf{u}). \qquad (6.25) $$
The coefficient matrices $A_i = \partial_u F_i$, $B_{ij}^d$, and $B_{ij}^c$, for $i, j \in \mathcal{D}$, have at least regularity $C^{\gamma-1}$ over the open set $O_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–87]. The corresponding symmetrized equations have symmetric dissipative terms issued from $B_{ij}^d$, $i, j \in \mathcal{D}$, and antisymmetric cohesive terms issued from $B_{ij}^c$, $i, j \in \mathcal{D}$ [87].

## 6.5. Extended thermodynamics

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

$$ \begin{aligned} \mathcal{E}^{\text{ex}} &= \mathcal{E} + \frac{1}{2}(\kappa - T\partial_T\kappa)|w|^2, & p^{\text{ex}} &= p - \frac{1}{2}\kappa|w|^2, \\ \mathcal{S}^{\text{ex}} &= \mathcal{S} - \frac{1}{2}\partial_T\kappa|w|^2, & g_k^{\text{ex}} &= g_k, & k \in \mathfrak{S}. \end{aligned} $$

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

$$ \Lambda = (\Lambda_{ij})_{i,j \in \mathfrak{S}}, \quad \Lambda_{ij} = \frac{\partial_{\rho_j} g_i}{T} = \Lambda_{ij}^{\text{ex}}. \qquad (6.26) $$

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

(H0) $\mathcal{E}^{\text{ex}}$, $\mathcal{P}^{\text{ex}}$, $\mathcal{S}^{\text{ex}}$ are $C^\gamma$ functions of $z$ defined over $O_z \subset (0, \infty)^n \times \mathbb{R}^d \times \mathbb{R}^d \times (0, \infty)$. The open set $O_z$ is nonempty and connected and $\kappa = \kappa(T)$ is a $C^{\gamma+1}$ function of temperature $T$. If $(T, \rho_1, \dots, \rho_n)^t \in O_z$, then $(\rho_1, \dots, \rho_n, 0, 0, T)^t \in O_z$, and if $(\rho_1, \dots, \rho_n, w, u, T)^t \in O_z$, then $(T, \rho_1, \dots, \rho_n)^t \in O_z$.

(H1) Letting $g_k^{\text{ex}} = g = \partial_{\rho_k} \mathcal{E}^{\text{ex}} - T \partial_{\rho_k} \mathcal{S}^{\text{ex}}$, we have the extended volumetric Gibbs' relation $Td\mathcal{S}^{\text{ex}} = d\mathcal{E}^{\text{ex}} - \sum_{i \in \mathfrak{S}} g_i^{\text{ex}} d\rho_i - \kappa w \cdot dw$ and the constraint $\mathcal{E}^{\text{ex}} + p^{\text{ex}} - T \mathcal{S}^{\text{ex}} = \sum_{i \in \mathfrak{S}} \rho_i g_i^{\text{ex}}$.

(H2) For any $(T, y_1, \dots, y_n) \in (0, \infty)^{1+n}$, such that $\sum_{i \in \mathfrak{S}} y_i = 1$, and any $u \in \mathbb{R}^d$, there exists a threshold density $\rho_m > 0$ such that $(\rho y_1, \dots, \rho y_n, 0, u, T)^t \in O_z$ for $0 < \rho < \rho_m$ and

$$ \lim_{\rho \to 0, w=0} \frac{1}{\rho}(\mathcal{E}^{\text{ex}} - \mathcal{E}^{\text{id}}) = 0, \quad \lim_{\rho \to 0, w=0} \frac{1}{\rho}(\mathcal{P}^{\text{ex}} - \mathcal{P}^{\text{id}}) = 0, \quad \lim_{\rho \to 0, w=0} \frac{1}{\rho}(\mathcal{S}^{\text{ex}} - \mathcal{S}^{\text{id}}) = 0. \quad (6.27) $$

Assumption (H0) describes the regularity class of extended thermodynamic functions, (H1) is associated with the extended Gibbs relation for volumetric quantities, (H2) and with the compatibility with ideal gases. These properties are natural extensions of those concerning the classical thermodynamics $(\mathcal{T}_0-\mathcal{T}_2)$. We now introduce that Property (H3) associated with the thermal stability condition.

(H3) The open set $O_z$ is increasing with temperature $T$, the specific heat is positive $\partial_T \mathcal{E}^{\text{ex}} > 0$, and the capillarity coefficient is positive $\kappa(T) > 0$ over $O_z$.
The thermodynamic stability of the augmented model is now analyzed in the following proposition.

**Proposition 6.1.** Assuming that $(H_0-H_3)$ hold, the Hessian matrix $\partial_{u,u}^2 S^{\text{ex}}$ is negative definite if, and only if, $\Lambda$ is positive definite, $\partial_T \mathcal{E}^{\text{ex}} > 0$, and $\kappa > 0$.

*Proof.* This is a consequence of the expression

$$ \langle \partial_{uu}^2 S^{\text{ex}} x, x \rangle = - \frac{1}{T^2 \partial_T \mathcal{E}^{\text{ex}}} \left( x_{\mathcal{E}} - \left( \partial_{\mathcal{E}} \mathcal{E}^{\text{ex}} - \frac{1}{2} |\mathbf{u}|^2 \right) \mathbf{x}_{\mathcal{E}} - \mathbf{u} \cdot \mathbf{x}_u - (\kappa - T \partial_T \kappa) \mathbf{w} \cdot \mathbf{x}_w \right)^2 \\ - \langle \Lambda \mathbf{x}_{\mathcal{E}}, \mathbf{x}_{\mathcal{E}} \rangle - \frac{\kappa}{T} |\mathbf{x}_w|^2 - \frac{1}{\rho T} |\mathbf{x}_u - \mathbf{u} \mathbf{x}_{\rho}|^2, \qquad (6.28) $$

where $\mathcal{E} = (\rho_1, \dots, \rho_n)^t$, $\mathbf{I} = (1, \dots, 1)^t$, and $x = (x_{\mathcal{E}}, x_w, x_u, x_{\mathcal{E}})^t$ denotes an arbitrary vector of $\mathbb{R}^n$ with $x_{\mathcal{E}} = (x_1, \dots, x_n)^t$. This expression then directly yields that stability is equivalent to $\partial_T \mathcal{E} > 0$, $\kappa > 0$, and $\Lambda$ is positive definite, keeping in mind that $T$ and $\rho$ are positive. $\square$

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^{\text{ex}}$ keeping in mind that thermal stability holds with $(\mathcal{T}_3)$ and $(H_3)$. Extra assumptions may also be introduced concerning the normal variable in the following form:

(H4) The map $(\rho_1, \dots, \rho_n, T) \mapsto (\rho, \frac{g_2 - g_1}{T}, \dots, \frac{g_n - g_1}{T}, T)^t$ is a $C^{\gamma-1}$ diffeomorphism from $O_z$ onto an open set $O_w$.

Assuming $(H_0-H_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_{uu}^2 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 $a_{il}$ the number of $l$th atom in the $i$th species, $\mathfrak{A} = \{1, \dots, n^a\}$ the set of atom indices, and $n^a \ge 1$ the number of atoms—or elements—in the mixture. The atomic vectors $a_l$, $l \in \mathfrak{A}$, are defined by $a_l = (a_{1l}, \dots, a_{nl})^t$. The assumptions for chemistry terms (CH) and transport coefficients (TRP) may also be written in the following form.

(CH) The atom conservation relations $\langle v_j, a_l \rangle = 0$ hold for $j \in \mathfrak{R}, l \in \mathfrak{A}$ and $m = \sum_{l \in \mathfrak{A}} \tilde{m}_l a_l$. The rate constants $K_j$, $j \in \mathfrak{R}$ are $C^\gamma$ positive functions of $T > 0$.

(TRP) The coefficients $v, \eta$, and the matrix $L$ are $C^\gamma$ functions over $O_z$. We have $\eta > 0$, $v \ge 0$, $v + \eta(1 - \frac{2}{d}) > 0$, $L$ is symmetric positive semi-definite, and $N(L) = \mathbb{R}(1, \dots, 1, 0)^t$.

Under the assumptions $(H_0-H_4)$ and (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 $κ$ 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 $ε$:

$$p = \bar{p} + \epsilon^2 \tilde{p}, \qquad (7.1)$$

with $\bar{p}$ denoting the zeroth-order thermodynamic pressure and $\tilde{p}$ the fluid perturbation. The classic thermodynamic properties may then be evaluated at the zeroth-order pressure $\bar{p}$. The zeroth order pressure tensor is then in the form $\bar{\mathcal{P}} = \bar{p}I - \frac{1}{2}\kappa|\nabla\rho|^2I + \kappa\nabla\rho\otimes\nabla\rho - \kappa\rho\Delta\rho I$, and substituting this expansion in the momentum equation, it is found at zeroth-order that

$$\nabla \cdot (\bar{p} I - \frac{1}{2}\kappa|\nabla\rho|^2 I + \kappa\nabla\rho\otimes\nabla\rho - \kappa\rho\Delta\rho I) = 0. \qquad (7.2)$$

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 $\bar{p}$ plays a role in the momentum equation, and from $p^{\text{ex}} = p - \frac{1}{2}\kappa|\nabla\rho|^2$ it is obtained that $p^{\text{ex}} = \bar{p} - \frac{1}{2}\kappa|\nabla\rho|^2 + \epsilon^2\bar{p}$. For standard fluids where $κ = 0$, we then recover that the ambient pressure $\bar{p}$ is spatially uniform [8] since $\nabla \cdot (\bar{p} I) = \nabla \bar{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 \mathcal{S}$, $\rho = \rho(t, y)$, $\bar{p} = \bar{p}(t, y)$, and $T = T(t, y)$. In this framework, the zeroth order pressure tensor $\bar{\mathcal{P}}$ is in the form $\bar{\mathcal{P}} = (\bar{p} - \frac{1}{2}\kappa\rho'^2 - \kappa\rho\rho'')I + \kappa\rho'^2 e_y \otimes e_y$, where $e_y$ denotes the unit vector normal to the interface, $\rho' = \partial_y\rho(t, y)$, and $\rho'' = \partial_{yy}\rho = \Delta\rho(t, y)$. From the normal component of the zeroth order momentum, we deduce that $\bar{p} + \frac{1}{2}\kappa\rho'^2 - \kappa\rho\rho''$ is a constant denoted by $p^{\infty} = \bar{p} + \frac{1}{2}\kappa\rho'^2 - \kappa\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

$$\bar{\mathcal{P}} = (p^{\infty} - \kappa\rho'^2)(I - e_y \otimes e_y) + p^{\infty}e_y \otimes e_y,$$

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} \kappa\rho'^2 dy$ may then be interpreted as surface tension [25, 27, 28].

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

$$\rho \partial_t h + \rho u \nabla \cdot h + \nabla \cdot q = \partial_t \bar{p} + u \cdot \nabla \bar{p}, \qquad (7.3)$$

with a term $u \cdot \nabla \bar{p}$ that is not anymore a priori negligible as for standard fluids. Note incidentally that there is a typographic error in reference [29] with $Q$ instead of $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 $\mathbf{u}$ are denoted by $\mathbf{u} = (\mathbf{u}, \mathbf{v})^T$. The solution is assumed to have the following self-similar structure [8, 29, 30, 76, 92, 93]

$$ 
\begin{aligned}
T &= T(t, y), & \rho &= \rho(t, y), & \bar{p} &= \bar{p}(t, y), \\
\mathbf{u} &= x\tilde{\mathbf{U}}(t, y), & \mathbf{v} &= \mathbf{v}(t, y), & \mathbf{y}_i &= \mathbf{y}_i(t, y), \\
\tilde{p} &= -\frac{1}{2}Jx^2 + \tilde{p}(t, y), & & &
\end{aligned}
$$

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 $J_i = J_i e_y$, $i \in \mathcal{S}$, and $q = q e_y$ and $J_i(t, y)$, $i \in \mathcal{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

$$ \partial_t \rho + 2 \delta \rho \tilde{\mathbf{U}} + \partial_y (\rho \mathbf{v}) = 0, \qquad (7.4) $$

$$ \rho \partial_t y_i + \rho v \partial_y y_i + \partial_y J_i = 0, \quad i \in \mathcal{S}, \qquad (7.5) $$

$$ \rho \partial_t \tilde{\mathbf{U}} + \rho \tilde{\mathbf{U}}^2 + \rho v \partial_y \tilde{\mathbf{U}} - J + \partial_y (\eta \partial_y \tilde{\mathbf{U}}) = 0, \qquad (7.6) $$

$$ \rho \partial_t h + \rho v \partial_y h - v \partial_y \bar{p} + \partial_y q = 0, \qquad (7.7) $$

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

$$ J_i = - \sum_{j \in \mathcal{S}} L_{ij} \partial_y \left( \frac{g_j}{T} \right) - L_{ie} \partial_y \left( \frac{-1}{T} \right), \quad i \in \mathcal{S}, \qquad (7.8) $$

$$ q = - \sum_{i \in \mathcal{S}} L_{ei} \partial_y \left( \frac{g_i}{T} \right) - L_{ee} \partial_y \left( \frac{-1}{T} \right), \qquad (7.9) $$

and it has been established in Section 7.1 that

$$ \bar{p} = p^{\infty} - \frac{1}{2} \kappa (\partial_y \rho)^2 + \kappa \rho \partial_y^2 \rho, \qquad (7.10) $$

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

$$ T(-\infty) = T^{-}, \qquad T(+\infty) = T^{+}, \qquad (7.11) $$

$$ y_i(-\infty) = y_i^{-}, \qquad y_i(+\infty) = y_i^{+}, \qquad (7.12) $$

$$ \tilde{u}(-\infty) = \alpha \sqrt{\rho^{+}/\rho^{-}}, \qquad \tilde{u}(+\infty) = \alpha, \qquad (7.13) $$
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

$$ \alpha = (J/\rho^+)^{1/2}. \qquad (7.14) $$

The stagnation plane is finally located at the origin with

$$ v(0) = 0. \qquad (7.15) $$

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

$$ 2\delta\rho\tilde{u} + \partial_y(\rho v) = 0, \qquad (7.16) $$

$$ \rho v \partial_y y_i + \partial_y J_i = 0, \quad i \in \mathcal{S}, \qquad (7.17) $$

$$ \rho \tilde{u}^2 + \rho v \partial_y \tilde{u} - J + \partial_y (\eta \partial_y \tilde{u}) = 0, \qquad (7.18) $$

$$ \rho v \partial_y h - v \partial_y \bar{p} + \partial_y q = 0, \qquad (7.19) $$

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

$$ \mu_i = \log x_i + \mu_i^{\text{sm}}, \qquad (7.20) $$

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

$$ \mu_i^{\text{sm}} = \frac{m_i g_i^{\text{id}\star}}{RT} + \log\left(\frac{RT}{(v-b)p^{\text{st}}m}\right) + \sum_{j \in \mathcal{S}} \frac{y_j m_i b_j}{m_j (v-b)} - \frac{m_i a b_i}{RT b(v+b)} \\ + \frac{a b_i - b \partial_{y_i} a}{b^2 RT} \log\left(1 + \frac{b}{v}\right), \qquad (7.21) $$

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

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

$$ L = \begin{pmatrix} \mathcal{D} & \mathcal{D}h \\ (\mathcal{D}h)^t & \lambda T^2 + \langle \mathcal{D}h, h \rangle \end{pmatrix}, \qquad (7.22) $$

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

$$ D_{ij} = \rho y_i y_j \frac{m}{R} D_{ij}, \quad h_i = h_i + RT \frac{\tilde{\chi}_i}{m_i}, \qquad (7.23) $$
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 = \partial_{y_i} h$, $i \in \mathfrak{S}$, where the derivative is thus taken at constant pressure, and $\tilde{\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–111] using

$$D = (\Delta + y \otimes y)^{-1} - \mathbb{I} \otimes \mathbb{I},$$

where $y$ denotes the mass fraction vector, $\mathbb{I} = (1, \dots, 1)^t$ the vector with unit components, $\Delta$ the Stefan-Maxwell matrix

$$\Delta_{kk} = \sum_{l \neq k} \frac{x_k x_l}{\mathcal{D}_{kl}^{\text{bin}}}, \quad \Delta_{kl} = - \frac{x_k x_l}{\mathcal{D}_{kl}^{\text{bin}}}, \quad k \neq l,$$

$x_k$ the mole fraction of the $k$th species, and $\mathcal{D}_{kl}^{\text{bin}}$, the binary coefficient for the species pair $(k, l)$, $k \neq l$. The binary diffusion coefficients $\mathcal{D}_{kl}^{\text{bin}}$, $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}^{\text{bin}} = \mathcal{D}_{kl}^{\text{bin id}} / \Upsilon_{kl}$, $k, l \in \mathfrak{S}$, $k \neq l$, where $\mathcal{D}_{kl}^{\text{bin id}}$ denotes the ideal gas binary diffusion coefficient and $\Upsilon_{ij}$ the steric factor

$$\Upsilon_{ij} = 1 + \sum_{k \in \mathfrak{S}} \frac{\pi \eta_k}{12} (8(\sigma_{ik}^3 + \sigma_{jk}^3) - 6(\sigma_{ik}^2 + \sigma_{jk}^2)\sigma_{ij} - 3(\sigma_{ik}^2 - \sigma_{jk}^2)^2 \sigma_{ij}^{-1} + \sigma_{ij}^3).$$

We have denoted here by $\sigma_{ij}$ the collision diameter of the species pair $(i, j)$, $\eta_k = \rho_k N/m_k$ the particle number density of the $k$th species, and by $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 $\tilde{\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}y$ [8, 107–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 = \partial_{y_i} h$, $i \in \mathfrak{S}$. We indeed have the identity

$$h_k = \partial_{y_k} h = \partial_{y_k} h - \partial_{\rho} h \frac{\partial_{y_k} \bar{p}}{\partial_{\rho} \bar{p}}, \qquad (7.24)$$

so that at a mechanical thermodynamic unstable point where $\partial_{\rho} \bar{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 $h$. To this aim, the derivative $\partial_{y_k} \bar{p}$ is rewritten in the form

$$\partial_{y_k} \bar{p} = \frac{\rho_m}{m_k} \partial_{\rho} \bar{p} + \mathcal{R}_k, \qquad (7.25)$$

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

$$\mathcal{R}_k = \left\{ \frac{RT}{(v-b)^2} + \frac{ma}{v(v+b)^2} \right\} \sum_{j \in \mathfrak{S}} y_j \left( \frac{b_k}{m_j} - \frac{b_j}{m_k} \right) + \frac{2m}{v(v+b)} \sum_{j,l \in \mathfrak{S}} y_j y_l \alpha_l \left( \frac{\alpha_j}{m_k} - \frac{\alpha_k}{m_j} \right).$$
The interest of this formulation (7.25) is that $\mathcal{R}_k$ vanishes for the pure species state $y_k = 1$ and $y_l = 0$, $l \in \mathfrak{S}, l \neq k$. The species enthalpies are then written in the modified form $\tilde{h}_k, k \in \mathfrak{S}$, with

$$ \tilde{h}_k = \partial_{y_k} h - \frac{\rho m}{m_k} \partial_{\rho} h - \frac{\mathcal{R}_k}{\mathfrak{f}(\partial_{\rho} \bar{p})} \partial_{\rho} h, $$

where $\mathfrak{f}$ is a smooth approximation of the function $\mathfrak{f}(x) = \max(m, x)$ for a positive constant $m > 0$ which depends on the mixture upon consideration. Such revised enthalpies are such that $\tilde{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–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₂H₆-N₂ 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 \mathcal{S}$, associated with the SRK equation of state (5.1) and (5.2), can be evaluated from the critical data [80]. Assuming that the ith species is chemically stable, the attractive parameter $\alpha_{c,i}$ at the critical temperature $T_{c,i}$ of Proposition 5.4 and the repulsive parameter $b_i$ are evaluated in the form

$$ \alpha_{c,i}^2 = a_{c,i} = 0.42748 \frac{R^2 T_{c,i}^2}{m_i^2 p_{c,i}}, \quad b_i = 0.08664 \frac{RT_{c,i}}{m_i p_{c,i}}, \qquad (8.1) $$

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

$$ s_i = 0.48508 + 1.5517\bar{\omega}_i - 0.151613\bar{\omega}_i^2, \qquad (8.2) $$

where $\bar{\omega}_i$ is the ith species acentric factor [113]. The attractive and repulsive parameters and the acentric factors that have been used for ethane C₂H₆ and nitrogen N₂ are presented in Table 1.

<table><caption>**Table 1.** Critical state and acentric factor of ethane and nitrogen from NIST.</caption><thead><tr><th>Species</th><th><i>T</i><sub>c,i</sub> (K)</th><th><i>p</i><sub>c,i</sub> (bar)</th><th><i>ω</i><sub>i</sub></th></tr></thead><tbody><tr><td>C<sub>2</sub>H<sub>6</sub></td><td>305.3</td><td>49.0</td><td>0.098</td></tr><tr><td>N<sub>2</sub></td><td>126.1</td><td>39.98</td><td>0.0372</td></tr></tbody></table>

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 ith species is a Lennard-Jones gas, however, it is possible to estimate [74,75] a pseudo-critical state and to deduce that

$$ a_i(T_{c,i}) = (5.55 \pm 0.12) \frac{N^2 \epsilon_i \sigma_i^3}{m_i^2}, \quad b_i = (0.855 \pm 0.018) \frac{N \sigma_i^3}{m_i}, \qquad (8.3) $$

where N is the Avogadro number, and $\sigma_i$ and $\epsilon_i$ are the molecular diameter and Lennard-Jones potential well depth of the ith 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–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_{\text{ref}}$ and temperature $T_{\text{ref}}$, for a given pure species with $n = 1$, $\mathcal{S} = \{1\}$, $y_1 = 1$, the corresponding specific volume(s) $v$ may be obtained by solving the cubic equation $p(T_{\text{ref}}, v, y_1) = p_{\text{ref}}$. Denoting by $v_{\text{inf}}$ and $v_{\text{sup}}$ the lowest and highest solutions, the fluid is liquid with volume $v_{\text{inf}}$ if $g(T_{\text{ref}}, v_{\text{inf}}, y_1) < g(T_{\text{ref}}, v_{\text{sup}}, y_1)$, gaseous with volume $v_{\text{sup}}$ if $g(T_{\text{ref}}, v_{\text{sup}}, y_1) < g(T_{\text{ref}}, v_{\text{inf}}, y_1)$, and both states coexist when $g(T_{\text{ref}}, v_{\text{inf}}, y_1) = g(T_{\text{ref}}, v_{\text{sup}}, y_1)$. The constant pressure heat capacity is then taken to be the partial derivative with respect to $T_{\text{ref}}$ of $h(T_{\text{ref}}, v_{\text{inf}}(T_{\text{ref}}, p_{\text{ref}}), y_1)$ or $h(T_{\text{ref}}, v_{\text{sup}}(T_{\text{ref}}, p_{\text{ref}}), y_1)$ 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].

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Constant pressure specific heat of pure nitrogen; — SRK; --- Ideal gas; + + + experimental measurements from Youglove [121] and Youglove and Ely [122].</figcaption></figure>

The saturation pressure $p_{\text{sat}}$ of nitrogen as a function of temperature $T$ is obtained by solving the nonlinear equation $g(T_{\text{ref}}, v_{\text{inf}}, y_1) = g(T_{\text{ref}}, v_{\text{sup}}, y_1)$, where $v_{\text{inf}}$ and $v_{\text{sup}}$ are the lowest and highest solutions of the cubic equation $p(T_{\text{ref}}, v, y_1) = p_{\text{ref}}$, and we then have $p_{\text{ref}} = p_{\text{sat}}(T_{\text{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.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Saturation pressure of pure nitrogen; ____ SRK; --- mechanical stability limits; +++ experimental measurements form Younglove and Ely [122].</figcaption></figure>

### 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 $\mathcal{S} = \{1, 2\} = \{C_2H_6, 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, y) = 0$ and the mechanical stability domain—of the mixture—by solving the equation $\partial_v p(T, v, y) = 0$. We have used a continuation in the variable $v$ for solving the system

$$ \left\{ \begin{array}{l} T = T_{\text{ref}}, \\ \det \Lambda(T, v, y), \\ y_1 + y_2 = 1, \end{array} \right. \qquad (8.4) $$

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 $\Lambda^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 $(y_{C_2H_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.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Stability limit domains of ethane-nitrogen mixtures at T = 200 K; ____ thermodynamic stability; --- mechanical stability limits.</figcaption></figure>

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

$$ egin{cases} T^\# = T^b = T_{\text{ref}}, \\ p(T^\#, v^\#, y^\#) = p(T^b, v^b, y^b), \\ g_1(T^\#, v^\#, y^\#) = g_1(T^b, v^b, y^b), \\ g_2(T^\#, v^\#, y^\#) = g_2(T^b, v^b, y^b), \\ y_1^\# + y_2^\# = y_1^b + y_2^b = 1. \end{cases} \quad (8.5) $$

The equilibrium curve between ethane and nitrogen at T = 220 K is presented in Figure 4 in the $(y_{C_2H_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.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Equilibrium between ethane and nitrogen at T = 220 K; ____ numerical simulation; + + + experimental measurements from Youglove [121] and Youglove and Ely [122].</figcaption></figure>
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_c$ of the system of equations [128]

$$ \left\{ \begin{array}{l} \partial_{x_1} \mathcal{G}(T_c, v_c, y_c) = 0, \\ \partial_{x_1}^2 \mathcal{G}(T_c, v_c, y_c) = 0, \\ y_{c1} + y_{c2} = 1. \end{array} \right. \qquad (8.6) $$

Only the relevant branch starting from the critical point of pure ethane given by $T_c = 305.3$ K and the $p_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 III according to the classification of Konynenburg and Scott [112, 113]. Binary $O_2-H_2O$ mixtures also have such a type III 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 $(y_{C_2H_6}, p)$ is also presented in Figure 6, and we notably observe a turning point with respect to the ethane mass fraction $y_{C_2H_6}$ that is easily taken into account by using continuation methods.

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Diagram of critical points for $C_2H_6-N_2$ mixtures in the $(T, p)$ plane; — numerical simulation; * * * experimental measurements from Eakin et al. [124]; △△△ experimental measurements from Stryjek et al. [125]; + + + experimental measurements from Wisotzki and Schneider [126].</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 6. Diagram of critical points for C₂H₆-N₂ mixtures in the (yC₂H₆, p) plane.</figcaption></figure>

## 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

$$T^{-} = 400 \text{ K}, \qquad T^{+} = 250 \text{ K}, \qquad (8.7)$$

$$\alpha = 100 \text{ s}^{-1}, \qquad p^{\infty} = 3 \text{ MPa.} \qquad (8.8)$$

The capillarity coefficient is varied to be

$$\kappa = 5.10^{-4}, 10^{-3}, 10^{-2}, 10^{-1} \text{ cm}^{7}\text{g}^{-1}\text{s}^{-2}.$$

For these pseudo-vaporizing interfaces, C₂H₆ is the only chemical species and no species equations are solved. The mass density ρ(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.

<figure><img src="image_9.png" /><figcaption>Figure 7. Mass density of ethane fronts for $p^{\infty} = 3$ MPa and various capillarity coefficients.</figcaption></figure>
The capillary parameter $\kappa$ has an important impact on the interface structure. The lower the value of $\kappa$, the sharper the interface, and when $\kappa \to 0$, the infinitely thin interface—discontinuous—solution is recovered. The corresponding local pressure $\bar{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 $\kappa$ [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_2H_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.

# References

1. E. A. Guggenheim, *Thermodynamics*, Amsterdam: North Holland, 1962.

2. N. Z. Shapiro, L. S. Shapley, Mass action law and the Gibbs free energy function, *SIAM J. Appl. Math.*, **13** (1965), 353–375. https://doi.org/10.1137/0113020
3. R. Aris, Prolegomena to the rational analysis of systems of chemical reactions, *Archiv. Rat. Mech. Anal.*, **19** (1965), 81–99. https://doi.org/10.1007/BF00282276

4. C. Truesdell, *Rational thermodynamics*, New York: McGraw-Hill, 1969.
https://doi.org/10.1007/978-1-4612-5206-1

5. F. J. Krambeck, The mathematical structure of chemical kinetics, *Arch. Rational Mech. Anal.*, **38** (1970), 317–347. https://doi.org/10.1007/BF00251527

6. J. Pousin, *Modélisation et analyse numérique de couches limites Réactives d'air*, Doctorat es Sciences, Ecole Polytechnique Fédérale de Lausanne, 1993.

7. P. Helluy, H. Mathis, Pressure laws and fast Legendre transform, *Math. Mod. Meth. Appl. S.*, **21** (2010), 745–775. https://doi.org/10.1142/S0218202511005209

8. V. Giovangigli, *Multicomponent flow modeling*, Boston: Birkhäuser, 1999.
https://doi.org/10.1007/978-1-4612-1580-6

9. V. Giovangigli, L. Matuszewski, Supercritical fluid thermodynamics from equations of state, *Phys. D*, **241** (2012), 649–670. https://doi.org/10.1016/j.physd.2011.12.002

10. K. O. Friedrichs, P. D. Lax, Systems of conservation laws with a convex extension, *P. Natl. Acad. Sci. USA*, **68** (1971), 1686–1688. https://doi.org/10.1073/pnas.68.8.1686

11. T. Ruggeri, Thermodynamics and symmetric hyperbolic systems, *Rend. Semin. Mat. U. Torino*, 1988, 167–183.

12. E. Godlevski, P. A. Raviart, *Numerical approximation of hyperbolic systems of conservation laws*, New York: Springer, **118** (1996). https://doi.org/10.1007/978-1-4612-0713-9

13. C. Dafermos, *Hyperbolic conservation laws in continuum physics*, Heidelberg: Springer, 2000.
https://doi.org/10.1007/978-3-662-22019-1

14. A. I. Vol'pert, S. I. Hudjaev, On the Cauchy problem for composite systems of nonlinear differential equations, *Math. USSR Sbornik*, **16** (1972), 517–544.
https://doi.org/10.1070/sm1972v016n04abeh001438

15. S. Kawashima, Y. Shizuta, On the normal form of the symmetric hyperbolic-parabolic systems associated with the conservation laws, *Tôhoku Math. J.*, **40** (1988), 449–464.
https://doi.org/10.2748/tmj/1178227986

16. V. Giovangigli, M. Massot, Asymptotic stability of equilibrium states for multicomponent reactive flows, *Math. Mod. Meth. Appl. S.*, **8** (1998), 251–297.
https://doi.org/10.1142/S0218202598000123

17. V. Giovangigli, M. Massot, Entropic structure of multicomponent reactive flows with partial equilibrium reduced chemistry, *Math. Method. Appl. Sci.*, **27** (2004), 739–768.
https://doi.org/10.1002/mma.429

18. V. Giovangigli, L. Matuszewski, Mathematical modeling of supercritical multicomponent reactive fluids, *Math. Mod. Meth. Appl. S.*, **23** (2013), 2193–2251.
https://doi.org/10.1142/S0218202513500309

19. V. Giovangigli, W. A. Yong, Asymptotic stability and relaxation for fast chemistry fluids, *Nonlinear Anal.*, **159** (2017), 208–263. https://doi.org/10.1016/j.na.2017.02.025
20. J. D. V. der Waals, Thermodynamische theorie der capillariteit in de onderstelling van continue Dichtheidsverandering, *J. Müller*, **20** (1979), 197–244. https://doi.org/10.1007/BF01011513

21. J. D. V. der Waals, Thermodynamisch theorie der kapillariät unter voraussetzung stetiger dichteanderung, *Z. Phys. Chem.*, **13** (1894), 657–725.

22. D. J. Korteweg, Sur la Forme que Prennent les Equations du Mouvement Fluide si l'on tient Compte de Forces Capillaires Causées par les Variations de Densité Considérables mais Continues et sur la Théorie de la Capillarité dans l'Hypothèse d'une Variations Continue de la Densité, *Arch. Neerl. Sci. Exactes*, **6** (1901), 1–20. Available from: https://archive.org/details/archivesnerland261901laha/page/n17/mode/2up.

23. J. E. Dunn, J. Serrin, On the thermomechanics of interstitial working, *Arch. Ration. Mech. An.*, **133** (1985), 95–133. https://doi.org/10.1007/BF00250907

24. P. G. de Gennes, Wetting: Statics and dynamics, *Rev. Mod. Phys.*, **57** (1985), 827–863. https://doi.org/10.1103/RevModPhys.57.827

25. J. S. Rowlinson, B. Widom, *Molecular theory of capillarity*, Courier Corporation, 2013.

26. D. M. Anderson, G. B. McFadden, A. A. Wheeler, Diffuse interface methods in fluid mechanics, *Annu. Rev. Fluid Mech.*, **30** (1998), 139–165. https://doi.org/10.1146/annurev.fluid.30.1.139

27. D. Jamet, *Diffuse interface models in fluid mechanics*, GdR CNRS Documentation, see Pmc. Polytechnique, 2001.

28. D. Jamet, O. Lebaigue, N. Coutris, J. M. Delhaye, The second gradient method for the direct numerical simulation of liquid-vapor flows with phase change, *J. Comput. Phys.*, **169** (2001), 624–651. https://doi.org/10.1006/jcph.2000.6692

29. P. Gaillard, V. Giovangigli, L. Matuszewski, A diffuse interface lox/hydrogen transcritical flame model, *Combust. Theor. Model.*, **20** (2016), 486–520. https://doi.org/10.1080/13647830.2016.1150518

30. P. Gaillard, V. Giovangigli, L. Matuszewski, Erratum for a diffuse interface lox/hydrogen transcritical flame model. Available from: https://hal.science/hal-04198003.

31. D. Nayigizente, S. Ducruix, T. Schmitt, Development of an interface thickening method for the direct numerical simulation of compressible liquid-vapor flows in the framework of the second gradient theory, *Phys. Fluids*, **33** (2021), 052119. https://doi.org/10.1063/5.0048715

32. Y. Le Calvez, *Modélisation mathématique et simulation numériques des mélanges fluides avec interfaces diffuses*, PhD thesis, (in preparation), 2024.

33. S. Gavrilyuk, S. Shugrin, Media with equations of state that depend on derivatives, *J. Appl. Mech. Techn. Ph.*, **37** (1996), 177–189. https://doi.org/10.1007/BF02382423

34. Y. Rocard, Equations hydrodynamiques avec termes capillaires, Théorie de la tension superficielle, *J. Phys. Radium*, **4** (1933), 533–548. https://doi.org/10.1051/jphysrad:01933004010053300

35. P. Barbante, A. Frezzotti, A comparison of models for the evaporation of a Lennard-Jones fluid, *Eur. J. Mech. B-Fluid.*, **64** (2017), 69–80. https://doi.org/10.1016/j.euromechflu.2017.01.020

36. V. Giovangigli, Kinetic derivation of diffuse-interface fluid models, *Phys. Rev. E*, **102** (2020), 012110. https://doi.org/10.1103/PhysRevE.102.012110
37. J. W. Cahn, J. E. Hilliard, Free energy of a non uniform system I, Interfacial free energy, *J. Chem. Phys.*, **28** (1958), 258–267. https://doi.org/10.1063/1.1744102

38. J. W. Cahn, J. E. Hilliard, Free energy of a non uniform system II, thermodynamic basis, *J. Chem. Phys.*, **30** (1959), 1121–1124. https://doi.org/10.1063/1.1730145

39. F. Falk, Cahn-Hilliard theory and irreversible thermodynamics, *J. Non-Equil. Thermody.*, **17** (1992), 53–65. https://doi.org/10.1515/jnet.1992.17.1.53

40. J. Kim, J. Lowengrub, Phase field modeling and simulation of three-phase flows, *Inter. Free Boundary*, **7** (2005), 435–466. https://doi.org/10.4171/IFB/132

41. H. W. Alt, The entropy principle for interfaces, fluids and solids, *Adv. Math. Sci. Appl.*, **19** (2009), 585–663.

42. H. Abels, H. Garcke, G. Grün, Thermodynamically consistent, frame indifferent diffuse interface models for incompressible two-phase flows with different densities, *Math. Mod. Meth. Appl. S.*, **22** (2012), 1150013. https://doi.org/10.1142/S0218202511500138

43. Z. Guo, P. Lin, A thermodynamically consistent phase-field model for two-phase flows with thermocapillarity effects, *J. Fluid Mech.*, **766** (2015), 226–271. https://doi.org/10.1017/jfm.2014.696

44. V. Giovangigli, Kinetic derivation of Cahn-Hilliard fluid models, *Phys. Rev. E*, **104** (2021), 054109. https://doi.org/10.1103/PhysRevE.104.054109

45. H. G. Lee, J. Yang, S. Kim, J. Kim, Modeling and simulation of droplet evaporation using a modified Cahn-Hilliard equation, *Appl. Math. Comput.*, **390** (2021), 125591. https://doi.org/10.1016/j.amc.2020.125591

46. C. Wang, S. M. Wise, A thermodynamically-consistent phase field crystal model of solidification with heat flux, *J. Math. Study*, **55** (2022), 337–357. https://doi.org/10.4208/jms.v55n4.22.01

47. F. Wang, H. Zhang, Y. Wu, B. Nestler, A thermodynamically consistent diffuse interface model for the wetting phenomenon of miscible and immiscible ternary fluids, *J. Fluid Mech.*, **970** (2023), A17. https://doi.org/10.1017/jfm.2023.561

48. H. Zhang, F. Wang, B. Nestler, Multi-component electro-hydro-thermodynamic model with phase-field method. I. Dielectric, *J. Comput. Phys.*, **505** (2024), 112907. https://doi.org/10.1016/j.jcp.2024.112907

49. A. Miranville, *The Cahn-Hilliard equation: Recent advances and applications*, CBMS-NSF Regional Conference Series in Applied Mathematics 95, Philadelphia: Society for Industrial and Applied Mathematics (SIAM), 2019. https://doi.org/10.1137/1.9781611975925

50. E. Bretin, S. Masnou, E. Oudet, Phase-field approximations of the Willmore functional and flow, *Numer. Math.*, **131** (2015), 115–171. https://doi.org/10.1007/s00211-014-0683-4

51. E. Benilov, The multicomponent diffuse-interface model and its application to water/air interfaces, *J. Fluid Mech.*, **954** (2023), A41. https://doi.org/10.1017/jfm.2022.1032

52. M. Dressler, B. Edwards, C. Öttinger, Macroscopic thermodynamics of flowing polymeric liquids, *Rheol. Acta*, **38** (1999), 117–136. https://doi.org/10.1007/s003970050162

53. D. Bruno, V. Giovangigli, Relaxation of internal temperature and volume viscosity, *Phys. Fluids*, **23** (2011), 093104. https://doi.org/10.1063/1.3640083
54. E. V. Kustova, E. A. Nagnibeda, On a correct description of a multi-temperature dissociating CO₂ flow, *Chem. Phys.*, **321** (2006), 293–310. https://doi.org/10.1016/j.chemphys.2005.08.026

55. R. H. Fowler, *Statistical mechanics*, Cambridge: Cambridge University Press, 1936.

56. J. H. Ferziger, H. G. Kaper, *Mathematical theory of transport processes in gases*, Amsterdam: North Holland, 1972.

57. J. Keizer, *Statistical thermodynamics of nonequilibrium processes*, New York: Springer-Verlag, 1987. https://doi.org/10.1007/978-1-4612-1054-2

58. K. Laasonen, S. Wonczak, R. Strey, A. Laaksonena, Molecular dynamics simulations of gas-liquid nucleation of Lennard-Jones fluid, *J. Chem. Phys.*, **113** (2000), 9741–9747. https://doi.org/10.1063/1.1322082

59. S. Chen, G. D. Doolen, Lattice Boltzmann method for fluid flows, *Annu. Rev. Fluid Mech.*, **30** (1998), 329–364. https://doi.org/10.1146/annurev.fluid.30.1.329

60. S. R. de Groot, P. Mazur, *Non-equilibrium thermodynamics*, Mineola: Dover publications, 1984.

61. C. Öttinger, *Beyond equilibrium thermodynamics*, Hoboken: John Wiley and Sons, 2005. https://doi.org/10.1002/0471727903

62. L. J. Gillespie, Equilibrium pressures of individual gases in mixtures and the mass-action law for gases, *J. Am. Chem. Soc.*, **47** (1925), 305–312. https://doi.org/10.1021/ja01679a003

63. M. Benedict, G. B. Webb, L. C. Rubin, An empirical equation for thermodynamic properties of light hydrocarbons and their mixtures, *J. Chem. Phys.*, **8** (1940), 334–345. https://doi.org/10.1063/1.1750658

64. J. A. Beattie, The computation of the thermodynamic properties of real gases and mixtures of real gases, *Chem. Rev.*, **18** (1948), 141–192. https://doi.org/10.1021/cr60137a010

65. O. Redlich, J. N. S. Kwong, On the thermodynamics of solutions. V. An equation of state. Fugacities of gaseous solutions, *Chem. Rev.*, **44** (1949), 233–244. https://doi.org/10.1021/cr60137a013

66. G. S. Soave, Equilibrium constants from a modified Redlich-Kwong equation of state, *Chem. Eng. Sci.*, **27** (1972), 157–172. https://doi.org/10.1016/0009-2509(72)80096-4

67. G. S. Soave, An effective modification of the Benedict-Webb-Rubin equation of state, *Fluid Phase Equilbr.*, **164** (1999), 157–172. https://doi.org/10.1016/S0378-3812(99)00252-6

68. D. Y. Peng, D. B. Robinson, A new two-constant equation of state, *Ind. Eng. Chem. Fundam.*, **15** (1976), 59–64. https://doi.org/10.1021/i160057a011

69. M. S. Graboski, T. E. Daubert, A modified Soave equation of state for phase equilibrium calculations. 3. Systems containing hydrogen, *Ind. Eng. Chem. Process. Dev. Des.*, **18** (1979), 300–306. https://doi.org/10.1021/i260070a022

70. E. D. Ozokwelu, J. H. Erbar, An improved Soave-Redlich-Kwong equation of state, *Chem. Eng. Commun.*, **52** (1987), 9–19. https://doi.org/10.1080/00986448708911854

71. K. G. Harstad, R. S. Miller, J. Bellan, Efficient high pressure state equations, *AICHE J.*, **43** (1997), 1605–1610. https://doi.org/10.1002/aic.690430624
72. A. Congiunti, C. Bruno, E. Giacomazzi, *Supercritical combustion properties*, 11th Aerospace Sciences Meeting and Exhibit, AIAA-2003-478, 2003. https://doi.org/10.2514/6.2003-478

73. P. Colonna, P. Silva, Dense gas thermodynamic properties of single and multicomponent fluids for fluid dynamics simulations, *J. Fluid Eng.*, **125** (2003), 414-427. https://doi.org/10.1115/1.1567306

74. W. A. C. Marín, U. E. G. Aconcha, J. D. O. Arango, Comparison of different cubic equations of state and combination rules for predicting residual chemical potential of binary and ternary Lennard-Jones mixtures: Solid-supercritical fluid phase equilibria, *Fluid Phase Equilibr.*, **234** (2005), 42–50. https://doi.org/10.1016/j.fluid.2005.05.014

75. W. A. C. Marín, J. D. O. Arango, U. E. G. Aconcha, C. P. S. Tavera, Thermodynamic derivative properties and densities for hyperbaric gas condensates: SRK equation of state predictions versus Monte Carlo data, *Fluid Phase Equilibr.*, **253** (2007), 147–154. https://doi.org/10.1016/j.fluid.2007.02.004

76. A. M. Saur, F. Behrendt, E. U. Franck, Calculation of high pressure counterflow diffusion flame up to 3000 bar, *Ber. Bunsenges Phys. Chem.*, **97** (1993), 900–908. https://doi.org/10.1002/bbpc.19930970710

77. H. Meng, V. Yang, A unified treatment of general fluid thermodynamics and its application to a preconditioning scheme, *J. Comput. Phys.*, **189** (2003), 277–304. https://doi.org/10.1016/S0021-9991(03)00211-0

78. J. C. Oefelein, Thermophysical characteristics of shear-coaxial LOX-H₂ flames at supercritical pressure, *P. Combust. Inst.*, **30** (2005), 2929–2937. https://doi.org/10.1016/j.proci.2004.08.212

79. G. Ribert, N. Zong, V. Yang, L. Pons, N. Darabiha, S. Candel, Counterflow diffusion flames of general fluids: Oxygen/hydrogen mixtures, *Combust. Flame*, **154** (2008), 319–330. https://doi.org/10.1016/j.combustflame.2008.04.023

80. V. Giovangigli, L. Matuszewski, F. Dupoirieux, Detailed modeling of planar transcritical H₂-O₂-N₂ flames, *Combust. Theor. Model.*, **15** (2011), 141–182. https://doi.org/10.1080/13647830.2010.527016

81. P. Gaillard, V. Giovangigli, L. Matuszewski, Nonmixing layers, *Phys. Rev. Fluids*, **1** (2016), 084001. https://doi.org/10.1103/PhysRevFluids.1.084001

82. S. Gavrilyuk, H. Gouin, Symmetric form of governing equations for capillary fluids, *arXiv Preprint*, 2008. https://doi.org/10.48550/arXiv.0802.1670

83. S. B. Gavage, R. Danchin, L. Mazet, D. Jamet, Structure of Korteweg models and stability of diffuse interfaces, *Interface. Free Bound.*, **7** (2005), 371–414. https://doi.org/10.4171/IFB/130

84. D. Bresch, F. Couderc, P. Noble, J. P. Vila, A generalization of the quantum Bohm identity: Hyperbolic CFL condition for Euler-Korteweg equations, *C. R. Math. Acad. Sci. Paris*, **354** (2008), 39–43. http://dx.doi.org/10.1016/j.crma.2015.09.020

85. D. Bresch, V. Giovangigli, E. Zatorska, Two-velocity hydrodynamics in fluid mechanics: Part I Well posedness for zero Mach number systems, *J. Math. Pure. Appl.*, **104** (2015), 762–800. https://doi.org/10.1016/j.matpur.2015.05.003

86. M. Kotschote, Dynamics of compressible non-isothermal fluids of new-newtonian Korteweg type, *SIAM J. Math. Anal.*, **44** (2012), 74–101. https://doi.org/10.1137/110821202
87. V. Giovangigli, Y. Le Calvez, F. Nabet, Symmetrization and local existence of strong solutions for diffuse interface fluid models, *J. Math. Fluid Mech.*, **25** (2023), 82. https://doi.org/10.1007/s00021-023-00825-4

88. M. R. Marcelin, Sur la mécanique des phénomènes irréversibles, *C. R. Acad. Sci. Paris*, 1910, 1052-1055.

89. M. R. Marcelin, *Contribution à l'étude de la cinétique physico-chimique*, Thèses de la Faculté des Sciences de Paris, 1914.

90. A. Ern, V. Giovangigli, The Kinetic equilibrium regime, *Phys. A*, **260** (1998), 49-72.
https://doi.org/10.1016/S0378-4371(98)00303-3

91. V. Giovangigli, *Solutions for models of chemically reacting mixtures*, Handbook of Mathematical Analysis in Mechanics of Viscous Fluids, Yoshikazu Giga and Antonín Novotný Editors, Springer-Verlag, 2018, 2979-3030. https://doi.org/10.1007/978-3-319-10151-4

92. L. Pons, N. Darabiha, S. Candel, G. Ribert, V. Yang, Mass transfer and combustion in transcritical non-premixed counterflows, *Combust. Theor. Model.*, **13** (2009), 57-81.
https://doi.org/10.1080/13647830802368821

93. V. Giovangigli, L. Matuszewski, Numerical simulation of transcritical strained laminar flames, *Combust. Flame*, **159** (2012), 2829-2840. https://doi.org/10.1016/j.combustflame.2012.05.011

94. R. J. Kee, F. M. Rupley, J. A. Miller, *Chemkin II: A Fortran chemical kinetics package for the analysis of gas phase chemical kinetics*, Livermore: Sandia National Lab. (SNL-CA), 1989.
https://doi.org/10.2172/5681118

95. D. G. Goodwin, R. L. Speth, H. K. Moffat, B. W. Weber, Cantera: An object-oriented software toolkit for chemical kinetics, thermodynamics, and transport processes, *Zenodo*, 2018.
https://doi.org/10.5281/zenodo.4527812

96. P. Deuflhard, *Newton methods for nonlinear problems*, Berlin: Springer, 2004.

97. H. B. Keller, Numerical Solution of bifurcation and nonlinear eigenvalue problems, *Appl. Bifurcat. Theory*, 1977, 359-384.

98. M. D. Smooke, Solution of burner stabilized premixed laminar flames by boundary value methods, *J. Comput. Phys.*, **48** (1982), 72-105. https://doi.org/10.1016/0021-9991(82)90036-5

99. M. D. Smooke, The computation of laminar flames, *P. Combust. Inst.*, **34** (2013), 65-98.
https://doi.org/10.1016/j.proci.2012.09.005

100. E. S. Oran, J. P. Boris, *Numerical simulation of reactive flows*, Cambridge: Cambridge University Press, 2001.

101. V. Giovangigli, M. D. Smooke, Adaptive continuation algorithms with application to combustion problems, *Appl. Numer. Math.*, **5** (1989), 305-331. https://doi.org/10.1016/0168-9274(89)90013-5

102. V. Giovangigli, N. Darabiha, *Vector computers and complex chemistry combustion*, Mathematical Modeling in Combustion and Related Topics, C. Brauner and C. S. Lainé Eds., Dordrecht: Springer Netherlands, **140** (1988), 491-503. https://doi.org/10.1007/978-94-009-2770-4

103. A. Ern, V. Giovangigli, Optimized transport algorithms for flame codes, *Combust. Sci. Tech.*, **118** (1996), 387-395. https://doi.org/10.1080/00102209608951988
104.L. Matuszewski, *Modélisation et simulation numérique des phénomènes de combustion en régime supercritique*, PhD Thesis, University Paris, 2011.

105.A. Ern, V. Giovangigli, Projected iterative algorithms with application to multicomponent transport, *Linear Algebra Appl.*, **250** (1997), 289–315. https://doi.org/10.1016/0024-3795(95)00502-1

106.V. Giovangigli, Multicomponent transport algorithms for partially ionized plasmas, *J. Comput. Phys.*, **229** (2010), 4117–4142. https://doi.org/10.1016/j.jcp.2010.02.001

107.V. Giovangigli, Mass conservation and singular multicomponent diffusion algorithms, *Impact Comput. Sci. Eng.*, **2** (1990), 73–97. https://doi.org/10.1016/0899-8248(90)90004-T

108.V. Giovangigli, Convergent iterative methods for multicomponent diffusion, *Impact Comput. Sci. Eng.*, **3** (1991), 244–276. https://doi.org/10.1016/0899-8248(91)90010-R

109.A. Ern, V. Giovangigli, *Multicomponent transport algorithms*, Lecture Notes in Physics, Monograph m24, Heidelberg: Springer-Verlag, 1994. https://doi.org/10.1007/978-3-540-48650-3

110.A. Ern, V. Giovangigli, The structure of transport linear systems in Dilute isotropic gas mixtures, *Phys. Rev. E*, **53** (1996), 485–492. https://doi.org/10.1103/PhysRevE.53.485

111.A. Ern, V. Giovangigli, Thermal diffusion effects in hydrogen-air and methane-air flames, *Combust. Theor. Model.*, **2** (1998), 349–372. https://doi.org/10.1088/1364-7830/2/4/001

112.P. H. Van Konynenburg, R. L. Scott, Critical lines and phase equilibria in binary van der Waals mixtures, *Philos. T. R. Soc.*, **298** (1980), 495–540. https://doi.org/10.1098/rsta.1980.0266

113.J. M. Prausnitz, R. N. Lichtenthaler, E. G. de Avezo, *Molecular thermodynamics of fluid-phase equilibria*, Pearson Education, 1999.

114.J. F. Ely, H. J. Hanley, Prediction of transport properties. 2. Thermal conductivity of pure fluids and mixtures, *Indus. Eng. Chem. Fundat.*, **22** (1983), 90–97. https://doi.org/10.1021/i100009a016

115.T. H. Chung, M. Ajlan, L. L. Lee, K. E. Starling, Generalized multiparameter correlation for nonpolar and polar fluid transport properties, *Ind. Eng. Chem. Res.*, **27** (1988), 671–679. https://doi.org/10.1021/ie00076a024

116.V. I. Kurochkin, S. F. Makarenko, G. A. Tirskii, Transport coefficients and the Onsager relations in the kinetic theory of dense gas mixtures, *J. Appl. Mech. Tech. Ph.*, **25** (1984), 218–225. https://doi.org/10.1007/BF00910464

117.R. J. Kee, F. M. Rupley, J. A. Miller, *The Chemkin thermodynamic data base*, SANDIA National Laboratories Report, 1987. https://doi.org/10.2172/7073290

118.M. W. Chase Jr, NIST-JANAF thermochemical tables, 4 Eds., *J. Phys. Chem. Ref. Data*, 1998.

119.B. J. McBride, M. J. Zehe, S. Gordon, *NASA Glenn coefficients for calculating thermodynamic properties of individual species*, National Aeronautics and Space Administration, John H. Glenn Research Center at Lewis Field, 2002.

120.H. Lin, Y. Y. Duan, Q. Min, Gradient theory modeling of surface tension for pure fluids and binary mixtures, *Fluid Phase Equilibr.*, **254** (2007), 75–90. https://doi.org/10.1016/j.fluid.2007.02.013

121.B. A. Younglove, Thermophysical properties of fluids I, Argon, ethylene, parahydrogen, nitrogen, nitrogen trifluoride, and oxygen, *J. Phys. Chem. Ref. Data*, **11** (1982).
122.B. A. Younglove, J. F. Ely, Thermophysical properties of fluids. II. Methane, ethane, propane, isobutane, and normal butane, *J. Phys. Chem. Ref. Data*, **16** (1987), 577–798.

123.M. K. Gupta, G. C. Gardner, M. J. Hegarty, A. J. Kidnay, Liquid-vapor equilibria for the $N_2^+CH_4^+C_2H_6$ System from 260 to 280 K, *J. Chem. Eng. Data*, **25** (1980), 313–318. https://doi.org/10.1021/je60087a016

124.B. E. Eakin, R. T. Ellington, D. C. Gami, *Physical-chemical properties of Ethane-Nitrogen mixtures*, Institute of Gas Technology, 1955.

125.R. Stryjek, P. S. Chappelear, R. Kobayashi, Low-temperature vapor-liquid equilibria of Nitrogen-Ethane system, *J. Chem. Eng. Data*, **19** (1974), 340–343. https://doi.org/10.1021/je60063a024

126.K. D. Wisotzki, G. M. Schneider, Fluid phase equilibria of the binary systems $N_2^+$ Ethane and $N_2^+$ Pentane Between 88 K and 313 K and Pressures up to 200 MPa, *Ber. Bunsenges. Phys. Chem.*, **89** (1985), 21–25. https://doi.org/10.1002/bbpc.19850890106

127.M. L. Japas, E. U. Franck, High pressure phase equilibria and PVT-data of the water-oxygen system including water-air to 673 K and 250 MPa, *Ber. Bunsenges. Phys. Chem.*, **89** (1985), 1268–1275. https://doi.org/10.1002/bbpc.19850891206

128.D. Y. Peng, D. B. Robinson, A rigorous method for predicting the critical properties of multicomponent systems from an equation of state, *AICHE J.*, **23** (1977), 137–144. https://doi.org/10.1002/aic.690230202

129.A. Mielke, A gradient structure for reaction-diffusion systems and for energy-drift-diffusion systems, *Nonlinearity*, **24** (2011), 1329–1346. https://doi.org/10.1088/0951-7715/24/4/016

130.M. Liero, A. Mielke, Gradient structures and geodesic convexity for reaction-diffusion systems, *Philos. T. R. Soc.*, **A371** (2013), 20120346. https://doi.org/10.1098/rsta.2012.0346

131.L. Dong, C. Wang, S. M. Wise, Z. Zhang, A positivity-preserving, energy stable scheme for a ternary Cahn-Hilliard system with the singular interfacial parameters, *J. Comput. Phys.*, **442** (2021), 110451. https://doi.org/10.1016/j.jcp.2021.110451

# 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 $\mathcal{E} = \{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 $y_1$. Since the mass fractions are assumed to be independent, the mass fraction of the single species $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, y_1)^t$, the thermodynamic is defined with $e(\zeta)$, $p(\zeta)$, and $s(\zeta)$, and Gibbs relation is in the form

$$ Tds = de + pdv - g_1dy_1. \qquad (A.1) $$

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

$$ v = y_1v_1. \qquad (A.2) $$
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, y_1) = y_1 e(T, v_1, 1)$, $p(T, v, y_1) = p(T, v_1, 1)$, and $s(T, v, y_1) = 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

$$e(\zeta) = y_1 e_1(T, v_1), \qquad (A.3)$$

$$p(\zeta) = p_1(T, v_1), \qquad (A.4)$$

$$s(\zeta) = y_1 s_1(T, v_1). \qquad (A.5)$$

From the definition $g_1 = \partial_{y_1} e - T \partial_{y_1} s$, we also deduce that

$$g_1 = \partial_{y_1} \left( y_1 e(T, v/y_1, 1) - T y_1 s(T, v/y_1, 1) \right) = e_1 - T s_1 + \frac{y_1 v}{y_1^2} \partial_v \left( e(T, v_1, 1) - T y_1 s(T, v_1, 1) \right),$$

so that

$$g_1 = e_1 - T s_1 + v_1 p_1. \qquad (A.6)$$

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 $Tds = \partial_T e \,dT + (p + \partial_v e) \,dv + (\partial_{y_1} e - g_1) \,dy_1$, we now deduce, after some algebra, that

$$y_1(Tds_1 - de_1 + pdv_1) + (Ts_1 - e_1 - v_1p_1 + g_1)dy_1 = 0,$$

and this yields the usual relation

$$Tds_1 = de_1 + pdv_1, \qquad (A.7)$$

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, y_1)$, but the link is easily made.

## B. Proof of Proposition 2.6

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

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

$$\tilde{\partial}_e \chi = \partial_T \chi \tilde{\partial}_e T, \qquad \partial_T \chi = \tilde{\partial}_e \chi \partial_T e, \qquad (B.1)$$

$$\tilde{\partial}_v \chi = \partial_v \chi + \partial_T \chi \tilde{\partial}_v T, \qquad \partial_v \chi = \tilde{\partial}_v \chi + \tilde{\partial}_e \chi \partial_v e, \qquad (B.2)$$

$$\tilde{\partial}_{y_k} \chi = \partial_{y_k} \chi + \partial_T \chi \tilde{\partial}_{y_k} T, \qquad \partial_{y_k} \chi = \tilde{\partial}_{y_k} \chi + \tilde{\partial}_e \chi \partial_{y_k} e, \qquad (B.3)$$

and letting $\chi = T$ in (B.2), we deduce that $\tilde{\partial}_v T = -(\tilde{\partial}_e T)\partial_v e$ and, finally, $T^2 \tilde{\partial}_{ee}^2 s = (\tilde{\partial}_e T)\partial_v e$. Similarly, we have $T^2 \tilde{\partial}_{ey_k}^2 s = -\tilde{\partial}_{y_k} T$, and using $\chi = T$ in (B.3), we deduce $\tilde{\partial}_{y_k} T = -(\tilde{\partial}_e T)\partial_{y_k} e$ so that $T^2 \tilde{\partial}_{ey_k}^2 s = (\tilde{\partial}_e T)\partial_{y_k} e$. Furthermore, we have $\tilde{\partial}_{vv}^2 s = \tilde{\partial}_v(\frac{p}{T})$, and using (B.2), we deduce that $\tilde{\partial}_v(\frac{p}{T}) = \partial_v(\frac{p}{T}) +$
$\partial_T(\frac{p}{T})\tilde{\partial}_v T$. From the compatibility conditions (2.8) and $\tilde{\partial}_v T = -(\tilde{\partial}_e T)\partial_v e$, we thus obtain $\tilde{\partial}_{vv}^2 s = \partial_v(\frac{p}{T}) - (\tilde{\partial}_e T)\frac{(\partial_v e)^2}{T^2}$. Similarly, $\tilde{\partial}_{vy_k}^2 s = \tilde{\partial}_{y_k}(\frac{p}{T})$, and thanks to (B.3), we have $\tilde{\partial}_{y_k}(\frac{p}{T}) = \partial_{y_k}(\frac{p}{T}) + \partial_T(\frac{p}{T})\tilde{\partial}_{y_k} T$. From the compatibility conditions (2.8) and $\tilde{\partial}_{y_k} T = -(\tilde{\partial}_e T)\partial_{y_k} e$, we get

$$\tilde{\partial}_{vy_k}^2 s = \partial_v(\frac{-g_k}{T}) - (\tilde{\partial}_e T)\frac{(\partial_v e)(\partial_{y_k} e)}{T^2}.$$

Since $\tilde{\partial}_{y_k y_l}^2 s = \tilde{\partial}_{y_l}(\frac{-g_k}{T})$, and using $\chi = \frac{-g_k}{T}$ in (B.3), we also obtain that

$$\tilde{\partial}_{y_k y_l}^2 s = \partial_{y_l}(\frac{-g_k}{T}) - (\tilde{\partial}_e T)\frac{(\partial_{y_k} e)(\partial_{y_l} e)}{T^2}.$$

Defining $\tilde{f}^e = (1, -\partial_v e, -\partial_{y_1} e, \dots, -\partial_{y_n} e)^t$, we have established that for any $x = (x_e, x_v, x_1, \dots, x_n)^t \in \mathbb{R}^{2+n}$

$$\langle (\tilde{\partial}_{\xi\xi}^2 s)x, x \rangle = -\frac{\tilde{\partial}_e T}{T^2} \langle \tilde{f}^e, x \rangle^2 + \frac{\partial_v p}{T} x_v^2 - 2 \sum_{k \in \Xi} \partial_v(\frac{g_k}{T}) x_v x_k - \sum_{k,l \in \Xi} \partial_{y_l}(\frac{g_k}{T}) x_k x_l. \quad (B.4)$$

We now prove that (i) implies (ii). Since $\xi = e f^e + v f^v + \sum_{i \in \Xi} y_i f^{vi}$ and $f^e$ are not proportional we deduce that $\langle (\tilde{\partial}_{\xi\xi}^2 s) f^e, f^e \rangle < 0$ so that $\tilde{\partial}_e T = -T^2 \langle (\tilde{\partial}_{\xi\xi}^2 s) f^e, f^e \rangle > 0$, and this yields $\partial_T e > 0$ since $\tilde{\partial}_e T \partial_T e = 1$. Letting $f = f^v + (\partial_v e) f^e$, we have $\langle \tilde{f}^e, f \rangle = 0$ and $f = f^v + (\partial_v e) f^e$, and since the mass fractions are positive, $\xi$ cannot be proportional. Thus, $\partial_v p/T = \langle (\tilde{\partial}_{\xi\xi}^2 s) f, f \rangle < 0$, and we have proved that $\partial_v p < 0$. In the same vein, assume that $x_1, \dots, x_n$ are arbitrary with $(x_1, \dots, x_n) \neq (0, \dots, 0)$, consider $f = \sum_{i \in \Xi} x_i f^i + x_e f^e$, and select $x_e$ such that $\langle \tilde{f}^e, f \rangle = 0$, that is, $x_e = \sum_{i \in \Xi} x_i \partial_{y_i} e$. Then, $f$ and $\xi$ cannot be proportional since the volume per unit mass is positive, so that $\langle (\tilde{\partial}_{\xi\xi}^2 s) f, f \rangle < 0$. This yields that $\sum_{k,l \in \Xi} \partial_{y_l}(\frac{g_k}{T}) x_k x_l = -\langle (\tilde{\partial}_{\xi\xi}^2 s) f, f \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 $\tilde{f}^v = (0, 1, -\frac{\partial_v g_1}{\partial_v p}, \dots, -\frac{\partial_v g_n}{\partial_v p})^t$ and using (B.4), we easily deduce that for $x = (x_e, x_v, x_1, \dots, x_n)^t$

$$\langle (\tilde{\partial}_{\xi\xi}^2 s)x, x \rangle = -\tilde{\partial}_e T \frac{\langle \tilde{f}^e, x \rangle^2}{T^2} + \partial_v p \frac{\langle \tilde{f}^v, x \rangle^2}{T} - \sum_{k,l \in \Xi} \left( \partial_{y_l} \left( \frac{g_k}{T} \right) + \frac{\partial_v g_k \partial_v g_l}{T \partial_v p} \right) x_k x_l. \quad (B.5)$$

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

$$\langle (\tilde{\partial}_{\xi\xi}^2 s)x, x \rangle = -\tilde{\partial}_e T \frac{\langle \tilde{f}^e, x \rangle^2}{T^2} + \partial_v p \frac{\langle \tilde{f}^v, x \rangle^2}{T} - \left\langle \left( \Lambda - \frac{\Lambda y \otimes \Lambda y}{\langle \Lambda y, y \rangle} \right) x_y, x_y \right\rangle, \quad (B.6)$$

where $x_y = (x_1, \dots, x_n)^t$. Since $\Lambda$ is positive definite, we first note that $\Lambda - \frac{\Lambda y \otimes \Lambda y}{\langle \Lambda y, y \rangle}$ is positive semi-definite with null space $\mathbb{R}y$ since for any $x_y = (x_1, \dots, x_n)^t$, we have

$$\left\langle \left( \Lambda - \frac{\Lambda y \otimes \Lambda y}{\langle \Lambda y, y \rangle} \right) x_y, x_y \right\rangle = \left\langle \Lambda(x_y - \frac{\langle \Lambda x_y, y \rangle}{\langle \Lambda y, y \rangle} y), (x_y - \frac{\langle \Lambda x_y, y \rangle}{\langle \Lambda y, y \rangle} y) \right\rangle. \quad (B.7)$$

Since $\partial_T e > 0$, $\partial_v p < 0$, and $\Lambda - \frac{\Lambda y \otimes \Lambda y}{\langle \Lambda y, y \rangle}$ is positive semi-definite with null space $\mathbb{R}y$, it is obtained from (B.6) that $\tilde{\partial}_{\xi\xi}^2 s$ is negative semi-definite. Moreover, any vector $x$ in the null space of $\tilde{\partial}_{\xi\xi}^2 s$ is
such that $\langle \hat{f}^e, x \rangle = 0$, $\langle \hat{f}^v, x \rangle = 0$, $(\Lambda - \frac{\Lambda y \otimes \Lambda y}{\langle \Lambda y, y \rangle})x_y = 0$. Since $x_y = (x_1, \dots, x_n)^t$ is in the null space of $\Lambda - \frac{\Lambda y \otimes \Lambda y}{\langle \Lambda y, y \rangle}$ we first deduce that $x_y = \lambda y$ for some $\lambda \in \mathbb{R}$. From the expression of $\hat{f}^v$, we first obtain that $x_v = \sum_{i \in \mathfrak{S}} x_i \partial_v g_i / \partial_v p$ so that $x_v = \lambda \sum_{i \in \mathfrak{S}} y_i \partial_v g_i / \partial_v p$. Using then the compatibility relations (2.8), we have $-\partial_v g_i = \partial_{y_i} p$ so that $x_v = -\lambda \sum_{i \in \mathfrak{S}} y_i \partial_{y_i} p / \partial_v p$ and, finally, $x_v = \lambda v$ because of the 0-homogeneity of $p$. In addition, from the expression of $\hat{f}^e$, we get that $x_e = x_v \partial_v e + \sum_{i \in \mathfrak{S}} x_i \partial_{y_i} e_i$ so that $x_e = \lambda v \partial_v e + \lambda \sum_{i \in \mathfrak{S}} y_i \partial_{y_i} e_i$, and $x_e = \lambda e$ because of the 1-homogeneity of $e$. We thus conclude that $x$ is proportional to $\xi$ so that $N(\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

$$ \langle \Lambda x_y, x_y \rangle = \frac{1}{\langle \Lambda y, y \rangle} \langle \Lambda x_y, y \rangle^2 + \langle \Lambda x_y, x_y \rangle, $$

where $x_y = (x_1, \dots, x_n)^t$. This identity may be obtained from (B.7) or by decomposing $x_y$ in the form $x_y = \alpha y + z$ with $\langle \Lambda y, z \rangle = 0$. $\square$

### 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 $\mathcal{P}$ and the heat flux $\mathcal{Q}$ are found in the form [44]

$$ \mathcal{P} = p^{\text{ex}} \mathbf{I} + \sum_{i,j \in \mathfrak{S}} \kappa_{ij} \nabla \rho_i \otimes \nabla \rho_j - \sum_{i,j \in \mathfrak{S}} \rho_i \nabla \cdot (\kappa_{ij} \nabla \rho_j) + \mathcal{P}^d, $$

$$ \mathcal{Q} = \sum_{i,j \in \mathfrak{S}} \kappa_{ij} \nabla \rho_j (\rho_i \nabla \cdot \mathbf{u} + \nabla \cdot \mathcal{J}_i) - \sum_{i,j \in \mathfrak{S}} \nabla \cdot (\kappa_{ij} \nabla \rho_j) \mathcal{J}_i + \mathbf{q}, $$

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

$$ \mathcal{P}^d = - \nu \nabla \cdot \mathbf{u} \mathbf{I} - \eta (\nabla \mathbf{u} + \nabla \mathbf{u}^t - \frac{2}{3} \nabla \cdot \mathbf{u} \mathbf{I}), $$

$$ \mathcal{J}_i = - \sum_{j \in \mathfrak{S}} L_{ij} \left( \nabla \left( \frac{g_j^{\text{ex}}}{T} \right) - \frac{\nabla \cdot \nabla \cdot \left( \sum_{l \in \mathfrak{S}} \kappa_{jl} \nabla \rho_l \right)}{T} \right) - L_{ie} \nabla \left( \frac{-1}{T} \right), $$

$$ \mathbf{q} = - \sum_{i \in \mathfrak{S}} L_{ei} \left( \nabla \left( \frac{g_i^{\text{ex}}}{T} \right) - \frac{\nabla \cdot \nabla \cdot \left( \sum_{l \in \mathfrak{S}} \kappa_{il} \nabla \rho_l \right)}{T} \right) - L_{ee} \nabla \left( \frac{-1}{T} \right), $$

where $\nu$ denotes the bulk viscosity, $\eta$ the shear viscosity, and $L = (L_{ij})_{i,j \in \mathfrak{S} \cup \{e\}}$ the mass and heat diffusion matrix. These coefficients $\nu$, $\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}^{\text{ex}} = \mathcal{F}(\rho_1, \dots, \rho_n, T) + \frac{1}{2} \sum_{i,j \in \mathcal{Z}} \kappa_{ij} \nabla \rho_i \cdot \nabla \rho_j$ is of Van der Waals type, and the corresponding thermodynamic properties are found in the form

$$ \mathcal{E}^{\text{ex}} = \mathcal{E} + \frac{1}{2} \sum_{i,j \in \mathcal{Z}} (\kappa_{ij} - T \partial_T \kappa_{ij}) \nabla \rho_i \cdot \nabla \rho_j, \qquad p^{\text{ex}} = p - \frac{1}{2} \sum_{i,j \in \mathcal{Z}} \kappa_{ij} \nabla \rho_i \cdot \nabla \rho_j, $$

$$ S^{\text{ex}} = S - \frac{1}{2} \sum_{i,j \in \mathcal{Z}} \partial_T \kappa_{ij} \nabla \rho_i \cdot \nabla \rho_j, \qquad g_i^{\text{ex}} = g_i, \quad i \in \mathcal{Z}, $$

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

$$ T dS^{\text{ex}} = d\mathcal{E}^{\text{ex}} - \sum_{i \in \mathcal{Z}} g_i^{\text{ex}} d\rho_i - \sum_{i,j \in \mathcal{Z}} \kappa_{ij} \nabla \rho_i \cdot d\nabla \rho_j. \qquad (C.1) $$

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:

$$ \kappa_{ij} = \kappa(T) \quad i, j \in \mathcal{Z}. $$

The governing equations may then be simplified with the relations $\sum_{i,j \in \mathcal{Z}} \kappa_{ij} \nabla \rho_i \otimes \nabla \rho_j = \kappa \nabla \rho \otimes \nabla \rho$, $\sum_{i,j \in \mathcal{Z}} \kappa_{ij} \nabla \rho_i \cdot \nabla \rho_j = \rho \nabla \cdot (\kappa \nabla \rho)$, and $\sum_{i,j \in \mathcal{Z}} \kappa_{ij} \nabla \rho_i \nabla \rho_j \cdot \mathbf{u} = \rho \nabla \rho \nabla \cdot \mathbf{u}$, as well as the thermodynamic properties with $\sum_{i,j \in \mathcal{Z}} \kappa_{ij} \nabla \rho_i \cdot \nabla \rho_j = \kappa |\nabla \rho|^2$. Various terms may further be eliminated from the heat flux with $\sum_{i,j \in \mathcal{Z}} \kappa_{ij} \nabla \rho_j \nabla \cdot \mathcal{J}_i = 0$ and $\sum_{i,j \in \mathcal{Z}} \nabla \cdot (\kappa_{ij} \nabla \rho_j) \mathcal{J}_i = 0$ since $\sum_{i \in \mathcal{Z}} \mathcal{J}_i = 0$. The Cahn-Hilliard type driving forces are also species independent, since $\sum_{i \in \mathcal{Z}} \kappa_{il} \nabla \rho_l = \kappa \nabla \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

$$ \mathcal{P} = p^{\text{ex}} \mathbf{I} + \kappa \nabla \rho \otimes \nabla \rho - \rho \nabla \cdot (\kappa \nabla \rho) \mathbf{I} + \mathcal{P}^{\text{d}}, \qquad \mathcal{Q} = \kappa \rho \nabla \rho \nabla \cdot \mathbf{u} + \mathbf{q}, \qquad (C.2) $$

$$ \mathcal{J}_i = - \sum_{j \in \mathcal{Z}} L_{ij} \nabla \left( \frac{g_j}{T} \right) - L_{ie} \nabla \left( \frac{-1}{T} \right), \qquad \mathbf{q} = - \sum_{i \in \mathcal{Z}} L_{ei} \nabla \left( \frac{g_i}{T} \right) - L_{ee} \nabla \left( \frac{-1}{T} \right). \qquad (C.3) $$

The thermodynamic properties also read

$$ \mathcal{E}^{\text{ex}} = \mathcal{E}(\rho_1, \dots, \rho_n, T) + \frac{1}{2}(\kappa - T \partial_T \kappa) |\nabla \rho|^2, \qquad p^{\text{ex}} = p(\rho_1, \dots, \rho_n, T) - \frac{1}{2} \kappa |\nabla \rho|^2, \qquad (C.4) $$

$$ S^{\text{ex}} = S(\rho_1, \dots, \rho_n, T) - \frac{1}{2} \partial_T \kappa |\nabla \rho|^2, \qquad g_i^{\text{ex}} = g_i(\rho_1, \dots, \rho_n, T), \quad i \in \mathcal{Z}, \qquad (C.5) $$

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

<figure><img src="image_10.png" /></figure>

AIMS Press

©2024 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0)