Generalized Procedure for Estimating the Fractions of Nonbonded

Farshad Tabasinejad , R. Gordon Moore , Sudarshan A. Mehta , Kees C. Van Fraassen ... S. B. Kiselev and J. F. Ely , S. P. Tan, H. Adidharma, and M. Ra...
0 downloads 0 Views 68KB Size
Ind. Eng. Chem. Res. 2004, 43, 203-208

203

GENERAL RESEARCH Generalized Procedure for Estimating the Fractions of Nonbonded Associating Molecules and Their Derivatives in Thermodynamic Perturbation Theory Sugata Pikatan Tan, Hertanto Adidharma, and Maciej Radosz* Department of Chemical and Petroleum Engineering, University of Wyoming, Laramie, Wyoming 82071-3295

A generalized calculation procedure for solving thermodynamic perturbation theory (TPT1) is found to give reliable numerical estimates of the fraction of nonbonded associating molecules, as well as their first- and second-order partial derivatives with respect to temperature, density, and mole fractions. This is in contrast to the conventional procedures that rely on specialized equations depending on the number and type of association sites and the number of components in the mixture. This procedure is applicable to SAFT-like equations of state and other models that utilize the TPT1 association term. Introduction (TPT1)1-4

Wertheim’s first-order perturbation theory underpins applied thermodynamic models for associating molecules. For example, all SAFT-like equations of state, such as that proposed by Chapman et al.,5 utilize an association term based on this theory. Mu¨ller and Gubbins6 have reviewed other versions and applications of the SAFT equation of state, so we do not repeat such a review here, except to say that the SAFT model has become a popular engineering tool. It is also worth noting that the significance of the TPT1 association term goes beyond SAFT-like equations of state. For example, this term has been demonstrated to work with cubic equations of state, such as SRK,7 and activity-coefficient models, such as NRTL,8 UNIQUAC,9 and UNIFAC.10 Even though the TPT1 association term is a simple but realistic approximation, it is hard to implement in process simulators because it calls for mass-action equations that have multiple analytical forms to estimate the fractions of nonbonded molecules. These equations depend on the association site assignment for each associating component.11 This means that, each time one changes a set of association sites, for example, in going from system to system or from one association concept to another, one must change the equations beforehand. The situation is actually much worse than it initially sounds because, in most engineering calculations, for example, phase equilibrium and energy balance calculations, one needs not only the fractions of molecules that are not bonded for each association site type for each component, but also their first- and secondorder partial derivatives, for example, with respect to mole fractions, temperature, and density. Therefore, if one wants to test a different set of bonding sites, which is often tempting, or to include a different set of * To whom correspondence should be addressed. E-mail: [email protected]. Tel.: 307-766-2500. Fax: 307-766-6777.

associating components, which is often needed, or both, one faces a time-consuming effort of rederiving the equations for the fractions of nonbonded molecules and their first- and second-order derivatives, which inhibits creative applications. The issue of selecting specialized association terms for special mixtures has been alluded to or discussed by others.11,12 Huang and Radosz,13 for example, provided a tabulated list of examples of such specialized association terms. The first serious attempt aimed at making TPT1 more efficient was that by Elliott,14 who simplified TPT1 for systems that can satisfy the following three criteria: each molecule must be represented by segments with only two types of association sites of the same quantity, the extent of association must be the same for all hydrogen-bonding segments on the same molecule, and the cross-association parameter must be approximated by the square-root rule. For a chemical system that roughly satisfy these criteria, Elliot14 proposed an efficient way to calculate the fractions of nonbonded molecules as well as their first-order partial derivatives. Given that many important chemical systems do not satisfy these criteria, however, a more general procedure for solving TPT1 is still needed. Michelsen and Hendriks15 showed that the calculation of the association contributions to the pressure and chemical potential from TPT1 can be simplified significantly by the optimization of a function that does not require the first derivatives of the fractions of nonbonded molecules. This approach works for generalizing the association calculation as long as the second derivatives of the Helmholtz energy are not needed. The generalization of this approach is not applicable to phase stability analysis using the tangent plane method16 and calculating heat capacity, spinodal curve, and critical point, which require the second derivatives. The purpose of this paper, therefore, is to develop a general calculation procedure for the fractions of non-

10.1021/ie034041q CCC: $27.50 © 2004 American Chemical Society Published on Web 11/22/2003

204

Ind. Eng. Chem. Res., Vol. 43, No. 1, 2004

bonded molecules and their derivatives up to the second order that is applicable to all associating systems that have proton donors and proton acceptors, regardless of the number and type of association sites and the number of components in the mixture.

association bonds, we set ljki ) 0, which makes the association strength ∆ljki ) 0 according to eq 3.

Helmholtz Energy Change Due to Association

The mass action, eq 2, can be rearranged to a form that facilitates general solutions

Because the assumptions and derivation are described elsewhere,17 we shall only cite the working equation needed to calculate the Helmholtz energy change due to association. The dimensionless Helmholtz energy due to this association is given by

Aassoc NkT

n

)

Generalized Approach to Estimating Xji and Its Derivatives

1

Xji )

n

1+F

s

∑ xk l)1,l*j ∑ Slk Xlk ∆ljki k)1

i ) 1, ..., n; j ) 1, ..., s (4)

s

∑ ∑ i)1 j)1 xi

Sji[ln(Xji) - 1/2Xji + 1/2]

(1)

where xi is the mole fraction of component i, Xji is the fraction of component (molecule) i that is not bonded at sites j, Sji is the number of association sites of type j in each molecule of component i, n is the number of components in the mixture, and s is the number of types of association sites. The fractions of component i not bonded at sites of type j, Xji, are calculated from a mass-action equation that relates the number of molecules i bonded at sites of type j subject to the probability of these nonbonded molecules to form bonds. This probability depends on the number of molecules not bonded at other site types l in any component k available to make association bonds with the sites of type j and the corresponding interaction strength between site type j in component i and site type l in component k. Note that l is never equal j, i.e., the associating sites are always of different types to reflect the concept of donor-acceptor association bonding. Thus

This form is similar to that in refs 6, 11, and 19, except that the number of a particular site type l in component k, Slk, is explicitly brought out in eqs 4. To develop a generalized approach to solving eqs 4, we have to include all possible combinations of Xji, including those where the association site type j is absent in component i. If we have s site types and n components in the mixture, there are s × n selfconsistent simultaneous eqs 4 that are solved iteratively for s × n values of Xji. Note that, if the association site type j is absent in component i, i.e., Sji ) 0 and ∆ljki ) 0, then eqs 4 give Xji ) 1, and this does not affect the Helmholtz energy calculated from eq 1. Attempts to reduce eqs 4 to specific expressions, for example, by eliminating terms containing ∆ljki ) 0 from the righthand-side denominator, thus eliminating all Xji ) 1 as unknowns, make the generalization less obvious. To solve eqs 4 for Xji, one can use a simple iterative procedure that starts with a set of approximate trial solutions (initial guesses) substituted for the right-hand side to obtain better estimates of the roots for the next iteration. This process is repeated until sufficiently accurate approximations of the roots are obtained. The iterations converge rapidly with arbitrary initial values of Xji between 0 and 1. The partial derivatives of Xji with respect to a variable ξ are calculated from eqs 4 as follows

∂Xji ∆ljki

is the associawhere F is the number density and tion strength between site type l in component k and site type j in component i, which can be approximated by18

[ ( ) ]

∆ljki ) Kljkigseg(σ) exp

ljki -1 kT

(3)

In eq 3, gseg(σ) is the segment pair correlation function evaluated at contact, Kljki is a measure of the volume available for bonding between a site of type l in component k and a site of type j in component i, and ljki is the well depth of the site-site interaction potential between a site of type l in component k and a site of type j in component i. If k and i are equal, we have a case of self-association; otherwise, we have a case of cross-association. The subscript-superscript pairs are symmetrical, that is, their sequence does not make a difference, because they correspond to the same interaction, for example, ljki ) jlik. For two sites with no

∂ξ

) -(Xji)2

n

∂ ∂ξ

(F

∑ k)1

s

Slk Xlk ∆ljki) ∑ l)1,l*j

xk

i ) 1, ..., n; j ) 1, ..., s (5)

The variable ξ can be density, temperature, or mole fractions, one at a time, while the other variables remain constant in the partial derivatives. If ξ is the mole fraction of component i (xi) in the mixture, the mole fractions of the other components, xk with k * i, remain constant as well as the temperature and density. The derivatives of Xji given by eq 5 lead to a set of linear simultaneous equations written in matrix form as

[∂ξ∂ [X ]] ) [Ψ ]

[Λpq]

j i

ξ p

(6)

where the coordinates of the elements of the matrix [Λpq] are

p ) n(j - 1) + i

i ) 1, ..., n; j ) 1, ..., s

q ) n(l - 1) + k

k ) 1, ..., n; l ) 1, ..., s

(7)

Ind. Eng. Chem. Res., Vol. 43, No. 1, 2004 205

elements are

ΨFF p

)

()

2

j 2 ∂Xi

Xji

j

2 ∂Xi

+

F ∂F

∂F

n

F

+ (Xji)2

s

∑ xk l)1,l*j ∑

Slk

k)1

ΨTT p

)

()

(

Xlk

s

∑ xk l)1,l*j ∑

Slk

k)1

Thus, the matrix [Λpq] has an order of (s × n) × (s × n), as shown in Figure 1. The matrix containing the derivatives, i.e., the second matrix on the left-hand side of eq 6, has an order of (s × n) × 1, as does the matrix [Ψξp] on the right-hand side of eq 6. The superscript of the elements of matrix [Ψξp] denotes the variable ξ, and the row p runs from 1 to (s × n). Note that the elements of matrix [Λpq] do not depend on ξ, whereas those of matrix [Ψξp] do. The elements of matrix [Λpq] are

Λpq )

{

FxkSlk ∆ljki(Xji)2 l * j

(8)

where the Kronecker delta δki equals 1 for k ) i and 0 otherwise. The elements of matrix [Ψξp] for the derivatives with respect to density, temperature, and mole fractions are

[( )

ΨFp ) -(Xji)2

1 1

F Xj i

n

-1 +F

∑ k)1

n

ΨTp

)

-(Xji)2F s

Ψxpm

)

-(Xji)2F



l)1,l*j

s

xk

∑ l)1,l*j

Slk Xlk

s

∑ xk l)1,l*j ∑

Slk

Xlk

k)1

(

Xlm

∆ljmi

+



]

∂∆ljki ∂F

∂∆ljki ∂T

n

Slm

xkSlk

Xlk

(9)

(10)

)

∂∆ljki

∂xm m ) 1, ..., n (11)

k)1

The second derivatives with respect to variables ξ and ω can be calculated in a similar fashion. We can still use a matrix equation similar to eq 6, with the same matrix [Λpq]

[

[Λpq]

]

∂2 [Xj] ) [Ψξω p ] ∂ξ∂ω i

)

( )( )

j j 2 ∂Xi ∂Xi

Xji ∂F

+

F ∂T

s

∑ xk l)1,l*j ∑

Slk

k)1

(

Xlk

( )( ) ∑[ ( ( ∑ (

Ψxpmxa )

2

∂F2

+2

∂F ∂F

∂Xji

∂2∆ljki

+2

∂T2

∂Xji

(13)

∂T ∂T

Slm

n

xkSlk

Xlk

(14)

-

∂2∆ljki ∂F∂T

+

∂Xlk ∂∆ljki

+

∂F ∂T

)

∂Xlk ∂∆ljki ∂T ∂F

(15)

-

s

(Xji)2F

Sla

)

∂Xlk ∂∆ljki

Xji ∂xm ∂xa

k)1

)]

∂Xlk ∂∆ljki

j

∂T

n

(Xji)2F

(

Xlk

1 ∂Xi

l)1,l*j

l)j

δki

ΨFT p

∂2∆ljki

-1 -

-

Xji ∂T

n

Figure 1. Structure of matrix [Λpq]; it has s n × n unit (identity) submatrices [I] as shown by shaded boxes; the matrix elements in the nonshaded region are Λpq ) FxkSlk ∆ljki(Xji)2.

F2 Xji

2

j 2 ∂Xi

(Xji)2F

[( ) 2 1

)

∂Xlm ∂∆ljmi l lj ∆mi + Xm + ∂xa ∂xa

∂Xla

∆ljai ∂xm ∂2∆ljki

∂xm∂xa

+

+

Xla

)

∂∆ljai ∂xm

∂Xlk ∂∆ljki ∂xm ∂xa

+

+

)]

∂Xlk ∂∆ljki ∂xa ∂xm

(16)

In eq 16, the subscripts m and a of the mole fractions run with m ) 1, ..., n, and a ) m, ..., n, so that there are n(1 + n)/2 combinations of m and a. For association sites that are absent, where Sji ) 0, lj ∆ki ) 0, and Xji ) 1, it is easy to show that the corresponding elements of matrix [Ψξp] and matrix [Ψξω p ] are zero, which means that the first- and second-order partial derivatives of Xji are also zero, which, in turn, means that these derivatives do not affect the corresponding partial derivatives of the residual Helmholtz energy and thus the generality of the procedure. Examples of code written in Fortran 90 are provided in Appendix A for solving eqs 4 with a simple iterative method and in Appendix B for constructing matrix [Λpq], whose elements are given by eq 8, and matrix [ΨFp]], whose elements are given by eq 9, and then solving eq 6 using matrix inversion. The values of ∆ljki and its derivatives must be provided as inputs.

(12)

By taking the derivatives of eq 5, we can determine the elements of matrix [Ψξω p ] on the right-hand side of eq 12. For some common second derivatives, these

Examples Illustrating the Use of This Procedure Figure 2 shows the molecules used in our examples, including their association sites. The hydrogen atom is

206

Ind. Eng. Chem. Res., Vol. 43, No. 1, 2004

Figure 2. (a) Water; (b) poly(ethylene glycol); (c) poly(ethylene oxide); (d) aliphatic alcohol; (e) hydrogen fluoride. Triangles are sites of type 1 (H) and diamonds are sites of type 2 (O) or 3 (F).

assumed to have one association site (type 1), the oxygen atom is assumed to have two association sites (type 2, one per each electron lone pair), and the fluorine atom in HF is assumed to have one association site (type 3). In fact, this procedure allows for hypothesizing other association-site assignments, for example, differentiating between the oxygen association sites in hydroxyl groups and ether groups. Pure Water. The self-association of pure water is due to bonding among sites shown in Figure 2a. We have two sites of type 1 and two sites of type 2 on each water molecule (S11 ) S21 ) 2). The site-site interactions that lead to bonding are characterized by two parameters, 12 12 11 and K11, which are used to estimate the association strength ∆12 11 from eq 3. Using a simple iterative procedure to solve eqs 4, we get X11 ) X21, which is equivalent to the common analytical solution for an association of type 4C in Huang and Radosz13

X11

)

-1 + x1 + 8F∆12 11 4F∆12 11

(17)

Water (1)/PEG (2) and Water (1)/PEO (2). For water (1)/polyethylene glycol (2) mixtures, the numbers of bonding sites are set to be S11 ) S21 ) S12 ) 2 and S22, which depends on the molecular weight of the PEG. As shown in Figure 2b, the repeating unit of PEG has two carbon atoms, four hydrogen atoms, and one oxygen atom, and hence, its molecular weight is 44. The degree of polymerization n is the number of repeating units. If the molecular weight of PEG is M, n ) (M - 18)/44, and S22 ) 2n + 2. All of these sites lead to the following 12 12 12 12 12 12 interaction parameters: 12 11, K11, 22, K22, 12, K12, 21, 12 and K21. The last four parameters are the crossassociation parameters, which are commonly estimated on the basis of the self-association parameters using mixing rules. X11, X12, X21, and X22 are obtained from eqs 4. The matrix [Λpq] is of order 4 × 4, and [Ψξp] is of order 4 × 1. Long PEG molecules are practically indistinguishable from poly(ethylene oxide) (PEO), which is shown in Figure 2c. For a system of water (1)/PEO (2), S12 ) 0 and S22 ) 2n, with n ) (M - 2)/44, where M is the molecular weight of PEO. Because of the absence 12 of hydrogen-type association sites in PEO, 12 22 ) 21 ) 1 0, so that we obtain X2 ) 1 and 0 for all of its derivatives. The solutions for these systems using our general numerical approach are equal to those obtained using an analytical approach derived from models II and

III proposed by Kraska,12 i.e., model II for water (1)/ PEO (2) and model III for water (1)/PEG (2). If we want to test the effect of the hydrogen association sites at both ends of PEG in water (1)/PEG (2) mixtures, we can turn 12 off the hydrogen sites by setting 12 22 ) 21 ) 0. This is equivalent to Kraska’s model II, similar to that in water (1)/PEO (2) mixture. Water (1)/PEO (2)/Alcohol (3). The ternary mixture of water (1)/PEO (2)/alcohol (3) has additional sites, S13 ) 1 and S23 ) 2, on the alcohol molecules as shown in Figure 2d, and consequently, X13 and X23. The matrix [Λpq] is therefore of order 6 × 6, and the matrix [Ψξp] is of order 6 × 1. Water (1)/HF (2). The water (1)/HF (2) mixture has an extra site of proton-acceptor type F in HF molecules as shown in Figure 2e, so that S11 ) S21 ) 2 and S12 ) S32 ) 1. Because of the absence of F in water and O in the HF molecule, we have S31 ) S22 ) 0. The interactions involving the absent sites are set equal to zero (13 11 ) 12 22 ) 0), as are the interactions among proton-acceptor sites (23 ki ) 0) for all components k and i. The matrix [Λpq] is therefore of order 6 × 6, and the matrix [Ψξp] is of order 6 × 1. This procedure results in X22 ) X31 ) 1 and values of 0 for their derivatives. If one assumes, rightly or wrongly, that the cross bonds between H-F 12 and H-O have the same bond strength, i.e., ∆13 12 ) ∆21, 1 2 1 3 then solving eqs 4, one obtains X1 ) X1 and X2 ) X2. In this case, eqs 4 reduces to those specialized eqs 15 and 16 of Galindo et al.20 To illustrate the CPU time needed, we calculate Xji and their derivatives of the first and second order with respect to density using double precision, with Compaq Visual FORTRAN 6.60 and Intel Pentium 4 CPU (1.80 GHz). For each example above, we need a CPU time of less than 10-4 s, most of which is used to calculate Xji. Application to Other, Non-Donor-Acceptor, Types of Association. If an association site is allowed to interact with sites of its own type, the restriction of l * j is relaxed in the mass-action equation, which leads to associations that are different from donor-acceptor associations, for example, hydrophobic associations among protein molecules.21 Our procedure is applicable to such types of association as follows: First, we remove the l * j restriction from the massaction equation and its derivatives. Second, we replace eq 8 by

Λpq ) FxkSlk ∆ljki(Xji)2 + δpq

(18)

where the Kronecker delta δpq equals 1 for p ) q and 0 otherwise. Conclusion A generalized procedure for solving TPT1 is found to give fast and reliable numerical solutions for the fraction of nonbonded molecules and their first- and second-order partial derivatives regardless of the number and type of association sites and the number of components in the mixture. This will facilitate the implementation of TPT1 for broader applications and its use in process simulators.

Ind. Eng. Chem. Res., Vol. 43, No. 1, 2004 207

Acknowledgment This work was funded by NSF Grant CTS-0244388. Nomenclature Roman Letters Aassoc ) Helmholtz energy due to association gseg ) segment pair correlation function k ) Boltzmann constant Kljki ) measure of the volume available for bonding between site l in component k and site j in component i n ) number of components; number of repeating units in Figure 2 N ) number of molecules p, q ) coordinates of the elements of matrix Λ s ) number of association-site types Sji ) number of association sites of type j in a molecule of component i T ) absolute temperature xi ) molar fraction of component i Xji ) molar fraction of component i not bonded at association sites of type j Greek Letters δkl ) Kronecker delta () 1 for k ) l and ) 0 for k * l) ∆ljki ) association strength between site l in molecule k and site j in molecule i ljki ) association energy between site l in molecule k and site j in molecule i Λ ) coefficient matrix on the left-hand side of derivative equations F ) number density σ ) contact distance between molecules ω , ξ ) variables (T, F, or xi) with respect to which the derivatives are taken Ψ ) matrix on the right-hand side of derivative equations

Appendix A. Fortran Code for Solving Eqs 4 Using Simple Iterative Method ! Initialization X_A ) 0.5d0 del ) 1.0d0 ! Iterations with tolerance of 10-9 do while (del > 1.0d-9) X_A_old ) X_A do i ) 1,nComp do j ) 1,nSite sum1 ) 0.0d0 do k ) 1,nComp sum2 ) 0.0d0 do l ) 1,nSite if(l))j)cycle sum2 ) sum2 + S(l,k)*X_A_old(l,k)*delta(l,k,j,i) end do sum1 ) sum1 + x(k)*sum2 end do X_A(j,i) ) 1.0d0/(1.0d0 + rho*sum1) dif(j,i) ) dabs((X_A(j,i)-X_A_old(j,i))/X_A(j,i)) end do end do del ) maxval(dif) end do

Appendix B. Fortran Code for Solving Eq 6 with Respect to Density Using Matrix Inversion The values of Xji are assumed to be calculated in Appendix A above.

!Calculate the elements of matrix [Λpq] p ) 0; q ) 0 do j ) 1,nSite do i ) 1,nComp p)p+1 do l ) 1,nSite do k ) 1,nComp q)q+1 if(l))j) then if(k))i) then m_lambda(p,q) ) 1.0d0 else m_lambda(p,q) ) 0.0d0 end if cycle end if m_lambda(p,q) ) x(k)*rho*S(l,k)*delta(l,k,j,i)& &*(X_A(j,i))**2 end do end do q)0 end do end do !Calculate the inverse of matrix Λ using some inverting procedure !(m_lambda is the input and output of subroutine “inverse”) order ) nComp*nSite CALL inverse(order,m_lambda) !Calculate the elements of matrix [Ψp] p)0 do j ) 1,nSite do i ) 1,nComp sum1 ) 0.0d0 do k ) 1,nComp sum2 ) 0.0d0 do l ) 1,nSite if(l))j) cycle sum2 ) sum2 + S(l,k)*X_A(l,k)*d_delta_d_rho(l,k,j,i) end do sum1 ) sum1 + x(k)*sum2 end do p)p+1 psi(p) ) -(X_A(j,i))**2*(rho*sum1+(1.0d0/X_A(j,i)-1.0d0)/rho) end do end do !Pre-multiply Ψ by Λ-1 dX_A ) 0.0d0 do p ) 1,order do q ) 1,order dX_A(p) ) dX_A(p) + m_lambda(p,q)*psi(q) end do end do !Collecting the derivatives in proper subscripts p)0 do j ) 1,nSite do i ) 1,nComp p ) p+1 d_X_A_d_rho(j,i))dX_A(p) end do end do

Literature Cited (1) Wertheim, M. S. Fluids with Highly Directional Attractive Forces. I. Statistical Thermodynamics. J. Stat. Phys. 1984, 35, 19. (2) Wertheim, M. S. Fluids with Highly Directional Attractive Forces. II. Thermodynamic Perturbation Theory and Integral Equations. J. Stat. Phys. 1984, 35, 35. (3) Wertheim, M. S. Fluids with Highly Directional Attractive Forces. III. Multiple Attraction Sites. J. Stat. Phys. 1986, 42, 459. (4) Wertheim, M. S. Fluids with Highly Directional Attractive Forces. IV. Equilibrium Polymerization. J. Stat. Phys. 1986, 42, 477. (5) Chapman, W. G.; Gubbins, K. E.; Jackson, G.; Radosz, M. New Reference Equation of State for Associating Liquids. Ind. Eng. Chem. Res. 1990, 29, 1709-1721. (6) Mu¨ller, E. A.; Gubbins, K. E. Molecular-Based Equations of State for Associating Fluids: A Review of SAFT and Related Approaches. Ind. Eng. Chem. Res. 2001, 40, 2193-2211.

208

Ind. Eng. Chem. Res., Vol. 43, No. 1, 2004

(7) Kontogeorgis, G. M.; Voutsas, E.; Yakoumis, I.; Tassios, D. P. An Equation of State for Associating Fluids. Ind. Eng. Chem. Res. 1996, 35, 4310. (8) Kontogeorgis, G. M.; Yakoumis, I. V.; Meijer, H.; Hendriks, E.; Moorwood, T. Multicomponent phase equilibrium calculations for water-methanol-alkane mixtures. Fluid Phase Equilib. 1999, 158-160, 201. (9) Fu, Y.-H.; Sandler, S. I.; Orbey, H. A Modified UNIQUAC Model That Includes Hydrogen Bonding. Ind. Eng. Chem. Res. 1995, 34, 4351-4363. (10) Mengarelli, A. C.; Brignole, E. A.; Bottini, S. B. Activity coefficients of associating mixtures by group contribution. Fluid Phase Equilib. 1999, 163, 195-207. (11) Mu¨ller, E. A.; Gubbins, K. E. Associating Fluids and Fluid Mixtures. In Equations of State for Fluids and Fluid Mixtures, Part II; Sengers, J. V., Kayser, R. F., Peters, C. J., White, H. J., Jr., Eds.; Elsevier: Amsterdam, 2000; Vol. V; pp 435-477. (12) Kraska, T. Analytic and Fast Numerical Solutions and Approximations for Cross-Association Models within Statistical Association Fluid Theory. Ind. Eng. Chem. Res. 1998, 37, 48894892. (13) Huang, S. H.; Radosz, M. Equation of State for Small, Large, Polydisperse, and Associating Molecules. Ind. Eng. Chem. Res. 1990, 29, 2284-2294. (14) Elliott, J. R., Jr. Efficient Implementation of Wertheim’s Theory for Multicomponent Mixtures of Polysegmented Species. Ind. Eng. Chem. Res. 1996, 35, 1624-1629.

(15) Michelsen, M. L.; Hendriks, E. M. Physical properties from association models. Fluid Phase Equilib. 2001, 180, 165-174. (16) Michelsen, M. L. The isothermal flash problem. Part I. Stability. Fluid Phase Equilib. 1982, 9, 1-19. (17) Chapman, W. G.; Jackson, G.; Gubbins, K. E. Phase equilibria of associating fluids. Chain molecules with multiple bonding sites. Mol. Phys. 1988, 65, 1057. (18) Jackson, G.; Chapman, W. G.; Gubbins, K. E. Phase equilibria of associating fluids. Spherical molecules with multiple bonding sites. Mol. Phys. 1988, 65, 1. (19) Huang, S. H.; Radosz, M. Equation of State for Small, Large, Polydisperse, and Associating Molecules: Extension to Fluid Mixtures. Ind. Eng. Chem. Res. 1991, 30, 1994-2005. (20) Galindo, A.; Whitehead, P. J.; Jackson, G.; Burgess, A. N. Predicting the Phase Equilibria of Mixtures of Hydrogen Flouride with Water, Difluoromethane (HFC-32), and 1,1,1,2-Tetrafluoroethane (HFC-134a) Using a Simplified SAFT Approach. J. Phys. Chem. 1997, 101, 2082-2091. (21) Jiang, J.; Prausnitz, J. M. Molecular Thermodynamics for Protein Precipitation with a Polyelectrolyte. J. Phys. Chem. B 1999, 103, 5560-5569.

Received for review August 4, 2003 Revised manuscript received October 6, 2003 Accepted October 10, 2003 IE034041Q