ethane mixtures

Fuel 252 (2019) 439–447

Contents lists available at ScienceDirect

Fuel journal homepage: www.elsevier.com/locate/fuel

Full Length Article

Application of Peng-Robinson equation of state for modelling the multiphase equilibrium properties in Athabasca bitumen/ethane mixtures

T

⁎

Ting Suna, Ali Takbiri-Borujenib, , Hossein Nouroziehc, Ming Gub a

China University of Petroleum, Beijing, China West Virginia University, United States c Computer Modelling Group Ltd, Canada b

A R T I C LE I N FO

A B S T R A C T

Keywords: Phase behavior Bitumen Peng-Robinson equation of state Solvent

Phase behavior of (bitumen + solvent) mixtures and their properties are essential for surface upgrading methods and enhanced oil recovery applications. The objective for this work is to systematically model the phase partitioning and component distribution in the (bitumen + solvent) mixtures. To perform equation of state modelling, compositional analysis of bitumen is utilized to characterize the bitumen sample and its components up to C100. The full range of components is utilized in Peng-Robinson equation of state (PR-EOS) to predict the equilibrium phase properties, such as composition and density, in vapor-liquid and liquid-liquid regions. Based on the results, PR-EOS accurately captures the composition of phases, component distribution in each phase, and phase densities. Thereafter, an algorithm is proposed to lump the components on the basis of two-phase envelope generation. The results indicate that Athabasca bitumen and its phase behavior in (bitumen + solvent) mixtures can be well represented in multi-phase regions by four pseudo-components.

1. Introduction Determining phase behavior and estimating the thermodynamic properties for complex fluid systems or ill-defined fluids, such as (bitumen + solvent) mixtures, present significant challenges. Despite the importance of phase behavior data for production, pipeline transportation, upgrading, and refining of such fluids, limited experimental data for bitumen-containing mixtures have been reported in the literature. This is due to the fact that experimental measurements for these systems are usually expensive and time-consuming. This issue is even more pronounced when multiple phases form at equilibrium conditions. Modelling is a credible alternative for studying the phase behavior and the interaction of (bitumen + solvent) systems. The thermodynamic modelling of (bitumen + solvent) systems is limited compared to simple hydrocarbon systems and conventional oil. This might be due to complex nature and lack of a known characterization method for heavy oil and bitumen. Although different characterization schemes, as stated by Kumar and Okuno [1], have been proposed for different reservoir fluids [2–6]; for example, gas condensates [7–10], volatile oils [11,12], near-critical fluids [13–16], the methods for heavy oil and bitumen characterizations are limited [1,17–18]. This led to variety of models in literature for pseudo-component determination in the phase behavior study of

⁎

(bitumen + solvent) mixtures. Mehrotra et al. [19] utilized Bishnoi et al. [20] characterization method to develop five pseudo-component in order to study the phase behavior and solubility of carbon dioxide and ethane in Athabasca and Peace River bitumens using Peng-Robinson equation of state (PR-EOS). The critical properties calculated from Kesler and Lee [21] correlation provided a good agreement with the experimental solubility data. Fu and Puttagunta [22] and Fu et al. [23] considered a pseudo-component in the vapour-liquid equilibrium calculations of (bitumen + solvent) systems for the sake of simplicity. They developed a temperature-dependent correlation for binary interaction parameters to represent the vapour-liquid data. Mehrotra and Svrcek [24,25] and Mehrotra et al. [26] utilized PR-EOS along with Kesler and Lee [21] correlation for calculating critical properties to estimate the nitrogen, carbon monoxide, methane, carbon dioxide, and ethane solubility in Athabasca, Cold Lake, and Wabasca bitumen at the temperatures up to 110 °C and at the pressures up to 10 MPa. The bitumen was characterized using three pseudo-components which represent distillable maltene, undistillable maltene, and asphaltene that could adequately match the experimental data. Frauenfeld et al. [27] estimated the methane, ethane, propane, and carbon dioxide solubility in Lloydminster Aberfeldy oil and a blended Cold Lake/Lloydminster oil using PR-EOS. The binary interaction parameters were adjusted to tune the model.

Corresponding author. E-mail address: [email protected] (A. Takbiri-Borujeni).

https://doi.org/10.1016/j.fuel.2019.04.106 Received 28 February 2019; Received in revised form 11 April 2019; Accepted 16 April 2019 Available online 01 May 2019 0016-2361/ © 2019 Elsevier Ltd. All rights reserved.

Fuel 252 (2019) 439–447

T. Sun, et al.

was then tuned to match solubility data of propane in each bitumen cut. The bitumen was also characterized using the boiling point or carbon number distribution and the generalized model was then proposed to obtain the solubility of propane in the bitumen samples. Most of the above-mentioned studies have focused on the modelling of the vapor-liquid equilibrium and solubility data. The phase behavior of (bitumen + solvent) mixtures is still challenging when the multiphase equilibrium occurs at specific conditions. The objective of this study is thus to develop a characterization method for heavy oil and bitumen. The method is applicable for both the multiple-phase equilibrium conditions, and the development of PVT models for the purpose of simulation study. In addition, it tests the applicability of equation of state for phase behavior modelling study of (bitumen + solvent) mixtures. The characterization is based on the simulated distillated data and physical properties of bitumen coupled with available correlations for the calculation of critical properties. WinProp software from Computer Modelling Group (CMG) version 2017 was used to model the phase behaviour of the (bitumen + solvent) mixtures. The phase behavior predictions were examined with the available literature data for bitumen/ethane mixtures including vapor-liquid (VLE) and liquid-liquid equilibrium (LLE) data. The compositional analysis, phase compositions, and phase densities are compared against the experimental data and the phase diagrams are generated. These information are extremely important for surface upgrading, production, pipeline transportation, and refining.

Kariznovi et al. [28] modelled the phase behaviour of (bitumen + solvent) mixtures using PR-EOS for four bitumens (Athabasca, Peace River, Wabasca, and Cold Lake) mixed with four different solvents (nitrogen, carbon dioxide, methane, and ethane). They used the lab data of Mehrotra and Svrcek [24,25,29] and Svrcek and Mehrotra [30,31] and considered three, five and six pseudo-components which were defined on the basis of different fractions of bitumen. A universal set of pseudo-components was also proposed applicable for different Alberta bitumens. Saber and Shaw [32] and Saber et al. [33] used a group contribution method to estimate the critical properties of bitumen pseudo-components for equation of state phase behavior modelling of (Athabasca bitumen + n-decane) mixtures. Diaz et al. [34] studied the phase behaviour for the mixtures of Athabasca bitumen, propane, and carbon dioxide using different characterization methods. They utilized simulated distillation (SimDis) data along with six characterization approaches to predict the saturation pressures and liquid-liquid phase boundaries of pseudo-ternary systems. Agrawal et al. [35] modelled the saturation pressure of live Peace River bitumen at different condensate concentrations with PR-EOS using the characterization method proposed by Diaz et al. [34] through the adjustment of binary interaction parameters. They fitted the saturation pressure and asphaltene onset data with temperature-dependent binary interaction parameters. The developed model reasonably predicted the saturation pressures of methane/bitumen, ethane/bitumen, carbon dioxide/bitumen, and n-pentane/bitumen. Li and Yang [36] modelled the phase behaviour of propane/n-butane/heavy-oil systems at the temperatures up to 123 °C and the pressures up to 5030.0 kPa with a volume-translated PR-EOS coupled with a modified alpha function. The experimental data was accurately matched with the binary interaction correlations and the volume-translated equation of state, while prediction accuracy was decreased near the pure solvent critical temperature. Badamchi-Zadeh et al. [37] proposed a bitumen characterization method with an equation of state model for the carbon dioxide, propane, and Athabasca bitumen mixtures. Four pseudo-components were defined for bitumen based on conversion of the SimDist data of asphaltene free bitumen to a normal boiling point distribution. Then, the binary interaction parameters were tuned to match the saturation pressure. Yang et al. [38] modelled the phase behavior of methane/ propane/heavy-oil (Lloydminster) mixture using the volume-translated PR-EOS along with a modified alpha function at the pressures up to 10 MPa and the temperature range of 22–50 °C. A new correlation for the binary interaction parameter was proposed in order to match the measured swelling factors and saturation pressures. Ghasemi and Whitson [39] proposed a systematic approach to model the phase behaviour of the mixtures containing light solvents (nitrogen, carbon dioxide, carbon monoxide, methane, and ethane) and Athabasca bitumen using cubic equations of state. Gamma distribution was used to fit the bitumen distillation data; then, single-carbon-number (SCN) fractions were defined, and Twu [40] correlation was used to calculate the critical properties of fractious. Finally, SCN fractions were lumped into five pseudo-fractions. Gao et al. [41] observed the multiphase behavior for (bitumen + nbutane) mixtures in their phase behavior studies at pressures up to 10 MPa and temperatures up to 160C. Eghbali and Dehghanpour [42] presented experimental and modelling study for multiphase behavior of (Clearwater bitumen + butane) and (Clearwater bitumen + CO2) mixtures. They calibrated PR-EOS to predict (bitumen + solvent) phase behavior. Kumar and Okuno [43] proposed a method for characterization of bitumen by improving the fluid characterization method which was developed based on perturbation from n-alkanes (PnA). The authors showed that LLE data needed positive interaction coefficients between (bitumen + solvent) while VLE was not sensitive to the interaction coefficients used. Azinfar et al. [44] fractionated Athabasca bitumen to four cuts by using modified vacuum distillation. The PR-EOS

2. Equilibrium calculation State of the system is determined by fugacity of each component in each phase at equilibrium conditions. The equilibrium condition for species i in phase I and II is presented as,

fiI = fiII

(1)

In terms of the fugacity coefficients, Eq. (1) can be re-written as,

yi φi ̂ P = x i φiIÎ P, I

(2)

where the partial fugacity coefficients are calculated using, ∞ − ∫v ⎛ ⎡ ∂∂nPi ⎤ ⎝ ⎣ ⎦T , v, nj ≠ i ⎜

ln φi ̂ =

RT

RT ⎞ ⎟ dv v

⎠

− ln Z ,

(3)

where P is pressure, T is temperature, v is volume, R is gas universal constant, and Z is compressibility factor. Equation of state (EOS) is typically utilized to express the relationship between the pressure, temperature and volume. Peng-Robinson (Robinson and Peng [45]) EOS, a cubic EOS, provides acceptable results for the hydrocarbon mixtures and petroleum systems and thus, it was employed in this study to calculate the equilibrium properties of bitumen containing mixtures. The general form of Peng-Robinson EOS is shown as,

P=

RT aλ (T ) , − v−b v (v + b) + b (v − b)

(4)

where,

a=

0.45724R2Tc2 Pc

b=

0.07780RTc , Pc

λ (T ) = [1 + κ (1 −

(5)

(6)

Tr1/2)]2 ,

κ = 0.379642 + 1.48503ω − 0.164423ω2 + 0.016666ω3

(7) (8)

where P is pressure, T is temperature, v is volume, a represents the attractive forces between molecules, b accounts for the volume of molecules, Tc, and Pc are the critical temperature and pressure. a and b are 440

Fuel 252 (2019) 439–447

T. Sun, et al.

calculated from the critical properties of the component for a single component system. However, for multicomponent mixtures, a and b are defined based on the individual components’ a and b using the van der Waals mixing rules:

am =

∑ xi Di ,

(9)

i

Di = bm =

ai

∑j xj (1 − δij )

∑ xi bi , i

Table 1 Bitumen properties [58].

aj , and

(10)

δij = 1 −

−1/6

μmix (T , P ) Tc,mix ⎞ = ⎜⎛ ⎟ μref (Tref , Pref ) ⎝ Tc,ref ⎠

1/2

⎛ MWmix ⎞ ⎝ MWref ⎠

⎜

⎟

⎛⎜ βmix ⎞⎟, ⎝ βref ⎠

(18)

(19)

The rotational coefficient (β) represents the deviation from simple corresponding states, and it is defined as,

(13)

d4 β = 1 + d3 (ρ / ρc )ref MW d5

(20)

where ρc, and ρ are the critical density and mixture density. The viscosity of a mixture calculated using the Pedersen model is dependent on the critical pressures, critical temperatures, molecular weights of each component, and the coefficients di (as shown in the above equations). A modified version of Pedersen correlation for methane viscosity gives a better correlation of the viscosity data [54].

(14)

where xi is the mole fraction of component i in the mixture and ci is defined as (15)

where Ωb is a constant specific to the EOS (for Peng-Robinson EOS, it is 0.077796), vs is the volume translation, Tc is the critical temperature, and Pc is the critical pressure. vs is also considered as a tuning parameter that commonly obtained through matching the lab volumetric data (e.g. densities, relative volumes, and etc.). Thus, the volume translation technique is available to modify the calculated volumetric and also to match the bitumen and equilibrium phases’ densities. In WinProp, three methods are available for the volume shift: constant volume shift, temperature dependent volume shift correlation, and linear temperature dependent volume shift. Using constant volume shift option, the densities of the equilibrium phases are inaccurately modeled with the EOS particularly when a wide range of temperatures is considered. Thus, it is essential to define temperature-dependent volume shifts. Pedersen et al. [51] proposed the following equation to calculate the temperature-dependent volume shifts for each component,

S = S0 + S1 (T − Tref )

2/3

⎛ Pc,mix ⎟⎞ ⎝ Pc,ref ⎠

⎜

MWmix = d1 (MW wd2 − MWnd2) + MWn

N

ci = Ωb vsi RTc / Pc

(17)

where Tc, and Pc are the critical temperature and pressure. The molecular weight (MW) is calculated using both weight and molar averaged values as,

where v is corrected molar volume and c in translation parameter. For multi-components systems, c is defined as,

i=1

]

In WinProp, Pedersen viscosity correlation [53] is used to calculate the viscosity of high viscous fluids such as heavy oil and bitumen. This correlation is based on the corresponding state method, and it calculates the viscosity of mixtures with respect to a reference substance at the same reduced temperature and pressure using the following relation,

cor

∑ x i ci ,

2/7

4. Viscosity calculation

EOS can also be utilized to determine the volumetric properties of petroleum/hydrocarbon mixtures. Peneloux et al. [50] proposed a volume translation technique to be applied to Peng-Robinson EOS in order to improve the volumetric calculations without impacting the other properties. In the volume shift technique, the volume is corrected by the following relation,

c=

11.76 57.00 21.61 9.62

Rackett’s compressibility factors are back-calculated from the specific gravity and critical properties with the assumption that the specific gravity is roughly equal to the saturated liquid density at 60°F.

(12)

3. Density calculation

= v − c,

SARA Fraction (wt%) Saturates Aromatics Resins Asphaltenes

[1 + (1 − Tr ) vs = (RTc / Pc ) ZRA

Oellrich et al. [48] showed that a value of 1.2 for n provides an acceptable match with the experimental results for the paraffin-paraffin interaction coefficients. This equation satisfactorily reproduces the methane interaction coefficients reported by Katz and Firoozabadi [49] (Mehra [46]). In this study, n was considered as a tuning parameter to fit the experimental VLE and LLE data. Worth to mention that, Diaz et al. [34] defined the binary interaction parameters as a function of components’ critical temperatures to match the saturation pressure data of (Athabasca bitumen + solvent) mixtures.

v cor

512.5 ± 6.9

Component volume shifts can be evaluated by finding the difference between the saturated liquid molar volume of a pure component obtained from the EOS and the saturated liquid molar volume calculated using a modified Rackett‘s Z-factor equation at any specific temperature [52],

(11)

where δij represents the binary interaction between the hydrocarbons in the mixture, which is typically defined as a function of the critical properties. In WinProp, the binary interaction parameters are related to critical volume of components through the following equation [46,47], 1/6 1/6 n ⎛ 2vci vcj ⎞ ⎜ v1/3 + v1/3 ⎟ cj ⎠ ⎝ ci

Molecular Weight (g/mol)

5. Bitumen characterization Before applying an EOS for a phase behavior study, the fluid samples should be thoroughly characterized. The characterization of heavy fractions such as heavy oil and bitumen is generally on the basis of the boiling point distribution from atmospheric and vacuum distillation. Crude oils and distillate fractions contains thousands of components with complex molecular structures that typically show a continuous boiling point distribution. Two commonly used approaches for oil characterization are the characterization of the crude oil or its fractions: (i) in terms of pseudo-components that are defined by hydrocarbon type, e.g., paraffinic, olefinic, naphthenic or aromatics, or an average boiling point (ii) the continuous mixture approach, in which the oil is represented by a continuous distribution function with respect to boiling point, molecular mass, and/or some other measurable properties. The continuous mixture approach is adequate when the oil consists of one type of hydrocarbon groups, such as paraffins or naphthenic;

(16) 441

Fuel 252 (2019) 439–447

T. Sun, et al.

3 Exp

Composition (wt%)

Model

2

1

0 5

15

25

35

45

55

65

75

85

95

Single Carbon Number Fig. 1. Component distribution in raw Athabasca bitumen; experimental data from Kariznovi et al. [57].

1000

Table 2 Pseudo-components properties for PR-EOS; MW, molecular weight; Pc, critical pressure; Tc, critical temperature; SG, specific gravity; w, acentric factor.

Boiling Point (°C)

800

600

MW (g/g-mol)

Pc (MPa)

Tc (°C)

SG

ω

C6-19 C20-34 C35-99 C100+

210.2 360.4 806.2 2034.0

2.234 1.482 0.833 0.542

461.9 587.4 746.4 959.5

0.892 0.952 1.021 1.153

0.575 0.877 1.289 1.590

400

however, heavy oil and bitumen are a complex mixture composed of different hydrocarbons types (e.g. saturates, aromatics, etc.). In this study, the pseudo-component approach is employed for the bitumen characterization. The pseudo-components are primarily defined, and their critical properties are estimated and along with the component’s composition are applied into the Peng-Robinson EOS. The multi-phase equilibrium is then calculated by a flash method combined with a phase stability check. In this study, WinProp software from Computer Modelling Group (CMG) was used to model the phase behaviour of the bitumen containing mixtures. The tangent-plane stability criterion is used in WinProp [55]. This method is on the basis of a stagewise procedure in which the number of phases gradually increases. The equilibrium equation (Eq. (3)) is solved with Quasi-Newton Successive Substitution (QNSS) [56]. The compositional analysis of Athabasca bitumen sample was already reported by Kariznovi et al. [57], and was used to define the pseudo-components and characterize the bitumen. Table 1 summarizes the bitumen properties [58]. Based on the lab distillation data, components were defined for the distillable fraction of Athabasca bitumen, which accounts for about 80 percent (in mass) of the fluid, and also matches the component distribution. Riazi [59] developed the following correlation to relate the component’s boiling point to corresponding carbon number,

200

0

0

20

40

60

80

100

Cumulative Distilled (wt%) Fig. 2. Boiling point temperature with respect to weight percent distilled for raw Athabasca bitumen (Δ) taken from Kariznovi et al. [57].and proposed characterization scheme (line).

2500

2000

Pressure (kPa)

Component

1500

1000

500

Tb = 1090 − exp(6.9955 − 0.11193Nc2/3), 0

0

200

400

600

(21)

where Tb is the boiling point in Kelvin, and Nc is carbon number. Riazi [59] showed that the data of Katz-Firoozabadi [49] is represented by this equation with an average absolute deviation (AAD) of 0.2%. Riazi [59] also utilized the physical properties of hydrocarbons (C6–C22) to develop correlations for the estimation of molecular weight as a function of normal boiling point, and the prediction of specific gravity as a function of molecular weight. The molecular weight of each single carbon number, as proposed by Riazi [59], can be estimated by the following equation,

800

Temperature (°C) Fig. 3. Pressure-temperature diagram for full characterization (solid line) and lumped schemes (dash line).

442

Fuel 252 (2019) 439–447

T. Sun, et al.

10

Equilibrium Pressure (MPa)

8

Pressure (MPa)

Fig. 4. The pressure-composition (P-x) diagram (a) and coexisting equilibrium phase compositions (b) for the mixture of (Athabasca bitumen + ethane) at ambient temperature; ○, Mehrotra and Svrcek [29] data; ■, Nourozieh et al. [65] data; solid line, equation of state using full characterization scheme; dashed line, equation of state using lumped scheme.

10

Single L LLE

6

4 VLE

2

25

50

75

LLE

6 4 VLE

2

(b)

(a)

0

0 0

8

0

100

25

50

75

100

Equilibrium Compositions (wt%)

100

Fig. 5. Effect of overall solvent composition on the co-existing phase composition for mixtures of (Athabasca bitumen + ethane) at 5 MPa (a) and 9 MPa (b) and ambient temperature;▲, solvent-enriched phase; □, asphalteneenriched phase (data from Nourozieh et al. [65]); solid line, equation of state using full characterization scheme; dashed line, equation of state using lumped scheme.

100

75

50

25

(a)

0 20

40

60

80

Overall Solvent Composition (wt%)

MW = [[6.97996 − ln(1080 − Tb)]/0.01964]3/2

100

Co-excisting Phase Composition (wt%)

Co-excisting Phase Composition (wt%)

Solvent Overall Composition (wt%)

75

50

25

(b)

0 20

40 60 80 Overall Solvent Composition (wt%)

calculating the critical properties of Alberta bitumen and they concluded that the Kesler and Lee [21] correlations provided the most acceptable results for bitumen/gas systems. Kariznovi et al. [28] and Diaz et al. [34] also examined different correlations to calculate critical properties and acentric factor of non-distillable fractions. They came into the conclusion that the phase boundaries are better predicted with Lee-Kesler correlation compared to other available methods. Thus, in this study, the Kesler and Lee [21] correlation was used to calculate the critical properties and acentric factor of the components. Kesler and Lee [21] proposed the following equations,

(22)

where Tb is the boiling point. The specific gravity of each component is also related to molecular weight by,

SG = 1.07 − exp[3.56073 − 2.93886 × MW 0.1]

(23)

where MW is the molecular weight. Eqs. (21)–(23) were used to calculate the boiling points, molecular weights, and specific gravities of SCN groups. The boiling points of the components along with their compositions reasonably represent the experimental distillation data of Athabasca bitumen. The un-distilled fraction of Athabasca bitumen was represented by a single component and assigned an average molecular weight based on Diaz et al.’s [34] work. Alboudwarej et al. [60] proposed the following equation for asphaltene fractions to relate the specific gravity and molecular weight,

SG = 670 × MW 0.0629

100

Tc = ⎛341.7 + 811SG + [0.76392 + 0.21132SG] Tb ⎝ (0.25939 − 1.81239SG ) × 105 ⎞ + /1.8 Tb ⎠ ⎜

⎟

(25)

(24)

Pc = exp[8.3634 − 0.0566/ SG − (0.24244 + 2.2898/ SG + 0.11857/SG2)

Eq. (24) along with the average molecular weight of asphaltene was used to estimate the specific gravity of non-distillable fraction in this study. Thus, on the basis of the developed characterization scheme, the distillable fraction of bitumen is represented by single carbon numbers, and the non-distillable fraction by one pseudo-component (composition of C100+ is 18.76 mass percent). The comparison of experimental and developed characterization scheme in terms of the component distribution and boiling point curve of Athabasca bitumen is presented in Figs. 1 and 2. The critical properties and acentric factor of components are typically correlated through the components’ molecular weight, specific gravity, and boiling point. Three commonly used correlations are available in WinProp to calculate these properties: Kesler and Lee [21], Twu [40], and Riazi and Daubert [61]. Svrcek and Mehrotra [30,31] and Mehrotra et al. [19] evaluated a number of these correlations for

× 10−3 (1.8 × Tb)+(1.4685 + 3.648/ SG + 0.47227/SG2) × 10−7 (1.8 × Tb)2 − (0.42019 + 1.6977/SG2) × 10−10 (1.8 × Tb)3]/14.7 (26)

(ln P − 5.92714 + (15.2518 − cr

ω=

6.09648 Tbr

15.6875 Tbr

+ 1.28862 ln Tbr − 0.169347Tbr6

− 13.4721 ln Tbr +

0.43577Tbr6

0.8

ω = −7.904 + 0.1352Kw − 0.007465K w2 + 8.359Tbr (1.408 − 0.01063Kw ) + for Tbr > 0.8 Tbr 443

)

) for T

br

⩽ (27)

(28)

Fuel 252 (2019) 439–447

T. Sun, et al.

10 L1 Exp.

(a)

L1 EOS

Composition (wt%)

8

L2 Exp. L2 EOS

6

4

2

0 5

15

25

35

45

55

65

75

85

95

Single Carbon Number 10

L1 Exp.

(b)

L1 EOS L2 Exp.

Composition (wt%)

8

L2 EOS

6

4

2

0 5

15

25

35

45

55

65

75

85

95

Single Carbon Number Fig. 6. Component distribution in ethane-free liquid samples, extract (black and green) and residue (blue and red), from LLE study of (Athabasca bitumen + ethane) mixtures at 5 MPa (a) and 9 MPa (b) and a constant overall ethane composition of 40 wt%; blue and black, experimental data from Kariznovi et al. [57]; red and green, predictions by PR-EOS using full characterization scheme. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Kw =

(1.8Tb)1/3 SG

rule [63]. This scheme calculates the number of components using the following equation,

(29)

where P is pressure in atm, T is temperature in K, Pcr and Tbr are reduced critical pressure and reduced boiling point temperature, respectively. As discussed earlier, the components’ critical volumes are needed in Eq. (12) to calculate the binary interaction parameters between the components. In this study, Twu’s correlation [40] was employed to calculate the critical volumes of the components using their boiling point, molecular weight, and specific gravity.

Ng = Int [1 + 3.3log10 (N − n)]

(30)

where, n and N are the minimum carbon number and the maximum one in the sample, respectively. WinProp has two options for lumping into hypothetical components: default and user defined lumping scheme. The default lumping is according to components log (k) values and the k-values are calculated using Wilson’s equations. The Lee-Kesler [64] mixing rule are then used to calculate the critical properties of the components,

6. Component lumping

vc =

Following the characterization of Athabasca bitumen, lumping or grouping of the components into fewer numbers of components is required because large number of pseudo-components is not commonly practical for modelling study. The accuracy of predictions could be reduced during lumping step while the calculation time is significantly improved. The adequate number of pseudo-components mainly depends on the mixture/system under study. Whitson [62] proposed a lumping approach that is developed originally on the basis of Sturges’

Tc =

ω=

1 8

1 8vc

i

(31)

j

∑ ∑ xi xk (vci1/3 + vck1/3 ) i

∑ xi ωi i

444

∑ ∑ xi xk (vci1/3 + vck1/3 )3

j

Tci Tcj (32)

(33)

Fuel 252 (2019) 439–447

80

80 60 60 40 40 20 20

0 0

20

40

60

80

0 100

10000000 1000000

Bitumen Viscosity (mPa.s)

100

Cumulative Carbon Number Distribution (L2)

Cumulative Carbon Number Ditsribution (L1)

T. Sun, et al.

100000 10000 1000 100 10 1 0

Single Carbon Number

50

100

150

200

Temperature (°C)

Fig. 7. Cumulative carbon number distribution in gas-free liquid samples from LLE study of (Athabasca bitumen + ethane) mixtures at 5 MPa (Δ, ▲) and 9 MPa (■, □) and a constant overall ethane composition of 40 wt%; dots, Kariznovi et al. [57] data; lines, predictions by PR-EOS using full characterization scheme.

Fig. 9. Comparison of calculated and measured bitumen viscosities; ■, Nourozieh et al. [65] data; solid line, corresponding state using full characterization scheme; dashed line, corresponding state using lumped scheme. 10

Equilibrium Pressure (MPa)

1030

Bitumen Density (kg/m3)

1005

980

955

930

8 LLE

6

4 VLE

2

905 0 0

880 0

50

100

150

Temperature (°C)

Zc RTc (0.2905 − 0.085ω) RTc = vc vc

∑ ∑ xi xk (vci1/3 + vck1/3 ) i

400

600

800

1000

Fig. 10. The coexisting equilibrium phase densities for the mixture of (Athabasca bitumen + ethane) at ambient temperature; ○, Mehrotra and Svrcek [29] data; ■, Nourozieh et al. [65] data; solid line, equation of state using full characterization scheme; dashed line, equation of state using lumped scheme.

Fig. 8. Comparison of calculated and measured bitumen densities; ▲, Nourozieh et al. [65] data; solid line, equation of state using full characterization scheme; dashed line, equation of state using lumped scheme.

Pc =

200

Density of Equilibrium Phases (kg/m3)

200

Tci Tcj

7. Results and discussions

j

(34)

As previously mentioned, the compositional analysis of Athabasca bitumen sample was used to characterize the oil. The compositions and properties of components were applied in PR-EOS. In order to fit the experimental data, the measured composition data of VLE and LLE was converted into k-values, and the exponent of binary interaction parameter equation was adjusted to match the k-values. The best fitted exponents (n in Eq. (12)) for un-lumped and lumped schemes are 0.4964 and 0.5034, respectively. Using these values, the pressurecomposition (p − x) diagrams for mixtures of (Athabasca bitumen + ethane) were generated at the temperature of 22 °C as shown in Fig. 4. The VLE data are from Mehrotra and Svrcek [29] study and the LLE data from Nourozieh et al. [65] study. As depicted in the figure, the LLE was observed at overall ethane compositions greater than 18 wt % and at pressures larger than 4 MPa. Two schemes (lumped and unlumped) accurately predict the VLE and LLE regions using the equation of state. The lumped scheme predicts the LLE separation in slightly lower overall ethane composition compared to un-lumped scheme.

where vc is the critical volume, Tc is the critical temperature, w is acentric factor, and Pc is the critical pressure. While the modelling results with full number of components might be acceptable, the number of pseudo-components can be reduced for computational efficiency. In this study, the two-phase envelope (pressure-temperature diagram) is generated with PR-EOS using full number of components. The defined pseudo-components are then lumped into fewer number of pseudocomponents using user defined method while the two-phase envelope is generated and then compared with full characterization scheme. The minimum number of pseudo-components resulting in the same twophase envelope as full characterization scheme is the optimum number of components. As depicted in Fig. 3, the Athabasca bitumen sample can be represented by four pseudo-components. The lumped pseudocomponents properties are listed in Table 2.

445

Fuel 252 (2019) 439–447

1000

Co-existing Phase Density (kg/m3)

Co-existing Phase Density (kg/m3)

T. Sun, et al.

800 600 400 200

(a) 0 20

40

60

80

Fig. 11. Effect of overall solvent composition on the co-existing phase density for (Athabasca bitumen + ethane) mixtures at 5 MPa (a) and 9 MPa (b) and ambient temperature; ▲, solventenriched phase; ■, asphaltene-enriched phase (experimental data from Nourozieh et al. [65]); solid line, equation of state using full characterization scheme; dashed line, equation of state using lumped scheme.

1000 800 600 400 200

(b) 0

100

20

40

AAD (kg/m3)

b

MAD (kg/m3)

Model

Data

AARD (%)

Un-lumped

Raw Bitumen Equilibrium

0.34 0.43

3.3 3.8

6.2 11.7

Lumped

Raw Bitumen Equilibrium

0.28 0.45

2.7 4.0

5.8 8.9

a b c

c

AARD: average absolute relative deviation. AAD: average absolute deviation. MAD: maximum absolute deviation.

Table 4 The deviations of models for the calculation of mixture viscosities. Model

Data

AARD (%)a

Un-lumped

Bitumen Equilibrium Data

14.6 21.3

Lumped

Bitumen Equilibrium Data

23.5 23.3

a

80

100

equation of state. Further investigation of results (Fig. 7) also shows that the non-distilled fractions of extract and residue are in close agreement with the lab values. The increase in overall ethane concentration or pressure leads in higher extraction of light components from the bitumen, and this results into an increase in the fraction of heavier components in the bitumen-enriched phase (L2) [65]. The proposed models (un-lumped and lumped schemes) have been examined for the solubility data over the temperature and pressure range of 20–200 °C and the pressure range of 0.1–10 MPa and they provide AARDs of 16.3% and 16.7%, respectively. The density and viscosity of bitumen is of importance for recovery processes and their accurate predictions provide better field design and optimization. Although the density and viscosity of bitumen are usually fitted and then predicted with available correlations, the gas-saturated and solvent-diluted bitumen properties are not accurately predicted with the correlations over the wide ranges of pressures, temperatures, and solvent compositions. This is even worse when the composition of original fluid is changing with temperature, pressure, and/or solvent composition due to precipitation or formation of a second liquid phase. Thus, it is desired to model the density and viscosity of fluid or equilibrium phases while the compositional calculation is done. In WinProp, the density and viscosity of equilibrium phases were simultaneously reported with compositional calculations. Using the temperature-dependent volume shift correlation, the density of raw bitumen is well predicted over the temperature range as shown in Fig. 8 using unlumped and lumped schemes. This indicates that the developed characterization and the defined components reasonably represent the bitumen. The viscosity of bitumen was also correlated with modified Pedersen’s correlation and the results are presented in Fig. 9. The coefficients of viscosity correlation are d1 = 0.000156, d2 = 2.309757, d3 = 0.005902, d4 = 2.216300, and d5 = 0.442047. The data is accurately modeled with this correlation over the studied temperature range. The density of equilibrium phases with dissolved solvent is also reasonably represented by equation of state in the VLE and LLE region as shown in Fig. 10. The effect of solvent overall composition on the equilibrium densities were evaluated and compared with experimental data. The results at two pressures 5 and 9 MPa are plotted in Fig. 11. The density of light liquid phase in LLE region is slightly overestimated at higher pressure (9 MPa). This indicates that the pressure dependency of density might be not captured well in the model. Tables 3 and 4 summarize the deviation of two models from experimental data.

Table 3 The deviations of models for the calculation of mixture densities. a

60

Overall Solvent Composition (wt%)

Overall Solvent Composition (wt%)

AARD: average absolute relative deviation.

Fig. 4 also demonstrates the solvent composition in the equilibrium phases at different pressures and at an overall ethane concentration of 40 wt%. As anticipated from the figure, the equilibrium compositions are reasonably predicted with two lumped and un-lumped schemes. Although the predictions show slightly deviation from experimental data at higher pressures, there is a good agreement between the measured data and predicted ones. As it was explained in the study by Nourozieh et al. [65], the solvent overall composition affects the solvent composition in equilibrium phases due to the nature of extraction process as well as the solvent affinity to extract light fractions/components. Thus, the equilibrium composition of solvent in two phases is plotted as a function of solvent overall composition at two different pressures in Fig. 5. The satisfactory predictions and similar results to experimental data were obtained for equilibrium compositions using two schemes. It is also worth to compare the component distribution in two liquid phases excluding the solvent. This indicates how the solvent extracts the light hydrocarbon components into the lighter liquid phase, and partitions the heavy constitutes into the heavy liquid phase. Fig. 6 shows the distribution of components in two flashed-off phases taken at the equilibrium pressure of 9 MPa and a constant overall ethane composition of 40 wt%. The results indicate that the Peng-Robinson EOS coupled with the developed characterization scheme provides a good agreement not only in the solvent composition in different phases but also in the component distribution in two phases. The same plot but at a constant pressure of 5 MPa is also shown in Fig. 6. Again, the component distribution in extract and residue are well represented by the

8. Conclusions The VLE and LLE data of (Athabasca bitumen + ethane) mixtures as well as the component distribution within each phase are well represented by PR-EOS and the proposed characterization scheme. The minimum number of components to represent the multi-phase 446

Fuel 252 (2019) 439–447

T. Sun, et al.

equilibrium was determined on the basis of two-phase envelope comparison. For the system under study, four components were found to be enough to accurately capture the equilibrium phase compositions and densities. The calibrated EOS with four components predicted multiphase equilibrium compositions, densities and viscosities within acceptable errors. The proposed models (un-lumped and lumped schemes) have been examined for the solubility data over the temperature and pressure range of 20–200 °C and the pressure range of 0.1–10 MPa and they provided AARDs of 16.3% and 16.7%, respectively. The density data were predicted with maximum AARD of 0.43%.

[30] Svrcek WY, Mehrotra AK. Gas solubility, viscosity and density measurements for Athabasca bitumen. J Can Petrol Technol 1982;21(4):31–8. 387. [31] Svrcek WY, Mehrotra AK. Properties of peace river bitumen saturated with field gas mixtures. J Can Petrol Technolol 1989;28(2):50–6. [32] Saber N, Shaw JM. On the phase behaviour of Athabasca vacuum residue+ n-decane. Fluid Phase Equilibr 2011;302(1–2):254–9. [33] Saber N, Zhang X, Zou XY, Shaw JM. Simulation of the phase behaviour of Athabasca vacuum residue+ n-alkane mixtures. Fluid Phase Equilibr 2012;313:25–31. [34] Diaz OC, Modaresghazani J, Satyro MA, Yarranton HW. Modeling the phase behavior of heavy oil and solvent mixtures. Fluid Phase Equilibr 2011;304(1–2):74–85. [35] Agrawal P, Schoeggl FF, Satyro MA, Taylor SD, Yarranton HW. Measurement and modeling of the phase behavior of solvent diluted bitumens. Fluid Phase Equilibr 2012;334:51–64. [36] Li H, Yang D. Phase behaviour of C3H8/n-C4H10/heavy-oil systems at high pressures and elevated temperatures. J Can Petrol Technol 2013;52(1):30–40. [37] Badamchi-Zadeh A, Yarranton, HW, Maini BB, Phase behavior and physical property modeling for vapex solvents: propane, carbon dioxide, and Athabasca bitumen. Paper SPE 165505 presented at SPE Heavy Oil Conference, Calgary, Alberta, Canada, 11–13 June 2013. [38] Yang P, Li H, Yang D, Phase behaviour of CH4-C3H8-heavy oil systems under reservoir conditions. Paper SPE 165502 presented at SPE Heavy Oil Conference, Calgary, Alberta, Canada, 11-13 June 2013. [39] Ghasemi M, Whitson CH, Phase behavior and viscosity modeling of Athabasca bitumen and light solvent mixtures. Paper SPE 165416 presented at SPE Heavy Oil Conference, Calgary, Alberta, Canada, 11-13 June 2013. [40] Twu CH. An internally consistent correlation for predicting the critical properties and molecular weights of petroleum and coal-tar liquids. Fluid Phase Equilibr 1984;16(2):137–50. [41] Gao J, Okuno R, Li HA. An experimental study of multiphase behavior for n-butane/ bitumen/water mixtures. SPE J 2017;22(03):783–98. [42] Eghbali S, Dehghanpour H. An experimental and modeling study of carbon dioxide/ bitumen and C 4/bitumen phase behavior at elevated temperatures using cold lake bitumen. SPE J 2018. [43] Kumar A, Okuno R. Reliable characterization of bitumen based on perturbation from n-alkanes for steam-solvent coinjection simulation. Fuel 2016;182:141–53. [44] Azinfar B, Haddadnia A, Zirrahi M, Hassanzadeh H, Abedi J. Phase behaviour of butane/bitumen fractions: experimental and modeling studies. Fuel 2018;220:47–59. [45] Robinson DB, Peng DY, The Characterization of the Heptanes and Heavier Fractions for the GPA Peng-Robinson Programs; Gas Processors Association, Research Report RR-28: Tulsa, Oklahoma, USA, 1978. [46] Mehra RK. The Computation of Multi-Phase Equilibrium in Compositional Reservoir Studies. University of Calgary; 1981. Ph.D. Thesis. [47] Li Y-K, Heavy Fraction Characterization and Hypothetical Component Selection for Oil and Gas Mixtures, Computer Modelling Group Research Report R12.04, May 1983. [48] Oellrich L, Ploecker U, Prausnitz JM, Knapp H. Equation-of-state methods for computing phase equilibria and enthalpies. Int Chem Eng 1981;21(1):1–15. [49] Katz DL, Firoozabadi A. Predicting phase behavior of condensate crude/oil systems using methane interaction coefficients. J Petrol Technol 1978;30(11):1649–55. [50] Peneloux A, Rauzy E, Freze R. A consistent correction for redlich-kwong-soave volumes. Fluid Phase Equilib 1982;8(1):7–23. [51] Pedersen KS, Milter J, Sørensen H. Cubic equations of state applied to ht/hp and highly aromatic fluids. SPE J 2004;9:186–92. [52] Kokal SL, Sayegh SG. Gas-saturated bitumen density predictions using the volumetranslated peng-robinson equation of state. J Can Petrol Technol 1990;29(05):31–9. [53] Pedersen KS, Fredenslund A, Christensen PL, Thomassen P. Viscosity of crude oils. Chem Eng Sci 1984;39(6):1011–6. [54] Pedersen KS, Fredenslund A. An improved corresponding states model for the prediction of oil and gas viscosities and thermal conductivities. Chem Eng Sci 1987;42(1):182–6. [55] Nghiem LX, Li YK. Computation of multiphase equilibrium phenomena of reservoir fluid. Fluid Phase Equilibr 1984;17(1):77–95. [56] Nghiem LX, Heidemann RA, General acceleration procedure for multiphase flash calculation with application to oil-gas-water systems. Paper presented at the 2nd European Symposium on Enhanced Oil Recovery, Paris, France, 8-10 November 1982. [57] Kariznovi M, Nourozieh H, Abedi J. Experimental determination of k-values and compositional analysis of liquid phases in the liquid-liquid equilibrium study of (Athabasca bitumen + ethane) systems. J Chem Eng Data 2013;58(6):1772–80. [58] Nourozieh H. Phase Partitioning and Thermo-physical Properties of Athabasca Bitumen/Solvent Mixtures PhD dissertation University of Calgary; 2013 [59] Riazi MR, Characterization and properties of petroleum fractions. 1st Edition, Philadelphia, PA, 2005. [60] Alboudwarej H, Akbarzadeh K, Beck J, Svrcek WY, Yarranton HW. Regular solution model for asphaltene precipitation from bitumens and solvents. AIChE J 2003;49(11):2948–56. [61] Riazi MR, Daubert TE. Prediction of the composition of petroleum fractions. Ind Eng Chem Proc DD 1980;19(2):289–94. [62] Whitson C. Characterizing hydrocarbon plus fractions. SPE J 1983;23(4):683–94. [63] Sturges HA. The choice of a class interval. J Am Stat Assoc 1926;21(153):65–6. [64] Lee BI, Kesler MG. A generalized thermodynamic correlation based on three-parameter corresponding states. AIChE J 1975;21(3):510–27. [65] Nourozieh H, Kariznovi M, Abedi J. Physical properties and extraction measurements for Athabasca bitumen + light hydrocarbon system: evaluation of pressure effect, solvent-to-bitumen ratio, and solvent type. J Chem Eng Data 2011;56(11):4261–7.

References [1] Kumar A, Okuno R. Characterization of reservoir fluids using an EOS based on perturbation from n-alkanes. Fluid Phase Equilibr 2013;358:250–71. [2] Whitson CH. Characterizing hydrocarbon plus fractions. SPE J. 1983;23:683–94. [3] Stamataki SK, Magoulas, KG, Characterization of Heavy Undefined Fractions, 2001 SPE International Symposium on Oilfield Chemistry, SPE 64996. [4] Montel F, Gouel PL, A new lumping scheme of analytical data for compositional studies. In: SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers, 1984. [5] Leibovici C, Stenby EH, Knudsen K. A consistent procedure for pseudo-component delumping. Fluid Phase Equilibr 1996;117(1–2):225–32. [6] Lomeland F, Harstad O. Simplifying the task of grouping fluid components in compositional reservoir simulation. SPE Comp App 1995;7(2):38–43. [7] Yarborough L. Application of the generalized equation of state to petroleum reservoir fluids, equation of state in engineering and research, in: K.C. Choa, R.L. Robinson Jr. (Eds.). Adv Chem Ser Am Chem Soc 1978;182:386–439. [8] Lawal AS, Van der Laan ET, Thambynayagam RKM, Four-parameter modification of the Lawal-Lake-Silberberg equation of state for calculating gas-condensate phase equilibria. In: SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers, 1985. [9] Pedersen KS, Thomassen P, Fredenslund A, AIChE Spring National Meeting New Orleans. LA, March 1988, pp. 6–10. [10] Guo TM, Du L, A New Three-Parameter Cubic Equation of State for Reservoir Fluids. III. Application to Gas Condensates. SPE-19374-MS, unsolicited, 1989. [11] Peneloux A, Jain C, Behar, E, Application of the Developed Redlich-Kwong Equation of State to Predict the Thermodynamic Properties of Condensate Gases. In: SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers, Jan 1979. [12] Whitson, CH, Torp SB, Evaluating constant volume depletion data. In: SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers, Jan 1981. [13] Zuo JY, Zhang D, Plus fraction characterization and PVT data regression for reservoir fluids near critical conditions. In: SPE Asia Pacific Oil and Gas Conference and Exhibition. Society of Petroleum Engineers, 2000. [14] Al-Meshari AA, McCain WD, An accurate set of correlations for calculating critical properties and acentric factors for single carbon number groups. In: SPE Technical Symposium of Saudi Arabia Section. Society of Petroleum Engineers, Jan 2006. [15] Hosein R, Dawe R. Characterization of trinidad gas condensates. Pet Sci Technol 2011;29(8):817–23. [16] Hosein R, Jagai T, Dawe RA. A method for predicting the phase behaviour of trinidad gas condensates. W Indian J Eng 2013;35(2):22–30. [17] Krejbjerg K, Pedersen KS, Controlling VLLE equilibrium with a cubic EOS in heavy oil modeling. In: Canadian International Petroleum Conference. Petroleum Society of Canada, Jan 2006. [18] Ghasemi M, Alavian SA, Whitson, CH, SPE Heavy Oil Conference and Exhibition held in Kuwait City. Kuwait, Dec 2011, pp. 12–14. [19] Mehrotra AK, Sarkar M, Svrcek WY. Bitumen density and gas solubility predictions using the Peng Robinson equation of state. AOSTRA J Res 1985;1(4):215–29. [20] Bishnoi PR, Heidemann RA, Shah MK, Calculation of the thermodynamic properties of bitumen system. In: The oil sands of Canada – Venezuela 1977, editors, David A. Redford and Alvin G. Winestock, Proceedings of the Canada-Venezuela Oil Sands Symposium 1977, Edmonton, Alberta, Canada, 248-254. [21] Kesler MG, Lee BI. Improve prediction of enthalpy of fractions. Hydrocarb Process 1976;55(3):153–8. [22] Fu C, Puttagunta R. Binary interaction coefficients of the modified soave-redlichkwong equation of state for bitumen-containing systems at in-situ conditions. J Can Petrol Technol 1986;25(5):38–43. [23] Fu CT, Puttagunta VR, Vilcsak G. Gas solubility of methane and ethane in Cold Lake bitumen at in situ conditions. J Can Petrol Technol 1988;27(4):79–85. [24] Mehrotra AK, Svrcek WY. Characterization of Athabasca bitumen for gas solubility calculations. J Can Petrol Technol 1988;27(6):107–10. [25] Mehrotra AK, Svrcek WY. Correlation and prediction of gas solubility in Cold Lake bitumen. Can J Chem Eng 1988;66(4):666–70. [26] Mehrotra AK, Patience GS, Svrcek WY. Calculation of gas solubility in Wabasca bitumen. J Can Petrol Technol 1989;28(3):81–3. [27] Frauenfeld TWJ, Kissel G, Zhou S, PVT and viscosity measurements for Lloydminster-Aberfeldy and Cold Lake blended oil systems. Paper SPE 2002-110 presented at Canadian International Petroleum Conference, Calgary, Alberta, Canada, 11-13 June 2002. [28] Kariznovi M, Nourozieh H, Abedi J. Bitumen characterization and pseudocomponents determination for equation of state modeling. Energ Fuel 2010;24(1):624–33. [29] Mehrotra AK, Svrcek WY. Viscosity, density and gas solubility data for oil sand bitumens. part I: Athabasca bitumen saturated with CO and C2H6. AOSTRA J. Res. 1985;1(4):263–8.

447