Ind. Eng. Chem. Res. 2000, 39, 4415-4421
4415
Phase Equilibrium Calculations for Semicontinuous Mixtures Subject to Gravitational Fields Roge´ rio O. Espo´ sito, Marcelo Castier,* and Frederico W. Tavares Escola de Quı´mica, Universidade Federal do Rio de Janeiro, C.P. 68542, Rio de Janeiro-RJ 21949-900, Brazil
A recently proposed algorithm for performing thermodynamic equilibrium calculations subject to the influence of the gravitational field is extended to semicontinuous mixtures. The characterization of the continuous families is done by a pseudocomponent discretization based on a convenient choice of orthogonal polynomials to sample the distribution function interval. Then, adequate quadrature points and weights are calculated. We present results for a family of hydrocarbons from C8 to C20 with 51 mol % of methane for one-phase, two-phase, three-phase, and near-critical systems. Our results show the validity of the procedure for simulations of deep reservoirs and high-pressure storage conditions, solving problems with a high number of independent variables. Introduction The characterization of oil reservoir fluids for phase equilibrium calculations is difficult because of the very large number of components, and the approach known as semicontinuous mixture thermodynamics1-5 has been increasingly used to solve these problems. In this approach, it is usual to characterize the lighter components individually and to group the heavier fractions into one or more continuous families. Much of the work done in this area prior to 1990 was reviewed by Cotterman and Prausnitz6 and by Rochocz.7 More recent work was reviewed by Lira-Galeana et al.8 and Rochocz et al.9 Most of the papers on algorithms for phase equilibrium calculations in semicontinuous mixtures are for saturation (bubble or dew) points or flashes at specified temperature and pressure. Relatively few papers focus on an important aspect for the oil industry, namely, the distribution of species in gravitational fields. A recent example is the paper of Lira-Galeana et al.,8 who used the thermodynamics of continuous mixtures to compute oil reservoir pressures and compositions at different levels once the composition at a given level was specified. The objective of this paper is to present a procedure for calculating the segregation of semicontinuous mixtures subject to a gravitational field, at specified temperature, density, and global mole fractions. Our goal is to determine the distribution of the species and the heights of the interfaces between the phases present in the system, which is important information in deep oil fields. Another application is in the determination of storage conditions in oil tanks. The procedure is based on a recently developed algorithm for the segregation of fluids in gravitational fields,10 which is based on the minimization of a modified Helmholtz free energy that includes a term to account for the effect of the gravitational field. Modeling of Semicontinuous Mixtures The procedure presented in this paper is not limited to any particular distribution function or to a single * Corresponding author. E-mail:
[email protected]. Fax: +55-21-542-6376. Phone: +55-21-562-7607.
continuous family. However, for simplicity, we will only consider one continuous family in each example. To characterize the continuous family, we will follow Shibata et al.11 and Rochocz et al.,9 choosing a bounded exponential function in which the single characterization variable I is the number of carbon atoms, assumed to vary continuously between two finite integration limits. For a continuous family k, the expression for this distribution function is
Fk(Ik) ) Ck exp(-DkIk)
(1)
The expression is normalized so that
∫ηφ Fk(Ik) dIk ) 1 k
k
(2)
Therefore, the mole fractions satisfy the following equation: nd
∑ i)1
nfam
xi +
∑ k)1
nd
xk )
∑ i)1
nfam
xi +
∑ xk ∫η k)1
φk k
Fk(Ik) dIk ) 1 (3)
where nd is the number of discrete componens, nfam is the number of continuous families, xk is the global mole fraction of family k, Ik is the characterization index of family k, and ηk and φk are the lower and upper bounds, respectively, of the kth family’s distribution function. Dropping the family index k, using the correlation of Katz and Firoozabadi12 for the average number of atoms of carbon in the molecules of the continuous family, CN, as function of the average molecular weight of the continuous fraction, MW+
CN )
MW+ + 4 14
(4)
and defining
∆ ) D(φ - η) it can be shown9,11 that
10.1021/ie000268z CCC: $19.00 © 2000 American Chemical Society Published on Web 09/19/2000
(5)
4416
Ind. Eng. Chem. Res., Vol. 39, No. 11, 2000
C)
(
(e
)
-Dη
D - e-Dφ)
CN - η 1 e-∆ ) φ-η ∆ (1 - e-∆)
(6)
energy, a modified Helmholtz (Aj*) free energy is obtained nc
(7)
where
Aj* ) -PjVj +
nijµij* ∑ i)1
(11)
where µij* is a modified chemical potential, given by
η ) Cm - 0.5
(8)
φ ) CM + 0.5
(9)
where Cm and CM are the minimum and maximum number, respectively, of carbon atoms in the molecules of the continuous family. Given Cm, CM, and MW+, eqs 8 and 9 are used to compute η and φ. Equation 4 is used to calculate CN, the nonlinear eq 7 is solved to determine ∆, and D is calculated in eq 5. Because of the complexity of the thermodynamic models used, it is often not possible to obtain analytical expressions for the integrals that appear in the semicontinuous formulation. We used the generalized GaussLaguerre method proposed by Shibata et al.,11 as implemented by Rochocz et al.9 The procedure for obtaining quadrature points and weights is based on the Wheeler algorithm,13 which allows the calculation of the two coefficients in the recurrence law of the orthogonal polynomials introduced by Sack and Donovan.14 These coefficients are computed in the routine orthog, a Fortran code presented by Press et al.15 Given the recurrence law, we used the Golub-Welsh16 algorithm (also implemented by Press et al.15 in their routine gaucof) to compute the quadrature points and weights that provide us the pseudocomponents that represent the continuous family, which are then treated in the rest of the simulation as discrete species. We used the Peng-Robinson17 equation of state with the van der Waals one-fluid mixing rules, i.e., quadratic and linear mixing rules for the a and b parameters, respectively. Pure-component properties for the discrete species were taken from Reid et al.18 For the critical properties and acentric factor of the pseudocomponents corresponding to the quadrature points, we used linear interpolations with the values obtained from the table of generalized single-carbon-number physical properties of Whitson.19 Because the emphasis of this paper is on the calculation procedure, we set the kij binary interaction parameters equal to zero, even though this choice may not yield reliable phase equilibrium predictions in systems containing compounds with large hydrocarbon chains. Formulation In the formulation of the gravitational segregation problem, we assume that the numerical treatment of the continuous fraction will require discretization. In this way, the problem is formulated as if the mixture contained only discrete components. We consider a system split in ny layers that are assumed to be separated by diathermal, rigid, and impermeable walls. Neglecting the effect of kinetic energy, the energy of a layer j (Ej) is
Ej ) Uj + Ep,j
(10)
where Uj and Ep,j are the internal and potential energies, respectively. From a Legendre transform of the
µij* ≡ µij + Wighj
(12)
The symbols Pj, Vj, nij and µij refer to properties at level j, being the pressure, volume, mole number, and chemical potential of each component, respectively. The symbols Wi and g denote the molar weight of the component i and the acceleration of gravity, respectively; hj is the height at the middle of layer j. For a system with a specified mole number of each species (N), total volume (V), temperature (T), and total height (h), thermodynamic equilibrium occurs at the minimum of A*, which is given by a sum over all ny discrete layers, i.e. ny
A* )
Aj* ∑ j)1
(13)
Solution Procedure The minimization procedure for the modified Helmholtz function is described in detail by Espo´sito et al.10 Here, we only emphasize the most important features. For convenience, we used A*/RT instead of A* as objective function. The minimization of A*/RT was carried out in two loops. In the outer loop, we used the Simplex method20 to iterate on the volume fraction of the phases. In the inner loop, we used the Murray21 method to find the minimum of A*/RT, given the configuration provided by the external loop. We used the independent yield factors θij as variables in the inner loop.22 They are defined as the ratio of the number of moles of component i in phase j, nij, to the total number of moles of component i, Ni, i.e.
θij )
nij Ni
(14)
The largest yield factor for a component was considered as a dependent variable. Therefore, there are [nc(ny - 1)] independent variables in the minimization. The first and second derivatives of A*/RT with respect to the independent variables are given by
∂(A*/RT) ) Ni(ln aij* - ln aiJ*) ∂θij
[
( ) ( )
∂ ln a*ij ∂2(A*/(RT)) ) NiNl (δjm - δjM) ∂θij∂θlm ∂nlj (δJm - δJM)
∂ ln a*iJ ∂nlJ
(15) -
T,Vj,nkj,k*1
T,Vj,nkJ,k*1
]
(16)
where the capital subscripts stand for the dependent layers (i.e., those having the largest yield factor for components i and l, respectively) and aij* is the gravitymodified activity of component i in the jth layer, given by
Ind. Eng. Chem. Res., Vol. 39, No. 11, 2000 4417
aij* ≡ aij exp
( ) Wighj RT
(17)
nqb
Initial Estimates for the Outer Loop. To generate good initial guesses for the volume fraction of the phases, we solved a series of phase equilibrium (PE) problems without the effect of the gravitational field to find the pressure that matches the specified global density within a desired precision. For a specified number of moles in the system, this problem can be formulated in terms of a volume-difference equation, which was solved by the Newton-Raphson method. The details of the solution procedure can be found in Espo´sito et al.10 Having matched the value of the specified global density, we obtained a good initial estimate for the volume fraction of each phase present in the system. Each phase was then split into layers that were supposed to have the same size. The total number of layers (ny) was therefore equal to the specified number of layers per phase times the number of phases (nf). Initial Estimates for the Inner Loop. To obtain initial estimates for the inner loop in the first iteration of the outer loop, we used the results from the PE calculation without gravity as follows. For the yield factors in a liquid phase, we equally distributed the amount of each component in that phase throughout its layers.To obtain an initial estimate for the yield factors in a vapor phase, we used the analytical solution for an ideal gas23
θij ) exp ∑ m
[
1
]
Wig(hj - hm) RT
Table 1. Carbon Number (Cnum) and Feed Mole Fraction (zi) for Pseudocomponentsa 2 Cnum
4 zi
Cnum
6 zi
Cnum
8 zi
Cnum
zi
pseudo 1 9.17 0.362 8.15 0.184 7.85 0.105 7.71 0.068 pseudo 2 15.90 0.128 10.84 0.200 9.29 0.162 8.61 0.121 pseudo 3 15.07 0.086 11.71 0.124 10.16 0.120 pseudo 4 19.2 0.020 14.74 0.065 12.21 0.087 pseudo 5 17.77 0.026 14.56 0.051 pseudo 6 19.93 0.008 16.92 0.027 pseudo 7 18.92 0.012 pseudo 8 20.18 0.004 a The global mole fraction of methane is 0.51. b n ) number of q quadrature of points.
(18)
where the index m runs over the layers of the vapor phase. After the first convergence of the inner loop, the volume fractions were updated in the outer loop. For the second and subsequent iterations of the outer loop, we used the converged values of the yield fractions from the previous iteration. Results and Discussion One-Phase Systems. We present results for a mixture composed of methane and a series of n-alkanes from C8 to C20 presented by Lira-Galeana et al.8 Our values of Cm and CM are then 8 and 20, respectively. We set MW+ equal to 149 kg/kgmol, according to Figure 1 of their paper. From specified data at a reference level, they computed composition and pressure at any other level of the fluid column, following a procedure similar to a bubble-point calculation. Here, we deal with different specifications. To provide a basis for comparison, our specified global density (569.5 kg/m3) was set in such a way that the pressure profile obtained in our simulation using eight quadrature points is very similar to that of Lira-Galeana et al.8 The global mole fraction of methane was set equal to 0.51. For the pseudocomponents, the weights of generalized Gauss-Laguerre quadrature were multiplied by the global mole fraction of the continuous family (equal to 0.49) and are shown in Table 1, together with their calculated carbon numbers. For convenience, the pseudocomponent critical properties and acentric factors (Table 2) were obtained by interpolation in the generalized
Figure 1. Pressure profile for the one-phase system.
single-carbon-number tables of Whitson,19 even though the original example of Lira-Galeana et al.8 referred to a series of n-alkanes. The temperature of the 500-mdeep reservoir was set to 394.4 K.8 As we set the number of layers per phase equal to 10, the phase equilibrium problem has (10 - 1) × 9 ) 81 variables when we deal with 8 pseudocomponents. Figures 1 and 2 show pressure and methane mole fraction profiles, respectively, obtained in runs using 2, 4, 6 and 8 quadrature points. One can note that the profiles obtained with 6, 4, and 2 quadrature points increasingly deviate from the 8-quadrature-point solution (used to match the reported data) for the same global density (569.5 kg/m3). Figures 3 and 4 show the mole fraction profiles of the pseudocomponents for 2 and 4 quadrature points, respectively. Results for 6 and 8 points are similar but were omitted for visual convenience. In these figures, the dotted lines represent the global mole fractions, which are included to help the visualization of the gravity action on the pseudocomponents, which tend to segregate to the bottom of the reservoir. Two-Phase Systems. In this example, we considered the same mixture of the one-phase reservoir (the continuous family being represented by 8 pseudocomponents), but we set global density equal to 500 kg/m3. Again, the reservoir was 500-m-deep, and the temperature was set equal to 394.4 K. Under these conditions, our calculations indicate the existence of a liquid and a vapor phase in the system. Their interface is at a depth of approximately 100 m, as one can see by the change
4418
Ind. Eng. Chem. Res., Vol. 39, No. 11, 2000
Table 2. Critical Properties and Acentric Factors for Pseudocomponents nqa
pseudo 1 pseudo 2 pseudo 3 pseudo 4 pseudo 5 pseudo 6 pseudo 7 pseudo 8 a
2
4
Tc (K)
Pc (106 Pa)
ω
606.9 738.4
2.62 1.69
0.354 0.579
6
Tc (K)
Pc (106 Pa)
ω
579.3 644.4 725 780.4
2.88 2.29 1.77 1.49
0.318 0.413 0.552 0.668
8
Tc (K)
Pc (106 Pa)
ω
570.9 609.7 662.2 719.3 764.2 789.2
2.96 2.60 2.15 1.81 1.56 1.44
0.307 0.359 0.444 0.541 0.632 0.688
Tc (K)
Pc (106 Pa)
ω
567.3 592.2 629.5 672.1 716.2 753.8 777.1 792.0
2.98 2.77 2.42 2.08 1.82 1.62 1.50 1.43
0.303 0.334 0.390 0.460 0.535 0.610 0.660 0.695
nq ) number of quadrature of points.
Figure 2. Methane mole fraction profile for the one-phase system.
Figure 3. Pseudocomponent mole fraction profiles for the onephase system (nq ) 2). Vertical dotted lines represent the global mole fractions for visual reference.
of slope in pressure profile in Figure 5. Figure 6 shows mole fraction profiles for methane and some pseudocomponents. One can see a methane-rich gas phase, as well as the gravity action mainly on the slope of this component profile. Once this is a two-phase system and we again set the number of layers per phase equal to
Figure 4. Pseudocomponent mole fraction profiles for the onephase system (nq ) 4). Vertical dotted lines represent the global mole fractions for visual reference.
Figure 5. Pressure profile for the two-phase system.
10, the algorithm dealt with (20 - 1) × 9 ) 171 variables in the calculation. Three-Phase Systems. In this example, we calculated the mole fraction and pressure profiles with the same global mole fraction ratio between methane and the continuous family (8 quadrature points) as in the first example. However, we added a considerable amount
Ind. Eng. Chem. Res., Vol. 39, No. 11, 2000 4419
Figure 6. Mole fraction profile for two-phase system.
Figure 7. Pressure profile for the three-phase system.
Table 3. Pseudocomponent Feed Composition for the Multiphase Reservoir feed mole fraction pseudo 1 pseudo 2 pseudo 3 pseudo 4 pseudo 5 pseudo 6 pseudo 7 pseudo 8 methane water
0.040 0.072 0.071 0.051 0.030 0.016 0.007 0.002 0.300 0.410
of water to allow the presence of three phases in the system. All global mole fractions are shown in Table 3. We set the temperature to 300 K. The specified global density was 400 kg/m3. The height of the reservoir was arbitrarily set equal to 10 m. We again chose 10 layers per phase as input for the simulation, and therefore, ny ) 30. The algorithm dealt with (30 - 1) × 10 ) 290 independent yield factors in this calculation. We found the water-oil and oil-gas interfaces at 0.625 and 5.360 m above the bottom of the system, respectively. In Figure 7, these interfaces correspond to the change of slope in the pressure profile. The densest phase is almost hydrocarbon-free (Figure 8). Figure 9 shows the mole fraction profiles for some of the pseudocomponents. These results show that the procedure is able to simulate strongly nonideal mixtures in storage tanks under high pressures. Systems Near the Critical Point. For a system near the critical point, some thermodynamic properties, such as the isothermal compressibility, tend to diverge. For such systems, small pressure variations in a fluid column even in short or moderate-height cells can lead to considerable variations in density and mole fractions profiles.24 To specify this example, we first used a modified9 version of the Hicks and Young25 algorithm and computed the critical conditions (Tc, Pc, and vc) for the n-alkanes-methane mixture (nq ) 8) presented in Table 1. We found Tc ) 615.11 K, Pc ) 106.98 × 105 Pa, and vc ) 342.87 × 10-3 m3/kgmol. The critical density (Fc) can then be computed as Fc ) (MW+ × 0.49 + 16 × 0.51)/ vc ) 236.7 kg/m3, where MW+ ) 149 and the numerator
Figure 8. Methane and water mole fraction profiles for the threephase system.
of the latter expression accounts for the mean molecular weight of the mixture. We made our segregation calculations with T ) Tc and global density set equal to 237.7 kg/m3. We also specified 30 layers per phase in a 0.30m-high cell, but only one phase was predicted for this system. Therefore, the total number of calculated independent yield factors was equal to (30 - 1) × 9 ) 261. The density profile is shown in Figure 10. One can see a considerable variation in this quantity, even for this small-height system. The relative mole fraction deviation profiles of methane and two pseudocomponents are shown in Figure 11. One can see a variation of approximately 5% in relation to the methane global composition throughout the cell. Also, one can see that the pseudocomponents tend to segregate to the bottom and that this effect is more pronounced for pseudocomponent 7 (Cnum ) 18.92) than for pseudocomponent 5 (Cnum ) 14.56). Although we have not found experimental data in the literature for a quantitative comparison, our results are in qualitative agreement with those of Chang et al.,26 who predicted the properties of binary
4420
Ind. Eng. Chem. Res., Vol. 39, No. 11, 2000
Figure 11. Relative mole fraction deviations near the critical point [(xi - zi)/zi]. Figure 9. Pseudocomponent mole fractions for the three-phase system.
the continuous families. In our examples, we successfully used a bounded exponential function that was discretized using a generalized Gauss-Laguerre quadrature. Our first example provided a semiquantitative comparison with data reported by Lira-Galeana et al.,8 as some of the specifications are not the same. By specifying a lower global density, one can obtain two phases for the same global composition. This was done in our second example. Water was added in the third example to obtain three phases. Finally, we computed the critical properties of the proposed mixture and obtained good results for a small-height cell. All examples show the validity of the procedure for simulating deep reservoirs and high-pressure storage conditions, solving problems with a high number of independent variables. Acknowledgment
Figure 10. Density profile for the semicontinuous mixture near its critical point.
mixtures near critical points using the Leung-Griffiths27 equation of state.
M.C. acknowledges the hospitality of Prof. Jose´ Tojo Sua´rez, from the University of Vigo, Spain, where part of this work was done. The authors thank Charlles R.A. Abreu for helping in the preparation of the figures. Financial support was provided by CNPq/Brazil, PRONEX (Grant 124/96), and ANP (Ageˆncia Nacional do Petro´leo). Symbol List
Conclusions The phase equilibrium problem in a system subject to the influence of a gravitational field, formulated and solved as a multiphase flash proposed by Espo´sito et al.,10 was extended to semicontinuous mixtures. The specified system variables were the temperature, volume, and global mole numbers. The phase equilibrium problem was solved by minimization of a modified form of the Helmholtz function that includes a term to account for the effect of the gravitational field. The minimization was carried out in two nested loops. In the outer loop, the phase volume fractions were updated. In the inner loop, given a set of phase volume fractions, the distribution of species among the layers was updated. The developed procedure is not limited to any thermodynamic model or to any characterization scheme of
aij ) traditional activity of component i in layer j aij* ) modified activity of component i in layer j A ) traditional Helmholtz free energy A* ) modified Helmholtz free energy CN ) average carbon number of the continuous family Cnum ) pseudocomponent carbon number Cm ) minimum carbon number of the continuous family CM ) maximum carbon number of the continuous family Ep ) potential energy g ) gravity acceleration hj ) height at the center of layer j nc ) number of components nd ) number of discrete components nfam ) number of continuous families nf ) number of phases present in the system nij ) number of moles of species i in layer j Ni ) total number of moles of species i in the system nq ) number of quadrature points
Ind. Eng. Chem. Res., Vol. 39, No. 11, 2000 4421 ny ) total number of layers Pj ) pressure of layer j R ) universal gas constant T ) temperature Tc ) critical temperature Uj ) internal energy of layer j vc ) critical molar volume Vj ) volume of layer j xi ) local mole fraction of component i Wi ) molecular weight of component i zi ) global mole fraction of component i Greek Letters δjm ) Kronecker delta function θij ) yield factor of component i in layer j (θij ) nij/Ni) µij ) traditional chemical potential of component i in layer j µij* ) modified chemical potential of component i in layer j η ) lower bound of family distribution function φ ) upper bound of family distribution function Fc ) critical density
Literature Cited (1) Bownman, J. R. Distilation of an Indefinite Number of Components. Ind. Eng. Chem. 1949, 41, 2004. (2) Edmister, W. C.; Buchanan, D. H. Applications of Integral Distillation Calculation Methodology. Chem. Eng. Prog. Symp. Ser. 1953, 6, 69. (3) Ra¨tzch, M. T.; Kehlen, H. Continuous Thermodynamics of Complex Mixtures. Fluid Phase Equilib. 1983, 14, 225. (4) Kehlen, H.; Ra¨tzsch, M. T. Separate Treatment of Paraffins and Aromatics in Complex Hydrocarbon Mixtures by Continuous Thermodynamics. Z. Phys. Chem. 1984, 265, 1049. (5) Kehlen, H.; Ra¨tzsch, M. T.; Bergmann, J. Continuous Thermodynamics of Multicomponent Systems. AIChE J. 1985, 31, 1136. (6) Cotterman, R. L.; Prausnitz, J. M. Application of Continuous Thermodynamics to Natural-Gas Mixtures. Rev.Inst. Fr. Pe´ t. 1990, 45, 633. (7) Rochocz, G. L. M.Sc. Thesis, COPPE, Federal University of Rio de Janeiro, Brazil, 1990. (8) Lira-Galeana, C.; Firoozabadi, A.; Prausnitz, J. M. Computation of Compositional Grading in Hydrocarbon Reservoirs. Application of Continuous Thermodynamics. Fluid Phase Equilib. 1994, 102, 143. (9) Rochocz, G. L.; Castier, M.; Sandler, S. I. Critical Point Calculations for Semi-Continuous Mixtures. Fluid Phase Equilib. 1997, 139, 137.
(10) Espo´sito, R. O.; Castier, M.; Tavares, F. W. Calculations of Thermodynamic Equilibrium in Systems Subject to Gravitational Fields. Chem. Eng. Sci. 2000, 55, 3495. (11) Shibata, S. K.; Sandler, S. I.; Behrens, R. A. Phase Equilibrium Calculations for Continuous and Semicontinuous mixtures. Chem. Eng. Sci. 1987, 42, 1977. (12) Katz, D. L.; Firoozabadi, A. Predicting Phase Behavior of Condensate Crude-Oil Systems Using Methane Interaction Coefficients. J. Pet. Technol. 1978, 20, 1649. (13) Wheeler, J. C. Modified Moments and Gaussian Quadratures. Rocky Mt. J. Math. 1974, 4, 287. (14) Sack, R. A.; Donovan, A. F. An Algorithm for Gaussian Quadrature Given Modified Moments. Numer. Math. 1972, 18, 465 (15) Press: W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes in Fortran, 2nd ed.; Cambridge University Press: Cambridge, U.K., 1992. (16) Golub, G. H.; Welsch, J. H. Calculations of Quadrature Rules. Math. Comput. 1969, 23, 221. (17) Peng, D. Y.; Robinson, D. B. A Simple Two-Constant Equation of State. Ind. Eng. Chem. Fundam. 1976, 15, 59. (18) Reid, R. C., Prausnitz, J. M.; Poling, B. E. The Properties of Gases and Liquids, 4th ed.; McGraw-Hill: New York, 1987. (19) Whitson, C. H. Characterizing Hydrocarbon Plus Fractions. Soc. Pet. Eng. J. 1983, August, 683. (20) Nelder, J. A.; Mead, R. A Simplex Method for Function Minimization. Comput. J. 1965, 7, 308. (21) Murray, W. Numerical methods for unconstrained optimization; Academic Press: New York, 1972 (22) Michelsen, M. L. The Isothermal Flash Problem. Part I. Stability. Fluid Phase Equilib. 1982, 9, 1. (23) Espo´sito, R. O. M.Sc. Thesis. School of Chemistry, Federal University of Rio de Janeiro, Brazil, 1999. (24) Sychev, V. V. Complex Thermodynamics Systems; Mir Publishers: Moscow, 1981. (25) Hicks, C. P.; Young, C. L. Theoretical Predictions of Phase Behavior at High Temperatures and Pressures for non-Polar Mixtures: 1. Computer Solution Techniques and Stability Tests. J. Chem. Soc., Faraday Trans. 2 1977, 73, 597. (26) Chang, R. F.; Levelt-Sengers, J. M. H.; Doiron, T.; Jones, J. Gravity-Induced Density and Concentration Profiles in Binary Mixtures Near Gas-Liquid Critical Lines. J. Chem. Phys. 1983, 79 (6), 3058. (27) Leung, S. S.; Griffiths, R. B. Thermodynamic Properties near the liquid-vapor critical line in mixtures of He3 and He.4 Phys. Rev. A 1973, 8, 2670.
Received for review February 22, 2000 Accepted July 22, 2000 IE000268Z
