Aspects To Be Considered for the Development of a Correlation

Sep 10, 2010 - The existing commercial software packages for equilibrium data correlation, or analogous tools included in chemical process simulators,...
0 downloads 10 Views 3MB Size
10100

Ind. Eng. Chem. Res. 2010, 49, 10100–10110

Aspects To Be Considered for the Development of a Correlation Algorithm for Condensed Phase Equilibrium Data of Ternary Systems A. Marcilla,* M. M. Olaya, M. D. Serrano, and J. A. Reyes-Labarta Chemical Engineering Department, UniVersity of Alicante, Apartado 99, Alicante 03080, Spain

The existing commercial software packages for equilibrium data correlation, or analogous tools included in chemical process simulators, have been widely used and their extensive capabilities are well-known. Nevertheless, and as far as LLE calculations are concerned, they have been designed for the simplest equilibrium behavior. For example, for ternary systems, only the correlation of type 1 and type 2 LLE data is possible. None of these applications allows for the correlation of type 0 LLE, type 3 LLLE, or type 4 LLSE systems, although these are interesting for many industrial purposes. To enable a possible extension of such software, a robust computation algorithm has been developed to correlate phase equilibrium data among condensed phases for all of these more complicated behaviors. Moreover, this algorithm includes some new strategies to avoid the main problems occurring through correlation procedures, which lead to convergence problems or inconsistent results. These strategies are explained in detail in the present work. 1. Introduction Liquid-liquid or solvent extraction is a common separation technique used in the chemical industry for a number of important applications in which distillation is impractical (mixtures containing sensitive materials, components with high or close to boiling temperatures, existence of azeotropes, etc.). The design of industrial extraction processes requires the correlation of experimental equilibrium data to obtain the parameters of a model that represents the phase equilibrium behavior of the system (data reduction). Significant savings in operating costs can be achieved by improving the accuracy of this data correlation, which depends on two factors: the capacity of the model used and the calculation algorithm followed. It is generally accepted that existing local composition models have some limitations in reproducing the equilibrium behavior among condensed phases for certain systems. In a previous paper,1 some limitations of the representative local composition model NRTL were analyzed and the necessary requirements for more capable models were deduced. With regard to the second factor mentioned above, there are still some important limitations in the equilibrium data correlation tools included in commercial software packages (i.e., DECHEMA Data Preparation Package,2 Chemcad3). These tools are usually designed for the simpler LL equilibrium cases (i.e., type 1 and type 2 LLE ternary systems4). For example, none of these applications allow for the equilibrium data correlation of singular type 0 (island type) LLE systems, which sometimes occur in mixtures with polymers or with the simultaneous presence of a weak acid and base.5,6 Neither can more complex systems be correlated, although these are important for certain industrial applications: type 3 LLLE, which are frequent in processes involving surfactants,7-10 or type 4 LLSE systems, used to improve separations by the salting-out effect.11-14 Many methods for equilibrium calculations have been developed in the past two decades. The existing procedures are classified in two main groups:15,16 the equation-solving approach (K-value method) and the Gibbs energy minimization method. The flash equation-solving approach is very common, and it is * To whom correspondence should be addressed. Tel.: (34) 965 903789. Fax: (34) 965 903826. E-mail: [email protected].

applied in many calculation algorithms developed for VLE:16 Cairns and Futzer, 1990;17 Seider and Widago, 1996;18 Boston and Brit, 1978.19 However, this method imposes a necessary but not sufficient condition of equilibrium and can give multiple solutions corresponding to local minima. Methods for Gibbs energy minimization are preferable because they unambiguously define the condition of equilibrium for a system. These methods have been applied through different techniques such as the tangent plane criterion,20,21 the equal area method,22 and the maximum area method.23 The first group of Gibbs energy minimization methods is mainly based on the Michelsen tangent plane stability approach.20,21 The stability analysis can be used for phase assignment, phase split initialization, and flash calculation validation.24 The original procedures do not guarantee that all stationary points have been found and sometimes fail to obtain the global minimum solution. Consequently, various improvements on the performance of the calculation methods consisting of global optimization strategies have been developed: that is, the deterministic branch and bound procedure by McDonald and Floudas;25-28 homotopy continuation methods by Sun and Seider29 and Jalali and Seider;30 interval methods by Hua et al.;31 and phase stability methods by Saber and Shaw.32 With regard to the second group, Eubank and Hall developed the equal area rule22 and applied it to binary two-phase equilibrium calculations. Shyu et al.33 applied it to ternary liquidphase systems using activity coefficient models, and Hanif et al.34 developed the procedures needed to calculate phase behavior of ternary mixtures with an EOS model. Shyu et al.35 also extended the equal area rule minimization method with the maximum partial area rule providing a necessary and sufficient condition for equilibria involving multicomponent systems. Eubank et al.23 proposed the alternative approach of the maximum area method. It ensured the Gibbs energy minimization for binary LL, VL, and VLL but failed for ternary multiphase equilibria.36,37 Some modifications have been carried out: first, the method was extended to a volume integral approach for three-phase ternary systems with no success38 and, recently, it has been validated and generalized for multicomponent and multiphase systems.39 As a consequence of the

10.1021/ie1010383  2010 American Chemical Society Published on Web 09/10/2010

Ind. Eng. Chem. Res., Vol. 49, No. 20, 2010

analysis of the extension of the above area methods, a new procedure based on the tangent plane intersection concept (TPI) was proposed by Hodges et al.,38 and different solution techniques for LL, VL, and VLLE calculation of some twoand three-phase binary and ternary systems were given. However, results showed that the proposed techniques were sensitive to the position of the overall mixture and the size of the multiphase region for multiphase systems.40 The extensive work dedicated to the design of equilibrium calculation strategies proves the great difficulty of this task. In fact, many problems can arise when solving the equilibrium condition whatever the method used. Because of this, the improvement and development of equilibrium algorithms for data computation are still necessary. Convergence lack is a frequent problem due to the nonlinearity of the equations; trivial solutions are also common, and convergence to the global minimum is not always guaranteed. Commercial equilibrium data correlation tools, such as those included in the chemical process simulation packages, do not allow finding out about their internal procedure details. They only allow changing some operational parameters (iteration limits, deviation definition, equilibrium condition, etc.) and, in the case of convergence deficiency during data correlation, which too often occurs, these programs do not permit finding out if this is due to the initial set of model parameters proposed, a lack of capability of the model itself, an algorithm error, or a problem of the numerical method used. However, this information would be valuable to improve the regression process and to guarantee that the result obtained is the best possibility that the chosen model can provide for the system under study. In a previous paper,41 a discussion was presented regarding some illustrative problems that arise using commercial software packages for the correlation of type 1 LLE ternary systems and also inconsistencies that can be found in the literature. These facts reveal the high complexity of the subject and encourage more effort investment to improve the procedures for equilibrium data correlation. The limitation in the data correlation success (of any property or equilibrium data) must not be the calculation tool used but the model selected. Reliable algorithms must be used in the data correlation tools to obtain the corresponding parameter values that reproduce the experimental system behavior as long as the model used allows.1 In the present paper, a robust algorithm for the simultaneous correlation of condensed equilibrium data has been developed. This algorithm allows the correlation of different types of ternary systems (types 0, 1, 2, 3 and 4) to enable a possible extension of the existing correlation tools and also includes some strategies based on the topology analysis of the Gibbs energy of mixing function (gM) that overcome frequent inconsistencies found during equilibrium data correlation. 2. Algorithm Description The correlation algorithm has been devised for the correlation of condensed equilibrium data (LLE, LLLE, LLSE, etc.) for different types of ternary systems at constant T and P. The general aspects of this algorithm and the new calculation strategies that this includes will be described in the next sections. 2.1. General Aspects of the Correlation Algorithm. a. Simultaneous Correlation of All the Equilibrium Regions of the System. A flexible model should be capable of representing, with a unique set of parameters, all of the equilibrium data of a system, at least at constant temperature and pressure, regardless of the aggregation state and the number of phases present. Simultaneous correlation parameters are very

10101

advantageous for design calculations, because they can be used with certainty in the entire composition range. This practice avoids serious problems of inconsistencies that frequently arise with some parameter values that have been obtained using specific equilibrium data in limited composition regions, that is, when they are included in chemical process simulation software. However, this is not a common practice when multiple equilibrium regions are present. For example, most papers dealing with experimental data determination where solid and liquid phases coexist consider only LLE regions when trying to fit experimental data with a model.42-46 The algorithm developed has been successfully applied to the simultaneous correlation at constant T and P of type 3 ternary systems47 with LLLE (tie-triangle) regions and more complex type 4 systems in which LLE, LSE, and LLSE regions are simultaneously present.48 Moreover, the same algorithm has been used for more complex systems with hydrated salts with the simultaneous presence of LSE, LShE, LLShE, LShSE, and LLE regions.49 To make possible the simultaneous correlation of so many diverse equilibrium regions, the equilibrium over the whole composition domain is calculated and adequately compared with the experimental data for each of the model parameter sets tested during the optimization process, as explained below. b. Models for the Gibbs Energy of Mixing. The correlation algorithm proposed in the present paper allows the use of any model to formulate the Gibbs energy of mixing. Some illustrative examples found in the next sections correspond to the use of the NRTL model, but any classic or even user-defined equation can be otherwise used. Thus, it is possible to use the proposed calculation algorithm as an equilibrium data correlation utility to evaluate new more capable models for the Gibbs energy of mixing, such as the one presented in a previous paper.1 The Complex Step Derivative Approximation Technique50 has been used to precisely generalize numerical first and second derivatives in the numerical method used to solve the equilibrium condition equations. This is very useful when user-defined equations are proposed to formulate the Gibbs energy of mixing because it avoids the necessity of programming extensive analytical expressions for each model. c. Equilibrium Criterion. The equilibrium criterion used is the tangent plane condition that, as discussed in a previous work,51 describes the equilibrium very precisely compared to the isoactivity criterion. This condition is based on the one proposed by Iglesias et al.52 for LLE, and its formulation has been generalized for any kind of equilibrium between condensed phases (LSE, LLLE, LLSE, etc.). The application of this criterion reduces to the simultaneous solution of a system of (c - 1)np nonlinear equations,47,49 where c is the number of components and np is the number of phases. The Newton-Raphson method is the numerical method used to solve this system of equations. The chosen reference state for all of the components of the system in order to calculate the gM function is the pure liquid component at the same T and P of the system. When the system contains only liquid phases, the Gibbs energy for the liquid mixtures gM is formulated by the selected model. However, if solid phases are present, the Gibbs energy for the anhydrous solid (gS) is calculated as the chemical potential change from the pure liquid to the pure solid at the system T and P conditions, using a thermodynamic cycle.53 In the case of hydrated solids (Sh) being present, their Gibbs energy of formation with respect to the reference state (hypothetical subcooled molten salt and liquid water) is calculated using a resemblance of salt hydration

10102

Ind. Eng. Chem. Res., Vol. 49, No. 20, 2010

to gas adsorption on solid surfaces (BET multilayer gas adsorption model) as proposed by Zeng and Voigt.54 d. Objective Function and Optimization Technique. The phase equilibrium data correlation procedure requires the comparison between the experimental and calculated equilibrium data for each given set of parameters tested during the optimization process. Because the equilibrium among all liquid and solid phases is simultaneously considered for the complete system, the objective function OF to be minimized is expressed as the sum of two contributions: LLE region (OFLL), and the regions where a solid (anhydrous or hydrated) is present in the equilibrium (OFLS), when this applies:

in the liquid-liquid equilibrium calculations for the LL tie-lines (or LLS and LLSh tie triangles), a method to constrain the composition zones of search is used to get rid of most of the problems arising from multiple solution roots and the need for smart composition initial guesses. The procedure takes advantage of the topological information given by the second derivative of the Gibbs energy of mixing function (gM) about the surface curvature changes. The stability criterion (including metastability) for a homogeneous phase with c components is that the Hessian matrix of the gM function (σ) be positive definite, which corresponds to the Gibbs energy surface being convex:56

OF ) OFLL + OFLS ) nLL

3

np

∑ ∑ ∑ [(x

i,n,p)exptl

- (xi,n,p)calcd]2 +

n)1 i)1 p)1 nLS

σ)

(1)

3

∑ ∑ [(x

i,n,L)exptl

- (xi,n,L)calcd]2

n)1 i)1

In eq 1, xi,n,p is the molar fraction of component i in phase p for the LL tie-line n and xi,n,L is the molar fraction of component i in the liquid phase L of the LS tie-line n; nLL and nLS denote the number of tie-lines in the LL and LS equilibrium regions, respectively; np is the number of liquid phases (np ) 2); subscripts exptl and calcd denote, respectively, the experimental and calculated equilibrium data. For the correlation of tie-triangle data, in terms of the OF, these are considered to be divided into three pairs of tie-lines (each with np ) 2). The minimization of the objective function OF, eq 1, is performed by applying the Simplex optimization method.55 It must be remarked that the algorithm has been developed for ternary systems, but it could be extended to four or more components. 2.2. Novel Strategies Included in the Algorithm Proposed. As commented before, several approaches have been proposed for the computation of solutions to the equilibrium problem when it is formulated as the minimization of the Gibbs energy of mixing. However, due to the complexity of the issue and the nonlinear nature of the equations used to formulate the gM function, it is not possible to guarantee convergence to the solution. The result obtained is highly dependent on both the supplied initial points and the procedure used, even when the model parameters are known and the goal is to calculate or predict the number and composition of the equilibrium phases, as widely discussed by McDonald and Floudas.26 Obviously, the difficulty and the incidence of problems, that is, convergence deficiencies or inconsistent results, are higher when the aim of the calculation is to estimate the values for the model parameters by correlation of the experimental equilibrium data. The main problems found when using correlation tools included in commercial chemical process simulation software deal with two kinds of problems: (a) The program output shows a nonconvergence failure and the only chance for the user is repetition with another optimization algorithm or another initial parameter set. If no solution is found after some trials, the user gets to a “desperate point of uncertainty” about what to do next. (b) Convergent results that properly fit partial experimental data provide nonconsistent equilibrium behavior outside the composition region of data used through correlation. In this section, some strategies devised to minimize the above problems, which have been included in the algorithm presented in this paper, are explained. a. Limited Composition Range for the Liquid-Liquid Equilibrium Root Search. To find the coexisting compositions

| | ∂2gM ∂2x1

∂2gM ∂x1∂x2

∂2gM ∂x2∂x1

∂2gM ∂2x2

>0

(2)

When a fold appears in the gM surface, convexity is lost and the system splits into two phases. The solution for σ ) 0 gives the inflection points curve of the gM surface (spinodal curve), which separates the intrinsically unstable and stable equilibrium regions. The plait point condition for ternary binodal curves was established by Reid and Model57 as

σ)

| | | | ∂2gM ∂2x1

∂2gM ∂x1∂x2

∂2gM ∂x2∂x1

∂2gM ∂2x2

) 0 and δ )

∂σ ∂x1

∂σ ∂x2

2 M ∂2gM ∂ g ∂x2∂x1 ∂2x2

)0

(3)

In the algorithm proposed, the solutions for σ ) 0 and for δ ) 0 are estimated in the entire composition range, for each of the model parameter sets tested during the optimization process. This can be conveniently done by evaluating σ and δ changes of sign in a grid of compositions to avoid any convergence dependence on the numerical solution method during the correlation process. The inflection points curve, together with the plait point location, is used to confine both the composition area where overall composition mixtures for phase splitting calculation are located (σ < 0) and the regions where stable coexisting equilibrium phases are searched (σ > 0). As an example, Figure 1 shows the inflection points curve (σ ) 0) and also the curve where δ ) 0 for a type 1 (a), and an island type (b) systems obtained using the NRTL model. The intersection points of both curves represent the plait points of the system, which allows an optimal distribution of unstable mixture compositions. The tangent line to the inflection points curve at the plait point can be used to define the two composition coexisting areas, as those shaded in Figure 1, where each of the coexisting LLE compositions, for each unstable overall composition mixture, will be restricted during the equilibrium calculations. This performance excludes trivial solutions, improves the convergence, and is less time-consuming because it avoids searching for equilibrium splitting in the composition range with no curvature changes. Besides, it possesses an additional advantage that is better understood in the context of the next section, which is that it makes possible the identification of all of the LL equilibrium regions of the system, in the entire composition range, for each of the model parameter sets. b. New Criteria for Matching the Calculated Tie-Lines to the Experimental Data for the Objective Function Comparison. Correlation data is an optimization process that

Ind. Eng. Chem. Res., Vol. 49, No. 20, 2010

Figure 1. Restriction of the composition areas for LLE calculations. NRTL parameters: (a) A12 ) -50, A21 ) -300, A13 ) 2059.9, A31 ) 246.37, A23 ) -50, A32 ) -1000, R ) 0.2, 25 °C (type 1); (b) A12 ) 450, A21 ) -1500, A13 ) 523.05, A31 ) -230.72, A23 ) 1278, A32 ) -417.09, R ) 0.2, 20 °C (island type).

Figure 2. Intermediate step of the LLE correlation with NRTL for the ternary type 0 (island) system water (1) + dimethyl sulfoxide (2) + tetrahydrofuran (3) at 20 °C (ref 60, p 396) during the regression process, still far from an appropriate fitting. Parameters as in Figure 1b.

can reach many intermediate steps where parameter sets very far from the solution must be evaluated. Because LLE phase splitting occurs in a restricted composition region, unlike LVE calculations in which it occurs in the whole composition diagram for conditions below the critical, for many sets of parameters evaluated during the iteration process, the splitting or heterogeneous region can be very far from the experimental data, or it can even not exist. When this is taken into account, the comparison criteria used for the objective function between experimental and calculated data, eq 1, are very significant. Some of the existing literature methods generally accept the use of some experimental composition data (i.e., setting the middle point of the experimental tie-line or one concentration in one phase equal to its experimental value58,59), but this may lead to uncertainties and convergence problems through the comparison process, for instance, if the parameter set reproduces a homogeneous region in the experimental LLE zone.

10103

Figure 3. Intermediate step of the LLE correlation with NRTL for the ternary type 1 system methanol (1) + diphenylamine (2) + cyclohexane (3) at T ) 25 °C (ref 60, p 129) during the regression process, still far from an appropriate fitting, which generates an additional nondesirable island type LL region. A12(K) ) -356.64, A21 ) -1576.1, A13 ) 250.07, A31 ) 698.19, A23 ) -1600, A32 ) -4000, R ) 0.2.

In the example of Figure 2, an intermediate step during the NRTL correlation of the type 0 (island type) system water (1) + dimethyl sulfoxide (2) + tetrahydrofuran (3) at 20 °C (ref 60, p 396) is represented, which corresponds with the one shown in Figure 1b. As can be observed, the result is still far from an appropriate fitting. The LLE island region calculated for that intermediate set of parameter values is smaller than the experimental one. Therefore, if overall composition mixtures for calculating phase splitting were selected with the experimental LL binodal curve taken into account, or any calculated composition in a phase was set to be equal to the experimental one, many convergence errors would be obtained. Generally, to overcome these situations, those parameter sets are denied or a penalty term is added to the comparison objective function, but these alternatives introduce discontinuities and do not guide the optimization procedure through a realistic path. Figure 3 shows a more complicated but possible situation, which can occur during the optimization. This example shows an intermediate step during the LLE correlation for the type 1 ternary system methanol (1) + diphenylamine (2) + cyclohexane (3) at T ) 25 °C (ref 60, p 129), where the binary parameters have been restricted to correspond with the real behavior of the system (1-3 partially miscible; 2-3 and 1-2 totally miscible), as explained under section 2.2.c. Despite this restriction, an intermediate parameter set has been found that reproduces one additional LLE region (island type) that does not appear in the experimental phase behavior of the system. If phase equilibrium is not evaluated during the optimization in the entire composition range, parameters sets that give a very good approximation to the experimental LLE data, but that also predict this type of nondesirable region, could be erroneously validated. These parameter sets could lead to serious inconsistencies, that is, when they are used in process simulation tools.

10104

Ind. Eng. Chem. Res., Vol. 49, No. 20, 2010

is being evaluated. As a consequence, the optimization path works realistically and without discontinuities. As previously discussed, the number of experimental and calculated heterogeneous regions can be different, as in the example shown in Figure 3. Because for each iteration parameter set the second step requires the evaluation of the phase equilibrium behavior of the system in the entire composition range, the detection of “additional” heterogeneous regions is possible. The OF to be minimized is the sum of the contribution of all calculated equilibrium regions and, consequently, is representative of the accuracy of the current parameter set that is being evaluated. For the intermediate iteration step presented in Figure 3, eq 1 can be written as follows: OF ) OFtype 1 region + OFtype 0 region ) nLL

3

np

∑ ∑ ∑ [(x

i,n,p)exptl

- (xi,n,p)calcd, type 1]2 +

n)1 i)1 p)1

Figure 4. Intermediate step of the LLE correlation with NRTL for the type 1 ternary system water (1) + acetic acid (2) + trichloroethene (3) at T ) 25 °C (ref 60, p 152). Experimental and calculated LL tie-lines together with interpolated data are represented. Parameters as in Figure 1a.

The correlation algorithm presented in this paper tries to overcome these problems. We propose that all of the heterogeneous equilibrium regions, in the entire composition range, be obtained and adequately evaluated in the objective function (eq 1) for each parameter set tested during the correlation process. In this sense, a novel matching criterion is devised to always accomplish a piece of calculated equilibrium data to compare with the experimental one, despite being far from the experimental heterogeneous region. That is, we propose to carry out a calculated tie-line assignment for each of the experimental data, which must also be representative of the global equilibrium behavior of the system. Consequently, an OF value (eq 1) that represents the true quality of the evaluated parameter set is always obtained. The proposed strategy consists of the steps below and has been schematized in Figure 4 for a simple type 1 system: 1. An adequate number of tie-lines must be uniformly distributed in each of the heterogeneous regions of the system by experimental data interpolation. In the example shown in Figure 4, 10 tie-lines have been interpolated that uniformly describe the LL experimental region. 2. A comparable set of tie-lines must be calculated using the corresponding model parameter set that is evaluated in each iteration. The σ ) 0 solutions (spinodal curves) allow identification of the number and location of the heterogeneous regions. The plait points location (simultaneous σ ) 0 and δ ) 0, eq 3) are used to uniformly distribute overall composition mixture points for LLE calculations. For example, for a type 1 binodal curve an optimal distribution direction for overall composition mixture points is the one defined between the middle point of the partially miscible binary and the calculated plait point. For a type 0 system, mixtures should be distributed between both plait points; for a type 2 system the optimum direction should be defined, for instance, by the middle points of the two partially miscible binaries. In Figure 4, 10 tie-lines have been uniformly calculated through the LLE region, using an intermediate NRTL parameter set during the correlation process, to be compared with the corresponding experimental ones. 3. The calculated and experimental (interpolated) tie-lines are matched to obtain the OF value for the current iteration (eq 1). The objective function obtained with this procedure will always be representative of the quality of the model parameter set that

nLL

3

(4)

np

∑ ∑ ∑ [0 - (x

2 i,n,p)calcd, type 0]

n)1 i)1 p)1

Equation 4 takes into account not only the type 1 region, which is somewhat deviated from the experimental one, but also the island type undesired region, which has no comparable experimental data. At this point, it is important to remark again that this routine, included in the phase equilibrium data correlation tool, avoids inconsistencies and improves convergence because it guides the optimization path properly. c. Empirical Restrictions for the Binary Parameters of the Model. A successful equilibrium data correlation requires a coherent behavior among the calculated and experimental data in all equilibrium regions. Partial or total miscibility behavior of the binary subsystems and all of the equilibrium ternary regions must be faithfully reproduced by the model, at least qualitatively. However, parameters accurately describing certain equilibrium regions, but that are inconsistent with the behavior of the system in the entire composition range, can be found in the literature. In a previous work, we showed some examples of correlations where the model parameter sets predicted incoherent miscibility behaviors of some binary subsystems, although the experimental binodal curves were accurately defined.61 For example, some type 2 ternary systems are correlated with NRTL binary parameters that give three partially miscible (type 3) binary subsystems (ref 60, p 71) or type 1 ternary systems are correlated with NRTL binary parameters that give two partially miscible (type 2) binary subsystems (ref 60, p 436). These inconsistencies occur because no restriction is imposed to the binary parameter values during the data correlation and the final parameters obtained as solution are not validated in the entire composition range. In our opinion, partial or total miscibility should be imposed to the binary constituent subsystems during the correlation procedure. The equations corresponding with these restrictions may be obtained by using the information given by the topology of the gM surface. A systematic analysis of the curvature changes in the gM binary function can result in empirical conditions or parameter value ranges that determine phase splitting (LL) or miscibility (L) behavior. These restrictions can avoid parameters that, despite accurately describing certain equilibrium region, are inconsistent with the global behavior of the system in the whole composition diagram. As an example, the restrictions for the NRTL equation deduced for the binary interaction Aij parameter

Ind. Eng. Chem. Res., Vol. 49, No. 20, 2010

values (K) in a previous paper more accurate definition:

48

10105

are presented below, with a

homogeneous binary mixture (total miscibility): Aij < F (Aji) (5) heterogeneous binary mixture (partial miscibility): Aij > F (Aji) (6) with F(Aji) ) -4.46564 × 10-8A3ji + 2.95745 × 10-4A2ji 1.20662Aji + 766.908 ∀Aji ∈ [-1500, 2600] F(Aji) ) -0.547Aji + 237.992 ∀Aji > 2600 F(Aji) ) -1.829Aji + 435.760 ∀Aji < -1500 For the singular case of type 0 (island type) systems, in which three binary subsystems are totally miscible but a ternary LL region exists, we also obtained the restrictions for NRTL parameters to reproduce this ternary behavior.62 The gM curve for one of the three binary subsystems must be much larger (dissimilar binary gM curve) than the other two binary gM curves. The following restrictions allow this behavior, 1 and 2 being the dissimilar binary subsystems: A12 + A21 0

(8)

A23 + A32 > 0

(9) M

Analogous relationships could be deduced for other g models and used to correlate condensed phase equilibrium data to guarantee the consistent behavior of the binary subsystems. d. Vector Method. Coexisting liquid phases will become identical in the plait point. Therefore, the surface is very flat in this region, and as a consequence the numerical method to solve the equilibrium condition is prone to fail in this zone, which complicates equilibrium calculation. Some papers have been published to allow an adequate representation of the equilibrium in this region.63,64 In these papers, a semitheoretical correction is added to a conventional expression for the excess Gibbs energy (i.e., van Laar or NRTL). This correction is significant in the critical region but not elsewhere. In a previous paper,51 we approached this problem by focusing on the calculation procedure close to the plait point and, as a conclusion, we proposed an extension of the vector method,23 which led to a good representation of the binodal curve in this region. The advantage of this strategy is that it can be used with any model without the necessity of modifications. This method has been implemented in the proposed algorithm to be used when LL calculations are required close to the plait point. This alternative guarantees the convergence of the LL calculation and can be used to obtain a very good definition of the binodal curve in the area around the plait point. e. Geometrical Method To Estimate Initial Values for LLL Phase Compositions. For type 3 ternary systems, a direct geometrical method based on the tangent plane intersection condition has been included in the proposed algorithm to calculate the tie-triangle. This method, described in detail in a previous work,47 makes use of the concept that phases to be in equilibrium must be contained in a common plane, which has to be tangent to its gM surface. Three composition zones, each including one of the liquid phases corresponding to the LLLE tie-triangle, are separated by a sequential intersection of the

Figure 5. Global scheme of the proposed correlation algorithm including the novel strategies presented.

surface with successive planes. This procedure determines the number of existing phases and locates good initial values very close to equilibrium compositions that can be used as initial estimations for the tie-triangle iterative calculations. It provides the advantage that results do not depend on the overall composition mixture, and the trivial solution cannot be found because separated regions containing the searched phases are always determined. 3. Highlights of the Proposed Correlation Algorithm The philosophy of the proposed equilibrium correlation algorithm is schematized in Figure 5. The first step is the experimental data treatment, according to section 2.2.b. Characteristic data are chosen to identify each type of experimental system: the partially miscible binary subsystem and plait point for a type 1 system, the two partially miscible binaries for a type 2 system, and the partially miscible binaries and tie-triangles for type 3 and 4 systems. Then interpolated equilibrium data are obtained to uniformly describe all of the existing equilibrium regions, as required by the proposed comparison procedure used in the objective function (section 2.2.b). The solid chemical potentials are calculated if solid phases are present in the system (section 2.1.c). The model to be used for the equilibrium data correlation is then selected. According to the system nature, parameter restrictions are imposed to accomplish coherent miscibility behavior of binary subsystems (section 2.2.c). Also, and as a first approximation, the parameters for the partially miscible binary subsystems are calculated using exclusively binary data and kept constant to fasten the optimization procedure. When an approximate fitting is obtained, the optimization method varies the whole set of parameters to accurately fit all of the equilibrium data.

10106

Ind. Eng. Chem. Res., Vol. 49, No. 20, 2010

Table 1. Compilation of the Published Correlation Results Obtained Using the Proposed Algorithm improvement

correlation reference

LL LL LL LL LL LL LL LL LL LL LL LL

not correlated not correlated not correlated not correlated not correlated not correlated not correlated not correlated inconsistent correlation inconsistent correlation inconsistent correlation inconsistent correlation

61 61 61 61 61 61 61 61 61 61 61 this work

2

LL

inconsistent correlation

this work

0 0 3 3 3 4 4 4 4 4

LL LL LLL, LL LLL, LL LLL, LL LL, LLS, LS LL, LLS, LS LL, LLS, LS LS (LL, LLSh, LS, LSh, LSSh)

impossible impossible impossible impossible impossible impossible impossible impossible impossible impossible

62 62 47 47 47 48 48 48 48 49

system

components 1 + 2 + 3

type

1 2 3 4 5 6 7 8 9 10 11 12

furfural + formic acid + water, 35 °C water + acetic acid + dichloromethane, 19 °C toluene + acetaldehyde + water, 17 °C 1-butanol + methanol + water, 60 °C 4-methyl-2-pentanone + acetonitrile + water, 30 °C propenoic acid + propanoic acid + water, 70 °C hexane + 2,5-dimethyl-THF + 1,2-ethanodiol, 120 °C trichloroethene + furfural + water, 20 °C propyl acetate + formic acid + water, 35 °C acetonitrile + 1-propanol + hexane, 25 °C 1-hexanol+ nitromethane + water, 40 °C hexane + 9-octadecenoic acid (cis), methyl ester + acetic acid, nitrile, 25 °C propanoic acid, nitrile + heptane + octane, 1,8-OXY, perfluoro, 25 °C water + dimethyl sulfoxide + tetrahydrofuran, 20 °C trichloroacetic acid + antipyrine + water, 30 °C 1-nonanol + nitromethane + water, 23 °C 1-hexanol + nitromethane + water, 21 °C lauryl alcohol + nitromethane + glycol, 22 °C water + 1-butanol + NaCl, 25 °C water + 3-pentanol + NaCl, 25 °C water + acetone + NaCl, 25 °C water + ethanol + NaCl, 25 °C water + 1-butanol + LiCl, 25 °C

1 1 1 1 1 2 2 2 1 1 2 2

13 14 15 16 17 18 19 20 21 22 23

For each set of parameters proposed, the topology of the Gibbs energy of mixing function is evaluated in the entire composition space. The gM surface and its derivative functions are analyzed: plait points are determined and convex regions identified and separated for LL conjugated phase composition searching (section 2.2.a). Next, all of the equilibrium regions are calculated. For type 3 and 4 systems specific procedures have been developed to obtain LLLE (section 2.2.e), LSE, LShE, LLSE, LSSE, and LSShE data.47-49 For all types of systems, overall composition mixtures are uniformly distributed through the unstable regions for the two-phase equilibrium calculations. The vector method is used to precisely define the critical region of the binodal curves (section 2.2.d) when necessary. Because all of the existing calculated equilibrium regions are obtained for the entire composition domain, according to the proposed matching criterion of the tie-line assignment among experimental and calculated data, the OF (section 2.1.d) represents the true quality of the evaluated parameter set (section 2.2.b). This procedure avoids discontinuities and always guides the optimization procedure through a realistic path. 4. Test Examples In this section we present and discuss the examples that we have studied to analyze the behavior of the algorithm suggested. The results obtained in previous papers47-49,61,62 dealing with specific problems included in the algorithm, together with some new correlations, are summarized in Table 1. The existing equilibrium regions and the type of system based on the Treybal classification are shown for each case, together with the improvement obtained and the corresponding references. LLE data correlation for eight ternary systems including types 1 and 2 that have not been previously fitted by other authors (systems 1-8 classified in Table 1 as “not correlated”) have been tackled in another paper.61 A possible explanation of the problems that could arise during their correlations would be the difficulties in avoiding inconsistent solutions, such as the unrealistic Gibbs energy of mixing surfaces obtained with very negative values of NRTL binary parameters65 or the false

equilibrium regions

comparison comparison comparison comparison comparison comparison comparison comparison comparison comparison

prediction of LL splitting for some homogeneous binary subsystems. The use of restrictions for the total or partial miscibility for the binary subsystems (eqs 5 and 6) has allowed the regression with the NRTL model using the proposed algorithm. Some other systems found in that paper61 (9-11) had been previously correlated, but the results presented in the literature were inconsistent with the experimental behavior of the system (“inconsistent correlation” in Table 1). In the same way as before, the use of restrictions for the total or partial miscibility for the binary subsystems (eqs 5 and 6) has avoided inconsistent results such as those corresponding to system 9, which is type 1 in the Treybal classification. In this case, the parameter values published60 caused a nonexistent, small, two-liquid region spreading from the homogeneous 1-2 binary subsystem, producing a false type 2 system. The same type of problem, but with a greater inconsistency region, occurred for systems 10 and 11. NRTL parameters published60 for the 2-3 binary subsystem in system 10 and those corresponding to the 1-2 binary subsystem in system 11 predicted LLE splitting when, actually, these two binary pairs are totally miscible for all of the compositions at the T and P of the system. Results obtained for these systems with the proposed algorithm, compiled in ref 61, are totally coherent with the experimental behavior. The “inconsistent correlation” group of Table 1 has been completed with the NRTL correlation results for two additional type 2 ternary systems, for which inconsistent parameter sets have also been found in the literature:60 • System 12: hexane (1) + 9-octadecenoic acid (cis), methyl ester (2) + acetic acid, nitrile (3) at 25 °C (ref 60, page 190, and ref 66, page 110). • System 13: propanoic acid, nitrile (1) + heptane (2) + octane, 1, 8-OXY, perfluoro (3) at 25 °C (ref 60, p 460). For system 12, NRTL binary interaction parameters published60 corresponded to a type 1 ternary system, whereas the system is actually a type 2 (Figure 6a). For system 13, parameters produced three partially miscible binary subsystems giving two LL regions, as can be observed in Figure 7a, the system being classified as type 2. Figure 7b shows the detail of

Ind. Eng. Chem. Res., Vol. 49, No. 20, 2010

10107

Figure 6. Experimental and calculated (NRTL) LLE tie-line data for system 12 (hexane (1) + 9-octadecenoic acid (cis), methyl ester (2) + acetic acid, nitrile (3) at 25 °C): (a) using the binary parameters published in the literature (refs 60, p 190, and 66, p 110); (b) correlation results with the devised program.

the gM curve for the 1-2 binary subsystem that proves the undesirable LL splitting region produced by the published parameters.60 Experimental and calculated results using the algorithm proposed in the present paper have been graphically represented in Figures 6b and 7c. The correlation of the experimental equilibrium data using the restriction for the parameters included in the proposed algorithm, complete miscibility condition for the 1-2 binary subsystem (eq 5) and partial miscibility condition for the two binary subsystems 1-3 and 2-3 (eq 6), managed to obtain a set of consistent parameters. Table 2 shows the NRTL binary parameters (Aij) previously published60 and those obtained with the presented algorithm, together with the OF (eq 1). Standard deviations, calculated using the following equation, have also been included:

dev ) 100 ×



nLL

3

np

∑ ∑ ∑ [(x

i,n,p)exptl

- (xi,n,p)calcd]2

Figure 7. Experimental and calculated (NRTL) LLE tie-line data for system 13 (propanoic acid, nitrile (1) + heptane (2) + octane, 1,8-OXY, perfluoro (3) systems at 25 °C): (a) using the binary parameters published in the literature (ref 60, p 460); (b) detail of gM curve for 1-2 binary subsystem; (c) correlation results with the devised program.

n)1 i)1 p)1

6nLL

(10)

Finally, the correlations for type 0, type 3, and type 4 systems presented in the literature47-49,62 have been included in Table

10108

Ind. Eng. Chem. Res., Vol. 49, No. 20, 2010

Table 2. Correlation Results Using NRTL (r ) 0.2) for Systems 12 and 13a NRTL binary parameters Aij (K) system (Treybal type)

A12

A21

A13

A31

A23

A32

OF

dev (%) (eq 10)

-344.39 -773.86

529.14 482.29

590.41 557.14

-628.17 -708.85

1533.3 1845.2

3.13 × 10-1 2.95 × 10-3

7.22 0.70

319.85 337.54

809.92 1109.7

626.55 572.02

723.4 756.60

1.58 × 10-2 4.03 × 10-3

1.62 0.82

12 (2)

(a) (b)

-77.716 938.24

13 (2)

(a) (b)

368.82 344.12

192.26 175.04

a Binary parameters Aij (K), objective function (OF), and standard deviation (dev) are presented for the following cases: (a) inconsistent parameters in the literature;60 (b) consistent results using the present algorithm.

1 (systems 14-23). In these cases no contrast has been possible with the results of the existing commercial tools, because they do not allow the equilibrium data correlation for this more complex type of systems (“impossible comparison”). Requirements for the total miscibility for the binary subsystems (eqs 5 and 6) are essential for the correlation of island type systems (systems 14 and 15). No fitting was found in the literature for this type of system because commercial regression tools usually need the existence of at least one partially miscible binary subsystem. It was verified that if the miscibility conditions were not used, the final solutions of the optimization process were generally a type 1 or type 2 ternary system and, therefore, not consistent with the behavior of the system. The application of these miscibility restrictions, together with the relationships for island systems (eqs 7-9), has allowed the correlation of these singular systems with the proposed algorithm.62 The simultaneous correlations of systems with different equilibrium regions with a unique set of parameters have also been possible with the application of this algorithm. Selected experimental data for systems 16-18 (type 3) were simultaneously correlated, that is, those that define the shape of the phase diagram: the three binary subsystems and the LLL tietriangle.47 Results obtained showed that it was qualitatively possible to fit all LLE and LLLE regions with the NRTL equation using the same set of parameters. With regard to type 4 systems, due to the importance of the salting-out effect, some efforts were found in the literature trying to carry out the simultaneous correlation of equilibrium data over the whole composition range, but usually poor or inconsistent results were obtained, or even SLE and LLE regions were fitted separately.48 The present algorithm has permitted the simultaneous correlation of LLSE data for type 4 systems (19-22),48 even when hydrates are formed in the solid phase, giving a considerably complex equilibrium diagram (system 23).49 5. Conclusions A robust algorithm to correlate equilibrium data for condensed phases (LL, LLL, LS, LLS) has been presented. The algorithm can be used not only for type 1 and 2 ternary systems but also for the correlation of type 3, type island, and type 4 ternary systems, which are not usually included in the capabilities of many of the existing commercial tools for data reduction. Equilibrium data in all of the regions of the system are simultaneously correlated. Therefore, the binary parameters obtained for the model used are consistent with the phase behavior of the system in the entire composition range. The algorithm includes some novel strategies that have been intended to avoid the well-known convergence faults and also to guarantee the consistency of the parameters obtained. The routines developed have been demonstrated to be very powerful in this sense. Obviously, the programmed algorithm has great room for improvement because its

implementation has been carried out with the thinking of solving specific situations, and due to the complex nature of the phase equilibrium matter, new problems for particular situations can arise. The suggested ideas should be beneficial in helping experts improve the existing commercial tools that still fail in more situations than expected. Nomenclature Rij ) binary interaction parameters (K) for components i and j Aij ) nonrandomness NRTL factor gM ) Gibbs energy of mixing (dimensionless) R ) gas constant (J K-1 mol-1) np ) number of phases T ) temperature (K) P ) pressure (Pa) c ) number of components dev ) standard deviation nLL ) number of liquid-liquid tie-lines nLS ) number of liquid-solid tie-lines LLE ) liquid-liquid equilibrium LLLE ) liquid-liquid-liquid equilibrium LSE ) liquid-solid equilibrium LLSE ) liquid-liquid-solid equilibrium LLShE ) liquid-liquid-hydrated solid equilibrium LShE ) liquid-hydrated solid equilibrium VLE ) vapor-liquid equilibrium Superscripts exptl ) experimental value calcd ) calculated value Subscripts i, j, h ) components i, j, h p ) phase p n ) tie-line number L ) liquid phase S ) solid phase

Acknowledgment We gratefully acknowledge financial support from the VicePresidency of Research (University of Alicante) and Generalitat Valenciana (GV/2007/125). Literature Cited (1) Marcilla, A.; Olaya, M. M.; Serrano, M. D.; Reyes-Labarta, J. A. Methods for improving models for condensed phase equilibrium calculations. Fluid Phase Equilib. 2010, 296, 15. (2) Dechema Data Preparation Package (DPP) DECHEMA e.V. (3) ChemCAD 5.1.5. 2000. Chemstations Inc. (4) Treybal, R. E. Extraccio´n en fase lı´quida; UTEHA: México, 1968. (5) Becker, F.; Richter, P. Non-aqueous ternary mixtures with ‘island’ miscibility gaps. Fluid Phase Equilib. 1989, 49, 157. (6) Patterson, D. Polymer compatibility with and without a solvent. Polym. Eng. Sci. 1982, 22, 64.

Ind. Eng. Chem. Res., Vol. 49, No. 20, 2010 (7) Stateva, R.; Cholakov, P.; Georgi, S.; Galushko, A. A.; Wakeham, W. A. A powerful algorithm for liquid-liquid-liquid equilibria predictions and calculations. Chem. Eng. Sci. 2000, 55, 2121. (8) Kahlweit, M.; Strey, R.; Busse, G. Microemulsions: a qualitative thermodynamic approach. J. Phys. Chem. 1990, 94, 3881. (9) Andersen, J. G.; Koak, N.; de Loos, T. W. Influence of pressure on the LLLE in water + n-alkyl polyoxyethylene ether + n-alkane systems. Fluid Phase Equilib. 1999, 163, 259. (10) Huan, Z.; Van Bochove, G. H.; de Loos, T. W. Three-liquid phase equilibria in water + benzene + caprolactam + (NH4)2SO4 mixtures. Thermodynamics 2003, 49, 745. (11) Boddu, V.; Krishnaiah, A.; Viswanath, D. S. Liquid-liquid equilibria of the benzene + water + acetic acid ternary system and solubility of benzene in water: effect of calcium chloride. J. Chem. Eng. Data 2001, 46, 1172. (12) Veljkovic´, V. B.; Lazic´, M. L.; Stankovic´, M. Z. Repeated use of yeast biomass in ethanol fermentation of juniper berry sugars. World J. Microbiol. Biotechnol. 2006, 22, 519. (13) Roy, B. C.; Awual, M. R.; Goto, M. Effect of inorganic salts on ternary equilibrium data of propionic acid-water-solvents systems. J. Appl. Sci. 2007, 7, 1053. (14) Dionysiou, D.; Tsianou, M.; Botsaris, G. Extractive crystallization for the production of calcium acetate and magnesium acetate from carbonate sources. Ind. Eng. Chem. Res. 2000, 39, 4192. (15) Carreo´n-Caldero´n, G.; Fernando Barraga´n-Aroche, J.; Bazu´a-Rueda, E. Algoritmo para equilibrio multifa´sico, VLL y LLL, de sistemas muticomponentes con ecuaciones de estado bicu´bicas. Tecnol. Cienc. 2002, 17, 112. (16) Wakeham, W. A.; Stateva, R. P. Numerical solution of the isothermal multiphase flash problem. ReV. Chem. Eng. 2004, 20, 1. (17) Cairns, B. P.; Futzer, I. A. Multicomponent three-phase azeotropic distillation. 2. Phase-stability and phase-splitting algorithms. Ind. Eng. Chem. Res. 1990, 29, 1364. (18) Seider, W. D.; Widago, S. Multiphase equilibria of reactive systems. Fluid Phase Equilib. 1996, 123, 283. (19) Boston, J. F.; Brit, H. I. A radically different formulation and solution of the single-stage flash problem. Comput. Chem. Eng. 1978, 2, 109. (20) Michelsen, M. L. The isothermal flash problem. Part I. Stability. Fluid Phase Equilib. 1982, 9, 1. (21) Michelsen, M. L. The isothermal flash problem. Part II. Phasesplit calculation. Fluid Phase Equilib. 1982, 9, 21. (22) Eubank, P. T.; Hall, K. R. An equal area rule and algorithm for determining phase compositions. AIChE J. 1995, 41, 924. (23) Eubank, P. T.; Elhassan, A. E.; Barrufet, M. A.; Whiting, W. B. Area method for prediction of fluid-phase equilibria. Ind. Eng. Chem. Res. 1992, 31, 942. (24) Nichita, V.; Gomez, S. Efficient location of multiple global minima for the phase stability problem. Chem. Eng. J. 2009, 152, 251. (25) McDonald, C. M.; Floudas, C. A. Global optimization and analysis for the gibbs free energy function using the UNIFAC, Wilson, and ASOG equations. Ind. End. Chem. Res. 1995, 34, 1674. (26) McDonald, C. M.; Floudas, C. A. Global optimization for the phase stability problem. AIchE J. 1995, 41, 1798. (27) McDonald, C. M.; Floudas, C. A. Global optimization for the phase and chemical equilibrium problem: application to the NRTL equation. Comput. Chem. Eng. 1995, 19, 1111. (28) McDonald, C. M.; Floudas, C. A. GLOPEQ: a new computational tool for the phase and chemical equilibrium problem. Comput. Chem. Eng. 1997, 21, 1. (29) Sun, A. C.; Seider, W. D. Homotopy-continuaton method for stability analysis in the global minimization of the Gibbs free energy. Fluid Phase Equilib. 1995, 103, 213. (30) Jalali, F.; Seader, J. D. Homotopy continuation method in multiphase multi-reaction equilibrium systems. Comput. Chem. Eng. 1999, 23, 1319. (31) Hua, J. Z.; Brenneke, J. F.; Stadherr, M. A. Reliable prediction of phase stability using an interval Newton method. Fluid Phase Equilib. 1996, 116, 52. (32) Saber, N.; Shaw, J. M. Rapid and robust phase behaviour stability analysis using global optimization. Fluid Phase Equilib. 2008, 264, 137. (33) Shyu, G. S.; Hanif, N. S. M.; Alvarado, J. F. J.; Hall, K. R.; Eubank, P. T. Equal area rule methods for ternary systems. Ind. Eng. Chem. Res. 1995, 34, 4562. (34) Hanif, N. S. M.; Shyu, G. S.; Hall, K. R.; Eubank, P. T. Calculation of multi-phase equilibria using the equal area rule with application to hydrocarbon/water mixtures. Fluid Phase Equilib. 1996, 126, 53.

10109

(35) Shyu, G. S.; Hanif, N. S. M.; Hall, K. R.; Eubank, P. T. Maximum partial area rule for phase equilibrium calculations. Ind. Eng. Chem. Res. 1996, 35, 4348. (36) Elhassan, A. E.; Tsvetkov, S. G.; Craven, R. J. B.; Stateva, R. P.; Wakeham, W. A. A rigorous mathematical proof of the area method for phase stability. Ind. Eng. Chem, Res 1998, 37, 1483. (37) Elhassan, A. E.; Lopez, A. A.; Craven, R. J. B. Solution of the multiphase problem for pure component, binary and ternary systems using the Area Method. J. Chem Soc., Faraday Trans. 1996, 92, 4419. (38) Hodges, D.; Pritchard, D. W.; Anwar, M. M.; Gregory, P. Calculating binary and ternary multiphase equilibria: extensions of the integral area method. Fluid Phase Equilib. 1997, 130, 101. (39) Balogh, J.; Craven, R. J. B.; Stateva, R. P. The area method for phase stability analysis revisited: further developments. Formulation in terms of the convex envelope of thermodynamic surfaces. Ind. Eng. Chem. Res. 2007, 46, 1611. (40) Hodges, D.; Pritchard, D. W.; Anwar, M. M. Calculating binary and ternary multiphase equilibria: the tangent plane intersection method. Fluid Phase Equilib. 1998, 152, 187. (41) Marcilla, A.; Olaya, M. M.; Serrano, M. D. Comments on liquid-liquid equilibrium data regression. J. Chem. Eng. Data 2007, 52, 2538. (42) Macedo, E. A.; Skovborg, P.; Rasmussen, P. Calculation of phase Equilibria for solutions of strong electrolytes in solvent-water mixtures. Chem. Eng. Sci. 1990, 45, 875. (43) Li, Z.; Tang, Y.; Li, Y. Salting effect in partially miscible systems of n-butanol-water and butanone-water. 1. Determination and correlation of liquid-liquid equilibrium data. Fluid Phase Equilib. 1995, 103, 143. (44) Tang, Y.; Li, Z.; Li, Y. Salting effect in partially miscible systems of n-butanol-water and butanone-water. 2. An extended Setschenow equation and its application. Fluid Phase Equilib. 1995, 105, 241. (45) Govindarajan, M.; Sabarathinam, P. L. Salt effect on liquid-liquid equilibrium of the methyl isobutyl ketone-acetic acid-water system at 35 °C. Fluid Phase Equilib. 1995, 108, 269. (46) Escudero, I.; Cabezas, J. L.; Coca, J. J. Liquid-liquid equilibria for 2,3-butanediol + water + 4-(1)-methylpropyl)phenol + toluene at 25 °C. J. Chem. Eng. Data 1996, 41, 2. (47) Marcilla, A.; Olaya, M. M.; Serrano, M. D.; Velasco, R.; ReyesLabarta, J. A. Gibbs energy based procedure for the correlation of type 3 ternary systems including a three-liquid phases region. Fluid Phase Equilib. 2009, 281, 87. (48) Olaya, M. M.; Marcilla, A.; Serrano, M. D.; Botella, A.; ReyesLabarta, J. A. Simultaneous correlation of liquid-liquid, liquid-solid, and liquid-liquid-solid equilibrium data for water + organic solvent + salt ternary systems. Anhydrous solid phase. Ind. Eng. Chem. Res. 2007, 46, 7030. (49) Marcilla, A.; Reyes-Labarta, J. A.; Olaya, M. M.; Serrano, M. D. Simultaneous correlation of liquid-liquid, liquid-solid, and liquid-liquidsolid equilibrium data for water + organic solvent + salt ternary systems: hydrated solid phase formation. Ind. Eng. Chem. Res. 2008, 47, 2100. (50) Lai, K. L.; Crassidis, J. L.; Cheng, Y.; Kim, J. AIAA-2005-5944, AIAA/GN&C Conference and Exhibit; American Institute of Aeronautics and Astronautics (AIAA): San Francisco, CA, 2005. (51) Olaya, M. M.; Ibarra, I.; Reyes-Labarta, J. A.; Serrano, M. D.; Marcilla, A. Computing liquid-liquid phase equilibria: an exercise for understanding the nature of false solutions and how to avoid them. Chem. Eng. Educ. 2007, 41, 218. (52) 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, 229. (53) Prausnitz, J. M.; Lichtenthaler, R. N.; Gomes de Acevedo, E. Molecular Thermodynamics of Fluid-Phase Equilibria, 3rd ed.; Prentice Hall PTR: Upper Saddle River, NJ, 1999. (54) Zeng, D.; Voigt, W. Phase diagram calculation of molten salt hydrates using the modified BET equation. CALPHAD 2003, 27, 243. (55) Himmelblau, D. M. Process Analysis Statistical Methods; Wiley: New York, 1968. (56) 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, 47. (57) Reid, R.; Modell, M. Thermodynamics and Its Application; Prentice Hall: Englewood Cliffs, NJ, 1983. (58) Renon, H.; Asselineau, L.; Cohen, G.; Raimbault, C. Calcul sur ordinateor des Equilibres Liquide-Vapeur et Liquide-Liquide; Technip: Paris, France, 1971.

10110

Ind. Eng. Chem. Res., Vol. 49, No. 20, 2010

(59) Sørensen, J. M. Correlation of Liquid-Liquid Equilibrium Data. Thesis, Instituttet for Kemiteknik Danmarks Tekniske Højskole, Denmark, 1980. (60) Sorensen, J. M.; Arlt, W. Liquid-Liquid Equilibrium Data Collection; Chemistry Data Series, Vol. V, Part 2; DECHEMA: Frankfurt, Germany, 1980. (61) Reyes-Labarta, J. A.; Olaya, M. M.; Velasco, R.; Serrano, M. D.; Marcilla, A. Correlation of the liquid-liquid equilibrium data for specific ternary systems with one or two partially miscible binary subsystems. Fluid Phase Equilib. 2009, 278, 9. (62) Olaya, M. M.; Reyes-Labarta, J. A.; Velasco, R.; Ibarra, I.; Marcilla, A. Modelling liquid-liquid equilibria for island type ternary systems. Fluid Phase Equilib. 2008, 265, 184. (63) de Pablo, J. J.; Prausnitz, J. M. Thermodynamics of liquidliquid equilibria including the critical region. AIChE J. 1988, 34, 1595.

(64) Lee, B. H.; Qin, Y.; Prausnitz, J. M. Thermodynamic representation of ternary liquid-liquid equilibria near-to and far-from the plait point. Fluid Phase Equilib. 2006, 240, 67. (65) Simoni, L. D.; Lin, Y.; Brennecke, J. F.; Stadtherr, M. A. Reliable computation of binary parameters in activity coeficient models for liquidliquid equilibrium. Fluid Phase Equilib. 2007, 255, 138. (66) Sørensen, J. M.; Arlt, W. Liquid-Liquid Equilibrium Data Collection; Chemistry Data Series, Vol. V, Part 1; DECHEMA: Frankfurt, Germany, 1979.

ReceiVed for reView May 6, 2010 ReVised manuscript receiVed August 15, 2010 Accepted August 19, 2010 IE1010383