A reduction method for phase equilibrium calculations with cubic equations of state

THERMODYNAMICS

A reduction method for phase equilibrium calculations with cubic equations of state

D. V. Nichita

Programa de Yacimientos Naturalmente Fracturados, Instituto Mexicano del Petroleo, Eje Central Lázaro Cárdenas 152, 07730 México D.F., México

ABSTRACT

In this work we propose a new reduction method for phase equilibrium calculations using a general form of cubic equations of state (CEOS). The energy term in the CEOS is a quadratic form, which is diagonalized by applying a linear transformation. The number of the reduction parameters is related to the rank of the matrix C with elements (1-C_ij), where C_ij denotes the binary interaction parameters (BIPs). The dimensionality of the problem depends only on the number of reduction parameters, and is independent of the number of components in the mixture.

Keywords: Equation of state; Binary interaction parameters; Reduction method; Linear transformation.

INTRODUCTION

For two-phase flash calculations, the algorithms used to ensure the convergence toward the solution iterates on nc (number of components in the mixture) independent variables (which can be mole fractions, number of moles or logarithms of equilibrium constants). In many petroleum and chemical engineering applications, which may require a considerable number of flash calculations, it is practically impossible to have an extended description of mixtures composition, because this imply solving large dimension problems. Usually, individual components are lumped into pseudocomponents to reduce the problem dimensionality. An alternative way of reducing the dimensionality of the problem is the use of the so-called reduction methods.

In any reduction method, a quadratic form is replaced by the sum of a small number of scalar products. This leads to a system of a reduced number of equations, usually much less than nc. The number of independent variables is not dependent on the number of components in the mixture.

The first reduced flash model was presented by Michelsen (1986). Michelsen showed that the phase equilibrium problem can be solved using only three independent variables if all BIPs in the CEOS are zero. Michelsen's three-equation flash is extremely efficient, but the restriction on BIPs may be unacceptable for many actual problems. Jensen and Fredenslund (1987) extended the Michelsen method, for nonzero BIPs of only one component in the mixture, by solving a system of only five equations. The reduction theorem introduced by Hendriks (1988) establish the conditions under which a reduced number of equations can be solved for flash calculations and phase stability. A set of independent variables for reduced flash using the spectral decomposition (SD) was given by Hendriks and van Bergen (1992). If the matrix has only m nonzero eigenvalues, then the number of independent variables is m+2. A different reduction method, requiring 2c+3 (where c is the number of components with nonzero BIPs) independent variables was presented by Nichita and Minescu (1998, 2004). For an exact solution, the 2c+3 method needs the same number of independent variables as SD. However, SD provides good approximation of the solution if terms corresponding to small nonzero eigenvalues are neglected. Care must be taken when some nonzero eigenvalues are removed, because the one-diagonal condition is no more fulfilled, and proper scaling is required. More recently, the hypersphere decomposition of has been used to reduce dimensionality, even for full-ranked matrices (Nichita, 2004b).

For any reduction method, the reduction parameters are defined as (Hendriks, 1988)

and

where M=m+1, and q_ai;a=1,m;i=1,nc are the elements of the reduction matrix (that depend on the reduction method used), with qMi = B_i.

The vector of the M reduction parameters is

The key factor that allow to solve a system of equations of a reduced dimensionality for the flash problem is that the fugacity coefficient depends at given pressure and temperature only on the reduction parameters, and not on composition,j=ji(p,T,Q).

The reduction method proved to be a useful tool for phase equilibrium calculation. We have solved different kinds of phase equilibrium problems: two-phase flash calculations (Nichita and Minescu, 1998, 2004, Nichita et al., 2003), phase stability analysis (Nichita et al., 2002), multiphase equilibrium calculations (Nichita et al., 2004), critical point calculation (Nichita, 2005, 2006), and phase envelope construction for mixtures with many components (Nichita, 2004a).

This paper presents a new reduction method and illustrates the diversity of approaches suitable for reducing the dimensionality of the phase equilibrium problems. The proposed method is based on the procedure suggested by Tisza (1977) for quadratic forms diagonalization by using linear transformations.

THE CUBIC EQUATION OF STATE

In this work, a general form of two-parameter CEOS is used. It incorporates the Soave-Redlich-Kwong (SRK, Soave, 1972) and Peng-Robinson (PR, Peng and Robinson, 1976) CEOS. However, it is worth noting that any EOS that observes the restrictions of the reduction theorem can be used.

For the SRK CEOS d₁=0 and d₂=1, and for the PR .

With, A=ap/R²T², B=bp/RT, and Z=pv/RT, the implicit form of the CEOS is obtained

The van der Waals mixing rules are used for the energy, A, and for the volume, B, coefficients of the CEOS

^where:

In Eq. (9) and (10), W_a, W_b, and m_i (w_i) are EOS dependent. Their particular values for the SRK and the PR CEOS can be found for example in Michelsen (1986).

The fugacity coefficients are given by

with:

and D = d1 - d2.

THE PROPOSED REDUCTION METHOD

The energy term in the CEOS given by the van der Waals mixing rules is a quadratic form

which can be written as

where ,

and

Obviously, u_ij=u_ji.

The key of the reduction method is to express the CEOS energy parameter A as

that is, to diagonalize the quadratic form (14).

Usually, the diagonalization is performed by spectral decomposition of the matrix (Hendriks and van Bergen, 1992). In this work we propose a different approach for diagonalization, by using a linear transformation (Tisza, 1977).

The transformation connecting mole fractions (via x_i) to Q_i is of the "triangular" form

where

Eq. (18) reads

and represents a linear transformation of the form

that connects the reduction parameters to mole fractions via x. The elements of the matrixT are

The matrix U is singular in many cases, i.e., its rank r<nc. If the rank of U is r=nc-s, the last s diagonal elements in Eq. (17) are vanishing, l_r+1=...l_nc=0, therefore we have

The coefficients l_i can be calculated by a step-by-step construction of the linear transformation in a straightforward manner. However, a simple procedure suggested by Tisza (1977) can be used, starting from the observation that the linear transformation (18) is unimodular, hence the determinant of U is invariant

with

Eq. (24) also holds for the principal minors D_k of D_nc

where

If the matrix U is not full-ranked, r<nc, then .

From Eq. (26), we have

The first m reduction parameters are

where m=r, and the last reduction parameter Q_M, for M=m+1 is given by Eq. (2).

The elements of the matrix U can be written as

The coefficient y_i in the CEOS is

and combining Eq. (30) with Eq. (31), we obtain

Because l_a=0 for a=r+1,...,nc the summation in Eq. (17) is only up to m=r, and the CEOS coefficient A is given by

Finally, at given pressure and temperature, the fugacity coefficients are function of the reduction parameters, and independent on composition

where the compressibility factor depends at given p andT only on the reduced variables, Z=Z(Q), if A from Eq. (6) is replaced in Eq. (5). The expression of the fugacity coefficients is exactly the same as used by Nichita et al. (2003, 2004); here l_a are those obtained from Eq. (28) instead of eigenvalues of U, and t_ai are replacing the corresponding eigenvectors. Note that for all BIPs zero, the Michelsen's three-equation flash is obtained, by putting l₁=nc, and .

The calculations proceed as follows: the minors D_k are calculated for increasing k, starting with D₁=u₁₁=1 (for convenience D₀=1), and l_k are calculated with Eq. (28); k is increased until D_k=0 (more precisely D_k<e). The last k giving a non-zero determinant is equal to the rank of U. Components must be properly ordered, the first c components being those with non-zero BIPs. If the BIPs between a component and two other consecutive components are equal (C_ij=C_i,j+1), and this leads to D_k=0 for some k<r, one of these BIPs is altered by a small perturbation (say 1%); this does not affect phase equilibrium results, but prevents computational problems.

The structure of the matrix U is doubly bordered, the rank depending on the number of components having non-zero BIPs. For hydrocarbon mixtures this matrix is generally rank-deficient, i.e. singular. The proposed method requires usually only the calculation of low order determinants, up to r. For systems with many components, this avoids matrix operations for large dimensions (calculation of eigenvalues and eigenvectors).

Implementation of the proposed reduction method requires only minor changes in the existing codes based on SD. The subroutine for eigenvalues and eigenvectors calculation is replaced with the subroutine based on the proposed method. None of the phase equilibrium routine needs any modification. Only different formal parameters are transferred to these routines, l_k instead of eigenvalues, and t_ij from Eq. (29) instead of eigenvectors.

APPLICATIONS

A variety of phase equilibrium problem we have studied before has been addressed using the proposed reduction method. Applications for mixtures with nc ranging from 6 to 52, described in previous papers on two-phase flash (Nichita et al., 2003, Nichita and Minescu, 2004), phase stability testing (Nichita et al., 2002), three-phase vapor-liquid-liquid flash (Nichita et al., 2004), critical point calculation (Nichita, 2005, 2006), phase diagram construction (Nichita, 2004a), were reworked.

For all situations, the results are reproduced using the proposed reduction method. Phase equilibrium calculations using the proposed reduction method and the SD reduction method requires approximately the same number of iterations, and thus the same computational effort.

The MY10 Mixture

The first exemplification is given here for the MY10 mixture, taken from Metcalfe and Yarborough (1979). The MY10 synthetic mixture contains 10 normal-alkanes, with composition and BIPs of methane (taken from Firoozabadi and Pan, 2002) with the remaining components given in Table 1.

The non-zero eigenvalues for SD and l_k for the proposed method are shown in Table 2. The rank of U is r=3, giving 5 independent variables for two-phase equilibrium calculations.

The 11-component MY10/CO₂ mixture is obtained by adding CO₂ in different proportions to the MY10 mixture. Table 3 gives the non-zero eigenvalues for SD and l_k for the proposed method in two cases:

a) CO₂ BIPs from Table 1 (r=5)

b) All CO₂ BIPs equals to 0.12 (r=4)

As can be seen from Tables 2 and 3, the number of non-zero l_k in our reduction method is equal to the number of eigenvalues in the SD method. Therefore the number of reduction parameters (M=r+1) and of independent variables for phase equilibrium calculations (for example M+1=r+2 for a two-phase flash) is the same for the two approaches.

For case a, calculations are detailed below step-by-step. The component ordering is: CO₂, C1, nC4 to nC14, C2, C3.

The determinant D for k=6 is

We start with D₁=1 (D₀ is set to 1 for convenience) and l₁=1. For k=2, D₂=0.177351 and; for k=3, D₃=6.405e-3 and; for k=4, D₄=-4.927e-6, giving, then D₅=2.97e-9 for k=5 and. For k=6 the determinant is zero, D₆=0, and all other determinants D_k, k>6 are also zero, thus l_k=0 for k=6, ,11, and we keep for l only the first five nonzero values listed in Table 3; The elements of the transformation matrixT are given by Eq. (22). This detailed calculation show clearly how easy to implement the proposed method is.

A 52 Component Mixture

The next example is for a natural occurring hydrocarbon mixture with many components. Sample C from Pedersen et al. (1985) is a heavy gas-condensate for which a detailed description of the C₇₊ fraction is available, including the paraffin-naphtene-aromatic distribution. The mixture has 52 components, critical properties and acentric factors being calculated as described in Nichita (2005). The BIPs of methane with other hydrocarbon components are assigned according to Katz and Firoozabadi (1978), and BIPs of carbon dioxide and nitrogen with the hydrocarbon components are C_CO2-j=0.12 and C_N2-j=0.1.

The rank of matrix U is r=5. The five nonzero eigenvalues and nonzero values of l_k obtained by the reduction method are given in Table 4.

Note that only determinants up to order five are computed to complete the reduction for this 52 component mixture.

CONCLUSIONS

A reduction method is proposed, based on diagonalization of quadratic forms by means of linear transformations. The number of reduction parameters is the same as for reduction method based on spectral decomposition.

The codes for phase equilibrium calculation with the spectral decomposition reduction method can be used for the proposed reduction method without any modification of the phase equilibrium routines. Some previously addressed problems were reworked for the new reduction method. Results are identical and the number of iterations required for phase equilibrium calculation, and thus the computational effort is almost the same as for the spectral decomposition method.

The proposed method requires usually only the calculation of low order determinants. For systems with many components, this avoids matrix operations for large dimensions (calculation of eigenvalues and eigenvectors).

ACKNOWLEDGEMENT

I gratefully acknowledge the support from Instituto Mexicano del Petroleo, Programa de Yacimientos Naturalmente Fracturados, under project no. D00084.

NOMECLATURE

Symbols

A

attractive parameter in the CEOS

A_i

component parameter in CEOS

A_ij

CEOS cross parameter

a

attractive parameter in the CEOS

a_ij component CEOS coefficient B

volume parameter in the CEOS

B_i

component CEOS coefficient

B

covolume in the CEOS

matrix with elements 1-C_ij

C_ij

binary interaction coefficients between components i and j

D_nc

determinant of the matrix U

D_k

principal minors of D_nc

M

number of reduction parameters

m

number of nonzero eigenvalues

nc

number of components

p

pressure

Q

vector of reduction parameters

Q_a

reduction parameters

q_ai

elements of the eigenvectors

R

universal gas constant

r rank of matrix T

temperature

T

transformation matrix

t_ij

elements of T

U

identic with

u_ij

elements of U

v

molar volume

x_i

mole fraction, component i

Z

compressibility factor

Greek Letters

d₁, d₂ constants in CEOS D d₁-d₂ j_i fugacity coefficients l

eigenvalues

y_i

component CEOS coefficient

x_i

given by Eq. (15)

W_a, W_b

coefficients in the CEOS

w

acentric factor

Subscripts

c

critical

i,j,k

component index

r

reduced

a

reduction parameter index

Received: December 12, 2005

Accepted: May 4, 2006

* To whom correspondence should be addressed

^* Current address: UMR 5150 CNRS, Laboratoire des Fluides Complexes, Phone: +(33) (5) 5940 7685, Fax: +(33) (5) 5940 7725, Université de Pau et des pays de l'Adour, BP 1155, 64013 Pau Cedex, France. E-mail: [e-mail]

• Firoozabadi, A. and Pan, H., Fast and Robust Algorithm for Compositional Modeling: Part I Stability Analysis Testing, SPE J., 7, 78-89 (2002).
• Hendriks, E. M., Reduction theorem for phase equilibria problems, Ind. Eng. Chem. Res., 27, 1728-1732 (1988).
• Hendriks, E. M. and van Bergen, A.R.D., Application of a Reduction Method to Phase Equilibria Calculations, Fluid Phase Equilibria, 74, 17-34 (1992).
• Jensen, B.H. and Fredenslund Aa., A Simplified Flash Procedure for Multicomponent Mixtures Containing Hydrocarbons and One Non-Hydrocarbon Using Two-Parameter Cubic Equations of State, Ind. Eng. Chem. Res., 26, 2129-2134 (1987).
• Katz, D.L. and Firoozabadi, A., Predicting Phase Behavior of Condensate/Crude-Oil Systems Using Methane Interaction Coefficients, J. Petrol. Technol., 30, 1649-1655 (1978).
• Metcalfe, R.S. and Yarborough, L., Effect of Phase Equilibria on the CO₂ Displacement Mechanism", Soc. Petr. Eng. J, 19, 242-252 (1979).
• Michelsen, M.L., Simplified Flash Calculations for Cubic Equation of State, Ind. Eng. Chem. Process Des. Dev., 25, 184-188 (1986).
• Nichita, D.V., 2004. Phase Diagram Construction for Mixtures with Many Components, Proc. The 16^th International Congress of Chemical and Process Engineering CHISA2004, Prague, August 22-26, 2004. submitted.
• Nichita, D.V., 2004. The Rank-Reduction Method Applied to Phase Equilibrium Calculation, IMP Report, December 2004, to be submitted.
• Nichita, D.V., Calculation of Critical Points Using a Reduction Method, Fluid Phase Equilibria, 228-229, 223-231 (2005).
• Nichita, D.V., A New Method for Critical Points Calculation from Cubic EOS, AIChE J., 52, 1220 1227 (2006).
• Nichita, D.V., Broseta, D. and de Hemptinne, J.-C., 2004. Multiphase Equilibrium Calculation Using Reduced Variables, paper SPE 89439, Proc. The SPE/DOE 14^th Symposium on Improved Oil Recovery, Tulsa, OK, April 17-21, 2004, submitted to Fluid Phase Equilibria.
• Nichita, D.V., Broseta, D., de Hemptinne, J.-C. and Lachet, V., 2003. Efficient Phase Equilibrium Calculation for Compositional Simulation, Proc. The 20^th European Symposium on Applied Thermodynamics, held on Oct. 9-12, 2003, Lahnstein, Germany, in press, Petroleum Science and Technology.
• Nichita, D.V., Gomez, S. and Luna, E., Phase Stability Analysis with Cubic Equations of State Using a Global Optimization Method, Fluid Phase Equilibria, 194-197, 411-437 (2002).
• Nichita, D.V. and Minescu, F., 1998. A New, Efficient Phase Equilibria Calculation Method in a Reduced Flash Context, The 13^th International Congress of Chemical and Process Engineering CHISA'98, August 23-28, 1998, Prague.
• Nichita, D.V. and Minescu, F., Efficient Phase Equilibrium Calculation in a Reduced Flash Context, Can. J. Chem. Eng., 82, no. 6 1225-1238 (2004).
• Pedersen, K.S., Thomassen, P. and Fredenslund, Aa., Thermodynamics of Petroleum Mixtures Containing Heavy Hydrocarbons. 3. Efficient Flash Calculation Procedures Using the SRK Equation of State, Ind. Eng. Chem. Process Des. Dev., 24, 948-954 (1985).
• Peng, D.Y. and Robinson, D.B., A New Two-Constant Equation of State, Ind. Eng. Chem. Fundamentals, 15, 59-64 (1976).
• Soave, G., Equilibrium Constants From a Modified Redlich-Kwong Equation of State, Chem. Eng. Sci., 27, 1197-1203 (1972).
• Tisza, L., Generalized Thermodynamics, The MIT Press (1977).

Publication Dates

• Publication in this collection
  13 Dec 2006
• Date of issue
  Sept 2006

History

• Accepted
  04 May 2006
• Received
  12 Dec 2005