Article Cite This: Ind. Eng. Chem. Res. 2019, 58, 10064−10072
pubs.acs.org/IECR
Prediction and Reasoning for the Occurrence of Lower Critical Solution Temperature in Aqueous Solution of Ionic Liquids Ya-Ruei Tsai and Shiang-Tai Lin* Department of Chemical Engineering, National Taiwan University, Taipei 10617, Taiwan
Downloaded via UNIV OF SOUTHERN INDIANA on July 27, 2019 at 08:48:14 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.
S Supporting Information *
ABSTRACT: In this work, the combined COSMO-SAC and Pitzer−Debye−Hückel model is used for the prediction of the liquid−liquid equilibrium for 252 binary ionic liquid solutions, including 234 LCST (lower critical solution temperature)negative mixtures and 18 LCST-positive mixtures. It is found that the combined model provides qualitative description of the mutual solubility of LCST-negative mixtures (AAD-x1 = 0.2451 and AAD-x2 = 0.04635) and LCST-positive mixtures (AAD-x1 = 0.3522 and AAD-x2 = 0.009418) when compared to experimental data. The model can also be used for the prediction of the LCST of such aqueous solutions. Though the ARD% of predicted lower critical solution temperatures is about 40%, the model usually provides correct relative LCST of different IL solutions. We found that the temperature derivative of the Gibbs free energy of mixing at the spinodal points provides very useful information regarding the occurrence of the critical solution behaviors. The fact that the model does not require any species-dependent parameters makes it a useful tool for screening and ranking of ILs in applications such as extraction and forward osmosis of seawater.
1. INTRODUCTION
upper-critical-solution temperature (UCST, Figure 1b), an hourglass shape LLE with both LCST and UCST (Figure 1c), and a closed-loop LLE with both LCST and UCST (Figure 1d). The occurrence of UCST is generally observed at high temperature as the effect of entropy of mixing is amplified by temperature. On the other hand, the occurrence of LCST behavior is believed to be an entropy-driven phenomenon4 and associated with the existence of hydrogen bonding interaction and the balance between hydrophilic and hydrophobic character.3,5 The closed-loop LLE can be considered the LCST behavior following common UCST behavior at higher temperature. Examples includes the aqueous solution of tetrahydrofuran (THF)6 or poly(ethylene oxide) (PEO).4 The hourglass-type LLE has been observed in some polymer
An ionic liquid (IL), a subclass of molten salts, refers to a class of ionic compounds that melt below 100 °C.1 Typical ionic liquids are composed of an organic cation and a suitably combined organic or inorganic anion to allow them to stay as a liquid near room temperature. Because of its high thermal and chemical stability, low volatility, and most importantly, the high tailorability by adjusting the combination of cations and anions, ILs are considered as promising candidates that may replace the use of conventional organic solvents.2 The mixture of ILs with other solvents may exhibit interesting phase behaviors because of their tunable, amphiphilic nature. For example, ILs may be completely miscible, partially miscible, or nearly immiscible with a solvent by changing the temperature.3 Typically, the phase behavior of such thermoresponsive materials can be categorized into four types, (1) a U-shaped LLE with a lower-critical-solution temperature (LCST, Figure 1a), (2) a bell-shaped LLE with an © 2019 American Chemical Society
Received: May 9, 2019 Accepted: May 16, 2019 Published: May 16, 2019 10064
DOI: 10.1021/acs.iecr.9b02551 Ind. Eng. Chem. Res. 2019, 58, 10064−10072
Article
Industrial & Engineering Chemistry Research
In a study of LCST of polymer mixtures, Melagraki et al.21 adopted the quantitative structure property relationship (QSPR) model to correlate LCSTs of polymers in organic solvents. However, perhaps due to the need of a large set of experimental data, no such model is available for IL mixtures. Zhao et al.22 performed molecular dynamics (MD) simulations on ionic liquids with amino-acid−based anions that showed LCST with water. They found that it is the competition of the hydrogen-bonding interaction between COO−, COOH, NH2, and water at different temperatures that caused the LCST type of phase separation, which agrees with the previous speculation that LCST is induced by the competition of different types of molecular interaction. However, in their work, the mutual solubility of ILs and water and the LCST are not quantified using the same force field. While MD simulations provide useful insights to the occurrence of the critical behaviors, the high computational demand makes it unsuitable for the screening of candidate IL solvents. In this work, we apply the predictive COSMO-SAC model, combined with the Pitzer−Debye−Hückel (PDH) model for long-range electrostatic interactions,20,23,24 for the prediction of LLE and LCST of IL aqueous solutions. This method requires input of only the molecular structures and the corresponding electron density distribution from quantumchemical calculations. This allows one to perform a priori prediction for the key thermodynamic properties, miscibility gaps, and lower critical solution temperatures. We show that this approach also offers a way to understand the occurrence of LCST as a competition between different molecular interactions.
Figure 1. Illustration of four typical types of LLE and the corresponding critical solution transitions.
mixtures. The miscible region or contraction in the hourglasstype LLE occurred when the dense polymer is mixed with highly expanded solvents, which results in negative entropies of mixing.7 The temperature-responsive LLE is of great interest as it can be used as a solvent for extraction of chemicals8 and in catalysis,3 sensors, and actuators.9,10 For example, Cai et al.11 proposed using the LCST in a forward osmosis process for desalination. Compared to conventional reverse osmosis (RO), this novel process replaces high-cost electricity in RO with lowcost waste heat to stimulate phase transition and can potentially be more economical for seawater desalination. The potential of customizing ILs with desired properties makes it a class of great chemicals for novel applications. At the same time, it is becoming a critical challenge to search for the right chemical structure or the right combination of cation and anion that produces the desired properties. There have been many efforts made on the modeling of LLE in ionic liquid mixtures. Paduszynski and Domanska used perturbed-chain statistical associating fluid theory (PC-SAFT) for the prediction of LLE of imidazolium and pyridinium ILs in different solvents.12 Simoni et al.13 utilized NRTL, electrolyteNRTL (eNRTL), and UNIQUAC to characterize LLE of ionic liquids. However, these models require the fitting of interaction parameters using experimental data, where a full prediction is not possible. Lei et al.14 and Paduszynski and Domanska15 extended the UNIFAC parameters to include ionic liquids. However, there are only a limited number of pair interaction parameters available, making it difficult for designing new ionic liquids. Palomar and co-workers16−18 have performed the prediction and analysis of several properties and phase behaviors of ionic liquids using the COSMO-RS model. Freire19 and Lee20 also used COSMO-RS and COSMO-SAC models, respectively, to predict LLE of IL solutions. Though they demonstrated the capability of performing priori prediction of LLE via COSMO-based models, the studied LLE were limited mostly to UCST-type behaviors. In addition, the behavior near critical solution temperature and the feasibility of extending their model into LCST mixtures remain unexamined.
2. THEORY AND MODEL 2.1. Liquid−liquid Equilibrium (LLE) Calculation. The LLE calculation is based on equality of fugacity of both phases. Since the reference states chosen for both phases are pure liquids, the equality of fugacity can be expressed as the equality of activity, as in eq 1. xiIγiI = xiIIγiII
(1)
2.2. Activity Coefficient Model: COSMO-SAC 2010 + PDH. In this work, we use the activity coefficient model COSMO-SAC 201025 combined with the PDH model for the prediction of activity coefficients, i.e. ln γi = ln γiPDH + ln γiCOSMO ‐ SAC
(2)
According to the PDH model, the activity coefficients for species i can be expressed as ln γi
PDH
ÄÅ ÉÑ Å 2 z i2 Ix − 2Ix3/2 ÑÑÑÑ 1000 ÅÅÅ 2z i ÑÑ A ⌀ÅÅÅ ln(1 + ρ Ix ) + ÅÅ ρ ms 1 + ρ Ix ÑÑÑÑ ÅÇ Ö (3)
=−
where γPDH is the activity coefficient from the PDH model, Ms i is the molecular weight of the species, zi is the charge of the charged species, Ix is the ionic strength of the solution (Ix = 1/ 2Σxiz21), and A⌀ is the Debye−Hückel model parameter ( A⌀ =
1 3
(
2πNads 1/2 2 (e /ϵskT )3/2 ), 1000
)
in which ds is the density of
solvent, Na is the Avogadro number, e is the charge of an electron, ϵs is the dielectric constant of a neutral solvent, k is the Boltzmann constant, and ρ is 10065
DOI: 10.1021/acs.iecr.9b02551 Ind. Eng. Chem. Res. 2019, 58, 10064−10072
Article
Industrial & Engineering Chemistry Research ρ=R
ΔW (σmt , σns) = c ES(σmt + σns)2 − chb(σmt , σns)(σmt − σns)2
2e 2Nads Ms ϵskT
(10)
(4)
where the electrostatic interaction parameter CES is treated as a temperature-dependent parameter, as shown in eq 11, where AES and BES in eq 11 are constant universal parameters. The value of the universal parameters can be found in Chen’s work.28 The hydrogen bonding interaction parameters, according to different interaction between surface, are defined in eq 12.
where R can be expressed as R = kra, in which a is the summation of each ion radii and kr is assumed to be 2.5 for all solvents.24,26 The COSMO-SAC 2010 model describes the nonideality with a combinatorial term and a residual term: ln γiCOSMO‐SAC = ln γicomb + ln γires
The Staverman−Guggenheim model combinatorial term, as shown in eq 6. ln γi
comb
(5) 27
C ES = AES +
is chosen for the
⌀ θ ⌀ z = ln i + qi ln i + li − i Σjxjlj xi 2 ⌀i xi
A it (σmt) Ai
(6)
The residue term is then determined from the summation of segment surface contributions, as expressed in eq 13 ln γires =
nhb,OH,OT
∑
∑ pit(σm)·[ln Γst(σm) − ln Γit(σm)] σm
t
where aeff is an effective surface area. Because of the use of asymmetric reference state in the PDH model (eq 3), the consistency of reference states of species needs to be taken care of. The reference states are chosen to be pure liquid for uncharged species and infinite dilution for charged species. The overall activity coefficient for uncharged species can be expressed as ln γi = ln γiCOSMO‐SAC +ln γiPDH
(14)
The activity coefficients of ionic species can be expressed as ln γi =ln(γiCOSMO ‐ SAC /γi∞,COSMO ‐ SAC) + ln γiPDH
(15)
γ∞,COSMO‑SAC i
(7)
where is the infinite dilution activity coefficients of ionic species. 2.3. Stability Limit. The spinodal curve indicates the boundary of thermodynamically stable phase. It can also be used to characterize the origin of different critical solution behaviors. The spinodal boundary of a binary mixture can be determined from the inflection point of the Gibbs free energy of mixing
ΣixiA i pit (σmt)
ij ∂ 2Δmix g yz jj zz = 0 jj 2 z z ∂ x k {T,P
(8)
where the subscript s represents the mixture solution. The segment activity coefficient Γtj (σm) of a segment of type t and charge density σm (j = i for pure liquid i and j = S for mixtures) can be determined using eq 9 l o o nhb,OH,OT Σσnpjs (σns)Γ sj (σns) ln Γtj(σmt) = −lnm oΣs o n ÄÅ É| ÅÅ −ΔW (σmt , σns) ÑÑÑo Å ÑÑo expÅÅÅ ÑÑ} o ÅÅÇ ÑÑÖo RT ~
Ai aeff
(13)
where the Ai is the total surface area of the species i and = ΣtAti (σtm is the total surface area of segments with the charge density σm. The σ-profile for a mixture can then be calculated as mole fraction average
ΣixiA i
(11)
(12)
Ati (σm)
pst (σm) =
T2
l o cOH − OH if s = OH, t = OH, and σmtσms < 0 o o o o o t s o o o cOT − OT if s = OT, t = OT, and σmσm < 0 chb(σmt , σms) = o m o o o cOH − OT if s = OH, t = OT, and σmtσms < 0 o o o o o 0 otherwise n
where φi is the normalized volume fraction, θi is the normalized surface-area fraction, li = z/2(ri − qi) − (ri − 1), z = 10 is the coordination number, xi is the mole fraction, ri and qi are the normalized volume and surface area parameter, where qi = Ai/q and ri = Vi/r, and Ai and Vi are the cavity surface and the cavity volume of a molecule, respectively, and both can be obtained from the COSMO calculation. The constants q = 79.53 Å2 and r = 66.69 Å3 are the surface area and the volume used for normalization. The residual term is used to describe the nonideality from molecular interactions. A σ profile, the distribution of screening charge density for a molecule when embedded in a conductor, is used to characterize the electronic nature of the molecule associated with surface level t (t = nhb, OH, or OT), where nhb stands for non-hydrogen bonding level, while the two hydrogen bonding (hb) levels, OH and OT, are introduced to differentiate the surfaces on hydroxyl group (AOH i (σ), e.g., OH on water and alcohol) and other types of hydrogen-bonding surfaces (AOH i (σ), e.g., O on ketones and NH2 in amines) pit (σm) =
BES
(16)
(
A phase is stable if (
(
∂ 2Δmix g ∂x 2
)
2
∂ Δmix g ∂x 2
)
> 0) and is unstable if
T ,P
< 0. From Figure 1, it can be seen that a weak
T ,P
criteria of LCST-positive mixture (i.e., a mixture that exhibits a LCST near the temperature under study) is that the mutual solubility (binodal curve) increases with decreasing temperature, whereas that of LCST-negative (i.e., a mixture that does not exhibit a LCST near the temperature under study and, instead, a UCST is expected) is the opposite. (Note that we use LCST-positive and LCST-negative to indicate whether a
(9)
The interaction energy between two surface segments ΔW can be calculated by considering electrostatic interaction and hydrogen-bonding interactions, as shown in eq 10 10066
DOI: 10.1021/acs.iecr.9b02551 Ind. Eng. Chem. Res. 2019, 58, 10064−10072
Article
Industrial & Engineering Chemistry Research lower-critical transition temperature is to be expected or not near the temperature where eq 16 is evaluated.) Since spinodal decomposition occurs within the mutual solubility, it is expected that the spinodal curve behave similarly to the binodal curve. Therefore, the temperature variation of the inflection point is different for LCST-positive and LCSTnegative systems. Figure 2 illustrates the Δmixg of a binary
Δmix g = RT Σxi ln xi + g ex(comb) + g ex(nhb) + g ex(OH) + g ex(OT) + g ex(PDH)
where gex(comb), gex(nhb), gex(OH), gex(OT), and gex(PDH) are contributions of excess Gibb’s free energy from the combinatorial term, non-hydrogen-bond level, OH level, and OT level segments in COSMO-SAC 2010 and from the PDH term, respectively. The first term, RTΣxi ln xi, is the entropic contribution from ideal mixing. We will show that the combination of eqs 17 and 18 provides good predictions to the occurrence of LCST and the dominating forces for its occurrence.
3. COMPUTATIONAL DETAILS 3.1. General Procedure of Phase Equilibrium Calculation. The general procedure of LLE flash calculation is as follows. At a given temperature, an initial guess of the feed composition of species i, xi, and the partition of each species in the two phases, Ki = xIIi /xIi , are required as the input. The compositions of species i in both phases in next step are solved by Adomian decomposition method29 according to the current Ki. For the mixtures considering full dissociation of the ions, the mean activity coefficient of anions and cations is applied so that the equal composition of anions and cations would be ensured through the loop. The updated Ki is obtained from the ratio of the activities of the two phases Ki = (xIIi γIIi )/(xIi γIi )), where the activity coefficients are calculated from COSMOSAC 2010 + PDH model. The equilibrium compositions are updated with the new Ki values and the procedure is repeated until Ki for every species reaches unity. The molecular structures, including those of organic solvents and ions, are generated by Avogadro,30 and the DFT calculation is conducted using Gaussian0931 with the combination of basis set and functional B3LYP/6-31G(d,p). The molecular geometry neutral ion pairs are not optimized in quantum-chemical calculations in this work. Instead, σ-profiles of ion pairs are generated through summation of σ-profiles of cations and anions. The COSMO surface area and volume of ion pairs follow the same manner. The temperature-dependent dielectric constants εs and density ds are all represented by the properties of nonionic solvents in this work. In this work, the LCST for aqueous solutions involving ionic liquids are determined by decreasing the temperature by 0.5 K and calculating the LLE until no phase separation occurs. The lowest temperature where phase separation occurs is regarded as the LCST for the IL−water mixture. The calculation of temperature derivative of the spinodal point is done by numerical derivatives. The inflection point x1S(T) at a given temperature T is determined by the second order derivative of Gibbs free energy of mixing over composition using eq 19, with Δx1 = 0.001. The composition, where the second-order derivative is zero (or closest to zero), is the inflection point at the given temperature.
Figure 2. Illustration of the temperature dependence near spinodal boundary. For LCST-positive mixtures, T2 > T1 > T0, while for LCSTnegative mixtures, T2 < T1 < T0. The spinodal and binodal composition of component 1 lean phase at temperature T1 are denoted as x1s(T1) and x1b(T1), respectively. The common tangent is used to indicate the binodal points at T0, while the dashed double arrow shows that we choose a higher temperature and a lower temperature compared with T1 at the composition of x1s(T1) for the numerical derivative. The dark colors of each isotherm indicates the unstable region.
mixture at three temperatures. There are two inflection points on each isotherm, indicating the stability limit of component 1rich and component 1-lean phases, respectively. In the case where T2 > T1 > T0, the two inflection points are closer to each other when temperature decreases. This indicates that the two phases become more miscible with decreasing temperature, and LCST may be expected. In the case where T2 < T1 < T0, the two phases become more miscible with increasing temperature, and UCST may be expected. Therefore, the temperature derivative of
(
∂ 2Δmix g ∂x 2
)
(18)
at the inflection point
T ,P
i 2 y ∂ jjjjjij ∂ Δmix g zyz zzzz z jjjjj z z ∂T jjk ∂x12 z{ zzz T , P k {x1s(T ) l o o>0 if T is close to UCST(LCST‐negative) m o o o (17) n
