Ind. Eng. Chem. Res. 2006, 45, 8449-8453
8449
Robust and Efficient Solution Procedures for Association Models Michael L. Michelsen* IVC-SEP, Institut for Kemiteknik, DTU, DK 2800 Lyngby, Denmark
Equations of state that incorporate the Wertheim association expression are more difficult to apply than conventional pressure explicit equations, because the association term is implicit and requires solution for an internal set of composition variables. In this work, we analyze the convergence behavior of different solution methods and demonstrate how a simple and efficient, yet globally convergent, procedure for the solution of the equation of state can be formulated. Introduction Association models such as statistical associating fluid theory (SAFT)1,2 and cubic plus association (CPA)3 have become very popular during the last 15 years, in particular in academia. Surprisingly, they do not yet seem to appear in their full form, i.e., including the association term, in commercial simulators. Representatives from one simulator manufacturer4 claim that this is due to a resistance from industry, because of the much higher level of complexity with the equation of state. Even in its nonassociating variants, the SAFT family of equations of state cannot be solved analytically for density. The addition of the association term adds a new set of internal variables, the fractions of association sites that are not occupied, and which are functions of the mixture density. A robust approach for finding the density would typically involve nested loops, where the fractions of free sites are calculated at a fixed density in an inner loop, with density being modified in the outer loop to match a specified pressure. Clearly, the complexity of this approach is much higher that that of, e.g., solving a cubic equation of state analytically, and it becomes important that the density iteration and, in particular, the inner multivariable iteration for determining the unbonded sites fraction is rapid, as well as extremely robust. The present paper examines computational approaches that can match these requirements.
components but not on the fractions of unbonded sites. The total number of components is C. At a maximum, the gradient of Q must equal zero, yielding
∂Q ∂XAi
1
XAi
C
nj
∆ A B XB ) 0 ∑j V∑ B
- 1 - ni
j
1
X Ai )
C
1+
∆A B XB ∑j (nj/V)∑ B i j
i
j
j
The notation, with one index for the molecule and another for the different sites on the molecule, is quite cumbersome, and we shall therefore simplify it as follows. Assume, for example, that we consider a three-component mixture where component 1 has two sites, component 2 has three sites, and component 3 has no sites. The sites on component 1 are indexed X1 and X2 and those on component 2 are indexed as X3, X4, and X5. The number of moles that host a given site k are denoted mk, and the total number of different sites is S. In our example, therefore, m1 ) m2 ) n1, m3 ) m4 ) m5 ) n2, and S ) 5. Using this notation, we may write
Q(X,m) )
1
S
S
∑k mk(ln Xk - Xk + 1) - 2V∑k ∑l mkmlXkXl∆kl
(4)
with
∂Q ∂Xk
(ln XA - XA + 1) ∑i ni∑ A
(2)
(3)
C
S
The computationally demanding calculation in the association models is to determine the fraction XAi of sites Ai that are not bonded to other sites. Here, Ai denotes the number of sites of type A on molecule i. Michelsen and Hendriks5 showed that the site fractions could be determined as the maximum of the function
i j
j
which is frequently presented in the form
Base Equations
Q(X, n) )
( )
) ni
) mk
( ) 1
Xk
S
- 1 - mk
∑l
ml ∆klXl ) 0 V
(5)
i
which leads to
i
1 2V
C
C
∑i ∑j ninj∑ ∑XA XB ∆A B A B
i j
i
i
j
(1)
j
Here, ni is the number of moles of component i, V the total volume, and ∆AiBj the association strength between A-sites on molecule i and B-sites on molecule j. The association strength is dependent on the temperature and molar density of the * To whom correspondence should be addressed. Tel.: (+45)45252865. Fax: (+45)45882258. E-mail:
[email protected].
1
Xk )
S
1+
∑l (ml/V)∆klXl
To simplify the notation further, we introduce
Klk ) Kkl )
10.1021/ie060029x CCC: $33.50 © 2006 American Chemical Society Published on Web 04/22/2006
mlmk∆lk V
(6)
8450
Ind. Eng. Chem. Res., Vol. 45, No. 25, 2006
The first step in the solution for X is to determine the elements of K. The objective function now becomes S
Q(X,m) )
∑k
mk(ln Xk - Xk + 1) -
1 2
S
S
∑k ∑l KklXkXl
(7)
with gradient vector
gk )
∂Q ∂Xk
) mk
( )
and Hessian matrix
Hkl ) -
1
Xk
S
-1 -
∑l KklXl
Xk2
(8)
δkl - Kkl
(9)
Note that, apart from the diagonal, the Hessian matrix is constant. Solution Procedures Successive Substitution. The form of eq 3 immediately suggests the use of a simple solution procedure using successiVe substitution. In terms of the variables used here, we obtain
mk
) fk(X(n)) ) X(n+1) k
(10)
S
mk +
∑l KklX(n)l
Successive substitution will be rapidly convergent when f is only weakly dependent on X, i.e., when the association constants Kkl are small. Convergence will be linear, with a rate determined by the numerically largest eigenvalue of the matrix J, with elements given by
Jkl )
∂fk
(mk +
∑l KklXk)2
1 f X1 ) 1 + k1X1 1+
-k1 (1 + k1X1)2
Clearly, ω ) 1/2 is the optimal choice when k1 is very large, with the price being that the damping will reduce the rate of convergence when K is small. A reasonable compromise is a value in the range of 0.2-0.3. Association Scheme 2B. Here, the associating molecule is assumed to have two sitessa “positive” site and a “negative” sitesand both can only associate with a site of the opposite polarity. The scheme is frequently used for alcohols. Thus, we have S ) 2, m1 ) m2 ) 1, K11 ) K22 ) 0, and K12 ) K21 ) k1, and the fraction of sites not bonded is determined from
X1 )
1 1 , X2 ) 1 + k 1X 2 1 + k1X1
Clearly, X1 ) X2 and we arrive at the same result as in the 1A case:
J)
2
x4k1 + 1
At the solution, the matrix J becomes
J11 ) λ1 )
) -k1X12 ) X1 - 1
(14)
2 1 + x4k1 + 1
The Jacobian, which is evaluated at the solution, now becomes
The general behavior of the method is best analyzed through the simple example of a single, pure component that can selfassociate. We consider, using the terminology of Huang and Radosz,2 three different association schemes, 1A, 2B, and 3B, where, for each, we analyze the convergence behavior in dependence of the magnitude of the association constant. Association Scheme 1A. In this scheme, the associating molecule is assumed to have a single site that can associate with a similar site on a different molecule. Thus, S ) 1 and m1 ) n1 ) 1 and K11 ) k1 ) ∆11/V. The scheme is typically used for organic acids. The fraction of sites not bonded is given by
X1 )
∑l
(13)
KklX(n) l
λd ) (1 - ω)λu + ω
(11)
S
+ ωXnk
S
where a value of ω ) 0 corresponds to the undamped case. The eigenvalue of the damped method is
X1 ) X 2 ) X )
mkKkl
)-
∂Xl
mk
X(n+1) ) (1 - ω) k mk +
() mk
convergence occurs. This is usually the case for a vapor phase where the molar volume V is large. For a liquid, however, k1 can be quite large (e.g., for acetic acid on the order of 105106) and, consequently, X1 becomes very small and the value of λ becomes close to -1, resulting in very slow convergence, where thousands of iterations are required. A simple remedy in this case is to use damped successive solution. We replace eq 10 by
(12)
The eigenvalue is always negative and convergence, therefore, oscillatory. When k1 is small, X1 will be close to 1 and rapid
(
0 X-1 X-1 0
)
(15)
with two eigenvalues of the same magnitude but opposite sign:
λ1 ) X - 1, λ2 ) 1 - X At first sight, this appears to destroy the speedup that damping provided. Proceeding as described for Scheme 1A with ω ) 0.25 yields an effective λ1 value that is closer to zero but unfortunately a value of the positive λ2 that is closer to 1 than in the undamped case also is attained. Therefore, the overall effect of damping seems to be a reduction in the rate of convergence. Fortunately, the increase in λ2 is inconsequential when X1 and X2 are identical, because the eigenvector associated with this eigenvalue is not “active”. However, there may be situations for multicomponent mixtures where sites on other molecules might be capable of interacting with one of the two sites but not the other and thereby create a situation where X1 and X2 differ, and, in such cases, damping might have a negative effect. We may also note that the 4C scheme, which is frequently used for water and for glycols, behaves similar to the 2B scheme, with the only difference being that two sites of each polarity is assumed.
Ind. Eng. Chem. Res., Vol. 45, No. 25, 2006 8451
Association Scheme 3B. This scheme has also been suggested for alcohols. Two identical sites are of one polarity and the third site, which can associate with the other two sites, has the opposite polarity. Therefore, we have S ) 3, m1 ) m2 ) m3 ) 1, and an association matrix,
The modification is performed as follows: We may write
0 0 k1 K ) 0 0 k1 k1 k1 0
where we have used eq 8. In the approximate Hessian, we drop the contribution gk (which is zero at the solution) and arrive at
( )
mk
1 mk
)
Xk2
Xk Xk
(
H ˆ kl ) -
resulting in
X 1 ) X2 )
1 1 , X3 ) 1 + k1X3 1 + 2k1X1
(16)
2
k1 + 1 + xk12 + 6k1 + 1 3k1 + 1 + xk12 + 6k1 + 1
The matrix J at the solution becomes
(
-k1X12 0 0 -k1X22 0 0 2 2 -k1X3 -k1X3 0
)
(17)
)
∑j KkjXj) δkl - Kkl
(21)
[
1
Xk
(mk +
]
∑l KkjXj) δkl - Xk-1KklXl
(22)
The sum of the off-diagonal elements of the kth row of E satisfies
|Ekl| ) Xk-1∑Xk-1KklXl < Xk-1(∑KklXl + KkkXk + mk) ∑ l*k l*k l (18)
2x2k1
, λ3 ) 0 (19) 3k1 + 1 + xk12 + 6k1 + 1
For k1 f ∞, |λ1| f 1/x2. Thus, undamped successive substitution converges more rapidly for the 3B scheme than for the schemes previously discussed. Second-Order Methods. Slow convergence of successive substitution, combined with the fact that the problem can be formulated as an unconstrained maximization, where global convergence can be guaranteed, makes second-order maximization-based methods an attractive alternative. This is in particular the case in the present situation, where the Hessian matrix is very simple and, therefore, inexpensive to evaluate. We have found it very advantageous to use the following iteration scheme:
H ˆ ∆X + g ) 0
(mk +
∑l KklXl + gk)
We show that H ˆ is negative definite as follows. Consider now the similarity transform
Ekl ) -
with eigenvalues
(λ1, λ2) ) (
Xk
Xk
S
(mk +
where D is the diagonal with diagonal elements Dkk ) Xk. E and H ˆ have identical eigenvalues, and the elements of E are
, k1 + 1 + xk12 + 6k1 + 1
X1 ) X 2 )
1
1
ˆD E ) D-1H
and
X3 )
)
(20)
where H ˆ is a modified Hessian matrix with the following properties: (i) It is negative definite for all X (ii) At the solution, H ˆ )H The first property implies that the direction ∆X will always be an ascent direction. Should the iteration fail in the sense that the objective function decreases when the full step is taken, a simple linesearch (e.g., by bisecting the step until an increase in Q is obtained) can be used as a remedy. The second property, on the other hand, ensures quadratic convergence.
That is, their sum is smaller than the magnitude of the diagonal element, and, therefore, by Gerschgorin’s theorem,6 all eigenvalues of H ˆ are negative. Note that this holds for any X. This is not the case with the original Hessian of eq 9. Take, for example, the 2B association scheme for a pure component and take as the initial estimate X1 ) X2 ) 1. Then, H becomes
(
-1 -k1 H ) -k 1 -1
)
with det(H) ) 1 - k12. If k1 > 1, one of the eigenvalues will be positive, and the matrix becomes indefinite. The second property of H ˆ , i.e., that it is equal to H at the solution, is evident because, at the solution, g is identical to zero. These properties suggest a very simple approach with guaranteed global convergence. Use an arbitrary initial estimate, and perform the following: (1) Calculate the step ∆X from eq 20. (2) Set Xnew ) Xold + ∆X. (3) Test that all elements of elements of Xnew are positive. A simple manner in which this can be obtained is by taking old old Xnew m ) max(Xm + ∆X, 0.20Xm )
(4) The value of the objective function is increased. (5) If any of the conditions are violated, set ∆X ) 1/2∆X and repeat from step 2. (6) Check for convergence. If not converged, set Xold ) Xnew and repeat from step 1. Practical Application The computational burden of solving the association equation is particularly high when molar density of the base equation of state must be determined. Typically, nested loops are used. The molar volume V (or the total volume V) is adjusted in an outer
8452
Ind. Eng. Chem. Res., Vol. 45, No. 25, 2006
loop, while the association equations must be solved for the matrix X corresponding to the assumed volume in the inner loop. When the inner loop has been converged, it is advantageous to calculate, in addition to the contribution to pressure from the association term,
Pasc ) -RT and, also, the derivatives calculated from
∂Pasc/∂V
X (∂Q ∂V )
(23)
and ∂X/∂V. The derivative is
∂g ∂X ∂g ∂g ) + X)0 ∂V ∂X V ∂V ∂V
( )
or
∂X ∂g + )0 ∂V ∂V
(24)
Here, H is already available in factored form from the solution for X. Finally,
( ) ()
∂Q ∂g ∂X -1 ∂P ) X+ 2 RT ∂V ∂V ∂V ∂V asc
2
T
The derivative of the association pressure is used in the outer loop to solve for the volume using a Newton-based method, and the volume derivative of X, which is obtained as a byproduct in the calculation of ∂Pasc/∂V, is used to create initial estimates for the inner loop solution of the association equations. When a correction ∆V has been determined in the outer loop at iteration n, we use, as an initial estimate for the inner loop,
X
(n+1)
)X
(n)
+ ∆V
(n)
1
10
102
103
104
105
106
103
2.97 3.13 3.58 3.93 5.20
2.97 3.21 3.62 4.04 5.16
2.98 3.18 3.61 4.02 5.09
3.00 3.19 3.49 3.93 5.01
3.19 3.06 3.48 3.91 4.96
3.48 3.48 3.75 3.89 5.09
3.93 3.91 3.89 4.35 5.43
104 105 106 107
k1/k2
1
10
102
103
104
105
106
103
5.44 5.90 6.01 6.56 6.86
5.75 5.78 6.32 6.76 6.96
5.26 5.99 6.65 6.83 7.29
5.26 5.99 6.53 6.89 7.59
5.84 5.88 6.57 6.93 7.84
6.02 6.05 6.13 6.89 7.99
6.04 6.23 6.66 6.94 7.96
104 105 106 107
by differentiation, with respect to V:
H
k1/k2
Table 2. Number of Iterations for the 3B-2B Scheme
g(X,V) ) 0
( )
Table 1. Number of Iterations for the 2B-2B Scheme
∂X (n) ∂V
( )
(25)
The tolerance for accepting an inner loop solution in the volume iteration is set fairly loose, and, as a consequence, only a single inner-loop iteration is necessary in most cases. This implies that the efficiency of the nested loop approach is only slightly lower than the more risky procedure of solving for V and X simultaneously. The safe procedure described earlier is only used in the few cases where a direct use of the Newton-based procedure fails to solve the association equations in a few iterations. Failures occur mostly when calculating liquid-phase properties for mixtures, where the association constants are very large. Here, the globally convergent procedure is used as a backup. Implementation of Safe Procedures We have implemented the back-up procedure in the following manner: (i) Initially, all elements of the X-vector are set to 0.2; (ii) A damping factor of ω ) 0.2 is used, and five steps of the substitution procedure are performed; and (iii) Subsequently, the second-order approach is used to converge the equations, with the modifications indicated earlier, i.e., (a) if an element of X is reduced by more than a factor of 5, the new value of that element is set to 0.2 times the previous value, and (b) a step that decreases the objective function is bisected until an increase in the objective function occurs.
We have tested the procedure as follows. We consider an equimolar binary mixture where both components can associate. The association schemes investigated are 2B-2B as well as 3B-2B. The reference association constants are called k1 and k2. Based on these, we create a sequence of “actual” association constants, which are calculated from
k11 ) r1k1, k22 ) r2k2, k12 ) r3xk11k22 where r1, r2, and r3 are evenly distributed random numbers in the range 0 < ri < 2. From this, we get the association matrices: For the 2B-2B association scheme:
(
0 k K ) 11 0 k12
(
k11 0 k12 0
0 k12 0 k22
k12 0 k22 0
)
For the 3B-2B association scheme:
0 0 K ) k11 0 k12
0 0 k11 0 k12
k11 k11 0 k21 0
0 0 k12 0 k22
k12 k12 0 k22 0
)
We assign values to the reference constants as shown in Tables 1 and 2. For each combination, 10 000 runs have been made, with random selection of the factors r1, r2, and r3. The equations are solved to machine accuracy (more than 10 significant digits). On the average, 3.9 iterations with the full second-order approach are required for the 2B-2B association scheme, whereas that for the 3B-2B association scheme is significantly higher (6.5). The reason for this is that, with high association constants and the 3B scheme, one of the elements of X becomes very small, and much smaller than is the case with the 2B scheme. A second indication of the 3B scheme being potentially more difficult is observed in the number of times that a stepsize reduction (bisection) is required, because the unmodified step leads to a decreased objective function. With the 2B-2B association scheme, this was never observed, whereas with the 3B-2B association scheme, this was observed, on the average, 0.67 times per case (i.e., in ∼10% of the iterations). The reason for the increase in the required number of bisections with the 3B scheme is probably that X3 becomes proportional to k1-1 (see eq 17), whereas with the 2B scheme, the site fraction only varies as k1-1/2. Consequently, the range of variation for the unknowns is orders of magnitude larger with the 3B association scheme.
Ind. Eng. Chem. Res., Vol. 45, No. 25, 2006 8453
The computation time for the 700 000 examples used in the two cases previously mentioned totaled 5 s, that is, ∼7 µs per calculation. A 3 GHz Pentium 4 computer with the Intel Fortran Compiler v. 9.0 was used for the runs. The damped successive substitution steps have a significant but modest effect. A reduction from five substitution steps to only a single substitution step increased the computation time by ∼30%. Xu et al.7 previously used interval arithmetic for computation of the phase stability with the SAFT equation of state, and, in their work, interval methods were also used for solving for the fractions of unbonded sites. The interval approach provides a mathematical guarantee that the correct solution is determined but carries a significant computational cost. Because the method used here is globally convergent and the maximum is unique, we obtain the same mathematical guarantee with a much simpler approach and at a cost that is orders of magnitude less. Solving for Volume A robust volume (or density) iteration can be implemented as follows. Use, as the independent variable, a reduced density ζ, which is the ratio of the actual density to that corresponding to the hard core volume. For CPA, we thus take ζ ) b/V. Therefore, the desired solution is in the range of 0 < ζ < 1. Set the initial limits as ζmin ) 0, ζmax ) 1. (1) Use, as the equation to solve, either F(ζ) ) P(ζ) - Pspec (SAFT) or F(ζ) ) (1 - ζ)(P(ζ) - Pspec) (CPA). (2) Assuming the hard core volume of the model to be b, take as initial estimates (i) ζ ) 0.5 (SAFT) or ζ ) 0.99 (CPA) if a liquid root is desired, or (ii) ζ ) 0 (SAFT) or ζ ) b/[b + (RT/P)] (CPA) if a vapor root is desired (3) Calculate, at step n, the current value of F and dF/dζ. An alternative to explicit calculation of dF/dζ is to approximate the derivative by ∆F/∆ζ, with ∆F ) F(n) - F(n-1) and ∆ζ ) ζ(n) - ζ(n-1). (4) If F > 0, set ζmax ) ζ; otherwise, set ζmin ) ζ. (5) Calculate a new value (ζnew) using Newton’s method:
ζnew ) ζ(n) -
F(k) dF/dζ(n)
(6) If ζmin < ζnew < ζmax, take ) ζnew. Otherwise, take ζ(n+1) ) (ζmin + ζmax)/2. (7) Continue until convergence is obtained. ζ(n+1)
Conclusion Procedures for solving the equation of state with association models have been developed. In particular, the use of a secondorder optimization method for solving the equations that describe the fractions of unbonded association sites ensures global convergence, combined with cost efficiency. Therefore, the complex and implicit nature of models based on the Wertheim association theory is not a hindrance for efficient and robust implementation. Nomenclature Ai ) type A site on molecule i b ) equation of state co-volume parameter in CPA Bj ) type B site on molecule j C ) number of components in mixture D ) diagonal matrix E ) transformed matrix fk ) function element
F ) target function for density solver gk ) element of gradient vector g ) gradient vector Hkl ) element of Hessian matrix H ) Hessian matrix i, j ) component indices Jkl ) element of Jacobian matrix k, l ) site indexes k1, k2 ) association constants Klm ) element of association matrix K ) association constant matrix ml ) amount of moles hosting site l n ) iteration counter ni, nj ) component molar amounts P ) pressure Q ) objective function r1, r2, r3 ) random factors, 0 < rl < 2 R ) gas constant S ) number of sites T ) temperature V ) molar volume V ) total volume XAi ) fraction of unbonded type A sites on molecule i X ) vector of unbonded site fractions Greek Symbols δkl ) Kronecker delta ∆AiBj ) association strength λ ) eigenvalue ω ) damping factor ζ ) packing fraction Literature Cited (1) 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. (2) Huang, S. H.; Radosz, M. Equation of state for small, large, polydisperse and associating molecules. Ind. Eng. Chem. Res. 1990, 29, 2284. (3) Kontogeorgis, G. M.; Voutsas, E. C.; Yakoumis, I. V.; Tassios, D. P. An equation of state for associating fluids. Ind. Eng. Chem. Res. 1996, 35, 4310. (4) Chen, C.-C.; Watanasiri, S.; Mathias, P.; De Leeuw, V. V. Industry perspective on the economic value of applied thermodynamics and the unmet needs of AspenTech clients. In Chemical Thermodynamics for Industry; Letcher, T. M., Ed.; Macmillan: London, England, 2004. (5) Michelsen, M. L.; Hendriks, E. M. Physical properties from association models. Fluid Phase Equilib. 2001, 180, 165. (6) Wilkinson, J. H. The Algebraic EigenValue Problem; Clarendon Press: Oxford, England, 1965. (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. 1998, 41, 938.
ReceiVed for reView January 9, 2006 ReVised manuscript receiVed February 3, 2006 Accepted February 3, 2006 IE060029X
