Solution Structure of Isoactivity Equations for ... - ACS Publications

Subscriber access provided by University of Massachusetts Amherst Libraries

Article

Solution Structure of Isoactivity Equations for Liquid–liquid Equilibrium Calculations Using the Nonrandom Two-Liquid Model Zheng LI, Kathryn A Mumford, Kathryn H. Smith, Jian Chen, Yong Wang, and Geoffrey W Stevens Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.5b04469 • Publication Date (Web): 18 Feb 2016 Downloaded from http://pubs.acs.org on February 20, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Industrial & Engineering Chemistry Research is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 14

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Solution Structure of Isoactivity Equations for Liquid–liquid Equilibrium Calculations Using the Nonrandom Two-Liquid Model Zheng Li,† Kathryn A. Mumford,† Kathryn H. Smith,† Jian Chen,§ Yong Wang,† Geoffrey W. Stevens*,† †Particulate Fluids Processing Centre, Department of Chemical and Biomolecular Engineering, The University of Melbourne, Victoria, 3010, Australia §State Key Laboratory of Chemical Engineering, Tsinghua University, Beijing, 10084, PR China Abstract: Calculation of liquid–liquid equilibrium (LLE) based on an activity coefficient model is a common and significant problem in chemical thermodynamics. This is usually carried out either by minimizing the Gibbs free energy of the system coupled with stability test (tangent plane distance criterion) or solving the isoactivity equations. While established, the stability test requires a very robust algorithm. By contrast, it is easy to solve the isoactivity equations, nevertheless, the solution strongly depends on the initial estimate. This study aims to investigate the structure of the solutions of isoactivity equations for LLE systems in order to provide a calculation strategy that is not so dependent on initial estimates while still maintaining its simplicity. A systematic study covering onephase and two-phase ternary, quaternary and quinary LLE systems and three-phase ternary systems has been carried out using the nonrandom two-liquid (NRTL) model and it has been found that the equation-solving approach can be applied to multiphase systems.

1. Introduction Calculation of liquid–liquid equilibrium (LLE) using an activity coefficient model is an important and common problem in chemical engineering and it provides rational design for separation processes such as liquid–liquid extraction. In general there are two approaches for calculation of LLE1-4: (1) minimization of the Gibbs free energy and (2) solving isoactivity equations (the K-value method). The minimum of the Gibbs free energy is a necessary and sufficient condition for stable phase equilibrium, however, the global minimum is often difficult to obtain for complex systems due to mathematical difficulty caused by the highly non-linear form of the objective function, additionally incorrect phase distributions may arise if local minima are obtained1-3, 5. As a result, it is essential to have a stability test to verify whether a result corresponds to the global minimum, which is usually carried out by the tangent plane distance (TPD) criterion6-8. Consequently minimizing the Gibbs free energy coupled with the stability test has become a popular procedure for calculating phase equilibrium and the procedure proceeds in a stepwise manner with composition estimates generated by stability test until a stable equilibrium state is found9-11. The stability test requires that F(x) (the distance between the Gibbs free energy surface and the tangent plane at composition z as a function of composition x) be non-negative for all composition x. While it is too time consuming to exhaustively test all the compositions, it is sufficient to test all the stationary points of F(x). Hence, the problem converts to locating all the stationary points (minima, maxima and saddle points) of F(x). A major challenge of this task is lack of certainty that all stationary points have been found and it requires a very robust algorithm for reliable results. The development of algorithms for stability test has attracted extensive investigations10, 12-16. The second approach, solving the isoactivity equations, is numerically easier than minimizing the Gibbs free energy, however, the results strongly depend on the initial estimate and poor estimates may lead to erroneous solutions that may correspond to local minima, maxima or saddle points of the Gibbs free energy1-3. Moreover, one estimate only produces one solution and there is no method to provide successive estimate as occurs with the stability test3. Consequently the equation-solving approach has been largely abandoned in favour of the first approach. Nevertheless, the simplicity and lower computational effort makes the equations-solving approach attractive. It would be a useful method if the solution obtained did not depend on the initial estimate. Towards this end, this study explores the solution structure of isoactivity equations for LLE systems including one-phase and twophase ternary, quaternary and quinary systems and three-phase ternary systems and makes use of this 1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 14

information to provide strategies for LLE calculations based on solving isoactivity equations using the nonrandom two-liquid (NRTL) model.

2. Solution structure of isoactivity equations 2.1 Isoactivity equations At equilibrium, the isoactivity equation of a LLE system with π phases and N components can be written as: x iφ γ iφ = x iφ + 1 γ iφ + 1 (i=1, 2, …, N; φ=1, 2, …, π-1)

(1)

where i is the component in the system, φ represents the number of phase, x is mole fraction of a component and γ is activity coefficient. The isoactivity equations must be constrained by mass balance. For a system with N components and π phases at given pressure and temperature, we can either specify mole fractions of (N-π) components in one phase or specify the total quantity of each component. The situation under the first type constraint has been discussed17-19, this paper investigates the second situation, i.e. with total amount of each component specified. The mass balance can be expressed as: π

∑ n iφ = n i

(φ=1, 2, …, N)

(2)

φ =1

where niφ represents moles of component i in phase φ and ni is the total moles of component i in the system. The thermodynamic model used in this study is the NRTL model20 although the same procedure could be applied to other thermodynamic models, the core NRTL equation is: N

ln γ i =

∑τ jiG ji x j j =1 N

∑ Gki xk k =1

N

N

+∑ j =1

x j Gij N

∑ Gkj xk

(τ ij −

k =1

∑ xlτ ljGlj l =1 N

∑ Gkj xk

)

(3)

k =1

with

τ ji = ( g ji − g ii ) / RT

(4)

G ji = exp( −α jiτ ji)

(5)

where gji is energy of interaction between a j-i pair of molecules, ߬ji (߬ji≠߬ij) and αji (αji=αij) are the energy parameter and non-randomness parameter respectively. In the following sections, the solution structure of the isoactivity equations using the NRTL model will be described. 2.2 Two-phase LLE systems 2.2.1 Two-phase ternary systems For a two-phase ternary system at given pressure and temperature, which perhaps is the most common case, the isoactivity equations (equation (1)) and the mass balance constraints (equation (2)) can be written as follows: x1I γ1I = x1II γ1II

(6a)

x 2I γ 2I = x 2II γ 2II

(6b) 2

ACS Paragon Plus Environment

Page 3 of 14

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


x 3I γ 3I = x 3II γ 3II

(6c)

n1I + n1II = n1

(7a)

n 2I + n 2II = n 2

(7b)

n3I + n3II = n3

(7c)

There are six equations with six variables (n1I, n2I, n3I, n1II, n2II, n3II) describing the equilibrium of a ternary LLE system. Activity coefficients as given in Eq.(3) are functions of mole fractions, which are functions of mole numbers. Inserting Eq. (3) into Eqs. (6a)-(6c) leads to complex non-linear equations that are difficult to obtain analytical solutions, hence we need to explore the numerical solutions as we usually do for complex functions. Three variables (n1II, n2II, n3II) can be eliminated by rearrangement of Eqs. (7a)-(7c) and inserting into Eqs. (6a)-(6c), resulting in three isoactivity equations with three variables. For convenience of further illustration, we use the system “Carbon tetrachloride (1) + 2propanol (2) + Water (3)” from Sørensen and Arlt1 (Page 10) as an example and specify the total amount of each component as: n1= 0.56972 (0.56812+0.00160); n2= 0.41035 (0.31166+0.09869); n3=1.01992 (0.12022+0.89970). This composition is an experimental two-phase composition and is listed as E in tables 1 and 2. Each of the three isoactivity equations represents a surface, made of points corresponding to values of variable sets (n1I, n2I and n3I) meeting the equation, in a three dimensional space with n1I, n2I and n3I as coordinates. The three surfaces are presented in Figure 1 using the mathematical software “Mathematica” (version 9.0.1) with equations (6a), (6b) and (6c) displayed in red, green and blue respectively.

Figure 1. Surfaces of solutions of isoactivity equations for the ternary system “Carbon tetrachloride (1) + 2propanol (2) + Water (3)” from Sørensen and Arlt1 (Page 10). (a), (b) and (c) are the solutions of Eqs. (6a), (6b) and (6c) respectively under constraint of Eqs. (7a)-(7c).

3 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 14

Each of the two surfaces may intersect with each other and the intersection lines are presented in Figures 2(a)-2(c). It can be seen that each pair of surfaces intersect on a straight line and two curves. Data points on the overlapping region of all the three groups of intersection lines satisfy all the three isoactivity equations. All the intersection lines are presented together in Figure 2(d). It is shown that the three straight lines exactly overlap, thus all the data points on the line are the solutions of the isoactivity equations. In fact, the straight line represents the symmetric solutions, i.e. niI = niII (i=1, 2, 3). These are obviously trivial solutions since they represent a single phase system. In addition, the intersection lines have two intersection points, S.1 and S.2. Compositions of the two points together with the experimental composition E are presented in Table 1. Interestingly S.1 and S.2 actually are the same with composition in two phases reversed. This solution (either S.1 or S.2) is named “exact solution”. In short, there is only one exact solution, which is the correct solution for stable phase equilibrium, and a number of symmetric solutions for this particular example.

Figure 2. Intersection lines (points) of the isoactivity equations for the ternary system “Carbon tetrachloride (1) + 2-propanol (2) + Water (3)” from Sørensen and Arlt1 (Page 10). Intersection lines of red and green surfaces, red and blue surfaces, green and blue surfaces are shown in (a), (b), and (c) respectively; intersection lines are presented together in (d), black and magenta points in (d) are solutions E, S.1 and S.2 in Table 1, the straight line in the centre represents symmetric solutions.

Figure 2(d) shows the dependence of the solution of isoactivity equations on the initial estimate. If the initial estimate is close to either S.1 or S.2, the correct solution would be found; otherwise, a symmetric solution is likely to be obtained. To investigate the effect of initial estimate, we solve the isoactivity equations multiple times (1000 in this example) with random estimations in the range of components’ total amounts. The equations are solved by the “fsolve” command in “MATLAB” (version R2011b). The command “fsolve” has three algorithms: trust-region-dogleg (default), trustregion-reflective and levenberg-marquardt21. 1000 solutions obtained using this method are presented in Figure 3. Solutions that are out of the range 0~ni ( i=1, 2 and 3) or with large residues (residue as 4 ACS Paragon Plus Environment

Page 5 of 14

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


defined in Eq.(8)) are excluded and assigned to be -0.1. 419 accurate solutions were found and a subset of the solutions is provided in Table 1. σ=

value of the left side of an equation −1 va lue of the right side of an equation

(8)

Table 1. A subset of solutions for isoactivity equations in moles (ternary system) Solutions E

n1 I 0.56812

n2 I 0.31166

n3I 0.12022

n1II 0.00160

n2II 0.09869

n3II 0.89970

G/RT -

S.1 S.2

0.55064 0.019078

0.31389 0.096462

0.12149 0.89844

0.019078 0.55064

0.096462 0.31389

0.89844 0.12149

-0.56729 -0.56729

S.3 S.4

0.19321 0.024096

0.13916 0.017356

0.34589 0.043137

0.37651 0.54562

0.27119 0.39299

0.67403 0.97678

-0.40987 -0.40987

S.5 0.067562 0.048663 0.12095 0.50216 Note E represents experimental two-phase compositions.

0.36169

0.89897

-0.40987

Figure 3. 1000 solutions of the isoactivity equations for the ternary system “Carbon tetrachloride (1) + 2propanol (2) + Water (3)” from Sørensen and Arlt1 (Page 10). (a): n1I, n1II; (b): n2I, n2II; (c): n3I, n3II. The solutions S.1 and S.2 occurred 82 and 93 times respectively and they are pointed out in the figures. Excluded solutions are assigned to be -0.1.

It can be seen that these solutions distribute in the range of 0~ni (i=1, 2 and 3) evenly and no obvious properties can be observed. To further understand the solution properties, the Gibbs free energy values of these solutions are calculated using the following equations and results are presented in Figure 4. 5 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Gj g = ×nj RT RT

(j=I, II)

g g ID g E = + = ∑ xi ln xi + ∑ xi ln γi RT RT RT i i

Page 6 of 14

(9) (10)

3

nφ = ∑ ni φ

(i=1, 2 and 3; φ=I, II)

(11)

i =1

G GI G II = + RT RT RT

(12)

R is the gas constant and T is absolute temperature, G is the Gibbs free energy and g is the molar Gibbs free energy, the superscripts “ID” and “E” represent ideal mixture property and excess property respectively. Interestingly the distributions of GI/RT and GII/RT are similar as if one is the reflection of another and the G/RT converges into two values: -0.40987 and -0.56729. Conversion of these solutions from moles into mole fractions (Table 2) reveals that all the accurate solutions except S.1 and S.2 (Table 1) are in fact the same solution, i.e. solution S.S in Table 2. The solution S.S is a symmetric solution and it represents a homogenous single phase, this is why the Gibbs free energy values for all of these solutions have the same value. The solutions S.1 and S.2, which are the same solution with the compositions in two phases reversed, are represented as S.A in Table 2. Solutions represented by S.S (S.3 onwards in Table 1) have larger Gibbs free energy values than the solutions S.1 and S.2, indicating that they are not stable solutions. Thus S.A (S.1 or S.2) is the correct solution. Solution A occurred 175 time (82 times for S.1 and 93 times for S.2) in 1000 solutions, meaning that the possibility of getting the correct solution for any random estimate is about 17.5% in this example.

Figure 4. Gibbs free energy of the 419 solutions for the ternary system “Carbon tetrachloride (1) + 2-propanol (2) + Water (3)” from Sørensen and Arlt1 (Page 10). (a): GI/RT; (b): GII/RT; (c): G/RT. The solutions S.1 and S.2 occurred 82 and 93 times respectively and they are pointed out in the figures.


Page 7 of 14

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 2. Solutions for isoactivity equations in mole fractions (ternary system) Solutions E

x1 I 0.56812

x2 I 0.31166

x3I 0.12022

x1II 0.00160

x2II 0.09869

x3II 0.89970

S.A S.S

0.55845 0.28486

0.31834 0.20518

0.12321 0.50996

0.018815 0.28486

0.095133 0.20518

0.88605 0.50996

Note E represents experimental two-phase compositions.

2.2.2 Two-phase quaternary and quinary systems For quaternary and quinary LLE systems, we are unable to plot the solution surfaces of isoactivity equations as done for ternary systems because there are too many variables. However, the solution structure can be explored by solving the equations multiple (e.g. 1000) times, followed by analysis of the solutions. The system of “n-Hexane (1) + n-Octane (2) + Benzene (3) + Sulfolane (4)” at 298.15 K from Chen et al.22 is taken as an example of quaternary system with initial components specified as: n1=0.313 (0.282+0.031), n2=0.183 (0.171+0.012), n3=0.950 (0.524+0.426) and n4=0.554 (0.023+0.531). Only one correct physical solution that corresponds to two phase-inversion solutions was found in 1000 solutions with random initial estimates. The correct solution S.A occurred 103 times and the symmetric solution S.S occurred 531 times. The solutions in moles and in mole fractions are given in Tables S1 and S2 respectively in Supporting Information (SI). Chen et al.22 also investigated the LLE of the quinary system “n-hexane (1) + n-octane (2) + benzene (3) + toluene (4) + sulfolane (5)” at 298.15 K, which is used as a calculation example in this study. The initial components are specified as: n1=0.186 (0.154+0.032), n2=0.150 (0.130+0.020), n3=0.245 (0.133+0.112) and n4=0.828 (0.494+0.334), n5=0.590 (0.088+0.502). Only one correct physical solution corresponding to two phase-inversion solutions was found in 1000 solutions with random initial estimates. The correct solution S.A occurred 260 times and the symmetric solution S.S occurred 360 times. The solutions in moles and in mole fractions are given in Tables S3 and S4 respectively in SI. Over 50 two-phase LLE systems including ternary, quaternary and quinary systems from a wide range of data sources1, 11, 22-31 have been tested with at least two different overall compositions for each system and it was found that the exact solution was always unique for all the tested systems. In other words, the seemingly multiple solutions are in fact the symmetric solutions which can be identified easily when expressed in mole fractions. The unique “exact solution” is the correct solution for a twophase LLE system. This is a simple and distinct property of the solution structure of two-phase LLE systems. Ten more calculation examples (5 ternary systems, 3 quaternary systems and 2 quinary systems) are provided in SI. 2.3 One-phase LLE systems The same ternary system “Carbon tetrachloride (1) + 2-propanol (2) + Water (3)” from Sørensen and Arlt1 (Page 10) was investigated for solution structure with components specified as: n1=0.56812, n2=0.31166, n3=0.12022. This overall composition is known to be a one-phase system, we however suppose it is a two-phase system and analyse its solution structure by plotting the solution surfaces in a three dimensional space. The results are shown in Figure 5. It can be seen that the three straight lines in the centre exactly overlap, data points on the line represent symmetric solutions, i.e. one-phase solution. No two-phase solution was found. Examination of a number of ternary, quaternary and quinary LLE systems from a range of literature shows that one-phase systems only have one onephase solution and do not have two-phase solutions.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 14

Figure 5. Solutions of the isoactivity equations for the ternary system “Carbom tetrachloride + 2-propanol (2) + H2O (3)” from Sørensen and Arlt1 (Page 10) with total amount specified as: n1=0.56812, n2=0.31166, n3=0.12022. The straight line in the centre represents one-phase solution and no two-phase solution was found.

2.4 Three-phase LLE systems While the majority of LLE systems have two phases, there are some systems that have three phases. Equilibrium data for a few three-phase ternary1, 32-36 and quaternary37-39 LLE systems have been determined. Unfortunately, none of these systems has available NRTL parameters in literature, thus direct investigations on solution structures of their isoactivity equations is impossible. Wasylkiewicz et al.10 found that the ternary system “Toluene (1) + 1-propanol (2) + Water (3)” exhibits three phases at 5 °C using the NRTL parameters regressed on the basis of equilibrium data at 25 °C from Sørensen and Arlt1 (Page 580). Despite the lack of experimental data supporing the existance of three phase region under this condition, Wasylkiewicz et al.10 found it a good example to perform their algorithm. Using the same system and the same NRTL parameters, we investigated the solution structure of isoactivity equations with total components specified as: n1= 0.50, n2= 0.30, n3=0.20. Firstly, we supposed this composition has two phases and explored its solution structure. Shown in Figure 6, the solution structure reveals three exact solutions (S.A, S.B and S.C in Table 3) plus phase-inversion solution for each of them, in addition to symmetric solutions (S.S in Table 3) as represented by the straight line in the centre. S.A, S.B and S.C occurred 206, 135 and 2 times respectively in 1000 solutions with random estimates. Secondly, we treated the system as a three-phase system and found one three-phase solution (S.D in Table 3). Comparison of the Gibbs free energy values indicates that the three-phase solution S.D is the stable phase equilibrium. Examination of three more three-phase compositions using the same system shows that three-phase systems always have more than two twophase solutions and one three-phase solution.

Figure 6. Solutions of the isoactivity equations for the ternary system “Toluene + 1-propanol (2) + Water (3)” from Sørensen and Arlt1 (Page 580) at 5 °C with total amount specified as: n1= 0.50, n2= 0.30, n3=0.20. Three exact solutions S.A, S.B and S.C are revealed.


Page 9 of 14

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 3. Two- and three-phase solutions in moles for a three-phase system Solutions n1I n2 I n3 I n1II n2II n3II S.A 0.37450 0.09555 0.02300 0.12550 0.20445 0.17700

n1III -

n2III -

n3III -

G/RT -0.38776

-

-

-

-0.38543 -0.38976

S.B S.C

0.48680 0.24294 0.08973 0.01320 0.05706 0.11027 0.49986 0.29600 0.11286 0.00014 0.00400 0.08714

S.D S.S

0.19949 0.19103 0.10138 0.00011 0.00323 0.07039 0.30040 0.10574 0.02823 -0.39004 Any point on the straight line in Figure 6 -0.37708

Note that the three-phase solution is related to the two-phase solutions: compositions of S.D in phase I, II and III is similar to phase II composition of S.A, phase II composition of S.C and phase I composition of S.A respectively. It was found that initial estimates based on combinations of compositions in different phases of two-phase solutions are significantly more efficient than random estimates for three-phase solutions. In this example, there are six compositions (two compositions in each two-phase solution), we thus have C(6,3) (C(6,3)=20) combinations of three-phase estimates. In this particular case, the estimate shown in eq. (13) produced solution S.D.    0.12550, 0.20445, 0.17700 0.00014, 0.00400, 0.08714 0.37450, 0.09555, 0.02300  (13) 4244444 3 14444 4244444 3 14444 4244444 3  14444 Phase II composition of S .C Phase I composition of S . A   Phase II composition of S . A

S.D in Table 3 has three compositions and represents three two-phase solutions: niI = niII , niI = niIII, and niII = niIII (i=1, 2, 3), on possibly one three-phase solution. This shows the relationship between two-phase solutions and a three-phase solution. A three-phase solution ideally would correspond to three two-phase solutions, each of them consists of two of the three compositions of the three-phase solution. However, there is a possibility that the six compositions of the three two-phase solutions may not be the same as (or close to) the three compositions of the three-phase solution. Furthermore, any three out of six or more compositions, which correspond to three or more two-phase solutions, may not be able to form a three-phase solution. On the other hand, one two-phase solution is insufficient to form a three-phase solution, therefore a system must have at least two two-phase solutions in order to have a three-phase solution. In short, the minimum number of two-phase solutions for a system to have a three-phase solution is two and the maximum number is not certain.

3. Strategies for LLE calculations 3.1 Summary of solution structure Solution structures of isoactivity equations for one-phase and two-phase ternary, quaternary and quinary systems and a three-phase ternary system have been systematically investigated. The main properties of solution structure can be summarized as follows: (i) Single-phase systems do not have two-phase solutions; (ii) Two-phase systems have a unique two-phase solution; (iii) Three-phase systems have multiple two-phase solutions and unique three-phase solution. All the tested systems satisfy the above solution structure properties, except for one example, which is the ternary system “Toluene (1) + 1-propanol (2) + Water (3)” at 5 °C using the NRTL parameters obtained from data at 25 °C from Sørensen and Arlt1 (Page 580) with compositions specified as: n1= 0.20, n2= 0.30, n3=0.50. In the study of Wasylkiewicz et al.10, one unstable and one stable two-phase solution were found. As shown in Figure 7, this system actually has three two-phase solutions, but no three-phase solution can be found. However, the extraplation of parameters from 25 °C to 5 °C is likely to have problems in thermodynamic consistency and there is no experimental data supporting the extrapolated parameters. Albeit whether an experimental three-phase equilibrium exists or not should not affect the phase equilinrium calculations, the thermodynamic consistency may influence 9 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 14

the solution structure considering that the NRTL model used in this study is derived based on thermodynamic principles. Nevertheless, more examinations of the solution structure of three-phase systems are warranted.

Figure 7. Solutions of the isoactivity equations for ternary system “Toluene + 1-isopropanol (2) + Water (3)” from Sørensen and Arlt1 (Page 580) at 5 °C with total amount specified as: n1= 0.20, n2= 0.30, n3=0.50. Three exact solutions S.A, S.B and S.C are revealed.

3.2 A general strategy for LLE calculations On the basis of the solution structures discussed above, a procedure for LLE calculations based on solving the isoactivity equations can be developed. For systems with unknown phase number, a general procedure is proposed as: (1). Solve two-phase isoactivity equations under mass balance constraints multiple times (e.g. 1000) using random estimates; (2). Exclude symmetric solutions; (3). (i) No solution found: single-phase system, end. (ii) One solution found: two-phase system, end. (iii) Multiple solutions found: three-phase system, go to (4). (4). Solve three-phase isoactivity equations under mass balance constraints multiple times with estimates generated based on two-phase solutions; (5). Exclude symmetric and two-phase solutions; (6) (i) The solution found, end. (ii) No solution found, repeat (4). The first step performs an exhaustive calculation to ensure all the “exact solutions” can be found. Step (3) determines the number of phases as well as the compositions for stable phase equilibrium if it is a one-phase or two-phase system, and provides rational estimates for solving three-phase solution if it is a three-phase system. If a system is known as a two-phase system, which is the most common case, only the first three steps of the procedure are required. In this case, solving the equations 20 times is recommended in step (1) of the procedure. If the correct solution does not occur in 20 solutions, the procedure can be repeated. The possibility of occurrence of the “exact solution” largely depends on the location of the solution in the space established by the amount of components: the possibility is high if the solution locates in the central area, while it is low if the solution locates at the edge of the space. Generally a solution can be found within five loops (in total solving the equations 100 times) and the calculation time is within 2 10 ACS Paragon Plus Environment

Page 11 of 14

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


seconds on an ordinary personal computer. The performance of the proposed calculation procedure may be improved if an initial estimate generating method that is more efficient than random estimates can be developed, particularly for the cases where the solution locates at the edge of the space. However, this requires investigations of a large number of LLE systems and comparing a number of possible strategies. Although interesting, it is beyond the scope of this paper that focuses on the solution structure of isoactivity equations. The proposed new calculation procedure is essentially different from minimizing the Gibbs free energy combined with the stability test. A brief comparison of these two approaches is listed in Table 8. Satisfaction of isoactivity equations is the necessary but not sufficient condition for stable phase equilibrium, however, when coupled with the solution structure analysis in all cases examined the correct solution has been obtained. The main advantage of the proposed procedure compared with the stability test approach is the ease of calculation because it does not involve minimizing challenging nonlinear functions and thus does not require robust numerical methods. Table 4. Comparison of the two procedures for LLE calculations Minimizing the Gibbs free energy coupled with TPD

Looks for solution in a stepwise way and repeats the procedure until the correct solution is found

Equation-solving coupled with solution structure Based on isoactivity equations and their solution structure Produces many candidate solutions and identify the correct one based on solution structure

Requires robust numerical methods

Simple numerical methods needed

Based on the global minimum of the Gibbs free energy

4. Conclusions Solution structures of isoactivity equations for one-phase and two-phase ternary, quaternary and quinary LLE systems and three-phase ternary LLE systems have been systematically investigated and distinct properties of solutions have been revealed. Based on this analysis, a simple calculation strategy using equations-solving approach, which mainly involves solving isoactivity equations multiple times (usually within 100 times), has been proposed for two-phase LLE systems. For a LLE system with unknown number of phases, a calculation strategy has also been proposed which starts from treating the system as a two-phase system. Additionally, two-phase solutions can be initial estimates for three-phase solution if the system indeed has three phases. This has provided more certainty to solving the isoactivity equations rather than relying on initial estimates and the proposed strategy for calculating LLE is essentially different from the stability test using TPD and has the advantage of not requiring robust numerical methods. Supporting Information Tables S1 and S2 listing solutions in moles and mole fractions respectively of the quaternary LLE system as discussed in section 2.2.2. Tables S3 and S4 listing solutions in moles and mole fractions respectively of the quinary LLE system as discussed in section 2.2.2. Ten more LLE calculation examples including five ternary systems, three quaternary systems and two quinary systems. This material is available free of charge via the Internet at http://pubs.acs.org.

Author Information Corresponding Author *E-mail: [email protected]. Fax: +61 3 83448824. Notes The authors declare no competing financial interest.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 14

Acknowledgement Support from the Australian Research Council and the Particulate Fluids Processing Centre is acknowledged. Prof. Fernando García-Sánchez at Instituto Mexicano del Petróleo is thanked for helpful discussions and for providing calculation examples using the stability test approach. Zheng Li is grateful for financial support from the China Scholarship Council.

References (1) Sørensen, J. M.; Arlt, W., Liquid-liquid equilibrium data collection-ternary systems. DECHEMA: Frankfurt, 1980. (2) Sørensen, J. M.; Magnussen, T.; Rasmussen, P.; Fredenslund, A., Liquid—liquid equilibrium data: Their retrieval, correlation and prediction Part II: Correlation. Fluid Phase Equilib. 1979, 3, (1), 47-82. (3) Nichita, D. V.; Gomez, S.; Luna, E., Multiphase equilibria calculation by direct minimization of Gibbs free energy with a global optimization method. Comput. Chem. Eng. 2002, 26, (12), 1703-1724. (4) Iglesias-Silva, G. A.; Bonilla-Petriciolet, A.; Eubank, P. T.; Holste, J. C.; Hall, K. R., An algebraic method that includes Gibbs minimization for performing phase equilibrium calculations for any number of components or phases. Fluid Phase Equilib. 2003, 210, (2), 229-245. (5) Gautam, R.; Seider, W. D., Computation of phase and chemical equilibrium: Part I. Local and constrained minima in Gibbs free energy. AIChE J. 1979, 25, (6), 991-999. (6) Baker, L. E.; Pierce, A. C.; Luks, K. D., Gibbs Energy Analysis of Phase Equilibria. Society of Petroleum Engineers Journal 1982, (October), 731-742. (7) Michelsen, M. L., The isothermal flash problem. Part I. Stability. Fluid Phase Equilib. 1982, 9, (1), 1-19. (8) Michelsen, M. L., The isothermal flash problem. Part II. Phase-split calculation. Fluid Phase Equilib. 1982, 9, (1), 21-40. (9) Nghiem, L. X.; Li, Y.-K., Computation of multiphase equilibrium phenomena with an equation of state. Fluid Phase Equilib. 1984, 17, (1), 77-95. (10) Wasylkiewicz, S. K.; Sridhar, L. N.; Doherty, M. F.; Malone, M. F., Global Stability Analysis and Calculation of Liquid−Liquid Equilibrium in Multicomponent Mixtures. Ind. Eng. Chem. Res. 1996, 35, (4), 1395-1408. (11) García-Flores, B. E.; Águila-Hernández, J.; García-Sánchez, F.; Aquino-Olivos, M. A., (Liquid–liquid) equilibria for ternary and quaternary systems of representative compounds of gasoline + methanol at 293.15K: Experimental data and correlation. Fluid Phase Equilib. 2013, 348, (0), 60-69. (12) Bonilla-Petriciolet, A.; Vázquez-Román, R.; Iglesias-Silva, G. A.; Hall, K. R., Performance of Stochastic Global Optimization Methods in the Calculation of Phase Stability Analyses for Nonreactive and Reactive Mixtures. Ind. Eng. Chem. Res. 2006, 45, (13), 4764-4772. (13) Nichita, D. V.; Gomez, S.; Luna, E., Phase stability analysis with cubic equations of state by using a global optimization method. Fluid Phase Equilib. 2002, 194–197, 411-437. (14) Rangaiah, G. P., Evaluation of genetic algorithms and simulated annealing for phase equilibrium and stability problems. Fluid Phase Equilib. 2001, 187–188, 83-109. (15) Hua, J. Z.; Brennecke, J. F.; Stadtherr, M. A., Enhanced Interval Analysis for Phase Stability: Cubic Equation of State Models. Ind. Eng. Chem. Res. 1998, 37, (4), 1519-1527. (16) Sun, A. C.; Seider, W. D., Homotopy-continuation method for stability analysis in the global minimization of the Gibbs free energy. Fluid Phase Equilib. 1995, 103, (2), 213-249. (17) Li, Z.; Mumford, K. A.; Shang, Y.; Smith, K. H.; Chen, J.; Wang, Y.; Stevens, G. W., Analysis of the Nonrandom Two-Liquid Model for Prediction of Liquid–liquid Equilibria. J. Chem. Eng. Data 2014, 59, (8), 2485-2489. (18) Marcilla, A.; Olaya, M. d. M.; Reyes-Labarta, J. A.; Wisniak, J., Comments on “Analysis of the Nonrandom Two-Liquid Model for Prediction of Liquid–Liquid Equilibria”. J. Chem. Eng. Data 2015, 60, (5), 1526-1529. (19) Li, Z.; Mumford, K. A.; Smith, K. H.; Chen, J.; Wang, Y.; Stevens, G. W., Reply to “Comments on ‘Analysis of the Nonrandom Two-Liquid Model for Prediction of Liquid–Liquid Equilibria’”. J. Chem. Eng. Data 2015, 60, (5), 1530-1531. (20) Renon, H.; Prausnitz, J. M., Local compositions in thermodynamic excess functions for liquid mixtures. AIChE J. 1968, 14, (1), 135-144.


Page 13 of 14

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(21) MathWorks. Documentation, fsolve. http://au.mathworks.com/help/optim/ug/fsolve.html;jsessionid=64086d01bf565d5d43faa4239df8, accessed on 4th Sep. 2015 (22) Chen, J.; Duan, L.-P.; Mi, J.-G.; Fei, W.-Y.; Li, Z.-C., Liquid–liquid equilibria of multi-component systems including n-hexane, n-octane, benzene, toluene, xylene and sulfolane at 298.15 K and atmospheric pressure. Fluid Phase Equilib. 2000, 173, (1), 109-119. (23) Sørensen, J. M.; Arlt, W., Liquid-liquid equilibrium data collection-ternary and quaternary systems. DECHEMA: Frankfurt, 1980. (24) Alvarez Gonzalez, J. R.; Macedo, E. A.; Soares, M. E.; Medina, A. G., Liquid-liquid equilibria for ternary systems of water-phenol and solvents: data and representation with models. Fluid Phase Equilib. 1986, 26, (3), 289-302. (25) Colombo, A.; Battilana, P.; Ragaini, V.; Bianchi, C. L.; Carvoli, G., Liquid−Liquid Equilibria of the Ternary Systems Water + Acetic Acid + Ethyl Acetate and Water + Acetic Acid + Isophorone (3,5,5-Trimethyl2-cyclohexen-1-one). J. Chem. Eng. Data 1999, 44, (1), 35-39. (26) Nakayama, T.; Sagara, H.; Arai, K.; Saito, S., High pressure liquid-liquid equilibria for the system of water, ethanol and 1,1-difluoroethane at 323.2 K. Fluid Phase Equilib. 1987, 38, (1–2), 109-127. (27) Zhang, H.; Gong, Y.; Li, C.; Zhang, L.; Zhu, C., Liquid−Liquid Equilibria for the Quaternary System Water (1) + Acrylic Acid (2) + Acetic Acid (3) + Cyclohexane (4) at (293.15, 303.15, and 313.15) K. J. Chem. Eng. Data 2011, 56, (5), 2332-2336. (28) Arce, A.; Blanco, M.; Soto, A., Determination and correlation of liquid–liquid equilibrium data for the quaternary system 1-octanol+2-methoxy-2-methylbutane+water+methanol at 25°C. Fluid Phase Equilib. 1999, 158–160, (0), 949-960. (29) Chen, J.; Mi, J.; Fei, W.; Li, Z., Liquid-liquid equilibria of quaternary and quinary systems including sulfolane at 298.15 K. J. Chem. Eng. Data 2001, 46, (1), 169-171. (30). Kumar, U. K. A.; Mohan, R., Quinary and Eight-Component Liquid–Liquid Equilibria of Mixtures of Alkanes, Aromatics, and Solvent (Furfural). J. Chem. Eng. Data 2013, 58, (8), 2194-2201. (31) Salem, A. B. S. H.; Hamad, E. Z.; Al-Naafa, M. A., Quaternary Liquid-Liquid Equilibrium of n-HeptaneToluene-o-Xylene-Propylene Carbonate. Ind. Eng. Chem. Res. 1994, 33, (3), 689-692. (32) Francis, A. W., Ternary Systems with Three Separate Binodal Curves. The Journal of Physical Chemistry 1956, 60, (1), 20-27. (33) Hou, H. Y.; Liu, S. T.; Geng, X. P., Three-liquid-phase equilibria of ternary system water-cyclohexanediethylene glycol monobutyl ether. Gaodeng Xuexiao Huaxue Xuebao/Chemical Journal of Chinese Universities 2008, 29, (11), 2249-2253. (34) Sinegubova, S. I., Critical phenomena and three-liquid-phase equilibrium in the ternary system nitromethane - water - n-heptane. High Temp. - High Pressures 1999, 31, (4), 363-366. (35) Negahban, S.; Willhite, G. P.; Walas, S. M.; Michnick, M. J., Three-liquid-phase equilibria of ternary and quaternary mixtures, water/n-decane/2-butyloxyethanol and water/n-octane/1-propanol/sodium chloride— experimental measurements and their correlation with the UNIQUAC model. Fluid Phase Equilib. 1986, 32, (1), 49-61. (36) Chakravarty, T.; Sivasubramanian, M. S.; Babu, S. V., Representation and calculation of three liquid-phase equilibria of two ternary systems. Fluid Phase Equilib. 1985, 19, (3), 221-231. (37) Xiang, Y.; An, X.; Shen, W., Three-liquid-phase equilibria for the quaternary system water + tert-butanol + n-decane + n-undecane in the tricritical region. J. Chem. Eng. Data 2007, 52, (3), 1113-1117. (38) Hou, H.; An, X.; Shen, W., Three-liquid-phases equilibria for the quaternary system water + formamide + cyclohexane + diethylene glycol monobutyl ether near the tricritical point. J. Chem. Eng. Data 2005, 50, (4), 1308-1312. (39) Bocko, P., Equilibrium of three liquid phases in water-acetonitrile-benzene-n-hexane mixtures at 65°C. Fluid Phase Equilib. 1980, 4, (1–2), 137-150.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 14

Table of Content Graphics


Solution Structure of Isoactivity Equations for ... - ACS Publications

Recommend Documents