Simulation of Refrigerant Phase Equilibria - American Chemical Society

Jun 1, 1997 - The complex domain trust region methods of Lucia and co-workers (Lucia ... trust region algorithms outperform Newton's method in singula...
0 downloads 0 Views 267KB Size
Download PDF
Ind. Eng. Chem. Res. 1997, 36, 2841-2848

2841

Simulation of Refrigerant Phase Equilibria Arthur S. Gow Department of Chemistry & Chemical Engineering, University of New Haven, 300 Orange Avenue, West Haven, Connecticut 06516-1999

Xinzhou Guo Department of Chemical Engineering, Clarkson University, Potsdam, New York 13699-5705

Delong Liu and Angelo Lucia* Department of Chemical Engineering, University of Rhode Island, Kingston, Rhode Island 02881-0805

Vapor-liquid equilibria for refrigerant mixtures modeled by an equation of state are studied. Phase behavior calculated by the Soave-Redlich-Kwong (SRK) equation with a single adjustable binary interaction parameter is compared with experimental data for binary refrigerant mixtures, two with a supercritical component and one that exhibits azeotropic behavior. It is shown that the SRK equation gives an adequate description of the phase envelope for binary refrigerant systems. The complex domain trust region methods of Lucia and co-workers (Lucia and Xu, 1992; Lucia et al., 1993) are applied to fixed vapor, isothermal flash model equations, with particular attention to root finding and root assignment at the equation of state (EOS) level of the calculations, and convergence in the retrograde and azeotropic regions of the phase diagram. Rules for assigning roots to the vapor and liquid phases in the case where all roots to the EOS are complex-valued are proposed and shown to yield correct results, even in retrograde regions. Convergence of the flash model equations is also studied. It is shown that the complex domain trust region algorithms outperform Newton’s method in singular regions of the phase diagram (i.e., at near azeotropic conditions and in the retrograde loop), primarily due to the eigenvalueeigenvector decomposition strategy given in Sridhar and Lucia (1995). A variety of geometric figures are used to illustrate salient points. 1. Introduction The modeling and simulation of vapor-liquid equilibria (VLE), liquid-liquid equilibria (LLE), and other multiphase phenomena are important in the analysis, design, operation, and control of many chemical processes. Oil recovery, high pressure and supercritical separations, and the production of environmentally safe refrigerants are examples of processes which require both accurate descriptive phase behavior models and reliable tools for solving the model equations. Models such as the Soave-Redlich-Kwong (SRK) equation (Soave, 1972) and the Peng-Robinson (PR) equation (Peng and Robinson, 1976) with conventional polynomial mixing rules and a single adjustable binary interaction parameter have been quite successful in correlating and extrapolating phase equilibria for mixtures composed of nonpolar or slightly polar components but have been observed to perform poorly when applied to mixtures containing highly polar or associating components and to systems at harsh processing conditions. Efforts to overcome mixing rule deficiencies have focused on the composition dependence of the equation of state (EOS) parameters. The most notable approach of this type is the development of mixing rules by combining an equation of state with an excess Gibbs free energy model which was pioneered by Huron and Vidal (1979) and refined by others (Dahl and Michelsen, 1990; Heidemann and Kokal, 1990; Wong and Sandler, 1992; Huang and Sandler, 1993) and is reviewed in detail elsewhere (Peng et al., 1995). While excess free* Author to whom correspondence should be addressed: tel, (401) 874-2814; fax, (401) 782-1180; e-mail, [email protected]. S0888-5885(97)00043-2 CCC: $14.00

energy-based mixing rules have been generally quite successful, there are a number of consistency issues which raise concern and are currently being addressed. Equation of state models are often used within flash equations to simulate all types of equilibrium flash processes, particularly at high pressure. Many equation-solving algorithms have been proposed for this purpose over the years, while more recent developments have concentrated on alleviating convergence problems (Michelsen, 1982a,b; Trangenstein, 1987; Whitson and Michelsen, 1989; Eubank et al., 1992, Sun and Seider, 1992; Schnepper et al., 1994; Heidemann and AbdelGhanni, 1994; McDonald and Floudas, 1995). These newer algorithms use numerical techniques such as tangent plane criteria, homotopy-continuation, interval arithmetic and mixed integer nonlinear programming (MINLP), as well as others. None of these papers, however, mentions the interrelated issue of finding correct roots to the underlying equation of state to any great extent. Most simply take for granted that the correct roots to the equation of state are readily available, but this is often not the case. A novel approach for solving phase equilibrium problems, which uses numerical methods that operate in the complex domain, was suggested by Lucia et al. (1990), developed by Lucia and Xu (1992) and Lucia et al. (1993), and applied to various models involving phase equilibria by Lucia and co-workers (Lucia and Taylor, 1992; Lucia and Wang, 1995; Taylor et al., 1996). They show that modeling in the complex domain provides smoothness across phase boundaries under very mild conditions, that real-valued solutions can be computed from complex-valued starting points, and that convergence to complex-valued solutions provides useful physi© 1997 American Chemical Society

2842 Ind. Eng. Chem. Res., Vol. 36, No. 7, 1997

cal information regarding phase existence. Numerical results clearly show that complex domain trust region methods are very reliable, often providing a way of solving problems which are difficult to solve using algorithms that are confined to the real domain. However, these papers do not address in detail issues such as root assignment when all roots to the EOS are complex-valued (which can happen when the flash temperature, pressure, and/or phase composition iterates become complex) or algorithmic behavior in singular regions of the phase diagram (i.e., near azeotropic points or in retrograde regions). These numerical issues, together with the correct modeling of refrigerant mixture vapor-liquid equilibria, are the focus of this paper. Accordingly, the paper is organized in the following way. Mathematical models for the equation of state and flash vessel are presented in section 2. In section 3, equation-solving and the complex domain methods of Lucia and co-workers are summarized. Particular attention is given to the root finding and root assignment tasks at the EOS level of the calculations, and rules for the assignment of compressibility factor roots to particular phases when all roots are complex-valued are proposed. Numerical results for a number of binary refrigerant mixtures are presented in section 4, with a focus on computational difficulties associated with EOS root finding and assignments and the presence of various singularities (trivial solutions, azeotropy, and retrograde behavior). General conclusions are given in section 5. 2. Mathematical Models In this section, models for the Soave-Redlich-Kwong (SRK) equation of state (EOS) for vapor and liquid phases and multiphase flash are presented. 2.1. Modeling Vapor and Liquid Phases. The SRK equation of state (Soave, 1972),

P)

aM(T) RT VM - bM VM(VM + bM)

(1)

is one of many simple cubic EOS that can be used to model the behavior of either vapor or liquid. Here, aM(T) and bM are fundamental mixture parameters whose composition dependence can be obtained from any number of mixing rules discussed earlier. In this paper, however, we apply the classical polynomial mixing rules

parameters are given by

ai(T) )

( x )]

[

ai(Tci) 1 + (0.480 + 1.574ωi - 0.176ωi2) 1 -

T Tci

2

(5) with

0.42748R2Tci2 ai(Tci) ) Pci

(6)

and

0.08664RTci Pci

bi )

(7)

where Tci, Pci, and Vci are the component i critical temperature, pressure, and molar volume, respectively, and ωi is the acentric factor of component i. The corresponding fugacity coefficient for component i is given by

(

ln φi ) ln

VM

VM - bM nc

2

)

bi +

VM - bM

zjaij(T) ∑ j)1

(

ln

)

VM - bM

+ RTbM VM VM - bM bM aM(T)bi ln - ln ZM (8) VM VM - bM RTb 2 M

[(

)

]

where φiL and φiV are the fugacity coefficients of species i in the liquid and vapor and where the subscript M denotes properties of the mixture. 2.2. Flash Model Equations. The usual equations that model a typical vapor-liquid equilibrium flash are the phase equilibrium equations

vi nc

∑i vi

) Ki

li nc

(9)

∑i li

nc nc

aM(T) )

∑i ∑j aij(T)zizj

(2)

fi ) li + vi

nc

bM )

∑i bizi

(3)

with the conventional corrected geometric mean combining rule for aij(T)

aij(T) ) (1 - kij)xai(T)aj(T)

and the component mass balance equations

(4)

where zi is either the liquid phase mole fraction of component i, xi, or the vapor phase mole fraction of component i, yi, and kij is an adjustable parameter which symbolizes the inability of the geometric mean combining rule in eq 4 to truly describe the molecular interactions between species i and j. The pure component

(10)

where fi, li, and vi are the ith component molar flows of the feed, liquid product, and vapor product, respectively, and where Ki is given by the ratio of φiL to φiV. Furthermore, depending on the flash specifications, an energy balance equation may also constitute part of the model. Here we consider the solution of isothermal, fixed vapor fraction (VT) flash problems, which are generalized bubble-point (or dew-point) pressure calculations, so instead of an energy balance, we use the specification equation nc

vi ) V ∑ i)1

(11)

Ind. Eng. Chem. Res., Vol. 36, No. 7, 1997 2843

Equations 9, 10, and 11 must be solved simultaneously using the equation of state for the unknown variables P, li, and vi. Note that the vapor-liquid phase diagram for an isothermal binary system may be generated by successively solving this problem over a range of feed compositions and by fixing the vapor to feed ratio (V/F) to zero so that the feed is in equilibrium with the vapor product (i.e., this gives the bubble-point pressure curve) and/or by fixing V/F to 1 so that the feed is in equilibrium with the liquid product (this gives the dew-point pressure curve). 3. Equation Solving In this section, mathematical preliminaries and computational difficulties associated with equation-solving methods in the real and complex domains are discussed and the complex domain algorithms of Lucia and coworkers are presented. 3.1. Mathematical Preliminaries. The flash model equations presented in the previous section can be written in the compact form

F(Z,p) ) 0

(12)

where F is a function defined on the complex domain and where Z and p denote complex-valued vectors of unknown variables and parameters of lengths n and m, respectively. Equation 12 can be recast in the form

Z ) G(Z,p)

(13)

where G is a nonlinear iterative map that depends on F and is therefore also defined on the complex domain. Moreover, it is well-known that solutions or roots to eq 12, say Z*, correspond to fixed points of G (i.e., Z such that Z* ) G(Z*,p)). For example, Newton’s method has the fixed-point form

Zk+1 ) Zk - [J(Zk,p)]-1F(Zk,p) ) GN(Zk,p) (14) where J(Z,p) denotes the Jacobian matrix of F evaluated at Z and p, and where k denotes an iteration counter. The subscript N denotes Newton’s method. Note that the roots, solutions, or fixed points, as well as the elements of the Jacobian matrix, can be real- or complex-valued. Other equation-solving methods (e.g., direct substitution and tearing algorithms) can also be written as nonlinear iterative maps, although the actual functionality of G is sometimes implicit. Equation-solving methods like Newton’s method can exhibit both periodic and chaotic behavior from real- and complex-valued starting points, even on seemingly simple nonlinear problems (see, for example, Cayley, 1879; Fatou, 1919; Julia, 1918; Curry et al., 1983). If iterates behave periodically or chaotically, then so does the two-norm (least-squares function) or any other norm of F. Stabilization procedures like line searching and trust region methods, which force norm reduction at each iteration, avoid periodic and chaotic behavior. However, they do so at the risk of terminating at singular points which often correspond to local minima in the two-norm of F in the real domain (Lucia and Xu, 1992). A singular point is any point of F at which the Jacobian matrix of F is rank-deficient (has a zero determinant or one or more zero eigenvalues). 3.2. Complex Domain Algorithms. The heart of the complex domain trust region methods of Lucia and co-workers (i.e., the single variable algorithm of Lucia and Xu (1992) and the multivariable method of Lucia

et al. (1993)) is the so-called “dogleg” strategy of Powell (1970), which is globally convergent to stationary points (i.e., local and global minima) of the least-squares function of F in the real domain. Local minima in the real domain correspond to singular points of F. This guaranteed convergence to either a solution or a singular point of F is achieved by using a linear combination of Newton and steepest descent steps, norm reduction in F at each iteration, and an iteratively adjusted trust region that forces the algorithm to take steepest descent steps when serious difficulties are encountered. Steepest descent guarantees reduction in the norm of F for small trust region sizes and thus global convergence. However, it is important to emphasize that global convergence here does not mean guaranteed convergence to a solution! It means convergence to either a solution or a singular point. This is a subtle and often misunderstood point. We refer the reader to the original paper by Powell (1970) for a more complete treatment of the basic steps of the algorithm, the “rules” for both selecting Newton, dogleg, and steepest descent steps as well as adjusting the size of the trust region iteratively, and the global convergence proofs. We use essentially the same rules in our complex domain trust region algorithms, which are really just complexified versions of the dogleg strategy with enhancements to avoid termination at singular points. The complex domain trust region algorithms work in the following way and represent a matriculation of ideas contained in Lucia et al. (p 649, 1990), Lucia and Xu (1992), Lucia et al. (1993), and Sridhar and Lucia (1995). The algorithm functions the same way as the original dogleg strategy, only in the complex domain, given a real- or complex-valued starting point and terminates when a real- or complex-valued solution is found. It is only when a singular point is encountered that our algorithms differ from that of Powell. In particular, an eigenvalue-eigenvector decomposition of the second derivative (or Hessian) matrix of the two-norm of F is performed at the singular point. Moreover, this Hessian matrix has at least one negative eigenvalue and, of course, an associated eigenvector because any singular point is a saddle point of F. The negative eigendirection and continuity guarantee that there is a point in the neighborhood of the singular point whose two-norm of F is less than that at the singular point. Once this point of lower norm is found, the algorithm returns to the complex domain dogleg strategy, terminating at either a solution or singular point, and this alternation between the eigenvalue-eigenvector decomposition and the complex domain dogleg repeats itself until a solution is found. 3.3. Finding Roots to an Equation of State. One of the most important aspects of solving phase equilibria modeled by a single equation of state is the assignment of a qualitatively correct compressibility root to each of the equilibrium phases. An incorrect root assignment (e.g., assigning a vapor-like root to a liquid phase, or vice versa) can lead to erroneous thermodynamic properties (fugacities, fugacity coefficients, K values, etc.) and is often computationally catastrophic, causing convergence to a trivial solution or even divergence. In either case, this is viewed as algorithmic failure. However, correct root assignments are complicated by the fact that some (or all) of the roots to the EOS can be complex-valued, particularly at high pressure, and by the nonlinearity and fractal nature of the basins of attraction for the roots. That is, widely different start-

2844 Ind. Eng. Chem. Res., Vol. 36, No. 7, 1997

ing points are not necessarily guaranteed to find different roots. For example, the standard compressibility factor initializations of ZM ) bMP/RT for the liquid phase and ZM ) 1 for the vapor phase, or any other rational choices, will not always lead to the physically meaningful liquid and vapor roots. One way to guarantee that correct root assignments are made is to compute all roots for each phase present in the model. Thus, given any phase type (i.e., liquid or vapor), an unambiguous decision regarding the compressibility of that phase can be made. Moreover, this is true regardless of whether the desired root is real- or complex-valued. The strategy we use to find all roots to the SRK (or any other) cubic EOS is straightforward and consists of finding one root numerically using a single variable complex domain trust region method, followed by deflation (or factorization) of the cubic to give a quadratic from which the remaining two roots can be calculated analytically by application of the quadratic formula. This removes the need for a “good” initial guess for finding the first root in either phase because regardless of the root found, the other two are easily computed using the quadratic formula. In practice, we actually initialize the numerical liquid compressibility calculations by setting ZM ) bMP/RT and initialize the vapor compressibility by setting ZM ) 1. 3.4. Root Assignments. Compressibility factor roots for vapor-liquid equilibria (VLE) are assigned in the following simple manner. ZM for the vapor is chosen to be the root with the largest positive real part of the three roots computed from T, P, and y. Conversely, ZM for the liquid is chosen to be the root with the smallest real part of those roots computed from T, P, and x. It is important to note that the assignment of distinctly different compressibility roots annihilates the existence of VLE trivial solutions. The rules are similar for liquid-liquid equilibria (LLE). In particular, ZM for the first liquid phase, liquid 1, is chosen to have the smallest real part of the roots computed using T, P, and x1 (the composition of liquid 1). Similarly, ZM for the second liquid phase (i.e., at T, P, and x2) is the root with the smallest positive real part. Clearly, other rules are easily developed for vapor-liquid-liquid equilibria (VLLE) and other multiphase systems. 4. Numerical Results Three representative binary systems composed from the pure components chlorodifluoromethane (R22), trifluoromethane (R23), methane (R50), ethane (R170), carbon dioxide (R744), and ethene (R1150) are used to illustrate the vapor-liquid equilibria generally observed for binary refrigerant mixtures. Note that all but the first are environmentally-friendly refrigerants since they contain no chlorine. The algorithms discussed in section 3 were implemented in Fortran 77, and all numerical results given in this section were computed to an accuracy of 10-8 in the two-norm of the model functions and done on an IBM 586 in double precision arithmetic. The compressibility calculations were initialized as described at the end of section 3.3, and the flash variables (i.e., x, y, L/F, and P) were initialized as follows unless stated otherwise. For V/F ) 0 and a feed composition of z, we set x ) z, L/F ) 1, and choose P and y arbitrarily based on assumed relative volatilities. On the other hand, for V/F ) 1, we set y ) z, L/F ) 0, and choose P and x arbitrarily based on assumed relative volatilities. 4.1. Roots to the Equation of State. The SRK EOS with a value of k12 ) 0 has been used to correlate

Figure 1. Iterative compressibilities for isothermal VLE for R23(1)-R22(2) at 343.15 K using the SRK equation.

the experimental isothermal VLE data of Roth et al. (1992) for R23(1)-R22(2) at 343.15 K, for which R23 is supercritical. Consider the tie line corresponding to a pressure of 41.6514 bar. We solved a single-stage flash problem with specified total vapor flow and temperature (a VT flash) in the complex domain to calculate the VLE represented by this tie line. Any value of V, 0 e V e 1, may be used; V ) 0 was chosen to make the problem equivalent to a bubble-point pressure calculation. All roots for the SRK equation for both the vapor and liquid phases (i.e., two sets of three in all) were determined at each flash iteration, from which the appropriate vapor and liquid compressibilities were chosen using the rules given in section 3.4. Figure 1 shows the vapor and liquid compressibilities for each of the six iterations needed to reach a convergence tolerance of 10-8 in the two-norm of the flash model equations. Note that the initial compressibility selected for the liquid phase is complex-valued while the one for the vapor is real but that the converged compressibilities for both phases are real-valued. Figures 2 and 3 show the level sets (or contours) of the two-norm of the SRK equation for the liquid and vapor phases for the second iteration of the flash (or bubble-point) calculations and give a good geometric illustration of the EOS root-finding task. It is important to note that the vapor composition and pressure iterates for this second iteration are complex-valued. This makes the coefficients of the EOS complex-valued. As a result, the SRK equation has no real-valued roots for either phase on the second iteration, which seems impossible at first glance. Note, however, that the critical compressibility predicted by the SRK equation (i.e., Zc ) 1/3) provides a clear delimiter between vaporlike and liquid-like roots. That is, any root, real- or complex-valued, whose real part is less than 1/3 can be considered liquid-like, while any root with a real part greater than 1/3 can be considered vapor-like. This is clearly illustrated in Figures 2 and 3 in which there are two liquid-like roots and one vapor-like root in each figure. Keep in mind that there are no real-valued roots in this case. It is also important to understand that we do not actually compare the real parts of the roots to 1/3 in determining whether they are liquid or vapor. We use the root assignment rules defined in section 3.4 exactly as stated and select the compressibility with the

Ind. Eng. Chem. Res., Vol. 36, No. 7, 1997 2845

Figure 4. Experimental and correlated isothermal VLE for R50(1)-R1150(2) at 223.15 K using the SRK equation.

Figure 2. Liquid compressibility on the second iteration for isothermal VLE for R23(1)-R22(2) at 343.15 K using the SRK equation.

Figure 3. Vapor compressibility on the second iteration for isothermal VLE for R23(1)-R22(2) at 343.15 K using the SRK equation.

smallest real part as the liquid root and the compressibility with the largest real part as the vapor root. The comparison with Zc ) 1/3 made here is simply intended to show that the rules used for selecting liquid and vapor compressibility roots give physically consistent results. 4.2. Retrograde VLE. Difficulties in either bubbleor dew-point calculations often occur because of the presence of a retrograde loop because vapor and liquid compressibilities rapidly approach each other in this region of the phase diagram. Consequently the trivial solution attracts more and more “good” initial values, and all numerical methods become sensitive to initial conditions because of rapidly changing physical properties. It is also very difficult to actually close the retrograde loop in the region between the points of maximum pressure and maximum composition (see, for example, Cui and Donohue, 1994). Figure 4 shows the experimental data of Sagara et al. (1972) and the complete phase diagram for R50(1)R1150(2) at 223.15 K generated using the SRK equation with k12 ) 0 and polynomial mixing rules. The flash

model equations were solved using the complex trust region algorithm, and the complex solutions that lie beneath the single phase behavior in the critical region in the upper right hand corner of the diagram are also shown in that figure. Note that we could close the loop! In contrast, Newton’s method in the real or complex domains experienced many difficulties in the retrograde region. To illustrate the usefulness of the eigenvalue-eigenvector decomposition of Sridhar and Lucia (1995) in improving reliability, consider a typical computation in the retrograde region using 1 mol of feed with a composition of z ) (z1,z2) ) (0.72,0.28) and flash vapor flow and temperature specifications of V ) 0 and T ) 223.15 K, respectively. Note that this feed composition is on the bubble-point (upper) curve near the point of maximum pressure and has a corresponding equilibrium vapor composition clearly in the retrograde loop. Starting from the real initial value Z0 ) (x1,x2,y1,y2,L,P) ) (0.1, 0.9, 0.2, 0.8, 1.0, 50.95), complex domain Newton’s method failed to converge (or diverge) after 200 iterations. The complex domain trust region, on the other hand, found the singular point given by Z18 ) (0.719996, 0.279996, 0.734338, 0.265672, 1.0, 59.0479 - 0.246036i) after 18 iterations, then a second singular point at Z22 ) (0.720915, 0.275076, 0.731994, 0.271545, 1.0, 57.7447 - 0.054522i) on iteration 22, followed by convergence to the real-valued solution Z* ) Z34 ) (0.72, 0.28, 0.76857, 0.23143, 1.0, 59.4203) on iteration 34. The principal eigenvectors of negative curvature and corresponding eigenvalues associated with each singular point are (el, λ)18 ) (0.0017 + 0.08656i, 0.00273 + 0.17410i, -0.00300 - 0.39576i, -0.02233 - 0.14476i, 0.16033 + 0.00446i, 0.54962 + 0.00557i, -0.00742) and (el, λ)22 ) (-0.00124 + 0.20813i, 0.00006 - 0.38620i, 0.00182 - 0.47461i, 0.00309 + 0.00002i, 0.23109 + 0.00046i, -0.39362 + 0.00144i, -1.06533). Note also that the liquid and vapor compositions for each singular point are nearly equal, indicating that the trivial solution is attracting the iterates in this region, but that the eigenvalue-eigenvector decomposition allows the iterates to re-enter the complex domain and eventually converge to the real-valued solution. Similar behavior was observed for many other feed compositions in this area of the phase diagram. For feed compositions greater than the maximum composition, single phase behavior exists. Thus, numerical difficulties can be expected with real domain algorithms, and normally Newton and other algorithms either converge to a trivial solution or simply diverge

2846 Ind. Eng. Chem. Res., Vol. 36, No. 7, 1997

Figure 5. Experimental and correlated VLE for R744(1)-R170(2) at 223.15 K using the SRK equation.

Figure 6. Experimental and correlated VLE for R744(1)-R170(2) at 293.15 K using the SRK equation.

in this region of the phase diagram. The complex domain trust region method, however, is capable of finding complex-valued two-phase solutions in the single phase region and therefore exhibits finite step termination and provides useful information regarding phase existence. As an example, consider flash specifications of V ) 1 and T ) 223.15 K and 1 mol of feed with composition of z ) (z1,z2) ) (0.78, 0.22). Starting from an initial value of Z0 ) (x1,x2,y1,y2,L,P) ) (0.68, 0.32, 0.78, 0.22, 1.0, 50.0), the complex domain trust region method found a succession of three singular points before it converged to the complex-valued solution, Z* ) (0.68627 - 0.056993i, 0.31373 + 0.056993i, 0.78 + 0i, 0.22 + 0i, 1.0, 57.9559 - 3.52565i), in 27 iterations. The three singular points found along the way were Z5 ) (x1,x2,y1,y2,L,P) ) (0.701968 - 0.00817i, 0.298032 + 0.00817i, 0.780617 + 0.000721i, 0.219383 - 0.000721i, 0.0, 53.9422 + 0.39193i), Z12 ) (0.657552 + 0.000823i, 0.342448 + 0.000823i, 0.778998 + 0.0017607i, 0.221002 - 0.0017607i, 0.0, 55.2027 - 1.565261i), and Z19 ) (0.683586 - 0.028753i, 0.316413 + 0.028753i, 0.779861 - 0.0003692i, 0.220139 - 0.0003692i, 0.0, 57.4291 1.29008i). Complex domain Newton’s method also found the solution in only nine iterations, while real domain Newton’s method found a trivial solution. We note that similar behavior was observed for all complex solutions shown in Figure 4 and that, in our opinion, the complex solution gives a clear indication that the specifications are in the single phase region. 4.3. Azeotropic Mixtures. Another stringent modeling and equation-solving test is provided by azeotropic mixtures. Figures 5 and 6 illustrate that the SRK equation with k12 ) 0.14 gives good quantitative representation of the VLE data for carbon dioxide (1)ethane (2) (R744(1)-R170(2)) at both 223.15 and 293.15 K (Fredenslund and Mollerup, 1974). Note that there is an azeotrope at 223.15 K at a carbon dioxide composition of 0.6190 and that it is predicted reasonably well by the SRK equation with polynomial mixing rules. At 293.15 K, the azeotrope has disappeared and the phase diagram is represented by a two-lobed figure. Each lobe comes to a cusp in the middle portion of the composition range and single phase behavior exists between the lobes. Moreover, the lobe on the right appears as a line for the scale used in Figure 6 because the calculated liquid and vapor compositions differ only in the third significant digit. Note again, however, that our complex domain equation-solver was able to close both lobes.

Numerical difficulties with Newton-like methods for this system tend to occur at the higher temperature near the cusps of the lobes. As an illustration, consider 1 mol of feed with composition z ) (z1,z2) ) (0.388, 0.612) and flash specifications of V ) 0 and T ) 293.15 K. Starting from the initial values given by Z0 ) (x1,x2,y1,y2,L,P) ) (0.30, 0.70, 0.33, 0.67, 1.0, 55.5), Newton’s method in the real domain finds the trivial solution while Newton’s method in the complex domain converges to a periodic orbit. In contrast, the complex domain trust region method finds the singular point, Z13 ) (x1,x2,y1,y2,L,P) ) (0.388, 0.612, 0.389588 + 0.000208i, 0.610412 - 0.000208i, 1.0, 55.3747 + 0.004759i), on iteration 13, then the singular point, Z24 ) (x1,x2,y1,y2,L,P) ) (0.388, 0.612, 0.393177 + 0.000550i, 0.606823 - 0.000550i, 1.0, 55.3836 + 0.062576i), on iteration 24, and then finally the real-valued nontrivial solution, Z* ) (x1,x2,y1,y2,L,P) ) (0.388, 0.612, 0.391495, 0.608505, 1.0, 55.3785), on iteration 40. Furthermore, the liquid and vapor compressibilities at the solution are 0.332 692 and 0.350 904 and give an indication of the root finding challenge at the EOS level of the calculations. In our opinion, these numerical results clearly show that the complex domain EOS root finding and root assignment rules proposed in this work, as well as the eigenvalue-eigenvector decomposition of Sridhar and Lucia (1995), provide significant improvements in simulation reliability capabilities. Moreover, while it is true that complex domain calculations are roughly 1.5 to 3 times as expensive as their real counterparts (see, for example, Taylor et al., 1996), it has been our experience that these real domain methods (e.g., real domain Newton’s method, real domain dogleg strategies, real domain continuation methods, etc.) often fail to converge in regions of the phase diagram that contain strong nonlinearities and/or singularities such as in the retrograde region. 5. Conclusions Several binary mixtures of refrigerants were shown to be modeled adequately by the SRK equation with conventional polynomial mixing rules due to the spherical nature of the molecules, even at high pressure. Isothermal, fixed vapor flow (VT) flash processes involving refrigerant binaries modeled by a single equation of state were simulated using the complex domain trust region methods of Lucia and co-workers. Numerical

Ind. Eng. Chem. Res., Vol. 36, No. 7, 1997 2847

results show that calculated results are in good agreement with experimental data. Computational difficulties associated with EOS root finding and assignment were studied, and new rules for assigning compressibility roots to phases in situations where all roots are complex-valued were proposed and shown to give correct decisions, even in retrograde regions of the phase diagram. Other commonly encountered computational problems associated with singularity (i.e., azeotropy, trivial solutions, and retrograde behavior) were also studied. It was shown that the eigenvalue-eigenvector decomposition strategy of Sridhar and Lucia (1995) provides a reliable means of moving from singular points to solutions and that complex domain trust region methods are capable of finding complex-valued twophase flash solutions in the single phase region of the phase diagram and providing a clear indication of phase existence. Acknowledgment The work of X.G., D.L., and A.L. was supported by the National Science Foundation under Grant No. CTS9409158. Nomenclature ai(T) ) component i equation of state energy-volume parameter (cm6 bar mol-2) ai(Tci) ) component i equation of state energy-volume parameter evaluated at component i critical temperature, Tci (cm6 bar mol-2) aij(T) ) i-j pair equation of state energy-volume parameter (cm6 bar mol-2) aM(T) ) liquid or vapor mixture equation of state energyvolume parameter (cm6 bar mol-2) bi ) component i equation of state covolume parameter (cm3 mol-1) bM ) liquid or vapor mixture equation of state covolume parameter (cm3 mol-1) e ) eigenvector of the Hessian matrix of the two-norm of F fi ) component i feed molar flow rate (mol s-1) F ) total feed molar flow rate (mol s-1) F ) vector function of model equations defined on the complex domain G ) nonlinear iterative map defined on the complex domain GN ) nonlinear iterative map associated with Newton’s method J ) Jacobian matrix of F kij ) combining rule i-j binary interaction parameter Ki ) component i equilibrium K-value li ) component i liquid molar flow rate from flash (mol s-1) l ) vector of liquid component molar flow rates from flash (mol s-1) L ) total liquid molar flow rate from flash (mol s-1) m ) length of vector p n ) length of vector Z nc ) number of components in the mixture p ) complex-valued parameter vector P ) pressure (bar) Pci ) component i critical pressure (bar) R ) gas constant (cm3 bar mol-1 K-1) T ) absolute temperature (K) Tci ) component i critical temperature (K) vi ) component i vapor molar flow rate from flash (mol s-1) v ) vector of vapor component molar flow rates from flash (mol s-1) V ) total overhead vapor molar flow rate from flash (mol s-1)

Vci ) component i critical molar volume (cm3 mol-1) VM ) liquid or vapor mixture molar volume (cm3 mol-1) xi ) component i liquid mole fraction x ) vector of mole fractions in liquid mixture x1 ) vector of mole fractions in liquid phase 1 x2 ) vector of mole fraction in liquid phase 2 yi ) component i vapor mole fraction y ) vector of mole fractions in vapor mixture zi ) component i liquid or vapor mole fraction (xi or yi) Z ) complex-valued vector of unknowns Zk ) Z evaluated at iteration k Z* ) solution ZM ) liquid or vapor mixture compressibility factor Greek Letters φi ) component i fugacity coefficient in the liquid or vapor mixture φiL ) component i fugacity coefficient in the liquid mixture φiV ) component i fugacity coefficient in the vapor mixture ωi ) component i acentric factor λ ) eigenvalue of the Hessian matrix of the two-norm of F Subscripts c ) component ci ) component i critical value i ) component i j ) component j ij ) i-j pair M ) liquid or vapor mixture N ) Newton’s method 12 ) 1-2 pair λ ) associated with the eigenvalue λ Superscripts k ) iteration counter L ) liquid phase value V ) vapor phase value * ) value at solution 1 ) liquid phase 1 2 ) liquid phase 2

Literature Cited Cayley, A. The Newton-Fourier Imaginary Problem. Am. J. Math. 1879, II, 97. Cui, Y.; Donohue, M. D. Hydrogen Bonding Interactions in Polymer Solutions. AIChE Annual Meeting, San Francisco, November 1994, Paper 195f. Curry, J. H.; Garnett, L.; Sullivan, D. On the Iteration of a Rational Function: Computer Experiments with Newton’s Method. Commun. Math. Phys. 1983, 91, 267. Dahl, S.; Michelsen, M. L. High-Pressure Vapor-Liquid Equilibrium with a UNIFAC-Based Equation of State. AIChE J. 1990, 36, 1829. Eubank, P. T.; Elhassen, A. E.; Barrufet, M. A.; Whiting, W. B. Area Method for Prediction of Fluid Phase Equilibria. Ind. Eng. Chem. Res. 1992, 31, 942. Fatou, P. Sur les Equations Fonctionelle. Bull. Soc. Math Fr. 1919, 47, 161. Fredenslund, Aa.; Mollerup, J. Measurement and Prediction of Equilibrium Ratios for Ethane + Carbon Dioxide System. J. Chem. Soc., Faraday Trans. 1 1974, 70, 1653. Heidemann, R. A.; Kokal, S. L. Combined Excess Free Energy Models and Equations of State. Fluid Phase Equilib. 1990, 56, 17. Heidemann, R. A.; Abdel-Ghanni, R. An Approach to Multiphase Equilibrium Calculations. AIChE Annual Meeting, San Francisco, November 1994, Paper 117e. Huang, H.; Sandler, S. I. Prediction of Vapor-Liquid Equilibria at High Pressures Using Activity Coefficient Parameters Obtained from Low Pressure Data: A Comparison of Two Equation of State Mixing Rules. Ind. Eng. Chem. Res. 1993, 32, 1498.

2848 Ind. Eng. Chem. Res., Vol. 36, No. 7, 1997 Huron, M. J.; Vidal, J. New Mixing Rules in Simple Equations of State for Representing Vapour-Liquid Equilibria of Strongly Non-Ideal Liquid Mixtures. Fluid Phase Equilib. 1979, 3, 255. Julia, G. Memoir sur l’Iteration des Function Rationelles. J. Math Pures Appl. 1918, 4, 47. Lucia, A.; Guo, X.; Richey, P. J.; Derebail, R. Simple Process Equations, Fixed-Point Methods, and Chaos. AIChE J. 1990, 36, 641. Lucia, A.; Taylor, R. Complex Iterative Solutions to Process Model Equations? Comput. Chem. Eng. 1992, 16S, 387. Lucia, A.; Xu, J. Global Minima in Root Finding. In Recent Advances in Global Optimization; Floudas, C. A., Pardalos, P. M., Eds.; Princeton University Press: Princeton, NJ, 1992. Lucia, A.; Guo, X.; Wang, X. Process Simulation in the Complex Domain. AIChE J. 1993, 39, 461. Lucia, A.; Wang, X. Complex Domain Process Dynamics. Ind. Eng. Chem. Res. 1995, 34, 202. McDonald, C. M.; Floudas, C. A. Global Optimization for the Phase Stability Problem. AIChE J. 1995, 41, 1798. Michelsen, M. L. The Isothermal Flash ProblemsPart I: Stability. Fluid Phase Equilib. 1982a, 9, 1. Michelsen, M. L. The Isothermal Flash ProblemsPart II: Phase Split Calculation. Fluid Phase Equilib. 1982b, 9, 21. Peng, C.-L.; Stein, F. P.; Gow, A. S. An Enthalpy-Based Cubic Equation of State Mixing Rule for Cross-Prediction of Excess Thermodynamic Properties of Hydrocarbon and Halogenated Refrigerant Mixtures. Fluid Phase Equilib. 1995, 108, 79. Peng, D.-Y.; Robinson, D. B. A New Two Constant Equation of State. Ind. Eng. Chem. Fundam. 1976, 15, 59. Powell, M. J. D. A Hybrid Method for Nonlinear Equations. In Numerical Methods for Nonlinear Algebraic Equations; Rabinowitz, P., Ed.; Gordon and Breech: New York, 1970. Roth, H.; Peters-Gerth, P.; Lucas, K. Experimental Vapor-Liquid Equilibria in the Systems R22-CO2, CS2-R22, R23-CO2, CS2R23 and Their Correlation by Equations of State. Fluid Phase Equilib. 1992, 73, 147.

Sagara, H.; Arai, Y.; Saito, S. Vapor-Liquid Equilibria of Binary and Ternary Systems Containing Hydrogen and Light Hydrocarbons. J. Chem. Eng. Jpn. 1972, 5, 339. Schnepper, C.; Stadtherr, M. A.; Brennecke, J. F. Robust Calculation of High Pressure Phase Equilibria Using Interval Methods. AIChE Annual Meeting, San Francisco, 1994, Paper 117f. Soave, G. Equilibrium Constants from a Modified Redlich-Kwong Equation of State. Chem. Eng. Sci. 1972, 27, 1197. Sridhar, L. N.; Lucia, A. Process Analysis in the Complex Domain. AIChE J. 1995, 41, 585. Stein, F. P.; Proust, P. C. Vapor-Liquid Equilibria of the Trifluoromethane-Trifluorochloromethane System. J. Chem. Eng. Data 1971, 16, 389. Sun, A. C.; Seider, W. D. Homotopy Continuation Algorithm for Global Optimization. In Recent Advances in Global Optimization; Floudas, C. A., Pardalos, P. M., Eds.; Princeton University Press: Princeton, NJ, 1992. Taylor, R.; Achuthan, K.; Lucia, A. Complex Domain Distillation Calculations. Comput. Chem. Eng. 1996, 20, 93. Trangenstein, J. A. Customized Minimization Techniques for Phase Equilibrium Computations in Reservoir Simulation. Chem. Eng. Sci. 1987, 42, 2847. Whitson, C. H.; Michelsen, M. L. The Negative Flash. Fluid Phase Equilib. 1989, 53, 51. Wong, D. S. H.; Sandler, S. I. A Theoretically Correct Mixing Rule for Cubic Equations of State. AIChE J. 1992, 38, 671.

Received for review January 15, 1997 Revised manuscript received April 16, 1997 Accepted April 17, 1997X IE970043X

X Abstract published in Advance ACS Abstracts, June 1, 1997.