A Model-Free Approach Data Treatment of Vapor ... - ACS Publications

2134

Ind. Eng. Chem. Res. 2001, 40, 2134-2148

A Model-Free Approach Data Treatment of Vapor-Liquid Equilibrium Data in Ternary Systems. 1. Theory and Numerical Procedures Elizabeth Lam,†,‡ Andre´ s Mejıá,† Hugo Segura,*,† Jaime Wisniak,§ and Sonia Loras| Department of Chemical Engineering, Universidad de Concepcio´ n, P.O. Box 160-C, Correo 3, Concepcio´ n, Chile, Department of Chemical Engineering, Ben-Gurion University of the Negev, Beer-Sheva, Israel, and Departamento de Ingenierıá Quı´mica, Facultad de Quı´mica, Universitat de Vale` ncia, 46100 Burjassot, Valencia, Spain

A new data treatment approach for ternary systems is presented. The method is a model-free technique based on Barker’s equation and integration of the Gibbs-Duhem relation. It yields numerical information for the excess Gibbs energy and activity coefficients of a ternary system and its constituent binaries. Depending on the availability of data, the method is applicable to completely miscible isothermal or isobaric systems in vapor-liquid equilibrium, at low pressures, and is useful for selecting a GE model able to correlate accurately the data. Physical constraints for the curvature of the bubble surface in the vicinity of multicomponent azeotropic concentrations are deduced from the displacement theory. The relation of these constraints with the modelfree approach is discussed. Introduction Experimental vapor-liquid equilibrium (VLE) data are needed for both theoretical and practical purposes and are particularly useful when correlated by means of excess Gibbs energy (GE) models and/or equations of state (EOS). For practical applications, GE models constitute accurate tools for correlating experimental data. In fact, in engineering practice GE models are effective relations for designing separation systems. Properly fitted GE models are useful also for predicting physical properties of mixtures, such as self-diffusion coefficients and interfacial tension. From a theoretical viewpoint, the GE function brings light on the molecular interactions that are present in the liquid phase. Nowadays, statistical mechanics and molecular simulation are able to generate semitheoretically based GE models that, in turn, are advantageous for providing physical explanations about phase behavior. However, independently of the theoretical background of GE models (generally constrained to very idealized molecular considerations), practice has shown that quantitative predictions, or the physical interpretation of phase behavior, are possible when correlation models are adequate and accurate VLE data are available. The procedure for assessing and correlating VLE data (data treatment) includes consistency analysis and data reduction of the same. Both subjects are standard sections in every publication reporting experimental work on phase equilibrium. Data treatment may be performed by regressing VLE data for which the validity of a particular GE correlation is assumed. The GE model parameters are then calculated by minimizing the * Corresponding author. E-mail: [email protected]. † Universidad de Concepcio ´ n. ‡ Present address: Department of Chemical Engineering, Universidad Cato´lica del Norte, Avenida Angamos 0610, Antofagasta, Chile. § Ben-Gurion University of the Negev. | Universitat de Vale ` ncia.

deviations between predicted and experimental values. A common regression technique is based on Barker’s equations1 C

P)

∑ i)1

xiγiPi0

(1a)

Φi

xiγiPi0 yi ) Φ iP

(1b)

In eq 1a,b, C is the number of components of the mixture, xi and yi are mole fractions of component i in the liquid and vapor phases, respectively, γi is the activity coefficient, P0i is the vapor pressure of pure species, and Φi is a fugacity correction of the phases defined as2

Φi )

[

φˆ i V ˜ Li exp (P - P0i ) φi RT

]

(2)

where φi is the fugacity coefficient and V ˜ Li is the liquidphase volume for component i. According to eq 1a,b, the vapor pressure P of a mixture can be calculated if the temperature T, the liquid-phase mole fractions, and the GE function are known. In particular, the GE function is used to obtain analytical relations for activity coefficients. To predict P and yi, it is necessary to solve eq 1a,b. The parameters of the GE function are then calculated by minimizing the deviations (P, yi) from experimental data. Another technique for calculating the parameters of the GE model consists of a direct fit of activity coefficients, as calculated from experimental VLE data. Regression of experimental activity coefficients requires the assumption of a valid GE model. The parameters of the GE model are calculated by minimizing the deviations between calculated and experimental values,

10.1021/ie0007522 CCC: $20.00 © 2001 American Chemical Society Published on Web 04/07/2001

Ind. Eng. Chem. Res., Vol. 40, No. 9, 2001 2135

yielding a procedure that avoids phase equilibrium calculations. In contrast to previous approaches, model-free techniques are methods that do not require the assumption of a particular GE model. In the model-free approach, the GE function may be obtained numerically from experimental data by integrating the Gibbs-Duhem relation. The main feature of model-free techniques is their ability to overcome modeling pitfalls because, when performing data treatment, the right selection of an adequate GE model is also a part of the problem. The application of this approach for isothermal VLE data of binary systems is well described by Mixon et al.,3 Van Ness,4 Sayegh and Vera,5 and Wisniak et al.6 Nowadays, many VLE determinations are carried out under isobaric conditions. A relevant reason for this practice is that isobaric data meet direct application in the synthesis and design of separation technologies, particularly distillation. In addition, pressure is easier to control in dynamic equilibrium stills. However, from a theoretical viewpoint, the treatment of isobaric VLE data is limited to the availability of excess enthalpy data for the liquid phase. This fact is a serious limitation of every reduction method because excess enthalpy data are scarce, generally limited to narrow temperature ranges and, in addition, difficult to predict. Model-free approaches do not constitute an exception to this limitation, and to the best of our knowledge, they have not been applied systematically to the treatment of isobaric data. Ternary systems play an important role for analyzing the capability of prediction of multicomponent systems from their pertinent binaries. Usually, such a capability is tested using multicomponent GE models with parameters previously fitted to binary data. Little attention has been given to systems that are difficult to predict. Considering that inaccurate predictions may be originated in model weaknesses or inaccurate binary data, model-free techniques are potentially useful for treating complex ternary systems. Mixon et al.3 outlined the application of model-free techniques to isothermal ternary systems, although using a weak numerical approach (point relaxation), from which the treatment of isobaric systems is not clearly deduced. In this work, we present a rigorous extension of the model-free approach to isothermal and isobaric VLE data of ternary systems. In addition, we develop the numerical strategy required for solving efficiently the problem.

3

dgE )

3

ln γ3) dx1 + (ln γ2 - ln γ3) dx2 +

The Gibbs-Duhem equation for the excess Gibbs energy of a ternary system G ˜ E is2

-

H ˜E RT2

dT +

V ˜E RT

3

dP -

xi d ln γi ) 0 ∑ i)1

(3)

In eq 3, T is the temperature, R is the gas constant, H ˜E and V ˜ E are the excess enthalpy and volume, respectively, for the liquid phase, and γ is the activity coefficient. In addition, the total differential of the excess energy can be written as

xi d ln γi ∑ i)1

(4)

˜ E/RT. Replacement of eq 3 in eq 4 yields where gE ) G

dgE ) -

H ˜E V ˜E dT + dP + (ln γ1 - ln γ3) dx1 + RT RT2 (ln γ2 - ln γ3) dx2 (5)

From eq 5 we get now the following derivatives:

( )

)-

( )

)-

∂gE ∂x1

x2

∂gE ∂x2

x1

( )

+

( )

+

H ˜ E ∂T RT2 ∂x1

x2

H ˜ E ∂T RT2 ∂x2

x1

( )

+ ln γ1 - ln γ3

( )

+ ln γ2 - ln γ3

V ˜ E ∂P RT ∂x1 V ˜ E ∂P RT ∂x2

x2

(6a)

x1

(6b)

Parts a and b eq 6, together with the relation gE ) ln γi, constitute a linear system which can be solved for ln γi. Thus, the activity coefficients for each component of a ternary system are given by 3 ∑i)1 xi

[( ) ∂gE ∂x1

ln γ1 ) gE - (x1 - 1)

[( ) E

∂g ∂x1

ln γ2 ) gE - x1

x2

[( ) ∂gE ∂x1

ln γ3 ) gE - x1

] [( )

x2

]

+ ψx1 + (1 - x2)

x2

∂gE ∂x2

+ ψx1 - x2

[( ) ∂g ∂x2

∂gE ∂x2

x1

]

ψ x2

E

] [( )

+ ψx1 - x2

+

x1

+

x1

]

]

ψ x2

+ ψx2

(7)

where ψxi is defined as

ψxi )

( )

H ˜ E ∂T RT2 ∂xi

xj*i

( )

V ˜ E ∂P RT ∂xi

(8)

xj*i

For a ternary system eqs 1 and 7 yield the following nonlinear partial differential equation on gE:

P)

Theory

3

ln γi dxi + ∑xi d ln γi ) (ln γ1 ∑ i)1 i)1

{

[( ) { [( ) ]} ] [( ) ]} [( ) ] [(

x1P10 ∂gE exp gE + (1 - x1) Φ1 ∂x1 E

∂g ∂x2

+

ψx1 + (1 - x2)

∂gE ∂x2

x1

E

x1

∂g ∂x1

x2

] [( ) {

+ ψx1 -

0

+ ψx2

x2

x2

x2P2 ∂gE exp gE - x1 Φ2 ∂x1

x1

+ ψx2

+

+ ψx1 - x2

+

x2

x3P30 exp gE Φ3

)

∂gE ∂x2

x1

+ ψx2

]}

(9)

Equation 9 may be solved numerically in gE (as ˜ E are known continuexplained below) if P, T, V ˜ E, and H ous functions of x1 and x2. A potential feature of this

2136

Ind. Eng. Chem. Res., Vol. 40, No. 9, 2001

approach is that no excess Gibbs energy model is required for performing data treatment, thus avoiding possible modeling pitfalls in the analysis of ternary systems. In addition, the required input variables, P, ˜ E, and x, can be measured directly, with known T, V ˜ E, H uncertainties. As usual in experimental measurements, the mentioned input data are discrete, and it requires smoothing by means of interpolating functions in order to apply model-free techniques.3,5,6

Pj,i )

[

[

1. A Simplified Mixed Approach. In this approach the x1 derivatives are approximated by a simple finite difference

( ) ∂gE ∂x1

gE(x1) - gE(x1 - R) ≈ R x2

(10)

]}

{

]

]} [

(i - 1)

E gj,i

E gj,i-1

R

]

( )

(11)

where R is a sufficiently small number. The reason for combining different numerical definitions of the derivatives will become clear later, once the description of this approach has been completed. Let us now divide x1 and x2 in a grid of N equally sized intervals given by

R)

1 N-1

(12)

The concentration range for mole fractions can be parametrized by defining

[

x2 ) R(j - 1)

1eieN 1ejeN+1-i

(13)

The way in which the range for j has been defined in eq 13 ensures that the physical condition x1 + x2 e 1 be fulfilled. Equations 10, 11, and 13 yield the following discrete representation for eq 9:

E - gj-1,i + 2R

]}

(14)

In addition, the mole fractions for the vapor phase (eq 1b) are given by

y1,j,i )

[

{

0 R(i - 1)P1,j,i E exp gj,i + (1 - R(i - 1)) Φ1,j,iPj,i

]

[

E E E E gj,i - gj,i-1 gj+1,i - gj-1,i + ψx1,j,i - R(j - 1) + ψx2,j,i R 2R

y2,j,i )

{

[

0 E E R(j - 1)P2,j,i gj,i - gj,i-1 E + exp gj,i - R(i - 1) Φ2,j,iPj,i R

]

[

E E gj+1,i - gj-1,i + ψx2,j,i 2R

y3,j,i ) 1 - y1,j,i - y2,j,i

]} ]}

(15)

where the subindices j and i indicate that the property is evaluated at the mole fractions given by eq 13. Equation 14 may now be written as the objective function

ϑj,i )

[

[

{ ]

0 R(i - 1)P1,j,i E exp gj,i + (1 - R(i - 1)) Φ1,j,i

[

E E E E gj,i - gj,i-1 gj+1,i - gj-1,i + ψx1,j,i - R(j - 1) + R 2R 0 R(j - 1)P2,j,i E ψx2,j,i + exp gj,i - R(i - 1) Φ2,j,i

]}

{

]

[

E E E E gj,i - gj,i-1 gj+1,i - gj-1,i + ψx1,j,i + (1 - R(j - 1)) + R 2R 0 (1 - R(i - 1) - R(j - 1))P3,j,i E ψx2,j,i + exp gj,i -R Φ3,j,i

]} [

(i - 1) x1 ) R(i - 1)

{

E gj+1,i

ψx2,j,i

while the x2 derivatives are approximated by the central difference

gE(x2 + R) - gE(x2 - R) ∂g ≈ ∂x2 x1 2R

[

+ ψx1,j,i - R(j - 1)

ψx1,j,i + (1 - R(j - 1))

E

[

E E E E gj,i - gj,i-1 gj+1,i - gj-1,i + ψx1,j,i - R(j - 1) + R 2R 0 R(j - 1)P2,j,i E exp gj,i - R(i - 1) ψx2,j,i + Φ2,j,i

E E E E gj,i - gj,i-1 gj+1,i - gj-1,i + ψx1,j,i + (1 - R(j - 1)) + R 2R 0 (1 - R(i - 1) - R(j - 1))P3,j,i E exp gj,i -R ψx2,j,i + Φ3,j,i

Numerical Procedure An obvious numerical solution for eq 9 is the finite difference approach, where the partial derivatives are replaced by numerical approximations based on difference equations. Many formulas for numerical derivatives are available;7 they differ fundamentally on their accuracy for approximating the derivative. In the sections that follow, two alternatives for integrating eq 9 are proposed. Each can be applied depending on the complexity of the ternary system under analysis.

{ ]


]

[

{

E E E E gj,i - gj,i-1 gj+1,i - gj-1,i + ψx1,j,i - R(j - 1) + R 2R

ψx2,j,i

]}

- Pj,i (16)

that must be zero to satisfy the equilibrium condition imposed by Barker’s relation. Inspection of eq 16 indicates that if P, T, ψx1, and ψx2 are known functions of the concentration (parametrized by specific j and i E E , gj-1,i , values, as in eq 13), then its unknowns are gj,i-1 E E gj,i, and gj+1,i.


purposes. In addition, each side of the triangle (A-B, B-C, and C-A) represents a binary system whose gE function, in turn, may be solved using a specific modelfree technique for binary systems. In the discretization procedure, given by eqs 10 and 11, the partial derivatives in x1 are treated by means of a simple finite difference. Such an approach affects every j file in Figure 2. Therefore, to keep the numerical consistency of the integration method, it is necessary that the binary system that appears in the side C-A (binary 1 + 3), be integrated using the finite difference scheme given by eq 10. Because j ) 1 for the binary in the C-A side, eq 16 yields

ϑ1,i )

{

0 R(i - 1)P11,i E exp g1,i + (1 - R(i - 1)) Φ11,i

[ {

E E g1,i - g1,i-1 + ψx1,1,i R

Figure 1. Graphical interpretation of eq 16 about the nodes j and i.

Figure 2. Complete set of nodal points for the application of eq 16 to a ternary system: (O) ternary nodes; (0) binary nodes; (9) pure-component nodes.

Figure 1 shows the unknowns of eq 16 for fixed values of j and i. It should be noted that fixed values of j and i constitute a discretization node for which the Barker’s relation is applied. Similarly, Figure 2 shows the discretization grid of the whole range of mole fractions of a ternary system, showing also the nodes j and i illustrated in Figure 1. Every vertex of the triangle in Figure 2 represents a pure component, for which the excess Gibbs energy is zero. Consequently, the conditions

vertex A: x1 ) 1 w i ) N and j ) 1 w

E gN,1

)0

E )0 vertex B: x2 ) 1 w i ) 1 and j ) N w g1,N E )0 vertex C: x3 ) 1 w i ) 1 and j ) 1 w g1,1

(17)

can be used as boundary conditions for integration

exp

E g1,i

[

]}

+

0 (1 - R(i - 1))P31,i Φ31,i

E E g1,i - g1,i-1 - R(i - 1) + ψx11,i R

]}

- P1,i (18)

Starting with the boundary condition of vertex C in E ) 0), we see that for i ) 2 eq 18 depends only eq 17 (g1,1 E on the single variable g1,2 . If this scheme is repeated for increasing values of i, then it will always be possible to solve ϑ1,i in the single variable g1,iE, provided that the previous node ϑ1,i-1 has been solved for gEi-1. The numerical approach for solving eq 18 is discussed in detail in Appendix A. The remaining binaries of component 2 (1 + 2 and 2 + 3, on the sides A-B and B-C of the triangle) must be solved using the discretization given by eq 11. Table 1 shows the pertinent discretized relations for these binaries. The solution strategy may be based on the approach of Mixon et al.3 or that of Wisniak et al.6 It is seen then that the vertexes and the sides of triangle A-B-C constitute nodes where the values of the excess Gibbs energy are known before solving the properties of a ternary system. Consequently, the ternary problem is reduced to the nodes inside the subtriangle D-E-F. Let us now establish a solution strategy for this subtriangle. As shown in Figure 2, the solution of nodes j and i for a column with a fixed value of i (a condition in which x1 is constant), requires knowledge of the gE values of the previous column (i - 1). In fact, when the gE values are known along column i - 1, then eq 16 applied to nodes j and i has less unknowns. E E E ,gj,i ,gj+1,i} ≡ 0 ϑj,i ) ϑj,i{gj-1,i

(19)

Equation 19 may be applied then to the whole column i by varying j (which is equivalent to varying x2) in the range 2 e j e N - i. For this case, eq 19 gives way to the following nonlinear system: E E E ,g2,i ,g3,i }≡0 ϑ2,i ) ϑ2,i{g1,i E E E ,g3,i ,g4,i }≡0 ϑ3,i ) ϑ3,i{g2,i · · · E E E ϑN-i,i ) ϑN-i,i{gN-i-1,i,gN-i,i ,gN-i+1,i }≡0

(20)

E In eq 20 g1,i corresponds to a gE value for which j ) 1 E is a previously and, consequently, x2 ) 0. Thus, g1,i

2138


Table 1. Discretized Relations for Binary Systems Binary 1 + 2

discretization

[

{

0 E E (1 - R(N - i))P1,N+1-i,i gN+2-i,i - gN-i,i E + ψx2N+1-i,i exp gN+1-i,i - R(N - i) ϑN+1-i,i ) Φ1,N+1-i,i 2R

{

[

]}

0 E E R(N - i)P2,N+1-i,i gN+2-i,i - gN-i,i E + exp gN+1-i,i + (1 - R(N - i)) + ψx2,N+1-i,i Φ2,N+1-i,i 2R

]}

- PN+1-i,i

range 2eieN-1

boundary conditions E )0 gN,1 E g1,N )0

Binary 2 + 3

discretization ϑj,1 )

{

[

]}

0 E E R(j - 1)P2,j,1 gj+1,1 - gj-1,1 E + ψx2,j,1 exp gj,1 + (1 - R(j - 1)) Φ2,j,1 2R

{

[

E exp gj,1 - R(j - 1)

+

0 (1 - R(j - 1))P3,j,1 Φ3,j,1

]}

E E gj+1,1 - gj-1,1 + ψx2,j,1 2R

- Pj,1

range 2ejeN-1 boundary conditions E g1,1 )0 E gN,1 )0

calculated excess Gibbs energy value of the binary E is a system 1 + 3 in x1 ) R(i - 1). Similarly, gN-i+1,i previously calculated Gibbs energy value of the binary system 1 + 2 in x1 ) R(i - 1). Consequently, eq 20 represents a nonlinear system of N - i - 1 equations having the same number of unknowns. Numerical details for the solution of eq 20 are discussed in Appendix B. Briefly, eq 20 can be solved efficiently for every fixed i value using a Newton technique and, provided that the gE values for the nodes of column i 1 are known, using a tridiagonal matrix algorithm. The general procedure for a ternary system is as follows: 1. Select N and determine the size of the discretization grid for the ternary system. Calculate R ) (N - 1)-1. Recommended values of R for ternary systems are in the range 1 × 10-2 < R < 1 × 10-4. 2. Solve each binary system using the discretization given in step 1 and a model-free technique applicable to binary systems. Once step 2 has been completed, we have the solution for column i ) 1 (see Figure 2), representing the binary 2 + 3. Column 1 is required for solving column 2. Now, let i ) 2. 3. Solve the nonlinear system in eq 25 for column i using the procedure given in Appendix B. It is recomE mended to initialize the gE values of column i using gj,i E ) gj,i-1 (2 < j < N - i).

4. Let i ) i + 1 and repeat from step 3 until i ) N 2. 5. Convergent activity coefficients can now be calculated in every discretization node using

{

[

{

]}

E E gj,i - gj,i-1 + ψx1,j,i + R E E gj+1,i - gj-1,i (1 - R(j - 1)) + ψx2,j,i 2R

]}

[

E γ2 ) exp gj,i - R(i - 1)

{

γ3 ) exp

E gj,i

]

E E gj,i - gj,i-1 + ψx1,j,i R E E gj+1,i - gj-1,i R(j - 1) + ψx2,j,i 2R

E + (1 - R(i - 1)) γ1 ) exp gj,i

[

[

[

] ]

E E gj,i - gj,i-1 - R(i - 1) + ψx1,j,i R E E gj+1,i - gj-1,i R(j - 1) + ψx2,j,i 2R

[

]}

(21)

The procedure presented here, as well as the possibility of using a Newton method based on a tridiagonal Jacobian matrix, is a consequence of the numerical approach for the mole fraction derivatives given by eqs 10 and 11. The main advantage of the mixed approach


is that it proceeds such as a marching procedure in x1 instead of simultaneously solving the whole grid. 2. A Second-Order Approach. The mixed approach described before is computationally efficient for obtaining numerical data of the excess energy of simple ternary systems. However, because of the reasons explained in Appendix B, it should not be used when the constituent binary 1 + 3 exhibits azeotropic behavior. In addition, we have detected that when the multicomponent system exhibits azeotropes at low x2 concentrations, then the mixed approach may fail in yielding the solution. The numerical instability of the mixed approach in azeotropic systems is attributed to the use of a simple finite difference for approximating the x1 derivatives in eq 10. In this section we develop a more robust approach that, on the one hand, enhances the numerical stability of the integration of azeotropic systems but, on the other hand, implies a major computational effort. As shown in refs 3, 5, and 6, the use of central differences for approximating mole fraction derivatives proceeds generally without numerical difficulties for every type of binary system. A possible exception may be the case of binary systems that exhibit positive deviation azeotropes, where the curvature of the bubble curve about the azeotropic concentration is constrained.6 Such a constraint is not always satisfied by VLE data, because of experimental error or because of the smoothing of the bubble curve with the mole fraction xi, using empirical functions. In the second-order approach, the numerical approximation of the x1 derivative in eq 10 is replaced by the central difference

( ) ∂gE ∂x1

x2

E

≈

E

g (x1 + R) - g (x1 - R) 2R

(22)

The numerical approximation of x2 derivatives remains as in eq 11, R is defined by eq 12, and mole fractions are parametrized according to eq 13. For a ternary system, we now get the following discrete representation of Barker’s equation:

Pj,i )

[

{


]

[

E E E E gj,i+1 - gj,i-1 gj+1,i - gj-1,i + ψx1,j,i - R(j - 1) + 2R 2R 0 R(j - 1)P2,j,i E exp gj,i - R(i - 1) ψx2,j,i + Φ2,j,i

[

]}

]

{

E E gj,i+1 - gj,i-1 + ψx1,j,i + (1 - R(j - 1)) R E E gj+1,i - gj-1,i + ψx2,j,i + 2R 0 (1 - R(i - 1) - R(j - 1))P3,j,i E exp gj,i - R(i - 1) Φ3,j,i

[

[

]

]} { [

E E E E gj,i+1 - gj,i-1 gj+1,i - gj-1,i + ψx1,j,i - R(j - 1) + 2R 2R

ψx2,j,i

]}

(23)

The mole fractions of the vapor phase are calculated from

y1,j,i )

[

{ ]

0 R(i - 1)P1,j,i E exp gj,i + (1 - R(i - 1)) Φ1,j,iPj,i

[


y2,j,i )

0 1)P2,j,i

R(j Φ2,j,iPj,i

]

{

[

E exp gj,i - R(i - 1)

[

ψx1,j,i + (1 - R(j - 1))

ψx2,j,i

E gj,i+1

2R

E gj,i-1

]}

+

E E gj+1,i - gj-1,i + ψx2,j,i 2R

y3,j,i ) 1 - y1,j,i - y2,j,i

]}

(24)

Equation 23 can be used to define the equilibrium function

ϑj,i )

[

{


]

[

E E E E gj,i+1 - gj,i-1 gj+1,i - gj-1,i + ψx1,j,i - R(j - 1) + 2R 2R 0 R(j - 1)P2,j,i E ψx2,j,i + exp gj,i - R(i - 1) Φ2,j,i

[

]}

]

{

E E gj,i+1 - gj,i-1 + ψx1,j,i + (1 - R(j - 1)) 2R E E gj+1,i - gj-1,i + ψx2,j,i + 2R 0 (1 - R(i - 1) - R(j - 1))P3,j,i E exp gj,i - R(i - 1) Φ3,j,i

[

[

]

]} { [


ψx2,j,i

]}

- Pj,i (25)

which must be zero at nodes j and i. When P, T, ψx1, and ψx2 are known functions of the concentration, then E E E E , gj,i+1 , gj,i-1 , gj+1,i , and the unknowns of eq 25 are gj,i E gj-1,i. Figure 3 depicts the grid applied to the whole range of mole fractions of a ternary system. The figure shows also a specific nodes j and i for this discretization scheme. Similar to Figure 2, every vertex of the triangle represents a pure component having zero excess Gibbs energy. In addition, every side of the main triangle A-B-C is a binary system, which can be discretized using eqs 11 and 22. Compared to the mixed approach, the binary 1 + 3 is the only one that requires a different calculation strategy. The discretized relations applicable to each binary are shown in Tables 1 and 2; the method suggested in ref 6 is recommended for performing the calculations. In this manner, ternary calculations are once again reduced to the nodes inside the subtriangle D-E-F. Figure 3 indicates that the solution of nodes j and i requires the solution of columns i - 1 and i + 1 and the files j + 1 and j - 1. No satisfactory and stable equation tearing procedure was found to reduce the dimension of the calculations for each node. For this reason, we recommend to solve simultaneously every unknown node in the subtriangle D-E-F using a Newton itera-

2140


[gE]k+1 ) [gE]k + [∆gE]k

(28)

where k is an iteration index, [∆gE] is the vector that solves the linear system

[J]k[∆gE]k ) -[ϑ]k

Figure 3. Complete set of nodal points for the application of a second-order approach, eq 25, to a ternary system: (O) ternary nodes; (0) binary nodes; (9) pure-component nodes. Table 2. Discretized Relations for the Binary System 1 + 3 in the Second-Order Approach Binary 1 + 3 discretization

ϑ1,i )

[

{

0 R(i - 1)P11,i E exp g1,i + (1 - R(i - 1)) Φ11,i

E g1,i+1

2R

E g1,i-1

{

]} [

+ ψx1,1,i

E g1,i - R(i - 1)

and J corresponds to the Jacobian matrix of eq 26 evaluated at [gE]k. Because of the dimension of the nonlinear system, the dimension of the Jacobian matrix is (N2 - 5N + 6)/2 × (N2 - 5N + 6)/2. For a discretization of N ) 101, the number of matrix elements of the Jacobian has an order of magnitude of 23 × 106, which cannot be handled by a personal computer. Figure 4 shows the typical characteristic of the Jacobian matrix in eq 29, and the nonzero elements of the Jacobian for N ) 7. As seen, the Jacobian is a regular sparse nonsymmetric matrix, a characteristic that may be used for solving the linear system given by eq 29. Appendix C presents the algorithm for building the Jacobian elements. It is recommended to solve the linear system in eq 29 using the routine DLSLXG (available from IMSL8). 3. Limiting Conditions for Ternary Azeotropes. Similar to those for binary systems, the curvature and shape of the bubble surface about ternary azeotropic concentrations is constrained by stability considerations. This section is devoted to treating these constraints in detail. According to Malesinski,9 for a ternary system 2

0 1))P31,i

(1 - R(i + Φ31,i

(29)

]}

E E gj,i+1 - gj,i-1 + ψx1,1,i 2R

exp

()

i j

∂P

∂xj

(yi - xi)G ˜ xL x ∑ i)1

)

yi∆V hi ∑ i)1

- P1,j

range 2eieN-1 boundary conditions gE1,1 ) 0 E )0 g1,N

tive scheme. For a given discretization N, such an approach produces a system of N - 3 × (N - 2)/2 nonlinear equations (similar to eq 25, applied to the same number of nodes), for which it is necessary to E . determine the value of gj,i Let us now consider the following vector arrangement for the functions

2

()

T

∂T

∂xj

)-

(yi - xi)G ˜ xL x ∑ i)1

i j

2

(31)

yi∆H hi ∑ i)1

T,xi*j,3

In these equations ∆V h i is the partial volume of vaporization (i.e., V h Vi - V h Li ) and ∆H ˜ i the partial heat of h Li ). G ˜ xLi,xj represents the second vaporization (i.e., H h Vi - H derivative of the Gibbs energy for the liquid phase at constant temperature and pressure. Because azeotropes correspond to stationary points (maximum, minimum, or saddle points) of the bubble surface, the derivatives in eqs 30 and 31 are zero at the azeotropic point and, consequently, the liquid and vapor phases have the same composition. In addition, as required by intrinsic stability considerations, the Hessian matrix of the Gibbs energy for the liquid phase

and the following vector arrangement for the unknowns L

H(G ˜ ))

According to Newton’s iterative scheme, we have

(30)

2

T,xi*j,3

[

G ˜ xL1x1 G ˜ xL1x2 ˜ xL2,x2 G ˜ xL1x2 G

]

(32)

must be positive definite when evaluated at the conditions of a ternary azeotrope. Equations 30 and 31 applied to an azeotropic point yield


[

][ ]

[ ][ [ ][

] ]

∆H h 1/T v11 v12 ∆H h 2/T v21 v22 QP ) ∆H h 3/T v31 v32 matrices:

-1

l11 l12 l21 l22 l31 l32

(35)

v11 v12 y1 - 1 y2 v v y2 - 1 × H(G ˜ L) V ) 21 22 ) y1 y1 v31 v32 y2 l11 l12 x 1 - 1 x2 x2 - 1 × H(G ˜ V) L ) l21 l22 ) x1 x1 l31 l32 x2 Figure 4. Sparse Jacobian matrix for the second-order approach, showing the position of nonzero elements (x).

[ {( ) }] ∂

∂PAz

∂xk

∂xj

T,xi*j,3

) PxAzk,xj )

T,xk*j,3 2

∑ i)1

([ ] ∂yi

Az T,xk*j,3

∂xk

-

)

2

[ {( ) }] ∂

∂xk

∂T

(33a)

yAz ˜i ∑ i ∆V i)1

P,xi*j,3

2

T -

∑ i)1

([ ] ∂yi

Az P,xk*j,3

∂xk

)

(33b)

∂yi ∂xk

P,xk*j,3

˜ yV1y2 - G ˜ yV1y1 G ˜ yV2y2}]/{∆V ˜ × Det[H(G ˜ V)]} (38b) 2G ˜ xL1x2 G

G ˜ yV1y1)G ˜ yV2y2} - G ˜ yV1y2{(G ˜ xL1x2)2 + G ˜ xL1x1 G ˜ xL2x2}]/{∆V ˜ × Det [H(G ˜ V)]} (38c) where Det[H(G ˜ V)] corresponds to the determinant of the

˜ xL1x2)2G ˜ yV1y1 - 2G ˜ xL1x1 G ˜ xL1x2 G ˜ yV1y2 + G ˜ xL1x1 TxAz1x1 ) -TAz[(G {(G ˜ yV1y2)2 + G ˜ yV2y2(G ˜ xL1x1 - G ˜ yV1y1)}]/{∆H ˜ × Det[H(G ˜ V)]} (39a)

T ) - eˆ 3,k+1 QTeˆ 2,i

TxAz2x2 ) -TAz[(G ˜ xL2x2)2G ˜ yV1y1 + (G ˜ xL1x2)2 G ˜ yV2y2 + G ˜ xL2x2 {(G ˜ yV1y2)2 - 2G ˜ xL1x2 G ˜ yV1y2 - G ˜ yV1y1 G ˜ yV2y2}]/{∆H ˜ × Det[H

T ) - eˆ 3,k+1 QPeˆ 2,i

k ) 1, 2; i ) 1, 2

(G ˜ V)]} (39b) (34)

where eˆ m,l is a Cartesian basis in the mth space with a T ) unitary element at the lth position (for example, eˆ 3,2 [0 1 0]). In addition, QT and QP are matrices defined as

[

Replacing eqs 34-37 in eq 33a,b yields, after some algebra, the following pressure and temperature derivatives of the bubble curve at the azeotropic condition:

PxAz1x2 ) [G ˜ xL1x2{G ˜ xL2x2 G ˜ yV1y1 + (G ˜ yV1y2)2 + ( G ˜ xL1x1 -

where δik is the Kro¨necker operator. Parts a and b of eq 33 establish a clear connection between the curvature of the bubble surface and the curvature of the Gibbs energy at the azeotrope composition. To the best of our knowledge, no expression has been deduced for the derivatives of vapor-phase mole fractions as a function of the liquid-phase concentration in ternary systems, at the conditions given by eq 33a,b. Using the displacement method,9 we get

T,xk*j,3

(37)

PxAz2x2 ) [(G ˜ xL2x2)2G ˜ yV1y1 + (G ˜ xL1x2)2 G ˜ yV2y2 + G ˜ xL2x2{(G ˜ yV1y2)2 -

yAz hi ∑ i ∆H i)1

[ ] [ ]

]

(38a)

- δik G ˜ xLi,xj

2

∂yi ∂xk

G ˜ yV1y2 G ˜ yV1y2

{(G ˜ yV1y2)2 + G ˜ yV2y2(G ˜ xL1x1 - G ˜ yV1y1)}]/{∆V ˜ × Det[H(G ˜ V)]}

) TxAzk,xj ) P,xk*j,3

[

G ˜ yV1y1 G ˜ yV1y2

PxAz1x1 ) [(G ˜ xL1x2)2G ˜ yV1y1 - 2G ˜ xL1x1 G ˜ xL1x2 G ˜ yV1y2 + G ˜ xL1x1

Az

∂xj

In eq 36, H(G ˜ V) stands for the Hessian matrix of the Gibbs energy of the vapor phase, evaluated at constant temperature and pressure

H(G ˜ V) )

δik G ˜ xLi,xj

(36)

∆V h 1 v11 v12 ∆V h 2 v21 v22 QT ) ∆V h 3 v31 v32

][ ] -1

l11 l12 l21 l22 l31 l32

where lij and vij are the elements of the following

˜ xL1x2G ˜ xL2x2 G ˜ yV1y1 + (G ˜ yV1y2)2 + (G ˜ xL1x1 TxAz1x2 ) -TAz[{G G ˜ yV1y1)G ˜ yV2y2} - G ˜ yV1y2{(G ˜ xL1x2)2 + G ˜ xL1x1 G ˜ xL2x2}]/{∆H ˜ × Det [H(G ˜ V)]} (39c) matrix defined in eq 37 (which must be positive as required by intrinsic stability). In addition, ∆V ˜ and ∆H ˜ are the vaporization volume and enthalpy, respectively, evaluated at the conditions of the azeotrope. Inspection of eqs 38 and 39 indicates that (a) eqs 38a and 39a ˜ xL1x2, (b) eqs 38b and 39b depend depend on G ˜ xL1x1 and G

2142


on G ˜ xL2x2 and G ˜ xL1x2, and (c) eqs 38c and 39c depend on L L ˜ x2x2, and G ˜ xL1x2. Hence, on the basis of the considG ˜ x1x1, G erations developed before, we can suggest the following strategy for solving eq 38a-c (or eq 39a-c). Assume that the azeotropic coordinates, the vaporization volume ∆V ˜ (or the vaporization enthalpy ∆H ˜ ), the elements of the Hessian matrix of the Gibbs energy for the vapor phase, and the curvature of the bubble curve about the azeotrope are known. Equation 38a (or eq 39a) may be solved in G ˜ xL1x1, yielding the following solutions that depend on G ˜ xL1x2:

G ˜ xL1x1 )

[

]

∆H ˜ Tx1x1 Tx1x2 + H(G ˜ V) Az Tx x Tx x T 1 2 2 2

-4

(43)

is positive definite for the isobaric ternary azeotrope. For every other case, the systems in eqs 38 and 39 present only complex solutions, indicating that there is no real Gibbs energy able to represent the azeotrope. For the Lewis-Randall reference state, the Gibbs energy of a liquid mixture is given by

G ˜L RT

3

E

(T,P) ) g (T,P) +

3

G ˜ Lk

∑ xk ln xk + k)1 ∑ xkRT(T,P)

(44)

k)1

eq 38a

{2G ˜ xL1x2 G ˜ yV1y2 - Det[H(G ˜ V)] ( xDet[H(G ˜ V)]

xG˜

V ˜ yV1y1 y2y2(G

+ 4∆V ˜ PxAz1x1) - (2G ˜ xL1x2 - G ˜ yV1y2)2}/{2G ˜ yV2y2} (40a)

G ˜ xL1x1 )

{ x ( eq 39a

2G ˜ xL1x2 G ˜ yV1y2 - Det[H(G ˜ V)] ( xDet[H(G ˜ V)]

G ˜ yV2y2 G ˜ yV1y1 - 4

)

}

∆H ˜ Az T - (2G ˜ xL1x2 - G ˜ yV1y2)2 / Az x1x1 T {2G ˜ yV2y2} (40b)

Similarly, eq 38b (or eq 39b) may be solved in G ˜ xL2x2, yielding the following solutions that depend also on G ˜ xL1x2:

G ˜ xL2x2 )

{2G ˜ xL1x2 G ˜ yV1y2 - Det[H(G ˜ V)] ( xDet[H(G ˜ V)] V ˜ yV2y2 y1y1(G

+ 4∆V ˜ PxAz2x2) - (2G ˜ xL1x2 - G ˜ yV1y2)2}/{2G ˜ yV1y1} (41a)

G ˜ xL2x2 )

eq 39b

2G ˜ xL1x2 G ˜ yV1y2 - Det[H(G ˜ V)] ( xDet[H(G ˜ V)]

x (

G ˜ yV1y1 G ˜ yV2y2 - 4

)

∆H ˜ Az T - (2G ˜ xL1x2 - G ˜ yV1y2)2/2G ˜ yV1y1 Az x2x2 T (41b)

Both solutions (eqs 40a and 41a or eqs 40b and 41b) may be substituted in eq 38c (or eq 39c), yielding a complex function on the single variable G ˜ xL1x2. Because of its complexity, we do not present it here; however, it may be treated and analytically solved with Mathematica 3.0 using the code presented in Appendix D. In addition, a study of the solution reveals that G ˜ xL1x2 yields real values by the strategy discussed here if the matrix defined as

[

]

Px x Px x ˜ V) 4∆V ˜ P 1 1 P 1 2 + H(G x1x2 x2x2

4∆V ˜

∂2P +G ˜ yV1y1 > 0 ∂x12

(45)

while the temperature curvature of isobaric azeotropes is constrained by

eq 38b

xG˜

where G ˜ Lk is the Gibbs energy of the pure component in the liquid phase at the same temperature and pressure of the mixture. From eq 44 it follows that if a given azeotropic condition does not produce real values for the liquid-phase Gibbs energy, then it will also be unable to produce real values for the excess Gibbs energy. Consequently, the model-free approaches discussed here cannot yield a solution for a ternary system, because no solution exists in the vicinity of azeotropic-like concentrations that violate the constraint discussed for the matrices in eqs 42 and 43. It is interesting to compare the similarity of eqs 42 and 43, valid for ternary azeotropes, with the constraining relations for binary azeotropes. The following relation constrains the pressure curvature of a binary isothermal azeotrope6

(42)

is positive definite (a useful and equivalent condition is that all of the eigenvalues of the matrix are positive10) for the isothermal ternary azeotrope or the matrix defined as

∆H ˜ ∂2T +G ˜ yV1y1 > 0 Az T ∂x12

-4

(46)

The Gibbs energy of the vapor phase is usually a monotonic convex function of the concentration, yielding thus positive concavity on concentration. One exception to such a behavior is those systems that exhibit gasgas immiscibility (systems type III in the topological classification of van Konynenburg and Scott11); however, this equilibrium phenomenon has not been observed in the subcritical range. Consequently, eqs 45 and 46 may be violated essentially by binary azeotropes that present positive deviation from ideal behavior and for which the pressure curvature is negative (maximum pressure azeotropes) and the temperature curvature is positive (minimum temperature azeotropes). Similar conclusions may be formulated for the constraining relations of ternary azeotropes (eqs 42 and 43). The Hessian matrix of the Gibbs energy for the vapor phase is generally positive definite for systems at low pressures, as follows from intrinsic stability considerations and from the fact that vapors are completely miscible in the subcritical range. Consequently, the matrices defined in eqs 42 and 43 could lose positive definition if the Hessian matrix of pressure derivatives is not positive definite or if the Hessian matrix of temperature derivatives is not negative definite, as happens with positive deviation and saddle azeotropes. Although the matrices in eqs 42 and 43 must be positive definite for physical ternary azeo-


tropes, the inherent errors of experimental determinations, or the fit of bubble surfaces by means of empirical functions, may violate the constraint for ternary positive and saddle azeotropes. In fact, inspection of the definition of the matrices in eqs 42 and 43 constitutes an excellent way for testing smoothing functions of ternary data about azeotropic concentrations. The satisfaction of curvature constraints is a previous requirement for the application of the model-free approaches discussed here. In general, we have observed that the topology and concentration of smoothed ternary azeotropes (which violate the curvature constraints) are difficult to predict from binary data when using standard excess models. Example In this section we show how the model-free approaches described above apply to ternary systems and the information that can be obtained from them. More advanced applications, including the treatment of experimental ternary data, will be considered in the second part of this work. Let us consider a fictitious ternary system whose vapor pressure is represented by the following empirical relation:

P/kPa ) 100x1 + 105x2 + 100x3 + 10x1x2 + 5x1x3 + 10x2x3 (47) Equation 47 may correspond to a smoothing of experimental data. It reduces to the pure-component vapor pressure of component i when xi ) 1 and to the vapor pressures of the constituent binaries when xi ) 0 (i ) 1-3). The contour plot of isobars for eq 47 is shown in Figure 5a, where it is seen that the system has binary azeotropes for the systems 1 + 2 (B) and 2 + 3 (A). No ternary azeotrope is observed. For integration purposes, fugacity corrections will be neglected (Φi ) 1). In addition, for isothermal systems, eq 8 becomes

ψxi ) -

( )

V ˜ E ∂P RT ∂xi

≈0

(48)

xj*i

which may be reasonably neglected because the ratio V ˜ E/RT is almost zero for many systems. As pointed out before, for azeotropic systems the second-order approach is preferable over the mixed approach, especially when the system 1 + 3 is azeotropic. In our example we will use both approaches so as to discuss their features when treating azeotropic ternaries. As in every numerical integration of partial differential equations, the first problem to solve is to establish the size of the discretization grid (i.e., the N value for eq 12). This problem has not been considered in previous studies about model-free approaches. Our results with ternary systems indicate that the optimal size of the discretization grid is sensitive to the nonidealities mixture and that it should not be selected arbitrarily. It is useful to explore the convergence of the method to the properties at the equimolar concentration of a ternary (x1 ) x2 ) x3 ) 1/3). The recommended control variable is the vapor-phase concentration, for which the experimental accuracy is known or can be estimated. The integration method should exhibit convergence errors in the range of the experimental determinations (i.e., less than 10-3). In addition and considering that pure-component vapor pressures and the vapor pressure of the system are given, then the

convergence of the concentration for the vapor phase also warrants the convergence of activity coefficients. Figure 5b shows the convergent values of vapor-phase mole fractions (within an error of 10-4), which are met for N g 250 in the case of the mixed approach and for N g 100 in the second-order approach. Such a behavior is to be expected because the mixed approach uses a less accurate approximation of the derivative (see eq 10). Figure 5b also shows that both approaches yield the same convergent values. Based on the size of the integration grid, the selection between the mixed or the second-order approaches is a tradeoff problem: the mixed approach runs faster than the second-order approach in the ratio 3:1 to 6:1. For N g 150, the second-order approach is very slow, although always stable. In fact, although the results of the mixed approach tend to the correct solution of vapor-phase mole fractions, numerical instability was observed for N > 290 because of the presence of binary azeotropes in the system. Once the optimal size of the grid and the approach have been determined, we can calculate the equilibrium properties. We will consider the results of the secondorder approach using N ) 101 (the mixed approach yields almost equivalent results for N ) 290). Figure 5c illustrates the calculated excess Gibbs energy surface, from which we deduce that the ternary system exhibits positive deviation from ideal behavior, as expected from an inspection of eq 52. Parts d and e of Figure 5 show the activity coefficients for each component, from which we can deduce the effect of the ternary mixture on the interactions that occur in the liquid phase. We observe that the values of the activity coefficients of each individual species are bounded by the activity coefficients of the binary mixture, indicating that the ternary mixture improves the affinity between components in the liquid phase. Finally, Figure 5f shows the equilibrium diagram of the ternary system including vapor-phase mole fractions. This simple application shows that model-free approaches are able to generate a complete set of equilibrium properties. Because the starting point of the approach is the Gibbs-Duhem equation, the properties so obtained constitute a set of consistent properties for the smoothed bubble surface and the approximations considered for the calculations (i.e., assumptions about vapor-phase corrections and the effect of excess volumes or enthalpies). The application of the method to isobaric systems is similar to the case of isothermal systems, although, in a isobaric case, a fit of the bubble temperature surface in liquid-phase mole fractions is required (a smoothing function similar to eq 47). Concluding Remarks In this work we have presented two numerical formulations of model-free techniques for treating ternary systems: the mixed and the second-order approaches. The techniques differ on the computational charge, and their range of applicability depends on the complexity of the ternary system to be treated. The mixed approach is recommended for nonazeotropic isothermal and isobaric systems. The second-order approach can be used in every type of miscible system, although it requires more computational effort. The model-free approaches presented here are capable of predicting consistent equilibrium properties starting from a fit of the bubble surface on liquid-phase mole

2144


Figure 5. (a) Isobars for the fictitious system predicted by eq 47. (A and B) binary azeotropic points. (b) Vapor-phase mole fraction convergence of the mixed and second-order approaches for the system predicted by eq 47. (c) Excess Gibbs energy surface calculated for the system in eq 47. Data generated using the second-order approach with N ) 101: (‚‚‚) ternary mixture; (s) binaries. (d) Activity coefficient plot for component 1. (e) Activity coefficient plot for component 2. (f) Activity coefficient plot for component 3. (g) Phase diagram calculated from the second-order approach for the system in eq 47.

fractions. These equilibrium properties, particularly the excess Gibbs energy function and vapor-phase mole

fractions, do not require the assumption of models, thus avoiding possible modeling pitfalls when treating mul-


ticomponent data. Compared to the usual treatment of data by means of GE models, the model-free approaches presented here are not efficient for performing VLE calculations. However, they are so for assessing data and for studying special features of multicomponent phase equilibrium surfaces. They are recommended for guiding the selection of models, for testing the reliability of ternary data, and for analyzing the possibility of predicting ternary data from binaries. It is important to establish the differences between the numerical approaches proposed here and that of Mixon et al.3 for ternaries. In our approach, model-free calculations are directly applicable to isothermal and isobaric data. However, one limitation may be the availability of data for excess volumes or enthalpies. From a numerical viewpoint, Mixon et al. proposed a point relaxation technique for solving the discrete equation of Barker. According to Mixon et al., such a strategy is recommended because of the difficulties in programming a block relaxation over the triangular concentration diagram. However, our approaches overcome this difficulty. When block relaxation calculations are performed, all of the nodes of the triangular grid are solved simultaneously, inducing a better numerical performance of the Newton method. According to the point relaxation approach, calculations are performed iteratively node per node in the discretization grid, a condition which may produce numerical instability when the problem requires thin grids. In addition, we have demonstrated that the values of the calculated equilibrium properties are sensitive to the size of the grid, whose optimal value may be determined by exploring the convergence of calculated properties. The time of calculation for the point relaxation technique increases substantially as the size of the grid increases, an inconvenient situation when treating complex ternary systems. Finally, the displacement theory was used to determine physical constraints for the curvature of the bubble surface in azeotropic concentration. These constraints affect specifically saddle and positive deviation azeotropes. No solution for the GE function may be achieved by means of model-free techniques when fitted or experimental azeotropes violate the constraints. Acknowledgment This work was partially financed by FONDECYT, Santiago, Chile (Projects 1990402 and 7990065), and by MEC, Valencia, Spain (Project PB96-0788). E.L. was sponsored by CONICYT, Chile. A.M. acknowledges a grant from SF, Concepioń, Chile. Appendix A: A Numerical Solution for Equation 18 As pointed out in the text, eq 18

ϑ1,i )

[

{

0 R(i - 1)P11,i E exp g1,i + (1 - R(i - 1)) Φ11,i

]} [

E E g1,i - g1,i-1 + ψx1,1,i R

+

{

0 (1 - R(i - 1))P31,i E exp g1,i Φ31,i

E E g1,i - g1,i-1 R(i - 1) + ψx1,1,i R

]}

- P1,i (A.1)

represents the binary system 1 + 3. It may be solved

E E provided that g1,i-1 is known. on the single variable g1,i Equation A.1 has two known values in the limits of pure E ) gE(x1)0) ) 0 and gE1,N ) gE(x1)1) ) components: g1,1 0. Therefore, when 2 e i e N - 1 is selected, it is possible to obtain gE values for midrange concentrations. Because eq A.1 is nonlinear, the Newton method is recommended. In such an iterative approach

E E [g1,i ]k+1 ) [g1,i ]k -

[

ϑ1,i

E dϑ1,i/dg1,i

]

(A.2)

k

where k is an iteration index, and the required derivative is given by

dϑ1,i E dg1,i

) [2 - i][ϑ1,i + P1,i] +

[

{

0 (i - 1)P11,i E exp g1,i + Φ11,i

]}

E E g1,i - g1,i-1 (1 - R(i - 1)) + ψx1,1,i R

(A.3)

The fugacity correction Φ is an implicit variable in eqs A.1-A.3 because it depends on vapor-phase mole fractions that, in turn, may be calculated from

y11,i )

{

0 R(i - 1)P11,i E exp g1,i + (1 - R(i - 1)) Φ11,iP1,i

[

E E g1,i - g1,i-1 + ψx1,1,i R

y31,i ) 1 - y11,i

]}

(A.4)

once the iterative scheme in eq A.2 converges. The recommended procedure for solving the binary 1 + 3 is the following: 1. Assume ideal gas behavior, for which Φl,1,i ) 1 (l ) 1 and 3). 2. Let i ) 2. 3. Solve eq A.1 using the iterative scheme indicated in eq A.2. 4. Calculate vapor-phase mole fractions as indicated in eq A.4. 5. Recalculate Φl,1,i ) 1 (l ) 1 and 3), and repeat from step 3 until no change is observed in vapor-phase mole fractions. Once step 5 converges, we obtain a numerical value for g1,Ei . 6. Let i ) i + 1, and repeat from step 3 until i ) N 2. When starting with i ) 2, eq A.3 is positive, independently of the thermodynamic behavior of the binary system 1 + 3. In addition, inspection of eq A.1 reveals that the sum ϑ1,i + P1,i is always positive. Consequently, as i increases above 2, eq A.3 may become negative and eventually be zero if a change of sign occurs. Such a behavior is undesirable for the Newton method, which may become numerically indefinite as follows from eq A.2, inducing a pitfall in the numerical procedure described here. Particularly, such a situation constrains the application of the mixed approach to cases in which eq A.3 remains positive, without exhibiting a change in sign. Let us now develop the criteria for establishing the requirements that allow application of the mixed approach. As follows from the physical significance of each of its terms, eq A.3 may be written as

2146

dϑ1,i


[

) (2 - i) E

dg1,i

]

Newton algorithm, according to which the set of unknowns is obtained from the iterative procedure

0 0 (1 - x1)γ3P31,i x1γ11,iP11,i + + Φ11,i Φ21,i 0 (i - 1)P11,i γ11,i (A.5) Φ11,i

Equation A.5 will remain positive, as required by the mixed approach, if 0 Φ11,i i - 1 γ31,iP31,i < x1 + (1 - x1) 0 γ11,iP11,iΦ21,i i - 2

∂[ϑ2,i‚‚‚ϑj,i‚‚‚ϑN-i,i] E E E ∂[g2,i ‚‚‚gj,i ‚‚‚gN-i,i ]

1 1

E ∆g2,i E ∆g3,i l ; ∆gEk ) E ∆gj,i l E ∆gN-i,i k

(A.7)

(A.8)

2

Thus, eq A.7 becomes

x1 + (1 - x1)R3,1 < 1

(A.9)

Equation A.9 is always satisfied if R3,1 < 1, implying that component 3 is the less volatile component. Therefore, the convergence of the mixed approach is assured provided that P3 < P1 and the binary system 1 + 3 is not azeotropic. Appendix B: A Numerical Solution for the Nonlinear System in Equation 20 Equation 19 is an abbreviation of eq 16, which represents the equilibrium condition for a discretization node. For a ternary system this equilibrium condition is expressed using Barker’s relation. Considering the numerical approximation of the derivatives of the gE function, to obtain activity coefficients, the ϑj,i function exhibits a nonlinear dependence on excess energy as follows:

ϑj,i )

E E E ϑj,i{gj-1,i ,gj,i ,gj+1,i}

≡0

(B.1)

When eq B.1 is applied to every node in the ith column of Figure 2, it is possible to obtain the system given by eq 20 E E E ,g2,i ,g3,i }≡0 ϑ2,i ) ϑ2,i{g1,i E E E ,g3,i ,g4,i }≡0 ϑ3,i ) ϑ3,i{g2,i E E E ,gN-i,i ,gN-i+1,i } ≡ 0 (B.2) ϑN-i,i ) ϑN-i,i{gN-i-1,i E E where g1,i and gN-i+1,i correspond to known excess Gibbs energies of the binary systems 1 + 3 and 1 + 2, respectively. Because the system in eq B.2 is nonlinear (see eq 16), an adequate candidate for its solution is the

[

|k∆gEk ) -Tk

(B.4)

[] []

where

In addition, in eq A.7 we can recognize the relative volatility R3,1 of component 3 in component 1 defined as

y3 x1 γ3P03Φ1 R3,1 ) ) x3 y1 γ P0Φ

j ) 2, N - i (B.3)

where k is an iteration index. In addition, the iterative E in eq B.3 are obtained from the Jacocorrections ∆gj,i bian matrix of the system in eq B.2 by solving

(A.6)

For a given discretization N, the ratio (i - 1)/(i - 2) yields the minimum value (N - 1)/(N - 2), which is approximately unitary if N is large. Therefore, eq A.6 is always satisfied if 0 γ31,iP31,i Φ11,i N - j, let i ) 2. 9. If j e N - 2, let j ) j + 1 and repeat from step 4. The derivatives needed in steps 4 and 5 may be easily calculated from eq 25. Appendix D. Mathematica 3.0 Codes for Solving Analytically Eqs 38c and 39c (1) Isothermal case

]}

0 (1 - R(i - 1) - R(j - 1))P3,j,i E Ω3 ) exp gj,i - R(i - 1) Φ3,j,i

[

]

[

]}

E E E E gj,i - gj,i-1 gj+1,i - gj-1,i + ψx1,j,i - R(j - 1) + ψx2,j,i R 2R (B.8)

The procedure is the following: 1. Assume ideal gas behavior, for which Φl,j,i ) 1 (l ) 1 and 3). E E E E 2. Give an initial estimate {g2,i , g3,i ‚‚‚gj,i ‚‚‚gN-i,i }. 3. Evaluate {ϑ2,i, ϑ3,i‚‚‚ϑj,i‚‚‚ϑN-i,i} using eq 16. 4. Evaluate the elements of the Jacobian matrix using eq B.7. 5. Solve the tridiagonal system in eq B.6. E estimate according to eq B.3. 6. Actualize the gj,i 7. Proceed from step 3 until ||∆gEk || e , where || || denotes a Euclidian norm and is the desired convergence error. 8. Calculate the vapor-phase mole fractions from eq 15. 9. With these vapor-phase mole fractions, recalculate the vapor-phase corrections Φl,j,i (l ) 1 and 3). 10. Proceed from step 3 until no change is observed in the calculated vapor-phase mole fractions. Appendix C: Computation of Nonzero Elements for the Jacobian Matrix in the Second-Order Approach In this appendix we present the algorithm for building the Jacobian matrix needed in eq 29. Pertinent elements of the matrix in question are called Jl,m: 1. Let k ) 1. 2. Let j ) 2. 3. Let i ) 2. 4. Let Jk,k ) ∂ϑj,i/∂gj,i. 5. If i < N - j, let

(2) Isobaric case

List of Symbols C ) number of components g ) dimensionless Gibbs energy (G ˜ /RT) G ) Gibbs energy H ) enthalpy N ) grid size n ) number of moles P ) absolute pressure R ) universal gas constant T ) absolute temperature V ) volume x, y ) compositions of the liquid and vapor phases Greek Letters R ) parameter defined in eq 12 ∆ ) vaporization property γ ) activity coefficient Φ ) fugacity correction, defined in eq 2 θ ) weighting parameter ψ ) function defined in eq 8 ϑ ) equilibrium function defined in eqs 16 and 25 φˆ ) fugacity coefficient ∂ ) partial derivative Subscripts P, T ) at constant pressure and temperature

2148


Superscripts ∼ ) molar property - ) partial property Az ) azeotropic composition 0 ) pure-component saturation property E ) excess property exptal ) experimental value calc ) calculated value L ) liquid phase V ) vapor phase T ) transposition operator

Literature Cited (1) Barker, J. A. Determination of Activity Coefficients from Total Pressure Measurements. Aust. J. Chem. 1953, 6, 207-210. (2) Van Ness, H. C.; Abbott, M. M. Classical Thermodynamics of Nonelectrolyte Solutions With Applications to Phase Equilibria; McGraw-Hill: New York, 1982. (3) Mixon, F. O.; Gumowski, B.; Carpenter, B. H. Computation of Vapor-Liquid Equilibrium Data from Solution Vapor Pressure Measurements. Ind. Eng. Chem. Fundam. 1965, 4, 455-459. (4) Van Ness, H. C. On the Integration of the Coexistence Equation for Binary Vapor-Liquid Equilibrium. AIChE J. 1970, 16, 18-22.

(5) Sayegh, S. G.; Vera, J. H. Model-Free Methods for Vapor Liquid Equilibria Calculations. Binary Systems. Chem. Eng. Sci. 1980, 35, 2247-2256. (6) Wisniak, J.; Apelblat, A.; Segura, H. Application of ModelFree Computation Techniques and Limiting Conditions for Azeotropy. An Additional Assessment of Experimental Data. Chem. Eng. Sci. 1997, 52, 4393-4402. (7) Carnahan, B.; Luther, H. A.; Wilkes, J. O. Applied Numerical Methods; Wiley: New York, 1969. (8) IMSL Library Reference Manual, 9th ed.; IMSL: Houston, TX, 1982. (9) Malesinski, W. Azeotropy and Other Theoretical Problems of Vapour-Liquid Equilibrium; Interscience Publishers: London, 1965. (10) Edgar, T. F.; Himmelblau, D. M. Optimization of Chemical Processes; McGraw-Hill: New York, 1988. (11) van Konynenburg, P.; Scott, R. L. Critical Lines and Phase Equilibria in Binary van der Waals Mixtures. Philos. Trans. R. Soc. London 1980, 298A, 495-539.

Received for review August 15, 2000 Revised manuscript received January 4, 2001 Accepted January 10, 2001 IE0007522

A Model-Free Approach Data Treatment of Vapor ... - ACS Publications

Recommend Documents