Ind. Eng. Chem. Res. 2001, 40, 6213-6220
6213
Continuous-Phase Equilibrium Problem: Quadrature Compositional Characterization and Asymptotic Convergence Joseph A. Labadie and Kraemer D. Luks* Department of Chemical Engineering, University of Tulsa, 600 South College Avenue, Tulsa, Oklahoma 74104-3189
In this paper, the authors examine the fluid-phase and solid-fluid-phase equilibrium problems in order to determine whether there is an inherent resistance to the convergence of quadraturebased computations in general for phase equilibrium problems. Our findings suggest that convergence can be expected for fluid-phase equilibria in general. For the less familiar solidfluid-phase equilibrium problem, convergence can depend on the postulated description of the solid phase. Introduction Continuous systems are by definition composed of so many components that their specific identity defies precise analysis. For example, it can be useful to view the C7+ portion of a crude oil as a continuous mixture. A live oil would then be a semicontinuous mixture of well-defined gaseous species plus a distributed C7+ fraction which may exceed 60 mol % of the system; in turn, a stock tank oil might be viewed as a single continuous system. With respect to such continuous systems, it has been shown1-8 that quadrature can be a very effective way of specifying the pseudocomponents of the C7+ portion of an oil, provided one can accurately represent the distribution of the many species as a continuous function of some parameter(s), say, the carbon number. Additionally, quadrature discretization can be shown to avoid certain conservation of mass inconsistencies that are inherent in the continuous formalism of phase equilibrium computations, including saturation computations.9 Convergence usually occurs at a relatively low order/ level of quadrature in the fluid-phase equilibrium problem when using an equation of state (EOS) for example. Such asymptotic convergence formally ensures that one’s computed result will be equivalent to that which would be obtained if many (n f ∞) pseudocomponents were employed. In earlier work,9 we showed that pseudocomponent characterization of continuous (crude oil) mixtures using finite Laguerre-Gauss quadrature is able to generate dew-point calculations that are asymptotically the same (“converged”) to a reasonable degree of precision for n ) 4 or more pseudocomponents. (The system addressed was a live oil with a C7+ fraction of about 30 mol % having a molecular weight of 178. A variant of the Redlich-Kwong EOS was used to model the thermodynamics.) This dew-point result is appealing because for such a system the molecular nature of the heaviest C7+ pseudocomponent is strongly dependent on * To whom correspondence should be directed. Phone: (918) 631-2974. Fax: (918) 631-3268. E-mail: kraemer-luks@ utulsa.edu.
the value of the last discretized component’s carbon number Cn. It appears that the effect of an increase in the molecular weight of the heaviest C7+ pseudocomponent on the computed dew-point pressure is offset by the accompanying decrease in its mole fraction. Equivalently, one could argue that the mathematical complexity of that dew-point problem was essentially captured when n g 4. In a recent paper,10 the authors presented computational results for the solid-fluid-phase equilibrium problem for an n-paraffinic C7+ model fluid whose pseudocomponent compositional characterization is determined by finite Laguerre-Gauss quadrature, with the assumption that the initial solid phase formed at the system’s crystal point was a “pure” pseudocomponent. The crystal-point temperature results did not display asymptotic convergence as the quadrature level increased, an unanticipated result. It was surmised that the highly nonlinear nature of the fluid phase’s Gibbs energy of mixing, particularly that manifested in the ideal part of the Gibbs energy of mixing, presents a serious challenge to obtaining a satisfactory asymptotic result. Furthermore, it was suggested that the convergence of fluid-phase equilibrium results with increasing quadrature level could likewise experience hindered convergence, although such an effect anecdotally appears to be of considerably lesser extent. In this present study, we examine the well-known EOS-based fluid-phase equilibrium problem in general, and the solid-fluid-phase equilibrium problem in particular, addressing specifically what effects Gibbs energy nonlinearities may have on asymptotic convergence of the quadrature problem. Fluid-Phase Equilibrium Problem The two-phase fluid-phase equilibrium problem formally requires one to find the solution to the set of equations balancing chemical potentials between the two phases R and β for each component i:
µiR(P,T,{xi}) ) µiβ(P,T,{yi}), i ) 1, ..., n
(1)
An equivalent statement of these mass flow equilibria
10.1021/ie0103512 CCC: $20.00 © 2001 American Chemical Society Published on Web 12/01/2001
6214
Ind. Eng. Chem. Res., Vol. 40, No. 26, 2001
Figure 1. SRK EOS prediction of the reduced Gibbs energy of mixing versus mole fraction of n-hexane for the binary mixture propane + n-hexane at 420 K and 20 bar (kij ) 0.0). The common tangent locating the equilibrium phase compositions is shown as a dashed line. The thin solid curve is the reduced ideal Gibbs energy of mixing.
equates the well-known fugacity function:
ˆfiR(P,T,{xi}) ) ˆfiβ(P,T,{yi}), i ) 1, ..., n
(2)
Equations 1 and 2 are a set of n equations in 2n unknowns (P, T, {xi}, {yi}). Closure of the problem is traditionally achieved by either fixing n of the variables (as in a saturation calculation) or coupling eq 2 with n mass balance equations and fixing T and P (as in a flash calculation):
system of n components with equal mole fractions, the ideal Gibbs energy of mixing would be -RT ln(n), which f -∞ as n f ∞. In our earlier paper on solid-fluid-phase equilibria,10 we demonstrated computationally that, as the level of quadrature increased, the ideal Gibbs energy of mixing function for a fluid phase grew; in contrast, the excess Gibbs energy function converged to five significant figures at the quadrature level n ) 5. Thus, divergence of the total Gibbs energy of mixing function as n f large (which correspondingly increases the dimensionality of the Gibbs energy function space) is unavoidable and raises the question as to why the solution to eq 1 (eq 2) should converge in the context of the quadrature problem as rapidly as it does, because the ideal Gibbs energy of mixing appears to be an essential component of the graphical construction of the solution (illustrated in Figure 1). On the other hand, gas-liquid dew-point calculations9 suggest that convergence of the fluid-phase equilibrium problem may be (convergence-wise) independent of the ideal Gibbs energy of mixing part of the total Gibbs energy of mixing function. In other words, convergence of the fluid-phase equilibrium problem is tantamount to convergence of the excess Gibbs energy function. Making this connection can be shown to be appropriate and will be demonstrated here. Equation 4 can be rewritten for an n-component system as
[( ∑ [ ( n
∆G
mix
xi ln ∑ i)1
(P,T,{xi}) ) RT
xi ln
i)1
Lxi + Vyi ) zi, i ) 1, ..., n
(3)
Equation 1 can be rigorously viewed as the solution geometrically obtained by constructing a common tangent hyperplane to the Gibbs energy surface at some T ) T0 and P ) P0, with the tangent points locating the phase compositions {xi} and {yi} consistent with a minimum Gibbs energy for the two-phase system. This geometry is illustrated in Figure 1 for the case n ) 2, where the reduced molar Gibbs energy of mixing ∆Gmix/ RT is chosen as the abscissa. The Gibbs energy of mixing is
xi(µi(P,T,{xi}) - µpure (P,T)) ∑ i i)1
(4)
It is often viewed as a combination of the ideal Gibbs energy of mixing and the excess Gibbs energy (of mixing):
∆Gmix(P,T,{xi}) ) ∆Gideal + Gexcess
(5)
The ideal Gibbs energy of mixing is
xi ln xi ∑ i)1
)]
(7)
φˆ iR(P,T,{xiR})xiR ) φˆ iβ(P,T,{xiβ})xiβ, i ) 1, ..., n
(8)
From eqs 6 and 7, the excess Gibbs energy can be expressed as n
[(
xi ln ∑ i)1
Gexcess(P,T,{xi}) ) RT
)]
φˆ i(P,T,{xi}) φpure (P,T) i
(9)
Tangent planes to the excess Gibbs energy curve at these same compositions xiR and xiβ are not common (collinear). Rather, these tangent planes at points R and β separately locate a pair of intercepts:
ln
[
] [
φˆ iR(P,T,{xi}) φpure (P,T) i
* ln
]
φˆ iβ(P,T,{xi}) φpure (P,T) i
(10)
Equivalently, the values for the fugacity coefficients are unequal:
n
∆Gideal(T,{xi}) ) RT
φpure (P,T) i
where φi is the fugacity coefficient fi/Pxi. The common tangent plane to the total Gibbs energy curve at points [phases] R and β therefore yields the intercept equalities
n
∆Gmix(P,T,{xi}) )
fpure (P,T) i
φˆ i(P,T,{xi})xi
n
) RT
)]
ˆfi(P,T,{xi})
(6)
It can easily be seen that this function will become large for a system of many components. For example, for a
φˆ iR(P,T,{xi}) * φˆ iβ(P,T,{xi})
(11)
This is a well-understood result. Figure 2 illustrates the noncollinearity of tangents to the excess Gibbs energy
Ind. Eng. Chem. Res., Vol. 40, No. 26, 2001 6215 Table 1. Liquid Mole Fraction L and the Compressibility Factors z for the Liquid and Vapor Phases for the Liquid-Vapor Flash Problem as a Function of the Number n of Finite Laguerre-Gauss Quadrature Pseudocomponents Used (T ) 350 K and P ) 10 bar)a z components n
L
liquid
vapor
bubble point P (bar)
2 3 4 5 6 7 8 9 10 11 12
0.589 74 0.588 31 0.586 78 0.585 51 0.584 64 0.584 13 0.583 86 0.583 75 0.583 70 0.583 69 0.583 69
0.117 50 0.117 96 0.118 17 0.118 33 0.118 44 0.118 51 0.118 55 0.118 57 0.118 57 0.118 58 0.118 58
0.972 43 0.971 77 0.971 20 0.970 79 0.970 57 0.970 46 0.970 42 0.970 42 0.970 42 0.970 42 0.970 43
72.725 72.800 72.890 72.953 72.990 73.009 73.017 73.022 73.021 73.022 73.022
a The mixture is 25 mol % CH , 12.5 mol % C H , 12.5 mol % 4 2 6 C3H8, and 50 mol % of the C7 continuous mixture described in appendices C and D. Also shown is the bubble-point pressure at T ) 350 K as a function of n.
Figure 2. SRK EOS prediction of the reduced excess Gibbs energy of mixing versus mole fraction of n-hexane for the binary mixture propane + n-hexane at 420 K and 20 bar (kij ) 0.0). Tangents at the equilibrium phase compositions are shown as dashed lines. The thin solid curve is the reduced ideal Gibbs energy of mixing.
curves at the equilibrium phase compositions for the two-phase system in Figure 1. The quadrature convergence (as n increases) of the excess Gibbs energy function for a continuous system can be related to the computational convergence of the fluid-phase equilibrium problem. For example, for finite Laguerre-Gauss quadrature, the quadrature efficiently performs the integration of a function f(X) by representing it as a sum:11
∫0
C
f(X) e-X dX )
n
wif(Xi) ∑ i)1
(A4)
As written here, the function f(X) depends on X, herein related to the carbon number I, and is a single-variable curve; in turn, the set of values {f(Xi)} are quadratureselected points on that curve. If the fugacity coefficients φi in eq 11 transparently depended only on the carbon number of species i, then it would follow that they would also fall on a single curve for a continuous system of interest and so would their ratios, which are the wellknown vaporization ratios, or K values:
φˆ iR φˆ iβ
)
xiβ ≡ KiβR, i ) 1, ..., n xiR
(12)
However, the fugacity functions are not quite this simple in nature. The value of the fugacity of a species i depends on the compositional makeup of the entire phase, which might be conceptually decomposed as (1) the mole fraction of the species i and (2) the composition of the “background” mixture in which species i resides. As the level of quadrature increases with n, the background description acquires additional detail. Convergence of a fluid-phase equilibrium result would appear to require that background effects likewise converge, resulting in convergence of the vaporization ratios to a
single curve of K ) K(I) at a relatively modest value of n. It will be demonstrated here that this convergence scenario occurs. Fluid-Phase Equilibrium Computations To demonstrate the convergence that a typical semicontinuous fluid-phase equilibrium calculation experiences, a system composed of 25 mol % methane, 12.5 mol % ethane, 12.5 mol % propane, and 50 mol % of a C7+ continuous mixture (described in and discretized according to appendices C and D) was flashed at 350 K and 10 bar, using the Soave-Redlich-Kwong (SRK) EOS as the fluid-phase model (appendix B). The results are shown in Table 1. One sees convergence of the liquid-phase fraction L to three significant figures at about n ) 7. The compressibility factors of the liquid and vapor phases zL and zV converge to four significant figures at about n ) 8. Separate calculations of the systems’ bubble-point pressures converge to three significant figures at about n ) 5. These results are consistent with the fluid-phase equilibrium results reported in earlier studies;9 other investigators who limited their quadrature levels to n e 4 were probably justified in claiming reasonably good results. Figure 3 is a plot of all of the equilibrium K values that this phase equilibrium problem produced for n ) 2-12. {Ki} appear to converge to a single problemspecific curve, depending only on the carbon number, which suggests that the variability of the background makeup (i.e., the composition {xj}, j ) 1, ..., n, j * i) is a secondary effect that does not significantly retard the quadrature convergence of the problem. The proposition that background effects may well be secondary, even negligible, has often been formalized in approximate statements such as the Lewis-Randall rule:
ˆfi(P,T,{xi}) ) f pure (P,T) xi i
(13)
A rough assessment of this background effect, convergence-wise, can be obtained by examining the convergence behavior of a function free of background effects, such as the pseudocritical temperature of the C7+ portion of the mixture, when calculated using the single equation (C1) for correlating Tc for the species.
6216
Ind. Eng. Chem. Res., Vol. 40, No. 26, 2001 Table 2. Crystal-Point Temperature for the Solid Solution Assumption versus the Pure Solid-Phase Assumption as a Function of the Number n of Finite Laguerre-Gauss Quadrature Pseudocomponents Used, for a C7+ Fluid Extending from Carbon Number 7 to 100 with a Molecular Weight of 200 (P ) 1 bar) crystal point T (K)
Figure 3. K values of the quadrature components for n ) 2-12 for the liquid-vapor flash problem at T ) 350 K and P ) 10 bar. Ki ) yi/xi. The mixture is 25 mol % CH4, 12.5 mol % C2H6, 12.5 mol % C3H8, and 50 mol % of a C7+ continuous mixture described in appendices C and D.
Pseudocritical Tc converges to four significant figures by n ) 5.10 It is recognized that the convergence performance will depend on the degree of nonlinearity of a problem; also, comparison of two different problems such as the bubble-point pressure and pseudocritical temperature is interpretation-wise risky. These concerns aside, it appears that compositional background effects retarded the convergence by about a count of 2 in the quadrature level employed. Solid-Fluid-Phase Equilibrium Problem Some aspects of the solid-fluid-phase equilibrium problem for the C7+ continuous mixture above were previously studied by us,10 employing the SRK EOS fluid-phase description coupled to a simple mathematical artifice (eq B2) for introducing the solid-phase fugacity function (because such classical equations of state do not directly admit a solid phase). The quadrature in appendix D and property correlations in appendix C were used. Of interest was the thermodynamic crystal point, which is similar to but not the same as the “cloud point”, which marks the visible precipitation of solids, often from a metastable dissolved state. The thermodynamic crystal point provides a computational upper bound to the temperature (hereafter, “crystalpoint temperature”) at which solid formation occurs. The crystal-point temperature problem was cast as the solution to the equations
ˆfil(P0,T,{xi}) ) fis°(P0,T), i ) n ˆfil(P0,T,{xi}) < fis°(P0,T), i ) 1, 2, ..., n - 1
(14)
where the heaviest pseudocomponent of the quadraturediscretized mixture has formed a trace of a (assumed) pure solid phase. (In every calculation performed, the heaviest species was thermodynamically the first to form a solid phase as the temperature decreased.) It was noteworthy that the computed crystal-point temperature did not achieve a clearly identifiable asymptotic crystal-point temperature as n f 12; rather, the tem-
components n
ideal solid solution
pure solid phase
2 3 4 5 6 7 8 9 10 11 12
328.084 354.091 371.959 381.538 385.749 387.445 388.083 388.301 388.367 388.384 388.387
327.506 353.741 371.798 381.430 385.642 387.301 387.846 387.894 387.712 387.424 387.088
perature appeared to pass through a maximum value of 387.894 K at n ) 9, after which it decreased as n increased. See the rightmost column of Table 2. Although the lack of convergence was not particularly pronounced for the case of B ) 100.5, it was much more pronounced for the problems B ) 40.5 and 1000.5. The fluid-phase equilibrium computations above suggest a reason for this apparent failure to converge. The fluidsolid-phase equilibrium K value for species n does not display a single problem-specific K-value curve (such as that in Figure 2), as was the case for the fluid-phase equilibrium flash and bubble-point calculations. The definition of this solid-fluid-phase equilibrium K value for species n would be
Knfs ) xn/sn
(15)
where sn is the mole fraction of the heaviest species n in the solid phase. Because sn was assumed to be unity in eq 14, the value of Kn is synonymous with xn, the liquid-phase mole fraction of the heaviest pseudocomponent from a given quadrature discretization. As can be seen from Table 5, xn becomes smaller as n increases, which would be expected from the finite Laguerre-Gauss quadrature for a fixed set of parameters (MW, A, B). The carbon number of the heaviest species for large n is relatively unchanged, approaching 100.5 as the quadrature level n f ∞. However, the mole fraction xn continues to decrease with increasing quadrature level. The result is that there is no clearly definable problem-specific K-value curve, such as that in Figure 2. Because one expects in nature a unique reproducible answer to the crystal-point temperature determination problem, the problem description in eq 14, appealing as it is, will not provide this answer in the context of the quadrature problem. Rather, it is suggested that a quadrature-based computation of the type used herein to determine the crystal-point temperature should admit the possibility of solid solution formation. Turner12 has pointed out in his review of n-paraffin phase equilibrium behavior that n-paraffins of neighboring carbon number will form solid solutions, the crystalline nature of which can be complicated and variable, even at carbon numbers of 30. Presumably, the potential for solid solution formation would increase as the carbon number increases. Admitting the possibility of a solid solution phase mathematically puts the crystal-point temperature
Ind. Eng. Chem. Res., Vol. 40, No. 26, 2001 6217 Table 3. Liquid-Phase Mole Fraction L, Compressibility Factors z, and Pseudocritical Temperatures T of the Liquid and Solid Solution Phases as a Function of the Number n of Finite Laguerre-Gauss Quadrature Pseudocomponents Used, for a C7+ Fluid Extending from Carbon Number 7 to 100 with a Molecular Weight of 200 (T ) 350 K and P ) 1 bar)a components n
L
z(liquid)
Tc,liq (K)
Tc,sol (K)
3 4 5 6 7 8 9 10 11 12
0.9925 0.9988 0.9965 0.9982 0.9977 0.9979 0.9979 0.9978 0.9979 0.9978
0.012 99 0.013 25 0.013 17 0.013 23 0.013 22 0.013 22 0.013 22 0.013 22 0.013 22 0.013 22
665.45 666.34 665.70 666.04 665.92 665.96 665.97 665.94 665.97 665.94
901.30 915.79 902.20 910.66 904.23 908.43 905.68 907.43 906.36 906.99
a The solid phase is assumed to be an ideal solution phase. No solid phase forms for n ) 2 at this temperature and pressure.
Figure 5. Composition profiles {si} of the solid solutions for n ) 4, 8, and 12 for the solid-liquid-phase equilibrium problem at T ) 350 K and P ) 1 bar, where the solid phase is an ideal solid solution. The mixture is the C7+ continuous mixture described in appendices C and D. Table 4. Ideal Solid Solution Phase Composition at Finite Laguerre-Gauss Quadrature Levels of n ) 4, 8, and 12 Pseudocomponents, for a C7+ Fluid Extending from Carbon Number 7 to 100 with a Molecular Weight of 200 (T ) 350 K and P ) 1 bar) components ) 4 carbon species number 1 2 3 4 5 6 7 8 9 10 11 12
Figure 4. K values of the quadrature components for n ) 2-12 for the solid-liquid equilibrium problem at T ) 350 K and P ) 1 bar, where the solid phase is an ideal solid solution. Ki ) si/xi. The mixture is the C7+ continuous mixture described in appendices C and D.
computation on the same footing as a bubble point or dew point in fluid-phase equilibria. All fluid-phase species in the dominant phase would be present in the infinitesimal solid phase. This description should lead to a quadrature-converged result, even without knowledge of the appropriate (correct) nonideal solid solution description for the homologous series of n-paraffins that composes the C7+ mixture described in appendix D. For example, one should be able to achieve reasonably rapid quadrature convergence of the crystal-point temperature computation with increasing quadrature level n, even with a relatively naı¨ve solid solution model. The overriding consideration with respect to computational convergence, in our opinion, is the admissibility of the solid solution formation (valid or otherwise). To demonstrate the need for admitting solid solution behavior into solid-fluid multicomponent computations, we performed the computation using the ideal solid solution model. Table 2 compares the crystal-point temperature as a function of the quadrature level for the cases of the formation of an ideal solid solution versus a pure solid. As can be seen, the crystal-point temperature calculation converges uniformly with increasing quadrature level n to a value of about 388.39
8.8654 19.2523 39.2984 72.0888
si
components ) 8
components ) 12
carbon number
carbon number
si
7.0598 9.4439 13.7093 19.8095 27.6598 37.1156 47.9363 59.7278 71.8604 83.3742 92.9342 98.9899
0.0052 0.0108 0.0140 0.0160 0.0210 0.0436 0.1866 0.5307 0.1375 0.0274 0.0059 0.0013
si
0.0203 7.5134 0.0092 0.0343 11.8438 0.0177 0.1357 19.6444 0.0219 0.8097 30.8978 0.0357 45.4573 0.1839 62.7762 0.6630 81.1619 0.0633 96.1163 0.0054
K. As pointed out earlier, the computations assuming the formation of a pure solid phase do not. Additional calculations were made for this mixture at a temperature of 350 K, which is below the crystalpoint temperature. These calculations are similar in nature to vapor-liquid flash calculations, with a finite amount of solid phase forming in the presence of a liquid phase. Table 3 lists the liquid fraction L and its associated compressibility factor z and pseudocritical temperature Tc as a function of quadrature level n. These satisfactorily converge as would be expected. A considerably more sensitive test of the quadrature convergence is the pseudocritical temperature of the solid solution phase, which appears to be converging in a damped oscillatory manner. The K values of the solid-fluid-phase equilibrium problem converge to a single curve (Figure 4) like those of the fluid-phase equilibrium problem (Figure 3). The composition of the solid-phase becomes distributed across the pseudocomponents as the quadrature level increases (Table 4 and Figure 5). Discussion Whereas the ideal solid solution model, for example, may in reality be a questionable representation of the laboratory solid-phase behavior, its effect on crystalpoint temperature convergence in the quadrature prob-
6218
Ind. Eng. Chem. Res., Vol. 40, No. 26, 2001
lem is significant and unequivocal. Our computations suggest that some degree of solid solution behavior must be included in quadrature-based solid-fluid-phase equilibrium descriptions in order to ensure computational convergence of the results as the quadrature level increases. Although operating under the assumption of pure solid-phase formation is often appealing from a conceptual viewpoint, the resultant mathematical inflexibility prevents computational convergence. Furthermore, if one assumes that one truly does have a “continuous” compositional system of the type candidate for quadrature characterization, the assumption of pure solid-phase formation is likely invalid, as Turner’s review12 suggests for the homologous series of n-paraffins. Convergence of the crystal-point temperature problem with an ideal solid solution model determining solidphase fugacity functions portends favorable convergence of the problem using (fugacity functions from) a realistic solid solution model. We anticipate that introduction of nonideality into the above ideal solid solution description will not significantly hinder the convergence of associated phase equilibrium results. For example, one can envision solid solutions of higher purity (i.e., narrower peaks than those shown in Figure 5) in equilibrium with the pseudoized fluid phase. The authors are currently studying problems in solid-fluid-phase equilibria using nonideal solid solution descriptions.
I is related to an n-paraffinic molecular weight through
MW ) 14.0I + 2.0
(A2)
The average carbon number for a given carbon number range (A, B) can be expressed in terms of the exponential decay parameter R through
〈I〉 ) [(RA + 1)e-RA - (RB + 1)e-RB]/R(e-RA - e-RB) (A3) In the limit B f ∞, 〈I〉 f A + 1/R. The higher the molecular weight MW of the continuous mixture, the smaller R is or, equivalently, the flatter the exponential distribution F(I) becomes. Briefly, finite Laguerre-Gauss quadrature approximates in an optimal way integrals of the type
∫0 f(X) e C
-X
n
dX )
wif(Xi) ∑ i)1
(A4)
The sets {wi} and {Xi} are determined by maximizing the polynomial order that f(X) can rigorously be. For a chosen n, the integral on the left-hand side is represented precisely to order 2n - 1.11 Appendix B. SRK EOS Formalism
Conclusions For a quadrature-characterized continuous mixture with increasing quadrature level n, the following apply: (i) Convergence of the phase equilibrium problem is concomitant with the convergence of the excess Gibbs energy function. (ii) Specifically, this functional convergence can be viewed as the establishment of a fixed equilibrium ratio curve K ) K(Cn), the points of which are represented by the set {Ki}, i ) 1, ..., n. (iii) Physically, it is understood that Ki depends primarily on the identity of the species carbon number i and secondarily on the nature of the background formed by the other species of the mixture. (iv) Convergence of the solid-fluid-phase equilibrium problem can be seriously hindered by assuming that the solid phase is pure; rather, adopting a solution model for the solid phase promotes convergence of the solidfluid-phase equilibrium results.
The continuous C7+ portion is discretized as a mixture of pseudocomponents using finite Laguerre-Gauss quadrature, as in appendix A. The resulting C7+ pseudocomponents comprise part or all of the fluid phase in a particular phase equilibrium problem. Fluidphase fugacities are determined from the SRK EOS model:16
P ) RT/(v - b) - a(T)/v(v + b) n
b) n
a)
In this study we use the simple exponential distribution
∑ ∑xixjaij i)1 j)1 n
xibi ∑ i)1
bi ) 0.08664 n
(A1)
to describe a continuous fluid, with its carbon number I limited to the range A ) 6.5 to B ) 100.5. Equation A1 qualitatively represents the nature of typical crude oil C7+ portions, with their composition functions exponentially decaying with the carbon number. (Du and Mansoori13,14 were the first, to our knowledge, to use the simple exponential distribution in phase equilibrium calculations, although Shibata et al.8 point out that Vogel et al.15 earlier presented experimental evidence supporting the use of the exponential form for C7+ mole fraction distributions in crude oils.) The carbon number
(B1a)
where
F(I) ) ae-aI/(e-aA - e-aB) ) 0, A g I, I g B
n
b)
Appendix A. Quadrature Formalism
xibi ∑ i)1
a(T) )
RTci Pci
n
∑ ∑xixjaij i)1 j)1
aij ) xaiiajj aii ) Riaci aci(Tci) ) 0.42747
R2Tci2 Pci
Ri1/2 ) 1 + mi(1 - TRi1/2) mi ) 0.480 + 1.574ωi - 0.176ωi2
(B1b)
Ind. Eng. Chem. Res., Vol. 40, No. 26, 2001 6219 Table 5. Compositions of the C7+ Continuous Mixture of n-Paraffins, Having Molecular Weight 200 with A ) 6.5 and B ) 100.5, as a Function of Quadrature Level n for the Finite Laguerre-Gauss Quadrature Prescription n)2 n)3 n)4
n)5
n)6
n)7
n)8
n)9
species
carbon number
mole fraction
1 2 1 2 3 1 2 3 4 1 2 3 4 5 1 2 3 4 5 6 1 2 3 4 5 6 7 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 9
10.970 800 32.547 169 9.640 376 23.797 524 53.612 087 8.865 414 19.252 258 39.298 370 72.088 773 8.346 969 16.350 261 31.251 578 53.908 798 84.098 287 7.981 271 14.353 291 26.007 525 43.203 216 65.840 619 90.635 746 7.714 218 12.915 609 22.337 180 36.035 742 53.868 196 74.771 253 94.135 151 7.513 400 11.843 770 19.644 440 30.897 814 45.457 347 62.776 195 81.161 861 96.116 327 7.358 606 11.022 028 17.599 978 27.048 380 39.233 652 53.807 895 69.927 331 85.676 895 97.317 634
8.529 85 × 10-1 1.470 15 × 10-1 7.059 15 × 10-1 2.827 13 × 10-1 1.137 27 × 10-2 5.858 49 × 10-1 3.671 02 × 10-1 4.605 57 × 10-2 9.933 02 × 10-4 4.883 91 × 10-1 4.088 02 × 10-1 9.554 48 × 10-2 7.101 49 × 10-3 1.603 76 × 10-4 4.103 35 × 10-1 4.201 01 × 10-1 1.466 21 × 10-1 2.151 13 × 10-2 1.384 16 × 10-3 4.732 20 × 10-5 3.480 43 × 10-1 4.124 77 × 10-1 1.904 08 × 10-1 4.343 65 × 10-2 5.244 00 × 10-3 3.707 06 × 10-4 2.063 55 × 10-5 2.980 73 × 10-1 3.943 11 × 10-1 2.234 64 × 10-1 6.961 66 × 10-2 1.287 25 × 10-2 1.518 56 × 10-3 1.330 16 × 10-4 1.132 81 × 10-5 2.576 36 × 10-1 3.710 01 × 10-1 2.457 45 × 10-1 9.652 17 × 10-2 2.432 30 × 10-2 4.172 88 × 10-3 5.335 52 × 10-4 6.010 84 × 10-5 7.150 35 × 10-6
Binary interaction parameters can be used for aij, i * j; in this solely computational study, wherein no comparisons are made with experimental data, no binary interaction parameters (corrections) are introduced. Pure solid-phase fugacities are determined using an artifice that relates them to pure liquid-phase fugacity functions: o
ln(fil /fiso)
)
(∆hfus i /RTti)
(
Tti -1 T
)
n ) 10
n ) 11
n ) 12
species
carbon number
mole fraction
1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 11 1 2 3 4 5 6 7 8 9 10 11 12
7.236 770 10.377 503 16.006 195 24.071 635 34.457 388 46.917 133 60.946 926 75.550 811 88.889 840 98.092 250 7.139 160 9.862 315 14.737 250 21.713 562 30.691 585 41.486 158 53.761 194 66.919 810 79.941 139 91.215 492 98.617 905 7.059 756 9.443 870 13.709 279 19.809 501 27.659 784 37.115 628 47.936 322 59.727 769 71.860 408 83.374 238 92.934 154 98.989 879
2.245 90 × 10-1 3.458 04 × 10-1 2.587 47 × 10-1 1.214 99 × 10-1 3.873 45 × 10-2 8.830 84 × 10-3 1.533 18 × 10-3 2.241 75 × 10-4 3.220 53 × 10-5 4.937 72 × 10-6 1.973 19 × 10-1 3.205 98 × 10-1 2.644 03 × 10-1 1.429 99 × 10-1 5.482 29 × 10-2 1.561 53 × 10-2 3.467 93 × 10-3 6.421 82 × 10-4 1.095 51 × 10-4 1.952 12 × 10-5 3.625 62 × 10-6 1.745 99 × 10-1 2.964 07 × 10-1 2.645 62 × 10-1 1.603 84 × 10-1 7.130 45 × 10-2 2.428 04 × 10-2 6.587 95 × 10-3 1.495 47 × 10-3 3.038 37 × 10-4 6.048 41 × 10-5 1.294 40 × 10-5 2.782 12 × 10-6
pseudocomponents are n-paraffins for the sake of calculation demonstrations. The following correlations were used for the pseudocomponent parameters:
Critical temperature (K):17 Tc ) 959.98 + exp(6.81536 - 0.211145nC2/3)
(C1)
Critical pressure (bar):10 (B2)
Equation B2 is a simplification of the form in work by Prausnitz et al.18 Tti is the triple-point temperature of is the heat fusion of pseudocomponent i, and ∆hfus i pseudocomponent i at the triple point. The pure liquid (hypothetical) fugacity fis° is calculated from the SRK EOS. Appendix C. Pseudocomponent Parameter Correlations To carry out phase equilibrium computations that include fluid and/or solid phases containing the quadrature-generated pseudocomponents, one must have the parameters Tc, Pc, ω, Tt, and ∆hfus as a function of pseudocomponent carbon. It was assumed that these
Pc ) 191.8823nC-0.9657
(C2)
Acentric factor:10 ω ) 1.85 - exp(0.68 - 0.05nC0.88)
(C3)
Triple-point temperature (K):18 Tt ) 374.5 + 0.02617[12.01115nC + 1.00797(2nC + 2)] - 20172/[12.01115nC + 1.00797(2nC + 2)] (C4) Enthalpy of fusion at the triple point temperature (cal/gmol):18 ∆hfus ) 0.1426[12.01115nC + 1.00797(2nC + 2)]Tt (C5)
6220
Ind. Eng. Chem. Res., Vol. 40, No. 26, 2001
Appendix D. Prototype Continuous Mixture and Its Discretization Quadrature conforming to the finite Laguerre-Gauss prescription in appendix A was used to discretize a prototype continuous C7+ portion of an oil extending from carbon number 7 to 100 with a molecular weight of 200. The carbon number bounds were specified as A ) 6.5 and B ) 100.5 (which determine the exponential decay parameter R in eq A1 using eqs A2 and A3). The upper limit of the finite integral in eq A4 is C ) R(B A), which for this prototype mixture is 0.130 834. A choice of n is made (depending on how precise one wishes the integral representation to be), and the resultant sets {wi} and {Xi} are translated into sets of mole fraction and carbon numbers. For example, the normalized {wi} are the mole fractions for the C7+ portion, which correspondingly have the carbon numbers {Ii}. This discrete description is used to perform the phase equilibrium computation. The precision with which one can calculate the sets of mole fractions and carbon numbers can be progressively challenging as n is increased. These sets must be calculated accurately in order to evaluate the convergence of a problem’s answer with increasing n. Using double precision in this study, the authors were able to generate reliable values of these sets for (A, B, MW) ) (6.5, 100.5, 200) for values of n to 12. These discretizations are given in Table 5. Nomenclature a(T) ) SRK parameter A ) lower carbon number bound b ) SRK parameter B ) upper carbon number bound C ) upper bound of the reduced carbon number variable f ) fugacity F ) distribution function ∆hfus ) enthalpy of fusion I ) carbon number mi ) SRK parameter MW ) molecular weight n ) number of quadrature pseudocomponents Cn ) carbon number P, P0 ) pressure R ) gas constant T ) temperature v ) molar volume wi ) quadrature weighting factor X ) quadrature variable xi ) mole fraction R ) exponential distribution decay parameter Ri ) SRK parameter ωi ) acentric factor Subscripts c ) critical point C ) carbon fus ) fusion i, j ) species l ) liquid phase R ) reduced s ) solid phase t ) triple point
Superscripts ° ) pure phase mix ) mixing property
Literature Cited (1) Cotterman, R. L.; Prausnitz, J. M. Flash Calculations for Continuous or Semicontinuous Mixtures Using an Equation of State. Ind. Eng. Chem. Process Des. Dev. 1985, 24, 434-443. (2) Cotterman, R. L.; Prausnitz, J. M. Continuous Thermodynamics for Phase-Equilibrium Calculations in Chemical Process Design. In Kinetic and Thermodynamic Lumping of Multicomponent Mixtures; Astarita, G., Sandler, S. I., Eds.; Elsevier Science Publishers BV: Amsterdam, The Netherlands, 1991; p 229. (3) Cotterman, R. L.; Bender, R.; Prausnitz, J. M. Phase Equilibria for Mixtures Containing Very Many Components. Development and Application of Continuous Thermodynamics for Chemical Process Design. Ind. Eng. Chem. Process Des. Dev. 1985, 24, 194-203. (4) Cotterman, R. L.; Chou, G. F.; Prausnitz, J. M. Comments on “Flash Calculations for Continuous or Semicontinuous Mixtures Using an Equation of State”. Ind. Eng. Chem. Process Des. Dev. 1986, 25, 840-841. (5) Cotterman, R. L.; Dimitrelis, D.; Prausnitz, J. M. Design of Supercritical Fluid-Extraction Processes Using Continuous Thermodynamics. In Supercritical Fluid Technology; Penninger, J. M. L., McHugh, M. A., Krukonis, V. J., Eds.; Elsevier Science Publishers BV: Amsterdam, The Netherlands, 1985; pp 107-120. (6) Behrens, R. A.; Sandler, S. I. Use of Semicontinuous Description to Model the C7+ Fraction in Equation of State Calculations. Presented at the SPE/DOE Fifth Symposium on Enhanced Oil Recovery of the Society of Petroleum Engineers, Tulsa, OK, Apr 20-23, 1986; SPE/DOE 14925. (7) Behrens, R. A.; Sandler, S. I. The Use of Semicontinuous Descriptions to Model the C7+ Fraction in Equations of State Calculations. SPE Reservoir Eng. 1988, 3, 1041-1047. (8) Shibata, S. K.; Sandler, S. I.; Behrens, R. A. Phase Equilibrium Calculations for Continuous and Semicontinuous Mixtures. Chem. Eng. Sci. 1987, 42, 1977-1988. (9) Luks, K. D.; Turek, E. A.; Kragas, T. K. Asymptotic Effects Using Semicontinuous vis-a`-vis Discrete Descriptions in Phase Equilibrium Computations. Ind. Eng. Chem. Res. 1990, 29, 21012106. (10) Labadie, J. A.; Luks, K. D. Convergence of the Continuous Fluid-Solid Equilibrium problem Using Quadrature Compositional Characterizations. Ind. Eng. Chem. Res. 2000, 39, 790796. (11) Lapidus, L. Digital Computation for Chemical Engineers; McGraw-Hill: New York, 1962; Chapter 2, pp 56-66. (12) Turner, W. R. Normal Alkanes. Ind. Eng. Chem. Prod. Res. Dev. 1971, 10 (3), 238-260. (13) Du, P. C.; Mansoori, G. A. Phase Equilibrium Computational Algorithms of Continuous Mixtures. Fluid Phase Equilib. 1986, 25, 840-841. (14) Du, P. C.; Mansoori, G. A. Phase Equilibrium of Multicomponent Mixtures: Continuous Mixture Free Energy Minimization and Phase Rule. Chem. Eng. Commun. 1987, 54, 139-148. (15) Vogel, J. L.; Turek, E. A.; Metcalfe, R. S.; Bergman, D. F. Applications of Equations of State to Calculate Reservoir Fluid Properties. Fluid Phase Equilib. 1983, 14, 103-116. (16) Soave, G. Equilibrium Constants From a Modified RedlichKwong Equation of State. Chem. Eng. Sci. 1972, 27, 1197-1203 (17) Prausnitz, J. M.; Lichtenthaler, R. N.; de Azevedo, E. G. Molecular Thermodynamics of Fluid-Phase Equilibria, 3rd ed.; Prentice-Hall: Englewood Cliffs, NJ, 1999; Chapter 11. (18) Tsonopoulos, C. Critical Constants of Normal Alkanes from Methane to Polyethylene. AIChE J. 1987, 33, 2080-2083. (19) Thomas, F. B.; Bennion, D. B.; Bennion, D. W.; Hunter, D. E. Experimental and Theoretical Studies of Solids Precipitation from Reservoir Fluids. J. Can. Pet. Technol. 1992, 31 (1), 22.
Received for review April 19, 2001 Revised manuscript received September 18, 2001 Accepted September 20, 2001 IE0103512
