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

Nataliya Shcherbakova, Vincent Gerbaud, Kevin Roger

► To cite this version:

Nataliya Shcherbakova, Vincent Gerbaud, Kevin Roger. Using the Intrinsic Geometry of Binodal Curves to Simplify the Computation of Ternary Liquid-Liquid Phase Diagrams. Entropy, 2023, 25 (9), pp.1329. 10.3390/e25091329. hal-04277815

HAL Id: hal-04277815

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

Submitted on 9 Nov 2023

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Article

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

Nataliya Shcherbakova *,† ID, Vincent Gerbaud † ID and Kevin Roger †

Laboratoire de Génie Chimique, Université de Toulouse, CNRS, INP, UPS, 31432 Toulouse, France;
vincent.gerbaud@ensiacet.fr (V.G.); kevin.roger@cnrs.fr (K.R.)

* Correspondence: nataliya.shcherbakova@ensiacet.fr

† These authors contributed equally to this work.

**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.

**Keywords:** LLE ternary diagrams; binodal curves; spinodal curves; differential path-following method; ordinary differential equations

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

check for
updates

**Citation:** 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

Academic Editor: Antony N. Beris

Received: 29 June 2023

Revised: 30 August 2023

Accepted: 11 September 2023

Published: 13 September 2023

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

**Copyright:** © 2023 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article
distributed under the terms and
conditions of the Creative Commons
Attribution (CC BY) license (https://
creativecommons.org/licenses/by/
4.0/).

## 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–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 $\mathbf{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, \mathbf{n})$. Being a homogeneous function of the first order with respect to $n_i$, $G$ can be expressed as

$$G(P, T, \mathbf{n}) = \sum_{i=1}^{N} n_i \mu_i(P, T, \mathbf{n}), \qquad (1)$$

where by definition, $\mu_i = \frac{\partial G}{\partial n_i}$ is the chemical potential of i-th component. The fundamental Gibbs equation

$$dG = -SdT + VdP + \sum_{i=1}^{N} \mu_i dn_i \qquad (2)$$

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 II. Then $G = G^I + G^{II}$. Equations (1) and (2) are valid for each phase, while the equilibrium condition $dG = 0$ reads $dG^I + dG^{II} = 0$. Moreover, since $n_{tot} = \sum_{i=1}^{N} n_i = \text{const}$ and $n_i = n_i^I + n_i^{II} = \text{const}$, it follows that $dn_i^I = -dn_i^{II}$. Then, the equilibrium assumption implies the equality between the pressure and the temperature in two phases: $T^I = T^{II}$, $P^I = P^{II}$, as well as the equality of chemical potentials of each component in both phases:

$$\mu_i^I(P, T, \mathbf{n}) = \mu_i^{II}(P, T, \mathbf{n}). \quad i = 1, \dots, N \qquad (3)$$

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_{tot}$ and $x_i = \frac{n_i}{n_{tot}}$ denote, respectively, the molar free Gibbs energy and the mole fractions of the components of the mixture. Since $\sum_{i=1}^{N} x_i = 1$, only $N-1$ of them can be used as independent variables. In what follows, we denote $\mathbf{x} = (x_1, \dots, x_{N-1})$, so that $x_N = 1 - x_1 - \dots - x_{N-1}$. In terms of molar variables, the first $N-1$ equilibrium conditions (3) are equivalent to the following relations:

$$\frac{\partial g^I(\mathbf{x}^I)}{\partial x_i^I} = \frac{\partial g^{II}(\mathbf{x}^{II})}{\partial x_i^{II}}, \quad i = 1, \dots, N-1, \qquad (4)$$

while the N-th condition (3) yields

$$g^I(\mathbf{x}^I) - g^{II}(\mathbf{x}^{II}) - \sum_{i=1}^{N-1} \frac{\partial g^I(\mathbf{x}^I)}{\partial x_i^I} (x_i^I - x_i^{II}) = 0. \qquad (5)$$
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^{II}), \quad g(x^I) - g(x^{II}) = g'(x^I)(x^{II} - x^I),$$

as it is shown in Figure 1. Here, $g' = \frac{\partial g}{\partial 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.

<figure><img src="image_4.png" /><figcaption>Figure 1. Two-phase separation in binary mixture: geometrical meaning of conditions (4) and (5).</figcaption></figure>

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

$$\Omega \stackrel{\text{def}}{=} \{ \mathbf{x} = (x_1, x_2) : x_i \in [0, 1], i = 1, 2 \text{ and } x_1 + x_2 \le 1 \}.$$

Consider a smooth surface

$$W \stackrel{\text{def}}{=} \{ (\mathbf{x}, z) \in \mathbb{R}^3 : z = g(\mathbf{x}), \mathbf{x} \in \Omega \}$$

associated with the graph of the function g : Ω → ℝ. The vector field

$$N(\mathbf{x}, z) = g'_{1}(\mathbf{x})\partial_{x_1} + g'_{2}(\mathbf{x})\partial_{x_2} - \partial_{z}, \qquad (6)$$

defines the normal direction to this surface. Here, ∂xᵢ are the coordinate vector fields in ℝ³ and g'ᵢ = ∂g(\mathbf{x}) / ∂xᵢ, i = 1, 2. Let P₁ = (x¹, g(x¹)) and P₂ = (x¹¹, g(x¹¹)) be two points in ℝ³ belonging to the surface W, and such that

$$g'_{1}(\mathbf{x}^{I}) = g'_{1}(\mathbf{x}^{II}), \quad g'_{2}(\mathbf{x}^{I}) = g'_{2}(\mathbf{x}^{II}). \qquad (7)$$

In view of (6), it means that the normals N(P₁) and N(P₂) are collinear.
Further, the vector $\overline{P_1P_2} = (x_1^{II} - x_1^I, x_2^{II} - x_2^I, g(\mathbf{x}^{II}) - g(\mathbf{x}^I))$ belongs to the tangent plane $T_{P_1}W$ attached to the surface W at $P_1$, which means that $\overline{P_1P_2} \perp N(P_1)$, i.e.,

$$g(\mathbf{x}^{II}) - g(\mathbf{x}^I) - g_1'(\mathbf{x}^I)(x_1^{II} - x_1^I) - g_2'(\mathbf{x}^I)(x_2^{II} - x_2^I) = 0. \quad (8)$$

In view of condition (7), the vector $-\overline{P_1P_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 $\Omega$ of the points $P_1$ and $P_2$ along the z-axis define the compositions of splitting liquid phases $\mathbf{x}^I$ and $\mathbf{x}^{II}$.

The projection of the bitangent segment $P_1P_2$ on $\Omega$ 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 $\mathbb{R}^3$, with the bitangent segments $P_1P_2$ being its generators. Projections of these directrices on $\Omega$ 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(\mathbf{x}) \stackrel{\text{def}}{=} \left[ \frac{\partial^2 g(\mathbf{x})}{\partial x_i \partial x_j} \right]_{i,j=1}^{N-1}$$

the Hessian associated to the function $g$. Then, the Gauss curvature of surface W takes the form [17]:

$$K(\mathbf{x}) = \frac{\det H(\mathbf{x})}{(g_1'(\mathbf{x})^2 + g_2'(\mathbf{x})^2 + 1)^2}$$

According to the sign of $K$, the composition domain $\Omega$ can be 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 $\{(\mathbf{x}, z) \in W : K(\mathbf{x}) = 0\}$. On the other hand, $K$ is a product of two principal curvatures of the surface: $K = k_1 \cdot 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 $\Omega$ 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 $\tilde{\Omega} \stackrel{\text{def}}{=} \{\mathbf{x} \in \Omega : K(\mathbf{x}) \ge 0\}$. It follows that the phase diagram can be described as the almost-Riemannian manifold $M = (\tilde{\Omega}, H)$ with a border, equipped with an almost-Riemannian metric associated with the Hessian $H(\mathbf{x})$. The sub-domain of $\tilde{\Omega}$ between the binodal and spinodal curves correspond to the metastable domain, while the remaining part delimited by the binodal curve and the boundary of $\Omega$ 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 $\tilde{\Omega}$. The plait point location can be found by solving equations due to Tompa [19]:

$$g_{222}''''g_{11}''^2 - 3g_{122}''''g_{11}''g_{12}'' + 3g_{112}''''g_{12}''^2 - g_{111}''''g_{12}''g_{22}'' = 0, \quad g_{111}''g_{22}'' - g_{12}''^2 = 0 \quad (9)$$
$$ \text{where } g'''_{ijk} = \frac{\partial^3 g}{\partial x_i \partial x_j \partial x_k} $$

Figure 2 illustrates the described geometrical concepts.

<figure><img src="image_5.png" /><figcaption>**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.</figcaption></figure>

### 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 Ω¹ and Ω¹¹ to each of the phases. Consider now an extended configuration space Σ = Ω¹ × Ω¹¹ of dimension four:

$$ \Sigma \stackrel{\text{def}}{=} \{ \mathbf{q} \in \mathbb{R}^4 : \mathbf{q} = (q_1, q_2) \text{ with } q_1 = \mathbf{x}^I \in \Omega^I, q_2 = \mathbf{x}^{II} \in \Omega^{II} \}. $$

Equations (7) and (8) define three smooth sub-manifolds in Σ associated with zero-level sets of three functions $F_i : \Sigma \to \mathbb{R}^1$ such that

$$ F_1(\mathbf{q}) = g'_1(q_1) - g'_1(q_2), \quad F_2(\mathbf{q}) = g'_2(q_1) - g'_2(q_2), \quad (10) $$

$$ F_3(\mathbf{q}) = g(q_2) - g(q_1) + (\nabla_x g(q_1)|_{q_2} - q_1), $$

where $\nabla_x g(q_1) = (g'_1(q_1), g'_2(q_1))$ and $(|)$ denote, respectively, the standard gradient and scalar product in $\mathbb{R}^2$. The intersection of these 3D sub-manifolds defines a one-dimensional sub-manifold $B \subset \Sigma$ such that

$$ B \stackrel{\text{def}}{=} \{ \mathbf{q} \in \Sigma : F_i(\mathbf{q}) = 0, i = 1, 2, 3 \}. \quad (11) $$

The orthogonal projections $\pi_i : \Sigma \to \Omega$ such that $\pi_I(\mathbf{q}) = q_1$ and $\pi_{II}(\mathbf{q}) = q_2$ define two branches $\pi_I(B)$ and $\pi_{II}(B)$ of the binodal curve. In what follows, the sub-manifold B will be referred to as the generalized binodal curve.
Let $V(\mathbf{q}) \in T_q B$ be the tangent vector at point $\mathbf{q} \in B$ of the generalized binodal $B$ defined by (11). By construction, $V \in \bigcap_{i=1}^3 \mathrm{Ker} \nabla_q F_i$, so that by definition,

$$ \nabla_q F_i(\mathbf{q}) \cdot V(\mathbf{q}) = 0, \quad 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 $\mathbf{q} \in B$:

$$ \begin{pmatrix} 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) \\ \Phi_1(\mathbf{q}) & \Psi_1(\mathbf{q}) & 0 & 0 \end{pmatrix} \cdot \begin{pmatrix} V_1 \\ V_2 \\ V_3 \\ V_4 \end{pmatrix} = 0, \quad (12) $$

where the functions $\Phi_1$ and $\Psi_1$ are defined by the relation

$$ \begin{pmatrix} \Phi_1 \\ \Psi_1 \end{pmatrix} = H(q_1) \cdot (q_1 - q_2). $$

Rewriting (12) in a more compact form yields

$$ \Phi_1 V_1 + \Psi_1 V_2 = 0, \quad H(q_2) \begin{pmatrix} V_3 \\ V_4 \end{pmatrix} = H(q_1) \begin{pmatrix} V_1 \\ V_2 \end{pmatrix}. \quad (13) $$

If $\det H(q_2) \neq 0$ and at least one of the functions $\Phi_1, \Psi_1$ is non-zero, one obtains the solution of (13) in the form

$$ V_1 = \Psi_1 \det H(q_2), \quad V_2 = -\Phi_1 \det H(q_2) $$

and

$$ \begin{pmatrix} V_3 \\ V_4 \end{pmatrix} = (H(q_2))^{-1} \cdot H(q_1) \cdot \begin{pmatrix} V_1 \\ V_2 \end{pmatrix}, $$

Performing all necessary simplifications, we obtain

$$ \begin{pmatrix} V_1 \\ V_2 \\ V_3 \\ V_4 \end{pmatrix} = \begin{pmatrix} \Psi_1 \det H(q_2) \\ -\Phi_1 \det H(q_2) \\ \Psi_2 \det H(q_1) \\ -\Phi_2 \det H(q_1) \end{pmatrix} \quad (14) $$

where, by definition,

$$ \begin{pmatrix} \Phi_2 \\ \Psi_2 \end{pmatrix} = H(q_2) \cdot (q_1 - q_2). \quad (15) $$

The above expressions define a smooth vector field $V \in T\Sigma$. It is well defined except for the singular points $\mathbf{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.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Four-dimensional lifting of the binodal curve. The complete configuration space of the binodal curve is $\Sigma = \Omega^I \times \Omega^{II}$. The vector field $V = (V_1, V_2, V_3, V_4)$ generates a smooth curve $B$ (blue) in $\Sigma$. The orthogonal projections of $B$ onto $\Omega$ 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.</figcaption></figure>

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 $\tilde{\Omega}$ associated with the Hessian matrix $H(\mathbf{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 $\dot{\mathbf{q}} = V(\mathbf{q})$ in $\Sigma$. 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 $\tilde{\Omega}$, 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 $\Sigma$, 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 $\Omega$, 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 $\Omega$.

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 $\Omega$. Indeed, being the solution of the equation $H(\mathbf{x}) = 0$, the spinodal curve is a solution to the differential equation associated with the vector field $S(\mathbf{x}) = \nabla_x H(\mathbf{x})^{\perp}$ in $\Omega$, which, in general, is regular at plait points. Here, the superscript $\perp$ denotes the orthogonal complement to the vector $\nabla_x H(\mathbf{x})$ in $\mathbb{R}^2$. In other words, the spinodal curve is the solution of the ODE equation $\dot{\mathbf{x}} = S(\mathbf{x})$, $\mathbf{x} \in \Omega$; the starting point for the integration can be detected by finding one inflexion point of the function $g$ on the boundary of $\Omega$ or inside $\Omega$ 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 UNIQUAC 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 = \frac{x_1}{N1} \ln x_1 + \frac{x_2}{N2} \ln x_2 + \frac{x_3}{N3} \ln x_3 + \chi_{12}x_1x_2 + \chi_{13}x_1x_3 + \chi_{23}x_2x_3 + \beta x_1x_2x_3. \quad (16)$$
Here $x_i$, $i = 1, 2, 3$ are the volume fractions of the components, and $N_2$ and $N_3$ are 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 $\chi_{ij}$ and $\beta$ 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 $\beta = 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 $\chi_{ij}$ 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$, $\chi_{ij}$, and $\beta$ 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} \left( \sum_{k=1}^{n_b} \frac{F_1^2(q_k) + F_2^2(q_k) + F_3^2(q_k)}{(1 - F_3^2(q_k))^2} \right), \qquad (17) $$

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)

$$ \sigma = 100 \left( \frac{1}{6n} \sum_{k=1}^{n_b} \sum_{i}^{3} \left( x_{i,k}^{I,FH} - x_{i,k}^{I,exp} \right)^2 + \left( x_{i,k}^{II,FH} - x_{i,k}^{II,exp} \right)^2 \right)^{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><img src="image_7.png" /><figcaption>**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.</figcaption></figure>

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

<table><caption>**Table 1.** Root mean square deviation $\sigma$ of analysed mixture.</caption><thead><tr><th>System</th><th>Model</th><th>T, C</th><th>σ</th></tr></thead><tbody><tr><td>water–DMSO–THF</td><td>FH (this work)</td><td>20°</td><td>3.53</td></tr><tr><td>water–DMSO–THF, [25]</td><td>NRTL</td><td>20°</td><td>5.83</td></tr><tr><td>water–DMSO–THF, [8]</td><td>UNIQUAC</td><td>20°</td><td>3.4</td></tr><tr><td>water–DMSO–THF, [24]</td><td>NRTL</td><td>20°</td><td>3.18, 3.09</td></tr><tr><td>water–DMSO–THF, [24]</td><td>UNIQUAC</td><td>20°</td><td>3.97, 3.4</td></tr><tr><td>water–acetone–phenol</td><td>FH (this work)</td><td>50°</td><td>0.88</td></tr><tr><td>water–acetone–phenol, [26]</td><td>NRTL</td><td>50°</td><td>0.81</td></tr><tr><td>water–acetone–phenol, [24]</td><td>NRTL</td><td>50°</td><td>1.13</td></tr><tr><td>water–acetone–phenol, [24]</td><td>UNIQUAC</td><td>50°</td><td>1.17</td></tr><tr><td>water–acetone–phenol</td><td>FH (this work)</td><td>56°</td><td>1.27</td></tr><tr><td>water–acetone–phenol, [23]</td><td>NRTL</td><td>56°</td><td>1.61</td></tr><tr><td>water–acetone–phenol, [23]</td><td>UNIQUAC</td><td>56°</td><td>1.52</td></tr><tr><td>water–acetone–hexadecane</td><td>FH (this work)</td><td>20°</td><td>2.92</td></tr></tbody></table>

### 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 $\chi_{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><img src="image_8.png" /><figcaption>**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.</figcaption></figure>

#### 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, Tables A1 and 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$, $\chi_{13}$, and $\chi_{23}$. Due to the significant variation of the solubility data of hexadecane in water in different sources, the term $\chi_{13}$ was replaced by the $\chi_{13} + \delta$, where the correction term $\delta$ 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.

<figure><img src="image_9.png" /><figcaption>**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.</figcaption></figure>
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_{tot}$ denote the molar free Gibbs energy. Since the mole fractions $x_i = \frac{n_i}{n_{tot}}$ verify $\sum_{i=1}^{N} x_i = 1$, only $N-1$ of them can be used as independent variables. In what follows, we denote $\mathbf{x} = (x_1, ..., x_{N-1})$, so that $n_N = n_{tot}(1 - x_1 - ... - x_{N-1})$. Since $dx_N = - \sum_{i=1}^{N-1} dx_i$, the Gibbs Equation (2) takes the form

$$dg = -sdT + v dP + \sum_{i=1}^{N-1} (\mu_i - \mu_N)dx_i, \quad (A1)$$

where $s = S/n_{tot}$ and $v = V/n_{tot}$ are the specific molar entropy and volume, correspondingly. Moreover, for all $i = 1, ..., N-1$ we obtain

$$\frac{\partial g}{\partial x_i} = \frac{1}{n} \frac{\partial G}{\partial x_i} = \frac{1}{n} \left( \frac{\partial G}{\partial n_i} \frac{\partial n_i}{\partial x_i} + \frac{\partial G}{\partial n_N} \frac{\partial n_N}{\partial x_i} \right) = \mu_i - \mu_N, \quad (A2)$$
and hence,

$$g = \frac{1}{n} \sum_{i=1}^{N} \mu_i n_i = \sum_{i=1}^{N-1} (\mu_i - \mu_N)x_i + \mu_N. \qquad (A3)$$

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:

$$\frac{\partial g^I(\mathbf{x}^I)}{\partial x_i^I} = \frac{\partial g^{II}(\mathbf{x}^{II})}{\partial x_i^{II}}, \quad i = 1, \dots, N-1, \qquad (A4)$$

while the N-th condition (3) yields

$$g^I(\mathbf{x}^I) - g^{II}(\mathbf{x}^{II}) - \sum_{i=1}^{N-1} \frac{\partial g^I(\mathbf{x}^I)}{\partial x_i^I} (x_i^I - x_i^{II}) = 0. \qquad (A5)$$

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

$$\begin{aligned} F_1 = & \chi_{12}(x_2^I - x_2^{II}) - \chi_{13}(2(x_1^I - x_1^{II}) + x_2^I - x_2^{II}) - \chi_{23}(x_2^I - x_2^{II}) \\ & + \beta(x_2^I - 2x_1^I x_2^I - (x_2^I)^2 - x_2^{II} + 2x_1^{II} x_2^{II} + (x_2^{II})^2) \\ & + \frac{1}{N_1} \ln\left(\frac{x_1^I}{x_1^{II}}\right) - \frac{1}{N_3} \ln\left(\frac{1 - x_1^I - x_2^I}{1 - x_1^I - x_2^{II}}\right) \\ F_2 = & \chi_{12}(x_1^I - x_1^{II}) - \chi_{13}(x_1^I - x_1^{II}) - \chi_{23}(x_1^I - x_1^{II} + 2x_1^I - 2x_2^{II}) \\ & + \beta(x_1^I - (x_1^I)^2 - x_1^{II} + (x_1^{II})^2 - 2x_1^I x_2^I + 2x_1^{II} x_2^{II}) \\ & + \frac{1}{N_2} \ln\left(\frac{x_2^I}{x_2^{II}}\right) - \frac{1}{N_3} \ln\left(\frac{1 - x_1^I - x_2^I}{1 - x_1^{II} - x_2^{II}}\right) \end{aligned} \qquad (A6)$$

$$\begin{aligned} F_3 = & \chi_{12}(x_1^I - x_1^{II})(x_2^I - x_2^{II}) - \chi_{13}(x_1^I - x_1^{II})(x_1^I - x_1^{II} + x_2^I - x_2^{II}) \\ & - \chi_{23}(x_2^I - x_2^{II})(x_1^I - x_1^{II} + x_2^I - x_2^{II}) \\ & + \beta(x_1^I x_2^I - 2(x_1^I)^2 x_2^I - x_1^{II} x_2^I + 2x_1^I x_1^{II} x_2^I - 2x_1^I (x_2^I)^2 + x_1^{II} (x_2^I)^2 \\ & - x_1^I x_2^{II} + (x_1^I)^2 x_2^{II} + x_1^{II} x_2^{II} - (x_1^{II})^2 x_2^{II} + 2x_1^I x_2^I x_2^{II} - x_1^{II} (x_2^{II})^2) \\ & + \frac{1}{N_1} \left(x_1^I - x_1^{II} - x_1^{II} \ln\left(\frac{x_1^I}{x_1^{II}}\right)\right) + \frac{1}{N_2} \left(x_2^I - x_2^{II} - x_2^{II} \ln\left(\frac{x_1^I}{x_1^{II}}\right)\right) \\ & - \frac{1}{N_3} \left(x_1^I - x_1^{II} + x_2^I - x_2^{II} - (1 - x_1^{II} - x_2^{II}) \ln\left(\frac{1 - x_1^I - x_2^I}{1 - x_1^{II} - x_2^{II}}\right)\right) \end{aligned}$$

The phase coexistence conditions (4) and (5) take the form $F_i(\mathbf{q}) = 0$ for $i = 1, 2, 3$ and $\mathbf{q} = (x_1^I, x_2^I, x_1^{II}, x_2^{II}) \in \Sigma$. Remarkably, for the FH model, these conditions are linear with respect to the interaction parameters $\chi_{ij}$, $\beta$ 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^{II}$ are available through the experimental procedure. The binary mixture of components 1–3 corresponds to the side $x_2 = 0$ of the boundary of $\Omega$. Setting $x_1 = x$, $x_3 = 1 - x$, the function $g$ can be rewritten in the form

$$g_1 = \frac{x}{N_1} \ln x + \frac{1-x}{N_3} \ln(1-x) + \chi_{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^{II})\chi_{13} + \frac{1}{N_3}\ln\left(\frac{1-x^I}{1-x^{II}}\right) - \ln\left(\frac{x^I}{x^{II}}\right) = 0, \qquad (A7)$$

$$x^I - x^{II} - \chi_{13}(x^I - x^{II})^2 + (1 - x^{II}) \ln\left(\frac{1 - x^I}{1 - x^{II}}\right) = 0. \qquad (A8)$$

If the values of $x^I$ and $x^{II}$ are known, the above equations form a system of two linear algebraic equations with respect to $\chi_{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$, $\chi_{13}$, and $\chi_{23}$ by solving a system of four linear algebraic equations in terms of the miscibility gap limits $x_1^I$, $x_1^{II}$ and $x_2^I$, $x_2^{II}$.

## Appendix C. Water–Acetone–Hexadecane Experimental Data

<table><caption>Table A1. Water–acetone–hexadecane tie-lines, mass fraction experimental data.</caption><thead><tr><th colspan="3">Phase I</th><th colspan="3">Phase II</th></tr><tr><th>Water</th><th>Acetone</th><th>Hexadecane</th><th>Water</th><th>Acetone</th><th>Hexadecane</th></tr></thead><tbody><tr><td>0.005026</td><td>0.088051</td><td>0.906923</td><td>0.1153287</td><td>0.88</td><td>0.0046713</td></tr><tr><td>0.002956</td><td>0.121763</td><td>0.875281</td><td>0.07389669</td><td>0.9115</td><td>0.014603313</td></tr><tr><td>0.004319</td><td>0.127764</td><td>0.867917</td><td>0.05061999</td><td>0.919</td><td>0.030380015</td></tr><tr><td>0.000775</td><td>0.153273</td><td>0.845952</td><td>0.02604692</td><td>0.9075</td><td>0.06645308</td></tr><tr><td>0.000313</td><td>0.284572</td><td>0.715115</td><td>0.01003813</td><td>0.865</td><td>0.124961871</td></tr><tr><td>0.002187</td><td>0.023547</td><td>0.974265</td><td>0.93897</td><td>0.06102</td><td>0.000000001</td></tr><tr><td>0.002148</td><td>0.029435</td><td>0.968417</td><td>0.87777</td><td>0.12223</td><td>0.000000001</td></tr><tr><td>0.001949</td><td>0.036801</td><td>0.96125</td><td>0.793004</td><td>0.20699</td><td>0.000000001</td></tr><tr><td>0.000308</td><td>0.048658</td><td>0.951034</td><td>0.67413</td><td>0.3259</td><td>0.000000001</td></tr><tr><td>0.003093</td><td>0.054403</td><td>0.942503</td><td>0.56358</td><td>0.4364</td><td>0.000000001</td></tr><tr><td>0.000276</td><td>0.065321</td><td>0.934403</td><td>0.49548</td><td>0.50452</td><td>0.000000001</td></tr><tr><td>0.005298</td><td>0.074969</td><td>0.919733</td><td>0.3655</td><td>0.63446</td><td>0.000000001</td></tr><tr><td>0.003066</td><td>0.082343</td><td>0.914591</td><td>0.2541</td><td>0.74585</td><td>0.000000001</td></tr><tr><td>0.002307</td><td>0.097121</td><td>0.900573</td><td>0.1747</td><td>0.82524</td><td>0.000000001</td></tr></tbody></table>
<table><caption>**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.</caption><thead><tr><th>Water</th><th>Acetone</th><th>Hexadecane</th></tr></thead><tbody><tr><td>0.3503</td><td>0.6496</td><td>8.98 × 10<sup>-5</sup></td></tr><tr><td>0.1989</td><td>0.7996</td><td>0.0015</td></tr><tr><td>0.1519</td><td>0.846</td><td>0.002</td></tr><tr><td>0.115</td><td>0.882</td><td>0.0033</td></tr><tr><td>0.10056065</td><td>0.89269389</td><td>0.00674546</td></tr><tr><td>0.06122825</td><td>0.91649876</td><td>0.02227299</td></tr><tr><td>0.04112357</td><td>0.9178071</td><td>0.04106932</td></tr><tr><td>0.02740997</td><td>0.90841079</td><td>0.06417925</td></tr><tr><td>0.01445356</td><td>0.88430202</td><td>0.10124442</td></tr><tr><td>0.0039</td><td>0.8456</td><td>0.1504</td></tr><tr><td>0.002</td><td>0.799</td><td>0.199</td></tr><tr><td>0</td><td>0.783</td><td>0.217</td></tr><tr><td>0</td><td>0.264</td><td>0.736</td></tr><tr><td>0.0005</td><td>0.1975</td><td>0.8021</td></tr><tr><td>0.0007</td><td>0.1589</td><td>0.8404</td></tr><tr><td>0.00092</td><td>0.1122</td><td>0.8869</td></tr><tr><td>0.00065</td><td>0.0603</td><td>0.9391</td></tr></tbody></table>

## References

1. Koningsveld, R.; Stockmayer, W.H.; Nies, E. *Polymer Phase Diagrams*; Oxford Univ. Press: Oxford, UK, 2001.

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. [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. [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 Equilibria* **2013**, *353*, 121–131. [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. [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.

7. Michelsen, M.L. The isothermal flash problem. Part I. Stability. *Fluid Phase Equilibria* **1982**, *9*, 1–19

8. Binous, H.; Bellagi, A. Calculation of ternary liquid-liquid equilibrium data using arc-length continuation. *Wiley Eng. Rep.* **2021**, *3*, e12296.

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. [CrossRef]

10. Levelt, S. *How Fluids Unmix*; Royal Netherlands Academy of Arts and Sciences: Amsterdam, The Netherlands, 2002

11. Arnold, V.I. *Singulatrities of Caustics and Wave Fronts*; Kluwer: Baltimore, MD, USA, 1990.

12. Varchenko, A.N. Evolutions of convex hulls and phase transitions in thermodynamics. *J. Math. Sci.* **1990**, *52*, 3305–3325. [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. [CrossRef]

14. Gaite, J.A. Phase transitions as catastrophes: The tricritical point. *Phys. Rev. A* **1990**, *41*, 5320–5324. [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. [CrossRef]

16. Treybal, R.E. *Liquid Extraction*, 2nd ed.; McGraw-Hill: New-York, NY, USA, 1963.

17. Do Carmo, M. *Differential Geometry of Curves and Surfaces*, 2nd ed.; Courier Dover Publications: Mineola, NY, USA, 2017.

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. [CrossRef]

19. Tompa, H. *Polymer Solutions*; Butterworths: London, UK, 1956.
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.

21. Gmehling, J.; Kolbe, B.; Kleiber, M.; Rarey, J. *Chemical Thermodynamics for Process Simulation*, 2nd ed.; Willey-VCH Verlag: Weinheim, Germany, 2019.

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.

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. [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. [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. [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. [CrossRef] [PubMed]

**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.