## Abstract

Global phase diagrams (GPD) of binary mixtures in phase equilibria modeling are analyzed. The mapping of the global equilibrium surface in the parameter space of the equation of state (EoS) model provides the most comprehensive system of criteria for predicting binary mixture phase behavior. One may obtain the relationships for azeotropic boundaries from the global phase diagram [A (azeotrope) and H (hetero-azeotrope)] regions. Analytical expressions to predict azeotrope and double azeotrope phenomena in terms of critical parameters of pure components were derived using global phase diagram. The problem estimations of phase behavior modeling under the uncertainty are formulated applying the Pareto-optimum parameter and different (crisp and fuzzy) convolution schemes. The Pareto-optimum parameters in the Redlich-Kwong equation of state used different conflicting data sets (simultaneous description of the phase equilibria and critical line data in binary mixtures, thermodynamically consistent description of the inhomogeneous data). Ionic liquids (ILs) are one of prospective new working media for different environmentally friendly technologies. Practically undetectable vapor pressure is considered the ILs as ideal solvents replacing conventional solvents in the frame of a “green chemistry.” Combination of ionic liquids with conventional natural and synthetic refrigerants promotes the increasing efficiency of absorption processes due to nonvolatile ionic liquids (absorbents).

### Keywords

- azeotrope
- phase equilibria
- global phase diagram
- equation of state
- the Pareto-optimum parameters
- azeotrope-breaking
- ionic liquids
- refrigerant blends

## 1. Introduction

The pioneering work of Van Konynenburg and Scott [1] demonstrated that the van der Waals one-fluid model has wide possibilities of qualitative reproducing the main types of phase diagrams of binary mixtures. The thermodynamic equilibrium mapping onto the space of parameters of an equation of state gives the possibility to obtain the criteria for the phase behavior of binary mixtures. The functions of temperature *T*, pressure *p*, chemical potential *μ*, and component concentration *x* are determined as the “field” variables that are the same for all phases coexisting in equilibrium of a substance *i*. The molar fraction is the “density” that is in principle of liquid *(l)* and gas *(g)* phases. Global phase diagrams (GPD) of binary mixtures represent boundaries between types of phase behavior in a dimensionless space of EOS parameters.

GPD also provides good visualization of the impact of model parameters of mixture components to the topology of phase behavior. The proposed types of phase behavior classification [1] are generally used to characterize the different types of phase behavior in binary mixtures.

Varchenko [2] has provided a more rigorous classification for conventional features of equilibrium surface and phase diagrams for binary mixtures with strict determination of the eight topologically different rearrangements.

Current topological analysis of equilibrium surfaces of binary fluid systems provides 26 singularities and 56 scenarios of phase behavior evolution depicted in *p*–*T* diagrams [3].

The various phase diagram classes and *p–T* projections of the main types of phase diagrams have been described in literature [1, 2, 3]. Global phase diagrams are a technique which can be used for the prediction of different phase behavior in the mixtures without vapor–liquid equilibrium calculations [4, 5, 6, 7, 8, 9, 10, 11].

Patel and Sunol [12] developed a robust automated routine for global phase diagram generation in binary systems. The approach uses any equation of state models, takes into account solid-phase existence, and provides type VI phase diagram generation. The generated data set includes calculations of the critical endpoints, quadruple points, critical azeotrope points, azeotrope endpoints, pure azeotrope points, critical line, liquid–liquid–vapor line, and azeotrope line.

Azeotrope-breaking is important for the successful distilling of industrially important mixtures. To simulate the mixture phase behavior, models based on the equation of state (EoS) presentation for thermodynamic properties are more preferable. The conventional methods of parameter identification use the statistical paradigm, which is based on maximum likelihood or a posteriori probability criteria and does not take into account uncertainties of vague nature. Decision-making process under various uncertainties requires mathematical methods, which include uncertainty evaluation a priori. Statistical methods interpret all variety of uncertainty types in the framework of the randomness concept. Nevertheless, there are ill-structured situations, which have not any strictly defined boundaries and cannot be accurately formulated.

The main challenge is to deal with such kind of expressions like “neighborhood” or “best fit,” which do not have strict boundaries, separating one class of objects from others. Generally, there are ambiguous verbal models, which can be treated as fuzzy formulated targets, depending on biased assessment of boundaries for approximations used. As a case study, we provide estimation of the optimal parameters of the Redlich-Kwong equation of state [13], retrieved from the different conflicting data sets resulting from the inconsistency problems arising in the modeling of phase equilibria. Such process reflects different types of uncertainties, including uncertainties of nonstatistical origin. The parameters of phase equilibria models are considered as alternatives, i.e., they allow meeting the targets and prescribed constraints. The parameter estimation problem of phase equilibrium modeling multicriteria approach is applied. To describe the uncertainties of such type, the “fuzzy set” approach introduced by L.A. Zadeh [14] is used.

The problem of optimum parameter estimation of thermodynamic and phase behavior under the uncertainty is a search of the Pareto set. The diverse computational methods of the Pareto-optimum parameter convolution crisp and fuzzy schemes to reduce a vector criterion into the scalar are provided.

Ionic liquids can be treated as “adjustable” working fluids given the fact that variations of different “R – ” groups and cation/anion ratio selection ensure meeting preferred trade-off solution between density, viscosity, melting point, and other physicochemical properties. A great role in the IL applications plays the azeotropic phenomena.

Considering a huge number of the “off the shelf” mixtures scheduled for destroying or recycling, it becomes key importance to develop the theoretically sound methods for reliable assessment of thermodynamic and phase behavior. Experimental retrieval of the azeotrope properties is time-consuming and a costly process. Availability of theoretical predictions for azeotropic phenomena would not only reduce this cost but also make efficient required experimental investigations.

The objective of this study is to present novel approach which encompasses global phase diagram technique and Pareto-optimal EoS parameter estimation for azeotropic and double azeotropic criteria evaluation in terms of critical parameters of pure components for the binary mixtures of refrigerants and ionic liquids.

This paper is organized as follows. In Section 2, the phase diagram classification and the one-fluid model of the equation of state for global phase diagrams are presented. Section 3 provides approach to sustainable property calculation based on a fuzzy set methodology providing a trade-off solution among different sources of data required for regression of model parameters. Section 4 discusses the results of azeotrope-breaking criteria for IL-refrigerant mixtures based at the approaches provided in previous sections.

## 2. Phase diagram classes and global phase diagram

Mapping of the global equilibrium surface into the parameter space of the equation of state model provides the most comprehensive set of criteria for prediction phase behavior of the binary mixture.

The impact of critical parameters of components on phase topology is provided via global phase diagrams. The diagrams are depicted in the space of the equation of state parameters, e.g., van der Waals *a* and *b* parameters. The specific points and lines of global phase diagrams, including tricritical points, double critical endpoints, azeotropic lines, etc., create the boundaries at the diagram and divide the model parameter space into the regions corresponding to the various types of the phase behavior. Generally, global phase diagram is expressed in dimensionless parameters which depend on the equation of state model. The global phase diagrams of such different models as the one-fluid EoS of binary Lennard-Jones fluid [6, 7] and the Redlich-Kwong model [5, 8, 9] are almost identical, in particular for the case of equal-sized molecules.

We consider here the cubic-type models:

where *R* is the universal gas constant and the EoS parameters *a* and *b* of mixture depend on the mole fractions *xi* and *xj* of the components *i* and *j* and on the corresponding parameters *aij* and *bij* for different pairs of interacting molecules:

The convenient set of dimensionless parameters for the Redlich-Kwong model is as follows [13]:

where.

where *Z1–Z2* plane at fixed values of Z_{3} and Z_{4} is a straight line. The values are equal to Ω_{a} = 0.42747, Ω_{b} = 0.08664, and Z_{c} = 0.333.

The combining rules for the binary interaction parameters are

where *kij* and *lij* are fitting coefficients in the Lorentz**-**Berthelot combining rule (*kij = lij* = 0).

The simplest boundary is a normal critical point when two fluid phases are becoming identical. Critical conditions are expressed in terms of the molar Gibbs energy derivatives in the following way:

Corresponding critical conditions for the composition-temperature-volume variables are

where A is the molar Helmholtz energy and.

*VC* and *TC* at a given concentration *x*.

Global phase diagrams generated for Redlich-Kwong model is presented in Figures 1 and 2.

Here C_{1} and C_{2} are critical points of components 1 and 2; C_{m} is hypothetic critical point beyond solidification line. The types of phase behavior indicated as Roman numbers are described in [1, 2, 3].

## 3. Uncertainties and conflicts in the parameter estimation

Accounting the effects of uncertainty in the fluid-phase equilibria, modeling became important in the recent decade. Generally, the uncertainty tried to be resolved via probability and random process theories. However, the probabilistic methods can lead to the unreliable estimation of parameters. There are three main complimentary models of uncertainty described in literature:

Convex models of uncertainty, developed by Ben-Haim and Elishakoff [15] as an extension of interval analysis

Game-theoretic models of uncertainty deriving from

*conflict*among the different goals (in our case, it is a conflict between thermodynamically inconsistent data)Verbal models of uncertainty deriving from

*vagueness*and initiated to be resolved by Zadeh [14]

Three models of uncertainty cannot exist one by one in parameter estimation, and mathematical tools should reveal this fact. The conventional single-criterion methods of parameter estimation are examples of lopsided vision of multicriteria decision-making problem where only one set of variables, strongly dependent on decision-maker experience, is recommended. A lack of single meaning of optimality concept, following from the basics of modern game theory, points to impossibility to construct a completely formal algorithm for parameter estimation and futility of eliminating the uncertainty in the search process. It is more correct to consider the trade-off or rational parameters, which adequately describe thermodynamic properties and phase behavior, generated by EoS.

The following sequence of decision-making steps in model parameter estimation is proposed:

Determination of the Pareto-optimal (or compromise or trade-off) set

**X**_{P}as the formal solution of the multicriteria problem to minimize a conflict source of uncertaintyInformal selection of convolution scheme to switch over a vector criterion

**K**into a scalar combination of the*Ki (Х)*Evaluation of the final decision vector

**X**_{opt}**∈ X**_{P}to minimize a vagueness source of uncertainty

### 3.1 Pareto set

The most important stage toward the multicriteria problem solving is an establishment of the Pareto domain *ХР*, i.e., such domain in model parameter space where it is not possible to reach a dominance of selected criteria over others. A geometrical visualization of the Pareto set (**AB** line) for the bi-criteria case in *2D* parameter space is given in Figures 3 and 4. For instance, the best least squares fit of the *p–T-x* phase equilibria data *X1* and *X2* could be interpreted as geometric and energetic parameters in the van der Waals EoS model.

To locate the Pareto set, we assume to compare two solutions, **X**^{*} and **X**_{0}, which hold the inequality:

It is obvious that **X**_{0} is preferable than **X**^{*}. Thus, all **X**^{*} vectors satisfying Eq. (8) may be excluded. It is worth to analyze that only those **X**^{*} vectors for which there is no such **X**_{0} when for all criteria the inequality (8) is satisfied. The set of all values **X**_{P} **= X**^{*} is the Pareto set, and the vector **X**_{P} is the unimprovable vector of the results, in case if from *Ki* **(X**_{0}**)** ≤ *Ki* **(X**_{P}**)** for any **i** it follows that **K (X**_{0}**) = K (X**_{P}**)**.

### 3.2 Crisp convolution schemes

Isolationistic and cooperative trends are usually considered for aggregation of local criteria vector into global scalar criterion. The isolationistic convolution schemes are additive, i.e., global criterion is represented as a weighted sum of local criteria and entropy as a sum of local criterion logarithms. If behavior of each criterion is complied with common decision to minimize some cooperative criterion, then a convolution scheme can be presented in the form

where *wi* are the weight coefficients, *KC* is a global trade-off criterion. If it will be possible to come to an agreement about preference (weight) for each criterion, then the final decision can be found as a solution of scalar nonlinear programming problem:

In terms of game-theoretical approach, such coalition is defined as the von Neumann coalition scheme [16].

If no concordance among decision-makers is concerning the weight choice, then arbitration network is preferable. Classical arbitration scheme was derived mathematically rigorously by Nash but very often criticized from common sense [16]:

All crisp convolution schemes under discussion try to reduce an uncertainty deriving from conflict among different criteria in the Pareto domain. The next step is extenuation of uncertainty driving from vagueness.

### 3.3 Fuzzy convolution scheme

Zadeh [14] put the theory of fuzzy sets forward with explicit reference to the vagueness of natural language, when describing quantitative or qualitative goals of the system. Here we assume that local criteria as well as different constraints in the ill-structured situation can be represented by fuzzy sets. A final decision is defined by the Bellman and Zadeh model [17] as the intersection of all fuzzy criteria and constraints and is represented by its membership function **μ(Х)** as follows:

This problem is reduced to the standard nonlinear programming problems: to find the values **Х** and *λ* that maximize *λ* that is subject to

Here, by way of illustration, the new approach was introduced to estimate the Redlich-Kwong EoS parameters for simultaneous description of the phase equilibria and critical line data in binary mixtures, thermodynamically consistent description of the inhomogeneous data, and other inconsistency problems arising in the modeling of phase equilibria.

### 3.4 Conflict between phase equilibria and critical line description in binary mixtures

Reliable models for thermodynamic and phase behavior description of binary mixtures are facing problems of adequate estimation of the binary interaction parameters from the restricted set of *VLE* data. In spite of availability of strict relationships to obtain from EoS model, the different derivatives of thermodynamic values, for instance, the prediction of critical lines from the EoS with parameters restored from *VLE* data, is questionable due to diverse sources of uncertainty in both the used models and experimental data.

For illustration, we consider thermodynamic and phase behavior of water-carbon dioxide system near the critical point of water. Experimental data on *p-ρ-T-x* properties have been taken from [18]; data on phase equilibria and critical curve have been extracted from [19, 20]. The critical lines have been calculated using algorithm developed in [21]. Phase equilibria calculations have been carried out with Michelsen and Mollerup method [22] for cubic EoS. We performed phase equilibria and critical line calculations for the RK EoS with binary interaction parameters *k12* and *l12* used in standard mixing rules

and restored from different crisp convolution schemes. For simplicity, the results of calculation for phase equilibria and critical line are presented in Figures 5 and 6 only for most selected compromise schemes. The challenge of conflict between parameters retrieved from different sections of thermodynamic surface remains important for practically all EoS including sophisticated multiparameter equation of state. In general, the conflict can be described as follows: if interaction parameters are retrieved from the one class of properties, the prediction of other properties is doubtful in spite of validity of thermodynamic relationships. The final solution will depend on the problem setting and decision-maker subjective experience.

## 4. Results and discussion

This study provides a novel approach for defining the quantifiable estimation of boundary states and specific points which define changes of phase behavior. The global phase diagram technique is applied for deriving the analytical expressions and Pareto-based approach with fuzzy convolution scheme used for adequate evaluation of model parameters.

There are no rigorous mathematical methods to construct adequate model without subjectively based solution due to a number of uncertainties and conflicts. To reduce the source of unavoidable uncertainty, the different compromise schemes of criteria convolution are considered. It is important to note that the selection of the convolution scheme is a subjective choice of the decision-maker.

Below we provide results and discussion of application of the abovementioned approach for identification of the azeotrope- and double azeotrope-breaking criteria for imidazolium (IL)-based ionic liquids and industrial refrigerant mixtures. IL doping leads to the breaking of azeotrope in binary refrigerant mixtures that open the way for the azeotrope refrigerant mixture separation technologies in order to remove the environmentally harmful substances.

### 4.1 Azeotrope-breaking criteria

Applying the cubic model of the equation of state, only the critical properties and acentric factor of the individual components in mixtures are sufficient to define the phase behavior in the interested sector of the parameters.

The degenerated critical azeotrope point is a boundary state separating azeotrope and non-azeotrope and produces the limit of the critical azeotrope at *xi* → *0* or at *xi* → *1*. The solution of the thermodynamic equation system for a degenerated critical azeotrope [5, 6, 7, 9]

gives the following relationship for dimensionless parameters *Zi*:

where the upper signs + or − correspond to *x2 = 0* and the lower signs - to *x2 = 1*, respectively.

To define the azeotropic states, the parameters *Zi* can be expressed via critical (pseudocritical) temperatures *Tci* and pressures *Pci* of pure components and empirical binary interaction parameters *k12* and *l12* (17). The unlike pair interaction parameters *Z2* and *Z4* (i.e., *a12* and *b12*) can be estimated solving simultaneously the system of Eqs. (17) for the given *Z1* and *Z3* (or the set of pure component constants *a11*, *b11*, *a22*, *b22*) that are determined from the critical parameters of the components:

The boundary between azeotrope and zeotrope states in *Z1–Z2* plane at fixed values of Z_{3} and Z_{4} is a straight line. The definition of azeotropic and zeotropic states for refrigerants is depicted in Figure 7.

### 4.2 Double azeotrope-breaking criterion

Double azeotropic phenomena can be observed in the global phase diagram only in the vicinity of the crossing point of two straight azeotropic borders in (*Z1*, *Z2*) plane (Figures 8 and 9):

The best results were demonstrated by application of the Schwarzentruber-Renon and Wang-Sandler mixing rule [10] with three adjustable parameters (Figure 8). The use of local mapping shows good results of double azeotropic description at 323.05 К using the standard mixing rule (Figure 8).

### 4.3 Applications of ionic liquids in azeotrope mixture distillation

The set of parameters for a given equation of state model univocally defines a global phase diagram and, consequently, evolution of phase behavior for binary mixture in a wide range of temperatures and pressures. Global phase diagram for model systems of components explores that binary mixture of industrial refrigerants with imidazolium (IL)-based ionic liquids does not experience azeotropic behavior in the practicable range of parameters.

Distribution of critical points for major components of refrigerants [25, 26, 27, 28, 29] and hypothetical critical parameters of a number of imidazole-based ionic liquids is depicted in Figure 10.

It shall be noted that at subcritical temperatures, IL may undergo thermal decomposition that causes uncertainty in assessment of critical parameters. Usually, the available values of critical points of ionic liquids are derived from low-temperature regions having extremely small saturation pressures. Computational procedures for finding model parameters are unstable. This can lead to errors in the predictions of phase equilibria at high temperatures.

Thanks to undetectable vapor pressure, ILs are considered as potential environmentally friendly candidates for replacing conventional solvents. The selective solubility properties of the ILs appeared for particular components of the mixtures and can be treated as extraction media for separation processes. Moreover, the increase of efficiency for absorption processes is promoted thanks to nonvolatile feature of ionic liquids acting as absorbents.

The information on breaking of azeotrope is crucial for the separation of different industrial mixtures and is intensively discussed in the literature. The proposed existing expressions for identification of azeotrope behavior generally are empirical.

To describe and predict the phase behavior of the mixture, the equation of state models is more suitable for calculation of the thermodynamic properties.

The calculations of phase equilibria for system imidazolium-based ionic liquid and refrigerants R134a and R1234yf are performed. The Redlich-Kwong one-fluid model equation of state was selected due to simplicity, few parameters to be estimated, existing robust computational algorithms for obtaining the derivations of different thermodynamic features, as well as a large amount of existing data on binary interacting parameters in the existing literature.

The phase behavior of imidazolium-based ionic liquids С_{8}H_{11}N_{3}F_{6}S_{2}O_{4} ([EMIm][Tf2N]) and C_{10}H_{19}N_{2}BF_{4} ([HMIm][BF_{4}]) with refrigerants R134а and R1234yf was evaluated using the parameters regressed at the low-pressure experimental data. It is noted that the variation of anionic group leads to the shift of critical point of ILs and obviously impact intermolecular interactions between ionic liquids and molecules of refrigerants. To this end the phase behavior pattern is also impacted. Variation in the *k12* interaction coefficient shifts the position of a specific point on the global phase diagram. For R1234yf-[HMIm][BF_{4}] system, the position of the specific point at different values *k12* = − 0.1, 0, + 0.1 demonstrates a tendency to transition from azeotrope to zeotrope state or vice versa. The binary interaction parameters *k12* and *l12* for R134a-Il blends were restored from experimental data provided by Ren et al. [28, 29] with Pareto-based method described in Section 3.

In case a specific point is located in the northern or southern quadrants of the diagrams depicted in Figure 11, the azeotropic phenomenon is expected to appear in the binary mixture. The pattern of specific point location for ionic liquid [HMIm][BF_{4}] + refrigerants R134a (R1234yf) systems is also provided in Figure 11.

It is considered that the azeotrope definitely appears in the R134a-R1234yf system. The literature search provides lack of experimental data for this system. The boundaries presented in Figure 11 for the R1234yf-[HMIm][BF_{4}] mixture practically coincide. The addition of the ionic liquid to azeotropic mixture leads to azeotrope-breaking that is demonstrated by a change of phase envelope for the R134a-R1234yf-[HMIm][BF_{4}] mixture in Figure 12.

Positive value of binary interaction coefficient *k12* can report on III-type of phase behavior for the system studied. The calculations were performed in MATLAB software on the base of the algorithms proposed by Michelsen-Mollerup [23].

## 5. Conclusions

The responses to the environmental challenges require development and application of new environmentally friendly working media. As one of the promising media for chemical and refrigeration industry, the mixtures of conventional refrigerants with ionic liquids are considered. This requires reliable data on thermophysical properties and phase behavior of the mixture.

Tremendously growing amount of data and development of the data science provides new basis for estimation of the model parameters which influence the description and further prediction of the different physical processes.

In this study we present new approach which does not require vapor–liquid equilibrium calculations for binary mixtures. This approach is based on synergetic combination of global phase diagram technique and Pareto-based regression to reduce the uncertainty level caused by different sources of the data.

The global phase diagram is used to determine the type of phase behavior and derive analytical expressions to predict azeotrope and double azeotrope phenomena in terms of critical parameters of pure components.

To restore the EoS model parameters under uncertainty, the Pareto-optimality concept with fuzzy convolution scheme is applied. This approach is quite general and can be applied to other mathematical models, which describe a wide spectrum of phenomena of thermodynamic and phase behavior.

The azeotrope and double azeotrope criteria were elaborated and assessed for the imidazolium (IL)-based ionic liquids and industrial refrigerant mixtures. It was shown that IL doping leads to the breaking of azeotrope in binary refrigerant mixtures that open the way for the azeotrope refrigerant mixture separation technologies in order to remove the environmentally harmful substances.

The azeotrope phenomena in the refrigerant-IL mixtures are discussed, and conclusion about the highly improbable azeotrope blend formation for these systems is given. Azeotrope-breaking in the R134a-R1234yf mixture at IL doping is considered as an illustration of zeotrope behavior in the refrigerant-IL mixtures. Global phase behavior of ionic liquid-industrial refrigerant mixture is analyzed, and possible types of phase behavior according to the classification scheme of Scott and Van Konynenburg are established.