Reliable Computation of All the Density Roots of the Statistical

Mar 23, 2006 - Prediction of multiphase equilibrium using the PC-SAFT equation of state and simultaneous testing of phase stability. Nayef M. Alsaifi ...
0 downloads 0 Views 269KB Size
Download PDF
Ind. Eng. Chem. Res. 2006, 45, 3303-3310

3303

Reliable Computation of All the Density Roots of the Statistical Associating Fluid Theory Equation of State through Global Fixed-Point Homotopy Naveed Aslam and Aydin K. Sunol* Department of Chemical Engineering, UniVersity of South Florida, 4202 E. Fowler AVenue, Tampa, Florida 33620

Strong molecular basis imbedded into equations of state such as statistical associating fluid theory (SAFT) enhanced their popularity for phase behavior prediction of complex fluids and their mixtures despite added computational challenges. The compressibility expressions for the SAFT equation of state are seventh-order and ninth-order polynomials for pure components and mixtures, respectively. Reliable computation of all the possible solutions of the SAFT compressibility equation and identification of the physically meaningful solutions to use is a challenging but necessary task. A global fixed-point homotopy continuation based method which can reliably compute all the possible solutions of the SAFT equation of state is described. The global fixedpoint homotopy continuation guarantees that all the solutions of a nonlinear equation can be located on a single homotopy path when it is forced to start from a single starting point. This method is very desirable for the SAFT equation of state since the number and nature of density roots for a given temperature and pressure is not known a priori. The method shows that it is possible to locate all the solutions of the SAFT compressibility equation on a single path if a starting point is selected from a criterion which minimizes the number of real roots of the global fixed-point homotopy function. Examples covering both pure components and mixtures are presented to illustrate the methodology. 1. Introduction Recent advances in computer power and statistical mechanics accelerated the incorporation of molecular principles into classical thermodynamics and utility of more complex but accurate molecular equations of state. The statistical associating fluid theory (SAFT) is an equation of state (EOS) which is based on molecular interactions. The SAFT-EOS incorporates molecular associations (H-bonding, etc.), chain flexibility, and repulsive as well as dispersive interactions. Thus, small, large, polydisperse, and associating molecules and their mixtures can be accurately modeled over the entire density range. The capabilities of SAFT-EOS models are described well by Chapman and Gubbins,1 Huang and Radosz,2 and Fu and Sandler.3 SAFT has the following three parameters for nonassociating fluids. These are typically regressed from pure liquid component densities and saturated vapor pressures: number of segments in the molecule (mi); molar volume of segment (υi∞); segment interaction energy (ui0). For associating species, two additional parameters are required to characterize the association between different sites. The most widely used form of SAFT-EOS was developed by Huang and Radosz.2 The compressibility expression in this form of SAFTEOS is a combination of ideal compressibility, compressibility due to hard sphere contribution, compressibility due to chain contribution, and compressibility due to dispersion contribution. Thus, the SAFT compressibility equation becomes highly nonlinear in terms of reduced density. A highly nonlinear compressibility expression yields a large number of real roots as pointed out by Koak and de Loos,4 Lucia and Luo,5 and Brennecke and Stadtherr.6 For pure components, the SAFT compressibility equation becomes a seventh-degree equation, and for mixtures, it becomes a ninth-degree equation (Lucia and Luo5). Since the exact number of real compressibility roots * To whom correspondence should be addressed. Phone: +1-813 974 3566. Fax: +1-813 974 3651. E-mail: [email protected].

at a given temperature and pressure is not known a priori, a reliable method is needed to compute all the possible roots of the SAFT compressibility equation. Lucia and Luo5 described an approach for solving the nonlinear SAFT compressibility equation. Their approach is based on local methods such as Newton and dogleg, required a large number of initial guesses, and often diverges or exhibit periodic behavior. The dogleg method provides no guarantee that all the possible solutions can be computed. Brennecke and Stadtherr6 described a reliable approach to compute all the possible solutions of the SAFT compressibility equation through interval analysis. Interval methods although very reliable are difficult to implement in process simulators and require a special computing environment. Continuation based methods are another class of global methods, which, if implemented properly, lead to computation of all the possible solutions of nonlinear equations. This paper describes a global approach based on homotopy continuation to compute all the possible solutions of the SAFT compressibility equation. Our approach uses a global fixed-point homotopy continuation algorithm to solve the compressibility equation. The global fixed-point homotopy continuation method allows for reliable computation of all the possible solutions of a nonlinear equation from a single starting point. In addition, it ensures that all the possible real solutions are located on a single path provided that the starting point selection follows a certain criterion. Other homotopy formulations such as Newton and Affine require multiple starting points and paths to compute all the possible solutions of nonlinear equations. 2. Theory The compressibility expression of the SAFT-EOS for nonassociating systems where Dij are the Chen and Kreglewski6

Z - Zid - Zhs - Zdisp - Zchain ) 0.0

(2.1)

Zid ) 1.0

(2.2)

10.1021/ie0489554 CCC: $33.50 © 2006 American Chemical Society Published on Web 03/23/2006

3304

Ind. Eng. Chem. Res., Vol. 45, No. 9, 2006

4η - 2η2 (1 - η)3

Zhs ) m 4

Zdisp ) m

(2.3)

[ ][ ]

9

u

∑ ∑ jDij kT i)1.0 j)1.0

i

η

j

(2.4)

τ

constants, k is Boltzmann’s constant, and τ ) πx2/6 is the packing fraction for closed-packed spheres in eq 2.4. For mixtures, the average segment number m, average segment energy u, and average reduced density (packing fraction) η are calculated from following equations. N

m)

η) N

u

)

kT

ximi ∑ i)1

(πNavF)

(2.5)

N

ximidi3 ∑ i)1.0

6

(2.6)

N

xixjmimjυij0(uij/kT) ∑ ∑ i)1.0 j)1.0 N

(2.7)

N

∑ ∑ xixjmimjυij0

i)1.0 i)1.0

uij ) (1 - kij)(uiuj)0.5

(2.8)

υij ) 1/2((υi0)1/3 + (υj0)1/3)3

(2.9)

(

ui ) ui0 1 + Zchain ) (1 - m)

[

e kT

)

(2.10)

(5/2)η - η (1 - η)(1 - (1/2)η) 2

]

compound

m

υ∞ (mL/mol)

u0/k (K)

e/k (K)

propane nitrogen cyclohexane carbondioxide ethane methane

2.69 1.0 3.97 1.0 1.0 1.0

13.45 19.45 13.52 19.81 31.07 21.57

193.03 123.53 236.41 274.37 301.38 190.29

34.0 3.0 10.0 52.0 16.0 1.0

Another important difference between newton and fixed-point homotopy is the range over which the homotopy parameter is varied. For newton homotopy, usually the parameter is within the range of 0-1. However, for global fixed-point homotopy, it is not restricted to the range of 0-1, rather it is allowed to attain any value. Although global newton homotopy provides excellent scaling of variables, it does not guarantee that all the real solutions of a given nonlinear equation are located on a single path starting from a single point. Multiple continuation paths that start from different starting points are required to compute all the real solutions through the global newton homotopy approach. This may not be a very desirable situation in case of the SAFT compressibility equation where the number of real roots at a given temperature and pressure is not known a priori. Kuno and Seader8,9 showed that it is possible to locate all the real roots of a nonlinear equation on a single path provided that the homotopy path is forced to start from a certain point. The starting point is computed through application of mathematical criteria which ensure that the homotopy path will not become unbounded. They demonstrated that for certain type of functions a starting point x0 can be selected by placing t f ∞ in eq 3.1. For t f ∞, eq 3.1 for fixed-point global homotopy becomes

F(x) - x + x0 ) 0 (2.11)

where d is the diameter of spheres and F is the molar density in eq 2.6. After placing expressions from equations 2.2-2.11 in eq 2.1, it becomes a nonlinear equation in terms of density (F) or reduced density (η). 3. Global Fixed-Point Homotopy Homotopy continuation methods are globally convergent methods for the solution of nonlinear equations; i.e., these methods can find the solution while moving from a known solution or a starting point to an unknown solution. The homotopy function consists of a linear combination of two functions G(x), an easy function, and F(x), a difficult function.

H(x,t) ) tF(x) + (1 - t)G(x)

Table 1. Pure Component Parameters for the SAFT Equation of State

(3.1)

where t is a homotopy parameter that is gradually varied from 0 to 1 as a path is tracked from G(x) (known solution) to F(x) (unknown solution). The function G(x) in eq 3.1 can take different forms which lead to different types of homotopies. Two main types are newton and fixed-point homotopy as described by Wayburn and Seader.7

Newton homotopy G(x) ) F(x) - F(x0)

(3.2)

Fixed-point homotopy G(x) ) x - x0

(3.3)

(3.4)

Thus, a starting point selected through solving eq 3.4 and which yields the minimum number roots of x0 will prevent the homotopy parameter from becoming unbounded and ensures that all real roots are located on a single path. 4. Result and Discussions Through global fixed-point homotopy, we are able to compute all the roots of the SAFT compressibility equation (eq 2.1) both for pure components and mixtures. The data for pure component SAFT parameters is presented in Table 1. We tested our approach with several combinations of temperature and pressure. To solve eq 2.1 through global fixed-point homotopy, a starting point criterion is applied. Figure 1 depicts the graph of eq 3.4 for pure nitrogen at 95.0 K and 2.0 bar, i.e., application of the starting point criterion on the SAFT compressibility equation. The general F(x) function in eq 3.4 is actually eq 2.1 for a specific case of SAFT compressibility. Thus, Figure 1 provides insights into the selection of a starting point for a fixedpoint global homotopy branch. As stated earlier, the starting point should be selected from a range which yields the minimum number of roots for eq 3.4. The three subgraphs of Figure 1 are included to illustrate the effect of range. As can be seen from Figure 1, below -800.0 or above 15.0, eq 3.4 for pure nitrogen has only one root, so we selected 15.0 as the starting point for this case. Figure 2 contains a homotopy plot for the SAFT compressibility equation for pure nitrogen at 95 K and 2.0 bar. Five real roots are located on a single path. Two additional complex roots are also located on a path starting from the

Ind. Eng. Chem. Res., Vol. 45, No. 9, 2006 3305

Figure 1. Application of starting point criterion to SAFT compressibility eq 2.1 for nitrogen at P ) 2.0 bar and T ) 95 K.

Figure 2. Homotopy path for SAFT compressibility eq 2.1 for nitrogen at 2.0 bar and 95 K.

complex domain. Table 2 contains the real roots for nitrogen. The complex roots are 0.3667 - 0.0108i and 0.3640 + 0.1033i, respectively. Figure 2 also contains the complex homotopy branches. Since the initial point is based on a criterion, the

homotopy path will not become unbounded. The homotopy branch is tracked with a pseudo-arc-length continuation algorithm as described in our earlier work13 and by Allgower and Georg.12 The selection of step size is very important since sharp

3306

Ind. Eng. Chem. Res., Vol. 45, No. 9, 2006

Figure 3. Application of starting point criterion to SAFT compressibility eq 2.1 for carbon dioxide at P ) 25.0 bar and T ) 278.0 K. Table 2. SAFT Compressibility Roots Obtained through Global Fixed-Point Continuation component

temperature (K)

pressure (bar)

reduced density

compressibility

nitrogen nitrogen nitrogen nitrogen carbon dioxide propane cyclohexane cyclohexane cyclohexane methane

95 80 125 142 278 250 340 380 460 150

2.0 2.0 30.0 80.0 25.0 3.0 0.64 2.07 10.08 10.0

0.006, 0.096, 0.365, 0.647, 0.918 0.008, 0.052, 0.406, 0.623, 0.931 0.264, 0.706, 0.883 0.273, 0.752. 0.852 0.044, 0.092, 0.323, 0.669, 0.905 0.007, 0.047, 0.352, 0.748, 0.851 0.001, 0.0472, 0.327 0.004, 0.054, 0.299 0.027, 0.058, 0.243 0.034, 0.070, 0.353, 0.654, 0.914

0.924, 0.063, 0.016, 0.009, 0.006 0.852, 0.138, 0.017, 0.011, 0.0078 0.263, 0.098, 0.0707 0.595, 0.216, 0.191 0.602, 0.288, 0.082,0.039, 0.029 0.856, 0.136, 0.018, 0.008, 0.0076 1.507, 0.0319, 0.004 1.091, 0.080, 0.014 0.650, 0.302, 0.072 0.622, 0.303, 0.060, 0.032, 0.023

turning points are encountered and if not selected correctly will cause a failure in branch tracking. For this case, the step size of 0.005 is used. For pure nitrogen at 95 K and 2.0 bar, the vapor compressibility is given by the largest root 0.924. On the other hand, liquid compressibility corresponds to a root which has the lowest chemical potential also pointed out by Koak and de Loos.4 That will correspond to a compressibility root of 0.006 for liquid nitrogen. Since the corresponding reduced density η is 0.918 for this root, which is greater than the nominal closed-packed reduced density of 0.74048, this is not a physically meaningful compressibility root for the liquid phase. Thus, the next compressibility root of 0.009 that corresponds to a reduced density of 0.647 will be the correct liquid-phase root under these conditions. Table 2 also contains the roots of the SAFT compressibility equation for nitrogen under other conditions of temperature and pressure. At 142 K and 80.0 bar, there is only one physically meaning full compressibility root, i.e., 0.595 corresponding to reduced density of 0.2733. The other two roots under this condition are corresponding to reduced densities of 0.752 and

0.852, respectively, which are greater than the nominal closedpacked reduced density of 0.74048. Thus, a compressibility of 0.595 describes nitrogen, which is in the supercritical state at this temperature. Figure 3 shows the graph of eq 3.4 for pure carbon dioxide at 25.0 bar and 278.0 K. The starting point for this case is selected as -28.0. Figure 4 contains the homotopy branch for pure carbon dioxide at 25.0 bar and 278.0 K. Five real roots are located on a single path starting from a single point. Table 2.0 contains all five real roots of the SAFT compressibility equation for pure carbon dioxide at 25.0 bar and 278.0 K. The vapor- and liquid-phase compressibility roots in this case are 0.602 and 0.039, respectively. Table 2 also contains the compressibility roots of the SAFTEOS for propane, cyclohexane, and methane through global fixed-point homotopy continuation. Five real and two complex roots are obtained for propane at 3.0 bar and 250 K. Only real roots are reported in Table 2.0 as complex roots do not carry any physical significance. The vapor-phase compressibility root for propane at 3.0 bar and 250 K is 0.856, while the liquidphase compressibility root is 0.018.

Ind. Eng. Chem. Res., Vol. 45, No. 9, 2006 3307

Figure 4. Homotopy path for SAFT compressibility eq 2.1 for carbon dioxide at 25.0 bar and 278 K.

Figure 5. Application of starting point criterion to SAFT compressibility eq 2.1 for a mixture of carbon dioxide and ethane (40% carbon dioxide and 60% ethane) at 17.8 bar and 263.15 K.

The SAFT compressibility roots for cyclohexane are computed at three different sets of temperature and pressure values. For 340 K and 0.64 bar, only three real roots are obtained. For

the other two conditions (380 K and 2.07 bar and 460 K and 10.08 bar), again, only three sets of real roots are obtained. Table 2 also contains all five real roots of pure methane at

3308

Ind. Eng. Chem. Res., Vol. 45, No. 9, 2006

Figure 6. Homotopy path for SAFT compressibility eq 2.1 for a mixture of carbon dioxide and ethane (40% carbon dioxide and 60% ethane) at 17.8 bar and 263.15 K.

Figure 7. Real and complex homotopy path for SAFT compressibility eq 2.1 for a mixture of carbon dioxide and ethane (40% carbon dioxide and 60% ethane) at 17.8 bar and 263.15 K.

150 K and 10.0 bar. The vapor- and liquid-phase compressibility roots for pure methane at these conditions are 0.622 and 0.032. We also calculated the compressibility roots for mixtures. Figure 5 contains the graph generated through application of starting point criterion to the SAFT compressibility equation (eq 2.1) for a mixture of carbon dioxide and ethane (40% carbon dioxide and 60% ethane) at 17.80 bar and 263.15 K. The binary

interaction parameter for this mixture is set at kij ) 0.15. Figure 6 contains a global fixed-point homotopy branch of the SAFT compressibility equation for this mixture. Five real roots are located on a single path, and four complex roots are also found. Only the real path is shown in Figure 6. Figure 7 contains both the real and complex paths of the compressibility roots. The complex roots for this system are 0.1465 - 0.1248i, 0.1430 +

Ind. Eng. Chem. Res., Vol. 45, No. 9, 2006 3309

Figure 8. Application of starting point criterion to SAFT compressibility eq 2.1 for a mixture of propane and cyclohexane (10% propane and 90% cyclohexane) at 1.0 bar and 215 K.

Figure 9. Homotopy path for SAFT compressibility eq 2.1 for a mixture of propane and cyclohexane (10% propane and 90% cyclohexane) at 1.0 bar and 215 K.

0.1248i, 0.3127 + 0.008i, and 0.1096 + 0.0001i. The step size used in the continuation algorithm is 0.00018. Figures 8 and 9 contain the graph of the application of starting

point criteria and the global fixed-point homotopy branch for a propane-cyclohexane (10% propane and 90% cyclohexane) mixture at 1.0 bar and 215 K. The binary interaction parameter

3310

Ind. Eng. Chem. Res., Vol. 45, No. 9, 2006

for this mixture is set at kij ) 0.1. The starting point of -8.5 is used in this case. Five real roots are located on a single path. The step size used for this example is 0.0025.

υi0 ) molar volume of segments. Z ) compressibility factor. Literature Cited

5. Conclusion All the roots of SAFT compressibility equations are reliably computed through global fixed-point homotopy continuation. Global fixed-point homotopy provides a robust approach to compute all the real roots of a single nonlinear equation without any prior information on the number and nature of roots. All the real roots of a nonlinear equation can be located on a single path if the homotopy path is forced to start from a starting point which is computed by a mathematical criterion. The mathematical criterion for selecting the starting point is based on transforming the general homotopy equation with t f ∞ and selecting a starting point from a range of values which yields the minimum number of real roots. Global fixed-point homotopy can also calculate complex roots of the nonlinear compressibility equation. The SAFT compressibility equation is solved through global fixed-point continuation both for pure components and mixtures at different temperatures and pressures. Nomenclature η ) reduced density µ ) chemical potential τ ) close-packed reduced density ) 0.74048 F ) molar density d ) temperature-dependent density diameter of a segment. e/k ) pure component constant (K) k ) Boltzmann constant m ) number of segments P ) pressure T ) temperature u/k ) temperature-dependent energy dispersion term u0 ) temperature-independent segment-segment energy

(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) Fu, Y.; Sandler, S. I. A Simplified SAFT Equation of State for Associating Compounds and Mixtures. Ind. Eng. Chem. Res. 1995, 34, 1897. (4) Koak, N.; de Loos, T. W.; Heidemann, R. A. Effect of the Power series Dispersion Term on the Pressure-Volume Behavior of Statistical Associating Fluid Theory. Ind. Eng. Chem. Res. 1999, 38, 1718. (5) Lucia, A.; Luo, Q. Binary Refrigerant-Oil Phase Equilibrium Using the Simplified SAFT Equation. AdV. EnViron. Res. 2001, 6, 123. (6) Xu, G.; Brennecke, J. F.; Stadtherr, M. A. Reliable Computation of Phase Stability and Equilibrium from SAFT Equation of State. Ind. Eng. Chem. Res. 2002, 41, 938. (7) Wayburn, T. L.; Seader, J. D. Homotopy Continuation Methods for Computer-Aided Process Design. Comput. Chem. Eng. 1987, 11, 7. (8) Seader, J. D.; Kuno, M.; Lin, S. A.; Johnson, K.; Wiskin, J. W. Mapped Continuation Methods for Computing All Solutions to General System of Nonlinear Equations. Comput. Chem. Eng. 1990, 14, 71. (9) Kuno, M.; Seader, J. D. Computing All Real Solutions to Systems of Nonlinear Equations with a Global Fixed-Point Homotopy. Ind. Eng. Chem. Res. 1988, 27, 1320. (10) Choi, S. H.; Book, N. L. Unreachable roots for global homotopy continuation methods. AIChE J. 1991, 37, 1093. (11) Gritton, K. S.; Seader, J. D.; Lin, W. J. Global Homotopy continuation procedure for seeking all roots of a nonlinear equation. Comput. Chem. Eng. 2001, 25, 1003. (12) Allgower, E. L.; Georg, K. Numerical continuation methods: An introduction; Springer: Berlin, 1990. (13) Aslam, N.; Sunol, A. K. Reliable Computation of Binary Homogeneous Azeotropes of Multi-Component mixtures at Higher Pressures through Equations of State. Chem. Eng. Sci. 2004, 59, 599.

ReceiVed for reView October 27, 2004 ReVised manuscript receiVed December 5, 2005 Accepted January 24, 2006 IE0489554