hydrocarbon systems at high pressures and elevated temperatures

Fluid Phase Equilibria 489 (2019) 48e66

Contents lists available at ScienceDirect

Fluid Phase Equilibria j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / fl u i d

A robust multiphase isenthalpic flash model for water/solvent/ hydrocarbon systems at high pressures and elevated temperatures Desheng Huang, Daoyong Yang* Petroleum Systems Engineering, Faculty of Engineering and Applied Science, University of Regina, Regina, S4S 0A2, Canada

a r t i c l e i n f o

a b s t r a c t

Article history: Received 11 October 2018 Received in revised form 4 February 2019 Accepted 5 February 2019 Available online 8 February 2019

Based on the assumption that only the presence of water and solvent is considered in the aqueous phase, a robust and pragmatic water-associated isenthalpic flash (WAIF) model has been developed to perform multiphase isenthalpic flash calculations for water/solvent/hydrocarbon mixtures at high pressures and elevated temperatures. More specifically, the phase boundary diagram in enthalpy-temperature (H-T) space is constructed first for a given water/solvent/hydrocarbon system at a specified pressure. Then, the phase transition enthalpies and temperatures are determined from the H-T diagram, and thus one-, two-, or three-phase isenthalpic flash calculations can be conducted directly without resorting to stability test for a given system at a specified enthalpy and pressure. In addition, a recursive method has been proposed to initialize the temperature so that the energy conservation equation and material balance equation can be solved simultaneously to overcome the convergence difficulties. The newly proposed isenthalpic flash model has been compared with the water-free isenthalpic flash (WFIF) model and conventional three-phase isenthalpic flash (CTIF) model through three case studies. The WAIF model is found to handle the narrow-boiling constraints effectively and accurately. Furthermore, as for Water/ C3H8/C16H34 mixture and Water/CO2/C16H34 mixture, a good agreement is achieved between the WFIF model and the WAIF model as the solubility of solvents (i.e., C3H8 and CO2) in the aqueous phase is quite small. With respect to Water/dimethyl ether (DME)/C16H34 mixture, however, the WAIF model is found to more accurately predict the phase behaviour of water/solvent/hydrocarbon mixtures than the WFIF model since the solubility of DME in the aqueous phase cannot be neglected. The case studies demonstrate that the newly developed water-associated isenthalpic flash model for water/solvent/hydrocarbon systems is robust, efficient, and accurate. © 2019 Elsevier B.V. All rights reserved.

Keywords: Isenthalpic flash Water/solvent/hydrocarbon systems Equation of state Water-associated model Mutual dissolution

1. Introduction The steam-based thermal recovery methods, such as steam assisted gravity drainage (SAGD) and cyclic steam stimulation (CSS), have been widely used to enhance oil recovery from heavy oil reservoirs [1,2]. Recently, co-injecting steam and solvent has been proved to be a promising technique in terms of oil production rate and ultimate oil recovery as it combines the advantages from both heat transfer of steam and mass transfer of solvent(s) to further reduce the viscosity of heavy oil [3e8]. Since water is inevitably closely associated with steam injection processes [9], three-phase vapour/hydrocarbon/aqueous equilibria of water/solvent/hydrocarbon mixtures are frequently encountered at the edge of a steam

* Corresponding author. E-mail address: [email protected] (D. Yang). https://doi.org/10.1016/j.fluid.2019.02.005 0378-3812/© 2019 Elsevier B.V. All rights reserved.

chamber [9,10]. It is still a challenging task to quantify such phase equilibria in the presence of water. In particular, isothermal flash shall be replaced with isenthalpic flash in steam-based EOR processes, though the prior knowledge of temperature lacks and narrow-boiling regions exist [11]. As such, it is of fundamental and practical importance to accurately quantify the phase behaviour near the steam chamber boundary for understanding the heavy oil recovery mechanisms and designing the steam-solvent co-injection processes [12]. In general, isothermal flash is commonly used to determine phase behaviour pertaining to conventional crudes as their isothermal conditions can be assumed [13,14]; however, thermalbased recovery processes cannot be treated to be isothermal due to the fact that high-temperature steam is injected into a lowtemperature heavy oil reservoir [15]. In this case, such a temperature variation voids the traditional assumption of constant temperature for conventional oil reservoirs during steam injection,

D. Huang, D. Yang / Fluid Phase Equilibria 489 (2019) 48e66

resulting in large errors associated with the isothermal flash [16]. As such, isenthalpic flash has been proposed to quantify the phase behaviour for the aforementioned steam injections where enthalpy can be assumed to be constant, though temperatures are varied within the matrix formation [17]. Physically, it is difficult to perform robust isenthalpic flash calculations due to the following two challenges: (1) number of the existing phases under flash conditions is unknown, and (2) the total enthalpy is significantly sensitive to temperature within the narrow boiling region [15e18]. The former increases the computational expenses, whereas the latter causes convergence problems for the energy conservation equation. Numerous efforts have been directed to tackle the aforementioned two challenges to certain extents. Michelsen [17] proposed a direct substitution multiphase isenthalpic flash method. The outer temperature iteration loop is used to satisfy the enthalpy constraint, while the inner loop is the traditional isothermal flash to determine the number of phases through stability analysis. By using a quasi-Newton iteration step to update K-values for the material balance equation and energy balance equation, Agarwal et al. [19] subsequently modified this method and proposed three schemes to perform isenthalpic flash calculations for multiphases and multicomponents. Both of these two methods found their convergence difficulties in the narrow boiling region. By simultaneously solving the matrix consisting of the Rachford-Rice equation, energy conservation equation, and stability equation, Gupta et al. [20] found that robustness of the isenthalpic flash calculation method has been improved for some cases; however, its application is still limited due to its inherent complexity. To overcome such a convergence problem, Zhu and Okuno [21] employed the condition number of the Jacobian matrix for the Rachford-Rice equation and the energy conservation equation to detect the presence of the narrow-boiling behaviour, while the Newton's method and bisection method are alternatively switched during the isenthalpic flash calculations, depending on the condition number of the Jacobian matrix. Recently, efforts have been made to improve the computational efficiency of isenthalpic flash calculations. By using the water-free isothermal flash proposed by Lapene et al. [22], Heidari et al. [18] modified the hybrid scheme proposed by Agarwal et al. [19] to drastically improve the computational speed by conducting pseudo-two-phase flash instead of conventional three-phase flash. The success of their method, nevertheless, is highly dependent on selection of appropriate initial K-values for water-containing mixtures. By using the new Rachford-Rice algorithm based on waterfree assumption [23], Li and Li [24] have improved the method proposed by Zhu and Okuno [21] to solve the inner loop isothermal flash calculation and outer loop energy conservation equation. However, all the aforementioned methods are not applicable to cases where a solvent is highly soluble in the aqueous phase at high pressures and elevated temperatures. In this paper, a robust and pragmatic technique has been developed to perform isenthalpic flash calculation for water/solvent/hydrocarbon systems at high pressures and elevated temperatures. The augmented water-free isothermal flash algorithm [25,26] is incorporated in the inner loop to solve the Rachford-Rice equation for phase fractions and compositions, while the outer loop is used to update the temperature by satisfying the enthalpy constraint. In addition, the phase boundary enthalpies are determined prior to isenthalpic flash for a given feed composition. Then, the isenthalpic flash calculation can be conducted straightly at one, two-, or three-phase region without performing stability test for a given feed at a specified enthalpy. Also, comparisons have been made to handle the narrow-boiling behaviour with respect to the WFIF model and the CTIF model under various conditions.

49

2. Mathematical formulations 2.1. Peng-Robinson equation of state (PR EOS) PR EOS [27] with the van der Waals’ mixing rule is a thermodynamic model used in this work because of its simplicity and computational accuracy [28e37]. The PR EOS is formulated as follows,

P¼

RT a  V  b VðV þ bÞ þ bðV  bÞ

(1)

a ¼ ac aðTr ; uÞ

ac ¼

b¼

(2a)

0:457235R2 Tc 2 Pc

(2b)

0:0777969RTc Pc

(2c)

where P is pressure, T is temperature, V is molar volume, a and b are constants at a given temperature, R is the universal gas constant, and aðTr ; uÞ is the alpha function which is a function of the reduced temperature Tr, and acentric factor u. The following alpha function for predicting the vapour pressure of pure substances is given by Peng and Robinson [27],

h



aðTr ; uÞ ¼ 1 þ 0:37464 þ 1:54226u  0:26992u2



1  T 0:5 r

i2 (3)

For a mixture system, a and b can be calculated with the van der Waals’ mixing rule, i.e.,

a¼

Nc X Nc X

 0:5  yi yj ai aj ð1  dij

(4a)

i¼1 j¼1

b¼

Nc X

yi bi

(4b)

i¼1

where the subscript i and j are the component index, Nc is the number of the components, yi and yj are the mole fractions of the ith and jth component, respectively, and dij is the binary interaction parameter (BIP) between the ith and jth component. The parameters a and b expressed in a dimensionless form are listed as,

A¼

aP R2 T 2

(5a)

B¼

bP RT

(5b)

The compressibility factor Z is calculated from the PR EOS,

    Z 3 þ ðB  1ÞZ 2 þ A  3B2  2B Z þ B3 þ B2  AB ¼ 0

(6)

The derivatives of various parameters associated with the PR EOS with respect to T are expressed as,

50

D. Huang, D. Yang / Fluid Phase Equilibria 489 (2019) 48e66

  vai 0:37464 þ 1:54226u  0:26992u2 ¼ 0:457235R2 T 1:5 ci vT   i. h  1 þ ð0:37464 þ 1:54226u0:26992u2 1  T 0:5 ri    Pci T 0:5 (7)   Nc X Nc  0:5  vaj vamk X va  ¼ þ aj i 1  dij ai 0:5xik xjk ai aj vT vT vT i¼1 j¼1

(8a)

vAmk P va 2A ¼ 2 2 mk  mk vT T R T vT

(8b)

vBmk B ¼  mk vT T

(8c)

where the subscript k is the phase index, and m is the mixture index. It is worthwhile to mention that Eq. (7) is obtained by deriving Eq. (2a) with respect to T, and alpha function is expressed as Eq. (3). 2.2. Three-phase isothermal flash For three-phase equilibrium of water/solvent/hydrocarbon mixtures, equal chemical potential of each component in each phase are to be satisfied,

fix ¼ fiy ¼ fiw ; i ¼ 1; 2; …; Nc

(9)

where fix, fiy, and fiw represent the fugacity of component i in hydrocarbon-rich liquid phase, vapour phase, and aqueous phase, respectively. The component and phase material balance constraints can be expressed as,

bx þ by þ bw ¼ 1

(10)

bx xi þ by yi þ bw wi ¼ zi

(11)

where bx , by , and bw are phase fractions in the hydrocarbon-rich liquid phase, vapour phase, and aqueous phase, respectively, xi, yi, and wi are mole fractions of component i in the hydrocarbon-rich liquid phase, vapour phase, and aqueous phase, respectively. Additionally, the mole fraction of each equilibrium phase must sum to unity, Nc X i¼1

xi ¼

Nc X

yi ¼

i¼1

Nc X

wi ¼ 1

(12)

i¼1

Introducing the equilibrium ratios (hydrocarbon-rich liquid phase is chosen as the reference phase),

8 yi fix > > > < Kiy ¼ xi ¼ fiy > w f > > : Kiw ¼ i ¼ ix xi fiw

(13)

where Kiy and Kiw are equilibrium ratios of component i in the vapour phase and aqueous phase with respect to the hydrocarbonrich liquid phase, respectively, fix , fiy , and fiw are fugacity coefficients of the ith component in the hydrocarbon-rich liquid phase, vapour phase, and aqueous phase, respectively.

Combining Eqs. (10)e(13) yields the following two RachfordRice equations for conventional full three-phase isothermal flash, which are commonly used to solve for phase fractions and compositions [38,39],

8   Nc Nc > X X zi Kiy  1 > > >   ¼ ðy  x Þ ¼ ¼0 RR y > i i > 1 þ by Kiy  1 þ bw ðKiw  1Þ > > i¼1 i¼1 > > > <

> > > Nc Nc > X X z ðK  1Þ > >  i iw > ðwi  xi Þ ¼ ¼0 RRw ¼ > > 1 þ b  1 þ bw ðKiw  1Þ K > iy y > i¼1 i¼1 : (14) By solving Eqs. (13) and (14), phase compositions can be determined from the material balance equations, 8

zi > >   xi ¼ > > 1 þ by Kiy  1 þ bw ðKiw  1Þ < yi ¼ Kiy xi > > > > : wi ¼ Kiw xi

(15)

Based on the water-free assumption that the aqueous phase consists of only pure water, Li and Li [23] simplified the work of Lapene et al. [22] and Okuno et al. [40] and then proposed an optimization method for the three-phase isothermal flash with water-free assumption. The objective function and constraints are expressed as follows [23],

8 Nc  i h   X > > * * < min : f by ¼ þ Kwz zi ln Kiy  Kwy  o n  i¼2 > > : subject to : K *  K b  min K *  z ; K *  K z ; is1 iy i iy i y wy wz wz (16) where

8 1  yw > * > > < K wy ¼ 1  x w > 1  z > w * > : K wz ¼ 1  xw

(17)

Eq. (17) is the “effective” equilibrium constants defined by Tang and Saha [41]. By combining Eqs. (16) and (17), phase compositions can then be obtained,

8 zi >  > ; isH2 O xi ¼  > > * > Kiy  K wy þ K *wz > > > > > < 1 xi ¼ ; i ¼ H2 O > K > ww > > > > y ¼ Kiy xi > > > i > : wi ¼ 1

(18)

where

8 Kwy > > Kww ¼ > > > yw > > > < T Pcw Kwy ¼ > Tcw P > > > > > > Pw > : yw ¼ sat P

(19)

where Pcw and Tcw are the critical pressure and temperature of water, P w sat denotes the saturation pressure of pure water that can

D. Huang, D. Yang / Fluid Phase Equilibria 489 (2019) 48e66

be obtained by using the empirical correlation developed by Bridgeman and Aldrich [42]. Considering that the solubility of some solvents in the aqueous phase cannot be negligible under reservoir conditions, Pang and Li [25,26] improved the water-free three-phase isothermal flash and proposed a new isothermal flash algorithm on the basis of the augmented water-free assumption. Based on the assumption that only a solvent (e.g., CO2) or dimethyl ether (DME) and water are present in the aqueous phase at high pressures and elevated temperatures [43], the objective function and constraints are modified as follows [25,26],

51

component i in phase k, xiNp is the mole fraction of component i in the reference phase, Ht is the total molar enthalpy of the system, Hk is the molar enthalpy of phase k, and Hspec is the specified molar enthalpy. The total molar enthalpy of Np phases is expressed as [17,45],

Ht ¼

Np X



dep bk H IGM þ Hk k



(24)

k¼1

More specifically, the molar ideal gas mixture enthalpy of phase k as well as the molar ideal gas enthalpy of component i can be

8 Nc i h h i   X     > > > zi ln 1  by 1  Kiy  bw  zw ln 1  by 1  Kwy  bw ð1  Kww Þ min : f by ¼ > > > > isH2 O; DME or CO2 > h h i >    i   > > < zCO2 ln 1  by 1  KCO2 y  bw 1  KCO2 w  zDME ln 1  by 1  KDMEy  bw ð1  KDMEw Þ 8    > by 1  Kiy þ bw  min 1  zi ; 1  Kiy zi ; isH2 O; DME or CO2 > > >  > > < > > > > subject to : by 1  Kwy  bw ð1  Kww Þ  min 1 zw ; 1  Kwy zw ; 1  Kww zw > > > b 1  KDMEy  bwð1  KDMEw Þ  min 1  zDME ; 1  KDMEy zDME ; 1  KDMEw z DME > > > > y : : by 1  KCO2 y  bw 1  KCO2 w  min 1  zCO2 ; 1  KCO2 y zCO2 ; 1  KCO2 w zCO2

(20)

calculated by using the following formula [46,47], By solving the above formula, the phase compositions can be expressed as,

8 z >  i  > ; isH2 O; DME or CO2 xi ¼ > > 1 þ b K > iy  1  bw y > > > > < zi   xi ¼ ; i ¼ H2 O; DME or CO2 1 þ b  1 þ bw ðKiw  1Þ K iy y > > > > > > yi ¼ Kiy xi > > > : wi ¼ Kiw xi ; i ¼ H2 O; DME or CO2

HIGM ¼ k

Nc X

xik HIG i

(25)

i¼1

 .  . 0 0 2 2 2 þ C 0P3i T 3  T 30 3 HIG i ¼ C P1i ðT  T0 Þ þ C P2i T  T 0  . 4 þ C 0P4i T 4  T 40 (26) (21)

where Ht is the total molar enthalpy, H IGM is the molar ideal gas k mixture enthalpy of phase k, H IG i is the molar ideal gas enthalpy of dep

component i, H k 2.3. Isenthalpic flash Isenthalpic flash calculations correspond to obtaining the temperature T, phase mole fraction bk , and phase composition xik which maximize the total entropy St of a system at a specified pressure P, enthalpy Hspec, and feed composition zi [17,44].

St ¼

Np X

bk Sk ; k ¼ 1; 2; …; Np

(22)

is the molar enthalpy departure, C 0P1i , C 0P2i , C 0P3i ,

and C 0P4i are ideal gas heat capacity coefficients of component i, and T0 is reference temperature, which is 273.15 K in this work. Then, the molar enthalpy departure for phase k is given by, dep

Hk

   

vA 21:5 Bmk ¼ RT 2 mk þ RTAmk vT " #   Z þ 1 þ 20:5 Bmk  ln k  þ RTðZk  1Þ Zk þ 1  20:5 Bmk

(27)

k¼1

The following Np e 1 material balance equation and energy conservation equation are to be satisfied,

gk ¼

NX p 1

  xik  xiNp ¼ 0

(23a)

i¼1

gNp ¼ Ht  Hspec ¼

Np X

bk Hk  Hspec ¼ 0

(23b)

k¼1

where St is the total molar entropy of the system, bk is the phase mole fraction of phase k, Sk is the molar entropy of phase k, zi is the overall mole fraction of component i, xik is the mole fraction of

3. Modified algorithms This section presents a modified isenthalpic flash algorithm based on the augmented water-free assumption, which is referred as water-associated isenthalpic flash (WAIF) model in this work. Two methods have been made available to perform the isenthalpic flash calculations in the literature, i.e., 1) the water-free isenthalpic flash (WFIF) model which is based on the water-free assumption that the solvent solubility in the aqueous phase is zero, and 2) the conventional three-phase isenthalpic flash (CTIF) model which considers all the solvent solubility in the aqueous phase. The former enhances the computational efficiency by solving the modified

52

D. Huang, D. Yang / Fluid Phase Equilibria 489 (2019) 48e66

Rachford-Rice monotonic objective function [24], but it is not suitable for the system with a high solvent solubility in the aqueous phase. The latter can predict all the solvent solubility in the aqueous phase; however, it is time-consuming because it needs numerous stability analyses to determine the number of phases of the system involved. The main objective of developing a WAIF model in this

study is to improve the computational efficiency for water/solvent/ hydrocarbon systems by avoiding stability analysis as those for the CTIF model while maintaining prediction accuracy for solvent solubility in the aqueous phase compared with the WFIF model. The new isenthalpic flash algorithm adopts a modified Rachford-Rice monotonic function i.e., Eq. (20) and the negative

Fig. 1. Flowchart of the WAIF algorithm.

D. Huang, D. Yang / Fluid Phase Equilibria 489 (2019) 48e66

53

Table 1 Physical properties of component H2O/C3H8/C16H34 [49]. Component

zi (mol%)

TC (K)

PC (bar)

u

C 0P1i J/(mol$K)

C 0P2i J/(mol$K2)

C 0P3i J/(mol$K3)

C 0P4i J/(mol$K4)

H2O

75

647.3

220.89

0.344

32.20

0.001907

1:06  105

 3:596  109

C3H8

15

369.8

42.46

0.152

4.22

0.306300

 1:60  104

3:215  108

C16H34

10

717.0

14.19

0.742

13.00

1.529000

 8:50  104

1:850  107

Table 2 BIPs between the components in the H2O/C3H8/C16H34 mixture [24]. Component

H2O

C3H8

C16H34

H2O C3H8 C16H34

0 0.6841 0.3583

0.6841 0 0

0.3583 0 0

recursive method is proposed to constrain the selection of initial temperature and overcome the problems associated with phase appearance and disappearance that are frequently encountered in a thermal process. The procedures for the new multiphase isenthalpic flash for water/solvent/hydrocarbon systems are briefly described below.

Table 3 Physical properties of component H2O/CO2/C16H34 [49]. Component

zi (mol%)

TC (K)

PC (bar)

u

C 0P1i J/(mol$K)

C 0P2i J/(mol$K2)

C 0P3i J/(mol$K3) 105

C 0P4i J/(mol$K4)  3:596  109

H2O

75

647.3

220.89

0.344

32.200

0.001907

1:06 

CO2

15

304.2

73.76

0.225

19.795

0.073430

 5:602  105

1:715  108

C16H34

10

717.0

14.19

0.742

13.000

1.529000

 8:50  104

1:850  107

The flowchart of the new algorithm is shown in Fig. 1.

Table 4 BIPs between the components in the H2O/CO2/C16H34 mixture [22]. Component

H2O

CO2

C16H34

H2O CO2 C16H34

0 0.1896 0.3583

0.1896 0 0.1250

0.3583 0.1250 0

flash concept [48] for calculating phase fractions and compositions. Therefore, there is no need to perform stability analysis to determine the number of phases for a given mixture. As for the WAIF model, it is the first task to determine the upper phase transition temperatures (TU) for a given feed composition consisting of water/ solvent/hydrocarbon mixtures at a specified pressure. More specifically, TU1 is the temperature at which the vapour phase just appears in the system, TU2 denotes the temperature at which the aqueous phase disappears, and TU3 represents the temperature at which the hydrocarbon-rich liquid phase leaves the system. Accordingly, phase fractions and compositions at the phase boundary can also be determined. Then, the enthalpies at upper phase transition temperatures can be calculated, which is referred to as phase transition enthalpies. Finally, the relation (i.e., H-T diagram) between phase transition enthalpies and temperatures can be constructed. Therefore, one-, two-, or three-phase isenthalpic flash calculations can be conducted directly without performing stability analysis by comparing a specified enthalpy of a mixture to the phase transition enthalpy. The energy conservation equation or material balance equation may not converge if a poor initial temperature T1 is assumed for a given three-phase isenthalpic flash. Li and Li [24] proposed a scheme to switch from three-phase to two-phase isenthalpic flash without the prior knowledge of the existing phases. In this paper, a

1) Specify P, Hspec, zi, together with thermodynamic properties such as critical temperature Tc, critical pressure, Pc, acentric factor, u, and BIPs. 2) Construct the H-T diagram at the specified P. In this step, detailed procedures are given to determine phase transition temperatures (i.e., appearance of a vapour phase at TU1, disappearance of liquid water at TU2, and disappearance of hydrocarbon liquid at TU3) and the corresponding phase transition enthalpies. (a) Appearance of a vapour phase at TU1. Perform three-phase isothermal flash by estimating a series of temperatures at the prespecified pressure until the vapour phase fraction approaches zero. Then estimated temperature is TU1, and the phase transition enthalpy at TU1 is calculated as,

HT U1 ¼ HHC bHC þ HA bA

(28a)

(b) Disappearance of liquid water at TU2. Perform three-phase isothermal flash by estimating a series of temperatures at the prespecified pressure until the aqueous phase fraction approaches zero. The estimated temperature is TU2, and the phase transition enthalpy at TU2 is calculated as,

HT U2 ¼ HHC bHC þ HV bV

(28b)

(c) Disappearance of liquid hydrocarbon at TU3. Perform twophase isothermal flash by estimating a series of temperatures at the specified pressure until the hydrogen liquid

Table 5 Physical properties of component H2O/DME/C16H34 [49]. Component

zi (mol%)

TC (K)

PC (bar)

u

C 0P1i J/(mol$K)

C 0P2i J/(mol$K2)

C 0P3i J/(mol$K3) 105

C 0P4i J/(mol$K4)  3:596  109

H2O

75

647.3

220.89

0.344

32.20

0.001907

1:060 

DME

15

400.05

52.92

0.200

17.02

0.179100

 5:234  105

 1:918  109

C16H34

10

717.0

14.19

0.742

13.00

1.529000

 8:500  104

1:850  107

54

D. Huang, D. Yang / Fluid Phase Equilibria 489 (2019) 48e66

Table 6 BIPs between the components in the H2O/DME/C16H34 mixture [50]. Component

H2O

DME

C16H34

H2O DME C16H34

0 0.1700 0.3583

0.1700 0 0.0150

0.3583 0.0150 0

phase fraction approaches zero. The estimated temperature is TU3, and the phase transition enthalpy at TU3 is calculated as,

HT U3 ¼ HV bV

(29c)

3) Determine the multiphase isenthalpic flash region from the H-T diagram by comparison with Hspec and obtain the lower temperature TL as well as the upper temperature TU according to the determined isenthalpic flash region. 4) Set the initial temperature Titer ¼ (TL þ TU)/2, the iteration number iter ¼ 1, and tolerance ε ¼ 106 .

Fig. 2. Comparison between phase transition enthalpies and temperatures obtained from the WAIF model and the WFIF model [24] at 80 bar for a mixture of 75 mol% water, 15 mol % C3H8, and 10 mol% C16H34.

Fig. 3. Comparison between the temperatures and phase fractions obtained from the WAIF model and the WFIF model in the three-phase region for the water/C3H8/C16H34 mixture at 80 bar and different enthalpies.

D. Huang, D. Yang / Fluid Phase Equilibria 489 (2019) 48e66

5) Perform isothermal flash at P and Titer at the specific phase region obtained from Step #3 and calculate phase fractions, bk , and compositions, xik. The method for initializing equilibrium ratios for the WAIF model is the same as that proposed by Pang and Li [25], while the modified Rachford-Rice equation (i.e., Eq. (20)) is solved by Newton with the line search method [25]. 6) If unphysical results, i.e., bk ;ð0; 1Þ occur, let T iterþ1 ¼ T iter þ DT if T iter > ðT L þ T U Þ=2 or T iterþ1 ¼ T iter  DT if T iter < ðT L þ T U Þ=2 and go to Step #5.

55

7) Calculate the residual of the energy conservation equation, g iter Np ,

giter Np ¼ Ht  Hspec

(30)



8) If g iter Np < ε, terminate the calculation; otherwise, continue to Step #9.

Fig. 4. Comparison of the temperatures and phase compositions that are obtained by using the WAIF model and the WFIF in three-phase region for the water/C3H8/C16H34 mixture at 80 bar and different enthalpies for (a) vapour phase; and (b) hydrocarbon-rich liquid phase.

56

D. Huang, D. Yang / Fluid Phase Equilibria 489 (2019) 48e66

Fig. 5. Mole fractions of C3H8 in the aqueous phase calculated by the WAIF model at 80 bar and different enthalpies in the three-phase region.

9) Update Titer. The secant method is used to update the temperature,

T iterþ1 ¼ T iter 

 iter  giter  T iter1 Np T iter1 g iter Np  g Np

(31a)

10) Check to see if T L < T iterþ1 < T U . If so, let iter ¼ iter þ1 and go to Step #11; otherwise, update Titer by using the Regula Falsi method [19],

T iterþ1 ¼ T iter 

 L  U g iter Np T  T g LNp  g U Np

(31b)

where g LNp and g U Np are the residual of the energy conservation equation at T L and T U , respectively.



11) If g iter1

< ε, terminate the calculation and output T iter1 , xik, Np and bk ; otherwise, go back to Step #5.

Fig. 6. The number of iterations required by the WAIF model and the WFIF model in the three-phase region for the water/C3H8/C16H34 mixture at 80 bar.

D. Huang, D. Yang / Fluid Phase Equilibria 489 (2019) 48e66

57

Fig. 7. Comparison between phase transition enthalpies and temperatures obtained from the WAIF model and the WFIF model at 80 bar for a mixture of 75 mol% water, 15 mol% CO2, and 10 mol% C16H34.

4. Results and discussion Three representative water/solvent/hydrocarbon mixtures with different solvent solubilities in the aqueous phase are selected as case studies to not only test the robustness, efficiency, and accuracy of the WAIF model, but also handle the narrow-boiling constraints, especially in the three-phase region. These three mixtures associated with the same water/solvent/hydrocarbon ratio are selected to

simulate a typical scenario which is frequently encountered in a solvent-steam co-injection process. More specifically, water represents the condensed steam, solvent represents the co-injected volatile components, and C16H34 represents the oil phase. The results calculated from the WAIF model are compared against with those obtained from the WFIF model and the CTIF model. The WFIF model developed by Li and Li [24] has been demonstrated as accurate, efficient, and robust for water/hydrocarbon mixtures with

Fig. 8. Comparison between the temperatures and phase fractions obtained from the WAIF model and the WFIF model in the three-phase region for the water/CO2/C16H34 mixture at 80 bar and different enthalpies.

58

D. Huang, D. Yang / Fluid Phase Equilibria 489 (2019) 48e66

water-insoluble solvents. Physical properties and BIPs for the mixtures are listed in Tables 1e6 [22,24,49,50]. Hereinafter, V represents the vapour phase, HC represents the hydrocarbon-rich liquid phase, and A represents the aqueous phase. 4.1. Case #1: water/C3H8/C16H34 mixture The water/C3H8/C16H34 mixture is tested at 80 bar and different

specified enthalpies. Li and Li [24] also used this mixture to compare the results obtained by the WFIF model with those of Zhu and Okuno [11]. Fig. 2 shows the comparison between phase transition enthalpies and temperatures, i.e., phase boundaries in HT space, obtained by the WAIF model and the WFIF model at 80 bar. As can be seen, the phase transition enthalpies and temperatures obtained by the two isenthalpic flash models are in an excellent agreement. The narrow-boiling behaviour occurs in the three-

Fig. 9. Comparison between the temperatures and phase compositions in the vapour phase that are obtained from the WAIF model and the WFIF model in the three-phase region for the water/CO2/C16H34 mixture at 80 bar and different enthalpies: (a) vapour phase compositions and (b) absolute deviations of vapour phase compositions.

D. Huang, D. Yang / Fluid Phase Equilibria 489 (2019) 48e66

phase region, which is indicated by a significant change in enthalpy with a small change in temperature. The phase transition temperatures at 80 bar are determined to be TU1 ¼ 491.04 K, TU2 ¼ 549.47 K, and TU3 ¼ 598.85 K. A number of specified enthalpies are selected from the three-phase region to evaluate the performance of the new isenthalpic flash model to handle the narrow-boiling constraints, as shown in Figs. 3e6. Therefore, phase transition enthalpies for the three-phase region are determined to be 14073 and 17599 J/mol from Fig. 2. Accordingly, TL ¼ 491.04 K, TU ¼ 549.47 K, and T1 ¼ (TL þ TU)/2 ¼ 520.255 K, respectively. Fig. 3 compares the temperatures and phase fractions of each phase calculated by these two isenthalpic flash models at 80 bar in the narrow-boiling region, i.e., the three-phase region. A wide range of specified enthalpies is selected from this region to test the robustness and accuracy of the WAIF model. The calculated temperatures and phase fractions by the two isenthalpic flash models agree well. An absolute average relative deviation (AARD) and a maximum absolute relative deviation (MARD) calculated between the two models for predicting the temperature are 0.001% and 0.003% as well as an AARD of 0.037% and a MARD of 0.057% for predicting the phase fractions, which indicates that the WAIF model is capable of dealing with narrow-boiling regions with high accuracy and robustness. Fig. 4 illustrates the comparison of the temperatures and mole fractions of individual components that are calculated by using the WFIF model and those calculated by using the WAIF model at 80 bar and different enthalpies in the threephase region for the vapour phase and hydrocarbon-rich liquid phase, respectively. A good agreement is observed between the calculated temperatures and mole fractions of individual components by these two methods. It is interesting to note that mole fractions of water in vapour phase and hydrocarbon-rich liquid phase both increase with an increment in temperature. Dissolved water in C16H34 acts as a low viscosity solvent that results in viscosity reduction of hydrocarbon-rich liquid phase, which is preferred in thermal recovery processes [51]. Fig. 5 shows the mole fractions of C3H8 in the aqueous phase calculated by the WAIF model. As can be seen, the largest mole fraction of C3H8 in the aqueous phase is less than 2:0  105 , which explains why the phase transition enthalpies, temperatures, phase

59

fractions and compositions obtained by the two isenthalpic flash models are in an excellent agreement. Based on the results from Case #1, the two isenthalpic flash models can also be applied to other water/hydrocarbon mixtures in which the solvent has a similar solubility of C3H8 in the aqueous phase. Fig. 6 shows the comparison of the number of iterations required by the two models in the narrow-boiling three-phase region. Compared with the WFIF model, the number of iterations required the WAIF model in the three-phase region, especially near the phase boundary between the three-phase and two-phase vapour-liquid regions, is reduced, indicating that the new isenthalpic flash model is efficient and robust in dealing with the narrow-boiling constraints. It is worthwhile noting that the initial temperature has a significant effect on the isenthalpic flash calculations. As such, a recursive method has been proposed in this study to accelerate the convergence speed within the narrow-boiling regions, resulting in a reduced number of iterations and thus less computational time. 4.2. Case #2: water/CO2/C16H34 mixture This mixture, which includes a slightly water-soluble solvent, CO2, is used to evaluate the performance of the two isenthalpic flash models. Multiphase isenthalpic flash calculations are conducted at 80 bar and different specified enthalpies. Fig. 7 shows the comparison between phase transition enthalpies and temperatures obtained by the WAIF model and those obtained by the WFIF model. The three-phase region temperature are located at 316.99e547.82 K, while the corresponding enthalpy is located at -39091-14903 J/mol. Compared to Case #1, there exists a relatively small discrepancy between the phase boundaries in the H-T space obtained from the two isenthalpic flash models. Fig. 8 illustrates the comparison between the temperatures and phase fractions of the vapour phase, hydrocarbon-rich liquid phase, and aqueous phase obtained from the two isenthalpic flash models at 80 bar in the three-phase region. As can be seen from Fig. 8, there is a good agreement between the temperatures and phase fractions obtained by the two isenthalpic flash models. The AARDs for predicting the temperature and phase fractions are 0.037% and 0.734% and the MARDs for predicting the temperature and phase fractions

Fig. 10. Mole fractions of CO2 in the aqueous phase calculated by the WAIF model at 80 bar and different enthalpies within the three-phase region.

60

D. Huang, D. Yang / Fluid Phase Equilibria 489 (2019) 48e66

Fig. 11. The number of iterations required by the WAIF model and the WFIF model in the three-phase region for the water/CO2/C16H34 mixture at 80 bar.

are 0.057% and 2.673%, respectively, implying that the two isenthalpic flash models both work well for the water/CO2/C16H34 mixture with a relatively small deviation. Fig. 9 compares the temperatures and phase compositions in the vapour phase that are obtained by the WFIF model and those obtained from the WAIF model for the water/CO2/C16H34 mixture at 80 bar within the threephase region. Obviously, an excellent agreement between the temperatures and phase compositions in the vapour phase is achieved.

Fig. 10 shows the dissolution of CO2 in the aqueous phase obtained by the WAIF model at 80 bar in the three-phase region. Compared to Fig. 5, the solubility of CO2 in the aqueous phase is a slightly higher than that of C3H8 in the aqueous phase, which explains why there exists a relatively small discrepancy between the phase boundaries in the H-T space obtained from the two isenthalpic flash models. It is obvious that the mole fraction of CO2 in the aqueous phase increases first, and then decreases with an increase in temperature. In addition, the largest mole fraction of CO2

Fig. 12. Comparison between phase transition enthalpies and temperatures obtained from the WAIF model and the WFIF model at 80 bar for a mixture of 75 mol% water, 15 mol% DME, and 10 mol% C16H34.

D. Huang, D. Yang / Fluid Phase Equilibria 489 (2019) 48e66

61

Fig. 13. Comparison between the temperatures and phase fractions obtained from the WAIF model and the WFIF model for the water/DME/C16H34 mixture at 80 bar and different enthalpies in the three-phase region.

in the aqueous phase is less than 5:0  103 . Fig. 11 depicts the number of iterations and temperatures obtained by the WAIF model and the WFIF model in the three-phase region for water/ CO2/C16H34 mixture at 80 bar. Again, the number of iterations of the WAIF model in which a recursive method is used for initializing the temperature for the inner loop isothermal flash is reduced compared with the WFIF model, illustrating the computational efficiency of the new isenthalpic flash model.

4.3. Case #3: water/DME/C16H34 mixture This mixture, which is comprised of water, heavy hydrocarbon C16H34 and DME, a partially water-soluble solvent due to its slight polarity [52], is used to evaluate the performance of the two isenthalpic flash models. The main objective of Case #3 is to examine the superiority of the WAIF model of dealing with narrow-boiling regions, i.e., the three-phase region, in terms of accuracy,

Fig. 14. Absolute deviations of phase fractions obtained by applying the WAIF model and the WFIF model in the three-phase region for the water/DME/C16H34 mixture at 80 bar.

62

D. Huang, D. Yang / Fluid Phase Equilibria 489 (2019) 48e66

robustness, and efficiency. Another purpose is to investigate the deviation of the results obtained from the WAIF model against with those from the WFIF model which uses water-free assumption and the CTIF model which considers the mutual solubility of components in all phases. Fig. 12 compares the phase transition enthalpies and temperatures obtained from the two isenthalpic flash models. The threephase region temperature is located at 513.95e549.03 K, while

the corresponding enthalpy is located at -9000-15500 J/mol. Compared to Cases #1 and #2, large discrepancies are observed between the phase boundaries in H-T space obtained from the WAIF model and the WFIF model due to that the solubility of DME in the aqueous phase is much higher than those of C3H8 and CO2 in the aqueous phase. More specifically, DME (i.e., CH3eOeCH3) is fairly soluble in water due to the presence of its oxygen atom which can interact with hydrogen atom in water to form hydrogen bonds [53].

Fig. 15. Comparison between the temperatures and phase compositions in the vapour phase that are obtained from the WAIF model and the WFIF in three-phase region for the water/DME/C16H34 mixture at 80 bar: (a) vapour phase compositions and (b) absolute deviations of vapour phase compositions.

D. Huang, D. Yang / Fluid Phase Equilibria 489 (2019) 48e66

As for C3H8 (i.e., CH3eCH2eCH3), the eCH3 group is hydrophobic and cannot form hydrogen bonds with water, thus its solubility in the aqueous phase is quite small [54,55]. With respect to CO2 (i.e., O]C]O), it has a slight negative charge near the oxygen atom and a slight positive charge near the carbon atom, resulting in polarized bonds [53]. CO2 is slightly soluble in water because water molecules are attracted to these polar areas and the bond between the carbon

63

and oxygen atoms is not as polar as the bond between the hydrogen and oxygen atoms [53]. Fig. 13 presents the comparison between the temperatures and phase fractions that are obtained from the WAIF model and those obtained by the WFIF model for the water/DME/C16H34 mixture at 80 bar and different enthalpies in the three-phase region. It is obvious that there exist large discrepancies between the

Fig. 16. Comparison between the temperatures and phase compositions in the vapour phase that are obtained from the WAIF model and the CTIF model in three-phase region for the water/DME/C16H34 mixture at 80 bar: (a) vapour phase compositions and (b) absolute deviations of vapour phase compositions.

64

D. Huang, D. Yang / Fluid Phase Equilibria 489 (2019) 48e66

temperatures and phases fractions that are obtained from the two isenthalpic flash models. This is due to the fact that the water-free isenthalpic flash model does not consider the DME solubility in the aqueous phase, i.e., the DME solubility in the aqueous phase is assumed as zero. In reality, however, the solubility of DME in the aqueous phase cannot be neglected. Fig. 14 plots the absolute deviations of phase fractions obtained by applying the WFIF model compared with the WAIF model in the three-phase region for the water/DME/C16H34 mixture at 80 bar. As can be seen, large absolute deviations are observed for the prediction of phase fractions by applying the WFIF model, especially for vapour phase fractions. The AARD and MARD for predicting phase fractions are 17.022% and 236.756%, respectively, implying that the WFIF model can cause significant deviations if applied to Case #3. Fig. 15 demonstrates the temperatures and phase compositions obtained from the WAIF model and the WFIF model for the water/ DME/C16H34 mixture at 80 bar in the three-phase region over a wide range of specified enthalpies for the vapour phase. As can be seen, large deviations are found between the calculated phase compositions by applying the two isenthalpic flash models. More specifically, the WFIF model yields an AARD of 2.128% and MARD of 4.597% for predicting phase compositions, indicating that the water-free assumption is not applied for Case #3. Fig. 16 compares the temperatures and phase compositions obtained from the WAIF model and the CTIF model in which the solubility of all components including the heavy hydrocarbon in the aqueous phase is taken into account. It should be noted that numerous stability analysis is required in the CTIF model to determine the number of phases for a given system at different enthalpies, which leads to more computational time compared with the WAIF model, though the CTIF model is able to predict the solubility of all the solvent in the aqueous phase. As can be seen from Fig. 16, an excellent agreement is achieved between the calculated phase compositions by these two isenthalpic flash methods. More specifically, the WAIF model is calculated to yield an AARD of 0.003% and MARD of 0.005% for predicting phase compositions, implying that the WAIF model is also capable of predicting the phase behaviour of Case #3 without

compromising prediction accuracy while improving computational efficiency. Fig. 17 plots the variation of DME solubility in the aqueous phase calculated by the WAIF model for the water/DME/C16H34 mixture at 80 bar in the three-phase region. As can be observed, the smallest mole fraction of DME in the aqueous phase is larger than 2:0  102 , which cannot be neglected. This explains why there exist large discrepancies between the phase transition enthalpies, temperatures, phase fractions, and compositions determined by the two isenthalpic flash models. In addition, the solubility of DME in the aqueous phase decreases with an increase in temperature, showing a different trend of solvent solubility in the aqueous phase compared to Cases #1 and #2. Fig. 18 shows the number of iterations used by the WAIF model and the CTIF model in the threephase narrow-boiling region for the water/DME/C16H34 mixture at 80 bar. Clearly, the number of iterations is significantly reduced by using the WAIF model, indicating the superiority of the new isenthalpic flash model proposed in this study to handle the narrow-boiling constraints. The reason why the number of iterations of the CTIF model is higher than that of the WAIF model is that the phase appearance and disappearance issue could happen for the CTIF model. It should be noted that the number of iterations is increased slightly when the specified enthalpies are near the phase boundary between the three-phase and two-phase vapour-hydrocarbon regions, whereas it is still 10 or less required by the WAIF model. Fig. 19 compares the computational time obtained by the WAIF model and the CTIF model for water/DME/C16H34 mixture at 80 bar in the three-phase region. The recursive method for the initial temperature estimation is applied in the WAIF model, which accelerates the energy conservation equation convergence speed. As can be seen from Fig. 19, the computational time is significantly reduced by using the WAIF model because a modified RachfordRice equation (i.e., Eq. (20)) is solved by the WAIF model instead of using two Rachford-Rice equations (i.e., Eq. (14)) used by the CTIF model, implying that the WAIF model is superior than the CTIF model in terms of the computational efficiency.

Fig. 17. Mole fractions of DME in the aqueous phase calculated by the WAIF model for the water/DME/C16H34 mixture at 80 bar in the three-phase region.

D. Huang, D. Yang / Fluid Phase Equilibria 489 (2019) 48e66

65

Fig. 18. The number of iterations required by the WAIF model and the CTIF model in the three-phase region for the water/DME/C16H34 mixture at 80 bar.

Fig. 19. Comparison between the computational time required by the WAIF model and the CTIF model in the three-phase region for the water/DME/C16H34 mixture at 80 bar.

5. Conclusions In this paper, practical computational techniques and formulations have been developed to determine phase boundaries in the enthalpy-temperature space (i.e., H-T diagram) together with water-associated isenthalpic flash calculations. The new isenthalpic flash model in this work can handle multiphase equilibria flash calculations at high pressures and elevated temperatures. By constructing an H-T phase boundary diagram, this new isenthalpic

flash model can be conducted directly without performing stability test for a given feed at a specified enthalpy and pressure. Compared to the traditional isothermal flash, the new isenthalpic flash model can be used to successfully deal with the convergence issues normally encountered in the narrow-boiling regions. In addition, both the WAIF model and the WFIF model work well for Cases #1 and #2 where the amount of solvent dissolved in the aqueous phase is quite small; however, large discrepancies are found for the WFIF model in Case #3 where the solubility of DME in the aqueous phase

66

D. Huang, D. Yang / Fluid Phase Equilibria 489 (2019) 48e66

cannot be neglected. Therefore, the newly developed WAIF model can be applied to calculate phase behaviour of water/solvent/hydrocarbon mixtures which are frequently encountered in solventsteam co-injection processes, especially for cases where the solvents are soluble in both the aqueous phase and the oil phase. Acknowledgements The authors acknowledge a Discovery Development Grant, a Discovery Grant and a Collaborative Research and Development (CRD) Grant from the Natural Sciences and Engineering Research Council (NSERC) of Canada to D. Yang and Combustion Solutions Inc., EHR Enhanced Hydrocarbon Recovery Inc., and Thermal Recovery Technologies Inc. for financial support. References [1] R.M. Butler, C.T. Yee, Progress in the in situ recovery of heavy oils and bitumen, J. Can. Petrol. Technol. 41 (1) (2002) 31e40. [2] A.M. Al Bahlani, T. Babadagli, SAGD laboratory experimental and numerical simulation studies: a review of current status and future issues, J. Petrol. Sci. Eng. 68 (3e4) (2009) 135e150. [3] K. Naderi, T. Babadagli, Solvent selection criteria and optimal application conditions for heavy-oil/bitumen recovery at elevated temperatures: a review and comparative analysis, J. Energy Resour. Technol. 138 (1) (2016) 0129041e012904-9. [4] D.W. Zhao, J. Wang, I.D. Gates, Optimized solvent-aided steam-flooding strategy for recovery of thin heavy oil reservoirs, Fuel 112 (2013) 50e59. [5] R.K. Jha, M. Kumar, I. Benson, E. Hanzlik, New insights into steam/solventcoinjection-process mechanism, SPE J. 18 (5) (2013) 867e877. [6] D. Ji, M. Dong, Z. Chen, Analysis of steam-solvent-bitumen phase behavior and solvent mass transfer for improving the performance of the ES-SAGD pocess, J. Petrol. Sci. Eng. 133 (2015) 826e837. [7] R.P. Leaute, B.S. Carey, Liquid addition to steam for enhancing recovery (LASER) of bitumen with CSS: results from the first pilot cycle, J. Can. Petrol. Technol. 46 (9) (2007) 22e30. [8] T. Kar, C. Ovalles, E. Rogel, J. Vien, B. Hascakir, The residual oil saturation determination for steam Assisted gravity drainage (SAGD) and solvent-SAGD, Fuel 172 (2016) 187e195. [9] W.R. Shu, K.J. Hartman, Effect of solvent on steam recovery of heavy oil, SPE Reservoir Eng. 3 (2) (1988) 457e465. [10] K.M. Brantferger, G.A. Pope, K. Sepehrnoori, Development of a thermodynamically consistent, fully implicit, equation-of-state, compositional steamflood simulator, in: Paper SPE 21253-MS, Presented at the SPE Symposium on Reservoir Simulation, Anaheim, CA, February 17-20, 1991. [11] D. Zhu, R. Okuno, Robust isenthalpic flash for multiphase water/hydrocarbon mixtures, SPE J. 20 (6) (2015) 1350e1365. [12] N. Jia, A. Memon, J. Gao, J.Y. Zuo, H. Zhao, H.J. Ng, H. Huang, Three-phase equilibrium study for heavy-oil/solvent/steam system at high temperatures, J. Can. Petrol. Technol. 50 (6) (2011) 68e79. [13] M.L. Michelsen, The isothermal flash problem. Part I. Stability, Fluid Phase Equilib. 9 (1) (1982) 1e19. [14] M.L. Michelsen, The isothermal flash problem. Part II. Phase-split calculation, Fluid Phase Equilib. 9 (1) (1982) 21e40. [15] D. Paterson, W. Yan, M.L. Michelsen, E.H. Stenby, Robust and efficient isenthalpic flash algorithms for thermal recovery of heavy oil, in: Paper SPE 179652 -MS, Presented at the SPE Improved Oil Recovery Conference, Tulsa, OK, April 11-13, 2016. [16] D. Zhu, R. Okuno, Multiphase isenthalpic flash integrated with stability analysis, Fluid Phase Equilib. 423 (2016) 203e219. [17] M.L. Michelsen, Multiphase isenthalpic and isentropic flash algorithms, Fluid Phase Equilib. 33 (1e2) (1987) 13e27. [18] M. Heidari, L.X. Nghiem, B.B. Maini, Improved isenthalpic multiphase flash calculations for thermal compositional simulators, in: Paper SPE 170029-MS, Presented at the SPE Heavy Oil Conference, Calgary, AB, June 10-12, 2014. [19] R.K. Agarwal, Y.K. Li, L.X. Nghiem, D.A. Coombe, Multiphase multicomponent isenthalpic flash calculations, J. Can. Petrol. Technol. 30 (3) (1991) 69e75. [20] A.K. Gupta, P.R. Bishnoi, N. Kalogerakis, Simultaneous multiphase isothermal/ isenthalpic flash and stability calculations for reacting/non-reacting systems, Gas Sep. Purif. 4 (4) (1990) 215e222. [21] D. Zhu, R. Okuno, A robust algorithm for isenthalpic flash of narrow-boiling fluids, Fluid Phase Equilib. 379 (2014) 26e51. [22] A. Lapene, D.V. Nichita, G. Debenest, M. Quintard, Three-phase free-water flash calculations using a new modified Rachford-Rice equation, Fluid Phase Equilib. 297 (1) (2010) 121e128. [23] R. Li, H. Li, New two-phase and three-phase Rachford-Rice algorithms based

on free-water assumption, Can. J. Chem. Eng. 96 (1) (2018) 390e403. [24] R. Li, H. Li, A robust three-phase isenthalpic flash algorithm based on freewater assumption, J. Energy Resour. Technol. 140 (3) (2018) 03290201e032902-10. [25] W. Pang, H. Li, An augmented free-water three-phase Rachford-Rice algorithm for CO2/hydrocarbons/water mixtures, Fluid Phase Equilib. 450 (2017) 86e98. [26] W. Pang, H. Li, Application of augmented free-water Rachford-Rice algorithm to water/hydrocarbons mixtures considering the dissolution of methane in the aqueous phase, Fluid Phase Equilib. 460 (2018) 75e84. [27] D.Y. Peng, D.B. Robinson, A new two-constant equation of state, Ind. Eng. Chem. Fundam. 15 (1) (1976) 59e64. [28] H. Li, S. Zheng, D. Yang, Enhanced swelling effect and viscosity reduction of solvent(s)/CO2/heavy-oil systems, SPE J. 18 (4) (2013) 695e707. [29] H. Li, D. Yang, Modified a function for Peng-Robinson equation of state to improve vapour pressure prediction of non-hydrocarbon and hydrocarbon compounds, Energy and Fuels 25 (1) (2011) 215e223. [30] H. Li, D. Yang, Phase behaviour of C3H8/n-C4H10/heavy-oil systems at high pressures and elevated temperatures, J. Can. Petrol. Technol. 52 (1) (2013) 30e40. [31] X. Li, H. Li, D. Yang, Determination of multiphase boundaries and swelling factors of solvent(s)-CO2-heavy oil systems at high pressures and elevated temperatures, Energy and Fuels 27 (3) (2013) 1293e1306. [32] Y. Shi, S. Zheng, D. Yang, Determination of individual diffusion coefficients of solvents-CO2-heavy oil systems with consideration of natural convection induced by swelling effect, Int. J. Heat Mass Transf. 107 (2017) 572e585. [33] Z. Chen, X. Li, D. Yang, Quantification of viscosity for solvents-heavy oil/ bitumen systems in the presence of water at high pressures and elevated temperatures, Ind. Eng. Chem. Res. 58 (2) (2019) 1044e1054. [34] Z. Chen, D. Yang, Correlations/estimation of equilibrium interfacial tension for methane/CO2-water/brine systems based on mutual solubility, Fluid Phase Equilib. 483 (2019) 197e208. [35] X. Dong, Y. Shi, D. Yang, Quantification of mutual mass transfer of CO2/N2-light oil systems under reservoir conditions, Ind. Chem. Eng. Res. 57 (48) (2018) 16495e16507. [36] Z. Chen, D. Yang, Quantification of phase behaviour of solvents-heavy oil systems in the presence of water at high pressures and elevated temperatures, Fuel 232 (2018) 803e816. [37] Z. Chen, D. Yang, Optimization of the reduced temperature associated with Peng-Robinson equation of state and soave-Redlich-Kwong equation of state to improve vapor pressure prediction for heavy hydrocarbon compounds, J. Chem. Eng. Data 62 (10) (2017) 3488e3500. [38] H.H. Rachford, J.D. Rice, Procedure for use of electronic digital computers in calculating flash vaporization hydrocarbon equilibrium, J. Petrol. Technol. 4 (10) (1952) 19e23. [39] Z. Li, A. Firoozabadi, General strategy for stability testing and phase-split calculation in two and three phases, SPE J. 17 (4) (2012) 1096e1107. [40] R. Okuno, R. Johns, K. Sepehrnoori, A new algorithm for Rachford-Rice for multiphase compositional simulation, SPE J. 15 (2) (2010) 313e325. [41] Y. Tang, S. Saha, An efficient method to calculate three-phase free-water flash for water-hydrocarbon systems, Ind. Eng. Chem. Res. 42 (1) (2003) 189e197. [42] O.C. Bridgeman, E.W. Aldrich, Vapor pressure Tables for water, J. Heat Tran. 86 (2) (1964) 279e286. [43] K. Sheng, R. Okuno, M. Wang, Dimethyl ether as an additive to steam for improved steam-assisted gravity drainage, SPE J. 23 (4) (2018) 1201e1222. [44] M.L. Michelsen, State function based flash specifications, Fluid Phase Equilib. 158 (1999) 617e626. [45] M.G. Kesler, B.I. Lee, Improve prediction of enthalpy of fractions, Hydrocarb. Process. 55 (1976) 153e158. [46] B.E. Poling, J.M. Prausnitz, J.P. O'connell, The Properties of Gases and Liquids, fifth ed., McGraw-Hill, New York, 2001. [47] C.A. Passut, R.P. Danner, Correlation of ideal gas enthalpy, heat capacity and entropy, Ind. Eng. Chem. Process Des. Dev. 11 (4) (1972) 543e546. [48] C.H. Whitson, M.L. Michelsen, The negative flash, Fluid Phase Equilib. 53 (1989) 51e71. [49] American Petroleum Institute (API), API Technical Data Book: Petroleum Refining, fourth ed., American Petroleum Institute, New York City, 1983. [50] K. Sheng, Analysis of Phase Behavior for Steam-Solvent Coinjection for Bitumen Recovery, M.A.Sc. Thesis, University of Alberta, Edmonton, AB, 2016. [51] C.A. Glandt, W.G. Chapman, The effect of water dissolution on oil viscosity, SPE Reservoir Eng. 10 (1) (1995) 59e64. [52] R.R. Ratnakar, B. Dindoruk, L. Wilson, Experimental investigation of DMEwater-crude oil phase behavior and PVT modeling for the application of DME-enhanced waterflooding, Fuel 182 (2016) 188e197. [53] J. Hine, P.K. Mookerjee, Structural effects on rates and equilibriums. XIX. Intrinsic hydrophilic character of organic compounds. Correlations in terms of structural contributions, J. Org. Chem. 40 (3) (1975) 292e298. [54] A. Chapoy, S. Mokraoui, A. Valtz, D. Richon, A.H. Mohammadi, B. Tohidi, Solubility measurement and modeling for the system propane-water from 277.62 to 368.16 K, Fluid Phase Equilib. 226 (2004) 213e220. [55] P. Flowers, K.H. Theopold, R. Langley, E.J. Neth, W.R. Robinson, Chemistry, Atoms First. OpenStax, Houston, 2016.