AN ANALYSIS OF THE UNIFIED FORMULATION FOR THE EQUILIBRIUM PROBLEM OF COMPOSITIONAL MULTIPHASE MIXTURES

. In this paper, we conduct a thorough mathematical analysis of the unified formulation advocated by Lauser et al . [ Adv. Water Res. 34 (2011) 957–966] for compositional multiphase flows in porous media. The interest of this formulation lies in its potential to automatically handle the appearance and disappearance of phases. However, its practical implementation turned out to be not always robust for realistic fugacity laws associated with cubic equations of state, as shown by Ben Gharbia and Flauraud [ Oil Gas Sci. Technol. 74 (2019) 43]. By focusing on the subproblem of phase equilibrium, we derive sufficient conditions for the existence of the corresponding system of equations. We trace back the difficulty of cubic laws to a deficiency of the Gibbs functions that comes into play due to the “unifying” feature of the new formulation. We propose a partial remedy for this problem by extending the domain of definition of these functions in a natural way. Besides, we highlight the crucial but seemingly unknown fact that the unified formulation encapsulates all the properties known to physicists on phase equilibrium, such as the tangent plane criterion and the minimization of the Gibbs energy of the mixture.


Motivation and objectives
In the numerical simulation of multicomponent (a.k.a.compositional) multiphase fluid flows, a delicate issue often arises in the handling of the appearance and disappearance of phases for various species, due to the laws of thermodynamic equilibrium.The traditional dynamic approach, known as variable-switching in reservoir simulations [10], considers only the unknowns and equations of the present phases.Albeit natural, it is awkward and even costly to implement, insofar as switching can occur all the time.Lauser et al. [17] proposed an alternative approach, called unified formulation, in which a fixed set of unknowns and equations is maintained during the calculations.This major theoretical advance is achieved by means of complementarity conditions, which allow distinct functioning regimes to be expressed in the same mathematical way, as is already the case in a wide range of areas such as mechanics, electronics or geology [1,31].Another key ingredient behind this "egalitarian" treatment of all regimes is the notion of extended partial fractions that must be assigned to species in all phases, including absent ones.
As a promise of more efficient simulations, the unified formulation has met with some success among numericists, as testified by the subsequent works by Ben Gharbia [5], Ben Gharbia and Jaffr [7], Masson et al. [20,21] and Beaude et al. [4].These are all based on simple fugacity coefficients, such as given by Henry's law.Another series of works at IFPEN [6,8,19,30] is focused on realistic fugacity coefficients given by cubic equations of state, such as Peng-Robinson's law.Although the latter investigations have demonstrated a clear superiority of the unified formulation over the variable-switching one regarding computational time in some cases, the outcome remains unclear in other cases with single-phase transition: the nonlinear solver for the (unified) algebraic system of equations may not converge at all, unlike its competitor.
There are two possible explanations for this observed lack of robustness from the unified formulation.To sketch them out in a precise manner, we need the following formal setup.After discretization in time and space of the continuous flow model using the unified formulation, the system of equations to be solved at each time-step takes the abstract form Λpq " 0, (1.1a) minppq, pqq " 0, (1.1b) where  P  Ă R ℓ is the unknown vector and Λ :  Ñ R ℓ´ ,  :  Ñ R  et  :  Ñ R  are continuously differentiable functions on the open domain .The componentwise action of the minimum function in (1.1b) is merely a convenient way of expressing the complementarity 0 ď pq K pq ě 0. For conciseness, let us put  pq " " Λpq minppq, pqq so that (1.1) becomes  pq " 0. We can then envision two scenarios that could cause the unified formulation to perform poorly: (1) System (1.1) is ill-posed for some data and thermodynamic laws.In other words, it may not have a unique solution or may not have a solution at all.An even worse situation is when some components of Λ -and therefore of  -are not well-defined over the whole domain of interest , so that (1.1) no longer makes sense.As will be seen later, this occurs for cubic equations of state frequently used in realistic simulations.(2) The numerical algorithm used to solve system (1.1) is not well suited to the semismooth nature of  .Indeed, the complementarity equations (1.1b) are not differentiable, which prevents the standard Newton method to be applied.A common remedy is the so-called Newton-min method [2,15].However, Newton-min may suffer from periodic oscillations for large time-steps, as evidenced by Ben Gharbia and Flauraud [6].
The first issue originates from physical modeling.It is the subject of this article, whose primary objective is to clarify the conditions on thermodynamic laws under which system (1.1) is well-behaved and to propose some improvements of the model so as to guarantee the existence of a solution.The second issue pertains to numerical methods.It requires a new method to be designed in order to replace Newton-min and was addressed in a previous paper [35].

Main results and outline of the paper
Valuable insights into the difficulty can be gained if, instead of the fully discretized flow model, we focus on an elementary phase equilibrium problem that lies at the core of the thermodynamic part.This is why we start by stating the phase equilibrium problem for multicomponent mixtures in Section 2, comparing the variable-switching formulation to the unified formulation (Sect.2.2) after recalling some preliminary notions in Section 2.1.Section 3 is devoted to the analysis of the unified formulation.We first revisit two thermodynamic properties in light of the new framework, namely, the principle of Gibbs energy minimization (Sect.3.1) and the tangent plane criterion (Sect.3.2).Although these properties are well-known in thermodynamics by virtue of various physical arguments, the point we would like to make here is that they are all mathematical consequences of the unified formulation.In (3.3), we introduce an important phasewise subproblem called local inversion of extended fugacities in Section 3.3.1.Sufficient conditions are worked out to ensure the existence and uniqueness of a solution to this extended fugacities inversion subproblem.In essence, we require strict convexity of the Gibbs functions in each phase, as well as invertibility and surjectivity of their gradient maps.These assumptions are also shown in Section 3.4 to guarantee the existence of a solution to the full phase equilibrium problem.
Given a fugacity or activity law from physics textbooks, there is no reason for the corresponding molar Gibbs function to fulfill the hypotheses of strict convexity and invertibility/surjectivity of the gradient map.In Section 4, we further investigate the question of strict convexity for some simple fugacity and activity models, namely, Henry's laws (Sect.4.1), Margules' law (Sect.4.2), and Van Laar's law (Sect.4.3).For each of these, we manage to determine the subregion in the space of parameters for which strict convexity holds.
A prominent category of fugacity laws widely used in realistic simulations of two-phase mixtures stems from cubic equations of state (EOS).As recalled in Section 5.1, the definition of the corresponding thermodynamic quantities involves solving a cubic equation which does not always have three real roots.After a careful study of the critical values (Sect.5.2) and the frontier between the 1-root and 3-root regions (Sect.5.3) for Peng-Robinson's law, one of the most advanced cubic EOS-based models, we explain the trouble with these laws regarding the domains of definition for different functions involved in (1.1).In a nutshell, since there are not always three real roots, the Gibbs functions and fugacity coefficients are not always well-defined simultaneously for both phases over the whole domain of generalized partial fractions.While this pathology is not detrimental to the variable-switching formulation, where only present phases are considered, it causes tremendous harm to the unified formulation, for which information relative to both phases must be permanently available.The uncovering of this difficulty in Section 6.1 prompts us to design an extension procedure for various thermodynamic functions, in an attempt to maintain a good behavior for the unified formulation, that is, to hope for the existence of a solution to (1.1).The basic idea, elaborated on in Section 6.2, is to extend the Gibbs functions by replacing the missing real root by the common real part of the two conjugate complex roots.This construction is supported by further calculations.

Preliminary notions
In this paper, we are concerned with the advantages and drawbacks of the unified formulation for the phase equilibrium problem at fixed pressure and temperature.To state this problem, one needs some prerequisites on the thermodynamics of multiphase multicomponent mixtures.

Species, phases and fractions
A multicomponent mixture is a physical system consisting of several chemically distinct components or species, e.g., hydrogen pH 2 q, water pH 2 Oq, carbon dioxide pCO 2 q, methane pCH 4 q . . .It can be thought of in a more abstract way by introducing the set of species  " tI, II, . . ., u,  ě 2, ( whose elements are labeled by Roman numerals.The total number of components  " || usually ranges from tens to hundreds.Each component  P  may be present under one or many phases.Intuitively, a phase is more or less a state of matter, e.g., gas pq, liquid pq, oil pq, solid pq . . .This notion may be subtler, though, at high pressure [12].Again, to lay down an abstract framework, we consider the set of all virtually possible phases  Although this choice somehow breaks the symmetry, it is commonly resorted to in practice.Finally, there is a third notion of fraction, called global fractions and denoted by   P r0, 1s, which quantifies the overall relative importance of each component  P  inside the mixture.Of course, we have (2.7) The vector  " `I , . . .,  ´1 ˘P Ω Ă R ´1 is called global composition of components.Again, because of the dependence (2.7),only the first  ´1 values in .Whenever a   appears in the text, it should be understood as   " 1 ´I ´. . .´´1 .The material balance of component  implies that Given the context Γ, the phasic fractions t  u PΓ and the partial fractions t   u p,qPˆP , it is straightforward to calculate the global composition t  u P by (2.8).The phase equilibrium problem takes exactly the opposite direction: given the global composition t  u P satisfying (2.7), is it possible to find the context Γ, the phasic fractions t  u PΓ and the partial fractions t   u p,qPˆP satisfying (2.3), (2.5) and (2.8) beside positivity?Obviously, we do not have enough equations yet.The missing ones will be supplied at the end of Section 2.1.3.
Remark 2.1.We have deliberately not specified whether the three kinds of fractions   ,    and   are molar, volumic or specific fractions.In fact, this does not matter.The mathematical structure of the problem remains the same and the theoretical development is similar in all cases.

Gibbs energy and chemical potential
The behavior of each phase  P P is governed by a single fundamental function   : Ω Ñ R known as the (intensive) Gibbs free energy of the phase.We require   to be as smooth as necessary in Ω and continuously extendable to BΩ.However, ∇  may blow up on BΩ.From   , we define  functions    : Ω Ñ R,  P , called chemical potentials by    pq "   pq `@∇  pq,   ´D (2.9) for  P Ω, where the vector   " p ,1 ,  ,2 , . . .,  ,´1 q P R ´1 is made up of Kronecker's symbols.The following statement gives two helpful identities between   and    .The first one (2.10a)relates the Gibbs energy to the potentials.The second one (2.10b)provides the gradient of the Gibbs energy from the potentials.(2.9).This cancels out   pq and the desired identity follows from   " p0, 0, . . ., 0q.

Lemma 2.2 (Connection between Gibbs energy and chemical potentials
Remark 2.3.In the above, we used the generic variable  to alleviate notations.Of course,   is to be evaluated at   , the composition of phase .As a matter of fact, the Gibbs function also depends on the pressure P  and the temperature T  of the phase [33].But since we work at fixed pressure and temperature, we purposely omit to write them down in order to concentrate on the dependency with respect to the fractions.

Fugacity, fugacity coefficient and quilibrium conditions
For a solid phase,    is a constant.For fluid phases such as gas, liquid and oil, the chemical potentials take the form in which Φ   is called the fugacity coefficient of component  in phase .Note, however, that it depends on the whole composition vector.As for the quantity it is known as the fugacity of component  in phase .Substituting the form (2.11a) into (2.10a),we obtain 12) The first sum ř  "I   ln   is the ideal part.The second sum, denoted by is the excess part.In this perspective, a fluid phase  is assimilated to a "perturbation" of the ideal gas.As will be done in Section 6, we shall act only on the excess part to modify the Gibbs function.
Owing to the regularity assumptions made on   and    , the functions Ψ  : Ω Ñ R and ln Φ   : Ω Ñ R are also as smooth as necessary, with Ψ  extendable by continuity to Ω but not the ln Φ   's.The very useful relations below between Ψ  and ln Φ   are similar to those between   and    .
Lemma 2.4 (Connection between excess energy and fugacity coefficients).For all  P Ω: ln Φ   pq " Ψ  pq `@∇Ψ  pq,   ´D , @ P ; (2.14a) pq " ln Φ   pq ´ln Φ   pq, @ P ztu; (2.14b) Proof.The proof is straightforward.For each identity from Lemma 2.2, we just have to separate the ideal part from the excess part.The ideal part vanishes trivially.
In a multicomponent mixture without any chemical reaction (also called non-reactive), the presence of two phases p, q P Γ ˆΓ implies that some equilibrium conditions must be achieved.According to thermodynamics, these conditions are the equalities across the two phases of pressure, temperature, and the chemical potentials corresponding to each component  P .In other words, the missing conditions for the phase equilibrium problem at fixed pressure and temperature are    p  q "    p  q, for all p, , q P  ˆΓ ˆΓ, or equivalently,    Φ   p  q "    Φ   p  q, for all p, , q P  ˆΓ ˆΓ. (2.15b) The fugacity coefficients Φ which differs from   by an additive affine function and a multiplicative constant.
A given family of positive real-valued functions tΦ   u p,qPˆP is said to be admissible if, for each  P P, there exists a Gibbs energy function   of which they are the fugacity coefficients.

Two mathematical formulations
Equipped with the preliminary notions of Section 1, we are now in a position to rigorously state the phase equilibrium problem in two different ways: the "traditional" one and the "modern" one.

Variable-switching formulation
The first formulation has the advantage of being "natural," insofar as it uses the variables that have been introduced so far.It also bears the name of natural variable formulation.

GIVEN
, P, tΦ   u p,qPˆP admissible,  P Ω, FIND Γ Ă P, t  u PΓ ą 0, t   u p,qPˆΓ ě 0 so as to satisfy ÿ PΓ      ´ " 0, @ P ; (2.17a) where  is a fixed phase of Γ. Obviously, equation (2.17b) is none other than (2.15b), but expressed in such a way to avoid redundancy.The material balances (2.17a), (2.17c) respectively match (2.8), (2.5).Note that (2.3) is not explicitly prescribed because it can be deduced from the existing equations by summing (2.17a) over  P , switching order, and invoking (2.17c).For a given context Γ, system (2.17) contains p `1q|Γ| equations and unknowns.It must of course be assumed that the physical properties of the species involved are such that the p|Γ| ´1q fugacity equalities (2.17b) are independent.
The price to be paid for naturality is that the context Γ is itself an unknown.To circumvent this difficulty, we start by making an "educated guess" for Γ.At every fixed Γ, we attempt to solve the algebraic equations (2.17): this is what physicists call a pP, Tq-flash.After exiting the flash, we check the positivity of   and the non-negativity of    , for  P Γ.Should one of these fractions have the wrong sign, we must change Γ by adding or deleting phases and go for another flash!The number of unknowns and equations for a flash strongly depends on the assumption currently made about the context Γ.There is a vast literature on numerical methods [22][23][24]36] for the flash problem (2.17) at fixed Γ.In addition to the classical and generic Newton-Raphson method [3,33], many special purpose algorithms have been dedicated to the flash problem.These are iterative methods based on various kinds of substitution [13], the most famous of them being the Rachford-Rice substitution [32].

Unified formulation
To avoid the annoyance of dynamically handling the context, Lauser et al. [17]  As a consequence, for each phase  P P, there are three possible regimes: -  ą 0 (phase  is present).This implies ř P    " 1 and by (2.18d),    "    .Hence, the extended fractions of a present phase coincide with the usual partial fractions.
-1 ´řP    ą 0. This entails   " 0 (phase  is absent) and    ‰    .The extended fractions of an absent phase differ from the usual partial fractions (see exception below).
-  " 0 and 1 ´řP    " 0. This corresponds to a transition point, where phase  starts appearing or disappearing.
It is legitimate to wonder about the origin of the sign condition 1 ´řP    ě 0. After all, it brings a new piece of information that was not included in the variable-switching formulation (2.17).As will be proven in Section 3.2, this condition ensures a stability property known as the tangent plane criterion by physicists.It can also be related to the minimization of the Gibbs energy of the mixture, as will be done in Section 3.1.
The ability of formulation (2.18) to deal with all possible configurations (arising from the presence or the absence of each phase) in a unified way is very attractive not only for convenience but also for computational efficiency.The context Γ no longer appears in the statement of the problem, but can be determined a posteriori by collecting those phases  for which   ą 0. As before, note that the phase balance (2.3) is not explicitly imposed because it can be recovered from the existing equations by summing (2.18a) over  P , permuting order and taking advantage of (2.18c).System (2.18) has p `1q equations and unknowns.It can be cast under the abstract form (1.1) with The existence of a solution to (2.18) can be guaranteed under some sufficient conditions on the Gibbs functions   , as elucidated in Section 3.4.

Properties of the unified formulation
The unified formulation enjoys many remarkable properties that seem to be unknown so far, at least to our knowledge.In particular, by postulating 1 ´řP    ě 0 from the beginning, it achieves a deep connection with some classical results in thermodynamics.

Connection with Gibbs energy minimization
We would like to better understand where this sign information comes from.In the literature, the condition 1 ´řP    ě 0 is customarily derived from a phase stability analysis [22] (see also [6] for a more recent presentation).However, this classical analysis suffers from a few limitations.First, it is restricted to two phases.Second, it is local: the Gibbs energy difference under study must be linearized via a first-order Taylor expansion, before minimizing.Third, the notion of extended fractions appears only at the end, in a very ad hoc way.
We propose a more direct connection between the unified formulation (2.18) and some Gibbs energy minimization problem expressed in terms of the extended fractions    , without any linearization.In this problem, the quantities 1 ´řP    will appear to be the Lagrange multipliers associated with the constraints   ě 0. Conversely, while not every critical point of the minimization problem pq is a solution of the unified formulation, some "natural" choice of critical points satisfies the unified formulation.

Towards a novel interpretation
In order to state the minimization problem, we need to introduce a new Gibbs function.For each phase  P P, let g  : R  `Ñ R be the extended molar Gibbs energy defined as using the renormalization (2.18d) to compute  P Ω from  " `I , . . .,   ˘P R  `zt0u.For normalized fractions, g  p I , . . .,   q "   pq.Thus,   lifts the intensive Gibbs function   to the domain of extended fractions, but it does not coincide with the usual extensive Gibbs function [33].The following Lemma summarizes its most useful properties.Lemma 3.1.For  P R  `zt0u and  P , we have ) Proof.The readers are referred to Lemma 2.3 from [34].The calculations involve the extensive Gibbs energy that we have not introduced here for conciseness, but are not difficult.
We can now consider the following minimization problem pq.
The last equation (3.4e) expresses the complementarity between each inequality constraint (3.3d) and its Lagrange multiplier at optimality.It can be rephrased as A set of values tp  ,   qu PP is said to be a critical point for problem (3.3) if there exists a set of values p, t  u P , t  u PP q such that the KKT optimality system (3.4) is satisfied.

From one formulation to the other
We first show that it is easy to go from the unified formulation to the minimization problem.
Theorem 3.2.Every solution `s   , s   ˘(PP of the unified formulation (2.18) is a critical point of the minimization problem (3.3), with where s   is the common value of the extended fugacity s    Φ   ps   q across all phases  P P.
Proof.Let `s   , s   ˘(PP be a solution of (2.18).The material balances (3.4c), (3.4d) are naturally met, as observed in Section 2.2.2.The equality of extended fugacities (2.18b) makes it possible to define s   " ´"ln ` ˘`1 ‰ in the way described in the theorem.This choice of s   trivially fulfills (3.4b) because of (3.2a).The choice of s   implies (3.4e) because of (2.18c).It remains to check (3.4a).To this end, we use Lemma 3.1 to write This completes the proof.
The reverse direction is more delicate.The main difficulty lies in the indetermination of the extended fractions for an absent phase.(1) If two phases p, q P P ˆP are both present, i.e., r   ą 0 and r   ą 0, then This implies that the complementarity condition (2.18c) holds for both phases and that the extended fugacity equalities (2.18b) hold between the two phases considered.
(2) If phase  is present and phase  is absent, i.e., r   ą 0 and r   " 0, then This is none other than the second part of (3.7).
To fully grasp the meaning of Theorem 3.3, it is capital to observe that when a critical point of (3.3) has a vanishing phase  P P for which r   " 0, the corresponding extended fractions r   cannot be uniquely determined.Indeed, r   plainly does not contribute to neither the objective function (3.3a) nor the constraint (3.3c) at fixed r   " 0. To put it another way, changing r   to any other vector R  `will provide another acceptable critical point.Thus, as soon as there is a critical point of (3.3) for which r   " 0, there are in fact an infinity of such critical points.Among this infinity of critical points, only those for which r    Φ   pr   q " r    Φ   pr   q for all  P , where  is present phase ´r   ą 0 ¯, will be also solutions of the unified formulation (2.18).Combining this with Theorem 3.2, we can interpret the unified formulation as a set of equations that is slightly "stronger" than that of the KKT system for the critical points.It is stronger in the sense that it helps selecting some special critical points -and hopefully just one -among the infinity of possible critical points that appear when one of the phases disappears.

A continuity principle
We even have heuristic arguments to claim the critical points selected by the unified formulation are "natural" ones.By this, we mean that the additional conditions (3.8) to be prescribed on the extended fractions of an absent phase  can be construed as the limit of a continuous process during which  was present before vanishing.To build up this process, let us reformulate the minimization problem pq or (3.3) as the bilevel or hierarchical problem min

Tangent plane criterion
Another set of properties can be established by looking at the geometric significance of the extended fugacity equalities (2.18b).Recall that Ω defined in (2.6b), is the domain of the partial fractions  renormalized from  by (2.18d B  B  ps   q for all  P tI, II, . . .,  ´1u.This completes the proof for (3.13c).
The first part of Theorem 3.4 reveals that, in general, there is no equality of chemical potentials if these were computed using the renormalized partial fractions.The second part of Theorem 3.4 is more interesting.Let us investigate this further by making an additional assumption on one of the phases.We recall the definition (3.12) of the linearized expansion     pq. ps   q ě    ps   q.Using (2.9) from Lemma 2.2, we can rewrite the previous inequality as   ps   q ´x∇    ps   q, s   y ě   ps   q ´x∇    ps   q, s   y. (3.17) On the other hand, taking the dot product of the equality of gradients (3.13c) with any  P Ω, we have x∇    ps   q, y " x∇    ps   q, y.
Adding together (3.17) and (3.18), we end up with   ps   q `x∇    ps   q,  ´s   y ě   ps   q `x∇    ps   q,  ´s   y which is the desired result (3.16).
This result, notoriously known as the tangent plane criterion, is usually derived by physicists from a local analysis of phase stability [22] (see also Sect.From this common tangent plane property, a purely geometric procedure can be devised in order to build a solution of the phase equilibrium problem (2.18).The construction involves the lower convex envelope of the function  Þ Ñ min PP   pq.More details will be given in Section 3.4.

Existence and uniqueness for a phasewise subproblem
The key step towards ensuring the existence of a solution to (2.18) is to study a phasewise subproblem that arises in two different forms.

Extended fugacities local inversion problem
To solve (2.18) in practice, Lauser et al. [17] advocated using the common values t  u P of extended fugacity across phases as main unknowns.This gives rise to a two-level algorithm.In the inner level, we solve  nonlinear systems of size  ˆ    Φ   p  q "   , @ P , one for each  P P.These local fugacity inversion problems express the extended fractions as implicit functions    pq of the extended fugacity vector  " `I , . . .  ˘P R  `.In the outer level, we solve one nonlinear system consisting of the  ` remaining equations min ´ , 1 ´řP    pq ¯" 0, @ P P, in the  ` unknowns `t  u PP ,   ( P ˘.This approach, the interest of which is to involve only "small" systems, was adopted by many authors [6,8,20,21]. Taking the logarithm of both sides of (3.21), writing    "      and proceeding as in the proof of Theorem 3.4, we can transform the inner system (3.21)into ln   `  p  q " ln   .
Thus, our ability to solve (3.23) for all reasonable inputs  P R  `relies on the existence of an unambiguous reciprocal function r∇  s ´1.

Extended fractions for a single-phase solution
The second situation occurs when the solution of (2.18) is single-phase, say, in phase .Hence, phase  would be entirely known.The ability to assign well-determined values to the extended fractions in an absent phase is an important feature of the unified formulation.System (3.24) has the same structure as (3.23).

Fundamental assumptions
The superiority of the unified formulation over the variable-switching formulation hinges upon the invertibility of (3.23) and (3.24), which cannot be taken for granted.To this end, additional assumptions need to be made.Below is the most natural set of assumptions.
Hypotheses 3.7.The gradient map ∇    : Ω Ñ R ´1 is surjective.Moreover, the Gibbs energy   : Ω Ñ R is strictly convex, that is, it satisfies one of the two conditions below, which are equivalent for a twice differentiable function: (a) For all p, q P Ω ˆΩ with  ‰ ,
We refer the reader to [9] for the notion of strict convexity and for the equivalence between the two conditions (a) and (b) for twice differentiable functions.Proof.Surjectivity provides existence of a solution  P Ω to ∇  pq "  for all  P R ´1 .Strict convexity enforces uniqueness of such a solution.
Hypotheses 3.7 is neither unrealistic nor unreachable, as shown by the following example.
Proof.The gradient ∇  : Ω Ñ R ´1 is given by This map is continuous over Ω.For any given  " `I , . . .,  ´1 ˘P R ´1 , the nonlinear system ∇  pq "  can be easily inverted and the only solution is ,  P tI, II, . . .,  ´1u.
This defines a unique continuous inverse map r∇  s ´1 : R ´1 Ñ Ω.From the expression (3.27) of the gradient, the Hessian matrix can be found to be where E is the matrix whose all entries are equal to 1.It follows that, for a generic  P R ´1 ,

Existence for the phase equilibrium problem in the unified formulation
Thanks to Hypotheses 3.7, a solution of (2.18) can also be worked out explicitly.Its construction is inspired by Gibbs' geometric one [12] for the two-phase binary p " 2q case.This settles the issue of existence under some minor technicalities.
Hypotheses 3.7 are taken for granted throughout this section.Additionally, we recall that the functions t  u PP are smooth (say, twice differentiable), take finite values on the boundary BΩ but their gradients blow up there, i.e., lim ÑBΩ ‖∇  pq‖ " `8.The latter is due to the presence of logarithms in the ideal parts of the Gibbs functions.The function  " min PP   (3.28) is continuous on s Ω but may not be differentiable.Let q  be the lower convex envelope of  on s Ω.By design, q  is a convex function.Like , q  is continuous.Here, we have a stronger property.
Lemma 3.10.The lower convex envelope q  is differentiable at all interior points  P Ω.
Proof.The lower convex envelope at  can be characterized as q pq " sup ℎď, ℎ affine ℎpq. (3.29) It is a convex function, which allows us to consider its subdifferential Bq pq at  P Ω.It is known that Bq pq is a nonempty and convex set [9].Let us distinguish two cases.
Case 1: q pq " pq.Let  P Bq pq.By definition of a subgradient, q pq ě q pq `x,  ´y for all  P s Ω.Let  P P such that pq "   pq.Combining the previous inequality with   pq ě pq ě q pq and q pq "   pq, we obtain   pq ě   pq `x,  ´y for all  P s Ω.This means that  P B  pq " t∇  pqu, which results in  " ∇  pq.Since the subdifferential Bq pq is reduced to a singleton, q is differentiable at . Case 2: q pq ă pq.Let  P Bq pq.Since pq ě q pq ě q pq `x,  ´y for all  P s Ω, the affine map ℎpq " q pq `x,  ´y is a legitimate "competitor" in the supremum (3.29).If the graph of ℎ does not intersect that of , namely, if ℎpq ă pq for all  P s Ω, we can find  ą 0 such that ℎpq ` ă pq for all  P s Ω, thanks to continuity of the functions and compactness of the domain.But then ℎ ` would be a better "candidate" in (3.29), as it would raise by  the value of q pq.To avoid this contradiction, there exists s   P s Ω such that ℎps   q "   ps   q " ps   q.
Let us investigate q ps   q.On the one hand, ℎps   q ď q ps   q ď ps   q.On the other hand, ℎps   q " ps   q as said above.Therefore, q ps   q " ps   q.By the same argument as in Case 1, we conclude that q  is differentiable at s   and ∇q ps   q " ∇  ps   q.From the inequality   pq ě pq ě q pq `x,  ´y and the equality   ps   q " q pq`x, s   ´y, we infer that s   achieves the minimum of the function  Þ Ñ   pq´q pq´x, ´y over s Ω.Since the latter function is strictly convex with unbounded derivatives on the boundary, the minimum cannot take place on BΩ.Thus, s   P Ω and minimality then entails  " ∇  ps   q.Hence, s   is a tangent point.At this stage, we have proved that to each  P Bq pq there corresponds a phase  P P and a point s   P Ω such that  " ∇  ps   q and q pq "   ps   q `x,  ´s   y (the last condition simply expresses that  belongs to the tangent hyperplane  s    .Assume that Bq pq contains two distinct elements  ‰ .By convexity, p1 ´q ` P Bq pq for all  P r0, 1s.To each  P r0, 1s there correspond a phase pq P P and point s  pq pq P Ω such that p1 ´q ` " ∇ pq ps  pq pqq.If necessary and up to a reparametrization, we can always take another  in this segment that is sufficiently close to  so that pq "  for all .Let   pq "  s pq   pq "   ps   pqq `xp1 ´q `,  ´s   pqy (3.30) be the value at  of the tangent map at s   pq.Since   pq " q pq for all  P r0, 1s, the derivative of   with respect to  must identically vanish.The calculation of this derivative leads to x ´,  ´s   pqy " 0. (3.31) Taking the difference of (3.31) between  " 0 and  " 1 leads to x∇  ps   p0qq ´∇  ps   p1qq, s   p0q ´s   p1qy " 0, which violates the strict convexity condition (3.25).Therefore, Bq pq is a singleton.
Thus, for  P Ω, it makes sense to speak about the gradient ∇q pq and the tangent hyperplane, defined as the graph of the linearized expansion   q pq " q pq `x∇q pq,  ´y.We introduce s Γpq " t P P | D s   P Ω,   ps   q "   q ps   q, ∇  ps   q " ∇q pqu (3.32) as the set of thoses phases whose Gibbs function   is tangent to the hyperplane   q .Lemma 3.11.For  P Ω, the following properties hold true: (1) s Γpq ‰ H. (2) For each phase  P s Γpq, the contact point s   is unique.(3) If  ď , then  P convts   u P s Γpq .
Proof.Existence of at least a contact point.The argument is similar to the proof of Lemma 3.10, with  " ∇q pq.
Uniqueness of the contact point in each phase.If the hyperplane   q  were tangent to   at two distinct points s   ‰ r   , then ∇q pq " ∇  ps   q " ∇  pr   q, and we would have x∇  ps   q ´∇  pr   q, s   ´r   y " 0, which violates the strict convexity condition (3.25).
Convex hull of the contact points.In the characterization (3.29) of q , the supremum is also a maximum reached at ℎ "   q  due to differentiability.The idea is to express this optimality by rotating the common tangent hyperplane of the contact points ts   u P s Γpq using the gradient vector  as parameter.Γpq.The equality of extended fugacities (2.18b) across phases is equivalent to that of ∇  ps   q and of ln s   μ ps   q, as was pointed out in Sections 3.2 and 3.3.On the one hand, it follows from (3.32) and (3.37b) that ∇  ps   q " ∇q pq for all  P P. On the other hand, if  P s Γpq, the common tangency q pq`x∇q pq, s   ´y "   ps   q implies that   ps   q "   ps   q ´x∇  ps   q, s   y " q pq ´x∇q pq, y.
Since s   " ř P s    " 1 after (3.37a), we have ln s   ` ps   q " q pq ´x∇q pq, y.If  P Pz s Γpq, by virtue of (3.37b), " expr  q ps   q ´ ps   qs.
(3.38) Therefore, ln s   ` ps   q "   q ps   q ´ ps   q " q pq ´x∇q pq, y.As a result, ln s   ` ps   q takes the same value for all  P P.
To prove the complementarity conditions (2.18c), we first notice from various definitions that s   ě 0 and s   `1 ´řP s    ˘" 0 for all  P P.Moreover, 1 ´řP s    " 0 for  P s Γpq.Hence, it just remains to prove that 1 ´řP s    ě 0 for  P Pz s Γpq.Starting from (3.38) and invoking the convexity of q , we have s   ď exprq ps   q ´ ps   qs ď exprps   q ´ ps   qs ď expp0q " 1, which is the desired result.
The assumption  ď  turns out to be true in practice: there are about two or three phases at most for tens to hundreds of components.For the two-phase binary case, namely, when  "  " 2, the previous solution can be proved to be unique ( [34], Thm.2.5).

Convexity analysis of simple laws
The goal of this section is to review some commonly used laws that satisfy Hypotheses 3.7 unconditionally or conditionally.Each law will be given by the excess Gibbs function Ψ  , which is connnected to the Gibbs function   by The subscript  stands for the phase to which the physical law under consideration applies.

Henry's law
In Section 3.3 (Prop.3.9), we saw that an ideal gas Ψ  " 0 fulfills Hypotheses 3.7.Next in the level of complexity is Henry's law [14] Ψ  pq " where   ( P are positive constants, each of them embodying a property of the corresponding species.The fugacity coefficients calculated by (2.14a)  Proof.Since Ψ  is affine with respect to  " `I , . . .,  ´1 ˘, its second derivatives all vanish.Therefore, the Hessian matrix ∇ 2   coincides with that of the Gibbs function of the ideal gas.But this matrix was shown to be positive definite in Proposition 3.9.We still have to check that the range of the gradient map ∇  pq " `ln `I  I ˘´ln `   ˘, . . ., ln `´1  ´1 ˘´ln `   ˘˘.
is equal to R ´1 .For a given  " `I , . . .,  ´1 ˘P R ´1 , the nonlinear system ∇  pq "  can be easily inverted and the only solution is "I expp  q{  ,  P tI, II, . . .,  ´1u.
This defines a unique continuous inverse map r∇  s ´1 : R ´1 Ñ Ω.

Margules' law
We now consider two laws dedicated to liquid binary mixtures p " 2q, namely, Margules' and Van Laar's [29].For liquids, physicists rather talk about activity coefficients instead of fugacity coefficients.This distinction is however anecdotical, since the mathematical structure of thermodynamic equilibria remains the same [27].
Since  ´1 "  The region indicated by (4.6) is colored in striped green in Figure 1.Its right-most point is located at p, q " p4, 0q, where it has a vertical tangent.
Proof.The first derivative of   is This is a continuous function over p0, 1q, with lim Ó0  1  pq " ´8 and lim Ò1  1  pq " `8.Thus,  1  has range in R.
Let us change the variable to  "  ´1 2 P `´1 2 , 1 2 ˘to work with the more symmetric function It remains to study  in order to determine the region of strict positivity pq ą 0. The technical details can be found in [16] or Proposition 3.2 of [34].

Van Laar's law
Van Laar's law is also a model for activity coefficients of a liquid [29].The excess Gibbs function associated with it is where p 12 ,  21 q P pR ˚q2 are two nonzero parameters.By (2.14a), the fugacity coefficients are ln Φ  To make sure that formulae (4.7) and (4.8) are well-defined over  P p0, 1q, the denominator  12  `21 p1 ´q must keep the same sign.This amounts to requiring that  12  21 ą 0. (4.9) Besides, the pair p 12 ,  21 q must be further restricted in order to comply with Hypotheses 3.7. where ℛ `" " 0 ă  ă 4 and || ă min ˆ; The region indicated by (4.10) is colored in yellow in Figure 2. It lies inside the cone  2 ă  2 that corresponds to condition (4.9).The origin p0, 0q must be excluded.
Proof.The first derivative of   is  ´2 " 1 2 p 12 `21 q `p 12 ´21 q remains to study  in order to determine the region of strict positivity pq ą 0. The technical details can be found in [16] or Proposition 3.3 in [34].

Cubic equation of states for two-phase mixtures
The fugacity laws investigated in Section 4 are simple and apply to a selected phase , regardless of the remaining ones.We are now going to examine a prominent category of laws for a two-phase (gas and liquid) mixture, in which the fugacity coefficients for both phases are computed in a "simultaneous" way.Throughout the rest of this paper, it is therefore assumed that The new labels  (gas) and  (liquid) are aimed at being more meaningful.To fix ideas, the presentation is done for Peng-Robinson's law [28].The philosophy is the same for other cubic laws.

Mixing rules and cubic equation
Each component  P  in a pure state is characterized by a pair of positive parameters   (cohesion term) and   (covolume).These are highly sophisticated functions of the pressure and the temperature, but at fixed pP, Tq can be considered as constants.This gives rise at fixed pP, Tq to a pair of positive dimensionless parameters As before, we shall never write down explicitly the dependency of ` ,   ˘on pP, Tq.
A multicomponent mixture is supposed to behave approximately as a fictitious pure component endowed with an averaged value for the pair p, q.The latter is computed from the ` ,   ˘'s and the current partial fractions by means of a mixing rule.More specifically, let  P s Ω be the partial fractions of some phase.There can be found [27] a wide variety of mixing rules  Þ Ñ ppq, pqq.We require mixing rules to be smooth and to satisfy the compatibility relation ``  ˘,  ` ˘˘" ` ,   ˘for all  P .We recall that   " p ,1 ,  ,2 , . . .,  ,´1 q was introduced in Section 2.1.2for  P  and consists of elementary Kronecker symbols  , .
The next step is to consider the cubic equation in the variable pq.This accounts for the name "cubic EOS."The dimensionless quantity pq, called compressibility factor, can be intuitively understood by noting that for a pure component, the cubic equation simply results from an algebraic transformation of the equation of state using Thus, for an ideal gas p "  " 0q, we have  " 1.
In the most favorable situation, there are three real roots, all greater than pq.These are then named In other words, the smallest root is associated with the liquid phase , while the largest one is associated with the gas phase .This is grounded on the physical fact that, at the same pressure and temperature, the gas phase occupies a larger volume than the liquid phase, which by (5.6) implies that   pq ą   pq.As for the intermediate root   pq, it does not have any physical meaning.

Gibbs functions and fugacity coefficients
Let  P t, u and assume that the real root   pq ą pq is well-defined.Then, the Peng-Robinson excess Gibbs energy is defined as From this, the fugacity coefficients can be deduced with the help of (2.14a).
Theorem Applying (2.14a) and using (5.8), we arrive at the desired result.

Two crucial questions
Formulae (5.8) and (5.9) are well-known to thermodynamicists.A delicate but less often clarified issue is to know which phase  P t, u they can be applied to, especially in the unfavorable situation when equation (5.3) has only one real root greater than pq.The simple-minded idea of taking   "   equal to this real root is of common practice in industrial codes, but is ill-advised since it gives rise to discontinuities in the Gibbs functions, as will be explained in Remark 6.1.
In fact, in the one-root scenario, two subcases have to be envisaged.If we manage to assign a "natural" phase label  "  or  to the real root, then the corresponding excess Gibbs energy Ψ  is defined by (5.8), leaving its counterpart in the other phase undefined.If we do not succeed in attributing a "logical" phase label to the real root, then Ψ  is undefined in both phases.This process is intuitive enough to describe with words, but raises two serious questions: (1) When does the cubic equation have three real roots greater than pq and when does it have only one real root greater than pq?(2) When and how can a "natural" phase label be assigned to the unique real root greater than pq?
The upcoming subsections answer to these questions by working on the generic form (5.4).Part of these issues was already addressed in [18] for Van der Waals' law.We are not aware of any similar work for Peng-Robinson's law.This is why we are taking this opportunity to undertake a rigorous study.

Assignment of phase labels to roots
Instead of the polynomial Υ , pq "  3 `p ´1q 2 `` ´2 ´3 Π , pq " ´1, (5.12) there is at least one root larger than .
( In physics textbooks [27,33], only decimal approximations (5.14d) of the critical point can be found, without any proof.The interest of Lemma 5.2 is to derive the exact values (5.14a)-(5.14c) of the critical point.
Proof.System (5.13) can be turned into a set of 3 polynomial equations in p  ,   ,   q.By eliminating   , we obtain the cubic equation  3   ´3 2  ´3  ´3 " 0 in   "   {  , whose only real root is From this   {  can be deduced exactly.Once this is done, we can get back to p  ,   ,   q.See Lemma 3.3 in [34] for more details.
(1) If { ą   {  « 0.170144420, the function Π , is decreasing over p, `8q and has only one zero greater than .(2) If { ă   {  « 0.170144420, the function Π , has two disctinct local extrema.In other words, there exist two distinct values   ă   in p, `8q such that Π 1 , p  q " Π 1 , p  q " 0.Then, Π , is decreasing on p,   q, increasing on p  ,   q and decreasing on p  , `8q.It may have one or three distinct zeros over p, `8q.Theorem 5.3 paves the way to a natural association of a root with a phase in the subcritical regime.Note that { plays the role of a temperature (up to a multiplicative constant).Definition 5.4 (Phase label assignment).The region 0 ă  ă p  {  q is said to be subcritical.In the subcritical region, a root  ą  of the cubic equation (5.4) is said to be associated with the liquid phase  if  ă   ; a root  ą  of the cubic equation (5.4) is said to be associated with the gas phase  if  ą   .
Let us elaborate on this definition before proving Theorem 5.3.If there is only one root  ą , this root cannot belong to p  ,   q.Therefore, either  P p,   q or  P p  , `8q.This way of assigning a phase label to  is most natural, since it extends by continuity the "topological" pattern observed in the case of three real roots.
The region  ą p  {  q is said to be supercritical.The graph of Π , no longer has two discernable branches.In this configuration, there is no natural way to associate  with a phase.Physically speaking, the critical threshold   {  corresponds to a critical temperature T  .Above the critical temperature, the distinction between gas and liquid phases no longer holds [12].Proof.The discriminant3 of the cubic equation (5.4) is ∆p, q " ´4 3 ´`8 It can be shown that ∆  pq ą 0 for  P p0,   q, ∆  p  q " 0 and ∆  pq ă 0 for  ą   .Therefore, if  ą   , only  0 pq exists.If  P p0,   q, there exist  0 pq ă   pq ă   pq.
The rest of the proof goes as follows.We first show that  0 pq ą 0.Then, we rule out the region t0 ă  ă   , 0 ă  ă  0 pqu which in fact belongs to the supercritical region.Next, in the region tp  {  q ă  ă  0 pq,  ˚ă u, where  ˚« 2.435425 is the ordinate at which the graphe of  0 p¨q enters the subcritical region, we show that the three real real roots cannot be all larger than .In conclusion, the only way for (5.4) to have three real roots, all greater than , is that  P p0,   q and  P p  pq,   pqq.This region is shown to be contained in the subcritical domain.The comprehensive discussion can be found in Theorem 3.6 from [34].

Domain extension for cubic EOS-based Gibbs functions
From Section 5.3, it appears that the cubic equation (5.3) may not always have three real roots greater than pq for all  P s Ω.As a consequence, the domain of definition for the functions Ψ  , Φ   for a given phase  may not always cover the whole simplex s Ω.This turns out to be detrimental to the unified formulation (2.18).

Difficulty of the unified formulation with cubic EOS
In a nutshell, the Gibbs energy functions   may grossly violate Hypotheses 3.7.To give a visual picture of the nature of the obstruction, let us consider the simplistic case of a two-phase binary p " 2q mixture, governed by Peng-Robinson's law combined with the mixing rule  Assume that `I ,  I ˘belongs to the one-root region labelled , while `II ,  II ˘belongs to the one-root region labelled , as illustrated in Figure 4.At  " 0, the curve  starts from `II ,  II ˘in the -root region.At some parameter value  "  5 P p0, 1q, it enters the three-root region.At some further value  "  7 P p 5 , 1q, it exits the three-root region.At  " 1, it finally meets `I ,  I ˘in the -root region.It is not difficult to realize that: -the quantities   pq, Ψ  pq,   pq are well-defined only for  P r0,  7 s;   ´7 ¯and   ps   q "  1  pq.Figure 5 depicts this situation.It could be argued that the same flaw of cubic EOS laws should cause the same prejudice to the variableswitching formulation of Section 2.2.1.But this is not true.In the variable-switching formulation, if the context is correctly guessed, we do not need to compute anything from the absent phase and the above problem is irrelevant.If the context is incorrectly alleged, the flash does not converge or may even crash, but there is an opportunity for us to make up for it by changing the context.The natural variable formulation does not seek to explore the regions where information is missing.The unified formulation has to do so, by its very vocation to treat all phases on an equal footing.Remark 6.1.From Figure 5, it can be seen that if we abruptly take   "   in the one-root regions  P r0,  5 q and  P p 7 , 1s, as often done by practitioners, then we will have   "   over these two intervals.As a consequence,   will exhibit a discontinuity at  5 and   a discontinuity at  7 .These unphysical discontinuities are not a favorable feature for numerical robustness.

A natural domain extension method
To enhance the performance of the unified formulation, it is essential that the domains of definition for the excess functions Ψ  can be extended to s Ω.Note that here we just want to extend the domains of definition of various functions.We do not strive to fulfill Hypotheses 3.7, since these assumptions may already be violated for the original unextended Gibbs functions.Even without strict convexity, a smooth extension of the Gibbs functions helps iterative methods [35] to remain well-defined everywhere.

Construction in the one-root region
When the cubic equation does not have three real roots, our idea is to use the common real part of the two complex conjugate roots, as a "surrogate" of the missing real root.Assume that   is the only real root greater than  of Peng-Robinson's cubic equation  3 `p ´1q 2 `` ´2 ´3 2 ˘ `` 2 `3 ´ ˘" 0, where the label  P t, u has been assigned in accordance with Definition 5.4.To alleviate notations, we do not explicitly indicate the dependency of ,  and  on .
Let  be the other phase, that is,  "  if  "  and  "  if  " .Since the sum of the three (complex) roots is equal to 1 ´, the two remaining conjugate roots share the common real part The following properties of   speak in favor of its enrollment as a substitute for   , which would have been subject the same constraints, had it existed.Lemma 6.2.Let p, q be a pair in the subcritical region 0 ă  ă p  {  q and assume that Peng-Robinson's cubic equation has only one real root   ą  that corresponds to phase .
Proof.The proof is similar to that of Theorem 5.1.
The gradient of   required by (6.5) can be computed by ∇  " ´1 2 p∇ `∇  q, in which ∇  comes from differentiating Peng-Robinson's cubic equation with respect to , i.e.,

Alteration in the three-root region
In the one-root region, the gradient ∇  " ´1 2 p∇ `∇  q remains bounded.If we start from the threeroot region and move towards the transition boundary where   disappears, the gradient ∇  does not remain bounded.Indeed,   also obeys (6.6) (just replace   by   ), and as   gets closer to being a double root, ∇  blows up.However, we need a finite gradient ∇  for the numerical resolution of system (2.18) by, say, the Newton method.
To achieve a smooth junction, we introduce a further approximation on a tiny portion of the three-root region.Assuming that we are in the three-root region, with  ă   ă   ă   , we introduce

𝜗 "
´   ´ P r0, 1s (6.7) as an indicator of the closeness to the transition boundary.Indeed, the cubic equation has double roots when  " 0 or  " 1.Let  P p0, 1{4q be a small threshold.
-If  P r2, 1 ´2s, we apply the usual formulas for the case of three real-roots.
-If  P p1 ´2, 1s, the two roots   and   are close to each other.We keep   but progressively replace   by   " 1 2 p1 ´ ´ q " 1 2 p  ` q whose gradient is bounded.Instead of   , we plug r   " r1 ´ pqs  ` pq  into (5.8)

Numerical validation of the extension method
Extensive numerical simulations are provided in [35], a companion paper to the present one, to demonstrate the relevance of the above extension method.In [35], we considered various systems of equations in the unified form (2.18) with a wide range of physical parameters and initial points.A careful comparison is carried out between two numerical methods used to solve these systems: the Newton-min method and the Non-Parametic Interior Point Method (NPIPM) that we designed on purpose to deal with such systems.
It turned out that very good results (nearly 100% convergence) can be achieved thanks to the combination of both ingredients, i.e., the extension of Gibbs functions and the NPIPM algorithm.A single ingredient alone is not enough in the following sense: unextended Gibbs functions always cause divergence of both numerical methods (Newton-min and NPIPM), but extended Gibbs functions combined with Newton-min does not bring significant improvement.

Conclusion
Beyond implementational advantages, the unified formulation has been shown to be able to recover all the properties known to physicists on phase equilibrium.Indeed, the complementary equations do encapsulate the tangent plane criterion (Thm.3.4), as the sign information is related to some stability condition.The unified Figure 9. Close-up comparison of the derivative of the extended Gibbs functions between  " 0.001 and  " 0.2 for Peng-Robinson' law with the indirect method.p I ,  I q " p0.275, 0.045q and p II ,  II q " p0.35, 0.04q.(A)  1   .(B)  1  .
formulation can also be regarded as a solution of some constrained minimization problem (Thms.3.2 and 3.3) in which the objective function is a Gibbs energy of the mixture.This solution is characterized by a set of equations that is slightly stronger than the KKT optimality conditions when a phase vanishes.The possibility of assigning well-defined values to the extended fractions of an absent phase can only be achieved if the Gibbs functions meet some restrictive requirements.Strictly convexity and surjectivity of the gradient over the whole domain of fractions (Hypotheses 3.7) are sufficient for this purpose.Remarkably, these assumptions also guarantee the existence of a solution to the phase equilibrium problem in the unified formulation (Thm.3.12).
Unfortunately, Hypotheses 3.7 are not satisfied by realistic Gibbs functions.The obligation of assigning well-defined fraction values to an absent phase then becomes a weakness that jeopardizes the whole unified approach.This is especially true for the Gibbs functions derived from cubic EOS, for which they are not even defined on the whole domain of fractions.The extension procedure of Section 6 is aimed at improving the robustness of the unified formulation.The corresponding numerical results will be the subject of another simulation-oriented article, where we also design an interior-point algorithm ( [34], Sect.5) in order to efficiently cope with complementarity conditions.
Despite its solid theoretical foundation, the current unified formulation is not able to support the phenomenon of phase separation, where the same phase is split into two or several distinct subphases due to the non-convexity of its Gibbs function.Note that the variable-switching formulation cannot do it either.Future works should tackle this question perhaps by combining the unified formulation with some judicious approaches advocated by [11,25].

Theorem 3 . 3 .
Let!´r  , r   ¯)PP be a critical point of the minimization problem (3.3).

Figure 1 .
Figure 1.Region of strict convexity for the parameters of Margules' law in the p, q-plane.

Figure 2 .
Figure 2. Region of strict convexity for the parameters of Van Laar's law.

Figure 3 .
Figure 3. Number of roots for Peng-Robinson's law in the p, q-quarter plane.
I `p1 ´q II .(6.1)For an arbitrary choice of the two pairs `I ,  I ˘and `II ,  II ˘in the subcritical region 0 ă  ă p  {  q, the parametrized curve  : r0, 1s Q  Þ Ñ ppq, pqq P `R˚˘2 is an arc of parabola.

Figure 4 .
Figure 4. Curve  defined by the mixing rule in the p, q-plane.

Figure 5 .
Figure 5.Typical situation where the fraction in the absent phase cannot be computed.

Figures 6
Figures 6 and 7 display a few examples of the extension method for Peng-Robinson's law in the binary case.Figures 8 and 9 provide a close-up comparison between two choices of .
labeled by Arabic numerals.The choice of P within a model is the (difficult) task of physicists:  should be large enough to take into account the appearance of new phases in models with time evolution, but not too large for computations to remain feasible.In reservoir simulations, the maximum number of possible phases  " |P| is commonly about 3.The relative importance of each phase  P P within the mixture is measured by the phasic fraction   P r0, 1s, such that ÿ   " 0 is said to be absent.Otherwise, it is present.The subset of present phases, namely,Γ " t P P |   ą 0u Ă P (2.4)is referred to as the context.For a present phase  P Γ, it is possible to compare the relative contribution of each component  P  within it by defining the partial fractions    P r0, 1s, such that P " t1, 2, . . .,  u,  ě 2, (2.2) whose elements are  " `I  , . . .,  ´1  ˘P Ω Ă R ´1 is called partial composition of phase .It consists of only the first  ´1 partial fractions, since the quantities  I  ,  II  , . . .,    are not independent, in view of (2.5).Whenever a    turns up in any formula, it should be interpreted as    " 1 ´I  ´. . .´´1  .The domain of   is the closure of ).For all  P Ω: Proof.Multiplying (2.9) by   , summing over  P  and noticing that ř P     " , we end up with (2.10a).To prove (2.10b), we subtract the last potential    pq "   pq `A∇  pq,   ´E from each    ,  P ztu, given by are given empirically or inferred from an equation of state.
ą 0 and r   " 0. It is no longer possible to divide (3.4b) by r   to retrieve information on the extended fugacities.Likewise, we now simply have   ě 0 from (3.4e).Equation (3.4a) for phase  leads to g  ´r   ¯`r .7)This implies that, in general, the complementarity condition (2.18c) does not hold for phase  and the extended fugacity equalities (2.18b) do not hold between  and .But (2.18c) is automatically met for phase  as soon as(2.18b)holds between  and .P , t r   u PP ˘be a solution of the KKT system (3.4).First, assume that r   ą 0 and r   ą 0. Dividing (3.4b) by r  , we obtain B   g  ´r   ¯`r   " 0 and B   g  ´r   ¯`r   " 0. From this, we deduce that B   g  ´r   ¯" B   g  ´r   ¯" ´r   .According to (3.2a) (Lem.3.1), this is equivalent to the equality of extended fugacities (2.18b), rewritten in the second part of (3.6).On the other hand, r   ą 0 implies r  ´r   ¯" 0. Combining this with (3.2b) (Lem.3.1), we infer that r   " r .Repeating the same reasoning for , we also get r   " r .Hence, r   " r   .This means that r  takes on the same value r  in all present phases.Let r Γ be set of  P P such that r   ą 0. Note that r Γ ‰ H because of (3.4c).Summing (3.4d) over  P  and  ´r   ¯´Bg  B  ´r   ¯ ´r   `1 ě 0.
p  qu PPztu , t r   p  qu PP .From the KKT conditions for (3.10) subject to (3.9b)-(3.9d), it follows that (see [34], Sect.2.3.2.4 for details) r    p  qΦpr   p  qq " r    p  qΦpr   p  qq for all  P .Now, we let   Ó 0. If all of the quantities involved in the above equality have finite limits, we clearly end up with (3.8).
PPztu  g  p  q ` g  p  q (3.10) for a fixed   ě 0. Assume that for each small enough   ą 0 there is a unique critical point, denoted by t r For an interior point   P Ω, we designate by     the tangent hyperplane to   at   .This tangent hyperplane, which exists thanks to the regularity assumptions on   , is the graph of the affine function     : R ´1 Ñ R defined as     pq "   p  q `x∇  p  q,  ´ y.   ˘PP ˘exists to the unified formulation(2.18), in which the renormalized fractions s   P Ω are computed from s   by (2.18d).We are going to learn as much as we can about it.Theorem 3.4.For any pair p, q P P ˆP of phases, present or absent:(1) The potentials in phase  are shifted from their counterparts in phase  by a same constant, i.e., ).The generic element of Ω ˆR is denoted by p, q.Let   " p, q P Ω ˆR |  "   pq ( (3.11) be the graph of the Gibbs energy function   : Ω Ñ R. s    and  s     are parallel.Put another way, ∇  ps   q " ∇  ps   q. (3.13c) Proof.For each phase  P P, let us define s   as in (3.13b), so that for  P  we have s  ps   q " ln s   `  ps   q. (3.15)From this, we deduce the translation property (3.13).Subtracting the last equality ln s   `  ps   q " ln s   μ  ps   q from (3.15) and recalling (2.10b) (Lem.2.2), we have B  B  ps   q " Theorem 3.5 (Tangent plane criterion).Assume that a phase  P P is present, i.e., s   ą 0.Then, for any other phase  P P, absent or present, From equality (3.13a), we have    ps   q "    ps   q ` ,   " ln s   ´ln s   .
s     pq ě  s    pq, for all  P R ´1 .(3.16) Thus, the tangent hyperplane  s     lies above or coincides with the tangent hyperplane  s    .Proof.
3.1).Theorem 3.5 testifies to the fact that this stability property is already encoded in the unified formulation via the sign of 1 ´řP s     " 1, the point  belongs to the interior of convpts   u P s Γ q.
.If phase  is "strictly" absent, namely, if 1 ´řP s    ą 0 and s   " 0, then the tangent hyperplane  s     will lie strictly above  s    .We now push one step further by looking at the case of several present phases.Let s Γ be the set of all  P P such that s   ą 0. Corollary 3.6 (Common tangent hyperplane).At a solution of the unified formulation satisfying s   P Ω for all  P P, the tangent hyperplanes t s    u P s Γ , are all the same.Moreover,  " `I , . . .,  ´1 ˘P intpconvpts   u P s Γ qq, (3.19) i.e., the global composition point belongs to the open convex hull spanned by the points ts   u P s Γ .Proof.Let p, q P s Γ ˆs Γ. Applying Theorem 3.5 twice and switching their roles, we have  s     pq ě  s    pq and  s    pq ě  s     pq, whence  s    pq "  s     pq for all  P Ω.Thus,  s    "  PΓ s After Theorem 3.4, the extended fractions in a vanishing phase  P Pztu satisfy ∇  ps   q " ∇  pq, Put another way, s   " 1 and s   " 0 for all  P Pztu.By (2.18a), rewritten as (3.20), we have s   " .Assume  P Ω.  ps   q "    pq, (3.24b)If the function ∇  were invertible, we could write s   " r∇  s ´1p∇  pqq.Then, it could deduced from (3.24b) that s   " expr   pq ´  ps   qs and s    " s   s    .
For  P R ´1 , we define r   pq " r∇  s ´1pq for  P s Γpq.Note that r   p∇q pqq " s   .The tangent map of   at r   pq reads  , the number of constraints does not exceed the space dimension and problem (3.34) keeps a chance of being feasible.The objective function (3.34a) is the value taken by this common tangent hyperplane at .∇q pq achieves a local optimum of(3.34).The first-order optimality conditions for (3.34) imply  ´∇ αpq ` `∇ αpq ´∇ β pq ˘`...` p∇ αpq ´∇ ω pqq " 0,(3.35)where,...,   are the Lagrange multiplies associated with the constraints (3.34b).Plugging  " ∇q pq into (3.35) and using ∇ αpq " r   pq, we end up with " p1 ´ ´. ..qs   ` s   `. ..` s   ` s   .(3.36)Since the coefficients in the right-hand side sum to 1, at least one of them must be positive.Up to a permuation of s Γpq, we can assume that 1 ´ ´. ..´ ą 0. Let us prove that the other coefficients are nonnegative.Suppose that   ă 0. The idea is now to rotate the common tangent hyperplane but to leave out the tangency constraint for , so that the new affine function becomes strictly lower that   , remains tangent to the other Gibbs functions in s Γpq and achieves a higher value at , which is contradiction.Let  " ∇q pq`, where  is orthogonal to the subspace spanned by s   ´s   , ..., s   ´s   .Since s   ´s   " ´∇ α ´∇ β ¯p∇q pqq and similarly for other phases, the | s Γpq|´2 constraints  αpq "  β pq, ...,  αpq "  π pq remain satisfied at first-order expansion.Let   pq "   ps   q be the values of the new hyperplane evaluated at  and s  Taking the dot product with  yields x∇  p∇q pqq, y "   x∇ s  pq pqq, y.Hence, it is possible to choose  so as to increase   and to decrease  s  , in other words such that  r pq   pq ą q pq and  r pq   ps   q ă   ps   q at first-order expansion.The affine function ℎ "  r pq   would then be a strictly better candidate than   q  in (3.29).This is impossible.The last property means that  is a convex combination of the contact points, that is, there exist t s   u P s   " .We are now ready to describe a solution.Theorem 3.12.Assume  ď ,  P Ω.Let ts   u P s Γpq , t s   u P s Γpq be defined as above and setThis procedure supplies us with a solution of (2.18).
r pq   pq "   pr   pqq `x,  ´r   pqy " x, y ´α pq, (3.33)where  α stands for the Legendre conjugate of   .Let s Γpq " t, , . . ., , u and consider the maximization problem max   r pq   pq (3.34a) subject to the | s Γpq| ´1 equality constraints  αpq "  β pq,  αpq "  γ pq, . . .,  αpq "  π pq,  αpq "  ω pq.(3.34b)The constraints (3.34b) are aimed at making the | s Γpq| functions (3.33) coincide with each other, so as to preserve common tangency.Since | s Γpq| ď  ď If  stays in a small enough neighborhood of ∇q pq, then  r pq remains below  (thanks to the compactness of s Ω) and ℎ "  r pq can be considered as a valid "candidate" in (3.29).Therefore, it is expected that  " r pq   .It is straigthforward to show that ∇  pq "  ´r   pq and ∇ s  pq " s   ´r   pq, so that ∇  p∇q pqq "  ´s   and ∇ s  pq pqq " s   ´s   .In view of (3.36),  ´s   "   ps   ´s   q `. . .` ps   ´s   q ` ps   ´s   q.
This is why this law is also called the constant coefficients law.Proposition 4.1.For all `I , . . .,   ˘P `R˚˘ , the Gibbs energy function   associated with Henry's law fulfills Hypotheses 3.7.
1, we simply write  instead of  I and .The excess function associated with Margules' law is Ψ  pq " p1 ´qr 12 p1 ´q `21 s, (4.4) where p 12 ,  21 q P pR ˚q2 are two nonzero parameters.By (2.14a), the fugacity coefficients are To meet Hypotheses 3.7, the pair p 12 ,  21 q must be restricted to some region of R 2 .Proposition 4.2.Let  "  12 `21 and  "  12 ´21 .Then, the Gibbs energy function   associated with Margules' law fulfills Hypotheses 3.7 if and only if The second derivative of   , multiplied by p1 ´q to get rid of singularities, is equal to 2Ünder assumption (4.9), this is a continuous function over p0, 1q, with lim Ó0  1  pq " ´8 and lim Ò1  1  pq " `8.Thus, 1 has range in R.
5.1.The Peng-Robinson fugacity coefficients are given by  and for any phase  P t, u in which   pq ą pq is well-defined.
´1˘, are both lesser than , Π , and Υ , have the same roots over p, `8q.
5.13)Critical points, also called triple points, are physically important.Here, this notion will help us divide the space of parameters into various subregions with physically distinct behaviors.