Modeling of Ionic Liquid Systems: Phase Equilibria and Physical Properties

1. Introduction

Ionic liquids (ILs) are a class of salts with a melting temperature below 100 °C, and the study of these compounds is considered priority by the U.S. Environmental Protection Agency. Due to their specific properties, which can be adjusted by changing either the cation or the anion, ILs have received great attention by the scientific community as potential replacements for volatile organic solvents (VOCs), and nowadays, ILs are starting to leave academic labs and find their way into a wide variety of industrial applications [1]. For example, ILs are used for the dispersion of nano-materials at IOLITEC, Air Products uses ILs instead of pressurized cylinders as a transport medium for reactive gases, ION Engineering is commercializing technology using ILs and amines for CO₂ capture and natural gas sweetening, and many others.

In order to apply these new compounds in different processes, the study of their physical properties, pure or mixed with other solvents, and phase equilibria (vapor-liquid, liquid-liquid, and solid-liquid equilibria) is crucial from a technological point of view. For example, density is fundamental to develop equations of state and it is also required for the design of different equipments, while viscosity is necessary for the design of processing units and pumping systems, and to study heat and mass transfer processes [2]. On the other hand, the refractive index can be used as a measure of the electronic polarizability of a molecule and can provide useful information when studying the forces between molecules or their behavior in solution [3].

As known, ILs also show an interesting potential to be used in separation processes and extraction media. Therefore, the knowledge of the mutual solubilities of molecular solvents and ILs prior to their industrial applications is also of primary importance. Moreover, many factors that control the phase behavior of these ionic salts with molecular solvents may be described from the phase equilibrium data.

However, as the number of possible ILs is enormous, this cannot be accomplished via experimental determination. Thus, it is very important to obtain models or empirical equations able to describe satisfactorily the experimental data.

In this chapter, a revision of the different equations applied for the modeling of physical properties of pure ILs and their mixtures, and phase equilibria of binary and ternary mixtures containing ILs, is presented and discussed. Future trends regarding the use of new models, namely equations of state accounting for association effects, are also focused.

2. Physical properties

2.1. Pure ionic liquids

Since ILs are relatively new compounds, experimental data on physical properties, such as density, viscosity, or refractive index of pure ILs and its mixtures with other solvents are required for the design of different equipment and processing units and very useful for developing accurate theoretical models.

Due to innumerable number of ILs that can be synthesized, experimental measurements are impractical for selection of a suitable IL for a specific application. Therefore, development of correlations and theoretical approaches allowing accurate modeling of IL-based systems is essential. This section shows the most common empirical equations used to correlate the temperature dependence of some of the physical properties of ILs.

For pure ILs, temperature dependence of physical properties such as density, speed of sound, or refractive index is very important for the successful and large-scale use of these compounds. Usually, this dependence is described using simple polynomial expressions, mainly equations of first, second and third order [4-6].

Several papers were also published concerning the experimental densities of pure ILs as a function of temperature and pressure [6-9]. The Tait equation [10] with four adjustable parameters is commonly used to fit these experimental data [6,8]. This equation is an integrated form of an empirical equation representative of the isothermal compressibility behavior versus pressure, and it can be expressed as:

$$\rho (T,p,C,B)=\frac{\rho (T,\mathrm{ }0.1\mathrm{ }\mathrm{M}\mathrm{P}\mathrm{a})}{1-C·Ln\left(\frac{B\left(T\right)+p}{B\left(T\right)+0.1\mathrm{ }\mathrm{M}\mathrm{p}\mathrm{a}}\right)}$$E1

where ρ (T, 0.1 MPa) represents the temperature dependence of density at 0.1 MPa. C is an adjustable parameter, and B(T) is commonly expressed as a second-order polynomial:

$$B\left(T\right)=\sum _{i=0}^2{B}_i·T^i$$E2

Regarding the variation of viscosity with temperature for pure ILs, a large number of empirical equations for correlating this property of pure fluids and mixtures can be found in literature [4,5,11-14]. The most commonly used equation is an Arrhenius-like law:

$$\mathrm{\eta }={\mathrm{\eta }}_{}·exp\left(-\frac{E_{\mathrm{a}}}{RT}\right)$$E3

where the viscosity at infinite temperature (η_∞) and the activation energy (E_a) are characteristic parameters generally adjusted from experimental data.

According to Seddon et al. [15], the Arrhenius law can generally be applied when the cation presents only a limited symmetry. If it is not the case, and especially in the presence of small and symmetrical cations with low molar mass, other equations such as the Vogel–Fulcher–Tamman (VFT) equation [16-18], or the modified VFT (mVFT) equation are recommended. This kind of expressions includes an additional adjustable temperature parameter (T₀) to the exponential term:

$$\mathrm{\eta }=A·exp\left(\frac{B}{T-T_0}\right)$$E4

$$\mathrm{\eta }=A·T^{0.5}·exp\left(\frac{B}{T-T_0}\right)$$E5

where A and B are also adjustable parameters.

Another empirical equation to correlate viscosity data with temperature was proposed by Litovitz [19]:

$$\mathrm{\eta }=A·exp\left(\frac{B}{R{T}^3}\right)$$E6

This equation is used at ambient pressures and has the advantage of containing only two fitting parameters. Comparing all these equations, in general, the best fits for the variation of viscosity with temperature for pure ILs are obtained with the VFT equation [5].

As reflected by Harris et al. [11], the general form of the pressure dependence of the shear viscosity of liquids is greater than exponential at moderate pressures and less than exponential at very high pressures. Taking this into account, these authors modified the VFT and Litovitz equations in order to include the temperature and pressure dependence of the viscosity for several pure imidazolium-based ILs. The new equations were defined as:

$$\mathrm{\eta }=exp\left(A+Bp+\frac{\left(C+Dp+E{p}^2\right)}{T^3}\right)$$E7

$$\mathrm{\eta }=exp\left(A´+B´+\frac{\left(C´+D´p+E´p^2\right)}{T-T_0}\right)$$E8

As for the density, a Tait-form equation is also used to correlate the pressure dependence of viscosity demonstrating good correlations [12,13].

$$\frac{\mathrm{\eta }(T,p,C,B)}{\mathrm{\eta }(T,0.1\mathrm{ }MPa)}=\mathrm{e}\mathrm{x}\mathrm{p}\left(C·Ln\left(\frac{B\left(T\right)+p}{B\left(T\right)+0.1\mathrm{ }\mathrm{M}\mathrm{p}\mathrm{a}}\right)\right)$$E9

where η(T, 0.1 MPa) represents the temperature dependence of viscosity at 0.1 MPa. C is an adjustable parameter, and B(T) is also commonly expressed as a second-order polynomial.

A hybrid Tait–Litovitz equation at elevated pressures (up to 126 MPa) was also presented in the literature to correlate the viscosity data for a series of room-temperature ILs [14]. The Litovitz equation is firstly used to correlate the data at ambient pressure, and then the Tait parameters are fitted for the higher pressures. This equation has the advantages of containing fewer fitting parameters than other models and simplicity of data analysis. The results show a good fit between the experimental data and those predicted by this equation.

2.2. Binary and ternary mixtures

In order to better understand the nature of ILs and design any future technological processes, detailed knowledge on the physical properties of ILs mixed with other solvents is required. During the last few years, the number of studies on thermophysical and thermodynamic properties of pure ILs and their mixtures with molecular solvents has increased significantly [4,20-24].

As for pure ILs, the dependence of the physical properties with temperature and composition is also correlated using empirical equations. In general, the change of density, speed of sound, refractive index and viscosity with composition is typically fitted to a polynomial expression although other more specific equations can be also found in literature. As example, the Connors and Wright equation [25] is employed to describe the variation of density with composition:

$$\rho ={\rho }_1-\left[1+\frac{b{x}_1}{\left(1-a{x}_1\right)}\right]\left(1-x_1\right)\left({\rho }_1-{\rho }_2\right)$$E10

where ρ₁ and ρ₂ are densities of pure compounds; x_i is the mole fraction of component i of the mixture; and a and b are fitting parameters. Although this equation was initially adopted only for the concentration dependencies of surface tension, Geppert-Rybczynska et al. [20] employed this equation for density getting a fit better than using any other simple polynomial.

As it is known, the above mentioned physical properties can be used to obtain the corresponding excess properties, which are generally fitted to a Redlich-Kister type equation [26]:

$$\delta Q=x_1\left(1-x_1\right)\sum _{\mathrm{i}=0}^M{A}_{\mathrm{i}}{\left({2x}_1-1\right)}^{\mathrm{i}}$$E11

where δQ is the excess property, x is the mole fraction, A_i are the adjustable parameters and M is the degree of the polynomial expansion.

An extended version of the Redlich-Kister equation, which takes into account the dependence on composition and temperature simultaneously, is also used to fit the excess properties [22]:

$$\delta Q=x_1\left(1-x_1\right)\frac{\sum _{\mathrm{i}=0}^M{A}_{\mathrm{i}}{\left({2x}_1-1\right)}^{\mathrm{i}}}{1+\sum _{j=1}^N{B}_j{\left(2x_1-1\right)}^j}$$E12

In order to take into account the influence of temperature on the excess properties, all the coefficients A_i and B_i are usually expressed as a second-order polynomial.

Density, refractive index, and viscosity data for ternary mixtures containing ILs are also common in the literature [27-29] and they are usually fitted to polynomial expressions. In these cases, weight fractions are often used instead of mole fractions, due to the large difference of molar mass between ILs and most organic solvents.

For the modeling of the excess properties, such as excess molar volumes, viscosity deviations, refractive index deviations or excess free energies of activation of viscous flow, the use of empirical equations is commonly adopted. Although the Redlich-Kister equation is also applied to correlate the excess properties for ternary systems containing ILs [29], the equations more widely employed are those proposed by Cibulka [30], Singh et al. [31] and Nagata and Sakura [32]. All the correlative models are capable of representing the behavior of the ternary mixtures with a higher or lesser degree of accuracy, although that developed by Cibulka usually leads to a better agreement with experimental data [28].

Cibulka equation:

$$Q_{123}^{\mathrm{E}}=Q_{12}^{\mathrm{E}}+Q_{13}^{\mathrm{E}}+Q_{23}^{\mathrm{E}}+x_1{x}_2{x}_3\left(A+B{x}_1+C{x}_2\right)$$E13

Singh et al. equation:

$${\mathrm{ }\mathrm{ }\mathrm{ }Q}_{123}^{\mathrm{E}}=Q_{12}^{\mathrm{E}}+Q_{13}^{\mathrm{E}}+Q_{23}^{\mathrm{E}}+{Ax}_1{x}_2{x}_3+B{x}_1\left(x_2-x_3\right)+C{x}_1^2{\left(x_2-x_3\right)}^2$$E14

Nagata and Sakura equation:

$$Q_{123}^{\mathrm{E}}=Q_{12}^{\mathrm{E}}+Q_{13}^{\mathrm{E}}+Q_{23}^{\mathrm{E}}+x_1{x}_2{x}_{3}A$$E15

where A, B and C are fit parameters and $Q_{\mathrm{i}\mathrm{j}}^{\mathrm{E}}$ is the contribution to the excess property of the constituent binary mixtures evaluated by Redlich-Kister equation.

3. Phase equilibria

Despite the large number of published articles and the broad fields of applications, there has not been a model explicitly derived for phase equilibria of ILs. This lack of models intended for systems containing ILs has forced researchers to use the equations available. But these equations were intended for ionic solutions, or not intended for ions at all. Thus, a model to account for a medium which is composed of ions, without a molecular solvent, is still needed. Nevertheless, the phase equilibria of ILs and their mixtures are being modeled in the literature. The models used will be described ahead. In general, those based on the excess Gibbs energy such as Wilson, NRTL or UNIQUAC were the first to appear. Lately, models with modifications for association (UNIQUAC ASM, NRTL1) and for electrolytes (PDH, e-NRTL) were also applied for this kind of systems in order to improve the results obtained with the initial models. Besides, Equations of State (EoS) have also been applied, especially for mixtures with gases and for broad ranges of pressure.

3.1. g^E-based models

Most of the g^E models available in the literature are for non-electrolytes. Thus, many authors have been using these models, or models for electrolyte solutions, for the phase equilibria of systems containing ILs. For example, for binary systems containing ILs it is common to use the model developed by Debye-Hückel (which was derived for small salt concentrations), although it is not recommendable for solutions at high ionic concentration.

The models for the correlation of these experimental data can be split into two main groups: i) ion-interaction Pitzer models and ii) local composition models.

i. The model developed by Pitzer has created a new generation of theories which use multiparameter regression. In the ion-interaction Pitzer model [33] the ion-interaction parameters are dependent on temperature and pressure, and it takes into account the Debye-Hückel constant.

The three-parameter Pitzer-ion interaction model has been successfully used for modeling vapor-liquid data of mixtures of ILs with water [34-37] or with alcohol [38] and has the following form for a binary 1:1 electrolyte solution:

ϕ -1 = f^ϕ + mB^ϕ + m²C^ϕ,E16

where

f^ϕ = -A_ϕ I^1/2 / (1 + bI^1/2),E17

A_ϕ = (1/3)(2πN_Ad_s)^1/2 (e² / 4πε₀εkT)^3/2,E18

B^ϕ = β ⁽⁰⁾ + β ⁽¹⁾ exp(-α₁I ^1/2)E19

In these equations, β ⁽⁰⁾, β ⁽¹⁾ and C^ϕ are ion-interaction parameters of the Pitzer model that are dependent on temperature and pressure, and A_ϕ is the Debye-Hückel constant for the osmotic coefficient on the molal scale. Also, N_A is the Avogadro number, e is the proton charge, ε₀ is the permittivity of vacuum, and k is the Boltzman constant. The term I is the ionic strength in molality, I = 1/2 Σm_iz_i², where m_i is the molality of i^th ion and z_i is the absolute value for i^th ionic charge. The remaining symbols have their usual meanings. The values for constants b and α were b = 1.2 kg^1/2 mol^-1/2 and α = 2 kg^1/2 mol^-1/2, respectively.

In the last years, the Extended Pitzer model of Archer [39,40], in which the third adjustable parameter in the Pitzer model is replaced by a two-parameter function depending on the ionic strength, has demonstrated its accuracy in modeling binary Vapor-Liquid Equilibria (VLE) of systems containing ILs. In this model, the equation for Β^ϕ is extended to:

B^ϕ = β ⁽⁰⁾ + β ⁽¹⁾ exp(-α₁I ^1/2) + β ⁽²⁾ exp(-α₂I ^1/2),E20

and a new equation is introduced:

C^ϕ = C ⁽⁰⁾ + C ⁽¹⁾ exp (-α₃I ^1/2)E21

In the previous equations, the ion interaction parameters of the extended Pitzer model of Archer are β ⁽⁰⁾, β ⁽¹⁾, β ⁽²⁾, C ⁽⁰⁾ and C ⁽¹⁾, dependent on temperature and pressure, and α₁, α₂, α₃, and b can be adjustable parameters or kept fixed at constant values; they are usually fixed to the values: α₁ = 2 kg^1/2 mol^-1/2, α₂ = 7 kg^1/2 mol^-1/2, α₃ = 1 kg^1/2 mol^-1/2 and b = 3.2 [41-44], although other values can be used [45]. This model is widely used in literature for binary mixtures containing ILs and water or alcohol or acetonitrile [46-55] with very satisfactory results.

ii. The local composition (LC) models give a better empirical description and have physical meaning for the correlation of osmotic and activity coefficients. There are several LC models reported in literature for the modeling of phase equilibria experimental data; such as UNIversal QUAsiChemical (UNIQUAC) [56], Non-Random Two Liquids (NRTL) [57], electrolyte NRTL (e-NRTL) [58], Non-Random Factor (NRF) [59], modified NRTL (MNRTL) [60], Mean Spherical Approximation NRTL (MSA-NRTL) [61] or Extended Wilson (EW) [62]. Furthermore, models with modifications for association such as UNIQUAC ASM or NRTL 1 are also applied for this kind of systems, although their use is less common [63,64].

In local composition models it is assumed that the activity coefficient is composed by two terms: a long-range contribution (LR) and a short-range contribution (SR):

$$\ln {\gamma }_i=\ln {\gamma }_i^{LR}+\ln {\gamma }_i^{SR}$$E22

The well-known models of UNIQUAC and NRTL calculate the activity coefficients as follows:

ii.i In the UNIQUAC model [56], the long-range contribution term is expressed as:

$$\text{ln }{\gamma }_{\text{i}}^{\text{LR,UNIQUAC}}=\text{ ln}\left(\frac{\varphi }{x_{\text{i}}}\right)+\frac{\text{z}}{\text{2}}{q}_{\text{i}}\ln \left(\frac{{\theta }_{\text{i}}}{{\varphi }_{\text{i}}}\right)+l_{\text{i}}\text{ - }\frac{{\varphi }_{\text{i}}}{x_{\text{i}}}{\displaystyle \sum _{\text{j}=\text{1}}^{\text{N}}{x}_{\text{j}}{l}_{\text{j}}}$$E23

where the term l_i is defined as a function of the external surface of the molecule and the bond segments:

$$l_{\text{i}}=\frac{\text{z}}{2}(r_{\text{i}}-q_{\text{i}})-(r_{\text{i}}-1)$$E24

and the coordination number, z, is set to 10.

The short-range contribution in this model is expressed as:

$$\text{ln }{\gamma }_{\text{i}}^{\text{SR,UNIQUAC}}=q_{\text{i}}\left[\text{1 - ln}\left({\displaystyle \sum _{\text{j}}{\theta }_{\text{j}}{\tau }_{\text{ji}}}\right)\text{ -}{\displaystyle \sum _{\text{j}=\text{1}}^{\text{N}}\frac{{\theta }_{\text{j}}{\tau }_{\text{ij}}}{{\displaystyle \sum _{\text{k}}{\theta }_{\text{k}}{\tau }_{\text{kj}}}}}\right]$$E25

The expression for the energy parameter, τ_ij, the volume fraction, ϕ_i, and the surface fraction, θ_i, are given for the next equations:

$${\tau }_{ji}=\exp \left(-\frac{u_{ji}-u_{ii}}{RT}\right)=\exp \left(\frac{\Delta u_{ij}}{RT}\right)$$E26

$${\theta }_{\text{i}}=\frac{q_{\text{i}}{x}_{\text{i}}}{{\displaystyle \sum _{\text{j}}{q}_{\text{j}}{x}_{\text{j}}}}{\varphi }_{\text{i}}=\frac{r_{\text{i}}{x}_{\text{i}}}{{\displaystyle \sum _{\text{j}}{r}_{\text{j}}{x}_{\text{j}}}}$$E27

where u_ij and u_ji are the energetic parameters, and r_i and q_i are the structural parameters corresponding to relative volume and surface area of the component i, respectively.

This model has been used for VLE of binary systems containing ILs [65]. The equation has also been used successfully for Liquid-Liquid Equilibria (LLE) of ternary systems including an IL [66,67], and even for Solid-Liquid Equilibria (SLE) of binary systems [68-70]. The main problem using the UNIQUAC model is the need for structural parameters r_i and q_i, the volume and surface area of the component. Despite for molecular components these parameters can be obtained easily (even from group-contribution data), in the case of ILs values for both anion and cation are calculated and summed. Different procedures have been proposed [66,68] and used successfully.

ii.ii In the NRTL model [57], the activity coefficients are calculated as follows:

$$\ln {\gamma }_i^{NRTL}=\left(\frac{{\displaystyle \sum _{j=1}^n{\tau }_{ji}{x}_j{G}_{ji}}}{{\displaystyle \sum _{k=1}^n{x}_k{G}_{ki}}}\right)+{\displaystyle \sum _{j=1}^n\left(\frac{x_j{G}_{ij}}{{\displaystyle \sum _{k=1}^n{x}_k{G}_{kj}}}\right)}\left({\tau }_{ij}-\frac{{\displaystyle \sum _{m=1}^n{\tau }_{mi}{x}_m{G}_{mi}}}{{\displaystyle \sum _{k=1}^n{x}_k{G}_{kj}}}\right)$$E28

where

$${\text{G}}_{\text{ij}}\text{= exp (-}{\alpha }_{\text{ij}}{\tau }_{\text{ij}}\text{)}$$E29

$${\tau }_{ij}=\frac{g_{ij}-g_{jj}}{RT}=\frac{\Delta g_{ij}}{RT}{\tau }_{ji}=\frac{g_{ji}-g_{ii}}{RT}=\frac{\Delta g_{ji}}{RT}$$E30

where x represents the mole fraction, g_ij is an energy parameter that characterizes the interaction of species i and j, R is the gas constant, T is the absolute temperature, and the parameter α_ij = α_ji is related to the nonrandomness in the mixture. Although α_ij can be adjusted, it can also be considered fixed in a value, usually between 0.2 and 0.47.

This is the local composition model most widely used in literature for binary and ternary systems containing ILs, regarding VLE [71-79] with alcohols or water, LLE of binary [80,81] and ternary systems [67,82-85] with ethanol or hydrocarbons or water and also SLE [68-70].

For the explanation of the following models, the Pitzer-Debye-Hückel (PDH) equation [33] has been used as the long-range term on a mole fraction scale as proposed by Chen et al. [86], assuming that the studied mixtures involve completely dissociated electrolytes. For the solvent, the formulation of this equation is:

$$ln{\gamma }_1^{PDH}=\frac{2A_{\varphi }{I}_x^{\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}{\sqrt{M_s}\left(1+\rho I_x^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right)}$$E31

in which ρ is the number density of ionic species and the term I is the ionic strength on a mole fraction basis, I = 1/2Σx_iz_i².

The short-range contribution calculated using different models, such as electrolyte NRTL (e-NRTL) [58], non-random factor (NRF) [59], modified NRTL (MNRTL) [60], mean spherical approximation NRTL (MSA-NRTL) [61] or Extended Wilson (EW) [62] are explained below.

ii.iii In the e-NRTL model [58], the short-range contribution for the activity coefficient is calculated as:

$$\ln {\gamma }_1^{e-NRTL}=2X_c^2\left[\frac{2{\tau }_{ca,m}\exp \left(-2\alpha {\tau }_{ca,m}\right)}{{\left(X_1+2X_c\exp \left(-\alpha {\tau }_{ca,m}\right)\right)}^2}+\frac{{\tau }_{m,ca}\exp \left(-\alpha {\tau }_{m,ca}\right)}{{\left(X_1\exp \left(-\alpha {\tau }_{m,ca}\right)+X_c\right)}^2}\right]$$E32

where α is the nonrandomness factor (usually set to 0.2) and X_i is the effective mole fraction of the component i, calculated as X_i = j_ix_i with j_i = z_i for ions and j_i = 1 for solvent. In this model, τ_ca,m and τ_m,ca are the adjustable parameters.

Examples of the correlation of VLE for binary mixtures [38,73-75] and ternary mixtures [73-75] containing ILs and ethanol or water can be found in literature. Nevertheless, the literature is scarce in examples for SLE [70] and LLE [80,81].

ii.iv The NRF model [59] calculates the short-range contribution for the activity coefficient of the solvent as:

$$\ln {\gamma }_1^{NRF}=X_c\left\{\begin{array}{l}\left(\frac{\nu X_c{\lambda }_E}{z_c{\nu }_c}\right)\left(1+\frac{\exp \left(-{\lambda }_E/Z\right)}{{\left(X_c\exp \left(-{\lambda }_E/Z\right)+X_1\right)}^2}\right)-\\ \frac{X_1{\lambda }_S\left(\left(\nu /z_c{\nu }_c\right)-2\exp \left(-{\lambda }_S/Z\right)\right)}{2X_c\exp \left(-{\lambda }_S/Z\right)+X_1}·\left(2-X_1\left(1+\frac{1}{2X_c\exp \left(-{\lambda }_S/Z\right)+X_1}\right)\right)\end{array}\right\}$$E33

where ν and ν_c are the total number of ions into which the salt dissociates and the number of cations in one mole of the salt, respectively, and Z is the coordination number (usually set to 8). In this model λ_E and λ_S are the adjustable parameters; its use is not common for describing the vapor-liquid equilibria of systems containing ILs [38,77].

ii.v The short-range contribution in the MNRTL model has been developed by Jaretum and Aly [60], and Sardroodi et al. [87] used it in the presented form for the first time:

$$ln{\gamma }_1{}^{MNRTL}=4X_c{}^2[{\tau }_{ca,m}{W}_{ca,m}/{\left(X_1+2W_{ca,m}{X}_c\right)}^2+{\tau }_{m,ca}\left(W_{m,ca}–1\right)/{\left(X_c+X_1{W}_{ca,m}\right)}^2]$$E34

where τ_ca,m and τ_m,ca are the parameters of the model, X_i is the effective mole fraction, and the next expression is assumed:

$$W_i=exp(-\alpha {\tau }_i+w_i)$$E35

where α is the nonrandomness factor (usually fixed to 0.2) and ω_ca,m and ω_m,ca are the adjustable parameters.

Among the local composition models, it is one of those giving lower deviations for the modeling of VLE of systems with ILs [34-37,47-54,77].

ii.vi The equation for calculating the short-range contribution for the activity coefficient of the solvent given by the MSA-NRTL model [61] is as follows:

$$\ln {\gamma }_1^{MSA-NRTL}=\frac{A}{B^2}+{\nu }_c{\nu }_a\left(C_c+C_a\right)+\nu {\tau }_{mc,ac}^2{x}_S^2{x}_S^i\left(D_c+D_a\right)$$E36

taking into account that:

$$A=({\nu }_c{p}_{cm}{\tau }_{cm}+{\nu }_a{p}_{am}{\tau }_{am})({\nu }_c{p}_{cm}+{\nu }_a{p}_{am})$$E37

$$B\text{}={\nu }_c{p}_{cm}+{\nu }_a{p}_{am}+x_S^1$$E38

$$C_c\text{}=p_{cm,ac}{\tau }_{mc,ac}/{\left({\nu }_a+x_S^i{p}_{mc,ac}\right)}^2$$E39

$$C_a\text{}=p_{ma,ca}{\tau }_{ma,ca}/{\left({\nu }_c+x_S^i{p}_{ma,ca}\right)}^2$$E40

$$D_c\text{}=\frac{{\nu }_c{p}_{mc,ac}}{{\left({\nu }_a+x_S^i{p}_{mc,ac}\right)}^2\left(1-\frac{{\nu }_a\alpha {\tau }_{mc,ac}}{{\nu }_a+x_S^i{p}_{mc,ac}}\right)}$$E41

$$D_a\text{}=\frac{{\nu }_a{p}_{ma,ca}}{{\left({\nu }_c+x_S^i{p}_{ma,ca}\right)}^2\left(1-\frac{{\nu }_c\alpha {\tau }_{ma,ca}}{{\nu }_c+x_S^i{p}_{ma,ca}}\right)}$$E42

being

$${\text{p}}_{\text{i}}\text{ = exp (-}\alpha {\tau }_{\text{i}}\text{)}$$E43

$${\tau }_{\text{ma,ca}}\text{ =}{\tau }_{\text{mc,ac}}\text{ +}{\tau }_{\text{am}}\text{ -}{\tau }_{\text{cm}}$$E44

$${\tau }_{mc,ac}={\tau }_{mc,ac}^1+{\tau }_{mc,ac}^2{X}_1$$E45

$$x^i =x_1/x_S=1/x_S-\nu$$E46

in which x₁ and x_S are the mole fraction of solvent and salt, respectively.

The correlation using this model it is not common for the treatment of the VLE data of systems containing ILs [77].

ii.vii The EW model was presented by Zhao et al. [62], and it describes the short-range contribution for the activity coefficient as:

$$\begin{array}{l}\ln {\gamma }_1^{EW}=-C\left(\ln \left(\frac{X_1+2X_c\exp \left(-E_{e1}/CRT\right)}{X_1+2X_c}\right)\right)+\\ 2X_c\left(\frac{\left(1-\exp \left(-E_{e1}/CRT\right)\right)X_1}{\left(X_1+2X_c\exp \left(-E_{e1}/CRT\right)\right)\left(X_1+2X_c\right)}+\frac{X_c\left(\exp \left(-E_{1e}/CRT\right)-1\right)}{\left(X_1\exp \left(-E_{1e}/CRT\right)+X_c\right)\left(X_1+X_c\right)}\right)\end{array}$$E47

where C is the coordination number (set to 10 [77]). In this model, E_e1 and E_1e are the adjustable parameters. This model has been checked for its use in VLE of aqueous systems with ILs [77].

Model	References
VLE	LLE	SLE
Pitzer	[34-38]
Extended Pitzer of Archer	[46-55]
UNIQUAC	[65]	[66,67]	[68-70]
NRTL	[71-79]	[67,80-85]	[70-72]
e-NRTL	[38,73-75]	[80,81]	[70]
NRF	[38,77]
MNRTL	[34-37,47-54,77]
MSA-NRTL	[77]
EW	[77]

Table 1.

Literature examples for the modeling of phase equilibria in systems containing ILs

A list of representative examples found in literature for the modeling of phase equilibria in systems containing ILs with the above mentioned models is presented in Table 1.

Among these correlation models, those which have demonstrated to give the best results in VLE are the Extended Pitzer model of Archer, the NRTL and the MNRTL models. Regarding LLE and SLE, there are less examples available, specially comparing different models. Nevertheless, it is clear that NRTL is, by far, the most used equation. It is also important to highlight that differences in performance among the models are small.

3.2. Equations of state

Equations of state (EoS) are powerful tools, which can be used to describe the properties of pure fluids or their mixtures. In the last 10 years, this kind of models has been widely applied to describe the properties of pure ILs, as well as to model the phase equilibrium (VLE and LLE) of mixtures containing them.

3.2.1. Peng-Robinson

The Peng-Robinson EoS was developed in 1976 by Peng and Robinson [88] and can be expressed as:

$$P=\frac{RT}{V-b}-\frac{a}{V\left(V+b\right)+b\left(V-b\right)}$$E48

with:

$$b=0.077796\frac{R{T}_c}{P_c}$$E49

$$a=0.457235\frac{R^2{T}_c^2}{P_c}\alpha$$E50

$$\alpha ={\left[1+m\left(1-\sqrt{\frac{T}{T_c}}\right)\right]}^2$$E51

$$m=0.37464+1.54226\omega -0.26992{\omega }^2$$E52

where T_c is the critical temperature, P_c is the critical pressure and ω is the acentric factor. For most ILs, the critical properties and the acentric factor are impossible to determine experimentally, because they start to decompose even before the temperature reaches the boiling point. For this reason, Valderrama and co-workers have applied an extended group contribution method, the modified Lydersen-Joback-Reid method, to determine the critical properties, boiling temperatures and acentric factors of several ILs [89-92]. This method only requires knowledge of the structure of the ILs and its molecular weight. Since there are no experimental data available for these properties, the accuracy of the method is verified by comparing calculated liquid densities of the ILs to experimental data available. The results show that the method is sufficiently accurate for several applications, with average absolute deviations (AAD) between calculated and experimental liquid densities in the range of 5 to 6%. The properties determined by this research group have been widely used by the scientific community when modeling the properties of ILs and its mixtures, using this or other equations.

For mixtures of fluids, mixing rules have to be applied to parameters a and b, which imply the use of one or more binary interaction parameters. For systems containing ILs, the Wong-Sandler mixing rules have been widely used [93-96], as well as the quadratic [97,98], van der Waals [94], [99,100] and Mathias-Klots-Prausnitz mixing rules [101,102].

The application of the Peng-Robinson EoS to systems with ILs has been mainly focusing on the VLE with CO₂ and other gases. For example, Shin et al. [103] modeled the high-pressure solubilities of CO₂ in ILs of the family [C_nmim][Tf₂N] (1-alkyl-3-methylimidazolium bis(trifluoromethylsulfonyl) imide) using this EoS and quadratic mixing rules with two temperature-dependent binary interaction parameters. They obtained AAD between calculated and experimental equilibrium pressures between 10 and 14%, for which they conclude that the model can satisfactorily predict the solubility of high-pressure CO₂ in this family of ILs over a wide range of pressures up to the supercritical region of CO₂. Also Álvarez and Aznar [93] used the Peng-Robinson EoS to model the VLE of binary systems composed of IL + supercritical CO₂ or CHF₃ and IL + hydrocarbons. In this work, the van der Waals and Wong-Sandler mixing rules were used, with UNIQUAC and NRTL models used as excess Gibbs energy models in the Wong-Sandler mixing rules. Their results present AAD in pressure between 2 and 24% for low pressures and between 7 and 52% for high pressures. The authors conclude that the EoS is not able to represent the data at high pressures; however, it performs well at low pressures. Later, Álvarez et al. [104] used the Peng-Robinson EoS with the Wong-Sandler mixing rules, once again with UNIQUAC and NRTL models for the excess Gibbs energy, but also using the COSMO-SAC model. They modeled the isobaric VLE of 1-ethyl-3-methylimidazolium ethylsulfate ([C₂mim][EtSO₄]) with propionaldehyde or valeraldehyde. When using UNIQUAC or NRTL models, they obtain AAD below 0.15% for both systems. Since the COSMO-SAC model is a predictive model, i.e., it does not require any adjustable binary interaction parameters, the authors regard the application of this model with the Peng-Robinson EoS and the Wong-Sandler mixing rules as pure predictive results. The AAD obtained in this way was 0.3 and 2% for each system, respectively. As a final example, Ren and Scurto [99] used the Peng-Robinson EoS with the van der Waals one fluid mixing rule to model the VLE and VLLE of imidazolium-based ILs and the refrigerant gas 1,1,1,2-tetrafluoroethane. The authors obtained AAD values between 0.4 and 7.4% for VLE and concluded that the model is able to represent the bubble-point data with excellent agreement. However, when using the interaction parameters regressed from the VLE data alone to predict compositions of the VLLE transition to LLE, only satisfactory results were obtained for all systems.

3.2.2. Soave-Redlich-Kwong

A Soave modification of the Redlich-Kwong EoS [105] has been frequently applied to systems with ILs. The Soave-Redlich-Kwogn (SRK) EoS was introduced in 1972 by Giorgio Soave, and can be expressed as:

$$P=\frac{RT}{V-b}-\frac{a\left(T\right)}{V\left(V+b\right)}$$E53

with:

$$a\left(T\right)=0.427480\frac{R^2{T}_c^2}{P_c}\alpha \left(T\right)$$E54

$$b=0.08664\frac{R{T}_c}{P_c}$$E55

In the works involving systems with ILs, the temperature dependent part of the a parameter has been modeled by the empirical form:

$$\alpha \left(T\right)=\sum _{k=0}^{\le 3}{{\beta }_k\left(\frac{T_c}{T}-\frac{T}{T_c}\right)}^k$$E56

Coefficients β_κ are usually determined to reproduce the vapor pressure of the pure compound. However, since there are no available experimental vapor pressure data for most ILs, the coefficients are treated as adjustable fitting parameters in the calculations [106]. The authors have found that for ILs, only one adjustable parameter, β₁, is sufficient (with β₀=1 and β₂= β₃=0). In the case of this EoS, modified van der Waals-Berthelot mixing rules have been used, with four binary interaction parameters.

Shiflett and Yokozeki [106] first used this EoS to correlate the solubility of CO₂ in ILs l-butyl-3-methylimidazolium hexafluorophosphate ([C₄mim][PF₆]) and l-butyl-3-methylimidazolium tetrafluoroborate ([C₄mim][BF₄]). The authors obtained an excellent fit between experimental and calculated solubility data, with standard deviations of 17 and 10.5 kPa for each system, respectively. They used the same model to correlate the VLE of CO₂ and other gases (like ammonia [107,108] or SO₂ [109]) in different ILs. Later, the same authors modeled the solubility of water in several different ILs using the same EoS [110]. The four binary interaction parameters were determined using binary VLE data from the literature, having obtained standard deviations lower than 0.6 kPa for all systems. The graphical results presented by the authors show a very good agreement between experimental and calculated solubilities of water in all ILs considered. The same research group has also used the SRK EoS to study the phase behavior of ternary mixtures CO₂/H₂/[C₄mim][PF₆] [111], CO₂/SO₂/1-butyl-3-methylimidazolium methyl sulfate ([C₄mim][MeSO4]) [112], CO₂/H₂S/[C₄mim][PF₆] [113] and CO₂/H₂S/[C₄mim][MeSO₄] [114]. For all the previous studies the authors obtained good agreement between calculated and experimental data. Very recently, Shiflett et al. [115] have also modeled the ternary system N₂O/CO₂/1-butyl-3-methylimidazolium acetate ([C₄mim][Ac]) with the purpose of understanding the separation of N₂O and CO₂ using ILs. They determined the binary interaction parameters using VLE data (either their own or from the literature) for the pairs N₂O/[C₄mim][Ac], CO₂/[C₄mim][Ac] and N₂O/CO₂. Unlike what happened in the previously mentioned ternary systems, the phase behavior prediction of the ternary system of N₂O/CO₂/[C₄mim][Ac] may not be guaranteed based on the binary interaction parameters alone. This happens because for systems containing mixtures with the chemical complex formation (hydrogen-bonding, charge-transfer complex, etc.) as in the case of the binary CO₂/[C₄mim][Ac] here present, the third component (in this case, N₂O) may interfere with the binary interactions of CO₂/[C₄mim][Ac], originating a decrease or an increase of the complex formation and changing the binary interaction parameters of this pair in the ternary system. Consequently, one of the binary interaction parameters of the pair CO₂/[C₄mim][Ac] was modified so as to fit the experimental ternary VLE data. This way, the model provides an excellent agreement between experimental and calculated data for the ternary system.

3.2.3. Statistical associating fluid theory

The original Statistical Associating Fluid Theory (SAFT) was developed in 1989 [116] based on Wertheim’s first-order thermodynamic perturbation theory [117-120]. Its main advantage over the traditional cubic EoS, is that it takes into account the structure of the molecule, similarly to group contribution models. It regards the molecules as chains of hard-spheres, which contain multiple association sites. Several variations of the SAFT EoS have been used to model the phase behavior of systems containing ILs: tPC-PSAFT [121-125], soft-SAFT [126-128], hetero-SAFT [129], PCP-SAFT [130] and PC-SAFT [131].

Kroon et al. [121] applied the tPC-PSAFT model to the phase behavior of IL + CO₂ systems and later Karakatsani et al. [122] applied it to the correlation of the solubility of CO₂, CO, O₂ and CHF₃ in the same ILs. Pure component parameters for the ILs were estimated from ILs experimental thermodynamic data (density, enthalpy, and entropy of dissolution of CO₂) and physicochemical data for the constituent ions (size, polarizability, number of electrons). The cross-association parameters were estimated from enthalpy and entropy data for the dissolution of CO₂ in the ILs. A temperature-dependent binary interaction parameter k_ij was adjusted in order to fit the model to experimental VLE data. The results of the model were in good agreement with experimental data. Later, Karakatsani and Economou [123] used the same parameters to model the VLE of the ternary system CO₂/acetone/[C₄mim][PF₆], also obtaining satisfactory agreement between calculated and experimental data. Additionally, Economou et al. [124] applied the same model to describe the VLE of the binary [C₈mim][BF₄]/benzene and the ternary CO₂/H₂O/[C₄mim][NO₃] (1-butyl-3-methylimidazolium nitrate) systems. For binary mixture calculations, a binary interaction parameter was fitted to the experimental data and the model was capable of accurately correlating the phase equilibria of the binary and of the ternary IL mixtures with polar solvents. In a more recent work, Karakatsani et al. [125] proved that the model can accurately predict the phase equilibrium of non-polar solvent/IL mixtures, without the use of any adjustable binary interaction parameter, by applying it to binary and ternary mixtures of ILs with organic solvents and water. They also concluded that in the case of aqueous solvents, the dissociation of IL has to be incorporated explicitly into the model in order to obtain a good correlation with the experimental data.

Andreu and Vega [126] used the soft-SAFT EoS to describe the solubility of CO₂ in ILs. They modeled the families of ILs [C_nmim][BF₄] and [C_nmim][PF₆] as Lennard-Jones chains with one associating site in each molecule. The chain length, size and energy parameters of the ILs were obtained by fitting the model predictions to available density data, obtaining AAD values lower than 0.2%. For the association parameters of ILs, values previously used for alkanols were adopted. The authors found that the model correlations and experimental data for VLE are in good agreement. They later used the same model to describe the solubility of hydrogen, CO₂ and xenon in ILs of the family [C_nmim][Tf₂N] [127]. In this case, the ILs were modeled as Lennard-Jones chains with three associating sites in each molecule. The pure component parameters for ILs were obtained as in the previous work, and a good description of experimental solubilities was obtained. Recently, Llovell et al. [128] re-parameterized the model for the [C_nmim][Tf₂N] family of ILs and modeled the solubilities of methanol and ethanol. Good agreement was found between the predictions of the model and the experimental data, with AAD values below 5% for methanol and around 10% for ethanol. They also modeled the VLE of IL/water systems, with the model providing good agreement with the experimental data, with AAD values between 6 and 12%. As for the LLE of mixtures of ILs with water, the authors were not able to obtain reasonable predictions, having obtained significant deviations in the water-rich phase. Quantitative agreement was achieved by using two adjustable parameters, which were adjusted only to the water-rich phase of the aqueous [C₄mim][Tf₂N] mixture and used in a predictive manner for the IL-rich phase and for the aqueous mixtures with [C₂mim][Tf₂N] and [C₆mim][Tf₂N].

Ji and Adidharma [129] used the heterosegmented SAFT (hetero-SAFT) to describe the solubility of CO₂ in the families of ILs [C_nmim][BF₄], [C_nmim][PF₆] and [C_nmim][Tf₂N]. The molecules of ILs were divided into groups representing the alkyl chain, the cation head and the anion. To account for the electrostatic/polar interactions between cation and anion, the spherical segments representing the cation head and the anion were assumed to have one association site each, which can only associate to each other. The parameters for the alkyl chains were obtained from those of the corresponding n-alkanes and the parameters for groups representing the cation head and the anion, including the two association parameters, were fitted to experimental IL density data. The model was capable of satisfactorily describing the solubility of CO₂ in the ILs studied.

Finally, Paduszynski et al. [130] used the Perturbed-Chain Polar Statistical Associating Fluid Theory (PCP-SAFT) to model the LLE of the IL 1-methyl-1-propylpiperidinium bis(trifluoromethylsulfonyl)imide [C₃mpip][NTf₂] with several alkan-1-ols. They modeled the IL as strongly associating molecules, with symbol A₁ representing a positive site which corresponds to the nitrogen atom on the cation and its proximity, and symbol B₁ representing a negative site which corresponds to the delocalized charge due to the oxygen molecules on the anion. They defined each type of associating site in an identical way; however, they only allowed A₁B₁ interactions to take place. Additionally, they assumed that each molecule has 5 positive sites of type A₁ and 5 negative sites of type B₁. The pure component parameters of the IL were determined by fitting to liquid density and total solubility parameter data. Self-association of IL and alkan-1-ols, as well as cross-association, were accounted for and one linearly temperature-dependent binary interaction parameter was needed in order to obtain qualitative agreement between calculated and experimental data. AAD values obtained were between 0.6 and 5%. On a different work, Paduszynski and Domanska [131] modeled the LLE of systems composed by piperidinium-based ILs ([C₃mpip][NTf₂] and [C₄mpip][NTf₂]) and several aliphatic hydrocarbons using the Perturbed-Chain Statistical Associating Fluid Theory (PC-SAFT) model. Pure component parameters for ILs were obtained as in the previous work and they used activity coefficients at infinite dilution of hydrocarbons in ILs reported in the literature to optimize the binary interaction parameter, which was again considered to be linearly temperature-dependent, with two adjustable parameters. In the calculation of the LLE of the systems, they obtained AAD values between 1.5 and 8%. Additionally, the authors decided to test their approach for the cross-associating systems [C₃mpip][NTf₂] + 1-pentanol and [C₃mpip][NTf₂] + water. Regarding the LLE of [C₃mpip][NTf₂] + 1-pentanol, the authors claim that the experimental data is well described by the model, although deviations between calculated and experimental compositions increase as the temperature increases and the model ends up over predicting the upper critical solution temperature by about 15 K. The AAD values obtained for this system were 52% for the IL-rich phase and 1.7% for the alcohol-rich phase. As for the system [C₃mpip][NTf₂] + water, the authors conclude that the PC-SAFT model is surprisingly good at describing the experimental data, but only when the IL-rich phase is considered alone. For the water-rich phase, the model predicts solubilities of IL much lower than those observed experimentally by several orders of magnitude. The authors justify this fact because the molecular model chosen for the IL is not appropriate for dilute solutions of IL in water, in which case the cation and anion of the IL are probably dissociated, which results on increased solubility.

3.2.4. Other EoS

Several other EoS have been used to model systems with ILs. For example, Tsioptsias et al. [132] used the Non-Random Hydrogen-Bonding (NRHB) model [133] to describe the phase behavior of binary systems containing ILs of the family [C_nmim][Tf₂N], obtaining good agreement between model correlations and experimental data. Wang et al. [134] used the square well chain fluid (SWCF) EoS [135] to model the solubilities of gases such as CO₂, C₃H₆, C₃H₈ and C₄H₁₀ in several ILs. Breure et al. [136] used a group contribution EoS to study the phase behavior of binary systems of ILs of the families [C_nmim][PF₆] and [C_nmim][BF₄] with CO₂, also obtaining good agreement between model predictions and experimental data. Very recently, Maia et al. [137] applied the Cubic Plus Association (CPA) EoS [138], which combines the SRK EoS with an advanced association term similar to that of the SAFT type models, to describe the VLE with CO₂ and the LLE with water of ILs [C₂mim][Tf₂N] and [C₄mim][Tf₂N]. Good agreement was obtained between calculated and experimental data, even though smaller AAD percentage values were obtained for the VLE than for the LLE.

4. Conclusions

Regarding the treatment of physical properties of pure ILs, temperature dependence of physical properties such as density, speed of sound, or refractive index is described using simple polynomial expressions, mainly equations of first, second and third order. For the viscosity, usually the VFT or mVFT equations are strongly recommended. For binary mixtures containing ILs, the dependence of the physical properties with temperature and composition is also correlated using empirical equations, and their excess properties are generally fitted to a Redlich-Kister type equation. For the fitting of the excess properties of ternary systems containing ILs, the most widely used equation is that proposed by Cibulka.

For the correlation of experimental data concerning phase equilibria of mixtures containing ILs, several g^E-based models have been applied in literature (Pitzer, Extended Pitzer model of Archer, UNIQUAC, NRTL, e-NRTL, NRF, MNRTL, MSA-NRTL, EW), being the NRTL model the one that unifies simplicity and satisfactory results for the treatment of vapor-liquid, liquid-liquid and solid-liquid equilibria.

The use of EoS for the modeling of phase equilibria involving ILs is frequent. Unlike what happens with g^E models, most of the literature with EoS involve VLE data, rather than LLE or SLE. Nevertheless, many authors have proved that excellent results can be obtained in data correlation or, for some cases, even prediction. The main difficulty with the application of EoS is the calculation of pure component parameters for ILs. Up to date, a general procedure has not yet been defined. Very recently, a complete review on the use of EoS with ILs, with special emphasis on the obtention of model parameters, has been published by Maia and co-workers [137]. The interested reader is directed to that work for further details.