Determination of Multiphase Boundaries and ... - ACS Publications

Article pubs.acs.org/EF

Determination of Multiphase Boundaries and Swelling Factors of Solvent(s)−CO2−Heavy Oil Systems at High Pressures and Elevated Temperatures Xiaoli Li, Huazhou Li, and Daoyong Yang* Petroleum Systems Engineering, Faculty of Engineering and Applied Science, University of Regina, Regina, Canada, S4S 0A2 ABSTRACT: A generalized methodology has been proposed and successfully applied to determine multiphase boundaries as well as swelling factors of solvent(s)−CO2−heavy oil systems at high pressures and elevated temperatures. Experimentally, twoand three-phase boundaries and swelling factors have been respectively measured by conducting PVT tests in the temperature range of 280.45 to 396.15 K. Theoretically, the Peng−Robinson equation of state (PR EOS) combined with the modified alpha function has been applied to describe phase behavior of the solvent(s)−CO2−heavy oil systems. More specifically, an exponential distribution function is used to split a heavy oil sample, whereas the logarithm-type lumping method is employed to group single carbon numbers (SCNs) into multiple carbon numbers (MCNs). The exponents associated with two binary interaction parameter (BIP) correlations are respectively tuned for the alkane solvent-pseudocomponent pair and CO2-pseudocomponent pair to match the measured saturation pressures. It is found that six pseudocomponents combined with the BIP correlation as a function of critical volume is sufficient to predict saturation pressure with an absolute average relative deviation (AARD) of 5.07%. In addition, the PR EOS model associated with the selected parameters is applied to predict three-phase boundaries for a C3H8−CO2−heavy oil mixture and a n-C4H10−CO2−heavy oil mixture yielding an overall AARD of 4.58%. As for swelling factors, the Peneloux et al. method provides the minimum AARD of 2.09% in comparison with 4.00% from the Jhaveri et al. method and 2.41% from the Twu et al. method, respectively.

1. INTRODUCTION Thermal technique has been considered as the most conventional method to recover heavy oil resources.1 It has been found, however, to be inefficient for thin heavy oil formations mainly due to the excessive heat losses and limited drainage height.2,3 Solvent-based recovery methods, for example, vapor extraction (VAPEX),4 cyclic solvent injection (CSI),5 steamalternating-solvent injection,6 and huff-n-puff (cyclic solvent injection),7 have been becoming the most attractive alternatives to enhance the heavy oil production in thin payzones. In general, solvents can be either nonhydrocarbon (e.g., CO2) or hydrocarbon (e.g., CH4, C2H6, C3H8, and n-C4H10) compounds.7−10 It has been experimentally demonstrated that dissolution of solvent(s) into heavy oil can significantly improve the heavy oil mobility,11 enhance the swelling factor,8,9 and decrease the equilibrium interfacial tension.12 Such beneficial effects are closely associated with the phase behavior of solvent(s)−CO2−heavy oil systems under in situ reservoir conditions. Therefore, it is of fundamental and practical importance to accurately model the phase behavior of solvent(s)−CO2−heavy oil systems at high pressures and elevated temperatures.13,14 Erroneous results in phase behavior evaluation would be produced if the plus fraction (e.g., C7+) is treated as a single pseudocomponent.15−18 Splitting the C7+ fraction into SCNs and lumping them needs to be implemented to achieve more accurate phase behavior prediction.19 Lohrenz-Bray-Clark (LBC) viscosity correlation was one of the earliest attempts to apply an exponential-type distribution for splitting C7+ fraction.20 Pedersen et al. proposed that there exists an approximate linear relationship between carbon number and © 2013 American Chemical Society

the logarithm of its corresponding mole fraction if the carbon number is larger than six.21,22 As for heavy oil with low API gravity, such a linear relationship starts at C11 or even larger carbon numbers.23−25 A more complex molar distribution model used to split the plus fraction is the three-parametergamma distribution model proposed by Whitson.26−28 Furthermore, a weight-based lumping method suggested that each lumped pseudocomponent contains the same weight amount.22 Danesh et al. recommended that summation of the mole fractions times the logarithm of the molecular weight should be the same.29 As for the number of lumped pseudocomponents, more than 4 (e.g., 5, 8, 9, 16, 18, 40) pseudocomponents are assumed to be appropriate for characterizing the heavy oil.30−32 In spite of the rapid growth in computing capability of modern computers, compositional reservoir simulation employing many pseudocomponents is still a time-consuming task. Because the number of (pseudo)components describing an oil sample can be reduced significantly without impairing accuracy of the predicted phase behavior results, it is desirable to characterize heavy oil with a small number of pseudocomponents.29 Limited experimental phase behavior data of solvent(s)− CO2−heavy oil systems at high pressures and elevated temperatures have been made available in the literature, mainly due to the difficulties in performing phase behavior experiments of such systems.33 Although several binary interaction parameter (BIP) correlations for describing the phase behavior Received: November 15, 2012 Revised: January 16, 2013 Published: February 15, 2013 1293

dx.doi.org/10.1021/ef301866e | Energy Fuels 2013, 27, 1293−1306

Energy & Fuels

Article

of pure components have been proposed,34−37 few attempts have been made to study the effect of BIP correlation on phase behavior prediction of solvent(s)−CO2−heavy oil systems in a wide range of pressures and temperatures. In this article, a generalized methodology has been proposed and successfully applied to determine multiple phase boundaries (i.e., liquid−vapor (L1V) and liquid−liquid−vapor (L1L2V) phase boundaries) and swelling factors of solvent(s)− CO2−heavy oil systems at high pressures and elevated temperatures. More specifically, the experimentally measured L1V phase boundaries, L1L2V phase boundaries and swelling factors of solvent(s)−CO2−heavy oil systems with 15 feeds in the temperature range of 280.45 to 396.15 K, are used to evaluate the prediction accuracy of the PR EOS model with the modified alpha function by treating a heavy oil sample as multipseudocomponents. The exponents in the two BIP correlations for the alkane solvent-pseudocomponent pair and CO2-pseudocomponent pair are respectively tuned to fit the experimental saturation pressure data. As a result, six pseudocomponents together with BIP expressed as a function of critical volume is selected to be applied in PR EOS model to predict L1L2V three-phase boundaries and the swelling factors of solvent(s)−CO2−heavy oil systems. In addition, three volume shift strategies, that is, Peneloux et al. method,38 Jhaveri et al. method,39 and Twu et al. method40 are examined to compare their efficacy in predicting the swelling factors.

Table 1. Compositional Analysis of the Lloydminster Heavy Oil

2. EXPERIMENTAL SECTION 2.1. Materials. The heavy oil sample is collected from the Lloydminster area in Saskatchewan, Canada, with the molecular weight and specific gravity of 482 g/mol and 0.9997 g/cm3, respectively. The compositional analysis result of the heavy oil is presented in Table 1. The solvents of CO2, C3H8 and nC4H10 have purities of 99.998%, 99.99% and 99.5% (Praxair, Canada), respectively. 2.2. Experimental Setup. A mercury-free DBR PVT system (PVT-0150−100−200−316−155, DBR, Canada) is employed to perform all of the experiments of solvent(s)− CO2−heavy oil systems. A schematic diagram of the setup can be found in Figure 1. A pressure-volume-temperature (PVT) cell, which is the core component of the PVT system equipped in the air-bath has an inner diameter of 3.177 cm and a total length of 20.320 cm, whereas it can sustain pressure up to 69 000.0 kPa over the temperature range of 238.15 to 473.15 K. The air-bath temperature is controlled by using a microprocessor-based controller in conjunction with a resistance temperature device (RTD) sensor. A floating piston in the PVT cell isolates the test fluid from the hydraulic oil, whereas a magnetic mixer is equipped at the bottom of the PVT cell, which is highly favorable to sufficiently stir the solvent(s) and heavy oil mixture by setting its rotation speed. The mixture can be pressurized by moving the isolation piston downward in the PVT cell through a high-pressure automatic positive-displacement pump (PMP-0500−1−10-MB-316-M4-C0, DBR, Canada), and vice versa. Meanwhile, a video-based digital cathetometer with a resolution of 0.002 cm allows direct and accurate measurements of the injection volume. A highprecision test gauge (2089, Ashcroft) connected to the inlet tubing of the PVT cell has an accuracy of 0.05% of full scale 13 790 kPa (2000 psia) for saturation pressure measurement, and full scale 6895 kPa (1000 psia) for three-phase boundary measurement.

carbon no.

wt %

mol %

carbon no.

wt %

mol %

C1 C2 C3 C4 C5 C6 C7 C8 C9 C10 C11 C12 C13 C14 C15 C16 C17 C18 C19 C20 C21 C22 C23 C24 C25 C26 C27 C28 C29 C30 C31

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.786 0.844 1.037 1.267 1.817 1.950 2.300 2.143 2.286 2.238 2.048 1.857 2.071 1.329 1.743 1.571 1.714 1.600 1.583 1.650 1.452 1.281 1.390

0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 2.730 2.640 2.960 3.320 4.390 4.380 4.830 4.220 4.240 3.920 3.400 2.930 3.110 1.910 2.390 2.070 2.170 1.950 1.850 1.860 1.580 1.350 1.420

C32 C33 C34 C35 C36 C37 C38 C39 C40 C41 C42 C43 C44 C45 C46 C47 C48 C49 C50 C51 C52 C53 C54 C55 C56 C57 C58 C59 C60 C61+

1.329 0.981 1.083 1.307 1.343 0.917 0.900 1.583 1.600 0.750 0.830 1.380 0.929 0.929 0.700 0.833 0.792 0.792 0.733 0.750 0.717 0.717 0.683 0.650 0.650 0.667 0.667 0.667 0.767 35.397

1.310 0.940 1.010 1.180 1.180 0.780 0.750 1.290 1.270 0.580 0.630 1.020 0.670 0.650 0.480 0.560 0.530 0.510 0.460 0.470 0.440 0.430 0.400 0.370 0.370 0.370 0.360 0.360 0.410 14.600

2.3. Experimental Procedure. In this study, the solvent(s)−CO2−heavy oil systems with 15 feeds have been measured, each of which is listed in Table 2 in terms of both weight percentage and mole percentage at various temperatures. The saturation pressures and swelling factors for Feeds #1−12 are measured, whereas the L 1 L 2 V three-phase boundaries for Feeds #13−15 are visually determined. The measurements of saturation pressure, volume, temperature, and the determination of the mixture’s composition have uncertainties of ±7 kPa, ± 1.59 × 10−2 cm3, ± 0.10 K, and ±3.0 wt % respectively, whereas the measurement of threephase boundary has an uncertainty of ±4 kPa. The procedure of determining saturation pressure and swelling factor is briefly described as follows. As for preparation work, the solvents are compressed to liquid phase in the transfer cylinder. The PVT cell and the tubing are cleaned with kerosene, flushed with the pure solvent, and finally evacuated. The temperature of the air bath is set to the prespecified value at least for 12 h prior to the experimental measurement. Subsequently, the solvent is first injected into the PVT cell at a pressure higher than its vapor pressure by using the syringe pump (500HP, Teledyne ISCO Inc., USA). If two solvents are involved, the solvent with a lower vapor pressure is first introduced, whereas its injected volume can be obtained from the digital cathetometer. After turning on the magnetic mixer, the heavy oil is then injected to the PVT cell. The volume of the heavy oil introduced can be determined by subtracting the 1294


Energy & Fuels

Article

Figure 1. Schematic diagram of the experimental setup for conducting PVT measurements for the solvents−CO2−heavy oil systems.

Table 2. Compositions of the Solvent(S)−CO2−Lloydminster Heavy Oil Systems with 15 Feeds at Experimental Temperatures Composition, wt %

Composition, mol %

feed no.

C3H8

n-C4H10

CO2

heavy oil

C3H8

n-C4H10

CO2

heavy oil

temperature, K

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

11.4 19.8 29.3 11.5 14.8 20.9 15.8 0.0 0.0 14.2 10.3 0.0 0.0 12.3 0.0

0.0 0.0 0.0 0.0 0.0 0.0 0.0 16.4 36.2 5.3 0.0 10.1 0.0 0.0 10.1

0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 11.0 7.0 60.4 35.1 54.2

88.6 80.2 70.7 88.5 85.2 79.1 84.2 83.6 63.8 80.5 78.7 82.9 39.6 52.6 35.7

58.5 73.0 81.9 58.8 65.6 74.4 67.3 0.0 0.0 55.4 36.2 0.0 0.0 23.6 0.0

0.0 0.0 0.0 0.0 0.0 0.0 0.0 62.0 82.5 15.8 0.0 34.3 0.0 0.0 11.8

0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 38.6 31.7 94.4 67.2 83.2

41.5 27.0 18.1 41.2 34.4 25.6 32.7 38.0 17.5 28.8 25.2 34.0 5.6 9.2 5.0

323.85 323.85 323.85 298.85 298.85 298.85 298.85, 323.15, 348.15, 396.15 298.85, 323.15, 348.15, 396.15 298.85, 323.15, 348.15, 396.15 298.85, 348.15, 396.15 280.45, 298.85, 318.75 298.85, 318.75, 347.65, 391.55 288.65, 298.65, 304.95 288.95, 298.55 288.65, 298.65, 308.55

relative high pressure at which L1L2 phase equilibrium can be observed. As for each minor withdrawn of pressure, it is necessary to take 12 h to obtain a stable cell pressure prior to turning off the mixer. Subsequently, the phase equilibrium state could be visually observed to be L1L2, L1L2V, or L1V equilibrium. As for Feeds #14−15, an amount of liquefied C3H8 or n-C4H10 has to be added into the PVT cell before heavy oil is injected. The determination procedure of threephase boundary is similar to that of Feed #13. The detailed experimental procedures together with the corresponding phase boundaries and swelling factors can also be found elsewhere.8−10

solvent volume from the total volume. To achieve sufficient dissolution, the mixer is kept rotating for 6 h prior to further measurement. The mixture is then depressurized starting from a liquid state at a withdrawal rate of 3 cm3/h. Finally, saturation pressure can be read from the measured diagram of the P−V relation by finding the inflection point, and the swelling factor can be determined by measuring the height of the fluid. As for determining three phase boundaries, it is found that the excessive CO2 is the essential condition to maintain a system at the three-phase equilibrium. Accordingly, a CO2− heavy oil mixture (i.e., Feed #13), a n-C4H10−CO2−heavy oil mixture (i.e., Feed #14), and a C3H8−CO2−heavy oil mixture (i.e., Feed #15) are experimentally studied by performing the following procedure. After completing all of the preparation of cleaning, evacuating, and preheating temperature, a large amount of liquefied CO2 is introduced to the PVT cell followed by turning on the magnetic mixer. Then, the heavy oil sample is injected to form Feed #13. Once the mixture has been rocked for 6 h, depressurization in PVT cell starts at a

3. MATHEMATICAL FORMULATION 3.1. Heavy Oil Characterization. 3.1.1. Splitting the Plus Fraction. As shown in Table 1, no light components (i.e., C1− C8) are found in the Lloydminster heavy oil. The plus fraction (C61+) of the heavy oil sample accounts for a large portion with 14.6 mol %, which is higher than that of any lighter 1295


Energy & Fuels

Article

ponents are performed to examine the effect of the number of pseudocomponents on EOS prediction accuracy with respect to two BIP correlations, whereas a lumping method based on concentration and molecular weight of SCN is applied. Such a lumping method requires that summation of the mole fractions times the logarithm of the molecular weight (ΣzilnMi) should be the same for each pseudocomponent.29 The remaining work is to assign properties to the lumped pseudocomponent based on properties of each SCN. Hong’s mixing rule is employed as follows,53

components. Because it is previously found that the component as heavy as C200 may also affect the phase behavior of heavy oil involved systems, it is necessary to further split the plus fraction into single carbon numbers (SCNs).41 The exponential distribution function proposed by Pedersen et al.23 is used in this study. It was found that, for heavy oil, the carbon number starting at C11 or higher exhibits an approximate linear relationship between carbon number and the logarithm of its corresponding mole fraction, that is,23 ln zi = A + B × CNi

(1)

u

where zi and CNi are mole fraction and carbon number of the component i, respectively. A and B are constants determined by fitting the experimental data. The molecular weight (M) of each extended SCN can be calculated from the following equation,24

Mi = 14CNi − 4

ηj =

(2)

CNmax

∑

zi (3)

i = CN+

p=

CN

M+ =

max ∑i = CN ziMi +

+

(5)

RT a − VM − b VM(VM + b) + b(VM − b)

a = acα(Tr , ω)

CN

max ∑i = CN zi

u

∑i = l ziMi

where ηj can be any property (e.g., Tb, Tc, Pc, vc, ZRA, M, γ, or ω) for the jth lumped pseudocomponent containing SCNs from l to u; zi and Mi are the mole fraction and molecular weight of the ith SCN. 3.2. PR EOS Model. 3.2.1. PR EOS. The PR EOS54 is one of the most widely used EOS in petroleum industry. It is formulated as follows,

It should be noted that eq 1 needs to satisfy the following mass balance equations:

z+ =

∑i = l ziMiηi

(4)

where z+ and M+ are mole fraction and molecular weight of the plus fraction. The Gamma-distribution model is also commonly used to split the plus fraction;26−28 nevertheless, the parameter defining the form of the distribution can be set to unity. As such, the Gamma-distribution model can be easily transformed into an exponential distribution function.24 To compare the above two widely used characterization methods, it is found that the exponential distribution function proposed by Pedersen et al. usually provides a better prediction.31 3.1.2. Determination of SCN Properties. Critical properties of any individual component need to be determined prior to performing the associated EOS calculations, though they are not usually available from compositional analysis reports. Numerous empirical correlations have been proposed to estimate the critical properties, most of which are a function of specific gravity and the boiling point.42−45 The correlations used in this study to compute the properties of each SCN are presented in Appendix. By assuming a constant Watson factor (Kw)46 for each SCN, specific gravity, γ, can be calculated.47 The boiling point temperature (TbR) is estimated from molecular weight and specific gravity by using the Soreide correlation48 as suggested by Whitson et al.47 The critical properties and acentric factor of SCN are then determined with the Kesler-Lee correlations42,49 as recommended by several researchers.14,32,50 Furthermore, the critical volume vc is calculated from the commonly used Twu correlations,51 which needs the critical properties of normalparaffin hydrocarbon as the precondition. Also, the Rackett parameter52 (ZRA), which is required input in volume shift calculations, is calculated by applying the corresponding equation in Appendix at 60 °F. 3.1.3. Lumping SCNs. It is necessary to group SCNs into appropriate numbers of pseudocomponents while preserving the prediction accuracy of the PR EOS model. In general, five to eight pseudocomponents are considered adequate for simulation purposes.47 In this study, two to eight pseudocom-

ac = 0.457235

(6a) (6b)

R2Tc2 pc

(6c)

RTc pc

(6d)

b = 0.0777969

where p is pressure, T is temperature, VM is molar volume. Recently, a modified alpha function, which can more accurately predict the vapor pressure of both nonhydrocarbon and light−heavy hydrocarbon pure substance, is proposed as follows,55 α = exp{(0.13280 − 0.05052ω + 0.25948ω 2)(1 − Tr) + 0.81769ln[1 + (0.31355 + 1.86745ω − 0.52604ω 2)(1 −

Tr )]2 }

(6e)

For mixtures, the parameters a and b are calculated with the following van der Waals’ mixing rule,56 a=

∑ ∑ yyi j (aiaj)0.5 (1 − δij) i

b=

j

∑ yj bj j

(7)

(8)

where yi and yj are the compositions of the ith and jth substance, respectively, and δij is BIP between the ith and the jth component. 3.2.2. BIP Correlations. In the literature, numerous efforts have been made to develop BIP correlations for mixtures.8,34,36,37,57−59 Two BIP correlations are widely used when an oil sample is treated as multipseudocomponents. This is attributed to the fact that the BIP matrix can be regenerated by adjusting only one parameter. BIP Correlation #134 is dependent on critical volumes, 1296


Energy & Fuels β ⎡ 1/3 1/3 ⎤ ⎢ 2 Vci Vcj ⎥ δij = 1 − ⎢ 1/3 ⎥ V + Vc1/3 j ⎥ ⎣⎢ ci ⎦

Article n

V* = V −

θ=

(9)

where V* is the corrected molar volume, V is the original molar volume calculated by the PR EOS, Ci is the correction factor for the ith component, and xi is the composition of the ith component. The following dimensionless shift parameter (s) is usually used to relate C with the b term in eq 6d of PR EOS,

s=

(10a)

s=1−

V2 V1(1 − S)

(13a)

d Me

(13b) 39

where d and e are 2.258 and 0.1823 for n-alkanes. In addition, the Rackett’s compressibility factor ZRA is involved in the following two correlations proposed by Peneloux et al.38 and Twu et al.40 to estimate the correction factor C, respectively.

(10b)

where Zc is the critical compressibility factor of a component. The exponent θ can be adjusted to fit experimental data. BIP Correlations #1 and #2 are applied in this study to predict the saturation pressures of solvent(s)−CO2−heavy oil systems with 12 feeds (Feeds #1−12 in Table 2), respectively. The BIPs between pseudocomponents of heavy oil are set to zero due to their similar properties, while the BIPs for CO2− C3H8 pair and CO2−n-C4H10 pair are set as 0.135 and 0.130, respectively.61 Moreover, the BIPs for C3H8−n-C4H10 pair can also be calculated by eqs 9 and 10, respectively. As for alkane solvent−heavy oil systems (i.e., Feeds #1−10 in Table 2), the exponent β or θ for pairs of alkane solvent and pseudocomponent is tuned to fit the measured saturation pressures. Once the exponent β or θ is optimized based on the measured saturation pressures for Feeds #1−10, they are then directly used to calculate the BIPs for solvent-pseudocomponent pairs (i.e., Feeds #11−12 in Table 2). Meanwhile, the exponents in the two BIP correlations for CO2-pseudocomponent pair are to be still tuned to match the experimental saturation pressure data. 3.2.3. Swelling Factor Calculation. Swelling factor is an important parameter for designing solvent-based processes for heavy oil recovery purpose. Swelling effect due to solvent dissolution in heavy oil generally can increase the in situ saturation of heavy oil in porous media resulting in an enhanced mobility of the swollen heavy oil. It is well recognized that a larger swelling factor is favorable to achieve a higher oil recovery during the solvent flooding process. The swelling factor can be calculated with the following equation,62 SF =

C b

Jhaveri et al. proposed the following equation to estimate s from molecular weight for heavy hydrocarbons,39

Zci + Zcj 2

(12)

i=1

where δij is the BIP between the ith and jth component, Vc is critical molar volume in m3/kmol, and β is an adjustable parameter. β = 1.2 for paraffin−paraffin interaction coefficients as recommended by Oellrich et al.60 BIP Correlation #236 is expressed as a function of critical temperature in the following form, ⎡ 2 T T ⎤θ ci cj ⎥ δij = 1 − ⎢ ⎢⎣ Tci + Tcj ⎥⎦

∑ xiCi

⎛ RT ⎞ C = 0.40768⎜⎜ c ⎟⎟(0.29441 − Z RA ) ⎝ pc ⎠

(14)

and ⎛ RT ⎞ C = 0.406501⎜⎜ c ⎟⎟(0.260484 − Z RA ) ⎝ pc ⎠

(15)

The three volume translation methods mentioned above are evaluated in terms of their accuracy for predicting the swelling factor. In this study, the absolute average relative deviation (AARD) is used as the objective function to optimize saturation pressure and swelling factor, AARD =

1 n

n

∑ i=1

Xicalcd − Xiexpt Xiexpt

(16)

where n is the number of data points, Xcalcd is the calculated i saturation pressure or swelling factor, and Xexpt is the measured i saturation pressure or swelling factor. As shown in Figure 2, a flowchart has been developed to illustrate the methodology adopted in this study to accurately model the phase behavior of solvent(s)−CO2−heavy oil mixtures by using the PR EOS model. The modeling procedure can be briefly described as follows: (1) Split the plus fraction of heavy oil and lump SCNs to a given number of pseudocomponents; (2) Assign physical and critical properties to each pseudocomponent based on the widely used characterization correlations; (3) Optimize the two exponents in eqs 9 and 10 for alkane solvent-pseudocomponent pair to match the measured saturation pressures of Feeds #1−10, respectively; (4) Optimize the two exponents in eqs 9 and 10 for CO2pseudocomponent pair to match the measured saturation pressures of Feeds #11−12, respectively. It should be noted that the optimized exponents in Step #3 are applied to calculate the BIPs for alkane solventpseudocomponent pair;

(11)

where SF is the swelling factor, V1 is the molar volume of the heavy oil at saturation temperature and atmospheric pressure (i.e., 101.3 kPa), V2 is the molar volume of solvent(s)− saturated heavy oil at saturation temperature and saturation pressure, and S is the solvent solubility in mole fraction. In this study, volume translation techniques38 are applied to improve the swelling factor prediction from PR EOS. It is noted that volume translation technique imposes no effect on the calculation of saturation pressures.38 The corrected molar volume can be readily calculated via the following equation,38 1297


Energy & Fuels

Article

Finally, the phase boundary at the fixed temperature is located on the pressure where the phase number changes.

4. RESULTS AND DISCUSSION 4.1. Heavy Oil Characterization. Once the molecular weight of the heaviest SCN becomes smaller than its previous one, the splitting procedure of plus fraction can be terminated.65 In this study, the plus fraction of heavy oil can be extended to the heaviest component C105+. The calculated molecular weight of C104 and C105+ are 1452.00 and 1457.76, respectively. As shown in Figure 3, the linear relationship

Figure 3. Molar distribution of the Lloydminster heavy oil.

between natural logarithm of SCN’s mole fraction and carbon number starts at C14. The two constants A and B in eq 1 regressed from the measured mole fractions of SCNs with carbon number between 14 and 60, are calculated to be −2.3985 and −0.0580 with a standard error of 0.0781 and 0.0020, respectively. As such, the mole fractions of C61−C104 can be easily determined by extrapolation. It is noted that the mole fraction and molecular weight of the heaviest component C105+ are determined by applying the mass balance equations (i.e., eqs 3 and 4), respectively. Figure 4 shows variation of the calculated critical temperature, boiling point temperature, acentric factor, and critical pressure with carbon number, respectively. The constant Watson factor Kw is calculated to be 12.04707. Consequently,

Figure 2. Flowchart for determining two- and three-phase boundaries as well as swelling factors of solvent(s)−CO2−heavy oil systems with 15 feeds.

(5) Steps #3 and 4 are performed for each number of pseudocomponents with consideration of the two BIP correlations, respectively; (6) Calculate the L1L2V three-phase boundaries for Feeds #13−15 with the optimized number of pseudocomponents and the optimized BIP correlation; (7) Evaluate three volume translation strategies in terms of the swelling factor prediction for Feeds #1−12 with the selected number of pseudocomponents and the optimized BIP correlation. In this study, the phase boundaries are computed by performing both Michelsen’s stability test 63 and flash calculation47 at various pressures and a fixed temperature. More specifically, multiple initial estimates are used to locate stationary points of the Gibbs tangent plan distance (TPD) function. If the Gibbs energy of a test mixture can be decreased by forming a new phase, the flash calculation should be further implemented. In this process, the multiple initialization corresponding to the multiple (pseudo)components in flash calculation can be obtained from the previous unstable solution associated with the stability test.47,64 The stability test could indicate that the mixture is in the state of either single-phase, or two-phase, and three-phase at one initiate pressure and a given temperature. Then, the pressure is increased by 1.0 kPa while keeping the temperature constant in each stability test iteration.

Figure 4. Variation of Tc, Pc, ω, and Tb of SCN fraction with carbon number for the Lloydminster heavy oil. 1298


Energy & Fuels

Article

Table 3. Physical and Critical Properties of Each Pseudocomponent in Different Lumping Schemes no. of PCsa

z mol %

M

pc, kPa

Tc, K

ω

γ

Tb, K

Vc, m3/kmol

ZRA

2

55.61 44.39 37.63 33.61 28.76 29.47 26.14 23.75 20.64 20.42 23.54 20.24 18.95 16.85 20.42 17.21 17.98 16.41 14.75 13.23 16.04 13.43 14.49 15.45 14.65 13.81 12.13 16.04 13.43 11.56 12.41 12.18 11.94 11.59 10.85

248.4 1005.5 204.4 375.6 1139.3 188.6 291.5 492.4 1248.9 171.2 245.6 356.0 577.7 1317.0 171.2 233.0 307.9 436.4 694.5 1396.5 162.1 212.4 259.8 336.9 467.1 741.5 1422.9 162.1 212.4 253.5 308.5 391.7 526.6 833.9 1450.6

1629.22 1060.48 1824.77 1242.19 1041.68 1921.40 1418.76 1092.30 1045.38 2048.14 1580.16 1258.00 1035.01 1053.21 2048.14 1633.05 1364.95 1136.15 995.86 1066.48 2123.91 1739.30 1515.98 1292.90 1103.37 989.31 1071.89 2123.91 1739.30 1539.29 1359.07 1190.43 1054.22 984.56 1078.27

752.84 1045.12 709.48 855.78 1080.74 691.15 797.28 920.58 1104.20 669.16 756.37 846.16 957.25 1116.75 669.16 744.08 811.45 894.83 997.24 1129.64 656.87 721.99 770.67 833.51 910.88 1010.49 1133.49 656.87 721.99 764.89 812.54 870.13 938.90 1033.57 1137.18

0.747 1.166 0.647 0.972 1.194 0.606 0.848 1.091 1.202 0.557 0.753 0.957 1.145 1.202 0.557 0.724 0.881 1.051 1.189 1.198 0.531 0.673 0.786 0.931 1.079 1.199 1.196 0.531 0.673 0.772 0.884 1.006 1.124 1.212 1.194

0.8487 1.0791 0.8213 0.9159 1.1114 0.8100 0.8766 0.9621 1.1346 0.7965 0.8504 0.9088 0.9903 1.1480 0.7965 0.8427 0.8857 0.9427 1.0236 1.1625 0.7889 0.8290 0.8594 0.9002 0.9543 1.0354 1.1670 0.7889 0.8290 0.8557 0.8863 0.9250 0.9752 1.0570 1.1716

579.66 865.64 533.68 688.69 897.89 514.44 626.64 756.13 917.59 491.56 582.96 678.75 792.99 927.49 491.56 569.84 641.79 729.94 831.43 937.02 478.89 546.41 598.19 665.33 746.55 843.56 939.70 478.89 546.41 592.01 642.97 704.18 775.16 863.95 942.12

0.9060 1.3586 0.7878 1.1649 1.3847 0.7376 1.0272 1.2877 1.3922 0.6789 0.9158 1.1501 1.3363 1.3936 0.6789 0.8813 1.0657 1.2507 1.3714 1.3934 0.6470 0.8191 0.9560 1.1215 1.2783 1.3789 1.3929 0.6470 0.8191 0.9400 1.0694 1.2050 1.3198 1.3887 1.3923

0.2465 0.3076 0.2489 0.2447 0.3220 0.2502 0.2438 0.2499 0.3350 0.2521 0.2454 0.2435 0.2562 0.3435 0.2521 0.2461 0.2433 0.2462 0.2670 0.3537 0.2531 0.2476 0.2445 0.2432 0.2478 0.2720 0.3572 0.2531 0.2476 0.2448 0.2432 0.2441 0.2516 0.2823 0.3609

3

4

5

6

7

8

a

PC is the abbreviation of pseudocomponent.

Table 4. Optimal Exponents in Two BIP Correlations for Each Lumping Scheme no. of PCs β θ

alkane solvent−PC pair in Systems #1−10 CO2−PC pair in Systems #11−12 alkane solvent−PC pair in Systems #1−10 CO2−PC pair in Systems #11−12

2

3

4

5

6

7

8

0.918 1.407 0.374 0.811

0.973 1.468 0.362 0.873

0.980 1.500 0.365 0.901

0.985 1.514 0.368 0.912

0.992 1.521 0.372 0.919

0.990 1.524 0.372 0.921

0.991 1.529 0.372 0.923

2−8 pseudocomponents are grouped with the Hong’s mixing rule (eq 5). The calculated physical and critical properties of each pseudocomponent in different lumping schemes (i.e., 2−8 pseudocomponents) are tabulated in Table 3. It should be noted that the value of ΣziInMi for each lumped group cannot be equal to one another perfectly due to the discontinuity of carbon numbers. 4.2. Performance of BIP Correlations. As for each lumping scheme, the exponent in either eq 9 or eq 10 is tuned twice, that is, one for alkane solvent-pseudocomponent pair based on the saturation pressure data of Feed #1−10, and the other for CO2−pseudocomponent pair based on the saturation pressure data of Feed #11−12. As a result, the optimal exponents in each lumping scheme are listed in Table 4. To account for the large physical and chemical difference between

CO2 and hydrocarbon pseudocomponent, it can be seen that the exponents for CO2−pseudocomponent pair are far larger than that for the alkane solvent−pseudocomponent pair in all lumping schemes. The AARD in predicting saturation pressure by applying the previously optimized exponents can be plotted as a function of the number of pseudocomponents (Figure 5). Part a of Figure 5 plots the calculation accuracy for alkane solvent(s)−heavy oil systems (i.e., Feeds #1−10), solvent(s)−CO2−heavy oil systems (i.e., Feeds #11−12), and all systems (i.e., Feeds #1−12), respectively. In general, the AARD decreases with the increased number of pseudocomponents. The minimum AARD of 5.06% is obtained with eight pseudocomponents for all the systems. As for BIP Correlation #2 (part b of Figure 5), the prediction accuracy for both Feeds #1−10 and overall Feeds 1299


Energy & Fuels

Article

As indicated in flowchart, the optimum scheme will be further used to calculate the swelling factors and three-phase boundaries. Although eight pseudocomponents with BIP Correlation #1 produces the most accurate results, AARD is improved only slightly when more than five pseudocomponents are applied (part a of Figure 5). With consideration of saving computing time without impairing the prediction accuracy, six pseudocomponents with BIP Correlation #1, which generates the total AARD of 5.07%, is selected to perform further calculations hereafter. 4.3. Evaluation of Volume Translation Methods. The molar volume (V1) of the heavy oil at saturation temperature and atmospheric pressure (i.e., 101.3 kPa) can be obtained by conducting single phase calculation, while the molar volume (V2) of solvent(s)-saturated heavy oil at saturation temperature and saturation pressure can be obtained when conducting the saturation pressure calculation. It should be noted that the involvement of correction factor imposes no impacts on phase equilibrium calculation results. As previously mentioned, the volume shifts used in EOS calculation are calculated by three strategies in this study, that is, Jhaveri et al. method,39 Peneloux et al. method,38 and Twu et al. method.40 The calculated volume shifts as well as AARDs of the calculated swelling factor can be found in Table 5. Overall, the Peneloux et al. correlation38 is found to be the most accurate for predicting swelling factors with an overall AARD of 2.09% in comparison to the other two methods. It is worthwhile noting that the Twu et al. method40 also has a strong prediction capacity with overall AARD of 2.41%. 4.4. Saturation Pressures and Swelling Factors. 4.4.1. C3H8−Heavy Oil System with Feeds #1−7. The exponent β in BIP Correlation #1 is found to be 0.992 for alkane solvent−pseudocomponent pair as shown in Table 4 when the heavy oil is characterized with six pseudocomponents. Figure 6 presents the measured and predicted saturation pressures and swelling factors as a function of mole fraction of C3H8 at the temperature of 323.85 and 298.85 K, respectively. As can be seen, the saturation pressure increases with an increase in the mole percentage of C3H8 at the same temperature, whereas it increases with temperature for a C3H8−heavy oil mixture. Similar findings can be observed for Feed #7 (Figure 7). Overall, the PR EOS model with the modified alpha function is able to accurately predict the saturation pressure with an AARD of 2.82% for the ten experimental data points in Figures 6 and 7. The prediction

Figure 5. Variation of AARD% versus the number of pseudocomponents with the optimal exponents in (a) BIP Correlation #1 and (b) BIP Correlation #2.

#1−12 is found not to be improved by increasing the number of pseudocomponents with an unexpectedly minimum AARD of 5.52% obtained from two pseudocomponents, though monotone decreasing trend is observed for Feeds #11−12. Overall, BIP Correlation #1 yields more accurate results in any lumping scheme than BIP Correlation #2 for predicting the saturation pressure of all feeds.

Table 5. AARD% of the Predicted Swelling Factors and the Calculated Volume Shifts with Three Volume Shift Strategies Volume Shift, s Jhaveri et al. PC

Solvent

AARD%

#1 #2 #3 #4 #5 #6 CO2 C3H8 n-C4H10 Feeds #1−10 Feeds #11−12 Overall Feeds #1−12

0.11575 0.16407 0.13002 0.09814 0.07719 0.07298 −0.09435 −0.07330 −0.05706 4.96 1.14 4.00 1300

39

Peneloux et al.38

Twu et al.40

0.22196 0.25326 0.26782 0.25289 0.14357 −0.31068 0.10905 0.09490 0.11324 2.00 2.38 2.09

0.04405 0.07526 0.08978 0.07489 −0.03411 −0.48705 −0.06853 −0.08264 −0.06435 1.95 3.79 2.41


Energy & Fuels

Article

Figure 8. Comparison of the measured and predicted saturation pressures and swelling factors for binary n-C4H10−heavy oil system with Feed #8 (62.0 mol % n-C4H10, 38.0 mol % heavy oil) and Feed #9 (82.5 mol % n-C4H10, 17.5 mol % heavy oil).

Figure 6. Comparison of the measured and predicted saturation pressures and swelling factors for binary C3H8−heavy oil system with Feeds #1−3 at 323.85 K and with Feeds #4−6 at 298.85 K, respectively.

14.57% is found at 396.15 K for Feed #9, which might be caused by experimental uncertainty. In addition, the dissolution of n-C4H10 in heavy oil can lead to a large swelling factor. 4.4.3. C3H8−n-C4H10−heavy oil system with Feed #10. Figure 9 shows the measured saturation pressures and swelling

Figure 7. Comparison of the measured and predicted saturation pressures and swelling factors for binary C3H8−heavy oil system with Feed #7 (67.3 mol % C3H8, 32.7 mol % heavy oil).

accuracy is slightly lower at a higher temperature close to the critical temperature of C3H8. In addition, it is found that the dissolution of C3H8 in heavy oil can lead to a large swelling factor indicating that the increment of in situ oil saturation is favorable to enhance oil recovery. In general, the PR EOS model in conjunction with the Peneloux et al. volume translation38 provides accurate prediction of the swelling factors yielding an AARD of 1.38% for the ten measured data points. 4.4.2. n-C4H10−Heavy Oil System with Feeds #8−9. Figure 8 depicts the measured saturation pressures and swelling factors for the n-C4H10−heavy oil system. Obviously, a higher saturation pressure and a larger swelling factor are measured for Feed #9 compared to Feed #8. This is because more nC4H10 is dissolved into heavy oil for Feed #9. It can also be seen that the saturation pressure increases with temperature for a given feed, implying that solubility of n-C4H10 in heavy oil at the same pressure is lower at the higher temperature. As for the eight saturation pressures measured for Feeds #8−9, its corresponding predicting accuracy by using the PR EOS model has an AARD of 6.80%. Furthermore, the swelling factors for n-C4H10−heavy oil systems are accurately predicted with an AARD of 2.78%. One exception with an AARD of

Figure 9. Comparison of the measured and predicted saturation pressures and swelling factors for ternary C3H8−n-C4H10−heavy oil system with Feed #10 (55.4 mol % C3H8, 15.8 mol % n-C4H10, 28.8 mol % heavy oil).

factors for the ternary C3H8−n-C4H10−heavy oil system. As can be seen, the results calculated from the model generally lies above the measured data in a high temperature range, implying the solubility of solvents in heavy oil is slightly underestimated at a given pressure. Although saturation pressure prediction for the ternary system is not as good as that for the aforementioned binary systems, it still yields a reasonable AARD of 8.70% for three data points. Furthermore, there exists a good agreement between the measured and predicted swelling factors with an AARD of only 1.98% by applying the Peneloux et al. volume translation strategy. It should be noting that the BIP between C3H8 and n-C4H10 is calculated to be 0.00087 by setting β = 1.2 for paraffin−paraffin interaction coefficients in BIP Correlation #1. 4.4.4. C3H8/n-C4H10−CO2−Heavy Oil Systems with Feeds #11−12. Ternary systems in the presence of CO2 are then 1301


Energy & Fuels

Article

studied. The BIPs between C3H8 and CO2, n-C4H10 and CO2 are set to 0.135 and 0.130, respectively.62 The optimized exponent β = 0.992 in BIP Correlation #1 is directly used for calculating BIP between alkane solvent−pseudocomponent pair, whereas the exponent β is tuned for the CO2 and each pseudocomponent pair. As shown in Table 4, the optimum β for CO2−pesudocomponent pair is found to be 1.521 when heavy oil is split into six pseudocomponents. As shown in Figure 10, the saturation pressure increases with an increase in

Figure 10. Comparison of the measured and predicted saturation pressures and swelling factors for ternary C3H8−CO2−heavy oil system with Feed #11 (36.2 mol % C3H8, 38.6 mol % CO2, 25.2 mol % heavy oil) and ternary n-C4H10−CO2−heavy oil system with Feed #12 (34.3 mol % n-C4H10, 31.7 mol % CO2, 34.0 mol % heavy oil).

temperature for both ternary systems. The predicted saturation pressures agree well with the measured values with an AARD of 4.67%. A larger deviation occurs in predicting the saturation pressure at the highest temperature for Feed #12 due to the fact that the system lies in the supercritical region. In addition, the swelling factor is found to be weakly dependent on temperature for a given mixture. Taking Feed #12 as an example, the swelling factor increases slightly from 1.315 to 1.364 when temperature is increased from 298.85 to 391.55 K. As a result, the swelling factor can be matched accurately by using the PR EOS model together with Peneloux et al. volume translation method with an AARD of 2.38% for seven data points. 4.5. L1L2V Three-Phase Boundaries. Experiments are also conducted to determine the three-phase boundaries for Feeds #13−15. In a typical L1L2V phase equilibrium, the L1 phase denotes the high-density liquid phase which consists of solvent(s) and heavy oil; the L2 phase denotes the lowerdensity liquid phase that mainly contains liquid solvents; and the V phase denotes the vapor phase composed of gaseous solvents. It is worthwhile noting that the optimized BIPs from Feeds #1−12 are employed directly to calculate L1L2V threephase boundary. It is possible to have a relatively high AARD between the measured and calculated phase boundaries, compared with the previously predicted saturation pressure. Experimentally, the three-phase boundaries of a binary CO2− heavy oil system with Feed #13 are observed at 288.65 and 298.65 K, while there is no three-phase equilibrium visually and qualitatively observed at 304.95 K. Theoretically, the model fails to locate three-phase boundary for Feed #13 in this study. As such, parts a and b of Figure 11 depict the measured L1L2V boundaries and the calculated pressure−composition (P−x) phase diagram for Feed #13 at 288.65 and 298.65 K,

Figure 11. Predicted P−x phase diagram for CO2−heavy oil system and the measured L1L2V phase boundaries for CO2−heavy oil system with Feed #13 (94.4 mol % CO2, 5.6 mol % heavy oil) at (a) 288.65 K and (b) 298.65 K.

respectively. It can be seen that at the given temperature, the measured upper three-phase boundary and lower three-phase boundary are very close to each other indicating an extremely narrow three phase zone. Besides, the measured boundaries are also close to the predicted transition pressure between L1−L1V phase boundary and L1L2−L1V phase boundary. As shown in Figure 12, although the predicted three-phase boundaries are narrower than the measured data for C3H8− CO2−heavy oil systems with Feed #14, there exists an excellent prediction for the lower boundary with an overall AARD of 2.34% for the four data points. Figure 13 compares the measured and predicted three-phase boundaries for n-C4H10− CO2−heavy oil systems with Feed #15. The lower boundary calculated from the model is overestimated, while the upper prediction is accurate with an overall AARD of 6.08% for six data points. In addition, it is found that the predicted threephase region is also smaller than the experimental measurements elsewhere.66

5. CONCLUSIONS A generalized methodology has been proposed to determine liquid−vapor (L1V) and liquid−liquid−vapor (L1L2V) phase boundaries as well as swelling factors of solvent(s)−CO2− heavy oil systems at high pressures and elevated temperatures. 1302


Energy & Fuels

■

Article

APPENDIX: CORRELATIONS OF SCN PROPERTIES

(1). Specific Gravity46,47

γ = 6.0108M 0.17947K w−1.18241

⎡ 0.16637γ (∑CNmax z M 0.82053) ⎤−0.84573 + i = CN+ i i ⎥ Kw = ⎢ ⎢ ⎥ z M + + ⎣ ⎦ (2). Boiling Point Temperature (Soreide)48

TbR = 1928.3 − (1.695 × 105)M −0.03522γ 3.266 × exp[−(4.922 × 10−3)M − 4.7685γ + (3.462 × 10−3)Mγ ] (3). Critical Temperature and Pressure (Kesler-Lee)42

Figure 12. Comparison of the measured and predicted L1L2V threephase boundaries for the ternary C3H8−CO2−heavy oil system with Feed #14 (23.6 mol % C3H8, 67.2 mol % CO2, 9.2 mol % heavy oil).

TcR = 341.7 + 811γ + (0.4244 + 0.1174γ )TbR −1 + (0.4669 − 3.2623γ ) × 105TbR

pcpsi = exp{8.3634 − 0.0566γ −1 − [(0.24244 + 2.2898γ −1 + 0.11857γ −2) × 10−3]TbR 2 + [(1.4685 + 3.648γ −1 + 0.47227γ −2) × 10−7]TbR 3 − [(0.42019 + 1.6977γ −2) × 10−10]TbR }

(4). Acentric Factor (Lee-Kesler)49

ω=

−1 6 − ln(pcpsi /14.7) + A1 + A 2Tbr + A3 ln Tbr + A4 ln Tbr −1 6 + A 7 ln Tbr + A8 ln Tbr A5 + A 6Tbr

if Tbr = TbR /TcR < 0.8

ω = −7.904 + 0.1352K w − 0.007465K w2 + 8.359Tbr Figure 13. Comparison of the measured and predicted L1L2V threephase boundaries for the ternary n-C4H10−CO2−heavy oil system with Feed #15 (11.8 mol % n-C4H10, 83.2 mol % CO2, 5.0 mol % heavy oil).

−1 + (1.408 − 0.01063K w )Tbr

if Tbr = TbR /TcR > 0.8

where A1 = −5.92714, A2 = 6.09648, A3 = 1.28862, A4 = −0.169347, A5 = 15.2518, A6 = 15.6875, A7 = −13.4721, and A8 = 0.43577.

The modified alpha function is successfully incorporated in PR EOS when the heavy oil is characterized as muti-pseudocomponents because the modified alpha function is widely suitable for both light and heavy hydrocarbon compounds. By separately tuning the BIPs for mixtures with and without nonhydrocarbon (i.e., CO2), a good performance of predicting saturation pressures and three phase boundaries has been obtained with AARDs of 5.07% and 4.58%, respectively. It is also found that the prediction accuracy increases with the increased number of pseudocomponents by employing the BIP correlation, which is a function of critical volume. However, six pseudocomponents is demonstrated to be sufficient to characterize the heavy oil, which is attributed to the fact that the heavier fractions of the heavy oil are scarcely involved in mass transfer process. Finally, the Peneloux et al. method for calculating volume shift is found to be the optimal method to predict swelling factors with an AARD of 2.09%, compared with the Jhaveri et al. method and Twu et al. method.

(5). Critical Volume (Twu)51

⎛ 1 + 2f ⎞2 v ⎟ vc = vcP⎜⎜ ⎟ ⎝ 1 − 2fv ⎠

⎡ 0.46659 ⎛ 3.01721 ⎞ ⎤ ⎟Δγv ⎥ ⎜ fv = Δγv ⎢ + − + 0.182421 0.5 0.5 ⎢⎣ TbR TbR ⎠ ⎥⎦ ⎝

Δγv = exp[4(γP2 − γ 2)] − 1 (6). Critical Properties of the Normal-Paraffin51

⎡ 0.533272 + 0.191017 × 10−3T ⎤−1 bP ⎢ ⎥ −7 2 ⎥ TcP = TbP⎢ + 0.779681 × 10 TbP− ⎢ ⎥ −10 3 −13 ⎢⎣ 0.284376 × 10 TbP + 95.9468(0.01TbP) ⎥⎦ 1303


Energy & Fuels

Article

⎡ ⎤2 ⎛ ⎞0.5 ⎢ 3.83354 + 1.19629⎜1 − TbP ⎟ ⎥ ⎢ ⎥ TcP ⎠ ⎝ ⎢ ⎥ ⎛ ⎞ T ⎢ + 34.8888⎜1 − bP ⎟+ ⎥ pcP = ⎢ ⎥ TcP ⎠ ⎝ ⎢ ⎥ 2 4⎥ ⎢ ⎛ ⎞ ⎛ ⎞ ⎢ 36.1952⎜1 − TbP ⎟ + 104.193⎜1 − TbP ⎟ ⎥ ⎢⎣ TcP ⎠ TcP ⎠ ⎥⎦ ⎝ ⎝

d, e = positive correlation coefficients defined in eq 13b f v = parameter defined in Twu critical volume correlation Kw = Watson factor l = lower boundary of single carbon number for a pseudocomponent Mi = molecular weight of the component i, g/mol M+ = molecular weight of the plus fraction MP = molecular weight of normal-paraffin hydrocarbons, lbm/lbm mol n = number of data points p = pressure, kPa pc = critical pressure, kPa pcP = critical pressure of normal-paraffin hydrocarbon, psia pcpsi = critical pressure, psia pST = saturation pressure, kPa R = universal gas constant s = volume shift S = solvent solubility in mole fraction SF = swelling factor T = temperature, K Tb = normal boiling point at 1 atm, K Tbr = reduced normal boiling point TbP = normal boiling point at 1 atm of normal-paraffin hydrocarbon, °R TbR = normal boiling point at 1 atm, °R Tc = critical temperature, K TcP = critical temperature of normal-paraffin hydrocarbon, °R TcR = critical temperature, °R Tr = reduced temperature u = upper boundary of single carbon number for a pseudocomponent vc = critical molar volume, ft3/lbm mol vcP = critical molar volume of normal-paraffin hydrocarbons, ft3/lbm mol V = original molar volume calculated by the PR EOS, m3/kmol Vc = critical molar volume, m3/kmol VM = molar volume, m3/kmol V1 = molar volume of the heavy oil at saturation temperature and atmospheric pressure, m3/kmol V2 = molar volume of solvent(s)-saturated heavy oil at saturation temperature and saturation pressure, m3/kmol V* = corrected molar volume, m3/kmol xi = liquid-phase composition vector Xcalcd = calculated objective parameter i Xexpt = measured objective parameter i yi = composition vector zi =mole fraction of the component i z+ = mole fraction of the plus fraction Zci = critical compressibility factor of the component i ZRA = Rackett parameter

⎡ ⎛ ⎛ T ⎞ vcP = ⎢1 − ⎜⎜0.419869 − 0.505839⎜1 − bP ⎟ ⎢ TcP ⎠ ⎝ ⎝ ⎣ −8 3 14 ⎞⎤ ⎛ ⎛ TbP ⎞ TbP ⎞ ⎟⎥ − 1.56436⎜1 − ⎟ − 9481.7⎜1 − ⎟ TcP ⎠ TcP ⎠ ⎟⎠⎥⎦ ⎝ ⎝

⎛ T ⎞ γP = 0.843593 − 0.128624⎜1 − bP ⎟ − 3.36159 TcP ⎠ ⎝ 3 12 ⎛ ⎛ TbP ⎞ TbP ⎞ ⎜1 − ⎟ − 13749.5⎜1 − ⎟ TcP ⎠ TcP ⎠ ⎝ ⎝

⎡ 5.71419 + 2.71579ln M − 0.28659(ln M )2 −⎤ P P ⎥ TbP = exp⎢ ⎢⎣ 39.8544(ln M )−1 − 0.122488(ln M )−2 ⎥⎦ P P − 24.7522ln MP + 35.3155(ln MP)2 MP ≈

TbP 10.44 − 0.0052TbP

(7). Rackett Parameter52 2/7

Z RA

■

⎛ Mpc ⎞1/[1 + (1 − Tr) =⎜ ⎟ ⎝ RTcγ ⎠

]

AUTHOR INFORMATION

Corresponding Author

*Phone: 1-306-337-2660, fax: 1-306-585-4855, e-mail: tony. [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors acknowledge a Discovery Grant and a CRD Grant awarded to Dr. Yang from the Natural Sciences and Engineering Research Council (NSERC) of Canada and EHR Enhanced Hydrocarbon Recovery Inc. for financial support.

■

Greek Letters

NOMENCLATURE a = attraction parameter in PR-EOS ac = PR EOS constant A = constant defined in eq 1 AARD = absolute average relative deviation b = van der Waals volume, m3/kmol B = constant defined in eq 1 Ci = volume correction factor for the ith component, m3/kmol CNi = carbon number of the ith component

α = α function in PR EOS β = adjustable exponent defined in eq 9 δij = BIP between component i and j γ = specific gravity, water = 1 γ+ = specific gravity of the plus fraction γP = specific gravity of normal-paraffin hydrocarbon, water = 1 Δγv = parameter in the Twu critical-volume correlation η = generic symbol for any property θ = adjustable exponent defined in eq 10

1304


Energy & Fuels

Article

ω = acentric factor

■

(20) Lohrenz, J.; Bray, B. G.; Clark, C. R. Calculating viscosities of reservoir fluids from their compositions. J. Pet. Technol. 1964, 16 (10), 1171−1176. (21) Pedersen, K. S.; Thomassen, P. Fredenslund, Aa. SRK-EOS calculation for crude oils. Fluid Phase Equilib. 1983, 14, 209−218. (22) Pedersen, K. S.; Thomassen, P. Fredenslund, Aa. Thermodynamics of petroleum mixtures containing heavy hydrocarbons. 1. Phase Envelope Calculations by Use of the Soave-Redlich-Kwong Equation of State. Ind. Eng. Chem. Process Des. Dev. 1984, 23 (1), 163−170. (23) Pedersen, K. S.; Thomassen, P.; Fredenslund, Aa. Properties of Oils and Natural Gases (Contributions in Petroleum Geology and Engineering) (Vol. 5). Gulf Publishing Co.: Houston, 1989. (24) Pedersen, K. S.; Christensen, P. L. Phase Behavior of Petroleum Reservoir Fluids; CRC Press: London, 2006, pp 88−89. (25) Krejbjerg. K.; Pedersen, K. S. Controlling VLLE equilibrium with a cubic EOS in heavy oil modeling. Proceedings of the Petroleum Society’s 7th Canadian International Petroleum Conference Canada (57th Annual Technical Meeting); Calgary, AB, Canada, June 13−15, 2006; Paper 2006−052. (26) Whitson, C. H. Characterizing hydrocarbon plus fractions. Proceedings of the European Offshore Petroleum Conference and Exhibition UK; London, UK, October 21−24, 1980; Paper SPE 12233. (27) Whitson, C. H. Effect of C7+ properties on equation-of-state predictions. Proceedings of the SPE Annual Technical Conference and Exhibition USA; New Orleans, LA, USA, September 26−29, 1982; Paper SPE 11200. (28) Whitson, C. H. Effect of C7+ Properties on equation-of-state predictions. SPE J. 1984, 24 (6), 685−696. (29) Danesh, A.; Xu, D.; Todd, A. C. A grouping method to optimize oil description for compositional simulation of gas injection processes. SPE Res. Eng. 1992, 7 (3), 343−348. (30) Rodrigurez, I.; Hamouda, A. A. An approach for characterization and lumping of plus fractions of heavy oil. Proceedings of the International Thermal Operations and Heavy Oil Symposium Canada; Calgary, AB, Canada, October 20−23, 2008; Paper SPE 117446. (31) Zuo, J. Y.; Zhang, D. Plus fraction characterization and PVT data regression for reservoir fluids near critical conditions. Proceedings of the SPE Asia Oil and Gas Conference and Exhibition Australia; Brisbane, QLD, Australia, October 16−18, 2000; Paper SPE 64520. (32) Diaz, O. C.; Modaresghazani, J.; Satyro, M. A.; Yarranton, H. W. Modeling the phase behavior of heavy oil and solvent mixtures. Fluid Phase Equilib. 2011, 304 (1−2), 74−85. (33) Jia, N.; Memon, A.; Gao, J.; Zuo, J.; Zhao, H.; Ng, H.-J. Threephase equilibrium study for heavy-oil/solvent/steam system at high temperatures. Proceedings of the Canadian Unconventional Resources and International Petroleum Conference Canada; Calgary, AB, Canada, October 19−21, 2010; Paper SPE 137453. (34) Chueh, P. L.; Prausnitz, J. M. Vapor-liquid equilibria at high pressures: Calculation of partial molar volumes in non polar liquid mixtures. AIChE J. 1967, 13 (6), 1099−1107. (35) Stryjek, R. Correlation and prediction of VLE data for n-alkane mixtures. Fluid Phase Equilib. 1990, 56, 141−152. (36) Gao, G. A simple correlation to evaluate binary interaction parameters of the Peng-Robinson equation of state: Binary Light Hydrocarbon Systems. Fluid Phase Equilib. 1992, 74 (15), 85−93. (37) Kordas, A.; Tsoutsouras, K.; Stamataki, S.; Tassios, D. A generalized correlation for the interaction coefficients of CO2hydrocarbon binary mixtures. Fluid Phase Equilib. 1994, 93, 141−166. (38) Peneloux, A.; Rauzy, E.; Freze, R. A consistent correction for Redlich-Kwong-Soave volumes. Fluid Phase Equilib. 1982, 8 (1), 7−23. (39) Jhaverl, B. S.; Youngren, G. K. Three-parameter modification of the Peng-Robinson equation of state to improve volumetric predictions. SPE Res. Eng. 1988, 3 (3), 1033−1040. (40) Twu, C. H.; Chan, H. S. Rigorously universal methodology of volume translation for cubic equation of state. Ind. Eng. Chem. Res. 2009, 48 (12), 5901−5906. (41) Pedersen, K. S.; Milter, J.; Sorensen, H. Cubic equations of state applied to HT/HP and highly aromatic fluids. SPE J. 2004, 9 (2), 186−192.

REFERENCES

(1) Zheng, S.; Yang, D. Pressure maintenance and improving oil recovery in terms of water-alternating-CO2 processes in thin heavy oil reservoirs. SPE Reservoir Evaluation & Engineering, Paper SPE 157719PA http://dx.doi.org/10.2118/157719-PA. (2) Fairfield, W. H.; White, P. D. Lloydminster fireflood performance, modifications promise good recoveries. Oil Gas J. 1982, 80 (6), 101−102. (3) Dyer, S. B.; Huang, S. S.; Farouq Ali, S. M.; Jha, K. N. Phase behaviour and scaled model studies of prototype saskatchewan heavy oils with carbon dioxide. J. Can. Pet. Technol. 1994, 33 (8), 42−48. (4) Upreti, S. R.; Lohi, A.; Kapadia, R. A.; El-haj, R. Vapor extraction of heavy oil and bitumen: A review. Energy Fuels 2007, 21 (3), 1562− 1574. (5) Ivory, J.; Chang, J.; Coates, R.; Forshner, K. Investigation of cyclic solvent injection process for heavy oil recovery. J. Can. Pet. Technol. 2010, 49 (9), 22−23. (6) Zhao, L. Steam alternating solvent process. SPE Res. Eval. Eng. 2007, 10 (2), 185−190. (7) Firouz, A. Q.; Torabi, F. Feasibility study of solvent-based huff-npuff method (cyclic solvent injection) to enhance heavy oil recovery. Proceedings of the SPE Heavy Oil Conference Canada; Calgary, AB, Canada, June 12−14, 2012; Paper SPE 157853. (8) Li, H.; Zheng, S.; Yang, D. Enhanced swelling effect and viscosity reduction of solvents-CO2-heavy oil systems. Proceedings of the SPE Heavy Oil Conference and Exhibition Kuwait; Kuwait City, Kuwait, December 12−14, 2011; Paper SPE 150168; SPE J., in press. (9) Li, H.; Yang, D. Phase behaviour of C3H8-n-C4H10-heavy oil systems at high pressures and elevated temperatures. J. Can. Pet. Technol., Paper SPE 157744-PA http://dx.doi.org/10.2118/157744PA. (10) Li, H.; Yang, D.; Li, X. Determination of three-phase boundaries of solvent(s)-CO2-heavy oil systems under reservoir conditions. Energy Fuels, 2013, 27 (1), 145−153. (11) Sayegh, S. G.; Maini, B. B. Laboratory evaluation of the CO2 huff-n-puff process for heavy oil reservoirs. J. Can. Pet. Technol. 1984, 23 (3), 29−36. (12) Li, H.; Yang, D.; Tontiwachwuthikul, P. Experimental and theoretical determination of equilibrium interfacial tension for the solvent(s)-CO2-heavy oil systems. Energy Fuels 2012, 26 (3), 1776− 1786. (13) Saber, N.; Shaw, J. M. On the phase behaviour of athabasca vacuum residue + n-decane. Fluid Phase Equilib. 2011, 302 (1−2), 254−259. (14) Kariznovi, M.; Nourozieh, H.; Abedi, J. Bitumen characterization and pseudocomponents determination for equation of state modeling. Energy Fuels 2010, 24 (1), 624−633. (15) Whitson, C. H. Critical properties estimation from an equation of state. Proceedings of the SPE/DOE 4th Symposium on Enhance Oil Recovery USA; Tulsa, OK, USA, April 15−18, 1984; Paper SPE 12634. (16) Behrens, R. A.; Sandler, S. I. The use of semicontinuous description to model the C7+ fraction in equation of state calculations. Proceedings of the SPE/DOE 5th Symposium on Enhance Oil Recovery USA; Tulsa, OK, USA, April 20−23, 1986; Paper SPE 14925. (17) Ahmed, T.; Cady, G. V.; Story, A. L. A generalized correlation for characterizing the hydrocarbon heavy fraction. Proceedings of the 60th SPE Annual Technical Conference and Exhibition USA; Las Vegas, NV, USA, September 22−25, 1985; Paper SPE 14266. (18) Li, Y.; Nghime, L. X.; Siu, A. Phase behavior computations for reservoir fluids: Effect of pseudo-components on phase diagrams and simulation results. Proceedings of the Petroleum Society of CIM Annual Meeting Canada; Calgary, AB, Canada, June 10−13, 1984; Paper CIM 84−35−19. (19) Naji, H. S. Feasible C7+ splitting methods an object-oriented approach. Int. J. Eng. Technol. 2010, 10 (1), 43−60. 1305


Energy & Fuels

Article

(42) Kesler, M. G.; Lee, B. I. Improve predictions of enthalpy of fractions. Hydro. Proc. 1976, 55, 153. (43) Roess, L. C. Determination of critical temperature and pressure of petroleum fractions. J. Inst. Pet. Tech. 1936, 22, 1270. (44) Cavett, R. H. Physical data for distillation calculations-vaporliquid equilibria. Proceeding of 27th American Petroleum Institute (API) Meeting USA; San Francisco, California, USA, 1962, 351−366. (45) Riazi, M. R.; Daubert, T. E. Simplify property prediction. Hydro. Proc. 1980, 115−116. (46) Watson, K. M.; Nelson, E. F. Improved methods for approximating critical and thermal properties of petroleum. Ind. Eng. Chem. 1933, 25 (8), 880−887. (47) Whitson, C. H.; Brule, M. R. Phase Behavior; Monograph Series, SPE: Richardson, TX, 2000. (48) Soreide, I. Improved phase behavior predictions of petroleum reservoir fluids from a cubic equation of state. Ph.D. Dissertation, Norwegian Institute of Technology (NTH), Trondheim, Morway, 1989. (49) Lee, B. I.; Kesler, M. G. A generalized thermodynamic correlation based on three-parameter corresponding states. AIChE J. 1975, 21 (3), 510−527. (50) Mehrotra, A. K.; Svrcek, W. Y.; Fu, C. T.; Puttagunta, R. Experiments to reconcile the data differences in the bitumen gas solubilities; UC/ARC Joint Report to AOSTRA; University of Calgary: Calgary, AB, 1985. (51) Twu, C. H. An internally consistent correlation for predicting the critical properties and molecular weights of petroleum and coal-tar liquids. Fluid Phase Equilib. 1984, 16 (2), 137−150. (52) Spencer, C. F.; Danner, R. P. Improved equation for prediction of saturated liquid density. J. Chem. Eng. Data 1972, 17 (2), 236−241. (53) Hong, K. C. Lumped-component characterization of crude oils for compositional simulation. Proceedings of the SPE/DOE Enhanced Oil Recovery Symposium USA; Tulsa, OK, USA, April 4−7, 1982; Paper SPE 10691. (54) Peng, D. Y.; Robinson, D. B. A new two-constant equation of state. Ind. Eng. Chem. Fundam. 1976, 15 (1), 59−64. (55) Li, H.; Yang, D. Modified α function for the Peng-Robinson equation of state to improve the vapor pressure prediction of nonhydrocarbon and hydrocarbon compounds. Energy Fuels 2011, 25 (1), 215−223. (56) Edmister, W. C.; Lee, B. I. Applied Hydrocarbon Thermodynamics (Vol. I, 2nd Ed.) Gulf Publishing Co.: Houston, 1984. (57) Graboski, M. S.; Daubert, T. E. A modified Soave equation of state for phase equilibrium calulations: 2. Systems Containing CO2, H2S, N2, and CO. Ind. Eng. Chem. Process Des. Dev. 1978, 17 (4), 448− 454. (58) Kato, K.; Nagahama, K.; Hirata, M. Generalized interaction parameters for the Peng-Robinson equation of state: Carbon dioxide n-paraffin binary systems. Fluid Phase Equilib. 1981, 7 (3−4), 219− 231. (59) Moysan, J. M.; Paradowski, H.; Vidal, J. Prediction of phase behaviour of gas-containing systems with cubic equations of state. Chem. Eng. Sci. 1986, 41 (8), 2069−2074. (60) Oellrich, L.; Plocker, U.; Prausnitz, J. M.; Knapp, H. Equationof-state methods for computing phase equilibria and enthalpies. Int. Chem. Eng. 1981, 21 (1), 1−15. (61) Computer Modelling Group Ltd. WINPROP Phase Property Program; Computer Modelling Group Ltd.: Calgary, AB, Canada, 2011. (62) Teja, A. S.; Sandler, A. I. A corresponding states equation for saturated liquid densities. II. Application to the calculation of swelling factors of CO2-crude oil systems. AIChE J. 1980, 26 (3), 341−345. (63) Michelsen, M. L. The isothermal flash problem: Part I. Stability. Fluid Phase Equilib. 1982, 9 (1), 1−19. (64) Firoozabadi, A. Thermodynamics of Hydrocarbon Reservoirs; McGraw-Hill: New York, NY, 1999. (65) Ali, A. A. New strategic method to tune equation-of-state to match experimental data for compositional simulation. Ph.D. Dissertation, Texas A&M University, College Station, TX, 2004.

(66) Trebble, M. A.; Salim, P. H.; Sigmund, P. M. A generalized approach to prediction of two and three phase CO2-hydrocarbon equilibria. Fluid Phase Equilib. 1993, 82, 111−118.

1306


Determination of Multiphase Boundaries and ... - ACS Publications

Recommend Documents