Next Article in Journal
Network Analysis of Depression Using Magnetoencephalogram Based on Polynomial Kernel Granger Causality
Previous Article in Journal
What Is It like to Be a Brain Organoid? Phenomenal Consciousness in a Biological Neural Network
Previous Article in Special Issue
Thermodynamic Entropy as a Noether Invariant from Contact Geometry
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Using the Intrinsic Geometry of Binodal Curves to Simplify the Computation of Ternary Liquid–Liquid Phase Diagrams

by
Nataliya Shcherbakova
*,†,
Vincent Gerbaud
and
Kevin Roger
Laboratoire de Génie Chimique, Université de Toulouse, CNRS, INP, UPS, 31432 Toulouse, France
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Entropy 2023, 25(9), 1329; https://doi.org/10.3390/e25091329
Submission received: 29 June 2023 / Revised: 30 August 2023 / Accepted: 11 September 2023 / Published: 13 September 2023
(This article belongs to the Special Issue Geometric Structure of Thermodynamics: Theory and Applications)

Abstract

:
Phase diagrams are powerful tools to understand the multi-scale behaviour of complex systems. Yet, their determination requires in practice both experiments and computations, which quickly becomes a daunting task. Here, we propose a geometrical approach to simplify the numerical computation of liquid–liquid ternary phase diagrams. We show that using the intrinsic geometry of the binodal curve, it is possible to formulate the problem as a simple set of ordinary differential equations in an extended 4D space. Consequently, if the thermodynamic potential, such as Gibbs free energy, is known from an experimental data set, the whole phase diagram, including the spinodal curve, can be easily computed. We showcase this approach on four ternary liquid–liquid diagrams, with different topological properties, using a modified Flory–Huggins model. We demonstrate that our method leads to similar or better results comparing those obtained with other methods, but with a much simpler procedure. Acknowledging and using the intrinsic geometry of phase diagrams thus appears as a promising way to further develop the computation of multiphase diagrams.

1. Introduction

Equilibrium phase diagrams testify how a complex meso/macroscopic behaviour emerges from a diversity of intermolecular interactions. Their knowledge thus guides us into unravelling the multi-scale relationship that underlies the properties of a given mixture, whether at equilibrium or during a dynamic process, for instance, liquid–liquid separation. In practice, gaining such knowledge requires either a straightforward computation from the prior knowledge of the thermodynamic potential (direct problem), or the determination of such a potential (inverse problem), usually the Gibbs free energy for liquid mixtures, from the experimental results. These direct and inverse problems are interrelated, as it is, in reality, impossible to obtain the expression of the thermodynamic potential without any prior experimental knowledge. Both experiments and computations are thus required to obtain phase diagrams. While this topic has received a lot of attention, due to its central role, it strikingly remains challenging nowadays to effectively compute phase diagrams from an experimental data set. Fast and reliable algorithms are needed to solve the computational problems arising such as the identification of model parameters from the experimental data, the localization of phase separation envelopes, and the detection of multi-phase regions. A vast literature points out the trial of treating even the common case of ternary systems that can undergo phase separation into two phases.
This observation is striking since the key features of the phase diagram of a common ternary mixture are well known and can be considered textbook material, for instance, binodal and spinodal curves, together with critical or plait points. These geometrical elements define the boundaries of stable, metastable, and unstable domains of the ternary mixture in the phase diagram. While the experimental determination of the spinodal curve is intrinsically challenging, many techniques have been devised to obtain the binodal curve, from which a reliable model can then be derived to predict the location of the spinodal curve [1]. Still, the nature of the binodal curve is often ignored by experimentalists, who usually record the location of phase separation or solubilization. While correct from a phenomenological standpoint, half the information is then lost as the binodal curve is actually the combination of two branches of conjugated compositions, or nodes.
From an adequate experimental data set, the first computational problem is an inverse one. The thermodynamic potential must be derived from the experimental data. A vast literature treats the techniques to solve the inverse problem using different optimization algorithms, including stochastic optimization [2], genetic algorithms [3], the ant colony optimization method [4], and others. The most popular optimization criteria include the least-square minimization on the distance between the experimental and model estimated points and the tangent plane distance function, although other formulations are implemented in stochastic optimization algorithms [2]. We refer the interested reader to the review paper [5] for the detailed study of existing numerical approaches for parameter estimation in the phase equilibria context. We also remark that most of the published literature deals with NRTL or UNIQUAC models.
The second computational problem is a direct one. Both the binodal and spinodal curves should be computed over the whole composition range from the expression of the previously extracted expression of the thermodynamic potential. In order to compute the binodal curve, one needs to solve the set of algebraic equations expressing the equality between the chemical potentials in each phase for every component of the multicomponent mixture. Different methods are currently used to address this problem. The simplest one consists in solving these algebraic equations by a kind of Newton–Raphson iterative procedure over the discrete mesh approximating the state space [6]. Clearly, the accuracy of the result, as well as the computational effort is directly related to the finiteness of the mesh that is used. The famous liquid–liquid multi-phase flash algorithm proposed by Michelsen [7] uses the phase stability analysis based on the minimization of the distance between the tangent plane and the Gibbs energy surface, and hence explores the geometrical properties of the potential surface. Recently, Binous and Bellagi [8] proposed to use the arc-length continuation method, reducing the computation of the binodals to the integration of differential-algebraic (DEA) systems of equations with high accuracy. It is worth also citing the homotopic methods, see [9] and references therein, which allow solving simultaneously direct and indirect problems.
Interestingly, both computation problems are challenging due to the intrinsic geometry of the equilibrium curves that we briefly recapitulate. Since the pioneering works of D.J. Korteweg [10], the mathematical description of phase diagrams can be set in terms of topological properties of the surface associated with the appropriate thermodynamic potential, for instance, the Gibbs free energy G. The problem is rather simple in two-parametric systems (one-component biphasic system with varying temperature or pressure, bi-component systems under isobaric–isothermal conditions, etc.). In this case, the analysis reduces to the detection of the bitangent lines and the inflection points of the graph of the potential. However, in the ternary case, this picture becomes much more complicated.
Indeed, the binodal curves result from the projections of the one-parametric families of bitangent lines to the surface, while the spinodal curves are the projections of the curves of zero Gauss curvature. Starting from the works of Arnold [11] and Varchenko [12], the types of possible singularities associated with different types of thermodynamic potentials have been the subject of many studies; see, for instance, refs. [13,14,15] and references therein. Most of these works treated the binary systems with varying pressure or temperature, though the developed methods can be generalized to any type of diagram.
In this work, we propose an algorithm able to treat both direct and inverse problems by the integration of ordinary differential equations (ODE). This huge simplification rests on the reformulation of the mathematical problems in a space adequate to the intrinsic geometry of ternary phase diagrams. More precisely, our key idea is to reformulate the definition of the binodal curves in an extended 4D space by associating a proper configuration space to each phase of the system. To this end, the phase coexistence conditions are rewritten in geometrical terms involving the notions of bitangent planes and conodal pairs of points. The binodal curves are then shown to be the projections of the integral curve of a certain vector field in the extended space. The same approach applied in the 2D composition plane yields the vector field generating the spinodal curves.
Using the described geometrical construction, the numerical computation of binodal and spinodal curves can be performed via the standard integration of a system of ODE equations by any conventional ODE solver. From the numerical point of view, the proposed method implements a kind of a differential path-following algorithm. This method is here applied to detect the binodal and spinodal curves of ternary liquid mixture diagrams of types 0, I, and II, defined according to Treybal’s [16] classification. The proposed approach can be implemented with a variety of models for the thermodynamic potential, in particular, excess Gibbs free energy models like NRTL or UNIQUAC. We chose here the expression based on the analytical Flory–Huggins model, modified with a ternary cross-term that accounts for the oversimplifying hypotheses underlying the standard Flory–Huggins model. Interaction parameters were found through a non-linear optimization procedure associated with a non-standard criterion, which accounts for the intrinsic geometry of the binodal curve. This choice of model is not exclusive to the computation of the binodal and spinodal curves via ODE integration, but it has several remarkable advantages. First of all, the Flory–Huggins model provides a good representation of mixtures composed of molecules of different lengths, for instance, when one or more compound is a polymer, as suggested by [1]. On the other hand, the linearity of the Flory–Huggins expression of the Gibbs free energy of mixing with respect to unknown parameters facilitates the resolution of the inverse LLE problem. Indeed, some of the unknown parameters can be computed directly in the binary case, provided the miscibility gap is known experimentally. This fact leads to an important simplification of the fitting procedure in the ternary case for the diagrams of type I and II.
The paper is organized as follows. In Section 2, we recall how the phase separation conditions of a biphasic system maintained at thermodynamic equilibrium are related to the Gibbs free energy surface topology. In Section 3, we present the main conceptual result of the paper. The problem is considered in the extended 4D space having the structure of product space of two copies of 2D configuration space associated with each phase. It is shown that the phase coexistence conditions define a smooth curve in this space, referred to as the generalized binodal curve. Each point of this curve projects on a tie-line of the phase diagram. The binodal curve computation is then reduced to the numerical integration of a system of four differential equations. Apart from the special case of zero (or “island”)-type diagrams, the starting point of such integration can be found by solving at most three binary problems on the boundaries of the composition triangle. In Section 4, the developed approach is applied to the analysis of the series of examples of ternary mixtures, modelled by using the Flory–Huggins equation with an additional triple interaction term. We show that the novel computation method performs nicely, and our fitting results are of similar to better accuracy relative to other authors who have used NRTL or UNIQUAC models with parameters carried out with much more powerful optimization algorithms.

2. Phase Separation: From Thermodynamics to Geometry

2.1. Phase Coexistence in Multicomponent Mixtures at Thermodynamic Equilibrium

Consider an N-component system of volume V characterized by temperature T, pressure P, and entropy S. Denote by n i the number of moles of the i-th component and set n = ( n 1 , , n N ) . By choosing P, T, and n as the coordinates of the thermodynamic state space, the physicochemical properties of the system as a whole can be described by the Gibbs free energy G ( P , T , n ) . Being a homogeneous function of the first order with respect to n i , G can be expressed as
G ( P , T , n ) = i = 1 N n i μ i ( P , T , n ) ,
where by definition, μ i = G n i is the chemical potential of i-th component. The fundamental Gibbs equation
d G = S d T + V d P + i = 1 N μ i d n i
describes the infinitesimal changes in the state of the system.
Consider an isolated system maintained at thermodynamic equilibrium without chemical reactions, and assume that its components coexist in two phases denoted by superscripts I and I I . Then G = G I + G I I . Equations (1) and (2) are valid for each phase, while the equilibrium condition d G = 0 reads d G I + d G I I = 0 . Moreover, since n t o t = i = 1 N n i = c o n s t and n i = n i I + n i I I = c o n s t , it follows that d n i I = d n i I I . Then, the equilibrium assumption implies the equality between the pressure and the temperature in two phases: T I = T I I , P I = P I I , as well as the equality of chemical potentials of each component in both phases:
μ i I ( P , T , n ) = μ i I I ( P , T , n ) . i = 1 , , N
In the remaining part of this paper, only the isobaric–isothermic conditions will be considered, and thus, the dependence of G on P and T will be neglected. It is also worth remarking that two different thermodynamic models might be used for the expressions of Gibbs energy and chemical potentials of two phases, and other types of diagrams (liquid–solid, solid–solid, etc.) can be modelled in this way. The computational method discussed in this paper can be easily adopted to this case, but for the sake of simplicity, it is assumed that both phases of the system are described by one single model of G.

2.2. Phase Coexistence Conditions in Partial Molar Variables

Most of the thermodynamic models of real mixtures, as well as the available data of phase separation, are given in terms of either mole, volume or mass fractions of the components. So, for practical purposes, it is more convenient to express the phase coexistence conditions with respect to one of these sets of variables, for example, the mole fractions. This choice makes no restriction on the computations presented below, but in the case of volume or mass fraction, the specific molar Gibbs energy should be replaced by Gibbs energy per unit of volume or mass.
Let g = G / n t o t and x i = n i n t o t denote, respectively, the molar free Gibbs energy and the mole fractions of the components of the mixture. Since i = 1 N x i = 1 , only N 1 of them can be used as independent variables. In what follows, we denote x = ( x 1 , , x N 1 ) , so that x N = 1 x 1 x N 1 . In terms of molar variables, the first N 1 equilibrium conditions (3) are equivalent to the following relations:
g I ( x I ) x i I = g I I ( x I I ) x i I I , i = 1 , , N 1 ,
while the N-th condition (3) yields
g I ( x I ) g I I ( x I I ) i = 1 N 1 g I ( x I ) x i I ( x i I x i I I ) = 0 .
The interested reader can find the detailed mathematical derivation of these conditions in Appendix A.
The phase coexistence conditions written in forms (4) and (5) admit a clear geometrical interpretation. Indeed, in the particular case N = 2 , Equations (4) and (5) reduce to the conditions of existence of a bitangent line to the graph of the molar free energy function g ( x ) , x [ 0 , 1 ] :
g ( x I ) = g ( x I I ) , g ( x I ) g ( x I I ) = g ( x I ) ( x I I x I ) ,
as it is shown in Figure 1. Here, g = g x . Other characteristics of the graph of g, like its convexity and the existence of inflexion points, are related to the material stability of the mixture. We will discuss them in a more general context in the next section.

2.3. Ternary Case: Phase Equilibrium Condition and Bitangent Planes Geometry

Let us now focus on a three-component liquid mixture whose components may coexist in one or two liquid phases. As in the binary case, in the ternary case, conditions (4) and (5) have a straightforward geometrical interpretation in terms of surface geometry.
Denote the composition domain of a ternary mixture by
Ω = def { x = ( x 1 , x 2 ) : x i [ 0 , 1 ] , i = 1 , 2 and x 1 + x 2 1 } .
Consider a smooth surface
W = def { ( x , z ) R 3 : z = g ( x ) , x Ω }
associated with the graph of the function g : Ω R . The vector field
N ( x , z ) = g 1 ( x ) x 1 + g 2 ( x ) x 2 z ,
defines the normal direction to this surface. Here, x i are the coordinate vector fields in R 3 and g i = g ( x ) x i , i = 1 , 2 . Let P 1 = ( x I , g ( x I ) ) and P 2 = ( x I I , g ( x I I ) ) be two points in R 3 belonging to the surface W, and such that
g 1 ( x I ) = g 1 ( x I I ) , g 2 ( x I ) = g 2 ( x I I ) .
In view of (6), it means that the normals N ( P 1 ) and N ( P 2 ) are collinear.
Further, the vector P 1 P 2 ¯ = ( x 1 I I x 1 I , x 2 I I x 2 I , g ( x I I ) g ( x I ) ) belongs to the tangent plane T P 1 W attached to the surface W at P 1 , which means that P 1 P 2 ¯ N ( P 1 ) , i.e.,
g ( x I I ) g ( x I ) g 1 ( x I ) ( x 1 I I x 1 I ) g 2 ( x I ) ( x 2 I I x 2 I ) = 0 .
In view of condition (7), the vector P 1 P 2 ¯ also belongs to T P 2 W . In other words, for N = 3 , conditions (4) and (5) mean that the surface W admits a bitangent plane passing through the points P 1 and P 2 . Such pairs of points on the surface W are called conodal. It is easy to see that the projections on Ω of the points P 1 and P 2 along the z-axis define the compositions of splitting liquid phases x I and x I I .
The projection of the bitangent segment P 1 P 2 on Ω is usually referred as node or tie-line. The curves on W formed by one-parametric families of conodal pairs define two directrices of a certain ruled surface in R 3 , with the bitangent segments P 1 P 2 being its generators. Projections of these directrices on Ω are called conodal or binodal curves of the phase diagram.

2.4. Differential Geometry of the Gibbs Energy Surface

As we have seen in the previous section, the geometry of the Gibbs energy surface W determines the topology of the underlying phase diagram. In fact, the most important properties of the phase diagram are encoded in the Gauss curvature of W. Denote by
H ( x ) = def 2 g ( x ) x i x j i , j = 1 N 1
the Hessian associated to the function g. Then, the Gauss curvature of surface W takes the form [17]:
K ( x ) = det H ( x ) ( g 1 ( x ) 2 + g 2 ( x ) 2 + 1 ) 2
According to the sign of K, the composition domain Ω can partitioned into elliptic ( K > 0 ), parabolic ( K = 0 ), and hyperbolic ( K < 0 ) sub-domains. The surface W is also subdivided into elliptic and hyperbolic sub-domains by the parabolic curve { ( x , z ) W : K ( x ) = 0 } . On the other hand, K is a product of two principal curvatures of the surface: K = k 1 · k 2 , so that the parabolic curve on W is the curve of zeros of one of the principal curvatures of the surface.
The vertical (i.e., along z-axis) projection of the parabolic curve on Ω defines the spinodal curve on the phase diagram. It bounds the elliptic sub-domain corresponding to the material stability domain of the mixture, which will be referred as Ω ˜ = def { x Ω : K ( x ) 0 } . It follows that the phase diagram can be described as the almost-Riemannian manifold M = ( Ω ˜ , H ) with a border, equipped with an almost-Riemannian metric associated with the Hessian H ( x ) . The sub-domain of Ω ˜ between the binodal and spinodal curves correspond to the metastable domain, while the remaining part delimited by the binodal curve and the boundary of Ω is the stable miscibility domain, where no phase separation occurs. In the binary case (see Figure 1), the binodal and spinodal curves reduce to two pairs of points that define the limits of stable and metastable domains.
The parabolic curves on smooth surfaces are very interesting geometrical objects. They may contain points where the curve has fourth-order contact with the tangent plane. Such special parabolic points (godrons in the terminology of [18]) correspond to the plait (critical) points of phase diagrams. It can be shown that special parabolic points give rise to two branches of conodal curves on W. In the underlying phase diagram, plait points are the only common points between spinodal and binodal curves. Besides plait points, binodal curves lie entirely in Ω ˜ . The plait point location can be found by solving equations due to Tompa [19]:
g 222 g 11 2 3 g 122 g 11 g 12 + 3 g 112 g 12 2 g 111 g 12 g 22 = 0 , g 11 g 22 g 12 2 = 0
where g i j k = 3 g x i x j x k
Figure 2 illustrates the described geometrical concepts.

3. Four-Dimensional Geometry of the Binodal Curve

So far, we have compared the 3D geometry of the surface W with the structure of the 2D phase diagram on Ω . But in order to better understand the intrinsic structure of the binodal curve, it would be more appropriate to associate a proper composition space Ω I and Ω I I to each of the phases. Consider now an extended configuration space Σ = Ω I × Ω I I of dimension four:
Σ = def { q R 4 : q = ( q 1 , q 2 ) with q 1 = x I Ω I , q 2 = x I I Ω I I } .
Equations (7) and (8) define three smooth sub-manifolds in Σ associated with zero-level sets of three functions F i : Σ R 1 such that
F 1 ( q ) = g 1 ( q 1 ) g 1 ( q 2 ) , F 2 ( q ) = g 2 ( q 1 ) g 2 ( q 2 ) ,
F 3 ( q ) = g ( q 2 ) g ( q 1 ) + x g ( q 1 ) | q 2 q 1 ,
where x g ( q 1 ) = ( g 1 ( q 1 ) , g 2 ( q 1 ) ) and ( | ) denote, respectively, the standard gradient and scalar product in R 2 . The intersection of these 3D sub-manifolds defines a one-dimensional sub-manifold B Σ such that
B = def { q Σ : F i ( q ) = 0 , i = 1 , 2 , 3 } .
The orthogonal projections π i : Σ Ω such that π I ( q ) = q 1 and π I I ( q ) = q 2 define two branches π I ( B ) and π I I ( B ) of the binodal curve. In what follows, the sub-manifold B will be referred to as the generalized binodal curve.
Let V ( q ) T q B be the tangent vector at point q B of the generalized binodal B defined by (11). By construction, V i = 1 3 K e r q F i , so that by definition,
q F i ( q ) · V ( q ) = 0 , i = 1 , 2 , 3 .
Computing the gradients of F i yields the following system of linear equations verified by components of the vector field V = ( V 1 , V 2 , V 3 , V 4 ) at q B :
g 11 ( q 1 ) g 12 ( q 1 ) g 11 ( q 2 ) g 12 ( q 2 ) g 12 ( q 1 ) g 22 ( q 1 ) g 12 ( q 2 ) g 22 ( q 2 ) Φ 1 ( q ) Ψ 1 ( q ) 0 0 · V 1 V 2 V 3 V 4 = 0 ,
where the functions Φ 1 and Ψ 1 are defined by the relation
Φ 1 Ψ 1 = H ( q 1 ) · q 1 q 2 .
Rewriting (12) in a more compact form yields
Φ 1 V 1 + Ψ 1 V 2 = 0 , H ( q 2 ) V 3 V 4 = H ( q 1 ) V 1 V 2 .
If det H ( q 2 ) 0 and at least one of the functions Φ 1 , Ψ 1 is non-zero, one obtains the solution of (13) in the form
V 1 = Ψ 1 det H ( q 2 ) , V 2 = Φ 1 det H ( q 2 )
and
V 3 V 4 = H ( q 2 ) 1 · H ( q 1 ) · V 1 V 2 ,
Performing all necessary simplifications, we obtain
V 1 V 2 V 3 V 4 = Ψ 1 det H ( q 2 ) Φ 1 det H ( q 2 ) Ψ 2 det H ( q 1 ) Φ 2 det H ( q 1 )
where, by definition,
Φ 2 Ψ 2 = H ( q 2 ) · q 1 q 2 .
The above expressions define a smooth vector field V T Σ . It is well defined except for the singular points q such that det H ( q 1 ) = det H ( q 2 ) = 0 , in other words, if q 1 and q 2 belong to the spinodal, i.e., if they coincide forming the plait point of the phase diagram. The described geometrical construction is shown in Figure 3.
Summing up, we showed that the binodal curve can be easily recovered from the projection of the integral curve of the vector field V. Moreover, the particular structure of formula (14), namely the properties (13) and (15), imply that the tie-lines are orthogonal to the binodal with respect the metric in Ω ˜ associated with the Hessian matrix H ( x ) .

Numerical Computation of Binodal and Spinodal Curves

Being the integral curve of the vector field V, the generalized binodal curve can be computed by solving the system of ordinary differential equations q ˙ = V ( q ) in Σ . Here, the dot notation stands for the derivative with respect to any suitable parameter. This result has an important practical application regarding its computation.
Using the vector field V, the numerical computation of binodal curves reduces to a simple ODE integration by any conventional solver with desired accuracy. The normalization of V allows one to avoid the eventual stiffness problem when approaching the border of Ω ˜ , so it is recommended to use the arc-length parameter for the integration. To start the computation, one needs to find an initial point in Σ , i.e., a starting tie-line of the binodal curve. This can be performed by analysing the profile of the W surface along the boundaries of the composition domain Ω , in other words, by solving at most three binary problems over the interval [ 0 , 1 ] . The only exception here are the closed-type 0 phase diagrams, for which the initial tie-line must be found inside Ω .
The same method can be applied to derive the differential equation of the spinodal curve. This case is even simpler because it is a problem in the 2D space Ω . Indeed, being the solution of the equation H ( x ) = 0 , the spinodal curve is a solution to the differential equation associated with the vector field S ( x ) = x H ( x ) in Ω , which, in general, is regular at plait points. Here, the superscript denotes the orthogonal complement to the vector x H ( x ) in R 2 . In other words, the spinodal curve is the solution of the ODE equation x ˙ = S ( x ) , x Ω ; the starting point for the integration can be detected by finding one inflexion point of the function g on the boundary of Ω or inside Ω for closed-type curves.
Once binodal and spinodal curves are computed, the plait point can be found by solving numerically Equation (9), taking the approximate common point between these curves as the initial guess for the non-linear solver in order to facilitate the converge procedure.
The implementation of the described computational method requires accessing the derivatives of the thermodynamic model defining the function g up to the order two. The finite-difference method may be not sufficient to obtain the necessary accuracy level for the Hessian computation. In fact, this is the main obstacle in the implementation of such types of algorithms in many other domains involving the thermodynamic models of real systems. However, we would like to stress the fact that the nowadays numerical technology allows one to easily overcome this inconvenience. For the academic use, the symbolic computation packages, like Maple or Mathematica, which we used in this paper, would be more than sufficient. To develop a stand-alone calculator of phase diagram, the automatic differentiation technique can be employed. A pilot version of a functional code of this kind was described in [20] for the computation of the univolatility curves of the residue curve maps.

4. Inverse Problem and Case Studies

Following [16], ternary LLE phase diagrams of biphasic systems can be roughly divided into three main classes depending on the number of the partially miscible binary pairs:
  • Type 0 or “island” type: The diagram is characterized by a closed heterogeneous domain inside the composition triangle, while all three binary pairs are miscible. The systems of this type exhibit two plait points.
  • Type I: One pair of components exhibits a miscibility gap on the border of the composition triangle. This type of diagram possesses one plait point where both liquid phases have the same composition. This is the most common type of phase diagrams (75% according to [21]).
  • Type II: This type is characterized by the presence of two partial miscibility gaps on the borders of the composition triangle. Such diagrams do not have plait points. They represent about 20% of known solutions ([21]).
Under variations in temperature or pressure, all these types of behaviour transform one into another, or split into two or more heterogeneous zones. Two heterogeneous zones can also melt, forming three phase domains. The analysis of this phenomenon goes beyond the aim of this paper. The three case studies presented below correspond to the three standard types of diagrams, using the Flory–Huggins model for the free energy function g.

4.1. The Flory–Huggins Model Equation

The choice of an appropriate thermodynamic model for function g depends on the particular application. For LLE diagrams, NRTL or UNIQIUAC models that describe a huge class of ternary systems, were studied by numerous authors, though the quality of the model parameter regression is not always satisfactory. As suggested by [1], for the polymer solutions, the Flory–Huggins (FH) model can be employed to describe the systems of type non-solvent–solvent–polymer, solvent–polymer–polymer, and even three polymers. The relative simplicity of the mathematical expression of constitutive equation is the main advantage of the FH model with respect to NRTL or UNIQUAC models used by other authors. Unlike NRTL or UNIQUAC, formulated in terms of mole fractions, the FH model uses partial volume fractions, assuming that the total volume of the systems is equal to the sum of partial volumes.
The classical Flory–Huggins model defines the free energy of mixing per unit of volume according to the expression
g = x 1 N 1 ln x 1 + x 2 N 2 ln x 2 + x 3 N 3 ln x 3 + χ 12 x 1 x 2 + χ 13 x 1 x 3 + χ 23 x 2 x 3 + β x 1 x 2 x 3 .
Here x i , i = 1 , 2 , 3 are the volume fractions of the components, and N 2 and N 3 and relative numbers of segments in the molecules of the compounds with respect to the first compound, which usually corresponds to water; thus, N 1 = 1 . The symbols χ i j and β stay for the binary and ternary interaction coefficients. In all computation below, x 3 was replaced by x 3 = 1 x 1 x 2 .
In most of the papers that use Equation (16), the ternary coefficient β = 0 , while some of the binary interaction coefficients may depend on the components of the composition vector x. However, the triple interaction term accounts for the possible linear dependency of χ i j on x, and allows one to relax the total volume conservation hypothesis when mixing the components. So, in this paper, to maintain the model equations as simple as possible, the six parameters N i , χ i j , and β are assumed to be constant. In other words, in this paper, we propose considering (16) as a formal mathematical expression without taking into account the exact physical meaning of its parameters, at least in the cases where constant volume assumption is difficult to justify.

4.2. Parameter Estimation Procedure and Case Studies

In order to calculate the LLE diagrams presented below, a set of n b tie-line measurements was used.
In order to fit the model, the six parameters of the FH model were computed via a non-linear optimization procedure, associated with the function
min p k = 1 n b F 1 2 ( q k ) + F 2 2 ( q k ) + F 3 2 ( q k ) ( 1 F 3 2 ( q k ) ) 2 ,
where F i are defined by (10), and p is the vector of unknown scalar parameters. The denominator term penalizes those pairs of points that do not belong to the same bitangent line. In our examples this criterion showed a better convergence to the experimental data comparing to the sole square term in the nominator of (17). The exact mathematical expressions of functions F i associated with the FH model are provided in Appendix B.1.
The described NLP minimization procedure was applied to compute two of the examples below: the type 0 diagram water–dimethyl sulfoxide (DMSO)–tetra hydrofuran (THF) (see in Figure 4) and water–phenol–acetone at 50 °C (Figure 5a). However, for the diagrams of types I and II, the computation can be simplified, provided the measurement of at least one of binary miscibility gaps is available. In fact, in this case, the linearity of the FH model with respect to the parameters allows one to compute the part of the parameters describing the binary pair straightforwardly, and thus reduces the number of unknowns of problem (17). The computation formulae are provided in Appendix B.1. Finally, the quality of the predicted model can be evaluated using the standard formula for the root mean square deviation (RMSD)
σ = 100 1 6 n k = 1 n b i 3 x i , k I , F H x i , k I , e x p 2 + x i , k I I , F H x i , k I I , e x p 2 1 / 2
used by other authors.
All numerical computations presented in this paper were realized by Mathematica 9 software [22]. We used the standard functions with the default choice of methods. Namely, the ODE integration was performed by NDSolve function, which implements the LSODA algorithm (combined Adams and BDF methods). The optimization problem was solved by the FindMinimum command with the quasi-Newton method chosen by default. The original unpublished data used for the LLE regression are included in Appendix C, together with the computational formulae for the Flory–Huggins model available in Appendix B.1.
Figure 4. Phase diagram of water–DMSO–THF in the volume fraction space. Here, black points and dashed lines correspond to the experimental points and tie-lines from [23]. The blue curve and the red dashed curve represent, respectively, the binodal and spinodal curves calculated using the Flory–Huggins model. White points indicate the end-points of the FH-calculated tie-lines (grey) used in the RMSD evaluation.
Figure 4. Phase diagram of water–DMSO–THF in the volume fraction space. Here, black points and dashed lines correspond to the experimental points and tie-lines from [23]. The blue curve and the red dashed curve represent, respectively, the binodal and spinodal curves calculated using the Flory–Huggins model. White points indicate the end-points of the FH-calculated tie-lines (grey) used in the RMSD evaluation.
Entropy 25 01329 g004

4.2.1. Type 0 Diagram: Water–DMSO–THF

The parametric regression of type 0 diagrams is known to be a difficult task. Figure 4 represents our result of the phase diagram prediction based on the experimental data available in [23] for the water–DMSO–THF system. These data were also studied in [8,24,25] using NRTL and UNIQUAC models. In Table 1, the computed RMSD value of the reconstructed FH model is compared with those obtained by other authors, showing a good quality of prediction for this mixture using the FH model.

4.2.2. Type I Diagram: Water–Phenol–Acetone

The ternary diagrams of the water–phenol–acetone system at 50 °C and 60 °C using NRTL and UNIQUAC models were studied in [24,26], both using the same data from [26]. Our result of the FH regression according the method described above is shown in Figure 5a). In Figure 5b, we show the phase diagram for 56 °C computed using the data from [23]. In this second case, the computation was simplified because [23] provides the measurement on the phenol–water miscibility gap. It allows one to compute explicitly the values of χ 13 and N 3 . The remaining four parameters were estimated via the general optimization procedure described above. Again, the comparison with the RMSD values published by other authors shows a very good quality of our result for 50 °C, and the best one for 56 °C.
Figure 5. Phase diagram of water–phenol–acetone in the volume fraction space at 50 °C (a) and at 56 °C (b). Black points and dashed lines correspond to the experimental points and tie-lines from [26] (left) and [23]. The blue curve and the red dashed curve represent, respectively, the binodal and spinodal curves calculated using the Flory–Huggins model. White points indicate the end-points of the FH-calculated tie-lines (grey) used in the RMSD evaluation.
Figure 5. Phase diagram of water–phenol–acetone in the volume fraction space at 50 °C (a) and at 56 °C (b). Black points and dashed lines correspond to the experimental points and tie-lines from [26] (left) and [23]. The blue curve and the red dashed curve represent, respectively, the binodal and spinodal curves calculated using the Flory–Huggins model. White points indicate the end-points of the FH-calculated tie-lines (grey) used in the RMSD evaluation.
Entropy 25 01329 g005

4.2.3. Type II Diagram: Water–Acetone–Hexadecane

For the water–acetone–hexadecane system, we used our own set of 14 tie-line measurements obtained with Raman’s spectroscopy. Another set of 17 measurements of the phase envelope, without taking into account the conodal pair correspondence, were obtained using the redisolution by adding acetone. More details on the experimental procedure are available in [27]; the data are included in Appendix C, Table A1 and Table A2. The phase envelope measurements allowed one to detect the acetone–hexadecane miscibility gap. For the miscibility gap of water–hexadecane, we used the data accessible in the literature. Knowing these two miscibility gaps, we first computed the four parameters N 2 , N 3 , χ 13 , and χ 23 . Due to the significant variation of the solubility data of hexadecane in water in different sources, the term χ 13 was replaced by the χ 13 + δ , where the correction term δ was adjusted in the next step of computation. After this important simplification, the number of unknown parameters was reduced to three. Applying the optimization protocol on the remaining set of parameters produced the result depicted in Figure 6.
This result shows an excellent prediction for shape of phase envelope, and a very good quality for the tie-lines’ directions. The RMSD value associated with this predictions equal to 2.92 .

5. Conclusions

This work showed that the challenging computation of liquid–liquid phase diagrams can be drastically simplified by taking into account the intrinsic geometrical structures associated with phase equilibrium conditions. The crucial mathematical idea is to work in an adequate extended space. Indeed, we show that the stability domain boundary defined by the 2D binodal curve is a projection of a higher-dimensional object, the 4D generalized binodal curve, which can be easily computed by solving a set of ordinary differential equations.
This mathematical viewpoint is employed to propose a new numerical algorithm for phase-diagram computation, which drastically reduces the computation effort and guarantees a high value of accuracy. Notably, the only iterative procedure is the one used to find the initial point for the integration. We tested this methodology on four ternary liquid phase diagrams of different topological types. We chose a modified Flory–Huggins model for the regression procedure applied to the available experimental data. The presented results show a good, sometimes excellent, accordance between the data and the calculated model, which opens a promising perspective for the further development of the method to study other types of multiphase diagrams. Furthermore, the developed approach can be used for any other type of ternary diagram and is valid for any chosen thermodynamic model.

Author Contributions

N.S.: conceptualization, methodology, software, validation, writing—original draft; V.G.: supervision, resources, writing—reviewing and editing; K.R.: validation, investigation, resources, data curation, writing—reviewing and editing. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The original experimental data are included in Appendix C.

Acknowledgments

The authors are grateful to Lison Raynal for the acquisition of the experimental data.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Phase Coexistence Conditions in Terms of Molar Fractions

Let g = G / n t o t denote the molar free Gibbs energy. Since the mole fractions x i = n i n t o t verify i = 1 N x i = 1 , only N 1 of them can be used as independent variables. In what follows, we denote x = ( x 1 , , x N 1 ) , so that n N = n t o t ( 1 x 1 . . . x N 1 ) . Since d x N = i = 1 N 1 d x i , the Gibbs Equation (2) takes the form
d g = s d T + v d P + i = 1 N 1 ( μ i μ N ) d x i ,
where s = S / n t o t and v = V / n t o t are the specific molar entropy and volume, correspondingly. Moreover, for all i = 1 , , N 1 we obtain
g x i = 1 n G x i = 1 n G n i n i x i + G n N n N x i = μ i μ N ,
and hence,
g = 1 n i = 1 N μ i n i = i = 1 N 1 ( μ i μ N ) x i + μ N .
Rewriting (A2), (A3) for each phase, the phase coexistence Equations (3) can be expressed in terms of mole fractions. Indeed, the first N 1 equilibrium conditions (3) are equivalent to the following relations:
g I ( x I ) x i I = g I I ( x I I ) x i I I , i = 1 , , N 1 ,
while the N-th condition (3) yields
g I ( x I ) g I I ( x I I ) i = 1 N 1 g I ( x I ) x i I ( x i I x i I I ) = 0 .

Appendix B. The Flory–Huggins Model Computational Formulae

Appendix B.1. Phase Coexistence Conditions Associated with Flory–Huggins Model and Linearity

In the case of the Flory-Huggins model (16) the mathematical expressions of functions F i take the forms
F 1 = χ 12 x 2 I x 2 I I χ 13 2 ( x 1 I x 1 I I ) + x 2 I x 2 I I χ 23 x 2 I x 2 I I
+ β x 2 I 2 x 1 I x 2 I ( x 2 I ) 2 x 2 I I + 2 x 1 I I x 2 I I + ( x 2 I I ) 2
+ 1 N 1 ln x 1 I x 1 I I 1 N 3 ln 1 x 1 I x 2 I 1 x 1 I x 2 I I
F 2 = χ 12 x 1 I x 1 I I χ 13 x 1 I x 1 I I χ 23 x 1 I x 1 I I + 2 x 2 I 2 x 2 I I +
+ β x 1 I ( x 1 I ) 2 x 1 I I + ( x 1 I I ) 2 2 x 1 I x 2 I + 2 x 1 I I x 2 I I
+ 1 N 2 ln x 2 I x 2 I I 1 N 3 ln 1 x 1 I x 2 I 1 x 1 I I x 2 I I
F 3 = χ 12 ( x 1 I x 1 I I ) ( x 2 I x 2 I I ) χ 13 ( x 1 I x 1 I I ) ( x 1 I x 1 I I + x 2 I x 2 I I )
χ 23 ( x 2 I x 2 I I ) ( x 1 I x 1 I I + x 2 I x 2 I I )
+ β x 1 I x 2 I 2 ( x 1 I ) 2 x 2 I x 1 I I x 2 I + 2 x 1 I x 1 I I x 2 I 2 x 1 I ( x 2 I ) 2 + x 1 I I ( x 2 I ) 2
x 1 I x 2 I I + ( x 1 I ) 2 x 2 I I + x 1 I I x 2 I I ( x 1 I I ) 2 x 2 I I + 2 x 1 I x 2 I x 2 I I x 1 I I ( x 2 I I ) 2
+ 1 N 1 x 1 I x 1 I I x 1 I I ln x 1 I x 1 I I + 1 N 2 x 2 I x 2 I I x 2 I I ln x 1 I x 1 I I
1 N 3 x 1 I x 1 I I + x 2 I x 2 I I ( 1 x 1 I I x 2 I I ) ln 1 x 1 I x 2 I 1 x 1 I I x 2 I I
The phase coexistence conditions (4) and (5) take the form F i ( q ) = 0 for i = 1 , 2 , 3 and q = ( x 1 I , x 2 I , x 1 I I , x 2 I I ) Σ . Remarkably, for the FH model, these conditions are linear with respect to the interaction parameters χ i j , β and to the 1 / N 2 and 1 / N 3 values. This fact permits the direct computation of these constants in the binary case, as it is explained in the next section.

Appendix B.2. Detection of the Parameters of the Fh Model from the Miscibility Gap Limits

Assume that components 1 and 3 are partially miscible and assume that the values x 1 I and x 1 I I are available through the experimental procedure. The binary mixture of components 1–3 corresponds to the side x 2 = 0 of the boundary of Ω . Setting x 1 = x , x 3 = 1 x , the function g can be rewritten in the form
g 1 = x N 1 ln x + 1 x N 3 ln ( 1 x ) + χ 13 x ( 1 x ) .
Using expressions (A6), the phase coexistence conditions (4) and (5) applied to the function g 1 then have the form
2 ( x I x I I ) χ 13 + 1 N 3 ln 1 x I 1 x I I ln x I x I I = 0 ,
x I x I I χ 13 ( x I x I I ) 2 + ( 1 x I I ) ln 1 x I 1 x I I
1 N 3 x I x I I x I I ln x I x I I = 0 .
If the values of x I and x I I are known, the above equations form a system of two linear algebraic equations with respect to χ 13 and 1 / N 3 , which can be easily computed explicitly.
Assume that both couples 1–3 and 2–3 are partially miscible. Then, the same procedure as above can be applied to find parameters N 2 , N 3 , χ 13 , and χ 23 by solving a system of four linear algebraic equations in terms of the miscibility gap limits x 1 I , x 1 I I and x 2 I , x 2 I I .

Appendix C. Water–Acetone–Hexadecane Experimental Data

Table A1. Water–acetone–hexadecane tie-lines, mass fraction experimental data.
Table A1. Water–acetone–hexadecane tie-lines, mass fraction experimental data.
Phase IPhase II
WaterAcetoneHexadecaneWaterAcetoneHexadecane
0.0050260.0880510.9069230.11532870.880.0046713
0.0029560.1217630.8752810.073896690.91150.014603313
0.0043190.1277640.8679170.050619990.9190.030380015
0.0007750.1532730.8459520.026046920.90750.06645308
0.0003130.2845720.7151150.010038130.8650.124961871
0.0021870.0235470.9742650.938970.061020.000000001
0.0021480.0294350.9684170.877770.122230.000000001
0.0019490.0368010.961250.7930040.206990.000000001
0.0003080.0486580.9510340.674130.32590.000000001
0.0030930.0544030.9425030.563580.43640.000000001
0.0002760.0653210.9344030.495480.504520.000000001
0.0052980.0749690.9197330.36550.634460.000000001
0.0030660.0823430.9145910.25410.745850.000000001
0.0023070.0971210.9005730.17470.825240.000000001
Table A2. Water–acetone–hexadecane phase envelope, mass fraction experimental data. The highlighted values correspond to the miscibility gap of the binary mixture acetone–hexadecane.
Table A2. Water–acetone–hexadecane phase envelope, mass fraction experimental data. The highlighted values correspond to the miscibility gap of the binary mixture acetone–hexadecane.
WaterAcetoneHexadecane
0.35030.64968.98 × 10−5
0.19890.79960.0015
0.15190.8460.002
0.1150.8820.0033
0.100560650.892693890.00674546
0.061228250.916498760.02227299
0.041123570.91780710.04106932
0.027409970.908410790.06417925
0.014453560.884302020.10124442
0.00390.84560.1504
0.0020.7990.199
00.7830.217
00.2640.736
0.00050.19750.8021
0.00070.15890.8404
0.000920.11220.8869
0.000650.06030.9391

References

  1. Koningsveld, R.; Stockmayer, W.H.; Nies, E. Polymer Phase Diagrams; Oxford Univ. Press: Oxford, UK, 2001. [Google Scholar]
  2. Ferrari, J.C.; Nagatani, G.; Corazza, F.C.; Oliveira, J.V.; Corazza, M.L. Application of stochastic algorithms for parameter estimation in the liquid–liquid phase equilibrium modeling. Fluid Phase Equilib. 2009, 280, 110–119. [Google Scholar] [CrossRef]
  3. Vatani, M.; Asghari, M.; Vakili-nejhaad, G. Application of Genetic Algorithm to Parameter Estimation in Liquid-liquid Phase Equilibrium Modeling. J. Math. Comp. Sci. 2012, 5, 60–66. [Google Scholar] [CrossRef]
  4. Fernandez-Vargas, J.A.; Bonilla-Petriciolet, A.; Segovia-Hernandez, J.G. An improved ant colony optimization method and its application for the thermodynamic modeling of phase equilibrium. Fluid Phase Equilib. 2013, 353, 121–131. [Google Scholar] [CrossRef]
  5. Zhang, H.; Bonilla-Petriciolet, A.; Rangaiah, G.P. A Review on Global Optimization Methods for Phase Equilibrium Modeling and Calculations. Open Thermodyn. J. 2011, 5, 71–92. [Google Scholar] [CrossRef]
  6. Aliena, F.W.; Smolders, C.A. Calculation of Liquid-Liquid Phase Separation in a Ternary System of a Polymer in a Mixture of a Solvent and a Nonsolvent. Macromolecules 1982, 15, 1491–1497. [Google Scholar]
  7. Michelsen, M.L. The isothermal flash problem. Part I. Stability. Fluid Phase Equilibria 1982, 9, 1–19. [Google Scholar]
  8. Binous, H.; Bellagi, A. Calculation of ternary liquid-liquid equilibrium data using arc-length continuation. Wiley Eng. Rep. 2021, 3, e12296. [Google Scholar]
  9. Bausa, J.; Marquardt, W. Quick and reliable phase stability test in VLLE flash calculations by homotopy continuation. Comp. Chem. Eng. 2000, 24, 2447–2456. [Google Scholar] [CrossRef]
  10. Levelt, S. How Fluids Unmix; Royal Netherlands Academy of Arts and Sciences: Amsterdam, The Netherlands, 2002. [Google Scholar]
  11. Arnold, V.I. Singulatrities of Caustics and Wave Fronts; Kluwer: Baltimore, MD, USA, 1990. [Google Scholar]
  12. Varchenko, A.N. Evolutions of convex hulls and phase transitions in thermodynamics. J. Math. Sci. 1990, 52, 3305–3325. [Google Scholar] [CrossRef]
  13. Aicardi, F.; Valentin, P.; Ferrand, E. On the classification of generic phenomena in one-parameter families of binary mixtures. Phys. Chem. Chem. Phys. 2002, 4, 884–895. [Google Scholar] [CrossRef]
  14. Gaite, J.A. Phase transitions as catastrophes: The tricritical point. Phys. Rev. A 1990, 41, 5320–5324. [Google Scholar] [CrossRef] [PubMed]
  15. Gaite, J.; Margalef-Roig, J.; Miret-Atrtes, S. Analysis of a three-component model phase diagram by catastrophe theory. Phys. Rev. B 1998, 57, 13527. [Google Scholar] [CrossRef]
  16. Treybal, R.E. Liquid Extraction, 2nd ed.; McGraw-Hill: New-York, NY, USA, 1963. [Google Scholar]
  17. Do Carmo, M. Differential Geometry of Curves and Surfaces, 2nd ed.; Courier Dover Publications: Mineola, NY, USA, 2017. [Google Scholar]
  18. Uribe-Vargas, R. A projective invariant for swallowtails and godrons, and global theorems on the flecnodal curve. Mosc. Math. J. 2006, 6, 731–768. [Google Scholar] [CrossRef]
  19. Tompa, H. Polymer Solutions; Butterworths: London, UK, 1956. [Google Scholar]
  20. Cots, O.; Shcherbakova, N.; Gergaud, J. SMITH: Differential homotopy and automatic differentiation for computing thermodynamic diagrams of complex mixtures. Comput. Aided Chem. Eng. 2021, 50, 1081–1086. [Google Scholar]
  21. Gmehling, J.; Kolbe, B.; Kleiber, M.; Rarey, J. Chemical Thermodynamics for Process Simulation, 2nd ed.; Willey-VCH Verlag: Weinheim, Germany, 2019. [Google Scholar]
  22. Wolfram Mathematica: The World’s Definitive System for Modern Technical Computing. Available online: https://www.wolfram.com/mathematica/ (accessed on 1 January 2020).
  23. Sørensen, J.M.; Arlt, W. Liquid-Liquid Equilibrium Data Collection, in der Reihe: Dechema Chemistry Data Series; DECHEMA: Frankfurt, Germany, 1980; Volume V. [Google Scholar]
  24. Zuber, A.; Raimundo, R.; Mafra, M.R.; Filho, L.C.; Oliveira, J.V.; Corazza, M.L. Thermodynamic modeling of ternary liquid-liquid systems with forming immiscibility islands. Braz. Arch. Biol. Technol. 2013, 56, 1034–1042. [Google Scholar] [CrossRef]
  25. Olaya, M.M.; Reyes-Labarta, J.A.; Velasco, R.; Ibarra, I.; Marcilla, A. Modelling liquid–liquid equilibria for island type ternary systems. Fluid Phase Equilibria 2008, 265, 184–191. [Google Scholar] [CrossRef]
  26. Mafra, M.R.; Krähenbühl, M.A. Liquid-Liquid Equilibrium of (Water + Acetone) with Cumene or r-Methylstyrene or Phenol at Temperatures of (323.15 and 333.15) K. J. Chem. Eng. Data 2006, 51, 753–756. [Google Scholar] [CrossRef]
  27. Roger, K.; Raynal, L.; Shcherbakova, N. Nanoprecipitation through solvent-shifting using rapid mixing: Dispelling the Ouzo boundary to reach large solute concentrations. J. Colloid Interface Sci. 2023, 650, 2049–2055. [Google Scholar] [CrossRef] [PubMed]
Figure 1. Two-phase separation in binary mixture: geometrical meaning of conditions (4) and (5).
Figure 1. Two-phase separation in binary mixture: geometrical meaning of conditions (4) and (5).
Entropy 25 01329 g001
Figure 2. Gibbs energy surface W vs. the corresponding phase diagram. The hyperbolic domain of surface W is delimited by the parabolic curve (black dashed curve). Its vertical projection on Ω defines the spinodal curve (black and white dashed curve), which bounds the unstable domain of the phase diagram. The one-parametric family of conodal pairs of points on W (co-nodal directrix curve, yellow) projects on the binodal curve of the phase diagram. It can touch the spinodal curve at a plait point. The binodal curve divides the elliptic domain on Ω into stable (white) and metastable (light blue) sub-domains.
Figure 2. Gibbs energy surface W vs. the corresponding phase diagram. The hyperbolic domain of surface W is delimited by the parabolic curve (black dashed curve). Its vertical projection on Ω defines the spinodal curve (black and white dashed curve), which bounds the unstable domain of the phase diagram. The one-parametric family of conodal pairs of points on W (co-nodal directrix curve, yellow) projects on the binodal curve of the phase diagram. It can touch the spinodal curve at a plait point. The binodal curve divides the elliptic domain on Ω into stable (white) and metastable (light blue) sub-domains.
Entropy 25 01329 g002
Figure 3. Four-dimensional lifting of the binodal curve. The complete configuration space of the binodal curve is Σ = Ω I × Ω I I . The vector field V = ( V 1 , V 2 , V 3 , V 4 ) generates a smooth curve B (blue) in Σ . The orthogonal projections of B onto Ω form the binodal curve (black) on the underlying phase diagram. The dotted lines indicate the coupled pairs of phases, and the dashed red curve indicates the location of the spinodal curve.
Figure 3. Four-dimensional lifting of the binodal curve. The complete configuration space of the binodal curve is Σ = Ω I × Ω I I . The vector field V = ( V 1 , V 2 , V 3 , V 4 ) generates a smooth curve B (blue) in Σ . The orthogonal projections of B onto Ω form the binodal curve (black) on the underlying phase diagram. The dotted lines indicate the coupled pairs of phases, and the dashed red curve indicates the location of the spinodal curve.
Entropy 25 01329 g003
Figure 6. Water–acetone–hexadecane system in the volume fraction space. Black points and dashed lines correspond to the experimental points and tie-lines. The blue curve and the red dashed curve represent, respectively, the binodal and spinodal curves calculated using the Flory–Huggins model. White points indicate the end-points of the FH-calculated tie-lines (grey) used in the RMSD evaluation.
Figure 6. Water–acetone–hexadecane system in the volume fraction space. Black points and dashed lines correspond to the experimental points and tie-lines. The blue curve and the red dashed curve represent, respectively, the binodal and spinodal curves calculated using the Flory–Huggins model. White points indicate the end-points of the FH-calculated tie-lines (grey) used in the RMSD evaluation.
Entropy 25 01329 g006
Table 1. Root mean square deviation σ of analysed mixture.
Table 1. Root mean square deviation σ of analysed mixture.
SystemModelT, C σ
water–DMSO–THFFH (this work)20°3.53
water–DMSO–THF, [25]NRTL20°5.83
water–DMSO–THF, [8]UNIQUAC20°3.4
water–DMSO–THF, [24]NRTL20°3.18, 3.09
water–DMSO–THF, [24]UNIQUAC20°3.97, 3.4
water–acetone–phenolFH (this work)50°0.88
water–acetone–phenol, [26]NRTL50°0.81
water–acetone–phenol, [24]NRTL50°1.13
water–acetone–phenol, [24]UNIQIAC50°1.17
water–acetone–phenolFH (this work)56°1.27
water–acetone–phenol, [23]NRTL56°1.61
water–acetone–phenol, [23]UNIQUAC56°1.52
water–acetone–hexadecaneFH (this work)20°2.92
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Shcherbakova, N.; Gerbaud, V.; Roger, K. Using the Intrinsic Geometry of Binodal Curves to Simplify the Computation of Ternary Liquid–Liquid Phase Diagrams. Entropy 2023, 25, 1329. https://doi.org/10.3390/e25091329

AMA Style

Shcherbakova N, Gerbaud V, Roger K. Using the Intrinsic Geometry of Binodal Curves to Simplify the Computation of Ternary Liquid–Liquid Phase Diagrams. Entropy. 2023; 25(9):1329. https://doi.org/10.3390/e25091329

Chicago/Turabian Style

Shcherbakova, Nataliya, Vincent Gerbaud, and Kevin Roger. 2023. "Using the Intrinsic Geometry of Binodal Curves to Simplify the Computation of Ternary Liquid–Liquid Phase Diagrams" Entropy 25, no. 9: 1329. https://doi.org/10.3390/e25091329

APA Style

Shcherbakova, N., Gerbaud, V., & Roger, K. (2023). Using the Intrinsic Geometry of Binodal Curves to Simplify the Computation of Ternary Liquid–Liquid Phase Diagrams. Entropy, 25(9), 1329. https://doi.org/10.3390/e25091329

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop