- Email: [email protected]
Prediction of multiphase equilibrium using the PC-SAFT equation of state and simultaneous testing of phase stability
Fluid Phase Equilibria 302 (2011) 169–178
Contents lists available at ScienceDirect
Fluid Phase Equilibria journal homepage: www.elsevier.com/locate/fluid
Prediction of multiphase equilibrium using the PC-SAFT equation of state and simultaneous testing of phase stability Nayef M. Alsaifi, Peter Englezos ∗ Department of Chemical and Biological Engineering, University of British Columbia, Vancouver, BC, Canada V6T 1Z3
a r t i c l e
i n f o
Article history: Received 29 May 2010 Received in revised form 1 September 2010 Accepted 3 September 2010 Available online 15 September 2010 Keywords: Multiphase calculations Stability analysis PC-SAFT Association Dipolar Quadrupole Hydrocarbon Water Ethylene glycol Trust-region Gauss–Newton method
a b s t r a c t This work employs the PC-SAFT equation of state to illustrate how to compute phase equilibria and stability of multiphase systems. The proposed computation scheme is based on the methodology of Gupta et al. (1991) but with modifications. Phase equilibria and stability are computed simultaneously and without the need to modify the exact stability equation. This approach does not encounter difficulty when the phase fraction approaches zero. A simple initialization strategy is utilized. In addition, this work tests the predictive capability of PC-SAFT in multiphase systems of water-hydrocarbon containing systems at different pressure–temperature conditions. The dipole and quadrupole interactions were included in the PC-SAFT equation and the binary interaction parameters (kij ) were set to zero. The work has demonstrated the ability of PC-SAFT to successfully predict vapour–liquid–liquid equilibrium of systems containing water, hydrocarbons and ethylene glycol with or without carbon dioxide. © 2010 Elsevier B.V. All rights reserved.
1. Introduction Computation of equilibrium and stability of multiphase systems plays an essential role in simulation and design of chemical processes such as distillation and extraction. The requirement for efficient and reliable multiphase computation methods is always of interest in terms of the effectiveness and reliability of process simulators. These methods are generally based on either the Gibbs free energy minimization approach or equation-solving techniques. The equation-solving techniques combine mass balances and equilibrium relations into a set of non-linear algebraic equations. Different approaches have been suggested for equation-solving techniques [1–3]. The Gibbs energy minimization approach, however, formulates the phase equilibrium problem into an optimization problem which can be solved by a global solver subject to material balance constraints. Various methodologies have been developed under this method by many researchers [4–10]. Irrespective of which method is adopted, the essence of phase equilibrium calculations is always to ensure the determination of the global minimum Gibbs free energy of a system at fixed isothermal and isobaric conditions.
∗ Corresponding author. Tel.: +1 604 822 6184; fax: +1 604 822 6003. E-mail address: [email protected] (P. Englezos). 0378-3812/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.fluid.2010.09.002
The main challenge in multiphase equilibrium calculations is that the number of the available phases in a system at equilibrium is not known a priori. A common way to overcome the encountered phase difficulty is to utilize stability criteria [11,12] to determine whether an obtained solution is thermodynamically stable. This can be handled either by assuming a fixed number of phases and checking the stability of the final solution [13] or by following a sequential procedure by adding a phase in the computation and testing the stability [5,14–16]. A novel development of stability and multiphase calculations which has some advantages over the previous approaches is the one proposed by Gupta et al. [1,2,17]. The methodology formulates the problem as a set of coupled non-linear equations that describe phase equilibrium, chemical equilibrium and the stability of the system with inequality constraints. Thus, it is an equation-solving technique type that combines mass balance and equilibrium relations as well as stability equations with inequality constraints. The advantages of the Gupta methodology are its simplicity and its ability to provide insight into the phase behaviour as it tracks the appearance or disappearance of the phases. In principle, the method presents the necessary and sufficient conditions for phase equilibrium and it offers a simple and clear stability criterion for multiphase reacting/non-reacting systems at a fixed pressure and temperature. Gupta et al. [1] modified the exact stability equation and replaced it with one that satisfies equivalently the behaviour of
170
N.M. Alsaifi, P. Englezos / Fluid Phase Equilibria 302 (2011) 169–178
appearance and disappearance of a phase. The reason was to avoid the problem of ill-conditioned Jacobian when the system of the non-linear equations is required to be solved simultaneously. The treatment of singular Jacobian was also taken into account with another modification. Furthermore, because Gupta et al. development requires two inequality constraints for each additional phase (stability variable and phase fraction, Equation (6)) to be greater than or equal to zero, it was necessary to readjust them during iterations when the two inequality constraints were violated. This is because the Newton-Raphson method cannot handle constraints when solving a system of non-linear equations. Therefore, the inequality constraints were implicitly included during the computation. In order to guarantee success of the previous modifications, a good initialization strategy of K-values was crucial. The K-values were initiated with the stability method suggested by Michelsen [12]. Although their procedure was successful for several systems, it encounters difficulty when the phase fraction approaches zero [18]. In this study, a different scheme is suggested to perform the simultaneous calculations. The suggested scheme does not encounter a difficulty when the phase fraction approaches zero. Furthermore, it is applied without modifying the stability equation. When the new scheme is applied to water-hydrocarbon containing systems, a simple initialization strategy is followed. Finally, when the system of interest could show more than two phases on the P–T continuum (say vapour–liquid–liquid) and the selected P–T condition is in one phase region (say Vapour phase), then the proposed scheme would easily assist in figuring out which liquid phase would appear first provided that the utilized model is able to qualitatively capture the phase behaviour of the studied system. For the scope of this paper, we restrict ourselves to single-stage multiphase equilibrium calculations. Chemical equilibrium calculations are not included. Moreover, due to the ongoing development of very sophisticated models and the requirement of evaluating these models for multiphase calculations, we are not going to utilize simple models such as cubic equations of state. Instead, perturbed-chain statistical association fluid theory (PC-SAFT) [19] is utilized including association, dipolar and quadrupolar terms. For the dipolar and quadrupolar terms, developments by Gross’ group are utilized in this work [20,21]. Practical systems such as waterhydrocarbon containing systems with ethylene glycol and carbon dioxide are tested and compared to experimental data. Binary interaction parameters are set to zero in all the calculations of this paper.
2. Multiphase and stability equations The appropriate set of equations which describe the stability and the isothermal–isobaric flash calculations of reacting and nonreacting systems are formulated by Gupta et al [1,2]. The details of their derivation are not repeated here. The reader could consult Gupta et al. [1] for more details about the basis of their development. The main governing equations of non-reacting systems at constant pressure and temperature are summarized here as follows.
Table 1 State of phase k in equations (5) and (6).
Case I Case II Case III
Stability variable ( k )
Phase fraction (˛k )
State of phase k
Zero Positive Zero
Positive Zero Zero
Stable Unstable Incipient
2.2. Mole fraction summation N
(kik ek − 1)xir = 0
(k = 1, . . . , ; k = / r)
(2)
i=1
xir =
nc
z
z )(1 + i=1 i
(
i
j j = 1 (kij e − 1)˛j ) j= / r
(i = 1, . . . , nc; j = 1, . . . , ; j = / r)
(3)
2.3. Phase equilibrium and stability relations xik = Kik xir ek ˛k k = 0
(i = 1, . . . , nc; k = 1, . . . , ; k = / r)
(k = 1, . . . , ; k = / r)
(4) (5)
Subject to ˛k ≥ 0
and
k ≥ 0
where ϕ kik = ir ϕik
(6)
(7)
In the above equations, z is feed composition, x is mole fraction, k is the stability variable, nc is the number of components and ϕik is the fugacity coefficient of component i in phase k. Furthermore, subscript r denotes a reference phase. The number of equations required to describe the phase equilibrium and stability of a specific system at a fixed temperature and pressure is equal to [(nc+2) − 1] for [(nc+2) − 1] unknowns, namely () ˛k , ( − 1) k and (nc·)xik . It should be emphasized that Equations (5) and (6) imply that either the phase fraction (˛k ) or the stability variable ( k ) or both must be zero. In case the phase fraction becomes zero, this indicates that the stability variable must be either a positive for an unstable phase or zero for an incipient phase. In a similar manner, if the stability variable becomes zero, then the phase fraction must be positive for a stable phase or zero for an incipient phase. These possible cases for phase (k) are summarized in Table 1. It should be noted that Equations (5) and (6) give the relationship between stability variables and phase fraction at the global minimum of Gibbs energy [2]. Finally, as demonstrated by Gupta et al. [1], at the point of minimum Gibbs free energy, the stability variables must satisfy the following relation which is written in terms of fugacity (fik ): k = ln
f ik
fir
(i = 1, . . . , nc; k = 1, . . . , ; k = / r)
(8)
3. Equation of state The SAFT model is one of the most commonly used theory based equations of state. The general framework of SAFT is written in terms of reduced Helmholtz energy as follows:
2.1. Phase fraction summation
˛k = 1
k=1
where is the number of phases and ˛k is the phase fraction.
(1)
ares = ahs + achain + adisp + aassoc + apolar
(9)
where res, hs, chain, disp, assoc and polar stand for residual, hard sphere, chain, dispersion, association and polar, respectively. The polar term represents the long-ranged electrostatic
N.M. Alsaifi, P. Englezos / Fluid Phase Equilibria 302 (2011) 169–178
171
Table 2 PC-SAFT parameters. Compound
M [g/mol]
m
[Å]
ε/k [K]
Methane Propane Butane Pentane Hexane Octane Carbon dioxide Water Ethylene glycol
16.043 44.096 58.123 72.146 86.177 114.231 44.010 18.015 62.068
1.000 2.002 2.332 2.690 3.058 3.818 1.513 1.226 3.404
3.704 3.618 3.709 3.773 3.798 3.837 3.187 2.792 2.844
150.030 208.110 222.880 231.200 236.770 242.780 163.330 203.000 178.370
interactions such as dipole–dipole (DD), dipole–quadrupole (DQ) and quadrupole–quadrupole (QQ) interactions. In this work, the dipole–dipole interactions are added for water and alcohols to increase the predictive capability of PC-SAFT as demonstrated in a previous study [23]. Furthermore, the quadrupole–quadrupole interactions are taken into account for carbon dioxide in order to improve the predictive capability of SAFT since the quadrupole–quadrupole interactions are major interactions in carbon dioxide. In this work, we mainly utilize perturbed-chain statistical association fluid theory (PC-SAFT) equation of state although other SAFT versions could be used. The simplified SAFT [36] is also used for testing some systems [Table 6]. PC-SAFT has been extensively described in the literature. The main terms of PCSAFT will not be repeated here. The reader could refer to the work of Gross and Sadowski [19] for details. Only the dipolar and quadrupolar terms will be given [20,21]. The dipolar and quadrupole terms are both developed based on the Pade approximate form of Rushbrooke et al. [22]:
apolar =
polar a2 . polar polar 1 − (a3 /a2 )
(10)
where polar superscript denotes dipolar (DD) or quadrupole (QQ) terms. The dipolar terms (DD) are defined as follows:
aDD 2 = −
i
aDD 3 = −
xi xj
j
(11)
42 2 ε εjj εkk xi xj xk ii 3 kB T kB T kB T i
×
3 3 εii εjj ii jj n n ∗2 ∗2 J DD j 2,ij kB T kB T 3 ,i ,j i ij
3 ii3 jj3 kk
ij ik jk
j
k
n,i n,j n,k ∗2 ∗2 ∗2 J DD i j k 3,ijk
(12)
In the above equations, εii is segment energy parameter of component i, ii is segment size parameter of component i, ni is number of dipole moments in component i and ∗2 is dimeni sionless squared dipole moment. J2,ij and J3,ijk are approximated integrals over the reference-fluid pair-correlation function and over three-body correlation functions. The approximations of these integrals are obtained by power series functions in which the constants are adjusted to simulation data. However the quadrupolar (QQ) terms are given by the following equations:
aQQ = − 2
3 2 4
i
j
xi xj
5 5 εii εjj ii jj n n Q ∗2 J kB T kB T 7 Q,i Q,j j 2,ij ij
(13)
εAB /k [K]
KAB
|Q| [DA]
[D]
Reference
1.85 1.7
[18] [18] [18] [18] [18] [18] [20] [22] [22]
4.4 0.072 0.135
= aQQ 3
3
2406.700 2152.400
3 3 4
i
xi xj
j
× nQ,i nQ,j Qj∗3 J3,ij +
×
ε 3/2 ε 3/2 15/2 15/2 jj ii jj ii kB T
kB T
42 3
ij7
3 3 4
2
i
j
xi xj xk
k
5 5 5 εii εjj εkk ii jj kk n n n Q ∗2 Qj∗2 Qk∗2 J3,ijk kB T kB T kB T 3 3 3 Q,i Q,j Q,k i ij ik jk
(14)
where Qj∗2 denotes the dimensionless square quadrupole moment and nQ,j is the number of quadrupole moments in a molecule j. To use the PC-SAFT equation with dipolar and quadrupolar terms in the multiphase calculations, the component parameters are required. Table 2 lists the component parameters which are used in the computations. Parameters are taken from different sources [19,21,23]. Two association sites are assumed for water while four association sites are used for ethylene glycol. 4. Gupta et al. (1991) algorithm As discussed in Section 2, the proposed methodology of Gupta (1991) for simultaneous stability and multiphase flash calculation provided a unified set of simultaneous non-linear equations with inequality constraints. Gupta et al. [1] have proposed to partition the system of equations into two groups to be solved in two nested loops. The first group consisted of Equations (1), (2) and (5). These equations were solved for phase fractions (˛k ) and stability variables ( k ) by the Newton-Raphson method in the inner loop. The second group consisted of Equations (3) and (4) in the outer loop to calculate mole fractions. The mole fractions were updated in the outer loop and accelerated by the use of the General Dominant Eigenvalue Method (GDEM) [24]. The K-values, which were initiated with Michelsen stability analysis [12], were re-calculated from Equation (7) by the use of a chosen model in the outer loop. As previously indicated, modifications were included in their work to avoid singular and ill-conditioned Jacobian. The stability equation (Equation (5)) was replaced by the following equation: ˛k k =0 (˛k+ k )
(15)
The equation describes the equivalent behaviour of Equation (5) in the appearance and disappearance of a phase. It improves the conditioning of Jacobian significantly. The singularities of the Jacobian, when both phase fraction and its corresponding stability variable become zero, were also avoided in the way that both would never become equal to zero during computation. Furthermore, during computations, the phase fractions become negative which of course violates inequality 6. In that case Gupta et al. forced the phase fraction to be zero. The corresponding stability variable was then re-calculated from Equation (2). Moreover, during computations, some of the appearing phases may have the same compositions and compressibility (i.e. triv-
172
N.M. Alsaifi, P. Englezos / Fluid Phase Equilibria 302 (2011) 169–178
ial solution). If that occurs, Gupta et al. [1] suggested adding the phases together to decrease the number of phases and the next computations were performed only for the distinct phases.
1.0 Vapour Water-Rich Liquid Hydrocarbon-Rich Liquid
Phase Fraction
0.8
5. The proposed algorithm In this section, the proposed algorithm is given and its advantages over Gupta’s methodology are discussed. Similar to Gupta et al. approach, Equations (1), (2) and (5) are solved in the inner loop but without modifying Equation (5) and with the inclusion of the inequality constraint (Equation (6)). This requires an iterative method other than Newton-Raphson method which cannot handle inequality constraints. A recent iterative method which could solve systems of nonlinear equalities and inequalities has been proposed by Macconi et al. [25]. The method will be referred as trust-region Gauss–Newton method. It is implemented efficiently using MATLAB in a solver called TRESNEI [26]. Unlike the Newton-Raphson method, the trust-region Gauss–Newton method is capable of handling inequality constraints and has global convergence properties. The use of TRESNEI will avoid some of the tedious Gupta et al work of controlling phase fraction and re-calculated stability variable. TRESNEI would help to obtain phase fraction and stability variables from solving directly Equations (1), (2), (5) and (6). Similar to Gupta et al. [1], the mole fraction vector is selected as the iterative variable in the flash computations. However, our updating procedure differs from that of Gupta et al. [1]. The updating procedure of Kinoshita and Takamatsu [27] is utilized in our approach. We found that the Kinoshita and Takamatsu updating procedure is relatively insensitive to the initial estimates of mole fraction. The combination of trust-region Gauss–Newton method and Kinoshita and Takamatsu makes it possible to avoid Gupta’s initialization strategies [28], particularly for water-hydrocarbon containing systems. In this study, the initial estimates of mole fraction are based on specifying a value of 0.99 to any potential component that may form liquid and to the lightest component in the vapour phase. This initialization is found to be effective for water-hydrocarbon containing systems. For example, the water-pentane system might form a maximum of three phases (vapour–liquid–liquid). It is adequate in our approach to start with a mole fraction value of 0.99 for pentane in the pentane rich liquidphase, a value of 0.99 for water in the water-rich liquid phase and a value of 0.99 for the lightest component in the vapour phase. Unlike the Gupta et al. algorithm, the current algorithm does not encounter a problem when the phase fraction approaches zero. Furthermore, an important feature of the current algorithm is its ability to handle more than one absent phase (i.e. an unstable phase). The feature would assist to figure out which absent phases would appear first when the P–T conditions are changed. This feature will be illustrated in details in Section 6.1.1.
0.6
0.4
0.2
0.0 400
410
420
430
440
450
460
Temperature / K Fig. 1. Phase composition of water–paraffin mixture at 24.1 bar predicted by dipolar PC-SAFT.
It should be noted that SAFT is not cubic and the method of finding density (or volume) roots differs from that of cubic equations of state. To determine the density roots, it is necessary to use a reliable iterative method for solving the equation of state for density. It is also important to use a small termination tolerance on the function value (ε = 10−14 ). TRESNEI works effectively in finding density roots. Besides TRESNEI, we tested trust-region dogleg method, which is a variant of the Powell dogleg method [29]. This method is also successful in finding the density roots of all the cases that have been studied. There are other reported reliable methods that could be used to find the density roots of SAFT such as the extended dogleg method [30], homotopy continuation [31], and interval analysis [32]. 6. Evaluation of the proposed algorithm In this section, several examples are considered for evaluating the proposed scheme of several systems. In all studied systems, the binary interaction parameter is set to zero. Of course, the PC-SAFT may not be able to predict accurately the behaviour of many systems if binary interaction parameters are not employed. However, it is interesting to assess SAFT as a predictive tool. 6.1. Water/octane/hexane/pentane/butane/propane mixture The first considered mixture is a synthetic system of water and paraffins. The system has been studied before by Heidemann [4], Peng and Robinson [33] and Erbar [34]. The overall mole composition is made up of 0.1667 propane, 0.1667 n-butane, 0.2
Table 3 Phase fractions and stability variables of water–paraffin mixture at 24.1 bar. T (K)
460.000a 454.513a 450.000a 441.531a 420.000a 406.981a 390.000b a b c
Water-rich liquid phase (LW )
Hydrocarbon-rich liquid phase (LHC )
Vapour phase (V)
˛LW
LW
˛LHC
LHC
˛V
V
0.0000 0.0000 0.0000 0.0000 0.1948 0.2418 0.2517
0.4816 0.3629 0.2382 0.0000 0.0000 0.0000 –
0.0000 0.0000 0.0483 0.1356 0.5163 0.7582 0.7483
0.0566 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
1.0000 1.0000 0.9517 0.8644 0.2889 0.0000 0.0000
– – – – – – 0.1849
Vapour is the reference phase. Water-rich liquid phase is the reference phase. Intel Pentium Dual E2200 2.2 GHz.
Mixture state
Iter #
CPUc time (s)
Vapour only LHC incipient Vapour-LHC LW incipient V-LW -LHC LW -LHC LW -LHC
13 13 13 13 13 12 12
10 11 10 10 9 17 7
6.1.1. One phase region (vapour only) At a temperature equal to 460 K, the vapour phase fraction is equal to 1 and the other two phase fractions of liquid phases are zero. This indicates that the two liquid phases are absent and we only have a vapour phase. Although one could delete a phase from the three phases during computation when the system shows only one phase (here vapour phase), the computation is conducted without deleting any phase. The stability variables of the two liquid phases at this temperature are positive (0.4816 for water-rich liquid phase and 0.05660 for hydrocarbon-rich liquid phase). The stability variable value of the hydrocarbon-rich liquid phase is much less than that of the water-rich liquid phase. This indicates that decreasing the temperature makes the hydrocarbon-rich liquid phase appear first. The change of phases and stability variables vs. iteration number are depicted in Fig. 2. The figure illustrates the effectiveness of the proposed scheme with the trust-region Gauss–Newton method in providing fast convergence. Computation is started with zero values for stability variables and phase fractions. As iteration proceeds, the phase fraction of the vapour (˛V ) and hydrocarbon-rich liquid (˛LHC ) appear first while the water-rich liquid phase (˛LW ) remains zero and its stability variable (LW ) becomes positive. At the third iteration, the hydrocarbon-rich liquid phase disappears and its stability variable (LHC ) attains a positive value. Both liquid stability variables converge to positive values while their phase fractions become zero indicating that they are unstable phases. The smaller magnitude of the stability variable of hydrocarbon-rich liquid phase confirms the behaviour observed before in Fig. 1 in which hydrocarbon-rich liquid phase would appear first if temperature is decreased. 6.1.2. Incipient hydrocarbon-rich liquid and vapour–liquid phase region As the temperature is decreased, the hydrocarbon-rich liquid phase is about to appear at 454.513 K where both its phase fraction and stability variable are zero as shown in Table 3. This demonstrates that at this condition (24.1 bar and 454.513 K) the hydrocarbon-rich liquid phase is incipient and a slight decrease in temperature brings the system to the two phase region; namely
1.0 0.8 0.6 0.4 0.2 0.0
(a)
0.6
(b)
173
αV
0.5 0.4
αL
0.3 0.2
θL
0.1 0.0 0.6
W
W
(c)
0.5
αL
0.4 0.3
θL
0.2 0.1
HC
HC
0.0 0
2
4
6
8
10
12
14
Iteration Number Fig. 2. Phase fraction (˛) and stability variable () vs. iteration number for water–paraffin mixture at 24.1 bar and 460 K for: (a) vapour (stable); (b) water-rich liquid phase (unstable); (c) hydrocarbon-rich liquid phase (unstable).
vapour and hydrocarbon-rich liquid phases. An example of the vapour and hydrocarbon-rich liquid phase is at a temperature of 450 K. The phase fraction of the water-rich liquid phase is still zero at this temperature indicating that it is absent and its stability variable is still positive (0.2382). Therefore, the water-rich liquid phase is still unstable at this condition but its stability value gets smaller. More insight into the behaviour of the phase fraction and stability variables during the iteration is shown in Figs. 3 and 4 at 454.513 and 450 K; respectively. The iteration represents the number of times the composition vectors are updated. As shown in Fig. 3, the phase fraction and stability variable converge to zero for hydrocarbon-rich liquid phase at the 5th iteration. In Fig. 4, the phase fraction of hydrocarbon-rich liquid phase starts develop-
Phase Fraction or Stability Variable
n-pentane, 0.0667 hexane, 0.1333 n-octane and 0.2667 water. The PC-SAFT with the proposed methodology is applied to study this mixture at fixed pressure (P = 24.1 bar) and a temperature range of 390–460 K. Both association and dipolar terms are included to take account of hydrogen bonding and dipole–dipole interactions of water molecules, respectively. The condensation curves are obtained and shown in Fig. 1. The general characteristic of the curves are similar to those obtained before [4,33,34] although different equations of state were utilized. As seen in Fig. 1, when the mixture is cooled from its vapour state, the formation of a hydrocarbon-rich liquid phase appears first. To demonstrate the proposed scheme, it is advantageous to show the computational details in this example. It should be noted that the system may have up to three phases in the selected temperature range (vapour only, vapour–liquid, liquid–liquid and vapour–liquid–liquid). An isobaric (24.1 bar) temperature path is followed by starting from 460 K (vapour only) to 390 K (liquid–liquid). We assume that the system consists of a maximum of three phases and the vapour is the reference phase. The reference phase could be switched to one of the liquid phases when the system does not show vapour phase as demonstrated below. Table 3 shows values for the phase fractions and stability variables for this mixture for seven temperatures from 390 to 460 K at 24.1 bar. The initial values of phase fractions and stability variables are arbitrary and set to zeros in this example. In the following, computation details of Table 3 are demonstrated for each region.
Phase Fraction or Stability Variable
N.M. Alsaifi, P. Englezos / Fluid Phase Equilibria 302 (2011) 169–178
1.0 0.8 0.6 0.4 0.2 0.0
(a)
αV
0.5
(b)
0.4
αL
0.3 0.2
θL
0.1
W
W
0.0 0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00
(c)
αL θL
0
2
4
6
8
10
12
HC
HC
14
Iteration Number Fig. 3. Phase fraction (˛) and stability variable () vs. iteration number for water–paraffin mixture at 24.1 bar and 454.513 K for: (a) vapour (stable); (b) waterrich liquid phase (unstable); (c) hydrocarbon-rich liquid phase (incipient).
N.M. Alsaifi, P. Englezos / Fluid Phase Equilibria 302 (2011) 169–178
1.0 0.8 0.6 0.4 0.2 0.0
(a)
αV
0.35 0.30 0.25 0.20 0.15 0.10 0.05 0.00
Phase Fraction or Stability Variable
Phase Fraction or Stability Variable
174
(b)
αL θL
0.18 0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00
(c)
W
W
αL
HC
θL
1.0 0.8 0.6 0.4 0.2 0.0
(a)
0.25
(b)
0.20
αV
0.15
αL
0.10
θL
0.05
W
0.00 0.6
(c)
0.5
αL
0.4
HC
0.3
HC
W
θL
0.2
HC
0.1 0.0 0
2
4
6
8
Iteration Number
10
12
14
0
2
4
6
8
10
12
14
Iteration Number
Fig. 4. Phase fraction (˛) and stability variable () vs. iteration number for water–paraffin mixture at 24.1 bar and 450 K for: (a) vapour (stable); (b) water-rich liquid phase (unstable); (c) hydrocarbon-rich liquid phase (stable).
Fig. 6. Phase fraction (˛) and stability variable () vs. iteration number for water–paraffin mixture at 24.1 bar and 420 K for: (a) vapour (stable); (b) water-rich liquid phase (stable); (c) hydrocarbon-rich liquid phase (stable).
ing after the 2nd iteration where its stability variable approaches zero and then the phase becomes present at the third iteration. The appearance of the hydrocarbon liquid phase is accompanied by a decrease in the vapour phase fraction which converges to 0.9517.
in the vapour–liquid–liquid region. An illustration of the characteristics of the vapour–liquid–liquid region is given in Table 3 and Fig. 6 at 420 K. As seen in Table 3, the liquid phase fractions are positive and their corresponding stability variables are zero indicating that they are present. In Fig. 6, the change of the stability variable and phase fraction with iteration is depicted. Both stability variables of the two liquid phases are zero while their phase fractions are positive showing that they are stable.
1.0 0.8 0.6 0.4 0.2 0.0
(a)
0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00
(b)
0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00
(c)
αV
αL θL
6.1.4. Liquid–liquid phase If we continue cooling the mixture, the vapour phase, which is the reference phase, disappears at 406.981 K indicating that the mixture exhibits two liquid phases; namely, water-rich and
Phase Fraction or Stability Variable
Phase Fraction or Stability Variable
6.1.3. Incipient water-rich liquid and vapour–liquid–liquid phase region The water-rich liquid phase is incipient at temperature equals to 441.531 K as shown in Table 3 where both stability variable and phase fraction have zero values. Both phase fraction and stability variable are zero (i.e. incipient) after the 3rd iteration as shown in Fig. 5. If the temperature is slightly decreased, the mixture coexists
W
W
αL
HC
θL
HC
1.0 0.8 0.6 0.4 0.2 0.0
(a)
0.30
(b)
0.25
αV
0.20
αL
0.15
θL
0.10 0.05
W
W
0.00 1.0
(c)
0.8 0.6
αL
0.4
θL
HC
HC
0.2 0.0
0
2
4
6
8
10
12
14
Iteration Number Fig. 5. Phase fraction (˛) and stability variable () vs. iteration number for water–paraffin mixture at 24.1 bar and 441.53 K for: (a) vapour (stable); (b) waterrich liquid phase (incipient); (c) hydrocarbon-rich liquid phase (stable).
0
2
4
6
8
10
12
14
Iteration Number Fig. 7. Phase fraction (˛) and stability variable () vs. iteration number for water–paraffin mixture at 24.1 bar and 406.981 K for: (a) vapour (incipient); (b) water-rich liquid phase (stable); (c) hydrocarbon-rich liquid phase (stable).
N.M. Alsaifi, P. Englezos / Fluid Phase Equilibria 302 (2011) 169–178
0.4
(a)
0.3
Phase Fraction or Stability Variable
hydrocarbon-rich liquid phases. At this point, both the two liquid phase fractions are positive and their corresponding stability variables are zero whereas the vapour phase fraction is zero. A fast convergence of phase fraction and the stability variables is achieved as shown in Fig. 7. It should be emphasized here, as previously indicated, that the proposed scheme does not encounter any computation difficulty when the phase fraction approaches zero even if that phase is the reference phase. The disappearance of the vapour phase at temperature equal or less than 406.981 K and 24.1 bar makes it necessary to switch the reference phase to one of the liquid phases if computation is required for temperature less than 406.981 K. If the reference phase is maintained to be the vapour phase, a decrease in temperature would cause the vapour phase fraction to be negative. A negative value is of course not a physical value. However, it is an indication of improper selection of the reference phase. The reference phase is switched to the water-rich liquid phase for the flash calculation at a temperature less than 406.981 K where the system exhibits liquid–liquid equilibrium. For example, at 390 K, the vapour phase is absent and its stability variable is positive (0.1849) as shown in Table 3. The characteristics of the phase fractions and stability variables with iterations are shown in Fig. 8. It is clear that the phase fraction of vapour is zero and its corresponding stability variable (see Fig. 8) becomes positive at the 2nd iteration. Finally, it should be noted that Equation (8) is verified in all the studied conditions.
175
0.2
αL
0.1 0.0 0.20
W
(b)
0.15
αV
0.10
θV
0.05 0.00 0.8
(c)
0.6
αL
0.4
θL
0.2
HC
HC
0.0 0
2
4
6
8
10
12
14
Iteration Number Fig. 8. Phase fraction (˛) and stability variable () vs. iteration number for water–paraffin mixture at 24.1 bar and 390 K for: (a) vapour (unstable); (b) waterrich liquid phase (stable); (c) hydrocarbon-rich liquid phase (stable).
0.7
The second mixture that was examined is the ethylene glycol–water–propane–methane quaternary. Dipole–dipole interactions are considered for both ethylene glycol and water. The dipole–dipole interactions are represented with two polar functional groups of the ethylene glycol. Hydrogen bonding interactions are also taken into account for water and ethylene glycol. Experimental data was reported for this mixture by Ng et al. [35] at two different P–T values. The system shows vapour–liquid–liquid behaviour according to the reported experimental P–T values. A comparison of the experimental and the calculated mole fraction of this system with the proposed scheme is shown in Table 4. As seen in the table, dipolar PC-SAFT predicts the mole fractions of this system well. The proposed scheme is also applied for this mixture at a range of pressures at a fixed temperature (283.15 K) to show the change of phase fractions of the three phases. The phase fraction vs. pressure is depicted in Fig. 9. As seen in the figure, as the pressure increases the mole phase fraction of vapour decreases while the phase fraction of hydrocarbon-rich liquid phase increases. How-
0.6
Mole Phase Fraction
6.2. Ethylene glycol/water/propane/methane
Vapour Water-Rich Liquid Hydrocarbon-Rich Liquid
0.5 0.4 0.3 0.2 0.1 0.0 10
20
30
40
50
60
70
80
90
Pressure / bar Fig. 9. Mole phase fraction vs. pressure of methane–propane–water–ethylene glycol at 283.15 K predicted by dipolar PC-SAFT.
Table 4 Equilibrium phase composition (mole fraction) of methane–propane–water–ethylene glycol system. T = 273.15 K
MEG Water Propane Methane
V (exp.)
Feed
V
0.092 0.317 0.284 0.307
3.620E−07 1.207E−04 0.1566 0.8433
T = 283.15 K
MEG Water Propane Methane
Feed
V
0.090 0.311 0.288 0.311
8.746E−07 2.266E−04 0.1961 0.8037
LHC
P = 70.4 bar LW 0.2236 0.7702 2.562E−03 3.719E−03
– – 0.1542 0.8458
2.204E−06 1.277E−04 0.5394 0.4605
3.400E−06 4.000E−05 0.5660 0.4340
V (exp.)
LHC
P = 69.0 bar
– – 0.1984 0.8016
4.497E−06 2.211E−04 0.5834 0.4163
LW (exp.)
LHC (exp.)
0.2243 0.7736 2.770E−04 1.790E−03 LW (exp.)
LHC (exp.)
LW
3.000E−06 7.200E−05 0.6018 0.3982
0.2235 0.7700 2.928E−03 3.583E−03
0.2244 0.7736 3.760E−04 1.580E−03
176
N.M. Alsaifi, P. Englezos / Fluid Phase Equilibria 302 (2011) 169–178
ever, the phase fraction of the water-rich liquid phase is incipient at temperature remains constant as the pressure increases. There is no experimental data to verify this behaviour but prediction of similar systems (No. 10, 11 and 12 in Table 6) which have been studied and compared to experimental at different pressures supports this behaviour. The proposed scheme performs the calculation of this system without any difficulty although the mole phase fraction was approaching zero for hydrocarbon-rich liquid phase around 15 bar and for vapour phase around 84 bar. In this example, the computations were performed for the vapour–liquid–liquid region but can be extended to other phase regions as done with water-paraffins system. The vapour phase was selected as the reference phase in this example. Furthermore, the number of iterations did not exceed 18 for each selected P–T values while the total CPU time is about 60 s.
0.4
αV
Phase Fraction or Stability Variable
0.3
6.3. Ethylene glycol/water/hexane/propane/CO2 /methane
0.2 0.1 0.0 0.6 0.5
αL
0.4 0.3
θL
0.2 0.1
6.4. Water-hydrocarbon and other systems The proposed algorithm was successfully tested for several water-hydrocarbon and other systems listed in Table 6. The listed systems exhibit various phase behaviours including vapour–liquid, liquid–liquid and vapour–liquid–liquid equilibrium. For the listed systems, two SAFT versions [19,36] have been utilized with the addition of hydrogen bonding and dipolar interactions. The reader
W
0.0 0.5 0.4 0.3
αL
0.2
θL
0.1
The experimental data of this six-component mixture were reported by Ng et al. [35] at two different temperatures (263.15 and 283.15 K) and a fixed pressure of 68 bar. The reported data show that this system exhibits vapour–liquid–liquid equilibrium at these conditions. To examine this system with the proposed scheme, the PC-SAFT is utilized with the inclusion of dipolar and hydrogen bonding interactions of ethylene glycol and water as well as quadrupole interactions for carbon dioxide. Fig. 10 represents the phase fraction and its corresponding stability variables vs. iteration. The stability variables are arbitrary initiated with non-zero values. Nonetheless, a fast convergence is achieved and Equations (5) and (6) are verified. At the 2nd iteration, the stability variables go to zero where the phase fractions of both liquids appear and proceed to converge at 10 iterations. The total CPU time is around 43 s. Table 5 illustrates the comparison between experimental and calculated mole fractions. It is apparent that the dipolar PC-SAFT compare well the experimental data although no adjustable parameter is introduced.
W
HC
HC
0.0 0
2
4
6
8
10
12
Iteration Number Fig. 10. Phase fraction (˛) and stability variable () vs. iteration number for ethylene glycol–water–hexane–carbon dioxide–methane mixture at 68 bar and 283.15 K for: (a) vapour (stable); (b) water-rich liquid phase (stable); (c) hydrocarbon-rich liquid phase (stable).
could refer to the PhD Dossier of Alsaifi [37] for the details of calculations. In this work, the binary interaction parameter was set to zero to assess the predictive capability of PC-SAFT in multiphase systems. The prediction was very satisfactory of systems of water–ethylene glycol–hydrocarbons based on the study of systems in Sections 6.2 and 6.3 and systems 10, 11 and 12 in Table 6. Finally, although we have employed PC-SAFT with the addition of association, dipolar and quadrupole terms, the total CPU time was not huge even for cross-association systems (e.g. water–ethylene glycol systems). As far as our knowledge goes, the only available study for computing phase stability and equilibrium using SAFT for cross association systems is the one that conducted by Xu et al. [32]. They described a reliable method to compute phase stability and equilibrium using interval analysis. Their work was limited to mainly VLE and binary systems. The study showed that the computing time for conducting stability and equilibrium calculations for cross-association systems was large in comparison to ours. For example, the total CPU time of 1-butanol-ethanol sys-
Table 5 Equilibrium phase composition (mole fraction) for ethylene glycol–water–hexane–propane–carbon dioxide–methane. T = 263.15 K
EG H2 O C6 H14 C3 H8 CO2 CH4
V (exp.)
Feed
V
0.114 0.392 0.198 0.049 0.049 0.198
1.120E−07 6.232E−05 0.0034 0.0180 0.0832 0.8953
T = 283.15 K
EG H2 O C6 H14 C3 H8 CO2 CH4
Feed
V
0.115 0.397 0.195 0.049 0.049 0.195
6.899E−07 2.270E−04 0.0066 0.0275 0.1044 0.8613
– – 0.0041 0.0217 0.1199 0.8543
LHC
2.175E−06 9.172E−05 0.4183 0.1025 0.0992 0.3800
P = 68 bar
LW (exp.)
LHC (exp.)
LW
5.400E−06 4.400E−05 0.4500 0.1100 0.0948 0.3452
0.2241 0.7704 2.974E−04 3.982E−04 0.00114 0.00373
V (exp.)
LHC
P = 68 bar LHC (exp.)
LW
– – 0.0071 0.0297 0.1306 0.8326
8.737E−06 2.725E−04 0.4640 0.1120 0.0989 0.3248
1.200E−05 9.600E−05 0.4800 0.1154 0.0870 0.3176
0.2233 0.7708 4.019E−04 5.022E−04 1.28E−03 3.68E−03
0.2235 0.7705 2.00E−05 3.00E−05 0.00448 0.00148 LW (exp.)
0.2234 0.7704 2.600E−06 4.800E−05 4.45E−03 1.67E−03
N.M. Alsaifi, P. Englezos / Fluid Phase Equilibria 302 (2011) 169–178
177
Table 6 Other studied systems.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
System
Model
Ethane–propane–butane–pentane–hexane (VLE) Methane–ethane–propane–pentane–heptane–decane (VLE) Methane–water–methylcyclohexane (VLLE) Methane–water–heptane (VLLE) Methane–water–tert-butyl methyl ether (VLLE) Methane–water–neohexane (VLLE) Methane–water–decane (VLLE) p-xylene–water–ethanol (VLLE) Propane–water–pentane–octane (VLLE) Methane–water–ethylene glycol–methylcyclohexane (VLLE) Methane–water–ethylene glycol–hexane (VLLE) Methane–water–ethylene glycol–octane (VLLE) Methanol–hexane (LLE) Methanol–heptane (LLE) Methanol–decane (LLE) Ethanol–tetradecane (LLE)
PC-SAFT PC-SAFT Dipolar PC-SAFT Dipolar PC-SAFT Dipolar PC-SAFT Dipolar PC-SAFT Dipolar PC-SAFT Dipolar PC-SAFT Dipolar PC-SAFT Dipolar PC-SAFT Dipolar PC-SAFT Dipolar PC-SAFT Dipolar Simplified-SAFT [36] Dipolar Simplified-SAFT [36] Dipolar Simplified-SAFT [36] Dipolar Simplified-SAFT [36]
tem was 1816.1 s at T = 343 K, P = 0.35 bar and zbutanol = 0.65 (feed composition) using the original SAFT [38] with a binary interaction parameter equals to 0.011. For the sake of comparison, we studied this system at the same conditions and the same SAFT version [38]. The total CPU time is found to be about 30 s which is significantly faster than that of Xu et al. [32]. The phase fraction of vapour, which is the reference phase, is 0.97082 while the liquid phase fraction is 0.029177 and its corresponding stability variable is 1.5955e−22 indicating that it is stable. Our algorithm shows a third and absent phase with stability variable equals to 2.4917. 7. Conclusions The procedure suggested originally by Gupta (1990) to perform multiphase equilibrium and stability calculations was modified. The new scheme combines the trust-region Gauss–Newton method [25,26] and the updating procedure of Kinoshita and Takamatsu [27]. Unlike the classical iterative method, the trust-region Gauss–Newton method is capable of solving non-linear equations with inequality constrains and it has global convergence properties. The computational scheme utilizes a simple initialization procedure. The current algorithm does not encounter problems when the phase fraction approaches zero. Furthermore, it could handle more than one absent phase easily. The proposed computation scheme was utilized and illustrated the predictive capability of the PC-SASFT equation of state for water–ethylene glycol–hydrocarbon systems with or without carbon dioxide. Hydrogen bonding and polar (dipolar and quadrupole) interactions were included in the PC-SAFT and computations were performed without the use of binary interaction parameters. The scheme is proved to be efficient for multiphase calculations for water-hydrocarbon containing systems. List of symbols a reduced Helmholtz free energy f fugacity J2,ij integral over the reference-fluid pair-correlation function J3,ijk integral over three-body correlation function k phase Boltzmann constant (J/K) kB kij binary interaction parameter m segment number nc number of components ni number of dipole moments in component i Q* dimensionless quadrupole moment r reference phase T absolute temperature (K)
x y z
liquid mole fraction vapour mole fraction feed composition
Greek letters ˛ phase fraction stability variable ϕ fugacity coefficient εii segment energy parameter of component i dipole moment i dipole moment for component i ∗2 dimensionless squared dipole moment i number of phases number density ii segment size parameter of component i (Å) Acknowledgments We are grateful to the King Fahd University of Petroleum & Minerals for its support of Nayef Alsaifi in funding his graduate scholarship. Financial support from the Natural Sciences and Engineering Research Council of Canada (NSERC) is also acknowledged. The authors wish also to give their sincere thanks to Benedetta Morini and Margherita Porcelli for making TRESNEI code available for public use. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21]
A.K. Gupta, P.R. Bishnoi, N. Kalogerakis, Fluid Phase Equilibr. 63 (1991) 65–89. A.K. Gupta, P.R. Bishnoi, N. Kalogerakis, Gas Sep. Purif. 4 (1990) 215–222. G. Han, G.P. Rangaiah, Comput. Chem. Eng. 22 (1998) 897–911. R.A. Heidemann, AIChE J. 20 (1974) 847–855. R. Gautam, W.D. Seider, AIChE J. 25 (1979) 991–999. M.E. Soares, A.G. Medina, C. McDermott, N. Ashton, Chem. Engng. Sci. 37 (1982) 521–528. M. Castier, P. Rasmussen, A. Fredenslund, Chem. Engng. Sci. 44 (1989) 237–248. W.X. Song, M.A. Larson, Chem. Engng. Sci. 46 (1991) 2513–2523. H. Pan, A. Firoozabadi, SPE Res. Eval. Engng. Feb. (1998) 36–42. A. Lucia, L. Padmanabhan, S. Venkataraman, Comput. Chem. Eng. 24 (2000) 2557–2569. L.E. Baker, A.C. Pierce, K.D. Luks, SPE/DOE Second Joint Symp. Enhanced Oil Recovery, Tulsa, Oklahoma, 1981. M.L. Michelsen, Fluid Phase Equilib. 9 (1982) 1–20. G.A. Iglesias-Silva, A. Bonilla-Petriciolet, P.T. Eubank, J.C. Holste, K.R. Hall, Fluid Phase Equilib. 210 (2003) 229–245. M.L. Michelsen, Fluid Phase Equilib. 9 (1982) 21–35. L.X. Nghiem, Y. Li, Fluid Phase Equilib. 17 (1984) 77–95. J.-S. Wu, P.R. Bishnoi, Comput. Chem. Eng. 10 (1986) 269–276. P.R. Bishnoi, A.K. Gupta, P. Englezos, N.E. Kalogerakis, Fluid Phase Equilib. 53 (1989) 97–104. R.M. Abdel-Ghani, Master Thesis, University of Calgary, Canada, 1995. J. Gross, G. Sadowski, Ind. Eng. Chem. Res. 41 (2002) 5510–5515. J. Gross, J. Vrabec, AIChE J. 52 (2006) 1194–1204. J. Gross, AIChE J. 51 (2005) 2556–2568.
178 [22] [23] [24] [25] [26]
N.M. Alsaifi, P. Englezos / Fluid Phase Equilibria 302 (2011) 169–178
G.S. Rushbrook, G. Stell, J.S. Hoye, Mol. Phys. 26 (1973) 1199–1215. N.M. Al-Saifi, E.Z. Hamad, P. Englezos, Fluid Phase Equilib. 271 (2008) 82–93. C.M. Crowe, M. Nishio, AIChE J. 21 (1975) 528–533. M. Macconi, B. Morini, M. Porcelli, Appl. Numer. Math. 59 (2009) 859–876. B. Morini, M. Porcelli, Comput. Optim. Appl., (2010) DOI: 10.1007/s10589-010rr9327-5. [27] M. Kinoshita, T. Takamatsu, Comput. Chem. Eng. 10 (1986) 353–360. [28] A.K. Gupta, Ph.D. Thesis, Dept. of Chem. and Petr. Eng., University of Calgary, Calgary, Canada, 1990. [29] P. Rabinowitz, Numerical Methods for Nonlinear Algebraic Equations, Ch.7, 1970.
[30] N. Aslam, A.K. Sunol, Ind. Eng. Chem. Res. 45 (2006) 3303–3310. [31] A. Lucia, Q. Luo, Adv. Environ. Res. 6 (2001) 123–134. [32] G. Xu, J.F. Brennecke, M.A. Stadtherr, Ind. Eng. Chem. Res. 41 (2002) 938–952. [33] D.-Y. Peng, D.B. Robinson, Canad. J. Chem. Eng. 54 (1976) 595–599. [34] J.H. Erbar, Proc. 52nd Ann. National. Gas Processors Assoc. Meeting, 62, 1973. [35] H.-J. Ng, C.-J. Chen, D.B. Robinson, GPA Research Report RR-92, 1986. [36] F. Yuan-H, S.I. Sandler, Ind. Eng. Chem. Res. 34 (1995) 1897–1909. [37] N.M. Al-Saifi, PhD Dossier, Department of Chem. & Bio Eng., UBC, 2009. [38] S.H. Huang, M. Radosz, Ind. Eng. Chem. Res. 30 (1991) 1994–2005.
Contents lists available at ScienceDirect
Fluid Phase Equilibria journal homepage: www.elsevier.com/locate/fluid
Prediction of multiphase equilibrium using the PC-SAFT equation of state and simultaneous testing of phase stability Nayef M. Alsaifi, Peter Englezos ∗ Department of Chemical and Biological Engineering, University of British Columbia, Vancouver, BC, Canada V6T 1Z3
a r t i c l e
i n f o
Article history: Received 29 May 2010 Received in revised form 1 September 2010 Accepted 3 September 2010 Available online 15 September 2010 Keywords: Multiphase calculations Stability analysis PC-SAFT Association Dipolar Quadrupole Hydrocarbon Water Ethylene glycol Trust-region Gauss–Newton method
a b s t r a c t This work employs the PC-SAFT equation of state to illustrate how to compute phase equilibria and stability of multiphase systems. The proposed computation scheme is based on the methodology of Gupta et al. (1991) but with modifications. Phase equilibria and stability are computed simultaneously and without the need to modify the exact stability equation. This approach does not encounter difficulty when the phase fraction approaches zero. A simple initialization strategy is utilized. In addition, this work tests the predictive capability of PC-SAFT in multiphase systems of water-hydrocarbon containing systems at different pressure–temperature conditions. The dipole and quadrupole interactions were included in the PC-SAFT equation and the binary interaction parameters (kij ) were set to zero. The work has demonstrated the ability of PC-SAFT to successfully predict vapour–liquid–liquid equilibrium of systems containing water, hydrocarbons and ethylene glycol with or without carbon dioxide. © 2010 Elsevier B.V. All rights reserved.
1. Introduction Computation of equilibrium and stability of multiphase systems plays an essential role in simulation and design of chemical processes such as distillation and extraction. The requirement for efficient and reliable multiphase computation methods is always of interest in terms of the effectiveness and reliability of process simulators. These methods are generally based on either the Gibbs free energy minimization approach or equation-solving techniques. The equation-solving techniques combine mass balances and equilibrium relations into a set of non-linear algebraic equations. Different approaches have been suggested for equation-solving techniques [1–3]. The Gibbs energy minimization approach, however, formulates the phase equilibrium problem into an optimization problem which can be solved by a global solver subject to material balance constraints. Various methodologies have been developed under this method by many researchers [4–10]. Irrespective of which method is adopted, the essence of phase equilibrium calculations is always to ensure the determination of the global minimum Gibbs free energy of a system at fixed isothermal and isobaric conditions.
∗ Corresponding author. Tel.: +1 604 822 6184; fax: +1 604 822 6003. E-mail address: [email protected] (P. Englezos). 0378-3812/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.fluid.2010.09.002
The main challenge in multiphase equilibrium calculations is that the number of the available phases in a system at equilibrium is not known a priori. A common way to overcome the encountered phase difficulty is to utilize stability criteria [11,12] to determine whether an obtained solution is thermodynamically stable. This can be handled either by assuming a fixed number of phases and checking the stability of the final solution [13] or by following a sequential procedure by adding a phase in the computation and testing the stability [5,14–16]. A novel development of stability and multiphase calculations which has some advantages over the previous approaches is the one proposed by Gupta et al. [1,2,17]. The methodology formulates the problem as a set of coupled non-linear equations that describe phase equilibrium, chemical equilibrium and the stability of the system with inequality constraints. Thus, it is an equation-solving technique type that combines mass balance and equilibrium relations as well as stability equations with inequality constraints. The advantages of the Gupta methodology are its simplicity and its ability to provide insight into the phase behaviour as it tracks the appearance or disappearance of the phases. In principle, the method presents the necessary and sufficient conditions for phase equilibrium and it offers a simple and clear stability criterion for multiphase reacting/non-reacting systems at a fixed pressure and temperature. Gupta et al. [1] modified the exact stability equation and replaced it with one that satisfies equivalently the behaviour of
170
N.M. Alsaifi, P. Englezos / Fluid Phase Equilibria 302 (2011) 169–178
appearance and disappearance of a phase. The reason was to avoid the problem of ill-conditioned Jacobian when the system of the non-linear equations is required to be solved simultaneously. The treatment of singular Jacobian was also taken into account with another modification. Furthermore, because Gupta et al. development requires two inequality constraints for each additional phase (stability variable and phase fraction, Equation (6)) to be greater than or equal to zero, it was necessary to readjust them during iterations when the two inequality constraints were violated. This is because the Newton-Raphson method cannot handle constraints when solving a system of non-linear equations. Therefore, the inequality constraints were implicitly included during the computation. In order to guarantee success of the previous modifications, a good initialization strategy of K-values was crucial. The K-values were initiated with the stability method suggested by Michelsen [12]. Although their procedure was successful for several systems, it encounters difficulty when the phase fraction approaches zero [18]. In this study, a different scheme is suggested to perform the simultaneous calculations. The suggested scheme does not encounter a difficulty when the phase fraction approaches zero. Furthermore, it is applied without modifying the stability equation. When the new scheme is applied to water-hydrocarbon containing systems, a simple initialization strategy is followed. Finally, when the system of interest could show more than two phases on the P–T continuum (say vapour–liquid–liquid) and the selected P–T condition is in one phase region (say Vapour phase), then the proposed scheme would easily assist in figuring out which liquid phase would appear first provided that the utilized model is able to qualitatively capture the phase behaviour of the studied system. For the scope of this paper, we restrict ourselves to single-stage multiphase equilibrium calculations. Chemical equilibrium calculations are not included. Moreover, due to the ongoing development of very sophisticated models and the requirement of evaluating these models for multiphase calculations, we are not going to utilize simple models such as cubic equations of state. Instead, perturbed-chain statistical association fluid theory (PC-SAFT) [19] is utilized including association, dipolar and quadrupolar terms. For the dipolar and quadrupolar terms, developments by Gross’ group are utilized in this work [20,21]. Practical systems such as waterhydrocarbon containing systems with ethylene glycol and carbon dioxide are tested and compared to experimental data. Binary interaction parameters are set to zero in all the calculations of this paper.
2. Multiphase and stability equations The appropriate set of equations which describe the stability and the isothermal–isobaric flash calculations of reacting and nonreacting systems are formulated by Gupta et al [1,2]. The details of their derivation are not repeated here. The reader could consult Gupta et al. [1] for more details about the basis of their development. The main governing equations of non-reacting systems at constant pressure and temperature are summarized here as follows.
Table 1 State of phase k in equations (5) and (6).
Case I Case II Case III
Stability variable ( k )
Phase fraction (˛k )
State of phase k
Zero Positive Zero
Positive Zero Zero
Stable Unstable Incipient
2.2. Mole fraction summation N
(kik ek − 1)xir = 0
(k = 1, . . . , ; k = / r)
(2)
i=1
xir =
nc
z
z )(1 + i=1 i
(
i
j j = 1 (kij e − 1)˛j ) j= / r
(i = 1, . . . , nc; j = 1, . . . , ; j = / r)
(3)
2.3. Phase equilibrium and stability relations xik = Kik xir ek ˛k k = 0
(i = 1, . . . , nc; k = 1, . . . , ; k = / r)
(k = 1, . . . , ; k = / r)
(4) (5)
Subject to ˛k ≥ 0
and
k ≥ 0
where ϕ kik = ir ϕik
(6)
(7)
In the above equations, z is feed composition, x is mole fraction, k is the stability variable, nc is the number of components and ϕik is the fugacity coefficient of component i in phase k. Furthermore, subscript r denotes a reference phase. The number of equations required to describe the phase equilibrium and stability of a specific system at a fixed temperature and pressure is equal to [(nc+2) − 1] for [(nc+2) − 1] unknowns, namely () ˛k , ( − 1) k and (nc·)xik . It should be emphasized that Equations (5) and (6) imply that either the phase fraction (˛k ) or the stability variable ( k ) or both must be zero. In case the phase fraction becomes zero, this indicates that the stability variable must be either a positive for an unstable phase or zero for an incipient phase. In a similar manner, if the stability variable becomes zero, then the phase fraction must be positive for a stable phase or zero for an incipient phase. These possible cases for phase (k) are summarized in Table 1. It should be noted that Equations (5) and (6) give the relationship between stability variables and phase fraction at the global minimum of Gibbs energy [2]. Finally, as demonstrated by Gupta et al. [1], at the point of minimum Gibbs free energy, the stability variables must satisfy the following relation which is written in terms of fugacity (fik ): k = ln
f ik
fir
(i = 1, . . . , nc; k = 1, . . . , ; k = / r)
(8)
3. Equation of state The SAFT model is one of the most commonly used theory based equations of state. The general framework of SAFT is written in terms of reduced Helmholtz energy as follows:
2.1. Phase fraction summation
˛k = 1
k=1
where is the number of phases and ˛k is the phase fraction.
(1)
ares = ahs + achain + adisp + aassoc + apolar
(9)
where res, hs, chain, disp, assoc and polar stand for residual, hard sphere, chain, dispersion, association and polar, respectively. The polar term represents the long-ranged electrostatic
N.M. Alsaifi, P. Englezos / Fluid Phase Equilibria 302 (2011) 169–178
171
Table 2 PC-SAFT parameters. Compound
M [g/mol]
m
[Å]
ε/k [K]
Methane Propane Butane Pentane Hexane Octane Carbon dioxide Water Ethylene glycol
16.043 44.096 58.123 72.146 86.177 114.231 44.010 18.015 62.068
1.000 2.002 2.332 2.690 3.058 3.818 1.513 1.226 3.404
3.704 3.618 3.709 3.773 3.798 3.837 3.187 2.792 2.844
150.030 208.110 222.880 231.200 236.770 242.780 163.330 203.000 178.370
interactions such as dipole–dipole (DD), dipole–quadrupole (DQ) and quadrupole–quadrupole (QQ) interactions. In this work, the dipole–dipole interactions are added for water and alcohols to increase the predictive capability of PC-SAFT as demonstrated in a previous study [23]. Furthermore, the quadrupole–quadrupole interactions are taken into account for carbon dioxide in order to improve the predictive capability of SAFT since the quadrupole–quadrupole interactions are major interactions in carbon dioxide. In this work, we mainly utilize perturbed-chain statistical association fluid theory (PC-SAFT) equation of state although other SAFT versions could be used. The simplified SAFT [36] is also used for testing some systems [Table 6]. PC-SAFT has been extensively described in the literature. The main terms of PCSAFT will not be repeated here. The reader could refer to the work of Gross and Sadowski [19] for details. Only the dipolar and quadrupolar terms will be given [20,21]. The dipolar and quadrupole terms are both developed based on the Pade approximate form of Rushbrooke et al. [22]:
apolar =
polar a2 . polar polar 1 − (a3 /a2 )
(10)
where polar superscript denotes dipolar (DD) or quadrupole (QQ) terms. The dipolar terms (DD) are defined as follows:
aDD 2 = −
i
aDD 3 = −
xi xj
j
(11)
42 2 ε εjj εkk xi xj xk ii 3 kB T kB T kB T i
×
3 3 εii εjj ii jj n n ∗2 ∗2 J DD j 2,ij kB T kB T 3 ,i ,j i ij
3 ii3 jj3 kk
ij ik jk
j
k
n,i n,j n,k ∗2 ∗2 ∗2 J DD i j k 3,ijk
(12)
In the above equations, εii is segment energy parameter of component i, ii is segment size parameter of component i, ni is number of dipole moments in component i and ∗2 is dimeni sionless squared dipole moment. J2,ij and J3,ijk are approximated integrals over the reference-fluid pair-correlation function and over three-body correlation functions. The approximations of these integrals are obtained by power series functions in which the constants are adjusted to simulation data. However the quadrupolar (QQ) terms are given by the following equations:
aQQ = − 2
3 2 4
i
j
xi xj
5 5 εii εjj ii jj n n Q ∗2 J kB T kB T 7 Q,i Q,j j 2,ij ij
(13)
εAB /k [K]
KAB
|Q| [DA]
[D]
Reference
1.85 1.7
[18] [18] [18] [18] [18] [18] [20] [22] [22]
4.4 0.072 0.135
= aQQ 3
3
2406.700 2152.400
3 3 4
i
xi xj
j
× nQ,i nQ,j Qj∗3 J3,ij +
×
ε 3/2 ε 3/2 15/2 15/2 jj ii jj ii kB T
kB T
42 3
ij7
3 3 4
2
i
j
xi xj xk
k
5 5 5 εii εjj εkk ii jj kk n n n Q ∗2 Qj∗2 Qk∗2 J3,ijk kB T kB T kB T 3 3 3 Q,i Q,j Q,k i ij ik jk
(14)
where Qj∗2 denotes the dimensionless square quadrupole moment and nQ,j is the number of quadrupole moments in a molecule j. To use the PC-SAFT equation with dipolar and quadrupolar terms in the multiphase calculations, the component parameters are required. Table 2 lists the component parameters which are used in the computations. Parameters are taken from different sources [19,21,23]. Two association sites are assumed for water while four association sites are used for ethylene glycol. 4. Gupta et al. (1991) algorithm As discussed in Section 2, the proposed methodology of Gupta (1991) for simultaneous stability and multiphase flash calculation provided a unified set of simultaneous non-linear equations with inequality constraints. Gupta et al. [1] have proposed to partition the system of equations into two groups to be solved in two nested loops. The first group consisted of Equations (1), (2) and (5). These equations were solved for phase fractions (˛k ) and stability variables ( k ) by the Newton-Raphson method in the inner loop. The second group consisted of Equations (3) and (4) in the outer loop to calculate mole fractions. The mole fractions were updated in the outer loop and accelerated by the use of the General Dominant Eigenvalue Method (GDEM) [24]. The K-values, which were initiated with Michelsen stability analysis [12], were re-calculated from Equation (7) by the use of a chosen model in the outer loop. As previously indicated, modifications were included in their work to avoid singular and ill-conditioned Jacobian. The stability equation (Equation (5)) was replaced by the following equation: ˛k k =0 (˛k+ k )
(15)
The equation describes the equivalent behaviour of Equation (5) in the appearance and disappearance of a phase. It improves the conditioning of Jacobian significantly. The singularities of the Jacobian, when both phase fraction and its corresponding stability variable become zero, were also avoided in the way that both would never become equal to zero during computation. Furthermore, during computations, the phase fractions become negative which of course violates inequality 6. In that case Gupta et al. forced the phase fraction to be zero. The corresponding stability variable was then re-calculated from Equation (2). Moreover, during computations, some of the appearing phases may have the same compositions and compressibility (i.e. triv-
172
N.M. Alsaifi, P. Englezos / Fluid Phase Equilibria 302 (2011) 169–178
ial solution). If that occurs, Gupta et al. [1] suggested adding the phases together to decrease the number of phases and the next computations were performed only for the distinct phases.
1.0 Vapour Water-Rich Liquid Hydrocarbon-Rich Liquid
Phase Fraction
0.8
5. The proposed algorithm In this section, the proposed algorithm is given and its advantages over Gupta’s methodology are discussed. Similar to Gupta et al. approach, Equations (1), (2) and (5) are solved in the inner loop but without modifying Equation (5) and with the inclusion of the inequality constraint (Equation (6)). This requires an iterative method other than Newton-Raphson method which cannot handle inequality constraints. A recent iterative method which could solve systems of nonlinear equalities and inequalities has been proposed by Macconi et al. [25]. The method will be referred as trust-region Gauss–Newton method. It is implemented efficiently using MATLAB in a solver called TRESNEI [26]. Unlike the Newton-Raphson method, the trust-region Gauss–Newton method is capable of handling inequality constraints and has global convergence properties. The use of TRESNEI will avoid some of the tedious Gupta et al work of controlling phase fraction and re-calculated stability variable. TRESNEI would help to obtain phase fraction and stability variables from solving directly Equations (1), (2), (5) and (6). Similar to Gupta et al. [1], the mole fraction vector is selected as the iterative variable in the flash computations. However, our updating procedure differs from that of Gupta et al. [1]. The updating procedure of Kinoshita and Takamatsu [27] is utilized in our approach. We found that the Kinoshita and Takamatsu updating procedure is relatively insensitive to the initial estimates of mole fraction. The combination of trust-region Gauss–Newton method and Kinoshita and Takamatsu makes it possible to avoid Gupta’s initialization strategies [28], particularly for water-hydrocarbon containing systems. In this study, the initial estimates of mole fraction are based on specifying a value of 0.99 to any potential component that may form liquid and to the lightest component in the vapour phase. This initialization is found to be effective for water-hydrocarbon containing systems. For example, the water-pentane system might form a maximum of three phases (vapour–liquid–liquid). It is adequate in our approach to start with a mole fraction value of 0.99 for pentane in the pentane rich liquidphase, a value of 0.99 for water in the water-rich liquid phase and a value of 0.99 for the lightest component in the vapour phase. Unlike the Gupta et al. algorithm, the current algorithm does not encounter a problem when the phase fraction approaches zero. Furthermore, an important feature of the current algorithm is its ability to handle more than one absent phase (i.e. an unstable phase). The feature would assist to figure out which absent phases would appear first when the P–T conditions are changed. This feature will be illustrated in details in Section 6.1.1.
0.6
0.4
0.2
0.0 400
410
420
430
440
450
460
Temperature / K Fig. 1. Phase composition of water–paraffin mixture at 24.1 bar predicted by dipolar PC-SAFT.
It should be noted that SAFT is not cubic and the method of finding density (or volume) roots differs from that of cubic equations of state. To determine the density roots, it is necessary to use a reliable iterative method for solving the equation of state for density. It is also important to use a small termination tolerance on the function value (ε = 10−14 ). TRESNEI works effectively in finding density roots. Besides TRESNEI, we tested trust-region dogleg method, which is a variant of the Powell dogleg method [29]. This method is also successful in finding the density roots of all the cases that have been studied. There are other reported reliable methods that could be used to find the density roots of SAFT such as the extended dogleg method [30], homotopy continuation [31], and interval analysis [32]. 6. Evaluation of the proposed algorithm In this section, several examples are considered for evaluating the proposed scheme of several systems. In all studied systems, the binary interaction parameter is set to zero. Of course, the PC-SAFT may not be able to predict accurately the behaviour of many systems if binary interaction parameters are not employed. However, it is interesting to assess SAFT as a predictive tool. 6.1. Water/octane/hexane/pentane/butane/propane mixture The first considered mixture is a synthetic system of water and paraffins. The system has been studied before by Heidemann [4], Peng and Robinson [33] and Erbar [34]. The overall mole composition is made up of 0.1667 propane, 0.1667 n-butane, 0.2
Table 3 Phase fractions and stability variables of water–paraffin mixture at 24.1 bar. T (K)
460.000a 454.513a 450.000a 441.531a 420.000a 406.981a 390.000b a b c
Water-rich liquid phase (LW )
Hydrocarbon-rich liquid phase (LHC )
Vapour phase (V)
˛LW
LW
˛LHC
LHC
˛V
V
0.0000 0.0000 0.0000 0.0000 0.1948 0.2418 0.2517
0.4816 0.3629 0.2382 0.0000 0.0000 0.0000 –
0.0000 0.0000 0.0483 0.1356 0.5163 0.7582 0.7483
0.0566 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
1.0000 1.0000 0.9517 0.8644 0.2889 0.0000 0.0000
– – – – – – 0.1849
Vapour is the reference phase. Water-rich liquid phase is the reference phase. Intel Pentium Dual E2200 2.2 GHz.
Mixture state
Iter #
CPUc time (s)
Vapour only LHC incipient Vapour-LHC LW incipient V-LW -LHC LW -LHC LW -LHC
13 13 13 13 13 12 12
10 11 10 10 9 17 7
6.1.1. One phase region (vapour only) At a temperature equal to 460 K, the vapour phase fraction is equal to 1 and the other two phase fractions of liquid phases are zero. This indicates that the two liquid phases are absent and we only have a vapour phase. Although one could delete a phase from the three phases during computation when the system shows only one phase (here vapour phase), the computation is conducted without deleting any phase. The stability variables of the two liquid phases at this temperature are positive (0.4816 for water-rich liquid phase and 0.05660 for hydrocarbon-rich liquid phase). The stability variable value of the hydrocarbon-rich liquid phase is much less than that of the water-rich liquid phase. This indicates that decreasing the temperature makes the hydrocarbon-rich liquid phase appear first. The change of phases and stability variables vs. iteration number are depicted in Fig. 2. The figure illustrates the effectiveness of the proposed scheme with the trust-region Gauss–Newton method in providing fast convergence. Computation is started with zero values for stability variables and phase fractions. As iteration proceeds, the phase fraction of the vapour (˛V ) and hydrocarbon-rich liquid (˛LHC ) appear first while the water-rich liquid phase (˛LW ) remains zero and its stability variable (LW ) becomes positive. At the third iteration, the hydrocarbon-rich liquid phase disappears and its stability variable (LHC ) attains a positive value. Both liquid stability variables converge to positive values while their phase fractions become zero indicating that they are unstable phases. The smaller magnitude of the stability variable of hydrocarbon-rich liquid phase confirms the behaviour observed before in Fig. 1 in which hydrocarbon-rich liquid phase would appear first if temperature is decreased. 6.1.2. Incipient hydrocarbon-rich liquid and vapour–liquid phase region As the temperature is decreased, the hydrocarbon-rich liquid phase is about to appear at 454.513 K where both its phase fraction and stability variable are zero as shown in Table 3. This demonstrates that at this condition (24.1 bar and 454.513 K) the hydrocarbon-rich liquid phase is incipient and a slight decrease in temperature brings the system to the two phase region; namely
1.0 0.8 0.6 0.4 0.2 0.0
(a)
0.6
(b)
173
αV
0.5 0.4
αL
0.3 0.2
θL
0.1 0.0 0.6
W
W
(c)
0.5
αL
0.4 0.3
θL
0.2 0.1
HC
HC
0.0 0
2
4
6
8
10
12
14
Iteration Number Fig. 2. Phase fraction (˛) and stability variable () vs. iteration number for water–paraffin mixture at 24.1 bar and 460 K for: (a) vapour (stable); (b) water-rich liquid phase (unstable); (c) hydrocarbon-rich liquid phase (unstable).
vapour and hydrocarbon-rich liquid phases. An example of the vapour and hydrocarbon-rich liquid phase is at a temperature of 450 K. The phase fraction of the water-rich liquid phase is still zero at this temperature indicating that it is absent and its stability variable is still positive (0.2382). Therefore, the water-rich liquid phase is still unstable at this condition but its stability value gets smaller. More insight into the behaviour of the phase fraction and stability variables during the iteration is shown in Figs. 3 and 4 at 454.513 and 450 K; respectively. The iteration represents the number of times the composition vectors are updated. As shown in Fig. 3, the phase fraction and stability variable converge to zero for hydrocarbon-rich liquid phase at the 5th iteration. In Fig. 4, the phase fraction of hydrocarbon-rich liquid phase starts develop-
Phase Fraction or Stability Variable
n-pentane, 0.0667 hexane, 0.1333 n-octane and 0.2667 water. The PC-SAFT with the proposed methodology is applied to study this mixture at fixed pressure (P = 24.1 bar) and a temperature range of 390–460 K. Both association and dipolar terms are included to take account of hydrogen bonding and dipole–dipole interactions of water molecules, respectively. The condensation curves are obtained and shown in Fig. 1. The general characteristic of the curves are similar to those obtained before [4,33,34] although different equations of state were utilized. As seen in Fig. 1, when the mixture is cooled from its vapour state, the formation of a hydrocarbon-rich liquid phase appears first. To demonstrate the proposed scheme, it is advantageous to show the computational details in this example. It should be noted that the system may have up to three phases in the selected temperature range (vapour only, vapour–liquid, liquid–liquid and vapour–liquid–liquid). An isobaric (24.1 bar) temperature path is followed by starting from 460 K (vapour only) to 390 K (liquid–liquid). We assume that the system consists of a maximum of three phases and the vapour is the reference phase. The reference phase could be switched to one of the liquid phases when the system does not show vapour phase as demonstrated below. Table 3 shows values for the phase fractions and stability variables for this mixture for seven temperatures from 390 to 460 K at 24.1 bar. The initial values of phase fractions and stability variables are arbitrary and set to zeros in this example. In the following, computation details of Table 3 are demonstrated for each region.
Phase Fraction or Stability Variable
N.M. Alsaifi, P. Englezos / Fluid Phase Equilibria 302 (2011) 169–178
1.0 0.8 0.6 0.4 0.2 0.0
(a)
αV
0.5
(b)
0.4
αL
0.3 0.2
θL
0.1
W
W
0.0 0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00
(c)
αL θL
0
2
4
6
8
10
12
HC
HC
14
Iteration Number Fig. 3. Phase fraction (˛) and stability variable () vs. iteration number for water–paraffin mixture at 24.1 bar and 454.513 K for: (a) vapour (stable); (b) waterrich liquid phase (unstable); (c) hydrocarbon-rich liquid phase (incipient).
N.M. Alsaifi, P. Englezos / Fluid Phase Equilibria 302 (2011) 169–178
1.0 0.8 0.6 0.4 0.2 0.0
(a)
αV
0.35 0.30 0.25 0.20 0.15 0.10 0.05 0.00
Phase Fraction or Stability Variable
Phase Fraction or Stability Variable
174
(b)
αL θL
0.18 0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00
(c)
W
W
αL
HC
θL
1.0 0.8 0.6 0.4 0.2 0.0
(a)
0.25
(b)
0.20
αV
0.15
αL
0.10
θL
0.05
W
0.00 0.6
(c)
0.5
αL
0.4
HC
0.3
HC
W
θL
0.2
HC
0.1 0.0 0
2
4
6
8
Iteration Number
10
12
14
0
2
4
6
8
10
12
14
Iteration Number
Fig. 4. Phase fraction (˛) and stability variable () vs. iteration number for water–paraffin mixture at 24.1 bar and 450 K for: (a) vapour (stable); (b) water-rich liquid phase (unstable); (c) hydrocarbon-rich liquid phase (stable).
Fig. 6. Phase fraction (˛) and stability variable () vs. iteration number for water–paraffin mixture at 24.1 bar and 420 K for: (a) vapour (stable); (b) water-rich liquid phase (stable); (c) hydrocarbon-rich liquid phase (stable).
ing after the 2nd iteration where its stability variable approaches zero and then the phase becomes present at the third iteration. The appearance of the hydrocarbon liquid phase is accompanied by a decrease in the vapour phase fraction which converges to 0.9517.
in the vapour–liquid–liquid region. An illustration of the characteristics of the vapour–liquid–liquid region is given in Table 3 and Fig. 6 at 420 K. As seen in Table 3, the liquid phase fractions are positive and their corresponding stability variables are zero indicating that they are present. In Fig. 6, the change of the stability variable and phase fraction with iteration is depicted. Both stability variables of the two liquid phases are zero while their phase fractions are positive showing that they are stable.
1.0 0.8 0.6 0.4 0.2 0.0
(a)
0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00
(b)
0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00
(c)
αV
αL θL
6.1.4. Liquid–liquid phase If we continue cooling the mixture, the vapour phase, which is the reference phase, disappears at 406.981 K indicating that the mixture exhibits two liquid phases; namely, water-rich and
Phase Fraction or Stability Variable
Phase Fraction or Stability Variable
6.1.3. Incipient water-rich liquid and vapour–liquid–liquid phase region The water-rich liquid phase is incipient at temperature equals to 441.531 K as shown in Table 3 where both stability variable and phase fraction have zero values. Both phase fraction and stability variable are zero (i.e. incipient) after the 3rd iteration as shown in Fig. 5. If the temperature is slightly decreased, the mixture coexists
W
W
αL
HC
θL
HC
1.0 0.8 0.6 0.4 0.2 0.0
(a)
0.30
(b)
0.25
αV
0.20
αL
0.15
θL
0.10 0.05
W
W
0.00 1.0
(c)
0.8 0.6
αL
0.4
θL
HC
HC
0.2 0.0
0
2
4
6
8
10
12
14
Iteration Number Fig. 5. Phase fraction (˛) and stability variable () vs. iteration number for water–paraffin mixture at 24.1 bar and 441.53 K for: (a) vapour (stable); (b) waterrich liquid phase (incipient); (c) hydrocarbon-rich liquid phase (stable).
0
2
4
6
8
10
12
14
Iteration Number Fig. 7. Phase fraction (˛) and stability variable () vs. iteration number for water–paraffin mixture at 24.1 bar and 406.981 K for: (a) vapour (incipient); (b) water-rich liquid phase (stable); (c) hydrocarbon-rich liquid phase (stable).
N.M. Alsaifi, P. Englezos / Fluid Phase Equilibria 302 (2011) 169–178
0.4
(a)
0.3
Phase Fraction or Stability Variable
hydrocarbon-rich liquid phases. At this point, both the two liquid phase fractions are positive and their corresponding stability variables are zero whereas the vapour phase fraction is zero. A fast convergence of phase fraction and the stability variables is achieved as shown in Fig. 7. It should be emphasized here, as previously indicated, that the proposed scheme does not encounter any computation difficulty when the phase fraction approaches zero even if that phase is the reference phase. The disappearance of the vapour phase at temperature equal or less than 406.981 K and 24.1 bar makes it necessary to switch the reference phase to one of the liquid phases if computation is required for temperature less than 406.981 K. If the reference phase is maintained to be the vapour phase, a decrease in temperature would cause the vapour phase fraction to be negative. A negative value is of course not a physical value. However, it is an indication of improper selection of the reference phase. The reference phase is switched to the water-rich liquid phase for the flash calculation at a temperature less than 406.981 K where the system exhibits liquid–liquid equilibrium. For example, at 390 K, the vapour phase is absent and its stability variable is positive (0.1849) as shown in Table 3. The characteristics of the phase fractions and stability variables with iterations are shown in Fig. 8. It is clear that the phase fraction of vapour is zero and its corresponding stability variable (see Fig. 8) becomes positive at the 2nd iteration. Finally, it should be noted that Equation (8) is verified in all the studied conditions.
175
0.2
αL
0.1 0.0 0.20
W
(b)
0.15
αV
0.10
θV
0.05 0.00 0.8
(c)
0.6
αL
0.4
θL
0.2
HC
HC
0.0 0
2
4
6
8
10
12
14
Iteration Number Fig. 8. Phase fraction (˛) and stability variable () vs. iteration number for water–paraffin mixture at 24.1 bar and 390 K for: (a) vapour (unstable); (b) waterrich liquid phase (stable); (c) hydrocarbon-rich liquid phase (stable).
0.7
The second mixture that was examined is the ethylene glycol–water–propane–methane quaternary. Dipole–dipole interactions are considered for both ethylene glycol and water. The dipole–dipole interactions are represented with two polar functional groups of the ethylene glycol. Hydrogen bonding interactions are also taken into account for water and ethylene glycol. Experimental data was reported for this mixture by Ng et al. [35] at two different P–T values. The system shows vapour–liquid–liquid behaviour according to the reported experimental P–T values. A comparison of the experimental and the calculated mole fraction of this system with the proposed scheme is shown in Table 4. As seen in the table, dipolar PC-SAFT predicts the mole fractions of this system well. The proposed scheme is also applied for this mixture at a range of pressures at a fixed temperature (283.15 K) to show the change of phase fractions of the three phases. The phase fraction vs. pressure is depicted in Fig. 9. As seen in the figure, as the pressure increases the mole phase fraction of vapour decreases while the phase fraction of hydrocarbon-rich liquid phase increases. How-
0.6
Mole Phase Fraction
6.2. Ethylene glycol/water/propane/methane
Vapour Water-Rich Liquid Hydrocarbon-Rich Liquid
0.5 0.4 0.3 0.2 0.1 0.0 10
20
30
40
50
60
70
80
90
Pressure / bar Fig. 9. Mole phase fraction vs. pressure of methane–propane–water–ethylene glycol at 283.15 K predicted by dipolar PC-SAFT.
Table 4 Equilibrium phase composition (mole fraction) of methane–propane–water–ethylene glycol system. T = 273.15 K
MEG Water Propane Methane
V (exp.)
Feed
V
0.092 0.317 0.284 0.307
3.620E−07 1.207E−04 0.1566 0.8433
T = 283.15 K
MEG Water Propane Methane
Feed
V
0.090 0.311 0.288 0.311
8.746E−07 2.266E−04 0.1961 0.8037
LHC
P = 70.4 bar LW 0.2236 0.7702 2.562E−03 3.719E−03
– – 0.1542 0.8458
2.204E−06 1.277E−04 0.5394 0.4605
3.400E−06 4.000E−05 0.5660 0.4340
V (exp.)
LHC
P = 69.0 bar
– – 0.1984 0.8016
4.497E−06 2.211E−04 0.5834 0.4163
LW (exp.)
LHC (exp.)
0.2243 0.7736 2.770E−04 1.790E−03 LW (exp.)
LHC (exp.)
LW
3.000E−06 7.200E−05 0.6018 0.3982
0.2235 0.7700 2.928E−03 3.583E−03
0.2244 0.7736 3.760E−04 1.580E−03
176
N.M. Alsaifi, P. Englezos / Fluid Phase Equilibria 302 (2011) 169–178
ever, the phase fraction of the water-rich liquid phase is incipient at temperature remains constant as the pressure increases. There is no experimental data to verify this behaviour but prediction of similar systems (No. 10, 11 and 12 in Table 6) which have been studied and compared to experimental at different pressures supports this behaviour. The proposed scheme performs the calculation of this system without any difficulty although the mole phase fraction was approaching zero for hydrocarbon-rich liquid phase around 15 bar and for vapour phase around 84 bar. In this example, the computations were performed for the vapour–liquid–liquid region but can be extended to other phase regions as done with water-paraffins system. The vapour phase was selected as the reference phase in this example. Furthermore, the number of iterations did not exceed 18 for each selected P–T values while the total CPU time is about 60 s.
0.4
αV
Phase Fraction or Stability Variable
0.3
6.3. Ethylene glycol/water/hexane/propane/CO2 /methane
0.2 0.1 0.0 0.6 0.5
αL
0.4 0.3
θL
0.2 0.1
6.4. Water-hydrocarbon and other systems The proposed algorithm was successfully tested for several water-hydrocarbon and other systems listed in Table 6. The listed systems exhibit various phase behaviours including vapour–liquid, liquid–liquid and vapour–liquid–liquid equilibrium. For the listed systems, two SAFT versions [19,36] have been utilized with the addition of hydrogen bonding and dipolar interactions. The reader
W
0.0 0.5 0.4 0.3
αL
0.2
θL
0.1
The experimental data of this six-component mixture were reported by Ng et al. [35] at two different temperatures (263.15 and 283.15 K) and a fixed pressure of 68 bar. The reported data show that this system exhibits vapour–liquid–liquid equilibrium at these conditions. To examine this system with the proposed scheme, the PC-SAFT is utilized with the inclusion of dipolar and hydrogen bonding interactions of ethylene glycol and water as well as quadrupole interactions for carbon dioxide. Fig. 10 represents the phase fraction and its corresponding stability variables vs. iteration. The stability variables are arbitrary initiated with non-zero values. Nonetheless, a fast convergence is achieved and Equations (5) and (6) are verified. At the 2nd iteration, the stability variables go to zero where the phase fractions of both liquids appear and proceed to converge at 10 iterations. The total CPU time is around 43 s. Table 5 illustrates the comparison between experimental and calculated mole fractions. It is apparent that the dipolar PC-SAFT compare well the experimental data although no adjustable parameter is introduced.
W
HC
HC
0.0 0
2
4
6
8
10
12
Iteration Number Fig. 10. Phase fraction (˛) and stability variable () vs. iteration number for ethylene glycol–water–hexane–carbon dioxide–methane mixture at 68 bar and 283.15 K for: (a) vapour (stable); (b) water-rich liquid phase (stable); (c) hydrocarbon-rich liquid phase (stable).
could refer to the PhD Dossier of Alsaifi [37] for the details of calculations. In this work, the binary interaction parameter was set to zero to assess the predictive capability of PC-SAFT in multiphase systems. The prediction was very satisfactory of systems of water–ethylene glycol–hydrocarbons based on the study of systems in Sections 6.2 and 6.3 and systems 10, 11 and 12 in Table 6. Finally, although we have employed PC-SAFT with the addition of association, dipolar and quadrupole terms, the total CPU time was not huge even for cross-association systems (e.g. water–ethylene glycol systems). As far as our knowledge goes, the only available study for computing phase stability and equilibrium using SAFT for cross association systems is the one that conducted by Xu et al. [32]. They described a reliable method to compute phase stability and equilibrium using interval analysis. Their work was limited to mainly VLE and binary systems. The study showed that the computing time for conducting stability and equilibrium calculations for cross-association systems was large in comparison to ours. For example, the total CPU time of 1-butanol-ethanol sys-
Table 5 Equilibrium phase composition (mole fraction) for ethylene glycol–water–hexane–propane–carbon dioxide–methane. T = 263.15 K
EG H2 O C6 H14 C3 H8 CO2 CH4
V (exp.)
Feed
V
0.114 0.392 0.198 0.049 0.049 0.198
1.120E−07 6.232E−05 0.0034 0.0180 0.0832 0.8953
T = 283.15 K
EG H2 O C6 H14 C3 H8 CO2 CH4
Feed
V
0.115 0.397 0.195 0.049 0.049 0.195
6.899E−07 2.270E−04 0.0066 0.0275 0.1044 0.8613
– – 0.0041 0.0217 0.1199 0.8543
LHC
2.175E−06 9.172E−05 0.4183 0.1025 0.0992 0.3800
P = 68 bar
LW (exp.)
LHC (exp.)
LW
5.400E−06 4.400E−05 0.4500 0.1100 0.0948 0.3452
0.2241 0.7704 2.974E−04 3.982E−04 0.00114 0.00373
V (exp.)
LHC
P = 68 bar LHC (exp.)
LW
– – 0.0071 0.0297 0.1306 0.8326
8.737E−06 2.725E−04 0.4640 0.1120 0.0989 0.3248
1.200E−05 9.600E−05 0.4800 0.1154 0.0870 0.3176
0.2233 0.7708 4.019E−04 5.022E−04 1.28E−03 3.68E−03
0.2235 0.7705 2.00E−05 3.00E−05 0.00448 0.00148 LW (exp.)
0.2234 0.7704 2.600E−06 4.800E−05 4.45E−03 1.67E−03
N.M. Alsaifi, P. Englezos / Fluid Phase Equilibria 302 (2011) 169–178
177
Table 6 Other studied systems.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
System
Model
Ethane–propane–butane–pentane–hexane (VLE) Methane–ethane–propane–pentane–heptane–decane (VLE) Methane–water–methylcyclohexane (VLLE) Methane–water–heptane (VLLE) Methane–water–tert-butyl methyl ether (VLLE) Methane–water–neohexane (VLLE) Methane–water–decane (VLLE) p-xylene–water–ethanol (VLLE) Propane–water–pentane–octane (VLLE) Methane–water–ethylene glycol–methylcyclohexane (VLLE) Methane–water–ethylene glycol–hexane (VLLE) Methane–water–ethylene glycol–octane (VLLE) Methanol–hexane (LLE) Methanol–heptane (LLE) Methanol–decane (LLE) Ethanol–tetradecane (LLE)
PC-SAFT PC-SAFT Dipolar PC-SAFT Dipolar PC-SAFT Dipolar PC-SAFT Dipolar PC-SAFT Dipolar PC-SAFT Dipolar PC-SAFT Dipolar PC-SAFT Dipolar PC-SAFT Dipolar PC-SAFT Dipolar PC-SAFT Dipolar Simplified-SAFT [36] Dipolar Simplified-SAFT [36] Dipolar Simplified-SAFT [36] Dipolar Simplified-SAFT [36]
tem was 1816.1 s at T = 343 K, P = 0.35 bar and zbutanol = 0.65 (feed composition) using the original SAFT [38] with a binary interaction parameter equals to 0.011. For the sake of comparison, we studied this system at the same conditions and the same SAFT version [38]. The total CPU time is found to be about 30 s which is significantly faster than that of Xu et al. [32]. The phase fraction of vapour, which is the reference phase, is 0.97082 while the liquid phase fraction is 0.029177 and its corresponding stability variable is 1.5955e−22 indicating that it is stable. Our algorithm shows a third and absent phase with stability variable equals to 2.4917. 7. Conclusions The procedure suggested originally by Gupta (1990) to perform multiphase equilibrium and stability calculations was modified. The new scheme combines the trust-region Gauss–Newton method [25,26] and the updating procedure of Kinoshita and Takamatsu [27]. Unlike the classical iterative method, the trust-region Gauss–Newton method is capable of solving non-linear equations with inequality constrains and it has global convergence properties. The computational scheme utilizes a simple initialization procedure. The current algorithm does not encounter problems when the phase fraction approaches zero. Furthermore, it could handle more than one absent phase easily. The proposed computation scheme was utilized and illustrated the predictive capability of the PC-SASFT equation of state for water–ethylene glycol–hydrocarbon systems with or without carbon dioxide. Hydrogen bonding and polar (dipolar and quadrupole) interactions were included in the PC-SAFT and computations were performed without the use of binary interaction parameters. The scheme is proved to be efficient for multiphase calculations for water-hydrocarbon containing systems. List of symbols a reduced Helmholtz free energy f fugacity J2,ij integral over the reference-fluid pair-correlation function J3,ijk integral over three-body correlation function k phase Boltzmann constant (J/K) kB kij binary interaction parameter m segment number nc number of components ni number of dipole moments in component i Q* dimensionless quadrupole moment r reference phase T absolute temperature (K)
x y z
liquid mole fraction vapour mole fraction feed composition
Greek letters ˛ phase fraction stability variable ϕ fugacity coefficient εii segment energy parameter of component i dipole moment i dipole moment for component i ∗2 dimensionless squared dipole moment i number of phases number density ii segment size parameter of component i (Å) Acknowledgments We are grateful to the King Fahd University of Petroleum & Minerals for its support of Nayef Alsaifi in funding his graduate scholarship. Financial support from the Natural Sciences and Engineering Research Council of Canada (NSERC) is also acknowledged. The authors wish also to give their sincere thanks to Benedetta Morini and Margherita Porcelli for making TRESNEI code available for public use. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21]
A.K. Gupta, P.R. Bishnoi, N. Kalogerakis, Fluid Phase Equilibr. 63 (1991) 65–89. A.K. Gupta, P.R. Bishnoi, N. Kalogerakis, Gas Sep. Purif. 4 (1990) 215–222. G. Han, G.P. Rangaiah, Comput. Chem. Eng. 22 (1998) 897–911. R.A. Heidemann, AIChE J. 20 (1974) 847–855. R. Gautam, W.D. Seider, AIChE J. 25 (1979) 991–999. M.E. Soares, A.G. Medina, C. McDermott, N. Ashton, Chem. Engng. Sci. 37 (1982) 521–528. M. Castier, P. Rasmussen, A. Fredenslund, Chem. Engng. Sci. 44 (1989) 237–248. W.X. Song, M.A. Larson, Chem. Engng. Sci. 46 (1991) 2513–2523. H. Pan, A. Firoozabadi, SPE Res. Eval. Engng. Feb. (1998) 36–42. A. Lucia, L. Padmanabhan, S. Venkataraman, Comput. Chem. Eng. 24 (2000) 2557–2569. L.E. Baker, A.C. Pierce, K.D. Luks, SPE/DOE Second Joint Symp. Enhanced Oil Recovery, Tulsa, Oklahoma, 1981. M.L. Michelsen, Fluid Phase Equilib. 9 (1982) 1–20. G.A. Iglesias-Silva, A. Bonilla-Petriciolet, P.T. Eubank, J.C. Holste, K.R. Hall, Fluid Phase Equilib. 210 (2003) 229–245. M.L. Michelsen, Fluid Phase Equilib. 9 (1982) 21–35. L.X. Nghiem, Y. Li, Fluid Phase Equilib. 17 (1984) 77–95. J.-S. Wu, P.R. Bishnoi, Comput. Chem. Eng. 10 (1986) 269–276. P.R. Bishnoi, A.K. Gupta, P. Englezos, N.E. Kalogerakis, Fluid Phase Equilib. 53 (1989) 97–104. R.M. Abdel-Ghani, Master Thesis, University of Calgary, Canada, 1995. J. Gross, G. Sadowski, Ind. Eng. Chem. Res. 41 (2002) 5510–5515. J. Gross, J. Vrabec, AIChE J. 52 (2006) 1194–1204. J. Gross, AIChE J. 51 (2005) 2556–2568.
178 [22] [23] [24] [25] [26]
N.M. Alsaifi, P. Englezos / Fluid Phase Equilibria 302 (2011) 169–178
G.S. Rushbrook, G. Stell, J.S. Hoye, Mol. Phys. 26 (1973) 1199–1215. N.M. Al-Saifi, E.Z. Hamad, P. Englezos, Fluid Phase Equilib. 271 (2008) 82–93. C.M. Crowe, M. Nishio, AIChE J. 21 (1975) 528–533. M. Macconi, B. Morini, M. Porcelli, Appl. Numer. Math. 59 (2009) 859–876. B. Morini, M. Porcelli, Comput. Optim. Appl., (2010) DOI: 10.1007/s10589-010rr9327-5. [27] M. Kinoshita, T. Takamatsu, Comput. Chem. Eng. 10 (1986) 353–360. [28] A.K. Gupta, Ph.D. Thesis, Dept. of Chem. and Petr. Eng., University of Calgary, Calgary, Canada, 1990. [29] P. Rabinowitz, Numerical Methods for Nonlinear Algebraic Equations, Ch.7, 1970.
[30] N. Aslam, A.K. Sunol, Ind. Eng. Chem. Res. 45 (2006) 3303–3310. [31] A. Lucia, Q. Luo, Adv. Environ. Res. 6 (2001) 123–134. [32] G. Xu, J.F. Brennecke, M.A. Stadtherr, Ind. Eng. Chem. Res. 41 (2002) 938–952. [33] D.-Y. Peng, D.B. Robinson, Canad. J. Chem. Eng. 54 (1976) 595–599. [34] J.H. Erbar, Proc. 52nd Ann. National. Gas Processors Assoc. Meeting, 62, 1973. [35] H.-J. Ng, C.-J. Chen, D.B. Robinson, GPA Research Report RR-92, 1986. [36] F. Yuan-H, S.I. Sandler, Ind. Eng. Chem. Res. 34 (1995) 1897–1909. [37] N.M. Al-Saifi, PhD Dossier, Department of Chem. & Bio Eng., UBC, 2009. [38] S.H. Huang, M. Radosz, Ind. Eng. Chem. Res. 30 (1991) 1994–2005.