488
Energy & Fuels 2008, 22, 488–495
Phase Envelope Construction for Mixtures with Many Components Dan Vladimir Nichita* Instituto Mexicano del Petroleo, Programa de Yacimientos Naturalmente Fracturados, Eje Central Lázaro Cárdenas 152, 07730 México D.F., México, and CNRS UMR 5150, Laboratoire des Fluides Complexes, UniVersité de Pau et des Pays de l’Adour, B.P. 1155, 64013 Pau Cedex, France ReceiVed May 17, 2007. ReVised Manuscript ReceiVed September 14, 2007
A reduction method for constructing vapor–liquid equilibrium phase envelopes with cubic equations of state is presented. The paper describes the calculation procedures for saturation (dewpoint/bubblepoint) pressures and temperatures, quality lines (for given mole fraction or volume fraction of one of the equilibrium phases), cricondentherm and cricondenbar points, the critical point, and the spinodal. The phase envelope construction is fully automatic. The saturation points are calculated throughout the critical region by stepping around the phase envelope in pressure or temperature increments. An extrapolation procedure taking advantage of the Jacobian matrix available from a previous step is used. Problem formulation in terms of reduced variables leads to simpler partial derivatives with respect to pressure and temperature. The proposed method is successfully tested for several representative hydrocarbon mixtures with various phase envelope shapes.
1. Introduction The construction of mixture phase envelopes is an important problem in phase equilibrium calculations, with many applications in chemical thermodynamics and hydrocarbon reservoir thermodynamics. It is very useful for process simulation and reactor design and for the study of hydrocarbon reservoir fluids. The most popular approaches for phase envelope construction are those of Michelsen,1,2 Li and Nghiem,3 and Ziervogel and Poling.4 The main difficulties in constructing phase envelopes are pointed out in the literature.1–3,5,6 Saturation points calculation is one of the basic phase equilibrium calculation problems. The number of solutions for a given specification is not known in advance; there may be multiple solutions or no solution at all. The fact that the solution is not unique imposes the use of carefully chosen initial estimates; otherwise, incorrect trivial or nontrivial solutions are found. Moreover, in the near critical region, it is not clear whether the saturation points are dewpoints (DP) or bubblepoints (BP). Complete phase envelope construction of multicomponent mixtures requires the resolution of nonlinear systems of equations of the order nc + 1 (nc is the number of components) for saturation pressures/temperatures and quality lines and evaluation of determinants of order nc for critical point and spinodal calculation. The computer time spent for linear algebra operations involved in these calculations is proportional at least to nc2. * Cooresponding author, currently with Université de Pau et des pays de l’Adour. Tel.: +33 5 5940 7685. Fax: +33 5 5940 7725. E-mail addresses:
[email protected];
[email protected]. (1) Michelsen, M. L. Fluid Phase Equilib. 1980, 4, 1–10. (2) Michelsen, M. L. Fluid Phase Equilib. 1984, 16, 57–76. (3) Li, Y-K.; Nghiem, L. X. The Development of a General Phase Envelope Construction Algorithm for Reservoir Fluid Studies. Presented at the 57th Annual Fall Technical Conference and Exhibition of the SPE of AIME, New Orleans, LA, September 26–29, 1982; paper SPE 11198. (4) Ziervogel, R. G.; Poling, B. E. Fluid Phase Equilib. 1983, 11, 127– 135. (5) Michelsen, M. L. Fluid Phase Equilib. 1985, 23, 181–192. (6) Nghiem, L. X.; Li, Y-K.; Heidemann, R. A. Fluid Phase Equilib. 1985, 21, 39–60.
For hydrocarbon mixtures, one may often need a detailed description of heavy fractions composition (including the paraffin–naphtene–aromatic distribution),7,8 leading to mixtures with many components. It is a common practice to lump several components into a pseudocomponent to reduce the overall number of components. However, if at some conditions, the detailed composition is again required, a delumping (or inverse lumping) procedure9,10 is used to recover it. Lumped and delumped phase envelopes differ (in some cases to a high extent) from the original phase envelope. This work is based on a different way of reducing the dimensionality of phase equilibrium calculations, namely the reduction method.11,12 The idea behind the application of any reduction method is to find a different formulation of the problem, based on a different set of independent variables, whose number is smaller, or preferably much smaller, than for the original set of variables. Several reduction procedures have been proposed in the literature.11–16 In terms of phase equilibrium calculation, this consists in replacing the traditional variables (mole numbers, mole fractions, partition coefficients) with some linear combinations of them. Consequently, the solution is no more sought in the compositional hyperspace (defined by mole fractions), but in the reduced hyperspace (defined by the new introduced variables). If for the reduction method the dimen(7) Nichita, D. V.; Goual, L.; Firoozabadi, A. SPE Prod. Facil. 2001, 16, 250–259. (8) Pedersen, K. S.; Thomassen, P.; Fredenslund, Aa. Ind. Eng. Chem. Process Des. DeV. 1985, 24, 948–954. (9) Nichita, D. V.; Leibovici, C. F. Fluid Phase Equilib. 2006, 245, 71–82. (10) Nichita, D. V.; Broseta, D.; Leibovici, C. F. Comput. Chem. Eng. 2006, 30, 1026–1037. (11) Michelsen, M. L. Ind. Eng. Chem. Process Des. DeV. 1986, 25, 184–188. (12) Hendriks, E. M. Ind. Eng. Chem. Res. 1988, 27, 1728–1732. (13) Jensen, B. H.; Fredenslund, Aa. Ind. Eng. Chem. Res. 1987, 26, 2129–2134. (14) Hendriks, E. M.; van Bergen, A. R. D. Fluid Phase Equilib. 1992, 74, 17–34. (15) Nichita, D. V.; Minescu, F. Can. J. Chem. Eng. 2004, 82, 1225– 1238. (16) Nichita, D. V. Braz. J. Chem. Eng. 2006, 23, 427–434.
10.1021/ef700251k CCC: $40.75 2008 American Chemical Society Published on Web 11/17/2007
Phase EnVelope Construction
Energy & Fuels, Vol. 22, No. 1, 2008 489
sionality of the problem (number of independent variables) is much smaller than for the classical methods, one can expect important decreases of computational effort. Besides, it seems that hypersurfaces in the reduction methods are smoother than those in conventional ones, and thus, robustness may be improved. Previously, we have successfully applied different reduction methods to solve various kinds of phase equilibrium problems: two-phase flash,15,17 phase stability analysis,18,19 multiphase flash,19 extrapolations,15,17 and critical points.20,21 In this work, we use for phase diagrams construction a reduction method applied for the general form of two-parameter cubic equations of state (EoS), including the Soave–Redlich–Kwong (SRK)22 EoS and the Peng–Robinson (PR)23 EoS. The reduction is effectively achieved by spectral decomposition of the matrix C with elements (1 – Cij), where Cij denotes the binary interaction parameters (BIPs). The dimensionality of the problem depends only on the number of reduction parameters, M ) m + 1 (where m is the number of nonzero or non-negligible eigenvalues of the matrix C ) and is independent of the number of components nc. The reduction method is extremely efficient for mixtures with many components and relatively few nonzero BIPs; this is the case for most practical applications involving hydrocarbon mixtures and classical contaminants (CO2, N2, H2S, H2O). We present a reduction procedure to calculate phase boundaries and key points on the phase envelope such as the critical point, cricondentherm (CMT), and cricondenbar (CMP). A recently developed reduction method is used for critical point calculation (and optional for spinodal construction). Here, we detail calculation of BP and DP pressures and temperatures, quality lines, CMT, and CMP, all calculation being based on the reduction method. The phase envelope construction is fully automatic; most features are adopted from previous algorithms,1,4 but in this work, all calculations are performed using reduction methods. The bubblepoints and dewpoints are calculated throughout the critical region by stepping around the phase envelope in pressure or temperature increments. An extrapolation procedure taking advantage of the Jacobian matrix available from a previous step is used. Specification of the independent variables is also done automatically. The proposed method is tested for several hydrocarbon mixtures having various phase envelope shapes, for nc values ranging from 6 to 52 and M ranging from 2 to 10. 2. Reduced Variables for Cubic EoS The general form of the two-parameter cubic EoS11,24 is used in this work RT a p) vm - b (vm + δ1b)(vm + δ2b)
(1)
The SRK EoS and PR EoS are particular cases of eq 1. For the PR EoS, δ1,2 ) 1 + 2 and for the SRK EoS, δ1 ) 0 and (17) Nichita, D. V.; Broseta, D.; de Hemptinne, J.-C.; Lachet, V. Pet. Sci. Technol. 2007, 25, 315–342. (18) Nichita, D. V.; Gomez, S.; Luna, E. Fluid Phase Equilib. 2002, 194–197, 411–437. (19) Nichita, D. V.; Broseta, D.; de Hemptinne, J.-C. Fluid Phase Equilib. 2006, 246, 15–27. (20) Nichita, D. V. Fluid Phase Equilib. 2005, 228–229, 223–231. (21) Nichita, D. V. AIChE J. 2006, 52, 1220–1227. (22) Soave, G. Chem. Eng. Sci. 1972, 27, 1197–1203. (23) Peng, D.-Y.; Robinson, D. B. Ind. Eng. Chem. Fundam. 1976, 15, 59–64. (24) Michelsen, M. L.; Mollerup, J. M. Thermodynamic Models: Fundamentals & Computational Aspects; Tie-Line Publications: Holte, Denmark, 2004.
δ2 ) 1. These two equations are currently the most widely used EoSs in chemical and petroleum industry simulations, due to their relatively simple form; high accuracy can be obtained by tuning the EoS parameters to available experimental data. The implicit form (in the compressibility factor, Z ) pvm/RT ) of the two-parameter cubic EoS is Z3 + [(δ1 + δ2 - 1)B - 1]Z2 + [A + δ1δ2B 2 (δ1 + δ2)B(B + 1)]Z - [AB + δ1δ2B (B + 1)] ) 0 (2)
is obtained using the classical notation A ) ap/(RT)2 and B ) bp/RT. The van der Waals mixing rules are used for the energy, A, and for the volume, B, parameters: nc
A)
nc
∑∑
nc
yiyjAij ;
B)
i)1 j)1
∑yB; j j
j)1
Aij ) Aji ) (AiAj)0.5(1 - Cij) ;
i, j ) 1, nc (3)
where Aj )
Ωa prj Trj2
Rj(T) ;
Bj )
Ωb prj ; Trj Rj(T) ) [1 + m(ωj)(1 - Trj0.5)]2 (4)
and the specific values of Ωa and Ωb and the expression of m as a function of the acentric factor ω for each EoS can be found in the literature.22,23 The fugacity coefficient of component i is given by the relation ln φi ) (Z - 1)
(
)(
)
Bi Z + δ1B A ψi Bi - ln(Z - B) 2 ln B ∆B A B Z + δ2B (5)
with ∆ ) δ1 - δ2 and nc
ψi )
∑A y
ij j
(6)
j)1
The key of the reduced phase equilibrium approach is to express EoS parameters (A and B) and ψi (eq 6) in terms of a small number of scalar products, i.e., QR, called “reduction parameters” (see the work of Hendriks12). In other words, all quantities involved in fugacity coefficient expressions must be linear forms or decomposable to linear forms. This can be done by using the spectral decomposition of the symmetric matrix C with elements (1 – Cij), i,j ) 1, nc, which is C ) SΛST (7) where Λ ) diag(λi); i ) 1, nc is the set of eigenvalues associated to the right-hand-side eigensystem S of C . If r ) rank C < nc, the matrix is singular and has only m ) r nonzero eigenvalues. In this case, the diagonal matix Λ ) diag(λi); i ) 1, r contains only the nonzero eigenvalues, S∈R(nc×r), and eq 2 reads m
1 - Cij )
∑λ q
R RiqRj
(8)
R)1
or, by multiplying eq 2 by (AiAj)1/2 m
Aij )
∑λ q
R RiqRj
R)1
where qRi ) qRi′Ai.
(9)
490 Energy & Fuels, Vol. 22, No. 1, 2008
Nichita
Let the reduction parameters be defined as nc
QR )
∑q
R ) 1, m
Riyi ;
(10a)
i)1
and nc
QM ≡ B )
∑By
(10b)
i i
i)1
where M ) m + 1, then the vector of reduction parameters is Q ) (Q1,Q2,...,Qm,B)T. The EoS coefficients A and ψi can be written as nc
A(Q) )
(
nc
m
∑∑yy ∑λ q
) ∑ (∑ ∑ ) ∑ (∑ ) ∑
R RiqRj
i j
i)1 j)1
m
R)1
m
nc
)
nc
qRiqRjyiyj )
λR
R)1 nc
i)1 j)1 2
m
qRiyi )
λR
R)1
λRQR2 (11)
R)1
i)1
and nc
ψi(Q) )
m
∑ ∑ yj
j)1
nc
m
λRqRiqRj )
R)1
∑
λRqRi
R)1
∑q
Rjyj )
j)1 m
∑λ q
R RiQR ;
R)1
i ) 1, nc (12)
Because A and ψi are functions only of the reduction parameters (and not on compositions), B is a reduction parameter itself, and the compressibility factor depends at given p and T only on the reduction parameters, Z ) Z[A(Q),B(Q)] ) Z(Q), the fugacity coefficient (eq 5) can be expressed as a function only of Q: ln φi(Q) ) (Z(Q) - 1)
(∑
qMi 1 - ln(Z(Q) - QM) × QM ∆QM
m
2
λRqRiQR -
R)1
)(
qMi m Z(Q) + δ1QM λRQ2R ln QM R)1 Z(Q) + δ2QM
∑
)
(13)
variables. It is worth mentioning that eq 14 holds for any alternative reduction method.15,16 For hydrocarbon mixtures containing the classical contaminants (carbon dioxide, nitrogen, hydrogen sulfide, water, etc.), usually the BIPs matrix has many zero entries, only a few components having nonzero BIPs with the remaining ones, that is m ) r
