6056
Ind. Eng. Chem. Res. 2006, 45, 6056-6062
GENERAL RESEARCH Efficient Solution of the Association Term Equations in the Statistical Associating Fluid Theory Equation of State Nikolaos M. P. Kakalis, Anthony I. Kakhu, and Constantinos C. Pantelides* Centre for Process Systems Engineering, Department of Chemical Engineering and Chemical Technology, Imperial College London, South Kensington Campus, London SW7 2AZ, United Kingdom
This paper proposes a hybrid successive substitution/modified Newton method for the solution of the nonlinear equations describing association effects in the SAFT equation of state. Extensive computational experiments are reported to show that the proposed method performs efficiently and reliably over wide ranges of conditions for multicomponent systems involving up to 57 distinct types of association sites. 1. Introduction The statistical associating fluid theory (SAFT) equation of state1,2 is particularly suited for the computation of physical properties of mixtures involving association between molecules, such as hydrogen bonding. The SAFT model assumes that each molecule may have zero or more association sites. Two sites on molecules of the same or different species will interact with a given energy once they are within a given range. Association interactions result in a contribution Aassoc to the mixture’s Helmholtz free energy of the form2
Aassoc
NC
)
RT
ni
xi ∑ ∑ i)1 a)1
(
ln Xai -
)
1 - Xai 2
(1)
where R is the ideal gas constant, T is the absolute temperature, NC is the number of components in the mixture, xi is the molar fraction of component i, and ni is the number of association sites on each molecule of component i. The quantity Xai denotes the fraction of sites a on molecules of component i that are not bonded. The Xai has been shown to satisfy the equations
1
Xai )
NC
1+F
nj
xj ∑ Xbj∆abij ∑ j)1 b)1
∀ i ) 1, ..., NC;
a ) 1, ..., ni (2)
where F is the mixture density and ∆abij is a quantity describing the association between site a on a molecule of species i and site b on a molecule of species j. In a typical evaluation of Helmholtz free energy, the fractions Xai are determined by solving the set of nonlinear algebraic equations 2; the values of both F and ∆abij are already known from earlier calculations. Because system 2 does not have a closed form analytical solution, it has to be solved numerically using an iterative method. Once the Xai are known, they can be used in eq 1 to determine Aassoc. The computation of the association contributions to most other physical properties of * To whom correspondence should be addressed. Tel: +44 20 7594 6622. Fax: +44 20 7564 6606. E-mail:
[email protected].
practical interest (e.g., fugacities) also requires the partial derivatives (“sensitivities”) of Xai with respect to the specified mixture composition, temperature, and volume. Overall, the overheads associated with the treatment of the association effects tend to make SAFT-based physical property evaluations much more expensive than those based on more widely used thermodynamic models. This is a particularly serious drawback if the EOS is used within a process model, in which case this numerical solution will need to be performed repeatedly. A typical dynamic process simulation may require hundreds or thousands of such solutions. Much of the work in the SAFT literature has employed successive substitution for the solution of eqs 2. As will be explained in more detail in Section 3.1, this numerical method offers robustness but lacks efficiency. This has motivated the consideration of special cases that allow the system 2 to be simplified to a form that can be solved more easily and efficiently.3-6 These approaches have also been extended to the computation of the sensitivities of Xai with respect to the external variables; however, the assumptions on which these simplifications are based restrict their applicability. An alternative is to employ an interval Newton/generalized bisection algorithm.7 This is an accurate and robust approach, but the additional computational cost is justified only if the SAFT equations are to be solved in the context of a more complex system that may have multiple solutions (e.g., for the determination of phase stability). Michelsen and Hendriks8 proposed an approach based on minimization of a state function (internal Helmholtz energy). This can calculate the various thermophysical property association contributions explicitly without any simplifying assumptions and without the need for the evaluation of the first partial derivatives of the solution of system 2. On the basis of the convexity properties of this function, the solution of the basic nonlinear system is achieved using an optimization-based approach that makes use of a trust region method that ensures quadratic convergence. This approach was also used by von Solms et al.,9 who presented a simplified formulation that is computationally faster than the original PC-SAFT EOS.10 Xu et al.7 and, later, Tan et al.11 showed that the first and second partial derivatives of Xai can be obtained at minimal cost via the solution of linear systems. The matrices involved
10.1021/ie051417m CCC: $33.50 © 2006 American Chemical Society Published on Web 07/21/2006
Ind. Eng. Chem. Res., Vol. 45, No. 17, 2006 6057
in these systems are computed analytically. Albeit robust and completely general, the methodology employs successive substitution for the solution of the nonlinear system. In a subsequent communication, Tan et al.12 presented computational time comparisons between their approach and the one by Michelsen and Hendriks8 for the calculation of first and second partial derivatives of the Helmholtz energy with respect to density. Their results indicated that the Michelsen and Hendriks method is faster as far as the computation of first derivatives is concerned; the two methods exhibit similar efficiency in the calculation of second derivatives. As is evident from the above brief review, much of the work to date has focused on the calculation of the nonlinear system sensitivities; however, the real computational bottleneck is the solution of the system itself. In this paper, we show that a hybrid approach, combining successive substitution with a modified Newton method, provides an efficient and robust way for the solution of the embedded nonlinear systems. Section 2 considers some aspects of the formulation of the association term of SAFT. Section 3 reviews the use of successive substitution and Newtontype methods for the solution of the association equations and describes the proposed hybrid solution algorithm for the efficient solution of the nonlinear system and the computation of its sensitivities. Finally, Section 4 presents computational results to demonstrate the efficiency and robustness of the suggested approach.
If we consider two sites a, a′ of the same type aj and two sites b, b′ of type bh, then we have ∆abij ) ∆a′b′ij, and we can therefore denote all these identical parameters as ∆ajbhij. Therefore, NSTj
∑
bh)1
Xbhj
∑
NSTj
∆abij )
b∈Sbhj
∑
bh)1
1
Xaji )
NC
1+F
(3)
NSTj
xj ∑ ∑ Xbhj∆abij ∑ j)1 bh)1 b∈S bhj
The inner double summation in the denominator can be written as NSTj
NSTj
Xbhj ∑ ∆abij ∑ ∑ Xbhj∆abij ) bh∑ )1 b∈S
bh)1 b∈Sbhj
bhj
(4)
b∈Sbhj
1
Xaji )
∑ ηbhjXbhj∆ajbhij
(5)
bh)1
NC
1+F
, ∀i ) 1, ..., NC;
NSTj
xj ∑ ηbhjXbhj∆ajbhij ∑ j)1 bh)1
aj ) 1, ..., NSTi (6)
NC which involves only ∑i)1 NSTi equations (and variables Xaji) NC instead of the ∑i)1 ni equations (and variables Xai) in 2. We can now also rewrite 1 as
Aassoc
NC
)
RT
[ ( NSTi
xi ∑ ηaji ln Xaji ∑ i)1 aj)1
1 - Xaji 2
)]
(7)
2.2. Standard Form of Nonlinear System. The nonlinear system 6 can be rewritten in the standard form F(X) ) 0, where NSTj
NC
The system of nonlinear eqs 2 needs to be solved each time the SAFT equation is used for the evaluation of properties of any process mixture involving association effects. We note that 2 treats each site a of each component i as a distinct entity; hence, the total number of variables Xai (and equations) in this NC ni. However, in most cases of practical interest, system is ∑i)1 the sites on a given molecule can be classified into a small number (usually fewer than three) of distinct types. All sites belonging to the same type are indistinguishable from each other as far as the SAFT EOS is concerned. This allows us to write formulation 2 in a more compact form. This observation was first made by Kraska6 for certain specific types of binary mixtures, but was later presented in a more general form by Gross et al.13 The formulation presented below is essentially identical to the one in the latter paper. 2.1. Reformulation of Association Equations. Let the number of distinct site types in component i be denoted by NSTi and the number of sites of component i that belong to a particular type aj ) 1, ..., NSTi be denoted by ηaji. Clearly, i ∑NST aj)1 ηaji ) ni. Let Saji denote the set of sites of component i that are of type aj. Since two sites a, a′ of the same type aj are indistinguishable, we would expect Xai ) Xa′i. We can therefore replace all variables Xai, a ∈ Saji by a single variable Xaji, and eq 2 can be rewritten as
∑
NSTj
∆ajbhij )
Overall, then, 2 can be rewritten in the more compact form
Faji ≡ Xaji + FXaji
2. The Association Term Equations
Xbhj
∑ ∑ (ηbhjXbhj∆ajbhij) - 1 ) 0; j)1 bh)1 xj
∀i ) 1, ..., NC; aj ) 1, ..., NSTi (8)
The Jacobian matrix of the above system is given by
∂Faji ∂Xaji
NC
)1+F
NSTj
xj ∑ (ηbhjXbhj∆ajbhij) + FXajixiηaji∆ajajii ∑ j)1 bh)1
∂Faji ) FXajixjηbhj∆ajbhij; ∀(bh, j) * (aj, i) ∂Xbhj
(9)
(10)
∀i ) 1, ..., NC; aj ) 1, ..., NSTi. We note that the first two terms of the right-hand side of eq 9 are simply the denominator of the fraction of the right-hand side of 6. Consequently, if the values of Xaji satisfy eq 6, then we can rewrite 9 in a shorter and more computationally efficient way as
|
∂Faji 1 + FX/ajixiηaji∆ajajii; ∀i ) 1, ..., NC; ) ∂Xaji / X/ aji aj ) 1, ..., NSTi (11) where the asterisk (*) indicates values computed at the solution of eq 6. The above Jacobian matrix is diagonally dominant and nonsingular and, hence, is always invertible. Therefore, as Xaji ∈ (0, 1) ∀i, aj, one can conclude that the system has a unique solution in a unit hypercube.7 3. Numerical Solution Methods for the Association Equations Most of the work reported in the literature on the practical use of the SAFT EOS is based on the use of successive substitution for the solution of the nonlinear system 8. This, unfortunately, leads to slow convergence or to solutions that often do not satisfy the equations very tightly. Although
6058
Ind. Eng. Chem. Res., Vol. 45, No. 17, 2006
convergence to tight tolerances may not be essential for standalone evaluations of physical properties, it is crucial for the use of EOS within the wider context of a process simulation or optimization calculation. This section briefly reviews the application of successive substitution and Newton-type methods to the solution of eq 8 and compares the performance of these methods via a simple example. It then proceeds to propose a hybrid approach that combines the two methods. In the interest of simplicity of notation, we will henceforth replace the symbols aj, bh used to denote site types with the simpler symbols a, b. 3.1. Successive Substitution. Successive substitution considers a nonlinear system of the matrix form X ) G(X). Beginning with an initial estimate (“guess”) of the solution vector, X(0), and a convergence tolerance, , the method computes a sequence of iterates X(k) via the formula
X(k+1) ) G(X(k)); k ) 0, 1, ...
(12)
until ||X(k+1) - X(k)|| < for some suitable norm ||‚||. An obvious form (cf. eq 6) of the function G(X) for the association system is
1
G(Xai) )
NSTj
NC
1+F
(13)
∑ ∑ ηbjXbj∆abij j)1 b)1 xj
A necessary condition for the convergence of the method to a solution X* (i.e., limkf∞ X(k) ) X*) is
F
(X*)] < 1 [∂G ∂X
(14)
where F(M) is the spectral radius of a matrix M (i.e., the magnitude of its largest eigenvalue). Whether this condition is satisfied depends on the (generally nonunique) choice of the function G(X). If this condition is not satisfied, then successive substitution will not converge, even if one starts arbitrarily close to the solution X*. Another potential problem with successive substitution is the speed of convergence. Even when it does converge, the method exhibits only first-order convergence, that is, ||X(k+1) - X*|| e R||X(k) - X*||, with 0 < R < 1. That means that tight convergence may require many iterations. On the other hand, successive substitution does not require matrix manipulations and, hence, has a relatively small cost per iteration. 3.2. Newton-Type Methods. In the standard Newton-type method, starting from an initial guess, X(0), the solution vector is iteratively corrected by a step δX computed as the solution of the following set of linear equations.
|
∂F δX(k) ) -F(X(k)); k ) 0, 1, ... ∂X X(k)
(15)
The above linear equations can be solved using, for instance, LU decomposition, and the corrections obtained are applied to the solution vector
X(k+1) ) X(k) + δX(k)
(16)
The iterations are repeated until a convergence criterion, such as ||X(k) - X(k-1)|| e or ||F(Xk)|| e , is satisfied for a given accuracy . Newton’s method can be shown to be locally convergent,14 that is, if ||X(0) - X*|| e for a sufficiently small number ,
Table 1. Computational Comparison: Temperature Variation conditions
successive substitution
Newton
T (K)
ζ
iterations
spectral radius
iterations
180 200 240 270 300 330 370 440 500
0.43 0.43 0.43 0.43 0.43 0.43 0.43 0.43 0.43
2315 1471 756 526 396 315 247 179 147
0.9875 0.9804 0.9621 0.9461 0.9291 0.9116 0.8888 0.8510 0.8216
13 12 11 10 10 9 9 8 8
Table 2. Computational Comparison: Density Variation conditions
successive substitution
Newton
T (K)
ζ
iterations
spectral radius
iterations
330 330 330 330 330 330 330
0.430 0.450 0.475 0.531 0.600 0.694 0.730
315 344 374 457 594 876 1033
0.9116 0.9188 0.9251 0.9382 0.9521 0.9673 0.9722
9 9 10 10 10 11 11
then, provided that the Jacobian ∂F(X*)/∂X is nonsingular, limkf∞X(k) ) X*. Thus, if we start sufficiently close to the solution, it will always converge. It is also locally quadratically convergent, that is, ||X(k+1) - X*|| e R||X(k) - X*||2, with 0 e R e 1. Thus, it converges very quickly once it gets in the neighborhood of the solution. In our case, the Jacobian matrix ∂F/∂X|X(k) is given by eqs 9 and 10. Modified Newton methods compute the step δX(k) by solving equations of the form
A(k)δX(k) ) -F(X(k)); k ) 0, 1, ...
(17)
where the matrix A(k) is an approximation to the true Jacobian ∂F/∂X|X(k). It can be shown (see, for instance, theorem 3.1 in Dennis and More15) that modified Newton methods for which limkf∞||A(k) - ∂F/∂X|X(k)|| ) 0 retain the superlinear convergence characteristics of the original Newton method. This is, indeed, the case if we replace eq 9 by the computationally more efficient 11 as the two expressions yield identical results as X(k) f X*. 3.3. Numerical Comparison. To better understand the differences in the efficiency of the solution approaches described in Sections 3.1 and 3.2 in the context of SAFT, we apply them to the solution of the association equations for a mixture of hydrogen fluoride with water. The SAFT models of the two molecules are illustrated in Figure 1a and b. Water molecules have four association sites of two distinct types. Hydrogen fluoride molecules have two sites, each of a different type. A complete set of SAFT model parameters for this system is given by Grice.16 We consider a mixture containing 68.28% HF (by mole) and perform a number of computational experiments by varying the temperature and packing fraction (ζ). The latter is closely related to the density of the mixture. For each experiment, we report the number of iterations required for the convergence of the successive substitution and modified Newton methods, as well as the spectral radius of the Jacobian matrix at the solution. All solutions are carried out to a convergence tolerance of ) 10-12. In Table 1, the packing fraction is kept constant and the temperature is varied, whereas in Table 2, temperature is constant and ζ is varied. Some of these tests correspond to extreme conditions, under which the system may cease to be in the fluid phase. These are deliberately chosen to examine the numerical behavior of the algorithm. We observe that at low
Ind. Eng. Chem. Res., Vol. 45, No. 17, 2006 6059 Table 3. Computational Comparison: Worst Case conditions
successive substitution
Newton
T (K)
ζ
iterations
spectral radius
iterations
180
0.73
8444
0.9966
15
temperatures or high densities, the spectral radius grows and approaches unity. As expected, this causes successive substitution to require very large numbers of iterations. On the other hand, the Newton method is much more efficient, requiring just one or two additional iterations. Table 3 shows results for the worst-case combination of low temperature and high density. In this case, successive substitution needs 8444 iterations, in contrast to Newton, which still requires only 15. Overall, these results indicate that successive substitution may not be appropriate for the solution of the SAFT association equations in the context of process modeling applications. A method that is often used for the acceleration of successive substitution iterations is successive overrelaxation (SOR) (see ref 14, section 7.4). In the current context, this would involve replacing the iterative scheme 12 by
X(k+1) ) ωG(X(k)) + (1 - ω)X(k); k ) 0, 1, ...
Figure 2. Effect of relaxation parameter on computational performance of SOR.
(18)
where ω ∈(0, 1] is a relaxation parameter. The rate of convergence of this scheme is determined by the spectral radius
[
∂G (X*) + (1 - ω)I ∂X
Fω
]
(19)
For ω ) 1, the SOR scheme, 18, reduces to the original successive substitution one, 12; however, smaller values of ω may result in smaller spectral radii and, therefore, faster convergence. In principle, this is an attractive possibility because the computational cost of 18 is only marginally higher than that of 12. To investigate the potential of SOR, we repeat three sets of computations (the first and last cases shown in Table 1 and also the case shown in Table 3) using SOR with a range of relaxation parameter values, ω. Figure 2 shows the acceleration factor achieved by the use of SOR expressed in terms of the ratio of the number of iterations required by simple successive substitution (as already given in Tables 1 and 3) divided by the number of iterations required by SOR. Acceleration factors of up to 75 can be achieved for ω ≈ 0.8. Nevertheless, two of the three cases considered still require more than 100 iterations to converge, which is much higher than the corresponding numbers of Newton iterations (shown in the last columns of Tables 1 and 3). Of course, since the cost per SOR iteration is smaller than the cost of a Newton iteration, it is possible that SOR will outperform Newton in some cases. However, this is unlikely to be the case in situations in which the initial guess is close to Figure 3. Hybrid algorithm for the solution of the association term equations and computation of the solution sensitivities.
Figure 1. H2O and HF modeling according to the SAFT formalism.
the solution (as usually happens, for example, in a dynamic process simulation context). 3.4. Hybrid Solution Method. The computational comparison presented in Section 3.3 illustrates the efficiency advantages of Newton-type methods over successive substitution. The main disadvantage of the former in the context of the application of interest to this paper is that superlinear convergence is obtained only in the neighborhood of the solution. On the other hand, successive substitution has the inherent characteristic of driving
6060
Ind. Eng. Chem. Res., Vol. 45, No. 17, 2006
the solution trajectory closer to the solution, and each iteration is relatively inexpensive when compared to a Newton-type iteration. The above considerations lead us to consider a hybrid solution approach that starts with a small number of successive substitution iterations aimed at bringing the solution point within the Newton region of convergence. The method then switches to a modified Newton iteration that ensures fast convergence to the solution within tight tolerances. A description of this hybrid algorithm is shown in Figure 3. The linear systems, 15, are solved using Choleski factorization combined with backward substitution. As already mentioned, the Jacobian of the system is always nonsingular, and therefore, its LU form always exists. On a more practical level, the use of analytical derivatives 10 and 11 improves both the reliability and the efficiency of the solution. 3.5. Evaluation of System Sensitivities. The evaluation of most thermophysical properties of practical interest requires the partial derivatives ∂Xai/∂Ψ, where Ψ ) {T, V, N}. The standard approach to deriving such sensitivities is by applying the implicit function theorem to the nonlinear system 8, written in the form F(X,Ψ) ) 0. By differentiating this with respect to Ψ at a solution point X*, we obtain
|( )
|
∂F ∂X* ∂F + )0 ∂X X* ∂Ψ ∂Ψ X*
(20)
from which we can obtain the required sensitivities of the solution X* with respect to Ψ by the solution of the linear system
|( )
|
∂F ∂F ∂X* )∂X X* ∂Ψ ∂Ψ X*
(21)
as implemented at step 9 of the algorithm of Figure 3. The matrix ∂F/∂X|X* is the Jacobian of the system 8 evaluated at the solution. The partial derivatives ∂F/∂Ψ can be evaluated by partial differentiation of 8
∂Fai ∂T ∂Fai
)
∂Fai
( )
1 1
V X/ ai
∂V
NSTj
NC
) FX/ai
∑ ∑ j)1 b)1
-1 +
xj
NC
FX/ai
(
ηbjX/bj NSTj
xj ∑ ∑ j)1 b)1
)
∂∆abij ∂T
(
ηbjX/bj
(22)
)
∂∆abij ∂V
(23)
)
∂Nk FX/ai
[
NSTk
∑ b)1
NC
(ηbkX/bk∆abik)
+
NSTj
xj ∑ ∑ j)1 b)1
(
ηbjX/bj
)]
|( ) |( X*
∂X* ∂F + ∂Ψ ∂X
X*
)
∂2X* ∂2F + ∂Ψ ∂Ψ′ ∂Ψ ∂Ψ′
(24)
∀k ) 1, ..., NC, i ) 1, ..., NC, and a )1, ..., NSTi. We note that the above sensitivity expressions are practically identical to those given by Xu et al.7 On the other hand, those derived by Tan et al.11 can be obtained immediately via the application of the implicit function theorem to the original eq 2 rather than the rearranged form 8. The extension of this approach to the calculation of the second-order sensitivities is straightforward. Subsequent differentiation of relation 20 at the solution point X* with respect to Ψ′ ) {T, V, N} yields
|
X*
)0
(25)
from which the second-order sensitivities ∂2Xai/∂Ψ ∂Ψ′ of the solution X* are derived as the solution of the linear system
∂F ∂X
|( X*
) [
∂2X* ∂2F )∂Ψ ∂Ψ′ ∂X ∂Ψ′
|( )
∂X* ∂2F + X* ∂Ψ ∂Ψ ∂Ψ′
|] X*
(26)
where the left-hand matrix is again the already available Jacobian of the nonlinear system and the partial derivatives ∂2F/ ∂X ∂Ψ′ and ∂2F/∂Ψ ∂Ψ′ are derived analytically in the same fashion as above. Note that at the last step of the procedure, 9b, the lower and upper triangular matrices L(l) and U(l) are the same as those used at the previous execution of step 6. Provided the convergence tolerance is sufficiently tight, we do not need to recalculate and refactorize the Jacobian at step 9. 4. Computational Experiments Any implementation of an equation of state to be used within a process modeling environment must be capable of performing not only efficiently but also reliably over wide ranges of conditions; even a single failure among several hundreds of successful evaluations may lead to failure of the entire process simulation. With this in mind, this section assesses the performance of the algorithm of Section 3.4 via extensive computational testing. 4.1. Testing Methodology. We start by writing the nonlinear system 8 in the form NC NSTj
F(Xai) ) Xai + Xai
∑ ∑ (Xbj∆˜ abij) - 1 ) 0 j)1 b)1
(27)
where we have combined all known quantities within the new parameters ∆ ˜ abij ≡ Fxjηbj∆abij. We now proceed to apply the algorithm to 27 over a wide range of values of ∆ ˜ . We consider values for the molar fractions xj in the range [10-5, 1.0]. In addition, examination of different molecules and mixtures under both normal and extreme conditions indicates that F∆ ∈ [10-5, 104]. Hence, ∆ ˜ ∈ [10-10, 104]. Our study considers mixtures involving from 1 to 19 components, each of which has three distinct types of association sites (NSTi ) 3), with a single site of each type (ηia ) 1, a ) 1, 2, 3). Thus, we solve systems with up to 57 distinct site types. All pairwise interactions are assumed to be possible. The ∆ ˜ abij values are created randomly in the range (10-10, 4 10 ) from a uniform [0, 1] distribution of random numbers r via the expression
∆ ˜ abij ) 1014r-10
∂∆abij ∂Nk
∂2F ∂X ∂Ψ′
(28)
Each test involves the solution of 10 000 systems, 27, with ∆ ˜ abij generated in the above manner, on an Intel Pentium 4 processor operating at 2.4 GHz with 256 MB of RAM. In the interest of brevity, the rest of this section reports results obtained for mixtures with 1, 6, 9, 12, 16, and 19 components; however, computations were also carried out with mixtures with all other numbers of components between 1 and 19. 4.2. Numerical Comparisons. Comparative computational results for successive substitution and modified Newton are given in Table 4, which reports average, minimum, and maximum numbers of iterations, as well as the standard deviations of the number of iterations over 10 000 computational
Ind. Eng. Chem. Res., Vol. 45, No. 17, 2006 6061 Table 4. Computational Comparison of Successive Substitution and Modified Newton Methods successive substitution iterations
modified Newton iterations
distinct site types
av
min
max
σ
CPUSS (ms)
av
min
max
σ
CPUNWT (ms)
CPUSS ÷ CPUNWT
3 18 27 36 48 57
316.3 505.6 748.8 1103.9 1656.7 2067.9
4 30 53 97 146 309
2858 2706 2635 2737 3258 3444
567.9 407.6 426.3 465.9 505.5 496.7
0.354 13.010 42.408 108.963 295.384 507.982
7.2 11.6 11.8 12.0 12.0 12.0
3 9 10 11 11 12
12 26 24 16 14 13
2.0 0.9 0.7 0.4 0.2 0.1
0.107 2.240 5.223 8.650 16.139 22.460
3.3 5.8 8.1 12.6 18.3 22.6
experiments. The average CPU time per system solution is also reported. As can be seen, the use of Newton-type iterations results in a very significant reduction in the average CPU time. This benefit increases as the size of the system increases, but even in the case of small-scale systems, a 3-fold increase in speed is achieved. We now proceed to consider the performance of the hybrid method described in Section 3.4. This depends on the number NSSITER of successive substitution iterations taken at step 2 of the algorithm. Table 5 shows computational statistics for NSSITER 4, 8, 10, and 15; additional runs with NSSITER )
Table 5. Computational Performance of Hybrid Method NSSITER ) 4 distinct site types 3 18 27 36 48 57
4.3 6.2 5.8 5.5 5.2 5.0
1 3 4 4 4 4
65 26 40 16 11 9
CPUHYB (ms)
2.6 1.7 1.3 0.9 0.6 0.6
0.068 1.291 2.793 4.367 7.636 10.287
Newton iterations av min max σ 3.7 5.1 4.7 4.4 4.1 3.8
1 2 3 3 3 3
NSSITER ) 10 distinct site types 3 18 27 36 48 57
Figure 4. Effect of number of initial successive substitution iterations on computational performance for a system with 57 distinct types of association sites.
Newton iterations av min max σ
NSSITER ) 8
Newton iterations av min max σ 3.6 4.8 4.4 4.0 3.6 3.4
1 1 2 3 3 2
13 22 15 11 8 7
2.6 1.6 1.2 1.0 0.8 0.7
23 23 29 11 9 9
CPUHYB (ms)
2.6 1.6 1.2 0.9 0.7 0.7
0.064 1.196 2.541 3.971 6.883 9.096
NSSITER ) 15
CPUHYB (ms) 0.065 1.186 2.509 3.909 6.676 8.871
Newton iterations av min max σ 3.2 4.2 3.8 3.5 3.2 2.9
1 1 2 2 2 2
9 15 12 8 8 8
2.5 1.5 1.2 0.9 0.8 0.7
CPUHYB (ms) 0.065 1.199 2.539 4.010 6.911 9.115
20 and NSSITER ) 25 were also performed but are not reported here. Additionally, Figure 4 shows the CPU time (averaged over 10 000 experiments) for a system with 57 distinct site types. These results clearly suggest that there is a nonmonotonic effect of NSSITER on overall computational performance. A few successive substitution iterations applied initially have a very beneficial effect because they bring the variables to the region of quadratic convergence of the modified Newton method with little cost. On the other hand, using too many of these iterations introduces an unnecessary cost because Newton-type iterations are much more effective in achieving tight convergence. The optimal number of initial successive substitution iterations seems to be around 10. Figure 5 shows the dramatic reduction in the number of Newton-type iterations achieved by introducing these relatively cheap iterations, each one of which essentially involves a single evaluation of the function 13. Finally, the proposed approach is robust as well as efficient: the hybrid method was tested with 190 000 computational experiments for each value of NSSITER, and every single one of these successfully converged to the solution. 5. Conclusions
Figure 5. Comparison of numbers of Newton iterations required by the hybrid method (NSSITER ) 10) and the modified Newton method for a system with 36 distinct site types.
This paper has proposed a hybrid algorithm for the solution of the nonlinear association term equations in the SAFT EOS. The method comprises an initial phase involving 10 successive substitution iterations, followed by a second phase involving a modified Newton method. The latter employs a Jacobian matrix approximation that is simpler to evaluate than the true Jacobian matrix but converges to it at the solution. The method is able to handle mixtures of any number of components with any number of associating sites per molecule. Extensive computational tests have been carried out to demonstrate that it performs both reliably and efficiently over wide
6062
Ind. Eng. Chem. Res., Vol. 45, No. 17, 2006
ranges of conditions. The exact partial derivatives of the system solution with respect to temperature, volume, and composition, needed for the calculation of physical properties, can also be computed via the solution of standard sensitivity equations at minimal additional cost. We note that the successive substitution part of the hybrid algorithm could be replaced with an SOR iterative scheme, which could lead to slightly fewer iterations being required during the first phase of the algorithm. A detailed investigation is required to determine the optimal number of such iterations. Acknowledgment The authors thank Professor George Jackson and Dr. Amparo Galindo for making available to us the SAFT parameters for the HF/water system and for numerous helpful comments and suggestions. We also thank an anonymous referee for suggesting to us to investigate the possible use of SOR. Our work was funded by the United Kingdom’s Engineering and Physical Sciences Research Council (EPSRC) under Platform Grant GR/ N08636 and the Alexander S. Onassis Public Benefit Foundation. Literature Cited (1) Chapman, W. G.; Gubbins, K. E.; Jackson, G.; Radosz, M. SAFT: Equation-of-state solution model for associating fluids. Fluid Phase Equilib. 1989, 52, 31. (2) 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. (3) Huang, S. H.; Radosz, M. Equation of state for small, large, polydisperse, and associating molecules. Ind. Eng. Chem. Res. 1990, 29, 2284. (4) 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.
(5) Elliott, J. R. Efficient implementation of Wertheim’s theory for multicomponent mixtures of polysegmented species. Ind. Eng. Chem. Res. 1996, 35, 1624. (6) Kraska, T. Analytic and fast numerical solutions and approximations for cross-association models within statistical association fluid theory. Ind. Eng. Chem. Res. 1998, 37, 4889. (7) Xu, G.; Brennecke, J. F.; Stadtherr, M. A. Reliable computation of phase stability and equilibrium from the SAFT equation of state. Ind. Eng. Chem. Res. 2002, 41, 938. (8) Michelsen, M. L.; Hendriks, E. M. Physical properties from association models. Fluid Phase Equilib. 2001, 180, 165. (9) von Solms, N.; Michelsen, M. L.; Kontogeorgis, G. M. Computational and physical performance of a modified PC-SAFT equation of state for highly asymmetric and associating mixtures. Ind. Eng. Chem. Res. 2003, 42, 1098. (10) Gross, J.; Sadowski, G. Application of the perturbed-chain SAFT equation of state to associating systems. Ind. Eng. Chem. Res. 2002, 41, 5510. (11) Tan, S. P.; Adidharma, H.; Radosz, M. Generalized procedure for estimating the fractions of nonbonded associating molecules and their derivatives in thermodynamic perturbation theory. Ind. Eng. Chem. Res. 2004, 43, 203. (12) Tan, S. P.; Adidharma, H.; Radosz, M. Reply to Comments on “Generalized procedure for estimating the fractions of nonbonded associating molecules and their derivatives in thermodynamic perturbation theory”. Ind. Eng. Chem. Res. 2004, 43, 6263. (13) Gross, J.; Spuhl, O.; Tumakaka F.; Sadowski, G. Modeling copolymer systems using the perturbed-chain SAFT equation of state. Ind. Eng. Chem. Res. 2003, 42, 1266. (14) Ortega, J. M.; Rheinboldt, W. C. IteratiVe Solution of Nonlinear Equations in SeVeral Variables; Academic Press Inc: New York, 1970. (15) Dennis, J. E.; More, J. J. Quasi-Newton methods, motivation and theory. SIAM ReV. 1977, 19, 46. (16) Grice, S. J. Predicting the phase equilibria of associating and reacting systems. Ph.D. Thesis, University of Sheffield, United Kingdom, 1999.
ReceiVed for reView December 20, 2005 ReVised manuscript receiVed May 17, 2006 Accepted May 22, 2006 IE051417M
