- Email: [email protected]
Critical point calculations for semi-continuous mixtures
FLUIDPHM[ EOUIIBBIIIA ELSEVIER
Fluid Phase Equilibria 139 (1997) 137-153
Critical point calculations for semi-continuous mixtures G e r a l d o L. R o c h o c z a,b, M a r c e l o C a s t i e r b,~, S t a n l e y I. S a n d l e r c,, Chemtech Sert,i~'os de Engenharia, Rua Nilo Pe~'anha 51, sala 2106, Rio de Janeiro, RJ, Brazil h Escola de Quhnica, Unit~ersidade Federal do Rio de Janeiro, C.P. 68542, Rio de Janeiro, RJ, 21949-900, Brazil ~ Center/or Molecular and Engineering Thermodynamics, Department of" Chemical Engineering Unit,ersi O, of Delaware Newark, DE 19716, USA
Abstract A typical problem in the oil and gas industry is the calculation of high pressure phase equilibria in systems containing a very large number of similar components. Treating the composition of these mixtures as a continuous distribution, instead of trying to characterize each individual component, has proven to be an elegant and efficient approach to solve this type of problem. In this paper, we rigorously formulate the critical point conditions for continuous and semi-continuous mixtures, but as is often the case with such formulations, the complexity of the expressions prevents an analytical solution of the equilibrium problem, and discretization is required in order to use realistic thermodynamic models. Once in discrete form, the critical point equations were solved using a modified form of the Hicks and Young algorithm and the thermodynamic stability of the solutions was guaranteed by the calculation of higher terms in the Helmholtz free energy expansion and by global phase stability analyses. The resulting procedure was employed to perform an investigation of the influence of the distribution parameters on the types of critical phase diagrams that can be obtained from the Peng-Robinson equation of state. © 1997 Elsevier Science B.V. Keywords: Critical state; Vapor-liquid equilibria; Method of calculation; Equation of state; Continuous thermodynamics;
Mixture
1. I n t r o d u c t i o n The characterization o f mixtures, such as oil reservoir fluids, coal polymer solutions, for physical property and equilibrium calculations very large number of components. Instead of trying to identify each these mixtures as an infinite set of components characterized by one
derived liquids and polydisperse may be difficult because of their individual component, modeling or more continuous variables, in
* Corresponding author. On leave at the Department of Chemical Engineering, University of Delaware, Newark, DE, USA. 0378-3812/97/$17.00 © 1997 Elsevier Science B.V. All rights reserved. PII S 0 3 7 8 - 3 8 1 2 ( 9 7 ) 0 0 2 1 2 - 4
138
G.L. Rochocz et al. / Fluid Phase Equilibria 139 (1997) 137-153
an approach known as the thermodynamics of continuous mixtures, has proven to be an elegant and efficient approach to solve this type of problem. A hybrid description, known as semi-continuous mixture thermodynamics, is to characterize some of the components individually, usually the lighter components, and to group the remaining species into one or more continuous families. In typical applications of the thermodynamics of continuous mixtures, summations over the mixture components are recast into integrals that often need to be numerically evaluated. The use of efficient integration methods generally results in computational implementations that compare favorably in terms of accuracy and speed with the traditional approach of lumping species into arbitrarily chosen pseudocomponents (Pedersen et al. [1-3]). The idea of using continuous distributions to characterize mixtures dates back to the work of Bowman [4] and Edmister and Buchanan [5], but the rigorous formulation of the thermodynamics of continuous mixtures based on functional analysis, found in the work of R~itzsch and Kehlen [6,7] and Kehlen et al. [8] renewed interest in this field. Work in this area has focused on the development of algorithms to solve the new phase equilibrium functional equations, the characterization of the continuous mixtures and the use of thermodynamic models in continuous form. Much of the work done upto 1990 was reviewed by Cotterman and Prausnitz [9] and Rochocz [10], and new applications of the continuous mixture approach have been published in recent years. Halpin and Quirke [11] used distribution functions for continuous mixtures that are written as an expansion in orthonormal functions, which allows the mass balance to be preserved in flash calculations. Haynes and Mathews [ 12] used the experimental TBP curve directly as distribution function to calculate VLE, and this approach was later incorporated into a reservoir simulator [13]. Luks et al. [14] studied the effect of varying the choice of pseudocomponents in each stage of a multiple-contact separation process. Kramarz and Wyczesany [15] formulated and solved the chemical equilibria of continuous mixtures. Lira-Galeana et al. [16] employed a continuous mixture approach to model the segregation of heavy components due to gravity in oil reservoirs. Despite the large amount of work done on the calculation of phase equilibria in continuous and semi-continuous systems, comparatively few papers have dealt with the determination of phase stability in these systems. The stability of polydisperse polymer systems was studied by Kehlen et al. [17] and Browarzik et al. [18], using segment fractions as independent variables in the formulation. Hu et al. [19] followed this approach to compute spinodal curves, critical points and other phenomena in polymer solutions. More recently, Browarzik and Kehlen [20] developed the theoretical framework for a continuous mixture stability theory that was applied to the van der Waals equation of state (EOS), using one-fluid mixing rules with the binary interaction parameters kij s e t equal to zero. In the van der Waals EOS, they further assumed that ~/a and b of each fraction were linear functions of a single characterization parameter M, chosen to be the number of carbon atoms in the components of the n-alkane mixture used in their procedure. Under these hypotheses, they determined spinodal curves and critical points for these mixtures. In this paper, we are first interested in formally extending the criticality conditions of Heidemann and Khalil [21] to continuous mixtures, using elements of functional analysis. We then investigate the effect of the distribution parameters on the types of critical loci that are obtained. We consider pseudobinary mixtures containing a light substance and a continuous family of components modeled by the Peng-Robinson [22] equation of state with the one-fluid mixing rules and non-zero binary interaction parameters. Under these hypotheses, the integrals that appear in the semi-continuous formulation do not have analytical solutions and we resort to numerical methods to obtain their
G.L. Rochocz et al. / Fluid Phase Equilibria 139 (1997) 137-153
139
values. The criticality conditions in the discretized variables are solved using a modified version [23,24] of the Hicks and Young algorithm [25].
2. Critical points in continuous mixtures The Helmholtz free energy of a continuous mixture is a state functional and can be written [8] as:
A=A(T,V,n,F(I))
(1)
where n is the total number of moles, F ( I ) is an intensive distribution function which in general is a function of several characterization variables, but for simplicity we consider a single characterization variable I. The development that follows is also applicable to semi-continuous mixtures, because the discrete components can be represented in the formulation by delta functions in the distribution function. From the normalization condition we have that:
fF(I)dI=
1
(2)
Following Heidemann and Khalil [21] and using the Kehlen et al. [8] notation for state functional partial derivatives, the Taylor expansion of the continuous mixture Helmholtz free energy function A at constant temperature and total volume is:
A-AO-fltx(I)AF(I)dI
= ~" f/f/+
--
OntOnt+ r,v
03A OnlOnl+Onl.
AF(I)dlAF(I+)dI +
AF(I)dlAF(I+)dl + AF(I*)dI* +...
(3)
T,V
where the chemical potential in a continuous mixture corresponds to the first partial derivative of A with respect to the distribution function:
/x(I)=
(OA) On1 r,v --
(4)
The system is locally stable if the first non-vanishing term on the right hand side of Eq. (3) is of even order and positive, and for any arbitrary perturbations AF(I)dI and AF(I+)dI +, this most commonly occurs when:
flf+ ( 02A ) AF(I)dlAF(I+)dI+>O OnzOnt+ r,v
(5)
At a critical point, both the second and third order terms of the Helmholtz free energy expansion
140
G.L. Rochocz et al. / Fluid Phase Equilibria 139 (1997) 137-153
must vanish, i.e.:
..z ) On1Ont+ r,v flfl fl ( 03A ) AF(1)dIAF(I+)dl+AF(I*)dI * = 0 ' * OntOnI+ On1" r,v
(7)
The corresponding expressions for the quadratic (q) and cubic (c) forms in a mixture containing only discrete components are:
q =
Onjc3n i j= I i=l
c=
iii( k=l j=l
i=l
T.V. A l ' l i A n j n,.,.j
O'A cgnkcgnjOni
: 0
"r,v. AniAnj An~ =
(8)
0
(9)
n,,,,.~
for any perturbations An i, An i and An k. The zero of the quadratic form in the discrete case can be detected by a zero value of the determinant of the matrix Q whose elements are:
(..z)
QiJ = OnjOni r,v,
(10)
tlt ~ t,i
or by obtaining a zero value for at least one of the eigenvalues of this matrix. For the continuous case, if a discretization is carried out, the discrete formulae for criticality and the numerical calculation procedures associated with them can be used, as we do in the applications which follow. It is, however, interesting to comment on the theoretical framework that formally supports the extension of the analysis for systems with infinitely many components. This framework comes from the fact that many results of linear algebra can be extended to spaces of infinite dimension, where matrix algebra is replaced by functional analysis of linear operators. Here, we will only mention the most relevant aspects of this extension and will refrain from giving proofs or providing the complete development of the formalism. For this we refer to the discussion of spectral theory in infinite dimensional Hilbert spaces found in Epstein [26]. This extension starts by considering the integral operator T:C[a,b] ~ C[a,b], where C[a,b] is the normed space of all continuously differentiable functions on [a,b], which is defined by Tf= g so that: (Tf)(s) =
g(s) = f bk(s,t)f(t)dt
(11)
where k(s,t) is a continuous function defined on [a,b] × [a,b], called the kernel of T. It can be shown [26] that the integral operator is a compact linear operator. An inner product can be defined in C[a,b] as:
(12)
G.L. Rochocz et al. / Fluid Phase Equilibria 139 (1997) 137-153
141
which defines a Hilbert space with the norm Hf]] = ( ( f , f ) , and this allows the definition of a norm for the integral operator. With the definition of the inner product, we can write the quadratic form for the integral operator T as:
( f ,Tf)= £b~bk( t,s)f( t)f(
,)dtds
(13)
The second functional derivative of the Helmholtz free energy ({(02A)}/{(0ni0ns+)})r.v = A"(I + I ÷) is a symmetric kernel function for the integral operator that appears in the quadratic form; in such a case the integral operator defined by this function is said to be hermitian. Epstein [26] demonstrates that if all the eigenvalues of a compact linear and hermitian operator are positive, the quadratic form ( f, Tf) is positive-definite. Higher-order terms can also be written using this formalism. For instance, we can introduce a third order kernel h(t,s,u) that defines a new integral operator F:C[a,b] ~ C[a,b] × C[a,b] as:
( F f )( t,s) = g( t,s) = ~bh( t,s,u)f( u)du
(14)
The new integral operator Ff(u) provides a function g(t,s) that can be considered to be a second order kernel g(t,s) for an integral operator T. It is therefore possible to define a cubic form: rb rb rb
(f, Vf,f}=(f,((Vf),f}}=((( Vf),f},f}= J. ],, Ja h(t,s,u)f(t)f(s)f(u)dtdsdu
(15)
At the stability limit, the quadratic form vanishes which means that at least one eigenvalue Ai is equal to zero, i.e., Ai=O
(16)
When the quadratic form becomes positive-semidefinite, it is possible to find a perturbation distribution function AF(I) so that:
flA"( I,I +) AF( I +)dI +=
0
(17)
and
flAFZ( l)dl
1
(18)
given by the normalized eigenfunction corresponding to the vanishing eigenvalue. Eqs. (17) and (18) are the continuous analogs of
QAn=O
(19)
and /7 c
E
= 1
(20)
i=1
that are found in the work of Heidemann and Khalil [21]. In fact, Eqs. (16), (17) and (9) represent a rigorous extension to continuous mixtures of the Heidemann and Khalil [21] formulation for the calculation of critical points.
G.L. Rochocz et al. / Fluid Phase Equilibria 139 (1997) 137-153
142
Although theoretically rigorous, the continuous criterion for critical points is cumbersome for numerical calculations. For discrete mixtures, one vanishing eigenvalue is detected by a zero determinant. In the continuous case, the integral operator can have infinite eigenvalues and it is very difficult to conclude when one of them is equal to zero. Furthermore, the integrals often do not have analytical solutions when realistic thermodynamic models are used, and must be numerically approximated using a method such as Gaussian quadrature. From a practical point of view, this corresponds to selecting a group a pseudocomponents and their amount in a given mixture, but based on a sound mathematical criterion. Once this is done, programs usually used for the calculation of phase equilibria, physical properties or, in our particular application, critical points, with mixtures of discrete components can be used. To calculate critical points, we used a modified version of the Hicks and Young [25] algorithm; additional details about the implementation of this method can be found elsewhere [23,24]. To test the stability of the calculated critical points, we numerically evaluated the fourth-order term of the A expansion using a four-point difference formula and the global stability test [27]. Only critical points that passed the stability tests are reported here. However, the global stability test was only used to verify the possibility of additional fluid phases; the possibility of forming solid phases was not tested in this work.
3. Modeling of semi-continuous mixtures The formulation of critical points that has been presented is not limited to any particular distribution function and discrete components can be included in the formulation as delta functions in the distribution function. However, it is usual to distinguish explicitly between discrete species and continuous families of components and we will follow this approach in the rest of this paper. As simple examples, consider the fact that the mole fraction of all species should sum to unity
Y'~xi+~xk=
~_,xi+ Y'~ x k
i=1
i=1
=1
k=l
Fk(Ik)dI k = 1
(21)
k
and that the expression to determine the parameter b of the Peng-Robinson equation of state when the one-fluid mixing rule is used is
b = Y'~ xib i + ~ i=1
xk
(I~.)bk(lk)dl k
(22)
k=l
where the indexes i and k denote a discrete component and a continuous family respectively, x denotes mole fractions, and rl~ and ~bk are the lower and upper limits for the characterization variable I~ of family k. Again, for simplicity, we assume that each family is characterized by a single variable Ik; also the distribution functions are normalized
f,7~kFk( Ik)dlk= l
(23)
Several procedures exist to characterize oil fluids and the discussion that follows can be adapted to accommodate different characterization procedures, but we will follow Shibata et al. [28] and choose a bounded exponential function in which the characterization variable I is the number of carbon
143
G.L. Rochocz et al. / Fluid Phase Equilibria 139 (1997) 137-153
atoms, assumed to vary continuously between two finite integration limits. Droping the family index k in order to simplify the notation, the expression for the distribution function is: F ( I ) = C e x p ( - DI)
(24)
From the normalization condition, Eq. (23), we obtain: D C = (e_tgn _ e_O6 )
(25)
A continuous family is characterized by three parameters: r/, 05 and D. The parameters ~7 and 05 are calculated from: 7~ = Cm-- 0.5
(26)
05 = C M + 0.5
(27)
where C m and C M are specified and correspond to the minimum and maximum number of atoms of carbon in the molecules of a continuous family. To determine D, we assume that the average molecular weight of the continuous fraction, MW +, is related to the average number of atoms of carbon in the molecules of the continuous family, C N, through the correlation of Katz and Firoozabadi [29]: MW++ 4 C N-
(28)
14
The average number of atoms of carbon in the molecules of the continuous family is given by
CN = f f lF( I)dI
(29)
Defining
(3o)
A = 0(05 -- r/) and substituting Eqs. (24) and (25) in Eq. (29), we obtain after some rearrangement: ( C N - - r/) _ 1
05-n
e -a
(1-e -j)
(31)
In summary, given Cm, C M and M W +, Eqs. (26) and (27) are used to compute r / a n d 05. Eq. (28) is used to calculate C N, the nonlinear Eq. (31) is solved to determine A, and D is calculated in Eq.
(30). Let us now consider, as an example, the calculation of the equation of state parameter b for a continuous family of components. Using the definition of A and introducing:
z = D ( I - rl)
(32)
G.L. Rochocz et al. / Fluid Phase Equilibria 139 (1997) 137-153
144
we obtain for this integral: 4> De-D/ e--- Dd~) b( Z ) d l -
f, (e-'-;--
1
(1 - e
a)
foae =b I = - - + ~ 7 D
dz=
,
(zj)
Y'~ Wjb I;=--D +~7
j= I
(33) where the last equality comes from the numerical approximation of the integral. In this calculation, Wj and zj are the weights and points of Gaussian quadrature respectively, and nq is the number of quadrature points. The extension of these expressions to multidimensional integrals, as required for instance to calculate the parameter a in the Peng-Robinson equation of state, presents no difficulty and is not presented here. The integral in Eq. (33) has two finite limits and could have been approximated using the Gauss-Legendre quadrature. However, the integrand has a decaying exponential and Shibata et al. [28] obtained satisfactory results in the calculation of bubble and dew lines using a generalized Gauss-Laguerre quadrature between finite limits. In their paper, tables for various values of J with two and three quadrature points can be found, determined by the solution of the nonlinear system of 2nq equations that results from satisfying the 0th to the ( 2 n q - 1)th moments of the distribution. However, this system of equations is ill-conditioned, resulting in loss of numerical accuracy; a better alternative to determine Wj and zj that we use here and that does not have these deficiencies is to apply the Wheeler [30] and then the Oolub-Welsch [31] algorithms, both implemented by Press et al. [32] in their routines 'orthog' and 'gaucof'. The Wheeler procedure requires the computation of 2nq modified moments that are definite integrals whose integrand is the product of the distribution function and Legendre polynomials; general expressions for these integrals can be determined without great difficulty for the bounded exponential distribution function. We tested our implementation of these computational procedures by reproducing the tables presented by Shibata et al. [28]. We used the Peng-Robinson [22] 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. [33] and the binary interaction parameters kij used for all examples are shown in Table 1. Following the suggestion of Robinson and Peng [34], if the acentric factor of a component or pseudocomponent was smaller or equal to that of n-decane (to < 0.489), the usual expression tbr K in the equation of state was employed, otherwise the K-value of this species was calculated from: K = 0.379642 + 1.48503o9 - 0.164423w 2 + 0.016666w 3 (34) Since our main interest was to explore the qualitative features of the critical diagrams we used simple interpolations with the values obtained from the table of generalized single-carbon-number Table 1 Binary interaction parameters (kii) between the discrete species and the pseudocomponents Discrete species kii CO, 0.10~' Ethane 0.01 n-Hexane 0.00
kij 0 between all pairs of pseudocomponents. "From Robinson and Peng [34]. =
G.L Rochocz et al. / Fluid Phase Equilibria 139 (1997) 137-153
145
physical properties of Whitson [35] to calculate the critical properties and acentric factor of the pseudocomponents corresponding to the quadrature points.
4. Results We first investigated the effect of the number of quadrature points (rtq) o n the calculation of critical points of a continuous mixture containing hydrocarbons from C 7 to C45 and with MW +-- 170. In Table 2, detailed results using from 1 to 4 quadrature points are shown together with only the mixture critical properties computed when 5 to 7 quadrature points were used. The results using a single quadrature point correspond to applying the correlations for the critical properties and acentric factor to just one pseudocomponent containing an average of 12.43 atoms of carbon. Interestingly, the
Table 2 Effect of the number of quadrature points on the calculated critical point of a continuous mixture ( r / = 6.5, q5 = 45.5, M W + = 170, A = 6.52, ki.j = 0 for all pairs of pseudocomponents) T~ (K)
Pc (bar)
o)
1
x
676.1 676.1
20.29 20.29
0.4669
12.43
1.0000
621.4 843.6 704.5
24.62 12.14 30.27
0.3776 0.8241
9.800 25.20
0.8294 0.1706
592.8 766.5 916.0 706.1
27.21 15.32 8.802 30.01
0.3349 0.6370 0.9403
8.637 17.96 34.83
0.6449 0.3266 0.0285
574.9 713.0 847.6 945.0 705.5
28.81 18.21 11.99 7.779 30.18
0.3119 0.5292 0.8366 0.9936
7.996 14.39 25.70 39.71
0.5004 0.4003 0.0910 0.0083
704.4
29.87
705.9
29.91
705.7
30.01
1 quadrature point Fraction 1 Mixture
2 quadrature points Fraction I Fraction 2 Mixture
3 quadrature points Fraction 1 Fraction 2 Fraction 3 Mixture
4 quadrature points Fraction Fraction Fraction Fraction Mixture
l 2 3 4
5 quadrature points Mixture
6 quadrature points Mixture
7 quadrature points Mixture
146
G.L. Rochocz et al. / Fluid Phase Equilibria 139 (1997) 137-153
results with only two quadrature points differ in mixture critical temperature and critical pressure by only 0.2% and 0.9% compared to the results for a mixture with seven quadrature points. We now consider the calculation of critical points in a pseudobinary system containing CO 2 and the continuous fraction characterized in Table 2, in order to further investigate the effect of the number of quadrature points, nq. For a given value of nq, the mole fraction of each pseudocomponent is kept constant on a CO2-free basis, throughout the computation of the critical locus. Fig. l a - c
400
(a)
~-.=17] j_= no-- jj I n =3 Ib
350
700
300 600
2so
g
f ~
200
500
u)
E a. 150
400 100
CO2 5O 0
\ I
I
~
300
400
500
.....
v
I
600
700
300
0.0
r
T
T
T
0.2
0.4
0.6
0.8
Temperature (K)
1.0
Mole fraction of CO 2
400
a
350
~-_ _
- -
(c)
nq=l __
n q=2
nq=3 _ _
300
n q=6
,250
2OO n
150
1 O0 5O 0 0.0
;
I
I
I
0.2
0.4
0.6
0.8
1.0
Mole fraction of CO 2
Fig. 1. Effect of the number of quadrature points (nq) on the critical curve of the system CO 2 +[continuous fraction] characterized in Table 2: (a) pressure-temperature projection; (b) temperature-mole fraction projection; (c) pressure-mole traction projection.
147
G.L Rochocz et al. / Fluid Phase Equilibria 139 (1997) 137-153
shows the critical projections for this system, calculated with one, two, three and six quadrature points. As would be expected, the results with a single point of quadrature are significantly different from those with larger number of pseudocomponents and will not be considered further. In the pressure-temperature ( P - T ) projection, a type III diagram in the classification of van Konynenburg and Scott [36] is obtained in all cases, with the branch of the critical line that starts at the continuous fraction passing first through a maximum in pressure and then through a minimum before proceeding to very high pressures. The branch of the critical line that starts at pure CO 2 (T~ = 304.1 K, Pc = 73.8 35
620
(a)
nq=2
\%.
(b)
--nq=3 34
//
~ ~ --..,..\ ~ ~ . / ........ -\\
33 A
.~ 32 O,.
31
30
/ /
/i
i
600
nq=6
\ \ \ ....\ \
llS////
580
\
,,\.\
560
\ \
R 54O
t
520
n-he×ariB
29 500
\
E
\.\\.\ \ :'\
I
I
i
I
520
540
560
580
\
\
\
500 600
620
0.0
0.2
Temperature (K)
0.4
0.6
0.8
1.0
Mole fraction of n-hexane
35 __-
nq=3 nq=6
34
(C)
nq=2
33
//
/ /I" / i
/
~.~ ",X \\~
//
~.k,
32 Q..
31
/i /.i i
\
\\
30
29 0.0
i
i
i
I
0.2
0.4
0.6
0.8
1.0
Mole fraction of n-hexane
Fig. 2. Effect of the number of quadrature points (nq) on the critical curve of the system n-hexane + [continuous fraction] ( 7 = 6 . 5 , ~b=45.5, MW + = 122): (a) pressure-temperature projection; (b) temperature-mole fraction projection; (c) pressure-mole fraction projection.
148
G.L. Rochocz et al. / Fluid Phase Equilibria 139 (1997) 137 153
bar) is very short, and shrinks as the number of quadrature points increases. Therefore, this branch is visible on the scale of the plots only when a single pseudocomponent is used. The critical projections with two, three or six quadrature points are similar, even at higher pressures not shown in the plots. The temperature-mole fraction (T-x) and pressure-mole fraction ( P - x ) projections show complex behavior in approximately the region 0.878 < Xco" < 0.889 (the exact limits depend o n F/q) and we found up to four critical points at some compositions. To test the effect of rtq when the discrete species and the continuous family are not as dissimilar as in the previous example, we consider a pseudobinary system composed of n-hexane and a continuous fraction containing hydrocarbons from C 7 to C45 with MW += 122, resulting in C N = 9. Fig. 2a-c shows the critical projections and the P - T projection is a type I system in the classification of van Konynenburg and Scott [36]. The results indicate that two quadrature points (pseudocomponents) are insufficient for a satisfactory representation of this system, while the results with three and six points are similar. We calculated critical points for pseudobinary mixtures using three quadrature points. Fig. 3 shows the effect of MW + on the critical diagrams of systems in which the discrete component is ethane and the continuous fraction contains hydrocarbons from C 7 to C45. However, there is an important difference between the curves that is not easily seen in the plot. The critical curves for the MW += 108 and MW += 136 mixtures are continuous and connect the critical points of the continuous fraction and of pure ethane; therefore this system is of type I. The critical loci for the systems of the heavier mixtures, MW += 164 and MW += 220, are discontinuous and these are type V systems. In these systems, the critical branches end at critical endpoints that are very close to the critical point of pure ethane. The transition from type I to V occurs through a tricritical point [37], but we made no attempt to determine the value of MW + at which this transition occurs. In systems of ethane with linear alkanes, n-octadecane (Tc = 748 K, Pc = 12 bar) is the first substance to give origin to type V behavior [37]. The values of MW += 136 and MW += 164 that bracket the transition correspond to
2OO
/ 180
/~
~\
\
/
\,
// /,/
160
\ k • -i- ~
/
/ / 140
/
~" •~
\
.~
//'
//
== ~, 1oo
x
MW*=220
......
MW*=136
\ \
///,/ ~ - ~
'
\
'\\
"\
60
40
MW*=164
---
\
•
80 i ~
MW*=108
\,
/
/ / 120
\\
- ----
\
ethane
\x\
\
\\
x
\
\\
\
x\
\ [
20 300
400
I -500
\\ r
600
700
-
800
Temperature (K)
Fig. 3. Effect of M W + on the critical curve of the system ethane + [continuous fraction] (~/= 6.5, ~b = 45.5).
G.L. Rochocz et al. / Fluid Phase Equilibria 139 (1997) 137 153
149
250
250
(a)
I
MW*=lO8
(b)
I
MW+=122
/~-'\ /
//
'
\
~
\L
l
MWf=136
r
MVV*=15o _
200
200
/
\\\
~ " 15o
~" 150
¢B
~- loo
~- 1oo
CO2
C02 50
50
0
200
i
r
I
]
I
I
I
I
250
300
350
400
450
500
550
600
650
250
J
I
r
i
I
i
r
I
300
350
400
450
500
550
600
650
Temperature (I0
700
Temperature (K)
Fig. 4. Effect o f M W + on the critical curve o f the s y s t e m carbon dioxide + [ c o n t i n u o u s fraction] ( 7 / = 6.5, (b = 45.5); (a) M W + = 108 and M W + = 122; (b) M W + = 136 and M W + = 150.
C N = 10 and C N = 12, but the critical properties of the corresponding continuous fractions are very different f r o m n-decane and n-dodecane respectively. Interestingly, in the series of n-alkanes, n-tetradecane (T~ = 695 K, Pc = 15.6 bar) and n-eicosane ( Tc = 767 K, Pc = 11.1 bar) have critical
250
250
(a)
/b)
Ca=25 ]
--
200
200
\
~" 150
~" 150 .o
#. 100
\
CO~ 50
0
200
i
I
250
300
350
\
CO2 50
I
I
I
~
I
400
450
500
550
600
Temperature (K)
\
lOO
650
200
250
300
350
400
450
500
550
600
650
Temperature (K)
Fig. 5. Effect o f the m a x i m u m n u m b e r o f c a r b o n s (C M) on the critical c u r v e o f the s y s t e m carbon dioxide + [ c o n t i n u o u s fraction] ( M W + = 122, -r/= 6.5); (a) C M = 25 (q5 = 25.5) and (b) C M = 20 ( ~ = 20.5).
150
G.L. Rochocz et al. / Fluid Phase Equilibria 139 (1997) 137 153
temperatures close to those of the continuous fractions with C N = l0 and C N = 12, but very different critical pressures. Fig. 4a,b shows the effect of MW + on the critical diagrams of systems in which the discrete component is CO 2 and the continuous fraction contains hydrocarbons from C v to C4.s. For MW += 108, there is a continuous vapor-liquid critical line between the critical points of CO 2 and of the continuous fraction, and there is a liquid-liquid critical line that proceeds to high pressures at approximately 250 K. The shape of this P - T projection is characteristic of type II systems. For MW += 122, the vapor-liquid critical line is interrupted at critical endpoints, and a critical diagram of type IV results. At MW += 136 a diagram of type IV is also calculated to occur, but for MW += 150 a diagram of type II! is obtained. The transition from type IV to type III critical behaviors occurs at a double critical endpoint [37] but we did not attempt to determine its thermodynamic coordinates and the value of MW + at which this occurs. Type IV behavior appeared here in a relatively narrow range of MW+; interestingly, CO 2 + n-tridecane is one of the few binary hydrocarbon systems that show type IV behavior, as such critical behavior also occurs only in a narrow range of heavy component molecular weights for real mixtures. Because the occurrence of type IV diagrams in the calculations are so sensitive to the parameters of the components in the mixture, it is interesting to use these parameters to investigate the effect of the minimum and maximum number of carbon assumed that are assumed to exist in a given distribution. The system CO 2 + [continuous fraction] with MW += 122 and hydrocarbons from C v to C45 (Fig. 4a) was selected as the standard for comparison. In Fig. 5a and b, we see the effect of changing the heaviest hydrocarbon from C45 to C2s and C2o. Fig. 5a (C M = 25) still retains the characteristics of a type IV diagram, while Fig. 5b (C M = 20) is of type II. In Fig. 6a and b, the effect of changing the lightest hydrocarbon from C 7 to C7.05 and C 7.J is shown. Conversely to what happens when the heavy end is changed, the topology of these critical diagrams is very sensitive to changes in the minimum carbon number, C m, employed to characterize the continuous fractions. 250
250
(a)
~
- -
cm=7°5
(b)
1
200
200
150
~" 15o
\
~- loo
/
\
C02 \,
5O
200
Cr.=7.1
lOO
C02
0
[--
I
I
[
I
I
I
250
300
350
400
450
500
Temperature (K)
50
\,
J
550
600
0
650
200
I
I
I
1
I
250
300
350
400
450
~
500
f - -
550
600
650
Temperature (K)
Fig. 6. Effect o f the m i n i m u m n u m b e r o f c a r b o n s ( C , , ) on the critical curve o f the s y s t e m carbon dioxide + [continuous fraction] ( M W + = 122, 4~ = 45.5); (a) C m = 7.05 (rl = 6.55) a n d (b) C m = 7.10 ( ~ / = 6.60).
G.L. Rochocz et aL / Fluid Phase Equilibria 139 (1997) 137-153
151
5. Conclusions The first objective of this paper was to show how the formulation of critical points of continuous mixtures fits into the general framework of integral operators and functional analysis. The formalism obtained in this context is readily generalized to the formulation of the conditions for critical transitions of higher order such as tricritical points, tetracritical points, etc. The second objective was to calculate critical points in systems containing a discrete component and a continuous fraction in order to study the effect of the distribution parameters on the critical projections. For this we used a generalized form of the Gauss-Laguerre quadrature to generate pseudocomponents and then used a modified version of the Hicks and Young [25] algorithm to perform the critical point calculations. The Peng-Robinson equation of state was used with a simple procedure to estimate the properties of the pseudocomponents. In our examples, three quadrature points were sufficient to provide acceptable representations of the systems studied. Despite the simple characterization procedure for the continuous fractions, by varying the light component and the average molecular weight of fraction, we could find five types of the van Konynenburg and Scott classification of critical projections in the P - T plane that cubic equations of state with one-fluid mixing rules are known to exhibit for binary mixtures. Using a type IV diagram as standard for comparison, the influence on the calculations of the mimimum and maximum number of carbons in the distribution was studied. Large variations of the maximum number of carbons could be made before the type of critical phase diagram changed. However, small variations in the minimum carbon number were sufficient to change the type of critical phase diagram obtained.
6. Notation a
A b c
C Cm CM CN D F
h(t,s,u) I Kij
k(s,t) MW + n n c r/d F/f
EOS attractive parameter Helmholtz free energy EOS covolume parameter Cubic form of the Helmholtz free energy expansion Auxiliary parameter of the distribution function Minimum number of atoms of carbon in the molecules of a continuous family Maximum number of atoms of carbon in the molecules of a continuous family Average number of carbons per molecule in a continuous mixture Parameter of the distribution function Intensive distribution function Third order kernel of integral operator T Characterization parameter of the continuous fraction Binary interaction parameter Kernel of the integral operator T Average molecular weight of a continuous fraction Total number of moles Number of components Number of discrete components Number of continuous families of components
152
G.L. Rochoc: et al. / Fluid Phase Equilibria 13911997) 137-153
Number of moles of component I Number of quadrature points Pressure Quadratic form of the Helmholtz free energy expansion Hessian matrix of the Helmholtz free energy expansion Temperature Total volume Weights of Gaussian quadrature Molar fraction of component i Points of Gaussian quadrature
n I
nq
P q
Q T V x i Zj
Greek letters
Parameter of the distribution function Parameter of the Peng-Robinson equation of state Eigenvalue Integral operator Lower limit of the distribution index Chemical potential Upper limit of the distribution index Acentric factor
A K
A T t-t
4, o)
Subscripts C
Critical property
i,j,k
Refer to components i,j,k respectively
Acknowledgements G.L.R. and M.C. acknowledge the comments and suggestions of Mfircio L.L. Paredes in the early stages of this work. The preparation of this manuscript was supported, in part, by Grant No. DE-FG02-85ER13436 from the U.S. Department of Energy, and Grant No. CTS-9521406 from the U.S. National Science Foundation, both to the University of Delaware. M.C. acknowledges the financial support of the Brazilian Ministry of Education (CAPES).
References [1] [2] [3] [4] [5] [6]
K.S. Pedersen, P. Thomassen, Aa. Fredenslund, Fluid Phase Equilibria 14 (1983) 21)9. K.S. Pedersen, P. Thomassen, Aa. Fredenslund, Ind. Eng. Chem. Process Des. Dev. 23 (1984) 163. K.S.P. Pedersen, P. Thomassen, Aa. Fredenslund, Ind. Eng. Chem. Process Des. Dev. 24 (1985) 948. J.R. Bownman, Ind. Eng. Chem. 41 (1949) 2004. W.C. Edmister, D,H. Buchanan, Chem. Eng. Prog. Symp. Ser. No. 6 (1953) 69. M.T. R~itzch, H. Kehlen, Fluid Phase Equilibria 14 (1983) 225.
G.L. Rochocz et al. / Fluid Phase Equilibria 139 (1997) 137-153
[7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37]
153
H. Kehlen, M.T. R~itzsch, Z. Phys. Chem. (Leipzig) 265 (1984) 1049. H. Kehlen, M.T. Rfitzsch, J. Bergmann, AIChE J. 31 (1985) 1136. R.L. Cotterman, J.M. Prausnitz, Revue de L'Institut Fran~ais du P6trole 45 (1990) 633. G.L. Rochocz, MSc Thesis, COPPE, Federal University of Rio de Janeiro, Brazil, 1990. T.J.P. Halpin, N. Quirke, SPE Reservoir Eng. 5 (1990) 617. H.W. Haynes Jr., M.A. Mathews, Ind. Eng. Chem. Res. 30 (1991) 1911. D.S. Chakravarty, M.A. Mathews, Ind. Eng. Chem. Res. 33 (1994) 1962. K.D. Luks, E.A. Turek, T.K. Kragas, Ind. Eng. Chem. Res. 32 (1993) 1767. J. Kramarz, A. Wyczesany, Chem. Eng. Sci. 48 (1993) 1665. C. Lira-Galeana, A. Firoozabadi, J.M. Prausnitz, Fluid Phase Equilibria 102 (1994) 143. H. Kehlen, M.T. Rfitzsch, J. Bergmann, J. Macromol. Sci. Chem. A24 (1987) 1. D. Browarzik, H. Kehlen, M.T. R~itzsch, J. Macromol. Sci. Chem. A27 (1990) 549. Y. Hu, X. Xing, D.T. Wu, J.M. Prausnitz, Fluid Phase Equilibria 104 (1995) 229. D. Browarzik, H. Kehlen, Fluid Phase Equilibria 123 (1996) 17. R.A. Heidemann, A.M. Khalil, AIChE J. 26 (1980) 769. D.-Y. Peng, D.B. Robinson, Ind. Eng. Chem. Fundam. 15 (1976) 59. M. Castier, S.I. Sandier, Critical point calculations with the Wong-Sandler mixing rule: I. Calculations with the van der Waals Equation of State, Chem. Eng. Sci. 52 (1997) 3393. M. Castier, S.I. Sandler, Critical point calculations with the Wong-Sandler mixing rule: II. Calculations with a modified Peng-Robinson Equation of State, Chem. Eng. Sci., in press. C.P. Hicks, C.L. Young, J. Chem. Soc., Faraday Trans. II 73 (1977) 597. B. Epstein, Linear Functional Analysis: An Introduction to Lebesque Integration and Infinite Dimensional Problems, Saunders, Philadelphia, PA, 1970. M.L. Michelsen, Fluid Phase Equilibria 8 (1982) 1. S.K. Shibata, S.I. Sandier, R.A. Behrens, Chem. Eng. Sci. 42 (1987) 1977. D.L. Katz, A. Firoozabadi, J. Pet. Technol. 20 (1978) 1649. J.C. Wheeler, Rocky Mountain J. Math. 4 (1974) 287. G.H. Golub, J.H. Welsch, Math. Comput., 23 (1969) 221-230 and A I - A I 0 . W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in FORTRAN, 2nd edn., Cambridge Univ. Press, Cambridge, UK, 1992. R.C. Reid, J.M. Prausnitz, B.E. Poling, The Properties of Gases and Liquids, 4th edn., McGraw-Hill, New York, 1987. D.B. Robinson, D.-Y. Peng, Research Report RR-28, Project 756, Gas Processors Association, 1978. C.H. Whitson, Soc. Petrol. Eng. J. (August, 1983) 683. P.H. van Konynenburg, R.L. Scott, Philos. Trans. R. Soc. London, Ser. A 298 (1980) 495. A. van Pelt, PhD Thesis, The Technical University of Delft, Holland, 1992.
Fluid Phase Equilibria 139 (1997) 137-153
Critical point calculations for semi-continuous mixtures G e r a l d o L. R o c h o c z a,b, M a r c e l o C a s t i e r b,~, S t a n l e y I. S a n d l e r c,, Chemtech Sert,i~'os de Engenharia, Rua Nilo Pe~'anha 51, sala 2106, Rio de Janeiro, RJ, Brazil h Escola de Quhnica, Unit~ersidade Federal do Rio de Janeiro, C.P. 68542, Rio de Janeiro, RJ, 21949-900, Brazil ~ Center/or Molecular and Engineering Thermodynamics, Department of" Chemical Engineering Unit,ersi O, of Delaware Newark, DE 19716, USA
Abstract A typical problem in the oil and gas industry is the calculation of high pressure phase equilibria in systems containing a very large number of similar components. Treating the composition of these mixtures as a continuous distribution, instead of trying to characterize each individual component, has proven to be an elegant and efficient approach to solve this type of problem. In this paper, we rigorously formulate the critical point conditions for continuous and semi-continuous mixtures, but as is often the case with such formulations, the complexity of the expressions prevents an analytical solution of the equilibrium problem, and discretization is required in order to use realistic thermodynamic models. Once in discrete form, the critical point equations were solved using a modified form of the Hicks and Young algorithm and the thermodynamic stability of the solutions was guaranteed by the calculation of higher terms in the Helmholtz free energy expansion and by global phase stability analyses. The resulting procedure was employed to perform an investigation of the influence of the distribution parameters on the types of critical phase diagrams that can be obtained from the Peng-Robinson equation of state. © 1997 Elsevier Science B.V. Keywords: Critical state; Vapor-liquid equilibria; Method of calculation; Equation of state; Continuous thermodynamics;
Mixture
1. I n t r o d u c t i o n The characterization o f mixtures, such as oil reservoir fluids, coal polymer solutions, for physical property and equilibrium calculations very large number of components. Instead of trying to identify each these mixtures as an infinite set of components characterized by one
derived liquids and polydisperse may be difficult because of their individual component, modeling or more continuous variables, in
* Corresponding author. On leave at the Department of Chemical Engineering, University of Delaware, Newark, DE, USA. 0378-3812/97/$17.00 © 1997 Elsevier Science B.V. All rights reserved. PII S 0 3 7 8 - 3 8 1 2 ( 9 7 ) 0 0 2 1 2 - 4
138
G.L. Rochocz et al. / Fluid Phase Equilibria 139 (1997) 137-153
an approach known as the thermodynamics of continuous mixtures, has proven to be an elegant and efficient approach to solve this type of problem. A hybrid description, known as semi-continuous mixture thermodynamics, is to characterize some of the components individually, usually the lighter components, and to group the remaining species into one or more continuous families. In typical applications of the thermodynamics of continuous mixtures, summations over the mixture components are recast into integrals that often need to be numerically evaluated. The use of efficient integration methods generally results in computational implementations that compare favorably in terms of accuracy and speed with the traditional approach of lumping species into arbitrarily chosen pseudocomponents (Pedersen et al. [1-3]). The idea of using continuous distributions to characterize mixtures dates back to the work of Bowman [4] and Edmister and Buchanan [5], but the rigorous formulation of the thermodynamics of continuous mixtures based on functional analysis, found in the work of R~itzsch and Kehlen [6,7] and Kehlen et al. [8] renewed interest in this field. Work in this area has focused on the development of algorithms to solve the new phase equilibrium functional equations, the characterization of the continuous mixtures and the use of thermodynamic models in continuous form. Much of the work done upto 1990 was reviewed by Cotterman and Prausnitz [9] and Rochocz [10], and new applications of the continuous mixture approach have been published in recent years. Halpin and Quirke [11] used distribution functions for continuous mixtures that are written as an expansion in orthonormal functions, which allows the mass balance to be preserved in flash calculations. Haynes and Mathews [ 12] used the experimental TBP curve directly as distribution function to calculate VLE, and this approach was later incorporated into a reservoir simulator [13]. Luks et al. [14] studied the effect of varying the choice of pseudocomponents in each stage of a multiple-contact separation process. Kramarz and Wyczesany [15] formulated and solved the chemical equilibria of continuous mixtures. Lira-Galeana et al. [16] employed a continuous mixture approach to model the segregation of heavy components due to gravity in oil reservoirs. Despite the large amount of work done on the calculation of phase equilibria in continuous and semi-continuous systems, comparatively few papers have dealt with the determination of phase stability in these systems. The stability of polydisperse polymer systems was studied by Kehlen et al. [17] and Browarzik et al. [18], using segment fractions as independent variables in the formulation. Hu et al. [19] followed this approach to compute spinodal curves, critical points and other phenomena in polymer solutions. More recently, Browarzik and Kehlen [20] developed the theoretical framework for a continuous mixture stability theory that was applied to the van der Waals equation of state (EOS), using one-fluid mixing rules with the binary interaction parameters kij s e t equal to zero. In the van der Waals EOS, they further assumed that ~/a and b of each fraction were linear functions of a single characterization parameter M, chosen to be the number of carbon atoms in the components of the n-alkane mixture used in their procedure. Under these hypotheses, they determined spinodal curves and critical points for these mixtures. In this paper, we are first interested in formally extending the criticality conditions of Heidemann and Khalil [21] to continuous mixtures, using elements of functional analysis. We then investigate the effect of the distribution parameters on the types of critical loci that are obtained. We consider pseudobinary mixtures containing a light substance and a continuous family of components modeled by the Peng-Robinson [22] equation of state with the one-fluid mixing rules and non-zero binary interaction parameters. Under these hypotheses, the integrals that appear in the semi-continuous formulation do not have analytical solutions and we resort to numerical methods to obtain their
G.L. Rochocz et al. / Fluid Phase Equilibria 139 (1997) 137-153
139
values. The criticality conditions in the discretized variables are solved using a modified version [23,24] of the Hicks and Young algorithm [25].
2. Critical points in continuous mixtures The Helmholtz free energy of a continuous mixture is a state functional and can be written [8] as:
A=A(T,V,n,F(I))
(1)
where n is the total number of moles, F ( I ) is an intensive distribution function which in general is a function of several characterization variables, but for simplicity we consider a single characterization variable I. The development that follows is also applicable to semi-continuous mixtures, because the discrete components can be represented in the formulation by delta functions in the distribution function. From the normalization condition we have that:
fF(I)dI=
1
(2)
Following Heidemann and Khalil [21] and using the Kehlen et al. [8] notation for state functional partial derivatives, the Taylor expansion of the continuous mixture Helmholtz free energy function A at constant temperature and total volume is:
A-AO-fltx(I)AF(I)dI
= ~" f/f/+
--
OntOnt+ r,v
03A OnlOnl+Onl.
AF(I)dlAF(I+)dI +
AF(I)dlAF(I+)dl + AF(I*)dI* +...
(3)
T,V
where the chemical potential in a continuous mixture corresponds to the first partial derivative of A with respect to the distribution function:
/x(I)=
(OA) On1 r,v --
(4)
The system is locally stable if the first non-vanishing term on the right hand side of Eq. (3) is of even order and positive, and for any arbitrary perturbations AF(I)dI and AF(I+)dI +, this most commonly occurs when:
flf+ ( 02A ) AF(I)dlAF(I+)dI+>O OnzOnt+ r,v
(5)
At a critical point, both the second and third order terms of the Helmholtz free energy expansion
140
G.L. Rochocz et al. / Fluid Phase Equilibria 139 (1997) 137-153
must vanish, i.e.:
..z ) On1Ont+ r,v flfl fl ( 03A ) AF(1)dIAF(I+)dl+AF(I*)dI * = 0 ' * OntOnI+ On1" r,v
(7)
The corresponding expressions for the quadratic (q) and cubic (c) forms in a mixture containing only discrete components are:
q =
Onjc3n i j= I i=l
c=
iii( k=l j=l
i=l
T.V. A l ' l i A n j n,.,.j
O'A cgnkcgnjOni
: 0
"r,v. AniAnj An~ =
(8)
0
(9)
n,,,,.~
for any perturbations An i, An i and An k. The zero of the quadratic form in the discrete case can be detected by a zero value of the determinant of the matrix Q whose elements are:
(..z)
QiJ = OnjOni r,v,
(10)
tlt ~ t,i
or by obtaining a zero value for at least one of the eigenvalues of this matrix. For the continuous case, if a discretization is carried out, the discrete formulae for criticality and the numerical calculation procedures associated with them can be used, as we do in the applications which follow. It is, however, interesting to comment on the theoretical framework that formally supports the extension of the analysis for systems with infinitely many components. This framework comes from the fact that many results of linear algebra can be extended to spaces of infinite dimension, where matrix algebra is replaced by functional analysis of linear operators. Here, we will only mention the most relevant aspects of this extension and will refrain from giving proofs or providing the complete development of the formalism. For this we refer to the discussion of spectral theory in infinite dimensional Hilbert spaces found in Epstein [26]. This extension starts by considering the integral operator T:C[a,b] ~ C[a,b], where C[a,b] is the normed space of all continuously differentiable functions on [a,b], which is defined by Tf= g so that: (Tf)(s) =
g(s) = f bk(s,t)f(t)dt
(11)
where k(s,t) is a continuous function defined on [a,b] × [a,b], called the kernel of T. It can be shown [26] that the integral operator is a compact linear operator. An inner product can be defined in C[a,b] as:
(12)
G.L. Rochocz et al. / Fluid Phase Equilibria 139 (1997) 137-153
141
which defines a Hilbert space with the norm Hf]] = ( ( f , f ) , and this allows the definition of a norm for the integral operator. With the definition of the inner product, we can write the quadratic form for the integral operator T as:
( f ,Tf)= £b~bk( t,s)f( t)f(
,)dtds
(13)
The second functional derivative of the Helmholtz free energy ({(02A)}/{(0ni0ns+)})r.v = A"(I + I ÷) is a symmetric kernel function for the integral operator that appears in the quadratic form; in such a case the integral operator defined by this function is said to be hermitian. Epstein [26] demonstrates that if all the eigenvalues of a compact linear and hermitian operator are positive, the quadratic form ( f, Tf) is positive-definite. Higher-order terms can also be written using this formalism. For instance, we can introduce a third order kernel h(t,s,u) that defines a new integral operator F:C[a,b] ~ C[a,b] × C[a,b] as:
( F f )( t,s) = g( t,s) = ~bh( t,s,u)f( u)du
(14)
The new integral operator Ff(u) provides a function g(t,s) that can be considered to be a second order kernel g(t,s) for an integral operator T. It is therefore possible to define a cubic form: rb rb rb
(f, Vf,f}=(f,((Vf),f}}=((( Vf),f},f}= J. ],, Ja h(t,s,u)f(t)f(s)f(u)dtdsdu
(15)
At the stability limit, the quadratic form vanishes which means that at least one eigenvalue Ai is equal to zero, i.e., Ai=O
(16)
When the quadratic form becomes positive-semidefinite, it is possible to find a perturbation distribution function AF(I) so that:
flA"( I,I +) AF( I +)dI +=
0
(17)
and
flAFZ( l)dl
1
(18)
given by the normalized eigenfunction corresponding to the vanishing eigenvalue. Eqs. (17) and (18) are the continuous analogs of
QAn=O
(19)
and /7 c
E
= 1
(20)
i=1
that are found in the work of Heidemann and Khalil [21]. In fact, Eqs. (16), (17) and (9) represent a rigorous extension to continuous mixtures of the Heidemann and Khalil [21] formulation for the calculation of critical points.
G.L. Rochocz et al. / Fluid Phase Equilibria 139 (1997) 137-153
142
Although theoretically rigorous, the continuous criterion for critical points is cumbersome for numerical calculations. For discrete mixtures, one vanishing eigenvalue is detected by a zero determinant. In the continuous case, the integral operator can have infinite eigenvalues and it is very difficult to conclude when one of them is equal to zero. Furthermore, the integrals often do not have analytical solutions when realistic thermodynamic models are used, and must be numerically approximated using a method such as Gaussian quadrature. From a practical point of view, this corresponds to selecting a group a pseudocomponents and their amount in a given mixture, but based on a sound mathematical criterion. Once this is done, programs usually used for the calculation of phase equilibria, physical properties or, in our particular application, critical points, with mixtures of discrete components can be used. To calculate critical points, we used a modified version of the Hicks and Young [25] algorithm; additional details about the implementation of this method can be found elsewhere [23,24]. To test the stability of the calculated critical points, we numerically evaluated the fourth-order term of the A expansion using a four-point difference formula and the global stability test [27]. Only critical points that passed the stability tests are reported here. However, the global stability test was only used to verify the possibility of additional fluid phases; the possibility of forming solid phases was not tested in this work.
3. Modeling of semi-continuous mixtures The formulation of critical points that has been presented is not limited to any particular distribution function and discrete components can be included in the formulation as delta functions in the distribution function. However, it is usual to distinguish explicitly between discrete species and continuous families of components and we will follow this approach in the rest of this paper. As simple examples, consider the fact that the mole fraction of all species should sum to unity
Y'~xi+~xk=
~_,xi+ Y'~ x k
i=1
i=1
=1
k=l
Fk(Ik)dI k = 1
(21)
k
and that the expression to determine the parameter b of the Peng-Robinson equation of state when the one-fluid mixing rule is used is
b = Y'~ xib i + ~ i=1
xk
(I~.)bk(lk)dl k
(22)
k=l
where the indexes i and k denote a discrete component and a continuous family respectively, x denotes mole fractions, and rl~ and ~bk are the lower and upper limits for the characterization variable I~ of family k. Again, for simplicity, we assume that each family is characterized by a single variable Ik; also the distribution functions are normalized
f,7~kFk( Ik)dlk= l
(23)
Several procedures exist to characterize oil fluids and the discussion that follows can be adapted to accommodate different characterization procedures, but we will follow Shibata et al. [28] and choose a bounded exponential function in which the characterization variable I is the number of carbon
143
G.L. Rochocz et al. / Fluid Phase Equilibria 139 (1997) 137-153
atoms, assumed to vary continuously between two finite integration limits. Droping the family index k in order to simplify the notation, the expression for the distribution function is: F ( I ) = C e x p ( - DI)
(24)
From the normalization condition, Eq. (23), we obtain: D C = (e_tgn _ e_O6 )
(25)
A continuous family is characterized by three parameters: r/, 05 and D. The parameters ~7 and 05 are calculated from: 7~ = Cm-- 0.5
(26)
05 = C M + 0.5
(27)
where C m and C M are specified and correspond to the minimum and maximum number of atoms of carbon in the molecules of a continuous family. To determine D, we assume that the average molecular weight of the continuous fraction, MW +, is related to the average number of atoms of carbon in the molecules of the continuous family, C N, through the correlation of Katz and Firoozabadi [29]: MW++ 4 C N-
(28)
14
The average number of atoms of carbon in the molecules of the continuous family is given by
CN = f f lF( I)dI
(29)
Defining
(3o)
A = 0(05 -- r/) and substituting Eqs. (24) and (25) in Eq. (29), we obtain after some rearrangement: ( C N - - r/) _ 1
05-n
e -a
(1-e -j)
(31)
In summary, given Cm, C M and M W +, Eqs. (26) and (27) are used to compute r / a n d 05. Eq. (28) is used to calculate C N, the nonlinear Eq. (31) is solved to determine A, and D is calculated in Eq.
(30). Let us now consider, as an example, the calculation of the equation of state parameter b for a continuous family of components. Using the definition of A and introducing:
z = D ( I - rl)
(32)
G.L. Rochocz et al. / Fluid Phase Equilibria 139 (1997) 137-153
144
we obtain for this integral: 4> De-D/ e--- Dd~) b( Z ) d l -
f, (e-'-;--
1
(1 - e
a)
foae =b I = - - + ~ 7 D
dz=
,
(zj)
Y'~ Wjb I;=--D +~7
j= I
(33) where the last equality comes from the numerical approximation of the integral. In this calculation, Wj and zj are the weights and points of Gaussian quadrature respectively, and nq is the number of quadrature points. The extension of these expressions to multidimensional integrals, as required for instance to calculate the parameter a in the Peng-Robinson equation of state, presents no difficulty and is not presented here. The integral in Eq. (33) has two finite limits and could have been approximated using the Gauss-Legendre quadrature. However, the integrand has a decaying exponential and Shibata et al. [28] obtained satisfactory results in the calculation of bubble and dew lines using a generalized Gauss-Laguerre quadrature between finite limits. In their paper, tables for various values of J with two and three quadrature points can be found, determined by the solution of the nonlinear system of 2nq equations that results from satisfying the 0th to the ( 2 n q - 1)th moments of the distribution. However, this system of equations is ill-conditioned, resulting in loss of numerical accuracy; a better alternative to determine Wj and zj that we use here and that does not have these deficiencies is to apply the Wheeler [30] and then the Oolub-Welsch [31] algorithms, both implemented by Press et al. [32] in their routines 'orthog' and 'gaucof'. The Wheeler procedure requires the computation of 2nq modified moments that are definite integrals whose integrand is the product of the distribution function and Legendre polynomials; general expressions for these integrals can be determined without great difficulty for the bounded exponential distribution function. We tested our implementation of these computational procedures by reproducing the tables presented by Shibata et al. [28]. We used the Peng-Robinson [22] 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. [33] and the binary interaction parameters kij used for all examples are shown in Table 1. Following the suggestion of Robinson and Peng [34], if the acentric factor of a component or pseudocomponent was smaller or equal to that of n-decane (to < 0.489), the usual expression tbr K in the equation of state was employed, otherwise the K-value of this species was calculated from: K = 0.379642 + 1.48503o9 - 0.164423w 2 + 0.016666w 3 (34) Since our main interest was to explore the qualitative features of the critical diagrams we used simple interpolations with the values obtained from the table of generalized single-carbon-number Table 1 Binary interaction parameters (kii) between the discrete species and the pseudocomponents Discrete species kii CO, 0.10~' Ethane 0.01 n-Hexane 0.00
kij 0 between all pairs of pseudocomponents. "From Robinson and Peng [34]. =
G.L Rochocz et al. / Fluid Phase Equilibria 139 (1997) 137-153
145
physical properties of Whitson [35] to calculate the critical properties and acentric factor of the pseudocomponents corresponding to the quadrature points.
4. Results We first investigated the effect of the number of quadrature points (rtq) o n the calculation of critical points of a continuous mixture containing hydrocarbons from C 7 to C45 and with MW +-- 170. In Table 2, detailed results using from 1 to 4 quadrature points are shown together with only the mixture critical properties computed when 5 to 7 quadrature points were used. The results using a single quadrature point correspond to applying the correlations for the critical properties and acentric factor to just one pseudocomponent containing an average of 12.43 atoms of carbon. Interestingly, the
Table 2 Effect of the number of quadrature points on the calculated critical point of a continuous mixture ( r / = 6.5, q5 = 45.5, M W + = 170, A = 6.52, ki.j = 0 for all pairs of pseudocomponents) T~ (K)
Pc (bar)
o)
1
x
676.1 676.1
20.29 20.29
0.4669
12.43
1.0000
621.4 843.6 704.5
24.62 12.14 30.27
0.3776 0.8241
9.800 25.20
0.8294 0.1706
592.8 766.5 916.0 706.1
27.21 15.32 8.802 30.01
0.3349 0.6370 0.9403
8.637 17.96 34.83
0.6449 0.3266 0.0285
574.9 713.0 847.6 945.0 705.5
28.81 18.21 11.99 7.779 30.18
0.3119 0.5292 0.8366 0.9936
7.996 14.39 25.70 39.71
0.5004 0.4003 0.0910 0.0083
704.4
29.87
705.9
29.91
705.7
30.01
1 quadrature point Fraction 1 Mixture
2 quadrature points Fraction I Fraction 2 Mixture
3 quadrature points Fraction 1 Fraction 2 Fraction 3 Mixture
4 quadrature points Fraction Fraction Fraction Fraction Mixture
l 2 3 4
5 quadrature points Mixture
6 quadrature points Mixture
7 quadrature points Mixture
146
G.L. Rochocz et al. / Fluid Phase Equilibria 139 (1997) 137-153
results with only two quadrature points differ in mixture critical temperature and critical pressure by only 0.2% and 0.9% compared to the results for a mixture with seven quadrature points. We now consider the calculation of critical points in a pseudobinary system containing CO 2 and the continuous fraction characterized in Table 2, in order to further investigate the effect of the number of quadrature points, nq. For a given value of nq, the mole fraction of each pseudocomponent is kept constant on a CO2-free basis, throughout the computation of the critical locus. Fig. l a - c
400
(a)
~-.=17] j_= no-- jj I n =3 Ib
350
700
300 600
2so
g
f ~
200
500
u)
E a. 150
400 100
CO2 5O 0
\ I
I
~
300
400
500
.....
v
I
600
700
300
0.0
r
T
T
T
0.2
0.4
0.6
0.8
Temperature (K)
1.0
Mole fraction of CO 2
400
a
350
~-_ _
- -
(c)
nq=l __
n q=2
nq=3 _ _
300
n q=6
,250
2OO n
150
1 O0 5O 0 0.0
;
I
I
I
0.2
0.4
0.6
0.8
1.0
Mole fraction of CO 2
Fig. 1. Effect of the number of quadrature points (nq) on the critical curve of the system CO 2 +[continuous fraction] characterized in Table 2: (a) pressure-temperature projection; (b) temperature-mole fraction projection; (c) pressure-mole traction projection.
147
G.L Rochocz et al. / Fluid Phase Equilibria 139 (1997) 137-153
shows the critical projections for this system, calculated with one, two, three and six quadrature points. As would be expected, the results with a single point of quadrature are significantly different from those with larger number of pseudocomponents and will not be considered further. In the pressure-temperature ( P - T ) projection, a type III diagram in the classification of van Konynenburg and Scott [36] is obtained in all cases, with the branch of the critical line that starts at the continuous fraction passing first through a maximum in pressure and then through a minimum before proceeding to very high pressures. The branch of the critical line that starts at pure CO 2 (T~ = 304.1 K, Pc = 73.8 35
620
(a)
nq=2
\%.
(b)
--nq=3 34
//
~ ~ --..,..\ ~ ~ . / ........ -\\
33 A
.~ 32 O,.
31
30
/ /
/i
i
600
nq=6
\ \ \ ....\ \
llS////
580
\
,,\.\
560
\ \
R 54O
t
520
n-he×ariB
29 500
\
E
\.\\.\ \ :'\
I
I
i
I
520
540
560
580
\
\
\
500 600
620
0.0
0.2
Temperature (K)
0.4
0.6
0.8
1.0
Mole fraction of n-hexane
35 __-
nq=3 nq=6
34
(C)
nq=2
33
//
/ /I" / i
/
~.~ ",X \\~
//
~.k,
32 Q..
31
/i /.i i
\
\\
30
29 0.0
i
i
i
I
0.2
0.4
0.6
0.8
1.0
Mole fraction of n-hexane
Fig. 2. Effect of the number of quadrature points (nq) on the critical curve of the system n-hexane + [continuous fraction] ( 7 = 6 . 5 , ~b=45.5, MW + = 122): (a) pressure-temperature projection; (b) temperature-mole fraction projection; (c) pressure-mole fraction projection.
148
G.L. Rochocz et al. / Fluid Phase Equilibria 139 (1997) 137 153
bar) is very short, and shrinks as the number of quadrature points increases. Therefore, this branch is visible on the scale of the plots only when a single pseudocomponent is used. The critical projections with two, three or six quadrature points are similar, even at higher pressures not shown in the plots. The temperature-mole fraction (T-x) and pressure-mole fraction ( P - x ) projections show complex behavior in approximately the region 0.878 < Xco" < 0.889 (the exact limits depend o n F/q) and we found up to four critical points at some compositions. To test the effect of rtq when the discrete species and the continuous family are not as dissimilar as in the previous example, we consider a pseudobinary system composed of n-hexane and a continuous fraction containing hydrocarbons from C 7 to C45 with MW += 122, resulting in C N = 9. Fig. 2a-c shows the critical projections and the P - T projection is a type I system in the classification of van Konynenburg and Scott [36]. The results indicate that two quadrature points (pseudocomponents) are insufficient for a satisfactory representation of this system, while the results with three and six points are similar. We calculated critical points for pseudobinary mixtures using three quadrature points. Fig. 3 shows the effect of MW + on the critical diagrams of systems in which the discrete component is ethane and the continuous fraction contains hydrocarbons from C 7 to C45. However, there is an important difference between the curves that is not easily seen in the plot. The critical curves for the MW += 108 and MW += 136 mixtures are continuous and connect the critical points of the continuous fraction and of pure ethane; therefore this system is of type I. The critical loci for the systems of the heavier mixtures, MW += 164 and MW += 220, are discontinuous and these are type V systems. In these systems, the critical branches end at critical endpoints that are very close to the critical point of pure ethane. The transition from type I to V occurs through a tricritical point [37], but we made no attempt to determine the value of MW + at which this transition occurs. In systems of ethane with linear alkanes, n-octadecane (Tc = 748 K, Pc = 12 bar) is the first substance to give origin to type V behavior [37]. The values of MW += 136 and MW += 164 that bracket the transition correspond to
2OO
/ 180
/~
~\
\
/
\,
// /,/
160
\ k • -i- ~
/
/ / 140
/
~" •~
\
.~
//'
//
== ~, 1oo
x
MW*=220
......
MW*=136
\ \
///,/ ~ - ~
'
\
'\\
"\
60
40
MW*=164
---
\
•
80 i ~
MW*=108
\,
/
/ / 120
\\
- ----
\
ethane
\x\
\
\\
x
\
\\
\
x\
\ [
20 300
400
I -500
\\ r
600
700
-
800
Temperature (K)
Fig. 3. Effect of M W + on the critical curve of the system ethane + [continuous fraction] (~/= 6.5, ~b = 45.5).
G.L. Rochocz et al. / Fluid Phase Equilibria 139 (1997) 137 153
149
250
250
(a)
I
MW*=lO8
(b)
I
MW+=122
/~-'\ /
//
'
\
~
\L
l
MWf=136
r
MVV*=15o _
200
200
/
\\\
~ " 15o
~" 150
¢B
~- loo
~- 1oo
CO2
C02 50
50
0
200
i
r
I
]
I
I
I
I
250
300
350
400
450
500
550
600
650
250
J
I
r
i
I
i
r
I
300
350
400
450
500
550
600
650
Temperature (I0
700
Temperature (K)
Fig. 4. Effect o f M W + on the critical curve o f the s y s t e m carbon dioxide + [ c o n t i n u o u s fraction] ( 7 / = 6.5, (b = 45.5); (a) M W + = 108 and M W + = 122; (b) M W + = 136 and M W + = 150.
C N = 10 and C N = 12, but the critical properties of the corresponding continuous fractions are very different f r o m n-decane and n-dodecane respectively. Interestingly, in the series of n-alkanes, n-tetradecane (T~ = 695 K, Pc = 15.6 bar) and n-eicosane ( Tc = 767 K, Pc = 11.1 bar) have critical
250
250
(a)
/b)
Ca=25 ]
--
200
200
\
~" 150
~" 150 .o
#. 100
\
CO~ 50
0
200
i
I
250
300
350
\
CO2 50
I
I
I
~
I
400
450
500
550
600
Temperature (K)
\
lOO
650
200
250
300
350
400
450
500
550
600
650
Temperature (K)
Fig. 5. Effect o f the m a x i m u m n u m b e r o f c a r b o n s (C M) on the critical c u r v e o f the s y s t e m carbon dioxide + [ c o n t i n u o u s fraction] ( M W + = 122, -r/= 6.5); (a) C M = 25 (q5 = 25.5) and (b) C M = 20 ( ~ = 20.5).
150
G.L. Rochocz et al. / Fluid Phase Equilibria 139 (1997) 137 153
temperatures close to those of the continuous fractions with C N = l0 and C N = 12, but very different critical pressures. Fig. 4a,b shows the effect of MW + on the critical diagrams of systems in which the discrete component is CO 2 and the continuous fraction contains hydrocarbons from C v to C4.s. For MW += 108, there is a continuous vapor-liquid critical line between the critical points of CO 2 and of the continuous fraction, and there is a liquid-liquid critical line that proceeds to high pressures at approximately 250 K. The shape of this P - T projection is characteristic of type II systems. For MW += 122, the vapor-liquid critical line is interrupted at critical endpoints, and a critical diagram of type IV results. At MW += 136 a diagram of type IV is also calculated to occur, but for MW += 150 a diagram of type II! is obtained. The transition from type IV to type III critical behaviors occurs at a double critical endpoint [37] but we did not attempt to determine its thermodynamic coordinates and the value of MW + at which this occurs. Type IV behavior appeared here in a relatively narrow range of MW+; interestingly, CO 2 + n-tridecane is one of the few binary hydrocarbon systems that show type IV behavior, as such critical behavior also occurs only in a narrow range of heavy component molecular weights for real mixtures. Because the occurrence of type IV diagrams in the calculations are so sensitive to the parameters of the components in the mixture, it is interesting to use these parameters to investigate the effect of the minimum and maximum number of carbon assumed that are assumed to exist in a given distribution. The system CO 2 + [continuous fraction] with MW += 122 and hydrocarbons from C v to C45 (Fig. 4a) was selected as the standard for comparison. In Fig. 5a and b, we see the effect of changing the heaviest hydrocarbon from C45 to C2s and C2o. Fig. 5a (C M = 25) still retains the characteristics of a type IV diagram, while Fig. 5b (C M = 20) is of type II. In Fig. 6a and b, the effect of changing the lightest hydrocarbon from C 7 to C7.05 and C 7.J is shown. Conversely to what happens when the heavy end is changed, the topology of these critical diagrams is very sensitive to changes in the minimum carbon number, C m, employed to characterize the continuous fractions. 250
250
(a)
~
- -
cm=7°5
(b)
1
200
200
150
~" 15o
\
~- loo
/
\
C02 \,
5O
200
Cr.=7.1
lOO
C02
0
[--
I
I
[
I
I
I
250
300
350
400
450
500
Temperature (K)
50
\,
J
550
600
0
650
200
I
I
I
1
I
250
300
350
400
450
~
500
f - -
550
600
650
Temperature (K)
Fig. 6. Effect o f the m i n i m u m n u m b e r o f c a r b o n s ( C , , ) on the critical curve o f the s y s t e m carbon dioxide + [continuous fraction] ( M W + = 122, 4~ = 45.5); (a) C m = 7.05 (rl = 6.55) a n d (b) C m = 7.10 ( ~ / = 6.60).
G.L. Rochocz et aL / Fluid Phase Equilibria 139 (1997) 137-153
151
5. Conclusions The first objective of this paper was to show how the formulation of critical points of continuous mixtures fits into the general framework of integral operators and functional analysis. The formalism obtained in this context is readily generalized to the formulation of the conditions for critical transitions of higher order such as tricritical points, tetracritical points, etc. The second objective was to calculate critical points in systems containing a discrete component and a continuous fraction in order to study the effect of the distribution parameters on the critical projections. For this we used a generalized form of the Gauss-Laguerre quadrature to generate pseudocomponents and then used a modified version of the Hicks and Young [25] algorithm to perform the critical point calculations. The Peng-Robinson equation of state was used with a simple procedure to estimate the properties of the pseudocomponents. In our examples, three quadrature points were sufficient to provide acceptable representations of the systems studied. Despite the simple characterization procedure for the continuous fractions, by varying the light component and the average molecular weight of fraction, we could find five types of the van Konynenburg and Scott classification of critical projections in the P - T plane that cubic equations of state with one-fluid mixing rules are known to exhibit for binary mixtures. Using a type IV diagram as standard for comparison, the influence on the calculations of the mimimum and maximum number of carbons in the distribution was studied. Large variations of the maximum number of carbons could be made before the type of critical phase diagram changed. However, small variations in the minimum carbon number were sufficient to change the type of critical phase diagram obtained.
6. Notation a
A b c
C Cm CM CN D F
h(t,s,u) I Kij
k(s,t) MW + n n c r/d F/f
EOS attractive parameter Helmholtz free energy EOS covolume parameter Cubic form of the Helmholtz free energy expansion Auxiliary parameter of the distribution function Minimum number of atoms of carbon in the molecules of a continuous family Maximum number of atoms of carbon in the molecules of a continuous family Average number of carbons per molecule in a continuous mixture Parameter of the distribution function Intensive distribution function Third order kernel of integral operator T Characterization parameter of the continuous fraction Binary interaction parameter Kernel of the integral operator T Average molecular weight of a continuous fraction Total number of moles Number of components Number of discrete components Number of continuous families of components
152
G.L. Rochoc: et al. / Fluid Phase Equilibria 13911997) 137-153
Number of moles of component I Number of quadrature points Pressure Quadratic form of the Helmholtz free energy expansion Hessian matrix of the Helmholtz free energy expansion Temperature Total volume Weights of Gaussian quadrature Molar fraction of component i Points of Gaussian quadrature
n I
nq
P q
Q T V x i Zj
Greek letters
Parameter of the distribution function Parameter of the Peng-Robinson equation of state Eigenvalue Integral operator Lower limit of the distribution index Chemical potential Upper limit of the distribution index Acentric factor
A K
A T t-t
4, o)
Subscripts C
Critical property
i,j,k
Refer to components i,j,k respectively
Acknowledgements G.L.R. and M.C. acknowledge the comments and suggestions of Mfircio L.L. Paredes in the early stages of this work. The preparation of this manuscript was supported, in part, by Grant No. DE-FG02-85ER13436 from the U.S. Department of Energy, and Grant No. CTS-9521406 from the U.S. National Science Foundation, both to the University of Delaware. M.C. acknowledges the financial support of the Brazilian Ministry of Education (CAPES).
References [1] [2] [3] [4] [5] [6]
K.S. Pedersen, P. Thomassen, Aa. Fredenslund, Fluid Phase Equilibria 14 (1983) 21)9. K.S. Pedersen, P. Thomassen, Aa. Fredenslund, Ind. Eng. Chem. Process Des. Dev. 23 (1984) 163. K.S.P. Pedersen, P. Thomassen, Aa. Fredenslund, Ind. Eng. Chem. Process Des. Dev. 24 (1985) 948. J.R. Bownman, Ind. Eng. Chem. 41 (1949) 2004. W.C. Edmister, D,H. Buchanan, Chem. Eng. Prog. Symp. Ser. No. 6 (1953) 69. M.T. R~itzch, H. Kehlen, Fluid Phase Equilibria 14 (1983) 225.
G.L. Rochocz et al. / Fluid Phase Equilibria 139 (1997) 137-153
[7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37]
153
H. Kehlen, M.T. R~itzsch, Z. Phys. Chem. (Leipzig) 265 (1984) 1049. H. Kehlen, M.T. Rfitzsch, J. Bergmann, AIChE J. 31 (1985) 1136. R.L. Cotterman, J.M. Prausnitz, Revue de L'Institut Fran~ais du P6trole 45 (1990) 633. G.L. Rochocz, MSc Thesis, COPPE, Federal University of Rio de Janeiro, Brazil, 1990. T.J.P. Halpin, N. Quirke, SPE Reservoir Eng. 5 (1990) 617. H.W. Haynes Jr., M.A. Mathews, Ind. Eng. Chem. Res. 30 (1991) 1911. D.S. Chakravarty, M.A. Mathews, Ind. Eng. Chem. Res. 33 (1994) 1962. K.D. Luks, E.A. Turek, T.K. Kragas, Ind. Eng. Chem. Res. 32 (1993) 1767. J. Kramarz, A. Wyczesany, Chem. Eng. Sci. 48 (1993) 1665. C. Lira-Galeana, A. Firoozabadi, J.M. Prausnitz, Fluid Phase Equilibria 102 (1994) 143. H. Kehlen, M.T. Rfitzsch, J. Bergmann, J. Macromol. Sci. Chem. A24 (1987) 1. D. Browarzik, H. Kehlen, M.T. R~itzsch, J. Macromol. Sci. Chem. A27 (1990) 549. Y. Hu, X. Xing, D.T. Wu, J.M. Prausnitz, Fluid Phase Equilibria 104 (1995) 229. D. Browarzik, H. Kehlen, Fluid Phase Equilibria 123 (1996) 17. R.A. Heidemann, A.M. Khalil, AIChE J. 26 (1980) 769. D.-Y. Peng, D.B. Robinson, Ind. Eng. Chem. Fundam. 15 (1976) 59. M. Castier, S.I. Sandier, Critical point calculations with the Wong-Sandler mixing rule: I. Calculations with the van der Waals Equation of State, Chem. Eng. Sci. 52 (1997) 3393. M. Castier, S.I. Sandler, Critical point calculations with the Wong-Sandler mixing rule: II. Calculations with a modified Peng-Robinson Equation of State, Chem. Eng. Sci., in press. C.P. Hicks, C.L. Young, J. Chem. Soc., Faraday Trans. II 73 (1977) 597. B. Epstein, Linear Functional Analysis: An Introduction to Lebesque Integration and Infinite Dimensional Problems, Saunders, Philadelphia, PA, 1970. M.L. Michelsen, Fluid Phase Equilibria 8 (1982) 1. S.K. Shibata, S.I. Sandier, R.A. Behrens, Chem. Eng. Sci. 42 (1987) 1977. D.L. Katz, A. Firoozabadi, J. Pet. Technol. 20 (1978) 1649. J.C. Wheeler, Rocky Mountain J. Math. 4 (1974) 287. G.H. Golub, J.H. Welsch, Math. Comput., 23 (1969) 221-230 and A I - A I 0 . W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in FORTRAN, 2nd edn., Cambridge Univ. Press, Cambridge, UK, 1992. R.C. Reid, J.M. Prausnitz, B.E. Poling, The Properties of Gases and Liquids, 4th edn., McGraw-Hill, New York, 1987. D.B. Robinson, D.-Y. Peng, Research Report RR-28, Project 756, Gas Processors Association, 1978. C.H. Whitson, Soc. Petrol. Eng. J. (August, 1983) 683. P.H. van Konynenburg, R.L. Scott, Philos. Trans. R. Soc. London, Ser. A 298 (1980) 495. A. van Pelt, PhD Thesis, The Technical University of Delft, Holland, 1992.