Loading [MathJax]/jax/element/mml/optable/BasicLatin.js
Research article

An innovative extended Bayesian analysis of the relationship between returns and different risk measures in South Africa

  • This study investigated the All Share Index (ALSI) returns and six different risk measures of the South African market for the sample period from 17 March 2000 to 17 March 2022. The risk measures analyzed were standard deviation (SD), absolute deviation (AD), lower semi absolute deviation (LSAD), lower semivariance (LSV), realized variance (RV) and the bias-adjusted realized variance (ARV). This study made an innovative contribution on a methodological and practical level, by being the first study to extend from the novel Bayesian approach by Jensen and Maheu (2018) to methods by Karabatsos (2017)—density regression, quantile regression and survival analysis. The extensions provided a full representation of the return distribution in relation to risk, through graphical analysis, producing novel insight into the risk-return topic. The most novel and innovative contribution of this study was the application of survival analysis which analyzed the "life" and "death" of the risk-return relationship. From the density regression, this study found that the chance of investors earning a superior return was substantial and that the probability of excess returns increased over time. From quantile regression, results revealed that returns have a negative relationship with the majority of the risk measures—SD, AD, LSAD and RV. However, a positive risk-return relationship was found by LSV and the ARV, with the latter having the steepest slope. Results were the most pronounced for the ARV, especially for the survival analysis. While ARV earned the highest returns, it had the shortest lifespan, which can be attributed to the volatile nature of the South African market. Thus, investors that seek short-term high-earning returns would examine ARV followed by LSV, whereas the remaining risk measures can be used for other purposes, such as diversification purposes or short selling.

    Citation: Nitesha Dwarika. An innovative extended Bayesian analysis of the relationship between returns and different risk measures in South Africa[J]. Quantitative Finance and Economics, 2022, 6(4): 570-603. doi: 10.3934/QFE.2022025

    Related Papers:

    [1] Yiping Meng . On the Rayleigh-Taylor instability for the two coupled fluids. AIMS Mathematics, 2024, 9(11): 32849-32871. doi: 10.3934/math.20241572
    [2] Osama Moaaz, Ahmed E. Abouelregal . Multi-fractional-differential operators for a thermo-elastic magnetic response in an unbounded solid with a spherical hole via the DPL model. AIMS Mathematics, 2023, 8(3): 5588-5615. doi: 10.3934/math.2023282
    [3] Ze Wang, Yan Zhang, Jincheng Shi, Baiping Ouyang . Spatial decay estimates for the Fochheimer equations interfacing with a Darcy equations. AIMS Mathematics, 2021, 6(11): 12632-12649. doi: 10.3934/math.2021728
    [4] Shu-Nan Li, Bing-Yang Cao . Anomalies of Lévy-based thermal transport from the Lévy-Fokker-Planck equation. AIMS Mathematics, 2021, 6(7): 6868-6881. doi: 10.3934/math.2021402
    [5] Ahmed G. Salem, Turki D. Alharbi, Abdulaziz H. Alharbi, Anwar Ali Aldhafeeri . Impact of a spherical interface on a concentrical spherical droplet. AIMS Mathematics, 2024, 9(10): 28400-28420. doi: 10.3934/math.20241378
    [6] Madeeha Tahir, Ayesha Naz, Muhammad Imran, Hasan Waqas, Ali Akgül, Hussein Shanak, Rabab Jarrar, Jihad Asad . Activation energy impact on unsteady Bio-convection nanomaterial flow over porous surface. AIMS Mathematics, 2022, 7(11): 19822-19845. doi: 10.3934/math.20221086
    [7] Xi Wang . Blow-up criterion for the 3-D inhomogeneous incompressible viscoelastic rate-type fluids with stress-diffusion. AIMS Mathematics, 2025, 10(1): 1826-1841. doi: 10.3934/math.2025084
    [8] Dario Bambusi, Simone Paleari . A couple of BO equations as a normal form for the interface problem. AIMS Mathematics, 2024, 9(8): 23012-23026. doi: 10.3934/math.20241118
    [9] Maxime Krier, Julia Orlik . Solvability of a fluid-structure interaction problem with semigroup theory. AIMS Mathematics, 2023, 8(12): 29490-29516. doi: 10.3934/math.20231510
    [10] Ahmed E. Abouelregal, Faisal Alsharif, Hashem Althagafi, Yazeed Alhassan . Fractional heat transfer DPL model incorporating an exponential Rabotnov kernel to study an infinite solid with a spherical cavity. AIMS Mathematics, 2024, 9(7): 18374-18402. doi: 10.3934/math.2024896
  • This study investigated the All Share Index (ALSI) returns and six different risk measures of the South African market for the sample period from 17 March 2000 to 17 March 2022. The risk measures analyzed were standard deviation (SD), absolute deviation (AD), lower semi absolute deviation (LSAD), lower semivariance (LSV), realized variance (RV) and the bias-adjusted realized variance (ARV). This study made an innovative contribution on a methodological and practical level, by being the first study to extend from the novel Bayesian approach by Jensen and Maheu (2018) to methods by Karabatsos (2017)—density regression, quantile regression and survival analysis. The extensions provided a full representation of the return distribution in relation to risk, through graphical analysis, producing novel insight into the risk-return topic. The most novel and innovative contribution of this study was the application of survival analysis which analyzed the "life" and "death" of the risk-return relationship. From the density regression, this study found that the chance of investors earning a superior return was substantial and that the probability of excess returns increased over time. From quantile regression, results revealed that returns have a negative relationship with the majority of the risk measures—SD, AD, LSAD and RV. However, a positive risk-return relationship was found by LSV and the ARV, with the latter having the steepest slope. Results were the most pronounced for the ARV, especially for the survival analysis. While ARV earned the highest returns, it had the shortest lifespan, which can be attributed to the volatile nature of the South African market. Thus, investors that seek short-term high-earning returns would examine ARV followed by LSV, whereas the remaining risk measures can be used for other purposes, such as diversification purposes or short selling.



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

    We initially discuss the mathematical structure and properties of classical thermodynamic functions in terms of intensive variables indispensable for fluid models. The complex situations of polymers [52], multitemperature flows [53,54], and links with statistical mechanics investigated by Fowler [55], Ferziger and Kaper [56], Keizer [57], Laasonen et al. [58], and Chen and Doolen [59], lay beyond the scope of this work. General physical formalisms including conservation laws and spatial transport like the thermodynamics of irreversible processes [60] or the GENERIC formalism of H. C. Öttinger [61] also lay beyond the scope of this work. A mathematical definition of thermodynamics is introduced and involves smoothness and homogeneity properties of thermodynamic functions, Gibbs' relation, and the compatibility with ideal gases at low density. We characterize thermal, mechanical and chemical thermodynamic stability. Thermal stability is generally guaranteed in practical applications and instabilities may then only be of mechanical or chemical type and determined by the spectrum of a single symmetric matrix Λ related to the species chemical potentials or Gibbs' functions [9,29,30]. The natural variable is denoted by ζ=(T,v,y1,,yn)t where T is the absolute temperature, v the volume per unit mass, y1,,yn the species mass fractions supposed to be independent, and n the number of species. The Gibbsian thermodynamic model is specified as e(ζ), p(ζ), and s(ζ), where e denotes the energy per unit mass, p the pressure, and s the entropy per unit mass, over an open set OζR2+n. We also address a different structure associated with the thermodynamic variable Z=(T,ρ1,,ρn)t, where ρ1,,ρn denotes the species masses per unit volume. Using the variable ζ=(T,v,y1,,yn)t leads to homogeneous thermodynamic functions—with singular entropy Hessians—and using the variable Z=(T,ρ1,,ρn)t leads to nonhomogeneous thermodynamic functions and nonsingular entropy Hessian matrices when stability holds. These two frameworks are useful and have been investigated in the literature generally for ideal gas mixtures [1,2,3,4,5,6,7,8,9]. We establish the mathematical equivalence between these thermodynamic formalisms under proper transformation rules. The thermodynamic functions are not defined for all states because of nonidealities and thermodynamic instabilities—associated with phase changes—are allowed. We also discuss the thermodynamic—part of the—entropic 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(ζ) with respect to (v,y1,,yn) and a quadratic expression of ppid as ρ goes to zero where pid denotes the ideal gas pressure. We determine the energy and entropy per unit mass and then establish that they constitute a thermodynamics according to the previous definition. We also also introduce rescaled versions of the stability matrix Λ that behave smoothly for vanishing mass fractions.

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

    We next investigate the mathematical structure of an extended thermodynamics associated with a simplified diffuse interface fluid model [20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51]. Diffuse interface models, like those of van der Waals or Cahn and Hilliard, are obtained in the presence of internal cohesive forces and involve density gradient squared terms in the free energies. Such diffuse interface models describe interphases—interfaces between phases—as regions with smooth variations of physical properties and have been used successfully in order to describe spinodal decomposition, droplets dynamics, three phase contact lines, and surface diffusion, as well as transcritical flames [20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51]. The simplified multicomponent flow model investigated in this paper is obtained by assuming that the capillarity coefficients are identical and derived in appendix C. In this situation, only the density gradient ρ 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,86,87] as well as for numerical simulations of diffuse interface fluids [84]. We also address the mathematical structure of chemical production terms for completeness even though chemical reactions will not be considered in the numerical simulations. The chemical production rates involve, in particular, Marcelin's expression for the chemical reaction rates of progress discussed by Marcelin [88,89], Keizer [57], and Giovangigli et al. [8,80,90,91].

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

    The thermodynamic formalism is validated by comparison with experimental measurements for mixtures of ethane C2H6 and nitrogen N2. 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 C2H6 and nitrogen N2 may be chemically unstable. These mixtures may split between an ethane-rich gaseous-like phase and an ethane-poor liquid-like phase. The two phase chemical equilibrium curves are found to be in very good agreement with experimental results of Gupta et al. [123]. We also investigate the critical points of ethane-nitrogen mixtures. Following the classification of Van Konynenburg and Scott [112,113], the curve of critical points is known to be of Type Ⅲ [124,125,126]. The numerical simulations confirm such a type Ⅲ behavior with an overall very satisfactory agreement in comparison with experimental results of Eakin at al. [124], Stryjek et al. [125], and Wisotzki and Schneider [126]. We finally investigate the structure of ethane liquid-vapor interfaces and study the influence of the capillarity coefficients on the flow structure.

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

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

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

    We investigate thermodynamics in terms of the natural variable ζ=(T,v,y1,,yn)t where T denotes the absolute temperature, v the volume per unit mass, i.e., v=1/ρ where ρ is the mass density, y1,,yn the species mass fractions, n1 the number of species, and S the species indexing set S={1,,n}. The thermodynamic model is specified as e(ζ), p(ζ), and s(ζ), defined on an open set Oζ(0,)2+n, where e denotes the energy per unit mass, p the pressure, and s the entropy per unit mass. The open set Oζ frequently differs from (0,)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 the derivation operator with respect to the variable ζ=(T,v,y1,,yn)t, by d the total differential operator, and for any λ>0, we define ζλ=(T,λv,λy1,,λyn)t. We denote by y the mass fraction vector y=(y1,,yn)t, 1I the vector with unit components 1I=(1,,1)t, , the Euclidean scalar product, and define  Σ ={y(0,)n; y,1I=1}. Physically, the mass fraction y must sum up to unity y Σ . However, as usual in the modeling of multicomponent flows, we assume that the mass fractions are independent and the mass constraint y,1I=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 y1 is usually formally eliminated as discussed in Appendix A. We denote by γN, γ2, the regularity class of thermodynamics functions.

    Definition 2.1. Let e, p, and s, be functions of the variable ζ defined over an open set Oζ(0,)2+n. These functions are said to define a thermodynamics if properties (T0T2) hold.

    (T0) 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(ζ).

    (T1) For any ζOζ, defining gk=ykeTyks, kS, we have Gibbs' relation:

    Tds=TedT+(p+ve)dv+kS(ykegk)dyk. (2.1)

    (T2) For any (T,y1,yn)t(0,)1+n, there exists vm so that v>vm implies (T,v,y1,yn)tOζ and

    lim (2.2)

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

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

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

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

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

    \begin{equation} T \mathrm{d} s = \mathrm{d} e + p \mathrm{d} v - \sum\limits_{k\in {\mathfrak S}} g_k \mathrm{d} \mathrm{y}_k, \end{equation} (2.3)

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

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

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

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

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

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

    e^\sharp - e^\flat = \displaystyle{\int}_{T^\flat}^{T^\sharp} \!\!\! \partial^{}_T e(\theta, v, \mathrm{y}_1,\ldots, \mathrm{y}_ n) \, d\theta = 0.

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

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

    \begin{equation} \Lambda_{kl} = \frac{ \partial^{}_{ \mathrm{y}_k} g_l}{T} = \frac{ \partial^{}_{ \mathrm{y}_l} g_k}{T}, \qquad k,l\in {\mathfrak S}, \end{equation} (2.4)

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

    \begin{equation} \widehat\Lambda = \Lambda - \frac{\Lambda \mathrm{y} {\otimes}\Lambda \mathrm{y}}{ \langle \Lambda \mathrm{y}, \mathrm{y}\rangle}. \end{equation} (2.5)

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

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

    \begin{equation} \Lambda \mathrm{y} = \frac{ v}{T} ( \partial^{}_{ \mathrm{y}_1} p,\ldots, \partial^{}_{ \mathrm{y}_ n} p)^t, \end{equation} (2.6)
    \begin{equation} \langle \Lambda \mathrm{y}, \mathrm{y}\rangle = - \frac{ v^2}{T} \partial^{}_ v p. \end{equation} (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^{}_{ \mathrm{y}_k} s = \partial^{}_{ \mathrm{y}_k} e - g_k , and this implies the compatibility relations

    \begin{equation} \partial^{}_T \bigl(\frac{p}{T}\bigr) = \frac{ \partial^{}_ v e}{T^2}, \qquad \partial^{}_T \bigl(\frac{- g_k}{T}\bigr) = \frac{ \partial^{}_{ \mathrm{y}_k} e}{T^2}, \qquad - \partial^{}_ v g_k = \partial^{}_{ \mathrm{y}_k} p. \end{equation} (2.8)

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

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

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

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

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

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

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

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

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

    \begin{equation} \mathrm{d} g = - s \mathrm{d} T + v \mathrm{d} p + \sum\limits_{k\in {\mathfrak S}} g_k \mathrm{d} \mathrm{y}_k. \end{equation} (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 :

    \begin{equation} \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, \end{equation} (2.10)
    \begin{equation} \widehat\partial^{}_ p \chi = \partial^{}_ v \chi \; \widehat\partial^{}_ p v, \qquad \partial^{}_ v \chi = \widehat\partial^{}_ p \chi \; \partial^{}_ v p, \end{equation} (2.11)
    \begin{equation} \widehat\partial^{}_{ \mathrm{y}_k}\! \chi = \partial^{}_{ \mathrm{y}_k}\! \chi + \partial^{}_ v \chi \; \widehat\partial^{}_{ \mathrm{y}_k}\! v, \qquad \partial^{}_{ \mathrm{y}_k}\! \chi = \widehat\partial^{}_{ \mathrm{y}_k}\! \chi + \widehat\partial^{}_ p \chi \partial^{}_{ \mathrm{y}_k} p. \end{equation} (2.12)

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

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

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

    \begin{equation} \widehat\partial^2_{ \mathrm{y} \mathrm{y}} g = T \Bigl( \Lambda - \frac{\Lambda \mathrm{y}{\otimes}\Lambda \mathrm{y}}{ \langle\Lambda \mathrm{y}, \mathrm{y}\rangle} \Bigr). \end{equation} (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^{}_{ \mathrm{y}_k} g = g_k , and this implies the compatibility relations

    \begin{equation} \widehat\partial^{}_T v = - \widehat\partial^{}_ p s, \qquad \widehat\partial^{}_T g_k = - \widehat\partial^{}_{ \mathrm{y}_k} s, \qquad \widehat\partial^{}_ p g_k = \widehat\partial^{}_{ \mathrm{y}_k} v. \end{equation} (2.14)

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

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

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

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

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

    \begin{equation} \mathsf{X}_i = \frac{ \mathrm{y}_im}{m_i}, \quad i\in {\mathfrak S}, \qquad \frac{\sum\nolimits_{i\in {\mathfrak S}} \mathrm{y}_i}{m} = \sum\limits_{i\in {\mathfrak S}}\frac{ \mathrm{y}_i}{m_i}, \end{equation} (2.15)

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

    \begin{equation} {{\rm{G}}} = m g, \qquad\quad {{\rm{G}}}_i = m_i g_i, \qquad i\in {\mathfrak S}. \end{equation} (2.16)

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

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

    \begin{equation} p^ {\rm{id}} = \frac{ R T}{ v} \sum\limits_{k\in {\mathfrak S}} \frac{ \mathrm{y}_k}{m_k}, \end{equation} (2.17)

    and the internal energy by

    \begin{equation} e^ {\rm{id}} = \sum\limits_{k\in {\mathfrak S}} \mathrm{y}_k^{} e_k^ {\rm{id}}, \qquad e_k^ {\rm{id}} = e^ \mathrm{st}_k + \displaystyle{\int}_{T^ \mathrm{st}}^T \!\! c_{ {{{\rm{v}}}} k}^ {\rm{id}}(\theta)\,d\theta, \end{equation} (2.18)

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

    \begin{equation} s^ {\rm{id}} = \sum\limits_{k\in {\mathfrak S}} \mathrm{y}_k^{} s_k^ {\rm{id}}, \qquad s_k^ {\rm{id}} = s^ \mathrm{st}_k + \displaystyle{\int}_{T^ \mathrm{st}}^T \!\! \frac{ c_{ {{{\rm{v}}}} k}^ {\rm{id}}(\theta)}{\theta} \,d\theta - \frac{ R}{m_k} \log \frac{ \mathrm{y}_k}{ v m_k\gamma^ \mathrm{st}}, \end{equation} (2.19)

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

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

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

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

    \mathcal{O}_ \xi^ {{\rm{id}}} = \{ \xi = ( \xi_ e, \xi_ v, \xi_1, \ldots, \xi_ n)^t \ {with}\ \xi_ v > 0, \, \xi_1 > 0, \ldots, \xi_ n > 0, \, \xi_ e > \sum\limits_{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^ \mathrm{st} - \int_0^{T^ \mathrm{st}} \!\! c_{ {{{\rm{v}}}} k}^ {\rm{id}}(\theta)\, d\theta , k\in {\mathfrak S} .

    Proof. The proof is straightforward [8].

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

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

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

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

    \begin{equation} \mathcal{E}(T, \rho_1,\ldots, \rho_ n) = (\sum\limits_{i\in {\mathfrak S}} \rho_i) \; e \Bigl( T, \frac{1}{\sum\nolimits_{i\in {\mathfrak S}} \rho_i}, \frac{ \rho_1}{\sum\nolimits_{i\in {\mathfrak S}} \rho_i}, \ldots, \frac{ \rho_ n}{\sum\nolimits_{i\in {\mathfrak S}} \rho_i} \Bigr), \end{equation} (3.1)
    \begin{equation} \mathcal{P}(T, \rho_1,\ldots, \rho_ n) = p \Bigl( T, \frac{1}{\sum\nolimits_{i\in {\mathfrak S}} \rho_i}, \frac{ \rho_1}{\sum\nolimits_{i\in {\mathfrak S}} \rho_i}, \ldots, \frac{ \rho_ n}{\sum\nolimits_{i\in {\mathfrak S}} \rho_i} \Bigr), \end{equation} (3.2)
    \begin{equation} \mathcal{S}(T, \rho_1,\ldots, \rho_ n) = (\sum\limits_{i\in {\mathfrak S}} \rho_i) \; s \Bigl( T, \frac{1}{\sum\nolimits_{i\in {\mathfrak S}} \rho_i}, \frac{ \rho_1}{\sum\nolimits_{i\in {\mathfrak S}} \rho_i}, \ldots, \frac{ \rho_ n}{\sum\nolimits_{i\in {\mathfrak S}} \rho_i} \Bigr), \end{equation} (3.3)

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

    \begin{align} & \mathcal{O}_ \mathcal{Z} = \{\; (T, \rho_1,\ldots, \rho_ n) \in (0,\infty)^{1+ n}; \ \Bigl( T, \frac{1}{\sum\nolimits_{i\in {\mathfrak S}} \rho_i}, \frac{ \rho_1}{\sum\nolimits_{i\in {\mathfrak S}} \rho_i}, \ldots, \frac{ \rho_ n}{\sum\nolimits_{i\in {\mathfrak S}} \rho_i} \Bigr)\in \mathcal{O}_ \zeta \;\}. \end{align}

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

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

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

    \begin{equation} T \mathrm{d} \mathcal{S} = {\mathfrak d}^{}\mspace{-2mu}_T \mathcal{E} \mathrm{d} T + \sum\limits_{k\in {\mathfrak S}} ( {\mathfrak d}^{}\mspace{-2mu}_{ \rho_k} \mathcal{E} - g_k) \mathrm{d} \rho_k, \end{equation} (3.4)

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

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

    \begin{equation} \lim\limits_{ \rho\to0} \tfrac{1}{ \rho} ( \mathcal{E}- \mathcal{E}^ {\rm{id}}) = 0, \qquad \lim\limits_{ \rho\to0} \tfrac{1}{ \rho} ( \mathcal{P}- \mathcal{P}^ {\rm{id}}) = 0, \qquad \lim\limits_{ \rho\to0} \tfrac{1}{ \rho} ( \mathcal{S}- \mathcal{S}^ {\rm{id}}) = 0. \end{equation} (3.5)

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

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

    \begin{equation} \mathrm{d} \mathcal{E} = \rho\, \partial^{}_T e\, \mathrm{d} T + \sum\limits_{k\in {\mathfrak S}} \partial^{}_{ \mathrm{y}_k}\! e\, \mathrm{d} \rho_k, \end{equation} (3.6)

    and, similarly,

    \begin{equation} \mathrm{d} \mathcal{S} = \rho\, \partial^{}_T s\, \mathrm{d} T + \sum\limits_{k\in {\mathfrak S}} \partial^{}_{ \mathrm{y}_k}\! s\, \mathrm{d} \rho_k. \end{equation} (3.7)

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

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

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

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

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

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

    \begin{equation} T \mathrm{d} \mathcal{S} = \mathrm{d} \mathcal{E} - \sum\limits_{k\in {\mathfrak S}} g_k \mathrm{d} \rho_k. \end{equation} (3.8)

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

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

    \begin{equation} \Lambda_{kl} = \frac{ \rho {\mathfrak d}^{}\mspace{-2mu}_{ \rho_k} g_l}{T} = \frac{ \rho {\mathfrak d}^{}\mspace{-2mu}_{ \rho_l} g_k}{T}, \end{equation} (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}}_{\textsf{0}}-{\mathcal {T}}_{\textsf{2}}) are satisfied and that the map \mathcal{Z} \mapsto \mathcal{X} is a C^ \gamma local diffeomorphism from an open set \mathcal{O}_ \mathcal{Z}' onto an open set \mathcal{O}_ \mathcal{X}' . Then, for any \mathcal{Z}\in \mathcal{O}_ \mathcal{Z}' , the following statements are equivalent:

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

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

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

    \begin{equation} \widetilde{\mathfrak d}^{}\mspace{-2mu}_ \mathcal{E}^{} \chi = {\mathfrak d}^{}\mspace{-2mu}_T^{} \chi \; \widetilde{\mathfrak d}^{}\mspace{-2mu}_ \mathcal{E}^{} T, \qquad {\mathfrak d}^{}\mspace{-2mu}_T^{} \chi = \widetilde{\mathfrak d}^{}\mspace{-2mu}_ \mathcal{E}^{} \chi \; {\mathfrak d}^{}\mspace{-2mu}_T^{} \mathcal{E}, \end{equation} (3.10)
    \begin{equation} \widetilde{\mathfrak d}^{}\mspace{-2mu}_{ \rho_k}\! \chi = {\mathfrak d}^{}\mspace{-2mu}_{ \rho_k}\! \chi + {\mathfrak d}^{}\mspace{-2mu}_T^{} \chi \; \widetilde{\mathfrak d}^{}\mspace{-2mu}_{ \rho_k}\!T, \qquad {\mathfrak d}^{}\mspace{-2mu}_{ \rho_k}\! \chi = \widetilde{\mathfrak d}^{}\mspace{-2mu}_{ \rho_k}\! \chi + \widetilde{\mathfrak d}^{}\mspace{-2mu}_ \mathcal{E}^{} \chi \; {\mathfrak d}^{}\mspace{-2mu}_{ \rho_k}\! \mathcal{E}. \end{equation} (3.11)

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

    \begin{equation} \widetilde{\mathfrak d}^2_{ \mathcal{E} \mathcal{E}} \mathcal{S} = - \frac{ \widetilde{\mathfrak d}^{}\mspace{-2mu}_ \mathcal{E}^{} T}{T^2}, \qquad \widetilde{\mathfrak d}^2_{ \mathcal{E} \rho_k} \mathcal{S} = \frac{ \widetilde{\mathfrak d}^{}\mspace{-2mu}_ \mathcal{E} T\, \partial^{}_{ \mathrm{y}_k} e }{T^2}, \quad k\in {\mathfrak S}, \end{equation} (3.12)
    \begin{equation} \widetilde{\mathfrak d}^2_{ \rho_k \rho_l} \mathcal{S} = - \frac{ \partial^{}_{ \mathrm{y}_k} g_l}{ \rho T} - \frac{ \widetilde{\mathfrak d}^{}\mspace{-2mu}_ \mathcal{E}^{} T\, \partial^{}_{ \mathrm{y}_k} e \; \partial^{}_{ \mathrm{y}_l} e}{T^2}, \quad k,l\in {\mathfrak S}. \end{equation} (3.13)

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

    \begin{equation} \bigl\langle ( \widetilde{\mathfrak d}^2_{ \mathcal{X} \mathcal{X}} \! \mathcal{S})\, \mathsf{X} , \mathsf{X} \bigr\rangle = - \widetilde{\mathfrak d}^{}\mspace{-2mu}_ \mathcal{E} T \frac{\langle \mathfrak{f}^ \mathcal{E}, \mathsf{X} \rangle^2}{T^2} - \frac{1}{ \rho} \bigl\langle \Lambda \mathsf{X}_ \mathrm{y}^{}, \mathsf{X}_ \mathrm{y}^{} \bigr\rangle, \end{equation} (3.14)

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

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

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

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

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

    \begin{equation} e (T, v, \mathrm{y}_1,\ldots, \mathrm{y}_ n) = v \; \mathcal{E} \Bigl( T, \frac{ \mathrm{y}_1}{ v}, \ldots, \frac{ \mathrm{y}_ n}{ v} \Bigr), \end{equation} (3.15)
    \begin{equation} p (T, v, \mathrm{y}_1,\ldots, \mathrm{y}_ n) = \mathcal{P} \Bigl( T, \frac{ \mathrm{y}_1}{ v}, \ldots, \frac{ \mathrm{y}_ n}{ v} \Bigr), \end{equation} (3.16)
    \begin{equation} s (T, v, \mathrm{y}_1,\ldots, \mathrm{y}_ n) = v \; \mathcal{S} \Bigl( T, \frac{ \mathrm{y}_1}{ v}, \ldots, \frac{ \mathrm{y}_ n}{ v} \Bigr), \end{equation} (3.17)

    and defined over the open set

    \begin{equation} \mathcal{O}_ \zeta = \{\; (T, v, \mathrm{y}_1,\ldots, \mathrm{y}_ n) \in (0,\infty)^{2+ n}; \ \Bigl( T, \frac{ \mathrm{y}_1}{ v}, \ldots, \frac{ \mathrm{y}_ n}{ v} \Bigr)\in \mathcal{O}_ \mathcal{Z} \;\}. \end{equation} (3.18)

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

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

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

    \mathrm{d} s = v {\mathfrak d}^{}\mspace{-2mu}_T \mathcal{S}\, \mathrm{d} T + \bigl( \mathcal{S}- \sum\limits_{k\in {\mathfrak S}} \rho_k {\mathfrak d}^{}\mspace{-2mu}_{ \rho_k} \mathcal{S}\bigr)\, \mathrm{d} v + \sum\limits_{k\in {\mathfrak S}} {\mathfrak d}^{}\mspace{-2mu}_{ \rho_k}\! \mathcal{S}\, \mathrm{d} \mathrm{y}_k,
    \mathrm{d} e = v {\mathfrak d}^{}\mspace{-2mu}_T \mathcal{E}\, \mathrm{d} T + \bigl( \mathcal{E}- \sum\limits_{k\in {\mathfrak S}} \rho_k {\mathfrak d}^{}\mspace{-2mu}_{ \rho_k} \mathcal{E}\bigr)\, \mathrm{d} v + \sum\limits_{k\in {\mathfrak S}} {\mathfrak d}^{}\mspace{-2mu}_{ \rho_k}\! \mathcal{E}\, \mathrm{d} \mathrm{y}_k.

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

    T( \mathcal{S} - \sum\limits_{k\in {\mathfrak S}} \rho_k {\mathfrak d}^{}\mspace{-2mu}_{ \rho_k} \mathcal{S}_k) - ( \mathcal{E} - \sum\limits_{k\in {\mathfrak S}} \rho_k {\mathfrak d}^{}\mspace{-2mu}_{ \rho_k} \mathcal{E}_k) = \mathcal{P}, \qquad T {\mathfrak d}^{}\mspace{-2mu}_{ \rho_k}\! \mathcal{S} - {\mathfrak d}^{}\mspace{-2mu}_{ \rho_k}\! \mathcal{E} = - g_k,

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

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

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

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

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

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

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

    \begin{equation} \mathcal{V} = - ( \widetilde{\mathfrak d}^{}\mspace{-2mu}_ \mathcal{X} \mathcal{S})^t = \frac{1}{T} \bigl( -1, g_1, \ldots, g_ n \bigl)^t. \end{equation} (3.19)

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

    \widetilde{\mathsf A}_{0} = \widetilde{\mathfrak d}^{}\mspace{-2mu}_ \mathcal{X} \mathcal{V} = \widetilde{\mathfrak d}^2_{ \mathcal{X} \mathcal{X}} \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,11,12,13,14,15,16,17,18,19].

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

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

    \begin{equation} \mathcal{W} = \Bigl( T, \rho, \frac{ g_2- g_1}{T}, \ldots, \frac{ g_ n- g_1}{T} \Bigl)^t, \end{equation} (3.20)

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

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

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

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

    \begin{equation} {\mathbb R} {1\kern -2.7pt \mathrm{I}}\oplus\Lambda( {1\kern -2.7pt \mathrm{I}}^\perp) = {\mathbb R}^ n. \end{equation} (3.21)

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

    {\mathfrak d}^{}\mspace{-2mu}_ \mathcal{Z} \mathcal{W} = \left( \begin{array}{cc} 1&0_{1, n}\\[3pt] { * _{n,1}}& {\mathcal J} \end{array} \right), \qquad {\mathcal J} = \left( \begin{array}{cccc} 1&1&\cdots&1\\[3pt] \partial_{\rho_1}(\frac{ g_2- g_1}{T})& \partial_{\rho_2}(\frac{ g_2- g_1}{T})& \cdots& \partial_{\rho_ n}(\frac{ g_2- g_1}{T})\\ \vdots&\vdots&&\vdots\\[3pt] \partial_{\rho_1}(\frac{ g_ n- g_1}{T})& \partial_{\rho_2}(\frac{ g_ n- g_1}{T})& \cdots& \partial_{\rho_ n}(\frac{ g_ n- g_1}{T})\\ \end{array} \right),

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

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

    {\mathcal J} \, - \sum\limits_{2 \le i \le n} \partial_{\rho_1}\Bigl(\frac{ g_i- g_1}{T}\Bigr) \, \mathsf{e}_i {\otimes} {1\kern -2.7pt \mathrm{I}} = \left( \begin{array}{cccc} 1&1&\cdots&1\\[4pt] 0_{ n-1,1}&{ {\mathcal J}'}\\ \end{array} \right),

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

    \langle \Lambda \mathcal{v}_i, \mathcal{v}_j \rangle = \Lambda_{i,j} - \Lambda_{1,j} - \Lambda_{i,1} + \Lambda_{1,1}, \qquad 2 \le i,j \le n,

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

    {\mathcal J}' = \rho \left( \begin{array}{ccc} \langle\Lambda \mathcal{v}_2, \mathcal{v}_2\rangle&\cdots& \langle\Lambda \mathcal{v}_2, \mathcal{v}_ n\rangle\\[3pt] \vdots&&\vdots\\ \langle\Lambda \mathcal{v}_ n, \mathcal{v}_2\rangle&\cdots& \langle\Lambda \mathcal{v}_ n, \mathcal{v}_ n\rangle \end{array} \right),

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

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

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

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

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

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

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

    \begin{equation*} \mathcal{W}^\sharp - \mathcal{W}^\flat = \displaystyle{\int}_0^1 {\mathfrak d}^{}\mspace{-2mu}_ \mathcal{Z} \mathcal{W} \bigl( \mathcal{Z}^\flat + \alpha( \mathcal{Z}^\sharp- \mathcal{Z}^\flat) \bigr) \, ( \mathcal{Z}^\sharp- \mathcal{Z}^\flat) \, d\alpha = 0, \end{equation*}

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

    \displaystyle{\int}_0^1 \Bigl\langle {\mathfrak d}^{}\mspace{-2mu}_ \mathcal{Z} \mathcal{W} \bigl( \mathcal{Z}^\flat + \alpha( \mathcal{Z}^\sharp- \mathcal{Z}^\flat) \bigr) \, ( \mathcal{Z}^\sharp- \mathcal{Z}^\flat) \, d\alpha, \mathcal{Z}^\sharp- \mathcal{Z}^\flat \Bigr\rangle = 0,

    and this yields

    \begin{equation} \sum\limits_{i,j\in \mathfrak{S}} \displaystyle{\int}_0^1 \frac{ \partial_{\rho_i} g_j}{T} \bigl( \mathcal{Z}^\flat + \alpha( \mathcal{Z}^\sharp- \mathcal{Z}^\flat) \bigr) d\alpha \; (\rho_i^\sharp - \rho_i^\flat) (\rho_j^\sharp - \rho_j^\flat) = 0. \end{equation} (3.22)

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

    \displaystyle{\int}_0^1 \Bigl\langle \Lambda \bigl( \mathcal{Z}^\flat + \alpha( \mathcal{Z}^\sharp- \mathcal{Z}^\flat) \bigr) d\alpha \,\, ( \varrho^\sharp- \varrho^\flat), \varrho^\sharp- \varrho^\flat\Bigr\rangle = 0,

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

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

    \begin{equation} \mathcal{P}^ {\rm{id}} = R T \sum\limits_{k\in {\mathfrak S}} \frac{ \rho_k}{m_k}, \end{equation} (3.23)

    and the energy per unit volume reads

    \begin{equation} \mathcal{E}^ {\rm{id}} = \sum\limits_{k\in {\mathfrak S}} \rho_k^{} e_k^ {\rm{id}}, \qquad e_k^ {\rm{id}} = e^ \mathrm{st}_k + \displaystyle{\int}_{T^ \mathrm{st}}^T \!\! c_{ {{{\rm{v}}}} k}^ {\rm{id}}(\theta)\,d\theta, \end{equation} (3.24)

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

    \begin{equation} \mathcal{S}^ {\rm{id}} = \sum\limits_{k\in {\mathfrak S}} \rho_k^{} s_k^ {\rm{id}}, \qquad s_k^ {\rm{id}} = s^ \mathrm{st}_k + \displaystyle{\int}_{T^ \mathrm{st}}^T \!\! \frac{ c_{ {{{\rm{v}}}} k}^ {\rm{id}}(\theta)}{\theta} \,d\theta - \frac{ R T}{m_k} \log \frac{ \rho_k}{m_k\gamma^ \mathrm{st}}, \end{equation} (3.25)

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

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

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

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

    \begin{equation} p = p^ {\rm{id}} + \phi. \end{equation} (4.1)

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

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

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

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

    \begin{equation} \vert \partial^ \beta_ \zeta\phi\vert \leqslant \frac{ \textsf{c}( \gamma,T)}{ v^{ \beta_ v+2}}, \qquad 0 \le \vert \beta\vert \le 1 + \gamma. \end{equation} (4.2)

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

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

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

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

    \begin{equation} e = e^ {\rm{id}} - \displaystyle{\int}_ v^\infty \!\!\! T^2 \partial^{}_T \Bigl( \frac{\phi}{T} \Bigr)(T, v', \mathrm{y}_1,\ldots, \mathrm{y}_ n) \, d v', \end{equation} (4.3)
    \begin{equation} s = s^ {\rm{id}} - \displaystyle{\int}_ v^\infty \!\!\! \partial^{}_T\phi(T, v', \mathrm{y}_1,\ldots, \mathrm{y}_ n) \, d v', \end{equation} (4.4)
    \begin{equation} g_k = g_k^ {\rm{id}} + \displaystyle{\int}_ v^\infty \!\!\! \partial^{}_{ \mathrm{y}_k}\phi(T, v', \mathrm{y}_1,\ldots, \mathrm{y}_ n) \, d v', \qquad k\in {\mathfrak S}. \end{equation} (4.5)

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

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

    e(T, v, \mathrm{y}_1, \ldots, \mathrm{y}_ n) - e(T, v'', \mathrm{y}_1, \ldots, \mathrm{y}_ n) = - \displaystyle{\int}_ v^{ v''} \!\!\! T^2 \partial^{}_T \bigl( \frac{\phi}{T} \bigr) (T, v', \mathrm{y}_1, \ldots, \mathrm{y}_ n) \, d v'.

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

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

    \begin{align*} & \bigl( T s - T s^ {\rm{id}} - e + e^ {\rm{id}} \bigr) (T, v, \mathrm{y}_1, \ldots, \mathrm{y}_ n) \\ & - \bigl( T s - T s^ {\rm{id}} - e + e^ {\rm{id}} \bigr) (T, v'', \mathrm{y}_1, \ldots, \mathrm{y}_ n) = - \int_ v^{ v''} \!\!\! \phi(T, v', \mathrm{y}_1, \ldots, \mathrm{y}_ n) \, d v'. \end{align*}

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

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

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

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

    \partial^{}_T e = \partial^{}_T e^ {\rm{id}} - T \displaystyle{\int}_ v^\infty \partial^2_T\phi\, d v', \qquad \partial^{}_T s = \partial^{}_T s^ {\rm{id}} - \displaystyle{\int}_ v^\infty \partial^2_T\phi\, d v',

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

    \partial^{}_ v e = \partial^{}_ v e^ {\rm{id}} + T \partial^{}_T\phi - \phi, \qquad \partial^{}_ v s = \partial^{}_ v s^ {\rm{id}} + \partial^{}_T\phi.

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

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

    \begin{equation*} \partial_T e = \partial_T e^ {\rm{id}} - T \displaystyle{\int}_ v^\infty \!\!\! \partial^2_T\phi\, d v', \end{equation*}

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

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

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

    \begin{equation} \Lambda_{kl} = \frac{ R}{ \mathrm{y}_k m_k} \delta_{kl} + \frac{1}{T} \displaystyle{\int}_ v^\infty \!\!\! \partial^2_{ \mathrm{y}_k \mathrm{y}_l} \phi \, d v', \qquad k,l\in {\mathfrak S}. \end{equation} (4.6)

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

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

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

    \varPi = {{{\rm{diag}}}} \bigl(\sqrt{m_1 \mathrm{y}_1}, \ldots, \sqrt{m_1 \mathrm{y}_ n} \bigr)/\sqrt{R}.

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

    \begin{equation*} \overline\Lambda = \varPi\, \Lambda \; \varPi, \qquad \widehat\Lambda_ \mathsf{r} = \Lambda \, \varPi^2, \qquad \widehat\Lambda_ \mathsf{l} = \varPi^2 \, \Lambda, \end{equation*}

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

    Proof. The proof is presented in [9].

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

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

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

    \begin{equation*} \overline\Lambda = \left( \begin{matrix} \overline\Lambda^{++}&0\\[2pt] 0&{\mathbb I} \end{matrix} \right), \qquad \widehat\Lambda_ \mathsf{r} = \left( \begin{matrix} \widehat\Lambda_ \mathsf{r}^{++}&0\\[5pt] \widehat\Lambda_ \mathsf{r}^{0+}&{\mathbb I} \end{matrix} \right), \qquad \widehat\Lambda_ \mathsf{l} = \left( \begin{matrix} \widehat\Lambda_ \mathsf{l}^{++}&\widehat\Lambda_ \mathsf{l}^{+0}\\[3pt] 0&{\mathbb I} \end{matrix} \right), \end{equation*}

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

    \varPi^{++} = {{{\rm{diag}}}}( \sqrt{m_1 \mathrm{y}_1}, \ldots, \sqrt{m_{ n^+} \mathrm{y}_{ n^+}})/\sqrt{R},

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

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

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

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

    Proof. We refer the reader to [9].

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

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

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

    \begin{equation} p = \sum\limits_{i\in {\mathfrak S}} \frac{ \mathrm{y}_i}{m_i} \frac{R T}{ v - b} - \frac{ a}{ v ( v + b)} \, , \end{equation} (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, \mathrm{y}_1, \ldots, \mathrm{y}_{ n}) and b (\mathrm{y}_1, \ldots, \mathrm{y}_{ n}) may be evaluated by using the Van der Waals mixing rules:

    \begin{equation} a = \sum\limits_{i,j\in {\mathfrak S}} \mathrm{y}_i \mathrm{y}_j \alpha_i \alpha_j \qquad b = \sum\limits_{i\in {\mathfrak S}} \mathrm{y}_i b_i, \end{equation} (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\geqslant2 .

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

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

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

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

    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

    \begin{equation} \phi = \Bigl( \sum\limits_{i\in {\mathfrak S}} \frac{RT \mathrm{y}_i}{m_i} \Bigr) \; \frac{ b}{ v( v - b)} - \frac{ a}{ v ( v + b)}. \end{equation} (5.3)

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

    \begin{equation} \mathcal{O}_ \zeta = \ \{\; \zeta\in(0,\infty)^{2+ n}; \; v > b \;\}, \end{equation} (5.4)

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

    Proof. The proof is presented in [9].

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

    \begin{equation} e = \sum\limits_{i\in {\mathfrak S}} \mathrm{y}_i e_i^{ {\rm{id}}} + \bigl( T \partial^{}_T a - a \bigr) \frac{1}{ b} \log \bigl( 1+ \frac{ b}{ v} \bigr) , \end{equation} (5.5)

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

    \begin{equation} s = \sum\limits_{i\in {\mathfrak S}} \mathrm{y}_i s_i^{ {\rm{id}}\star} - \sum\limits_{i\in {\mathfrak S}} \frac{ \mathrm{y}_i R}{m_i} \log \Bigl( \frac{ \mathrm{y}_i R T}{m_i ( v- b) p^ \mathrm{st}} \Bigr) + \frac{ \partial^{}_T a}{ b} \log \bigl( 1 + \frac{ b}{ v} \bigr), \end{equation} (5.6)

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

    \begin{align} \Lambda_{ij} = & \frac{R\delta_{ij}}{m_i \mathrm{y}_i} + \frac{R}{ v - b} \bigl( \frac{ b_{i}}{m_j} + \frac{ b_{j}}{m_i} \bigr) + \sum\limits_{k\in {\mathfrak S}} \frac{ \mathrm{y}_k}{m_k} \frac{R}{( v - b)^2} b_{i} b_{j} - \frac{2}{T} \frac{ \alpha_i \alpha_j}{ b} \log \bigl( 1+ \frac{ b}{ v} \bigr) \\[7pt] & + \frac{2}{T} \sum\limits_{k\in {\mathfrak S}} \mathrm{y}_k \bigl( \alpha_i \alpha_k b_{j} + \alpha_j \alpha_k b_{i} \bigr) \Bigl( \frac{1}{ b^2} \log \bigl( 1+ \frac{ b}{ v}\bigr) - \frac{ 1}{ b ( v + b) } \Bigr) \\[7pt] & + \frac{1}{T} a b_{i} b_{j} \Bigl( - \frac{2}{ b^3} \log \bigl( 1+ \frac{ b}{ v} \bigr) + \frac{ 2}{ b^2 ( v + b)} + \frac{1}{ b ( v + b)^2 } \Bigr), \qquad i,j\in {\mathfrak S}. \end{align} (5.7)

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

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

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

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

    Proposition 5.4. Assume that \alpha_i is in the form

    \begin{equation} \frac{ \alpha_i}{ \alpha_{ {{{\rm{c}}}},i}} = 1 + A \bigl( {\rm{s}}_{i} ( 1 - \sqrt{T}/\sqrt{T_{ {{{\rm{c}}}},i}}) \bigr), \end{equation} (5.8)

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

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

    \frac{ \partial^{}_T \alpha_i}{ \alpha_{ {{{\rm{c}}}},i}} = - \frac{ {\rm{s}}_{i}}{2(TT_{ {{{\rm{c}}}},i})^{1/2}} A' \bigl( {\rm{s}}_{i} ( 1 - \sqrt{T}/\sqrt{T_{ {{{\rm{c}}}},i}}) \bigr),
    \frac{ \partial_{TT}^2 \alpha_i}{ \alpha_{ {{{\rm{c}}}},i}} = \frac{ {\rm{s}}_{i}}{ 4T^{3/2}_{}T_{ {{{\rm{c}}}},i}^{1/2}} A' \bigl( {\rm{s}}_{i} ( 1 - \sqrt{T}/\sqrt{T_{ {{{\rm{c}}}},i}}) \bigr) + \frac{ {\rm{s}}_{i}^2}{4TT_{ {{{\rm{c}}}},i}} A'' \bigl( {\rm{s}}_{i} ( 1 - \sqrt{T}/\sqrt{T_{ {{{\rm{c}}}},i}}) \bigr).

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

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

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

    \begin{equation} p = \frac{ r T}{ v - b} - \frac{ a}{ v ( v+ b)}. \end{equation} (5.9)

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

    \begin{equation} T^\star = \frac{ a^\star}{ r} \frac{( v - b)^2 (2 v+ b)}{ v^2 ( v+ b)^2}, \end{equation} (5.10)

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

    \begin{equation} \partial^2_{T v} p \; \partial^{}_ v T^\star + \partial^2_{ v v} p = 0, \qquad \partial^2_{T v} p \; \partial^{}_{ \mathrm{y}_k} \! T^\star + \partial^2_{ v\, \mathrm{y}_k} p = 0, \quad k\in {\mathfrak S}. \end{equation} (5.11)

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

    Proof. The proof is presented in [9].

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

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

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

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

    \begin{equation} T^{\star\star}( \mathrm{y}) = T^\star(\theta_ {{{\rm{c}}}} b, \mathrm{y}), \qquad p^{\star\star}( \mathrm{y}) = p^\star \bigl( \theta_ {{{\rm{c}}}} b, \mathrm{y} \bigr) = p \bigl( T^\star(\theta_ {{{\rm{c}}}} b, \mathrm{y}), \theta_ {{{\rm{c}}}} b, \mathrm{y} \bigr), \end{equation} (5.12)

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

    Proof. We refer the reader to [9].

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

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

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

    \begin{equation} \mathcal{F}^ {\rm{ex}} = \mathcal{F}^ {} + {\tfrac{1}{2}} \varkappa \vert \boldsymbol{\nabla}\mspace{-3mu}\rho\vert^2, \end{equation} (6.1)

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

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

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

    \begin{equation} \mathcal{E}^ {\rm{ex}} = \mathcal{E}^ {} + {\tfrac{1}{2}}( \varkappa - T \partial_T \varkappa) \vert {\boldsymbol\nabla}\rho\vert^2, \qquad p^ {\rm{ex}} = p^ {} - {\tfrac{1}{2}} \varkappa \vert {\boldsymbol\nabla}\rho\vert^2. \end{equation} (6.2)

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

    \begin{equation} T \mathrm{d} \mathcal{S}^ {\rm{ex}} = \mathrm{d} \mathcal{E}^ {\rm{ex}} - \sum\limits_{i\in \mathfrak{S}} g^ {}_i \mathrm{d} \rho_i - \varkappa \boldsymbol{\nabla}\mspace{-3mu}\rho {\boldsymbol{\cdot}} \mathrm{d} \boldsymbol{\nabla}\mspace{-3mu}\rho, \end{equation} (6.3)

    where \mathrm{d} is the differentiation operator.

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

    \begin{align} & \partial_t^{}\rho_i + \boldsymbol{\nabla} {\boldsymbol{\cdot}}(\rho_i {\boldsymbol{u}}) + \boldsymbol{\nabla} {\boldsymbol{\cdot}} \boldsymbol{\mathcal{J}}_i = m_i \omega_i, \qquad \quad i\in \mathfrak{S} = \{1,\ldots, n\}, \end{align} (6.4)
    \begin{align} [3pt] & \partial_t^{}(\rho {\boldsymbol{u}}) + \boldsymbol{\nabla} {\boldsymbol{\cdot}}(\rho {\boldsymbol{u}}{\otimes} {\boldsymbol{u}}) + \boldsymbol{\nabla} {\boldsymbol{\cdot}} \boldsymbol{\mathcal{P}} = 0, \end{align} (6.5)
    \begin{align} [3pt] & \partial_t^{}\bigl( \mathcal{E}^ {\rm{ex}} + {\tfrac{1}{2}}\rho\vert {\boldsymbol{u}}\vert^2 \bigr) + \boldsymbol{\nabla} {\boldsymbol{\cdot}} \bigl( {\boldsymbol{u}} ( \mathcal{E}^ {\rm{ex}} + {\tfrac{1}{2}}\rho\vert {\boldsymbol{u}}\vert^2) \bigr) + \boldsymbol{\nabla} {\boldsymbol{\cdot}}( \boldsymbol{\mathcal{Q}} + \boldsymbol{\mathcal{P}} {\boldsymbol{\cdot}} {\boldsymbol{u}}) = 0, \end{align} (6.6)

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

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

    \begin{equation} \partial_t^{} \rho + \boldsymbol{\nabla} {\boldsymbol{\cdot}}(\rho {\boldsymbol{u}}) = 0. \end{equation} (6.7)

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

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

    \begin{align} \boldsymbol{\mathcal{P}} & = p^ {\rm{ex}} \boldsymbol{I} + \; \varkappa \boldsymbol{\nabla}\mspace{-3mu}\rho {\otimes} \boldsymbol{\nabla}\mspace{-3mu}\rho - \rho \boldsymbol{\nabla} {\boldsymbol{\cdot}}( \varkappa \boldsymbol{\nabla}\mspace{-3mu}\rho) \boldsymbol{I} + \boldsymbol{\mathcal{P}}^{{\rm{d}}}, \end{align} (6.8)
    \begin{align} \boldsymbol{\mathcal{Q}} & = \varkappa \rho \boldsymbol{\nabla} {\boldsymbol{\cdot}} {\boldsymbol{u}} \, \boldsymbol{\nabla}\mspace{-3mu}\rho + \boldsymbol{q}, \end{align} (6.9)

    with the dissipative fluxes given by

    \begin{align} \boldsymbol{\mathcal{P}}^{{\rm{d}}} = & - \mathfrak{v} \boldsymbol{\nabla} {\boldsymbol{\cdot}} {\boldsymbol{u}} \, \boldsymbol{I} - \eta \bigl( \boldsymbol{\nabla} {\boldsymbol{u}} + \boldsymbol{\nabla} {\boldsymbol{u}}^t - {\tfrac{2}{3}} \boldsymbol{\nabla} {\boldsymbol{\cdot}} {\boldsymbol{u}} \, \boldsymbol{I} \bigr), \end{align} (6.10)
    \begin{align} \boldsymbol{\mathcal{J}}_i = & - \sum\limits_{j\in \mathfrak{S}} L_{ij} \boldsymbol{\nabla}\Bigl( \frac{g_j^ {}}{T} \Bigr) - L_{i \mathrm{e}} \boldsymbol{\nabla}\Bigl( \frac{-1}{T} \Bigr), \qquad i\in \mathfrak{S}, \end{align} (6.11)
    \begin{align} \boldsymbol{q} = & - \sum\limits_{i\in \mathfrak{S}} L_{ \mathrm{e} i} \boldsymbol{\nabla}\Bigl( \frac{g_i^ {}}{T} \Bigr) - L_{ \mathrm{e} \mathrm{e}} \boldsymbol{\nabla}\Bigl( \frac{-1}{T} \Bigr), \end{align} (6.12)

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

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

    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

    \begin{equation} \sum\limits_{i\in \mathfrak{S}} \nu_{ij}^ {{{\rm{f}}}} \mathfrak{M}_i \rightleftarrows \sum\limits_{i\in \mathfrak{S}} \nu_{ij}^ {{{\rm{b}}}} \mathfrak{M}_i, \qquad j\in {\mathfrak R}, \end{equation} (6.13)

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

    \omega_i = \sum\limits_{j\in {\mathfrak R}} \nu_{ij} \tau_j,

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

    \nu_j = \left( \begin{array}{c} \nu_{1j}\\ \vdots\\ \nu_{ n j} \end{array} \right), \qquad \nu_j^ {{{\rm{f}}}} = \left( \begin{array}{c} \nu_{1j}^ {{{\rm{f}}}}\\ \vdots\\ \nu_{ n j}^ {{{\rm{f}}}} \end{array} \right), \qquad \nu_j^ {{{\rm{b}}}} = \left( \begin{array}{c} \nu_{1j}^ {{{\rm{b}}}}\\ \vdots\\ \nu_{ n j}^ {{{\rm{b}}}} \end{array} \right), \qquad \mu = \left( \begin{array}{c} \mu_1\\ \vdots\\ \mu_ n \end{array} \right),

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

    \begin{equation} \tau_j = \mathcal{K}_j \bigl( \exp\langle\nu_j^ {{{\rm{f}}}},\mu\rangle - \exp\langle\nu_j^ {{{\rm{b}}}},\mu\rangle \bigr), \qquad j\in {\mathfrak R}, \end{equation} (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 \mathcal{S}^ {\rm{ex}} is then obtained in the form

    \begin{align} \partial_t^{} \mathcal{S}^ {\rm{ex}} + \boldsymbol{\nabla} {\boldsymbol{\cdot}}( {\boldsymbol{u}} \mathcal{S}^ {\rm{ex}}) & + \boldsymbol{\nabla} {\boldsymbol{\cdot}} \Bigl( \frac{ \boldsymbol{\mathcal{Q}}}{T} - \frac{ \varkappa\rho \boldsymbol{\nabla}\mspace{-3mu}\rho \boldsymbol{\nabla} {\boldsymbol{\cdot}} {\boldsymbol{u}}}{T} - \sum\limits_{i\in \mathfrak{S}} \frac{ g_i}{T} \boldsymbol{\mathcal{J}}_i \Bigr) \; = \\ & - \frac{1}{T} \Bigl( \boldsymbol{\mathcal{P}} - p \boldsymbol{I} - ( \varkappa \boldsymbol{\nabla}\mspace{-3mu}\rho {\otimes} \boldsymbol{\nabla}\mspace{-3mu}\rho - \rho \boldsymbol{\nabla} {\boldsymbol{\cdot}}( \varkappa \boldsymbol{\nabla}\mspace{-3mu}\rho) \boldsymbol{I}) \Bigr) \boldsymbol{:} \boldsymbol{\nabla} {\boldsymbol{u}}\\ & \; - \; \Bigl( \boldsymbol{\mathcal{Q}} - \varkappa\rho \boldsymbol{\nabla}\mspace{-3mu}\rho \boldsymbol{\nabla} {\boldsymbol{\cdot}} {\boldsymbol{u}} \Bigr) {\boldsymbol{\cdot}} \boldsymbol{\nabla}\Bigl(\frac{-1}{T}\Bigr) \; - \; \sum\limits_{i\in \mathfrak{S}} \boldsymbol{\mathcal{J}}_i {\boldsymbol{\cdot}} {\boldsymbol\nabla} \Bigl(\frac{ g_i}{T}\Bigr) \; - \; \sum\limits_{i\in \mathfrak{S}} \frac{g_i m_i \omega_i}{T}. \end{align} (6.15)

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

    \begin{align*} \partial_t^{} \mathcal{S}^ {\rm{ex}} & + \boldsymbol{\nabla} {\boldsymbol{\cdot}}( {\boldsymbol{u}} \mathcal{S}^ {\rm{ex}}) + \boldsymbol{\nabla} {\boldsymbol{\cdot}} \Bigl( \frac{ \boldsymbol{q}}{T} - \sum\limits_{i\in \mathfrak{S}} \frac{ g_i}{T} \boldsymbol{\mathcal{J}}_i \Bigr) = \frac{ \mathfrak{v}}{T} ( \boldsymbol{\nabla} {\boldsymbol{\cdot}} {\boldsymbol{u}})^2 + \frac{\eta}{2T} \vert \boldsymbol{\mathsf S}\vert^2 \\ & + \!\! \sum\limits_{i,j\in \mathfrak{S}\cup\{ \mathrm{e}\}} \!\!\!\! L_{ij} \boldsymbol{\nabla}\Bigl(\frac{g_i}{T}\Bigr) {\boldsymbol{\cdot}} \boldsymbol{\nabla}\Bigl(\frac{g_j}{T}\Bigr) + \sum\limits_{j\in {\mathfrak R}} R \mathcal{K}_j \bigl( \langle\nu_j^ {{{\rm{f}}}},\mu\rangle - \langle\nu_j^ {{{\rm{b}}}},\mu\rangle \bigr) \bigl( \exp\langle\nu_j^ {{{\rm{f}}}},\mu\rangle - \exp\langle\nu_j^ {{{\rm{b}}}},\mu\rangle \bigr), \end{align*}

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

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

    \begin{equation} {\boldsymbol{w}} = {\boldsymbol\nabla}\rho. \end{equation} (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 {\boldsymbol{w}} is obtained by applying the differential operator {\boldsymbol\nabla} to (6.7) and reads

    \begin{equation} \partial_t {\boldsymbol{w}} + \sum\limits_{i\in \mathcal{D}} \partial_i^{} \bigl( {\boldsymbol{w}} \, u_i + \rho {\boldsymbol\nabla} u_i \bigr) = 0, \end{equation} (6.17)

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

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

    \begin{equation} \mathsf{u} = \bigl( \rho_1, \ldots, \rho_ n, {\boldsymbol{w}}, \rho {\boldsymbol{u}}, \mathcal{E} + {\tfrac{1}{2}}\rho\vert {\boldsymbol{u}}\vert^2 \bigr)^t, \qquad \mathsf{z} = \bigl( \rho_1, \ldots, \rho_ n, {\boldsymbol{w}}, {\boldsymbol{u}}, T \bigr)^t, \end{equation} (6.18)

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

    \begin{equation} \partial_t \mathsf{u} + \sum\limits_{i\in \mathcal{D}} \partial_i^{} \bigl( {\mathsf F}_{i}+ {\mathsf F}_{i}^{ {\rm{c}}}+ {\mathsf F}_{i}^{ {\rm{d}}}\bigr) = \Omega, \end{equation} (6.19)

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

    \begin{align} {\mathsf F}_{i}& = \Bigl( \rho_1 u_i, \ldots, \rho_ n u_i, {\boldsymbol{w}} u_i, \rho {\boldsymbol{u}} u_i + (p^ {\rm{ex}} + \varkappa\vert {\boldsymbol{w}}\vert^2) \boldsymbol{\mathsf{e}}_i, ( \mathcal{E}^ {\rm{ex}} + {\tfrac{1}{2}}\rho\vert {\boldsymbol{u}}\vert^2 + p^ {\rm{ex}} + \varkappa\vert {\boldsymbol{w}}\vert^2) u_i \Bigr)^t, \qquad i\in \mathcal{D} \end{align} (6.20)
    \begin{align} {\mathsf F}_{i}^{ {\rm{c}}}& = \Bigl( 0, \ldots, 0, \rho {\boldsymbol\nabla} u_i, - \rho {\boldsymbol\nabla}( \varkappa w_i), \rho \varkappa {\boldsymbol{w}} {\boldsymbol{\cdot}} {\boldsymbol\nabla} u_i - \rho {\boldsymbol{u}} {\boldsymbol{\cdot}} {\boldsymbol\nabla}( \varkappa w_i) \Bigr)^t, \end{align} (6.21)
    \begin{align} {\mathsf F}_{i}^{ {\rm{d}}}& = \Bigl( \mathcal{J}_{1i}, \ldots, \mathcal{J}_{ n i}, 0_{ d,1}, \, \boldsymbol{\mathcal{P}}^{{\rm{d}}}\!\!\!_i, \, q_i + \sum\limits_{j\in \mathcal{D}} \mathcal{P}^{{\rm{d}}}_{ij} u_j \Bigr)^t, \end{align} (6.22)
    \begin{align} \Omega & = \Bigl( m_1\omega_1 \ldots, m_ n\omega_ n, 0,0,0\Bigr)^t, \end{align} (6.23)

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

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

    \begin{equation} {\mathsf F}_{i}^{ {\rm{d}}} = - \sum\limits_{j\in \mathcal{D}} {\mathsf B}^{ {\rm{d}}}_{ij} \partial_j^{} \mathsf{u}, \qquad {\mathsf F}_{i}^{ {\rm{c}}} = - \sum\limits_{j\in \mathcal{D}} {\mathsf B}^{ {\rm{c}}}_{ij} \partial_j^{} \mathsf{u}, \qquad i\in \mathcal{D}, \end{equation} (6.24)

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

    \begin{equation} \partial_t^{} \mathsf{u} + \sum\limits_{i\in \mathcal{D}} {\mathsf A}_{i}( \mathsf{u}) \partial_i^{} \mathsf{u} - \sum\limits_{i,j\in \mathcal{D}} \partial_i^{}\bigl( {\mathsf B}^{ {\rm{d}}}_{ij}( \mathsf{u}) \partial_j^{} \mathsf{u} \bigr) - \sum\limits_{i,j\in \mathcal{D}} \partial_i^{}\bigl( {\mathsf B}^{ {\rm{c}}}_{ij}( \mathsf{u}) \partial_j^{} \mathsf{u} \bigr) = \Omega( \mathsf{u}). \end{equation} (6.25)

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

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

    \begin{align*} \mathcal{E}^ {\rm{ex}} & = \mathcal{E}^ {} + {\tfrac{1}{2}} ( \varkappa-T \partial_T \varkappa) \vert {\boldsymbol{w}}\vert^2, \qquad p^ {\rm{ex}} = p^ {} - {\tfrac{1}{2}} \varkappa \vert {\boldsymbol{w}}\vert^2, \\[3pt] \mathcal{S}^ {\rm{ex}} & = \mathcal{S}^ {} - {\tfrac{1}{2}} \partial_T \varkappa \vert {\boldsymbol{w}}\vert^2, \qquad\quad\quad\quad g_k^ {\rm{ex}} = g^ {}_k, \qquad k\in \mathfrak{S}. \end{align*}

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

    \begin{equation} \Lambda = (\Lambda_{ij})_{i,j\in \mathfrak{S}}, \qquad \Lambda_{ij} = \frac{ \partial_{\rho_j} g_i}{T} = \Lambda^ {\rm{ex}}_{ij}. \end{equation} (6.26)

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

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

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

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

    \begin{equation} \lim\limits_{ \rho\to0, {\boldsymbol{w}} = 0} \tfrac{1}{ \rho} ( \mathcal{E}^ {\rm{ex}}- \mathcal{E}^ {\rm{id}}) = 0, \qquad \lim\limits_{ \rho\to0, {\boldsymbol{w}} = 0} \tfrac{1}{ \rho} ( \mathcal{P}^ {\rm{ex}}- \mathcal{P}^ {\rm{id}}) = 0, \qquad \lim\limits_{ \rho\to0, {\boldsymbol{w}} = 0} \tfrac{1}{ \rho} ( \mathcal{S}^ {\rm{ex}}- \mathcal{S}^ {\rm{id}}) = 0. \end{equation} (6.27)

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

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

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

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

    Proof. This is a consequence of the expression

    \begin{align} \langle \partial^2_{ \mathsf{u} \mathsf{u}} \mathcal{S}^ {\rm{ex}} \, \mathsf{X} , \, \mathsf{X} \rangle = & - \frac{1}{T^2 \partial_T \mathcal{E}^ {\rm{ex}}} \Bigl( \mathsf{X}_ \mathcal{E} - ( \partial_\varrho \mathcal{E}^ {\rm{ex}} - {\tfrac{1}{2}}\vert {\boldsymbol{u}}\vert^2 {1\kern -2.7pt \mathrm{I}}) \mathsf{X}_\varrho - {\boldsymbol{u}} {\boldsymbol{\cdot}} \mathsf{X}_ {\boldsymbol{u}} - ( \varkappa-T \partial_T \varkappa) {\boldsymbol{w}} {\boldsymbol{\cdot}} \mathsf{X}_ {\boldsymbol{w}} \Bigr)^2 \\[3pt] & - \bigl\langle\Lambda \mathsf{X}_\varrho, \mathsf{X}_\varrho \bigr\rangle - \frac{ \varkappa}{T} \, \bigl\vert \mathsf{X}_ {\boldsymbol{w}}\bigr\vert^2 - \frac{1}{\rho T} \bigl\vert \mathsf{X}_ {\boldsymbol{u}} - {\boldsymbol{u}} \mathsf{X}_\rho \bigr\vert^2, \end{align} (6.28)

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

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

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

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

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

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

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

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

    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.

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

    \begin{equation} p^ {} = \overline p^ {} + \epsilon^ 2 \widetilde p^ {}, \end{equation} (7.1)

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

    \begin{equation} \boldsymbol{\nabla} {\boldsymbol{\cdot}}\Bigl( \overline p^ {} \boldsymbol{I} - {\tfrac{1}{2}} \varkappa \vert \boldsymbol{\nabla}\mspace{-3mu}\rho\vert^ 2 \boldsymbol{I} + \varkappa \boldsymbol{\nabla}\mspace{-3mu}\rho {\otimes} \boldsymbol{\nabla}\mspace{-3mu}\rho - \varkappa \rho \Delta\rho \boldsymbol{I} \Bigr) = 0. \end{equation} (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 \widetilde p plays a role in the momentum equation, and from p^ {\rm{ex}} = p^ {} - {\tfrac{1}{2}} \varkappa \vert \boldsymbol{\nabla}\mspace{-3mu}\rho\vert^2 it is obtained that p^ {\rm{ex}} = \overline p^ {} - {\tfrac{1}{2}} \varkappa \vert \boldsymbol{\nabla}\mspace{-3mu}\rho\vert^2 + \epsilon^ 2 \widetilde p . For standard fluids where \varkappa = 0 , we then recover that the ambiant pressure \overline p^ {} is spatially uniform [8] since \boldsymbol{\nabla} {\boldsymbol{\cdot}}(\overline p^ {} \boldsymbol{I}) = \boldsymbol{\nabla} \overline p^ {} = 0 . When capillary effects are included, the small Mach number relation (7.2) is not practical without further simplifying assumptions.

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

    \overline {\boldsymbol{\mathcal{P}}} = \bigl( p^\infty - \varkappa \rho^{\prime \, 2} \bigr) \bigl( \boldsymbol{I} - \mathsf{e}_ y{\otimes} \mathsf{e}_ y \bigr) + p^\infty \mathsf{e}_ y{\otimes} \mathsf{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} \varkappa \rho^{\prime \, 2} \, d y may then be interpreted as surface tension [25,27,28].

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

    \begin{equation} \rho \partial_t^{} h + \rho {\boldsymbol{u}} \boldsymbol{\nabla} {\boldsymbol{\cdot}} h + \boldsymbol{\nabla} {\boldsymbol{\cdot}} \boldsymbol{q} = \partial_t^{} \overline p^ {} + {\boldsymbol{u}} {\boldsymbol{\cdot}} \boldsymbol{\nabla} \overline p^ {}, \end{equation} (7.3)

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

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

    \begin{align*} T & = T(t, y), & \rho & = \rho(t, y), & \overline p^ {} & = \overline p^ {}(t, y),\\[3pt] {\rm{U}} & = x \, \widetilde{\rm{U}}(t, y), & {\rm{V}} & = {\rm{V}}(t, y), & \mathrm{y}_i & = \mathrm{y}_i(t, y),\\[3pt] {\widetilde p}^ {} & = - {\tfrac{1}{2}} J x^2 + {\widehat p}^ {}(t, y), \end{align*}

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

    \begin{align} \partial_t \rho + 2^\delta \rho {\widetilde {\rm{U}}} + \partial_ y^{}(\rho {\rm{V}}) & = 0, \end{align} (7.4)
    \begin{align} \rho \partial_t \mathrm{y}_i + \rho {\rm{V}} \partial_ y^{} \mathrm{y}_i + \partial_ y^{} \mathcal{J}_i & = 0, \qquad i\in \mathfrak{S}, \end{align} (7.5)
    \begin{align} \rho \partial_t \widetilde{\rm{U}} + \rho \widetilde{\rm{U}}^2 + \rho {\rm{V}} \partial_ y^{} \widetilde{\rm{U}} - J + \partial_ y^{}(\eta \partial_ y^{} \widetilde{\rm{U}}) & = 0, \end{align} (7.6)
    \begin{align} \rho \partial_t h + \rho v \partial_ y^{} h - {\rm{V}} \partial_ y \overline p + \partial_ y^{} q & = 0, \end{align} (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

    \begin{align} \mathcal{J}_i & = - \sum\limits_{j\in \mathfrak{S}} L_{ij} \partial_ y^{} \Bigl( \frac{g_j^ {}}{T} \Bigr) - L_{i \mathrm{e}} \partial_ y^{} \Bigl( \frac{-1}{T} \Bigr), \qquad i\in \mathfrak{S}, \end{align} (7.8)
    \begin{align} q & = - \sum\limits_{i\in \mathfrak{S}} L_{ \mathrm{e} i} \partial_ y^{} \Bigl( \frac{g_i^ {}}{T} \Bigr) - L_{ \mathrm{e} \mathrm{e}} \partial_ y^{} \Bigl( \frac{-1}{T} \Bigr), \end{align} (7.9)

    and it has been established in Section 7.1 that

    \begin{equation} \overline p^ {} = p^\infty - {\tfrac{1}{2}} \varkappa ( \partial_ y^{}\rho)^ 2 + \varkappa \rho \, \partial_ y^2\rho, \end{equation} (7.10)

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

    \begin{align} T(-\infty) & = T^ -, \qquad & T(+\infty) & = T^ +, \end{align} (7.11)
    \begin{align} \mathrm{y}_i(-\infty) & = \mathrm{y}_i^ -, \qquad & \mathrm{y}_i(+\infty) & = \mathrm{y}_i^ +, \end{align} (7.12)
    \begin{align} \widetilde u(-\infty) & = \alpha \sqrt{\rho^ +/\rho^ -}, \qquad & \widetilde u(+\infty) & = \alpha, \end{align} (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

    \begin{equation} \alpha = (J / \rho^ +)^{1/2}. \end{equation} (7.14)

    The stagnation plane is finally located at the origin with

    \begin{equation} {\rm{V}}(0) = 0. \end{equation} (7.15)

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

    \begin{align} 2^\delta \rho {\widetilde {\rm{U}}} + \partial_ y^{}(\rho {\rm{V}}) & = 0, \end{align} (7.16)
    \begin{align} \rho {\rm{V}} \partial_ y^{} \mathrm{y}_i + \partial_ y^{} \mathcal{J}_i & = 0, \qquad i\in \mathfrak{S}, \end{align} (7.17)
    \begin{align} \rho \widetilde{\rm{U}}^2 + \rho {\rm{V}} \partial_ y^{} \widetilde{\rm{U}} - J + \partial_ y^{}(\eta \partial_ y^{} \widetilde{\rm{U}}) & = 0, \end{align} (7.18)
    \begin{align} \rho v \partial_ y^{} h - {\rm{V}} \partial_ y \overline p + \partial_ y^{} q & = 0, \end{align} (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.

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

    \begin{equation} {\mu_{{i}}} = \log \mathsf{X}_i + \mu_i^ \mathrm{sm}, \end{equation} (7.20)

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

    \begin{eqnarray} {\mu_{{i}}}^{ \mathrm{sm}} & = & \frac{m_i g_i^{ {\rm{id}}\star}}{RT} + \log \Bigl( \frac{RT} {( v - b) p^ \mathrm{st} m} \Bigr) + \sum\limits_{j\in \mathfrak{S}} \frac{ \mathrm{y}_j}{m_j} \frac{m_i b_i}{ v - b } - \frac{m_i a b_i} {RT b ( v + b )} \\ && \qquad + \frac{ a b_i - b \partial_{ \mathrm{y}_i} a}{ b^2} \frac{m_i}{RT} \log \bigl( 1+ \frac{ b}{ v} \bigr), \end{eqnarray} (7.21)

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

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

    \begin{equation} L = \left( \begin{array}{cc} \mathcal{D}& \mathcal{D} \mathcal{h}\\[5pt] ( \mathcal{D} \mathcal{h})^t&\ \ \lambda T^2 + \langle \mathcal{D} \mathcal{h}, \mathcal{h}\rangle\\ \end{array} \right), \end{equation} (7.22)

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

    \begin{equation} \mathcal{D}_{ij} = \rho \mathrm{y}_i \mathrm{y}_j \frac{m}{R}D_{ij}, \qquad \mathcal{h}_i = h_i + RT\frac{\widetilde\chi_i}{m_i}, \end{equation} (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 = \widehat\partial^{}_{ \mathrm{y}_i} h , i\in \mathfrak{S} , where the derivative is thus taken at constant pressure, and \widetilde\chi_i , i\in \mathfrak{S} , the reduced thermal diffusion ratios [18,80]. The multicomponent diffusion coefficients D_{ij} are obtained with an approximate inversion [8,105,106,107,108,109,110,111] using

    D = (\Delta + \mathrm{y}{\otimes} \mathrm{y})^{-1} - {1\kern -2.7pt \mathrm{I}}{\otimes} {1\kern -2.7pt \mathrm{I}},

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

    \Delta_{kk} = \sum\limits_{l\ne k} \frac{ \mathsf{X}_k \mathsf{X}_l}{{\mathcal D}^ {{{\rm{bin}}}}_{kl} }, \qquad\quad \Delta_{kl} = - \frac{ \mathsf{X}_k \mathsf{X}_l}{{\mathcal D}^ {{{\rm{bin}}}}_{kl}}, \quad \ k \ne l,

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

    \Upsilon_{ij} = 1 + \sum\limits_{k\in \mathfrak{S}} \frac{\pi {\mathfrak n}_k}{12} \Bigl( 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 \Bigr).

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

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

    \begin{equation} h_k = \widehat\partial^{}_{ \mathrm{y}_k}h = \partial^{}_{ \mathrm{y}_k}h - \partial^{}_\rho h \frac{ \partial^{}_{ \mathrm{y}_k} \overline p^ {}} { \partial^{}_\rho \overline p^ {}}, \end{equation} (7.24)

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

    \begin{equation} \partial^{}_{ \mathrm{y}_k} \overline p^ {} = \frac{\rho m}{m_k} \partial^{}_\rho \overline p^ {} + \mathcal{R}_k, \end{equation} (7.25)

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

    \mathcal{R}_k = \Bigl\{ \frac{RT}{( v-b)^2} + \frac{ma}{ v( v+b)^2} \Bigr\} \! \sum\limits_{j\in \mathfrak{S}} \mathrm{y}_j \Bigl( \frac{b_k}{m_j} - \frac{b_j}{m_k} \Bigr) + \frac{2m}{ v( v+b)} \sum\limits_{j,l\in \mathfrak{S}} \! \mathrm{y}_j \mathrm{y}_l \alpha_l \Bigl( \frac{\alpha_j}{m_k} - \frac{\alpha_k}{m_j} \Bigr).

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

    \widetilde h_k = \partial^{}_{ \mathrm{y}_k} h - \frac{\rho m}{m_k} \partial^{}_\rho h - \frac{ \mathcal{R}_k}{ \mathfrak{f}( \partial^{}_\rho \overline p^ {})} \partial^{}_\rho h,

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

    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.

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

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

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

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

    \begin{equation} \alpha_{ {{{\rm{c}}}},i}^2 = a_{ {{{\rm{c}}}},i} = 0.42748 \frac{ R^2 T_{ {{{\rm{c}}}},i}^2}{m_i^2 p_{ {{{\rm{c}}}},i}}, \qquad b_{i} = 0.08664 \frac{R T_{ {{{\rm{c}}}},i}}{m_i p_{ {{{\rm{c}}}},i}}, \end{equation} (8.1)

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

    \begin{equation} {\rm{s}}_{i} = 0.48508+1.5517\varpi_{{i}} - 0.151613\varpi_{{i}}^2, \end{equation} (8.2)

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

    Table 1.  Critical state and acentric factor of ethane and nitrogen from NIST.
    Species T_{ {{{\rm{c}}}}, i} (K) p_{ {{{\rm{c}}}}, i} (bar) \varpi_{{i}}
    C _2 H _6 305.3 49.0 0.098
    N _2 126.1 39.98 0.0372

     | Show Table
    DownLoad: CSV

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

    \begin{equation} a_{i} \bigl(T_{ {{{\rm{c}}}},i}\bigr) = \bigl(5.55 \pm 0.12 \bigr) \frac{ {\rm{N}}^2 \epsilon_{{i}} \sigma_{{i}}^3}{m_i^2}, \qquad b_{i} = \bigl( 0.855 \pm 0.018 \bigr) \frac{ {\rm{N}} \sigma_{{i}}^3}{m_i}, \end{equation} (8.3)

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

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

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

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

    Figure 1.  Constant pressure specific heat of pure nitrogen; ____ SRK; _ _ _ Ideal gas; +++ experimental measurements from Younglove [121] and Younglove and Ely [122].

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

    Figure 2.  Saturation pressure of pure nitrogen; ____ SRK; _ _ _ mechanical stability limits; +++ experimental measurements form Younglove and Ely [122].

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

    \begin{equation} \left\{ \begin{array}{l} T = T_ {{{\rm{ref}}}},\\[3pt] \det\Lambda(T, v, \mathrm{y}),\\[3pt] \mathrm{y}_1 + \mathrm{y}_2 = 1, \end{array} \right. \end{equation} (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 \widehat\Lambda^ \mathsf{r} of Proposition 4.2 that is more convenient than \Lambda . The thermodynamic stability domain and the mechanical stability domain are presented at temperature T = 200 K in the plane (\mathrm{y}_{{\rm C}_2{\rm H}_6}, v) on Figure 3. The unstable zones are inside the curves and the mechanical unstable domain is naturally interior to the thermodynamic stability domain.

    Figure 3.  Stability limit domains of ethane-nitrogen mixtures at T = 200 K; ____ thermodynamic stability; _ _ _ mechanical stability limits.

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

    \begin{equation} \left\{ \begin{array}{l} T^\sharp = T^\flat = T_ {{{\rm{ref}}}},\\[3pt] p(T^\sharp, v^\sharp, \mathrm{y}^\sharp) = p(T^\flat, v^\flat, \mathrm{y}^\flat),\\[3pt] g_1(T^\sharp, v^\sharp, \mathrm{y}^\sharp) = g_1(T^\flat, v^\flat, \mathrm{y}^\flat),\\[3pt] g_2(T^\sharp, v^\sharp, \mathrm{y}^\sharp) = g_2(T^\flat, v^\flat, \mathrm{y}^\flat),\\[3pt] \mathrm{y}_1^\sharp + \mathrm{y}_2^\sharp = \mathrm{y}_1^\flat + \mathrm{y}_2^\flat = 1. \end{array} \right. \end{equation} (8.5)

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

    Figure 4.  Equilibrium between ethane and nitrogen at T = 220 K; ____ numerical simulation; +++ experimental measurements from Younglove [121] and Younglove and Ely [122].

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

    \begin{equation} \left\{ \begin{array}{l} \partial_{ \mathsf{X}_1} \mathcal{G}(T_ {{{\rm{c}}}}, v_ {{{\rm{c}}}}, \mathrm{y}_ {{{\rm{c}}}}) = 0,\\[3pt] \partial_{ \mathsf{X}_1}^2 \mathcal{G}(T_ {{{\rm{c}}}}, v_ {{{\rm{c}}}}, \mathrm{y}_ {{{\rm{c}}}}) = 0,\\[3pt] \mathrm{y}_{ {{{\rm{c}}}}1} + \mathrm{y}_{ {{{\rm{c}}}}2} = 1. \end{array} \right. \end{equation} (8.6)

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

    Figure 5.  Diagram of critical points for C _2 H _6 -N _2 mixtures in the (T, p) plane; ____ numerical simulation; \ast\ast\ast experimental measurements from Eakin et al. [124]; \triangle\triangle\triangle experimental measurements from Stryjek et al. [125]; +++ experimental measurements from Wisotzki and Schneider [126].
    Figure 6.  Diagram of critical points for C _2 H _6 -N _2 mixtures in the (\mathrm{y}_{{\rm C}_2{\rm H}_6}, p) plane.

    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

    \begin{align} T^ - & = 400\ {\rm K}, \qquad & T^ + & = 250\ {\rm K}, \end{align} (8.7)
    \begin{align} \alpha & = 100\ {\rm s}^{-1}, \qquad & p^\infty & = 3 \ {\rm MPa}. \end{align} (8.8)

    The capillarity coefficient is varied to be

    \varkappa = 5.10^{-4}, \ \; 10^{-3}, \ \; 10^{-2}, \ \; 10^{-1} \ \ {\rm cm}^7 {\rm g}^{-1} {\rm s}^{-2}.

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

    Figure 7.  Mass density of ethane fronts for p^\infty = 3 MPa and various capillarity coefficients.

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

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

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

    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.

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

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

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

    \begin{equation} T d s = d e + p d v - g_1 d \mathrm{y}_1. \end{equation} (A.1)

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

    \begin{equation} v = \mathrm{y}_1 v_1. \end{equation} (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, \mathrm{y}_1) = \mathrm{y}_1 e(T, v_1, 1) , p(T, v, \mathrm{y}_1) = p(T, v_1, 1) , and s(T, v, \mathrm{y}_1) = \mathrm{y}_1 s(T, v_1, 1) . We may thus define the densities e_1(T, v_1) = e(T, v_1, 1) , p_1(T, v_1) = p(T, v_1, 1) , and s_1(T, v_1) = s(T, v_1, 1) , in such a way that

    \begin{align} e( \zeta) = \mathrm{y}_1 e_1(T, v_1), \end{align} (A.3)
    \begin{align} [3pt] p( \zeta) = p_1(T, v_1), \end{align} (A.4)
    \begin{align} [3pt] s( \zeta) = \mathrm{y}_1 s_1(T, v_1). \end{align} (A.5)

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

    g_1 = \partial^{}_{ \mathrm{y}_1}\ \Bigl( \mathrm{y}_1 e(T, v/ \mathrm{y}_1,1) - T \mathrm{y}_1 s(T, v/ \mathrm{y}_1,1) \Bigr) = e_1 - T s_1 + \frac{ \mathrm{y}_1 v}{ \mathrm{y}_1^2} \partial^{}_ v \Bigl( e(T, v_1,1) - T \mathrm{y}_1 s(T, v_1,1) \Bigr),

    so that

    \begin{equation} g_1 = e_1 - T s_1 + v_1 p_1. \end{equation} (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 T \mathrm{d} s = \partial^{}_T e \; \mathrm{d} T + (p + \partial^{}_ v e) \, \mathrm{d} v + (\partial^{}_{ \mathrm{y}_1}\! e - g_1) \, \mathrm{d} \mathrm{y}_1 , we now deduce, after some algebra, that

    \mathrm{y}_1 ( T \mathrm{d} s_1 - \mathrm{d} e_1 + p \mathrm{d} v_1 ) + (T s_1 - e_1 - v_1 p_1 + g_1) \mathrm{d} \mathrm{y}_1 = 0,

    and this yields the usual relation

    \begin{equation} T \mathrm{d} s_1 = \mathrm{d} e_1 + p \mathrm{d} v_1, \end{equation} (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, \mathrm{y}_1) , but the link is easily made.

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

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

    \begin{equation} \widetilde\partial^{}_ e \chi = \partial^{}_T \chi \; \widetilde\partial^{}_ e T, \qquad \partial^{}_T \chi = \widetilde\partial^{}_ e \chi \; \partial^{}_T e, \end{equation} (B.1)
    \begin{equation} \widetilde\partial^{}_ v \chi = \partial^{}_ v \chi + \partial^{}_T \chi \; \widetilde\partial^{}_ v T, \qquad \partial^{}_ v \chi = \widetilde\partial^{}_ v \chi + \widetilde\partial^{}_ e \chi \; \partial^{}_ v e, \end{equation} (B.3)
    \begin{equation} \widetilde\partial^{}_{ \mathrm{y}_k}\! \chi = \partial^{}_{ \mathrm{y}_k}\! \chi + \partial^{}_T \chi \; \widetilde\partial^{}_{ \mathrm{y}_k}\!T, \qquad \partial^{}_{ \mathrm{y}_k}\! \chi = \widetilde\partial^{}_{ \mathrm{y}_k}\! \chi + \widetilde\partial^{}_ e \chi \; \partial^{}_{ \mathrm{y}_k}\! e, \end{equation} (B.3)

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

    \widetilde\partial^2_{ v \mathrm{y}_k}s = \partial^{}_ v\bigl(\frac{- g_k}{T}\bigr) - ( \widetilde\partial^{}_ e T) \frac{( \partial^{}_ v e)( \partial^{}_{ \mathrm{y}_k} e)}{T^2}.

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

    \widetilde\partial^2_{ \mathrm{y}_k \mathrm{y}_l}s = \partial^{}_{ \mathrm{y}_l}\bigl(\frac{- g_k}{T}\bigr) - ( \widetilde\partial^{}_ e T) \frac{( \partial^{}_{ \mathrm{y}_k} e) ( \partial^{}_{ \mathrm{y}_l} e)}{T^2}.

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

    \begin{equation} \bigl\langle ( \widetilde\partial^2_{ \xi \xi} s) \mathsf{X}, \mathsf{X} \bigr\rangle = - \frac{ \widetilde\partial^{}_ e T}{T^2}\, \langle \, \mathsf{\widehat f}^e, \mathsf{X} \rangle^2 + \frac{ \partial^{}_ v p}{T}\, \mathsf{X}_ v^2 - 2 \sum\limits_{k\in {\mathfrak S}} \partial^{}_ v \bigl( \frac{ g_k}{T} \bigr)\, \mathsf{X}_ v \mathsf{X}_k - \sum\limits_{k,l\in {\mathfrak S}} \partial^{}_{ \mathrm{y}_l} \bigl( \frac{ g_k}{T} \bigr)\, \mathsf{X}_k \mathsf{X}_l. \end{equation} (B.4)

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

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

    \begin{equation} \bigl\langle ( \widetilde\partial_{ \xi \xi}^2 s) \mathsf{X}, \mathsf{X} \bigr\rangle = - \widetilde\partial^{}_ e T \frac{\langle \, \mathsf{\widehat f}^e, \mathsf{X} \rangle^2}{T^2} + \partial^{}_ v p \frac{\langle \, \mathsf{\widehat f}^ v, \mathsf{X} \rangle^2}{T} - \sum\limits_{k,l\in {\mathfrak S}} \Bigl( \partial^{}_{ \mathrm{y}_l} \bigl( \frac{ g_k}{T} \bigr)\, + \frac{ \partial^{}_{ v} g_k\, \partial^{}_{ v} g_l}{ T \partial^{}_{ v} p} \Bigr) \mathsf{X}_k \mathsf{X}_l. \end{equation} (B.5)

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

    \begin{equation} \bigl\langle ( \widetilde\partial_{ \xi \xi}^2 s) \mathsf{X}, \mathsf{X} \bigr\rangle = - \widetilde\partial^{}_ e T \frac{\langle \, \mathsf{\widehat f}^ e, \mathsf{X} \rangle^2}{T^2} + \partial^{}_ v p \frac{\langle \, \mathsf{\widehat f}^ v, \mathsf{X} \rangle^2}{T} - \Bigl\langle \bigl(\Lambda - \frac{\Lambda \mathrm{y}{\otimes}\Lambda \mathrm{y}}{ \langle \Lambda \mathrm{y}, \mathrm{y}\rangle} \bigr) \mathsf{X}_ \mathrm{y}^{}, \mathsf{X}_ \mathrm{y}^{} \Bigr\rangle, \end{equation} (B.6)

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

    \begin{equation} \Bigl\langle \bigl( \Lambda - \frac{\Lambda \mathrm{y}{\otimes}\Lambda \mathrm{y}}{ \langle\Lambda \mathrm{y}, \mathrm{y}\rangle} \bigr) \mathsf{X}_ \mathrm{y}, \mathsf{X}_ \mathrm{y}\Bigr\rangle = \Bigl\langle \Lambda \bigl( \mathsf{X}_ \mathrm{y} - \tfrac{\langle \Lambda \mathsf{X}_ \mathrm{y}, \mathrm{y} \rangle}{ \langle\Lambda \mathrm{y}, \mathrm{y}\rangle} \mathrm{y}\bigr), \; \bigl( \mathsf{X}_ \mathrm{y} - \tfrac{\langle \Lambda \mathsf{X}_ \mathrm{y}, \mathrm{y} \rangle}{ \langle\Lambda \mathrm{y}, \mathrm{y}\rangle} \mathrm{y} \bigr) \Bigr\rangle. \end{equation} (B.7)

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

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

    \bigl\langle \Lambda \mathsf{X}_ \mathrm{y}, \mathsf{X}_ \mathrm{y} \rangle = \frac{1}{\langle\Lambda \mathrm{y}, \mathrm{y}\rangle} \langle\Lambda \mathsf{X}_ \mathrm{y}, \mathrm{y}\rangle^2 + \bigl\langle \widehat\Lambda \mathsf{X}_ \mathrm{y}, \mathsf{X}_ \mathrm{y} \bigr\rangle,

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

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

    \begin{align*} \boldsymbol{\mathcal{P}} & = p^ {\rm{ex}} \boldsymbol{I} + \sum\limits_{i,j\in \mathfrak{S}} \varkappa_{ij} \boldsymbol{\nabla}\mspace{-3mu}\rho_i {\otimes} \boldsymbol{\nabla}\mspace{-3mu}\rho_j - \sum\limits_{i,j\in \mathfrak{S}} \rho_i \boldsymbol{\nabla} {\boldsymbol{\cdot}}( \varkappa_{ij} \boldsymbol{\nabla}\mspace{-3mu}\rho_j) + \boldsymbol{\mathcal{P}}^{{\rm{d}}}, \\ \boldsymbol{\mathcal{Q}} & = \sum\limits_{i,j\in \mathfrak{S}} \varkappa_{ij} \boldsymbol{\nabla}\mspace{-3mu}\rho_j \bigl( \rho_i \boldsymbol{\nabla} {\boldsymbol{\cdot}} {\boldsymbol{u}} + \boldsymbol{\nabla} {\boldsymbol{\cdot}} \boldsymbol{\mathcal{J}}_i \bigr) - \sum\limits_{i,j\in \mathfrak{S}} \boldsymbol{\nabla} {\boldsymbol{\cdot}}( \varkappa_{ij} \boldsymbol{\nabla}\mspace{-3mu}\rho_j) \boldsymbol{\mathcal{J}}_i + \boldsymbol{q}, \end{align*}

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

    \begin{align*} \boldsymbol{\mathcal{P}}^{{\rm{d}}} = & - \mathfrak{v} \boldsymbol{\nabla} {\boldsymbol{\cdot}} {\boldsymbol{u}} \, \boldsymbol{I} - \eta \bigl( \boldsymbol{\nabla} {\boldsymbol{u}} + \boldsymbol{\nabla} {\boldsymbol{u}}^t - {\tfrac{2}{3}} \boldsymbol{\nabla} {\boldsymbol{\cdot}} {\boldsymbol{u}} \, \boldsymbol{I} \bigr),\\[3pt] \boldsymbol{\mathcal{J}}_i = & - \sum\limits_{j\in \mathfrak{S}} L_{ij} \Bigl( \boldsymbol{\nabla} \Bigl( \frac{ g_j^ {\rm{ex}}}{T} \Bigr) - \frac{ {\boldsymbol\nabla}\, {\boldsymbol\nabla} {\boldsymbol{\cdot}} (\sum\nolimits_{l\in \mathfrak{S}} \varkappa_{jl} \boldsymbol{\nabla}\mspace{-3mu}\rho_l)}{T} \Bigr) - L_{i \mathrm{e}} {\boldsymbol\nabla}\Bigl( \frac{-1}{T} \Bigr),\\ \boldsymbol{q} = & - \sum\limits_{i\in \mathfrak{S}} L_{ \mathrm{e} i} \Bigl( {\boldsymbol\nabla} \Bigl(\frac{ g_i^ {\rm{ex}}}{T} \Bigr) - \frac{ {\boldsymbol\nabla}\, {\boldsymbol\nabla} {\boldsymbol{\cdot}} (\sum\nolimits_{l\in \mathfrak{S}} \varkappa_{il} \boldsymbol{\nabla}\mspace{-3mu}\rho_l)}{T} \Bigr) - L_{ \mathrm{e} \mathrm{e}} {\boldsymbol\nabla}\Bigl( \frac{-1}{T} \Bigr), \end{align*}

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

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

    \begin{align*} & \mathcal{E}^ {\rm{ex}} = \mathcal{E}^ {} + {\tfrac{1}{2}} \sum\limits_{i,j\in \mathfrak{S}} ( \varkappa_{ij}-T \partial_T \varkappa_{ij}) \boldsymbol{\nabla}\mspace{-3mu}\rho_i {\boldsymbol{\cdot}} \boldsymbol{\nabla}\mspace{-3mu}\rho_j, \qquad p^ {\rm{ex}} = p^ {} - {\tfrac{1}{2}} \sum\limits_{i,j\in \mathfrak{S}} \varkappa_{ij} \boldsymbol{\nabla}\mspace{-3mu}\rho_i {\boldsymbol{\cdot}} \boldsymbol{\nabla}\mspace{-3mu}\rho_j, \\[3pt] & \mathcal{S}^ {\rm{ex}} = \mathcal{S}^ {} - {\tfrac{1}{2}} \sum\limits_{i,j\in \mathfrak{S}} \partial_T \varkappa_{ij} \boldsymbol{\nabla}\mspace{-3mu}\rho_i {\boldsymbol{\cdot}} \boldsymbol{\nabla}\mspace{-3mu}\rho_j, \qquad\qquad g_i^ {\rm{ex}}\ = g_i^ {}, \quad i\in \mathfrak{S}, \end{align*}

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

    \begin{equation} T \mathrm{d} \mathcal{S}^ {\rm{ex}} = \mathrm{d} \mathcal{E}^ {\rm{ex}} - \sum\limits_{i\in \mathfrak{S}} g_i^ {\rm{ex}} \mathrm{d} \rho_i - \sum\limits_{i,j\in \mathfrak{S}} \varkappa_{ij} {\boldsymbol\nabla}\rho_i {\boldsymbol{\cdot}} \mathrm{d} {\boldsymbol\nabla}\rho_j. \end{equation} (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:

    \begin{align*} \varkappa_{ij} = \varkappa(T) \qquad \quad i,j \in \mathfrak{S}. \end{align*}

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

    \begin{align} & \boldsymbol{\mathcal{P}} = p^ {\rm{ex}} \boldsymbol{I} + \varkappa \boldsymbol{\nabla}\mspace{-3mu}\rho {\otimes} \boldsymbol{\nabla}\mspace{-3mu}\rho - \rho \boldsymbol{\nabla} {\boldsymbol{\cdot}}( \varkappa \boldsymbol{\nabla}\mspace{-3mu}\rho) \boldsymbol{I} + \boldsymbol{\mathcal{P}}^{{\rm{d}}}, \qquad \boldsymbol{\mathcal{Q}} = \varkappa \rho \boldsymbol{\nabla}\mspace{-3mu}\rho \boldsymbol{\nabla} {\boldsymbol{\cdot}} {\boldsymbol{u}} + \boldsymbol{q}, \end{align} (C.2)
    \begin{align} [3pt] & \boldsymbol{\mathcal{J}}_i = - \sum\limits_{j\in \mathfrak{S}} L_{ij} \boldsymbol{\nabla} \Bigl( \frac{ g_j}{T} \Bigr) - L_{i \mathrm{e}} {\boldsymbol\nabla}\Bigl( \frac{-1}{T} \Bigr), \qquad\quad\ \boldsymbol{q} = - \sum\limits_{i\in \mathfrak{S}} L_{ \mathrm{e} i} {\boldsymbol\nabla} \Bigl( \frac{ g_i}{T} \Bigr) - L_{ \mathrm{e} \mathrm{e}} {\boldsymbol\nabla}\Bigl( \frac{-1}{T} \Bigr). \end{align} (C.3)

    The thermodynamic properties also read

    \begin{align} & \mathcal{E}^ {\rm{ex}} = \mathcal{E}^ {}(\rho_1,\ldots,\rho_ n,T) + {\tfrac{1}{2}} ( \varkappa-T \partial_T \varkappa) \vert \boldsymbol{\nabla}\mspace{-3mu}\rho\vert^2, \qquad p^ {\rm{ex}} = p^ {}(\rho_1,\ldots,\rho_ n,T) - {\tfrac{1}{2}} \varkappa\vert \boldsymbol{\nabla}\mspace{-3mu}\rho\vert^2, \end{align} (C.4)
    \begin{align} [3pt] & \mathcal{S}^ {\rm{ex}} = \mathcal{S}^ {}(\rho_1,\ldots,\rho_ n,T) - {\tfrac{1}{2}} \partial_T \varkappa \vert \boldsymbol{\nabla}\mspace{-3mu}\rho\vert^2, \qquad\ g_i^ {\rm{ex}}\ = g_i^ {}(\rho_1,\ldots,\rho_ n,T), \qquad i\in \mathfrak{S}, \end{align} (C.5)

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



    [1] Alles L, Murray L (2017) Asset pricing and downside risk in the Australian share market. Appl Econ 49: 4336–4350. http://dx.doi.org/10.1080/00036846.2017.1282143 doi: 10.1080/00036846.2017.1282143
    [2] Basher SA, Sadorsky P (2006) Oil price risk and emerging stock markets. Glob Financ J 17: 224–251. http://dx.doi.org/10.1016/j.gfj.2006.04.001 doi: 10.1016/j.gfj.2006.04.001
    [3] Bayes T (1763) An Essay towards Solving a Problem in the Doctrine of Chances. by the Late Rev. Mr. Bayes, F. R. S. Communicated by Mr. Price, in a Letter to John Canton, A. M. F. R. S. Phil Trans R Soc 53: 370–418. https://doi.org/10.1098/rstl.1763.0053 doi: 10.1098/rstl.1763.0053
    [4] Brooks C (2014) Introductory Econometrics for Finance. Cambridge University Press: New York.
    [5] Devore JL (2012) Probability and Statistics for Engineering and the Sciences. California Polytechnic State University: San Luis Obispo.
    [6] Emmert-Streib F, Dehmer M (2019) Introduction to Survival Analysis in Practice. Mach Learn Know Extr 1: 1013–1038. https://doi.org/10.3390/make1030058 doi: 10.3390/make1030058
    [7] Fassas AP, Siriopoulos C (2020) Implied volatility indices—A review. Q Rev Econ Financ 79: 303–329. https://doi.org/10.1016/j.qref.2020.07.004 doi: 10.1016/j.qref.2020.07.004
    [8] Floros C, Gkillas K, Konstantatos C, Tsagkanos A (2020) Realized Measures to Explain Volatility Changes over Time. J Risk Financ Manag 13: 1–19. https://doi.org/10.3390/jrfm13060125 doi: 10.3390/jrfm13060125
    [9] Gao G, Bu Z, Liu L, et al. (2019) A Survival Analysis Method for Stock Market Prediction. 2015 International Conference on Behavioral, Economic and Socio-cultural Computing (BESC), 116–122. https://doi.org/10.1109/BESC.2015.7365968 doi: 10.1109/BESC.2015.7365968
    [10] Gepp A, Kumar, K (2015) Predicting Financial Distress: A Comparison of Survival Analysis and Decision Tree Technique. Procedia Comput Sci 54: 396–404. https://doi.org/10.1016/j.procs.2015.06.046 doi: 10.1016/j.procs.2015.06.046
    [11] Griffin JE, Kalli M, Steel M (2018) Discussion of "Nonparametric Bayesian Inference in Applications": Bayesian nonparametric methods in econometrics. Stat Method Appl 27: 207–218. https://doi.org/10.1007/s10260-017-0384-0 doi: 10.1007/s10260-017-0384-0
    [12] Hansen PR, Lunde A (2006) Realized variance and market microstructure noise. J Bus Econ Stat 24: 127–161. https://doi.org/10.1198/073500106000000071 doi: 10.1198/073500106000000071
    [13] Hoque ME, Low S (2020) Industry Risk Factors and Stock Returns of Malaysian Oil and Gas Industry: A New Look with Mean Semi-Variance Asset Pricing Framework. Mathematics 8: 1–28. https://doi.org/10.3390/math8101732 doi: 10.3390/math8101732
    [14] Jensen MJ, Maheu JM (2018) Risk, Return and Volatility Feedback: A Bayesian Nonparametric Analysis. J Risk Financ Manag 11: 1–29. https://doi.org/10.3390/jrfm11030052 doi: 10.3390/jrfm11030052
    [15] Kalli M, Griffin J, Walker S (2011) Slice sampling mixture models. Stat Comput 21: 93–105. http://dx.doi.org/10.1007/s11222-009-9150-y doi: 10.1007/s11222-009-9150-y
    [16] Karabatsos G (2017) A menu-driven software package of Bayesian nonparametric (and parametric) mixed models for regression analysis and density estimation. Behav Res Methods 49: 335–362. https://doi.org/10.3758/s13428-016-0711-7 doi: 10.3758/s13428-016-0711-7
    [17] Markowitz H (1952) Portfolio selection. T J Financ 7: 77–91. https://doi.org/10.2307/2975974 doi: 10.2307/2975974
    [18] Moghaddam MD, Liu J, Serota RA (2020) Implied and realized volatility: A study of distributions and the distribution of difference. Int J Financ Econ 26: 2581–2594. https://doi.org/10.1002/ijfe.1922 doi: 10.1002/ijfe.1922
    [19] National Treasury (2021) Economic Overview. In: Chapter 2. Budget Review. Republic of South Africa. Available from: http://www.treasury.gov.za/documents/National%20Budget/2021/review/Chapter%202.pdf.
    [20] Papaspiliopoulos O (2008) A note on posterior sampling from Dirichlet mixture models. Available from: http://www2.warwick.ac.uk/fac/sci/statistics/crism/research/2008/08-20wv2.pdf.
    [21] Sehgal S, Pandey A (2018) Predicting Financial Crisis by Examining Risk-Return Relationship. Theor Econ Lett 8: 48–71. https://doi.org/10.4236/tel.2018.81003 doi: 10.4236/tel.2018.81003
    [22] Sharma S (2017) Markov Chain Monte Carlo Methods for Bayesian Data Analysis in Astronomy. Annu Rev Astron Astrophys 55: 213–259. https://doi.org/10.48550/arXiv.1706.01629 doi: 10.48550/arXiv.1706.01629
    [23] Stankovic JZ, Petrocvic E, Dencic-Mihajlov K (2020) Effects of Applying Different Risk Measures on the Optimal Portfolio Selection: The Case of the Belgrade Stock Exchange. Facta Universitatis Series: Economics and Organization 17: 17–26. https://doi.org/10.22190/FUEO191016002S doi: 10.22190/FUEO191016002S
    [24] Steyn JP, Theart L (2019) Are South African equity investors rewarded for taking on more risk? J Econ Financ Sci 12: 1–10. https://doi.org/10.4102/jef.v12i1.448 doi: 10.4102/jef.v12i1.448
    [25] Trichilli Y, Abbes MB, Masmoudi A (2020) Islamic and conventional portfolios optimization under investor sentiment states: Bayesian vs Markowitz portfolio analysis. Res Int Bus Financ 51: 1–21. https://doi.org/10.1016/j.ribaf.2019.101071 doi: 10.1016/j.ribaf.2019.101071
    [26] van Ravenzwaaij D, Cassey P, Brown SD (2018) A simple introduction to Markov Chain Monte-Carlo sampling. Psychon Bull Rev 25: 143–154. https://doi.org/10.3758/s13423-016-1015-8 doi: 10.3758/s13423-016-1015-8
    [27] Walker SG (2007) Sampling the dirichlet mixture model with slices. Commun Stat-Simul C 36: 45–54. http://dx.doi.org/10.1080/03610910601096262 doi: 10.1080/03610910601096262
    [28] Yildiz ME, Erzurumlu YO, Kurtulus B (2022) Comparative analyses of mean-variance and mean-semivariance approaches on global and local single factor market model for developed and emerging markets. Int J Emerg Mark 17: 325–350. https://doi.org/10.1108/IJOEM-01-2020-0110 doi: 10.1108/IJOEM-01-2020-0110
  • Reader Comments
  • © 2022 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(1463) PDF downloads(108) Cited by(0)

Figures and Tables

Figures(38)  /  Tables(1)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog