A Group Contribution Pseudocomponent Method for Phase

Aug 18, 2015 - We present a group contribution pseudocomponent (PC) method for mixtures of petroleum fluids and a solvent based on the characterizatio...
7 downloads 19 Views 3MB Size
Article pubs.acs.org/IECR

A Group Contribution Pseudocomponent Method for Phase Equilibrium Modeling of Mixtures of Petroleum Fluids and a Solvent Ping He* and Ahmed F. Ghoniem Department of Mechanical Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, United States ABSTRACT: We present a group contribution pseudocomponent (PC) method for mixtures of petroleum fluids and a solvent based on the characterization method (CM) for petroleum fractions and the predictive Peng−Robinson 1978 (PPR78) model. The objective of this study is to model the phase equilibrium of crude oil and a solvent consisting of one or multiple pure components. We develop a structural analysis of crude oil to bridge the CM and the PPR78 model, which predicts the temperature-dependent binary interaction parameters (BIP) among the PCs and pure components. The predicted results using our method agree well with experimental data for mixtures of oil sand bitumen and propane, CO2, or water. The approach is applicable to phase equilibrium predictions of mixtures of crude oil and water over a wide range of temperatures (273−644 K or higher as reasonable extrapolations) and pressures (4 kPa−26 MPa). most widely used model. Riazi and Daubert16,17 developed a systematic framework of correlation methods to estimate the thermodynamic properties as well as the physical and structural properties of petroleum fractions, such as the paraffins, naphthenes, and aromatics (PNA) compositions, refractive index, sulfur content, carbon and hydrogen weight ratio, etc. Compared to the real molecule approach, the petroleum fraction approach may produce larger deviations in the estimated properties. To overcome this, Peng18 introduced a molecular-type homologous series (MTHS) matrix to describe crude oils by grouping different molecule types into lumps, and Gomez-Prado et al.19 developed a modified MTHS model for heavy crude oils. Reiter et al.20 proposed to replace the PCs with real components using a property matching algorithm. In addition to the critical properties provided by these thermodynamic models, phase equilibrium calculations require the binary interaction parameters (BIP), which are critical to most polar molecules and obtained from experimental data. Gao et al.21 proposed a correlation to estimate the BIPs in the Peng−Robinson (PR) EoS between two light hydrocarbons, in which good matches with experiments were achieved for vapor−liquid equilibrium of light hydrocarbon mixtures. Diaz et al.22 used a similar approach to estimate the BIPs in the PR-EoS with a fitted parameter from experimental phase equilibrium data. Both correlations have the limitation that they cannot handle highly elevated temperatures and pressures, such as the condition of supercritical water. Jaubert’s group23−36 has been developing a group contribution method called the predictive Peng−Robinson 1978 (PPR78) model to estimate temperature-dependent BIPs in the PR-EoS between any two species such as hydrocarbons, mercaptans, H2, CO2, H2S, water, N2, alkenes, CO, He, Ar, etc. The PPR78 model has been built upon a theoretical relation between the mixing rules and the

1. INTRODUCTION Phase equilibrium modeling of crude oil is essential in several petroleum industrial processes, e.g., enhanced oil recovery using gas injection of sub- or supercritical CO2,1,2 upgrading of crude oil using supercritical water,3−7 etc. To model the thermodynamic properties of whole crude, two approaches have commonly been used: (1) the real molecule approach and (2) the characterization method (CM). In early studies of the real molecule approach, a synthetic crude oil was constructed using limited numbers of hydrocarbons8−10 by matching the phase behavior such as the bubble point of the actual crude. Using a synthetic crude oil greatly reduces the complexity of describing the crude composition and its properties, but the choice of representative species remains a major difficulty especially for the heavy fractions. To avoid choosing limited molecule types and gain higher accuracy, Klein’s group11,12 developed a molecular-level model to describe heavy crude oils using 5000+ unique molecules with corresponding probability density functions (PDFs) and reaction networks. Klein’s model is considered as the state-of-art method, because the PDFs of 5000+ molecules have been optimized by minimizing the difference between the computed and measured values for chemical and thermodynamic properties. Nevertheless, using a large number of molecules increases the computing time of the phase equilibrium calculations, for instance, time for a flash calculation is proportional to the square of the number of molecules. The petroleum characterization method (CM) approach has been developed since the early 20th century because of its extreme importance for refinery analysis. In the CM approach, each petroleum distillation fraction is treated as a pseudocomponent (PC), also called pseudospecies. Watson and Nelson (1933)13 analyzed the critical and thermal properties of petroleum distillation fractions and defined a Watson characterization factor for crude oils to describe its paraffinic, naphthenic, or aromatic properties. Since then, several CMs have been proposed to estimate the critical properties, molecular weights, molar enthalpies, etc. for crude PCs. Among these, the Lee−Kesler correlations14,15 have been the © 2015 American Chemical Society

Received: Revised: Accepted: Published: 8809

July 10, 2015 August 11, 2015 August 18, 2015 August 18, 2015 DOI: 10.1021/acs.iecr.5b02516 Ind. Eng. Chem. Res. 2015, 54, 8809−8820

Article

Industrial & Engineering Chemistry Research

Figure 1. Flowchart of the group contribution pseudocomponent method developed in this study.

knowledge, this work is the first to model the phase behavior of crude oil and water under elevated temperatures and pressures. Compared to previous work,44 (1) we use direct structural analysis to estimate the occurrence of functional groups in PCs instead of correlations; (2) we expand the model to predict the phase behaviors of mixtures of petroleum fluids and a solvent instead of the phase behavior of the petroleum fluids alone; and (3) we validate our approach by comparing the predictions with experimental data. The rest of the paper is organized as follows. The method and submodels are presented in Section 2 using Athabasca bitumen as an example to demonstrate the application for the method. In Section 3, the phase equilibrium results of mixtures of bitumen and propane, CO2 or water are compared with experimental data to validate our approach. Finally, the conclusions are presented in Section 4.

excess Gibbs energy models;37,38 thus, it overcomes the wellknown difficulty in BIP correlations for different species and temperatures. The model covers a wide range of pressures and temperatures. The PPR78 model could also be easily modified for another cubic equation,39 such as the Soave−Redlich− Kwong EoS. Hajiw et al.40 proposed a cubic-plus-association (CPA) EoS using the BIPs calculated from the PPR78 model, which modified only the BIPs of water and increased the accuracy of phase equilibrium calculations of aqueous mixtures. Gros et al.41 developed a group-contribution-associating equation of state (GCA-EoS) to calculate phase equilibria of hydrocarbons, alcohols, and water, and Ferreira et al.42 and Sánchez et al.43 extended the model to mixtures containing acids, esters, ketones, and aromatics with water and alcohols. The GCA-EoS is a SAFT-type (Statistical Associating Fluid Theory) EoS, whose framework is different from the PPR78 model and CPA-EoS but with the same phase equilibrium capabilities. A combination of phase equilibrium models (PPR78, CPAEoS, or GCA-EoS) and thermodynamic property models for crude oils (either the real molecule approach or the CM approach) can be utilized to model phase behavior of a mixture of crude oil and a pure component such as CO2 or water. Jaubert and co-workers,44 the authors of the PPR78 model, have recently published an example for how to do that, in which they used a group contribution method for critical properties and acentric factor to estimate the occurrences of functional groups in crude PCs. In this work, we propose an alternative approach based on the development of a direct structural analysis of crude PCs to link the petroleum CM with the PPR78 model. Both Jaubert et al. and our proposed method can be categorized as the group contribution pseudocomponent methods for phase equilibrium modeling of petroleum fluids. The objective of this work is to develop a generalized method to model the phase equilibrium of mixtures of crude oil and solvent such as CO2 and water over a wide range of temperatures and pressures. Compared to the Gao et al.21 and Diaz et al.22 correlations, the current model is capable of predicting phase equilibrium of crude oil mixtures over a wide range of temperatures and pressures, especially in the range of supercritical water. However, the current model is not capable of analyzing asphaltene precipitation. To the best of our

2. THERMODYNAMIC MODELS OF PETROLEUM FLUIDS As a complex mixture of hydrocarbons and other species, crude oil is often modeled using PCs based on the true boiling points in the distillation analysis for refinery processes. A PC represents a petroleum distillation fraction within a boiling range and is treated as a single species in the analysis. In this study, a group contribution pseudocomponent method for crude oil is developed to model the phase equilibrium of crude oil and a pure component. Figure 1 shows the flowchart describing our modeling approach. Namely, a crude sample is processed through distillation, and each distillation fraction is treated as a PC. For each PC, the specific gravity, SG, is measured or estimated using an average Watson factor if available in the crude assay. The critical properties, Tc and Pc, acentric factor, ω, and molecular weight, M, are estimated from the boiling point, Tb, and specific gravity, SG, using one of the CMs, e.g., the Lee−Kesler correlations (cf. Section 2.2). The PNA analysis is performed using a PIONA analyzer to get the compositions of paraffins, naphthenes, and aromatics or using the Riazi−Daubert correlations17 calculated from Tb and SG. The C/H ratios and sulfur content are measured using an elemental analyzer or estimated using the Riazi−Daubert correlations calculated from SG and M. The aromaticity, i.e., 8810

DOI: 10.1021/acs.iecr.5b02516 Ind. Eng. Chem. Res. 2015, 54, 8809−8820

Article

Industrial & Engineering Chemistry Research the ratio of aromatic carbon number over the total carbon number, is measured using the NMR spectroscopy45,46 or estimated based on NMR measurements available in the literature for the same or similar crude sample. Tb, SG, M, C/H ratios, sulfur content, and aromaticity are used in structural analysis to estimate the occurrences of the functional groups in the paraffins, naphthenes, and aromatics, the results of which are then averaged using the PNA composition to obtain the occurrences of functional groups in each PC (cf. Section 2.3). Using the thermodynamic properties Tc, Pc, ω, and occurrences of functional groups of crude PCs and pure components, the PPR78 model predicts temperature-dependent BIPs among PCs and pure components. Flash calculations for phase equilibrium of mixtures of crude and pure components can then be performed using Tc, Pc, ω, and BIPs. The rest of this section is organized in three parts: (1) The thermodynamic models for mixtures of crude oil and a solvent, (2) the characterization method to estimate Tc, Pc, ω, and M of crude PCs, and (3) the structural analysis to obtain the BIPs among PCs and pure components. 2.1. Thermodynamics Models for Mixtures of Crude Oil and a Solvent. 2.1.1. Peng−Robinson EoS with Mixing Rules. The cubic Peng−Robinson equation of state (PR-EoS) with standard van der Waals one-fluid mixing rules47 is used in this study p=

am =

am(T ) RT − 2 V − bm V + 2bmV − bm 2

∑ ∑ xixj(1 − kij) i

bm =

Ωb =

X=

bi am (Z − 1) − 8 bmRT bm 2 )bmρ ⎤⎡⎢ 2 ∑j xj(1 − kij) aiaj b ⎤ − i⎥ ⎥ 2 )bmρ ⎦⎢⎣ am bm ⎥⎦

ln(ϕi)̂ = −ln[Z(1 − bmρ)] + ⎡ 1 + (1 + ln⎢ ⎣ 1 + (1 −

(4)

where Z = PV/RT is the compressibility and ρ = 1/V is the molar concentration. The fugacity criteria for phase equilibrium is f1,̂ i = f2,̂ i

x1, iϕ1,̂ i = x 2, iϕ1,̂ i

or

(5)

where “1” denotes phase 1 and “2” phase 2 in the equilibrium state. In a binary system, the equilibrium composition is calculated using eq 5 and the unity condition (∑x1,i = ∑x2,i = 1). For a multicomponent system (≥3), there is an infinite number of equilibrium states defined by the composition. In this case, species conservations are used in order to compute a unique solution. Species conservation is expressed as follows zi = cx1, i + (1 − c)x 2, i

(6)

where zi is the average mole fraction in the mixture and c is the mole fraction of phase 1. 2.1.3. Predictive Peng−Robinson 1978 Model. The BIPs, kij, in the one-fluid mixing rules are critical to the phase equilibrium calculations. Instead of using application-specific BIP models, in this study, a group contribution method named the predictive Peng−Robinson 1978 (PPR78) model23−36 is used to calculate temperature-dependent BIPs among crude PCs, propane, CO2, water, etc. The PPR78 model data (Akl and Bkl) were correlated from large body of phase equilibrium experimental data of the listed species and covered a wide range of temperatures (273700 K) and pressures (4 kPa40 MPa). In PPR78, a BIP is calculated using the follow equation

2

Fij

kij =

RTc , i

−2 −

⎛ ⎜ ⎝ 2

(3)

ai bi



aj bj

⎞2 ⎟ ⎠

aiaj bibj

(7)

where

where Ωa =

≈ 0.253076587

2.1.2. Phase Equilibrium Criteria. The fugacity coefficient is calculated using the Peng−Robinson EoS and the one-fluid mixing rules48

(2)

Pc , i

6 2 −8

ωi3

where R is the universal gas constant, V is molar volume of the mixture, xi is molar fraction of species i, and kij is the temperature-dependent binary interaction parameter (BIP) describing the interaction between two components in the mixture and obtained using the PPR78 model. The values of the BIPs depend on the molecular properties of the components and capture their molecular interaction. A positive kij indicates reduced attraction and leads to repulsion. A negative kij indicates attraction. Positive BIPs favor phase separations, while negative BIPs favor mixing.48 BIPs are critical to the calculations of the phase equilibrium of mixtures. Calculating these values for pseudocomponents is a major challenge, and one such approach is described in this paper. ai and bi are calculated from the critical temperature, Tc, critical pressure, Pc, and acentric factor, ω, of each species:

bi = Ωb

3

mi = 0.379642 + 1.48503ωi − 0.164423ωi2 + 0.0166667

∑ xibi

⎞⎤ ⎟⎥ ⎟⎥ ⎠⎦

6 2 +8 − 3

for ωi > 0.491, mi is used in the following equation for heavy hydrocarbons49

aiaj

T Tc , i

3

mi = 0.374640 + 1.54226ωi − 0.26992ωi 2

j

⎛ R2Tc , i 2 ⎡ ⎢1 + m ⎜1 − ai = Ωa i⎜ Pc , i ⎢⎣ ⎝

−1 +

and for ωi ≤ 0.491

(1)

i

X ≈ 0.0777960739 X+3

Ng

Fij =

8(5X + 1) ≈ 0.457235529 49 − 37X

Ng

⎛ 298.15 ⎞(Bkl / Akl − 1) ⎜ ⎟ T ⎠

∑ ∑ (αik − αjk)(αil − αjl)Akl ⎝ k=1 l=1

(8) 8811

DOI: 10.1021/acs.iecr.5b02516 Ind. Eng. Chem. Res. 2015, 54, 8809−8820

Article

Industrial & Engineering Chemistry Research

certain functional groups in the PCs as represented by αik, assuming that the corresponding constants Akl and Bkl are known for these mixtures within the temperature and pressure range of interest. Therefore, and in order to estimate the BIPs for crude PCs and their mixtures with pure components, we develop a structural analysis to estimate the occurrences of the following functional groups whose data are available in PPR78 (a subset of the 24 functional groups defined in PPR78): (1) alkanic group (CH3, CH2, CH, C), (2) cycloalkanic groups (CH2,cyclic, CH/Ccyclic), (3) aromatic groups (CHaro, Caro, Cfused aro) (4) sulfhydryl group (SH) Each PC is assumed to be a mixture of three types of species: paraffins, naphthenes, and aromatics (PNA). The PNA analysis measures or estimates the fractions of paraffinic, naphthenic, and aromatic species in each PC, and the structural analysis is performed to estimate the occurrence of functional groups in paraffinic, naphthenic, and aromatic species in that PC. The molecular weights of paraffinic, naphthenic, and aromatic species of the same PC are assumed to be equal. Thus, the molar fractions of PNA in each PC are the same as the mass fractions. Let yPi, yNi, and yAi denote such fractions of paraffinic, naphthenic, and aromatic species of PC i and αkPi, αkNi, and αkAi denote the occurrence of functional group k in the paraffinic, naphthenic, and aromatic species of PC i, respectively. The occurrence of functional group k in PC i is averaged as follows αki = αkPi·yPi + αkNi·yNi + αkAi·yAi (9)

Ng = 24 is the number of functional groups (see Section 2.3 for definition of these groups) defined in PPR78, Akl and Bkl are constants correlated from experimental data, and αik is the occurrence of group k in molecule i. The occurrences are calculated using the structural analysis which is described in Section 2.3. Although the PPR78 model includes the interaction parameters for the groups included in this work, some of them have to be refitted. This is due to the fact that the estimated chemical structures of the PCs used in this study are far beyond the original correlation database of the PPR78 model. For example, the largest aromatic ring number of PCs used in this study is close to 30, while 3 is the largest number in the original database of the PPR78 model. Thus, in this study, the estimated BIPs between the PCs and water are corrected as shown in Section 3.4. 2.2. Characterization of Athabasca Bitumen. The Athabasca oil sand bitumen is selected as an example to demonstrate the application of our apporach because of the availability of distillation and phase behavior data in the literature. The boiling points Tb, specific gravities SG, and mole fractions of the PCs of Athabasca bitumen created by Diaz et al.22 from distillation measurements are used in the following analysis. Note that Diaz et al. artificially extended the boiling points beyond 617 °C to account for asphaltenes. Lee-Kesler correlations14,15 are used to calculate Tc, Pc, ω, and M, which are compared with the results of Diaz et al.22 as shown in Figure 2.

The occurrences are then used in eq 8 to calculate the BIPs among PCs and pure components. The rest of this subsection is used to describe the structural analysis developed in the current study and is organized as follows: (1) estimations for the PNA compositions, C/H ratios and sulfur contents in the crude oil (taken here as Athabasca oil sand bitumen), (2) the structural analysis for paraffins, (3) the structural analysis for naphthenes, and (4) the structural analysis for aromatics. 2.3.1. PNA Compositions, C/H Ratios, and Sulfur Contents in Bitumen. The PNA compositions yP, yN, and yA (cf. Figure 3), C/H weight ratios CH, and sulfur contents %S (cf. Figure 4a) in the bitumen PCs (shown in Table 1) are estimated using interpolations and extrapolations from the crude assay of Athabasca bitumen published by ExxonMobil.50 (Note that for simplicity, the sulfur content is assumed to exist only in the aromatic molecules, which does not affect our analysis because the occurrence of functional groups is averaged using the PNA compositions in the end of the structural analysis.) By neglecting oxygen, nitrogen, metals, etc. in bitumen, the carbon number of each PC is estimated as follows (cf. Figure 4b)

Figure 2. Properties of crude PCs created by Diaz et al.22 (crosses) from distillation data and calculated using Lee−Kesler correlations (circles), (a) critical temperatures, (b) critical pressures, (c) acentric factors, and (d) molecular weights.

NC =

In the current study, Lee−Kesler results of Tc, Pc, and ω are used for bitumen PCs. Note that Lee-Kesler correlations are known to predict smaller molecular weights for the heavy fractions of crude oil. As a result, and to satisfy the mass distributions in the bitumen sample, Diaz et al. results of M are used in the current study. The selected pseudocomponents and their corresponding properties are listed in Table 1. 2.3. Structural Analysis. As described in Section 2.1.3, the PPR78 model estimates temperature-dependent BIPs using group contributions of 24 functional groups, which are necessary for the phase equilibrium calculations. According to eqs 7 and 8, these parameters depend on the occurrence of

M(1 − %S /100) 12.011(1 + 1/CH )

(10)

2.3.2. Structural Analysis for Paraffins. Paraffins consist of normal and branched alkanes. Thus, four types of functional groups exist in paraffins: CH3, CH2, CH, and C. We assume that no highly branched alkanes exist in Athabasca bitumen, and hence the C group is neglected. The total paraffin weight fraction in the current bitumen sample is ε ⎩

molecules. We further assume that the average carbon number per aromatic ring is CperRing = 4 −

2 arctan[0.1(CA − 6)] π

(28)

(35)

Thus, the aromatic ring number (cf. Figure 7b) is calculated using Raro =

CA − 6 +1 CperRing

The numbers of CH2,cyclic and CHcyclic groups are assumed to be

(29)

Accordingly, the fused aromatic carbon number per molecule

(30)

and the nonfused aromatic carbon number per molecule is NC,nonfused aro = CA − NC,fused aro

(36)

total NCH,cyclic = 0.4Ccyclic

(37)

The numbers of alkanic functional groups are assumed as follows

is NC,fusedaro = 2(Raro − 1)

total NCH2,cyclic = 0.6Ccyclic

⎡ ⎤ π NCH3 = Calkyl ⎢1 − arctan(Calkyl − 1)⎥ ⎣ ⎦ 1.5

(38)

NCH = 0.1Calkyl

(39)

(31)

NCH2 = Calkyl − NCH3 − NCH

The nonfused aromatic carbons are categorized into two groups: CHaro and Caro, the latter has attached groups, e.g., an alkyal chain or a naphthenic ring. The number of the Caro group is assumed to be the same as the aromatic ring number Raro, which means that on average each aromatic ring connects to an alkyal chain or half of a naphthenic ring. Thus, we have the following relations: NC,aro = Raro

(32)

NCH, aro = NC,nonfusedaro − NC, aro

(33)

(40)

where the alkyl carbon number Calkyl = C − The occurrences of functional groups in the aromatics calculated using above equations are shown in Figure 8. sat

Ctatal cyclic.

The above groups are the aromatic functional groups in the aromatic species. The alkyl chains and naphthenic rings attached to the aromatic core are analyzed in the following. The naphthenic carbon number Ctatal cyclic is assumed to be 10% of total saturated carbons Csat in the aromatics when Csat > 5, and sat Ctatal ≤ 5 as shown in the following equation: cyclic = 0 when C total Ccyclic = 0.1C sat × Heaviside(C sat − 6, 1)

(34) Figure 8. Occurrences of functional groups in the aromatic part of bitumen PCs. Note that the jump at 617 °C is related to data extrapolation for asphaltenes.

where Csat = NC − CA and the Heaviside function is used to generate a smooth transition in the region of 5 ≤ Csat ≤ 7 8815

DOI: 10.1021/acs.iecr.5b02516 Ind. Eng. Chem. Res. 2015, 54, 8809−8820

Article

Industrial & Engineering Chemistry Research 2.3.5. Functional Groups in Bitumen. The data in Section 2.3.22.3.4 are now combined to estimate the average occurrences of functional groups for bitumen PCs calculated using eq 9, and the results are shown in Figure 9. These results are used next to estimate the BIPs for the bitumen, and for its mixtures with three different pure components.

Figure 10. BIPs among selected PCs calculated in the current study. Diaz et al.22 estimated the following constant BIPs for the same PCs: 0.029 (Malt_ABP[285]C and Asph_6), 0.023 (Malt_ABP[323]C and Asph_6), 0.008 (Malt_ABP[359]C and Asph_4), 0.014 (Malt_ABP[359]C and Asph_1), 0.008 (Malt_ABP[395]C and Asph_3), and 0.0004 (Asph_4 and Asph_6). Figure 9. Occurrences of functional groups in Athabasca bitumen PCs. Note that the jump at 617 °C is related to data extrapolation for asphaltenes.

Figures 12. In addition, we note that the current model is not capable to analyze asphaltene precipitation. 3.3. Mixing of Bitumen and CO2. BIPs of PCs and CO2 are shown in Figure 13. The BIP of a PC and CO2 decreases with temperature below ∼350 K and increases slightly above ∼350 K. The current study only validates the BIP of a bitumen PC and CO2 for temperature below 369 K. Predictions of saturated CO2 compositions in the bitumen show good agreement with the experimental data55,56 as shown in Figure 14. Similar results to ours are predicted by Diaz et al.22 for mixing bitumen with propane or CO2. In the following, we show that our approach can model crude mixtures with water under elevated temperatures. 3.4. Mixing of Bitumen and Water. The BIPs of bitumen PCs and water are the most interesting and important for the applications of this study. The BIPs directly calculated using the PPR78 model decrease with temperature, reaching a minimum at a value close to the critical temperature of water (the temperature at which the BIP of a PC and water reaches a minimum rises slowly as the molecular weight of the PC increases), and the same phenomenon is found for CO2 as shown in Figure 13. At higher temperatures, the calculated BIPs increase rapidly. This is not physically correct, because using the original BIPs of PCs and water calculated using the PPR78 model we find that the solubility of water in a PC is decreasing with temperature after a certain point. It is contrary to the fact that the solubility of water in a PC is expected to increase monotonically with temperature before full miscibility is achieved. Thus, the following guidelines are used to correct the BIP of a PC and water near and beyond the temperature where the BIP reaches a minimum: (1) BIPs of PCs and water should decrease monotonically toward zero or negative values; (2) At higher temperatures, e.g., above 500 K, the BIP of a heavier PC and water should be higher than that of a lighter one; (3) Corrections should make the phase equilibrium predictions close to experimental data.

3. RESULTS AND VALIDATIONS For the purpose of model validation, mutli-component flash calculations of phase equilibrium have been performed for mixtures of bitumen and propane, CO2, or water, to compare with experimental measurements54−57 of saturate compositions in bitumen. BIPs of bitumen PCs and of bitumen PCs and pure components are calculated using the PPR78 model (cf. eqs 7 and 8) based on the occurrence of functional groups (cf. Figure 9). In each case, the calculated BIPs are shown, and then the phase equilibrium results are compared with the experiments. 3.1. BIPs of Bitumen Pseudocomponents. The calculated BIPs among PCs are in the range of −0.035 to 0.004 over the temperature range of 300700 K, as shown for 6 selected pairs in Figure 10. The BIP between the lightest and heaviest PCs has the largest negative value indicating that their attraction is the strongest among the different PCs. Their absolute value decreases with temperature indicating less attraction. The absolute value of the BIP decreases as the boiling points of two PCs get closer. The BIPs among similar PCs are close to 0 indicating that similar PCs have negligible deviation from the combination rule. The same trend and similar values were predicted in Diaz et al.’s22 work, although they estimated positive BIPs of the same bitumen PCs. Negative BIPs of maltenes and asphaltenes were used in the study of asphaltene precipitation in the crude oil,53 similar to what is predicted in our study. However, because the BIPs of the crude PCs are small in magnitude, they do not influence the phase behaviors of the mixtures of crude oil and pure components. 3.2. Mixing of Bitumen and Propane. The calculated BIPs between PCs and propane over 300600 K are shown in Figure 11. Generally, the BIP of a PC and propane decreases when temperature increases and when the PC becomes heavier. Predictions of saturated propane compositions in the bitumen show good agreement with the experimental data54 as shown in 8816

DOI: 10.1021/acs.iecr.5b02516 Ind. Eng. Chem. Res. 2015, 54, 8809−8820

Article

Industrial & Engineering Chemistry Research

Figure 11. BIPs between PCs and propane calculated in the current study.

Figure 12. Saturate propane composition in bitumen at different temperatures and pressures (experimental data are from the literature54).

Figure 14. Saturate CO2 composition in Athabasca bitumen at different temperatures and pressures (experimental data are from the literature40,41).

We apply the above guidelines to choose an extrapolation starting temperature T0 for each PC (cf. Table 2), i.e., (i) when temperature is lower than T0, the BIP calculated using the original PPR78 data are used, and (ii) when temperature is higher than T0, the BIP is extrapolated using the value and slope at T0. Figure 15 shows that the corrected BIPs and the BIPs of maltenes and water are in the region of smaller BIP values compared to those of asphaltenes and water. This correction is critical in predicting results consistent with the phase equilibrium experimental data shown in Figure 16.

Our results using the corrected BIPs based on the PPR78 predictions and estimated functional groups show good agreement with the experimental data57 for the solubility of water in bitumen fractions as shown in Figure 16. We note that the self-association behavior of water may cause relatively larger deviations in the model predictions than in the cases of propane or CO2. Modeling the association behavior, the CPA-EoS40 and GCA-EoS43 could potentially increase the accuracy in phase behavior predictions of crude and water. Nevertheless, our results show that the deviation from the experimental data of water solubility in bitumen is within 2 wt % and decreases when

Figure 13. BIPs between PCs and CO2 calculated in the current study. Note that these results are consistent with the BIPs of alkanes and CO2 in Péneloux and co-workers’ work,37 in which they validated the BIP values with experiments up to 663 K. 8817

DOI: 10.1021/acs.iecr.5b02516 Ind. Eng. Chem. Res. 2015, 54, 8809−8820

Article

Industrial & Engineering Chemistry Research Table 2. Correction Starting Temperature T0 for BIPs of PCs and Water pseudocomponent

T0 (K)

pseudocomponent

T0 (K)

pseudocomponent

T0 (K)

Malt_ABP[285]C Malt_ABP[323]C Malt_ABP[359]C Malt_ABP[395]C Malt_ABP[432]C Malt_ABP[456]C

610 607 599 587 572 560

Malt_ABP[469]C Malt_ABP[482]C Malt_ABP[493]C Malt_ABP[525]C Malt_ABP[574]C Malt_ABP[618]C

555 550 546 532 508 489

Asph_1 Asph_2 Asph_3 Asph_4 Asph_5 Asph_6

499 486 472 458 442 432

assumptions in two aspects: (1) the paraffinic branch ratios and (2) aromaticity. Since fractions of paraffins are small in the Athabasca bitumen, the assumptions of paraffins are not analyzed here. We denote Case 1 as the case presented in Section 2 and former subsections of Section 3. Case 2 is obtained by increasing the paraffinic branch ratios based on Case 1. The corresponding assumptions are modified as follows: (1) Methyl branch number in naphthenes is greatly increased compared to eq 18 NMe − b = 0.2Cakyl

(41)

(2) Naphthenic carbon number in aromatics is doubled Figure 15. BIPs of PCs and water calculated in the current study. The curves of PCs are shown alternately in solid lines and dashed lines according to the orders of the PCs listed in Tables 1 and 2. Note that the arrows show the curves from light fractions toward heavy fractions, i.e., from the PC, Malt_ABP[285]C, toward the PC, Asph_6. Specifically, the two blocks of curves beyond 500 K are of maltenes (lower) and asphaltenes (upper), respectively.

total Ccyclic = 0.2C sat × Heaviside(C sat − 6, 1)

(42)

(3) Methyl branch number (the same as CH number) in aromatics is doubled as well NCH = 0.2Calkyl

(43)

Case 3 is obtained by only increasing the aromatic carbon number CA using ⎧ 6, Tb < 218°C ⎪ CA = ⎨ ⎪ 1.3 −5 ⎩ 6 + 1.2 × 10 M(Tb − 218) , Tb ≥ 218°C (44)

CA in Case 3 increases 19.68% on average compared to Case 1. Functional groups of Cases 2 and 3 are estimated using the same structural analysis in Case 1 except for the operations mentioned above. BIPs among the PCs and pure components are calculated using the PPR78 model and the occurrence of functional groups estimated in the current study. Phase equilibrium calculations show that the sensitivity of the structural assumptions is small for the mixing of bitumen with (1) propane, (2) CO2, and (3) water. Here we present the results of CO2 in bitumen in Figure 17. Case 3 of a large aromaticity has relatively larger deviation from Case 1 than Case 2 of a large paraffinic branch ratio. Although the deviation of Case 3 is small as shown in Figure 17, these results indicate that it is important to use the correct aromatic carbon number based on experimental measurements, e.g., the NMR spectroscopy, in our proposed method.

Figure 16. Saturate water composition in bitumen at different temperatures and pressures (experimental data are from the literature57). The lower part of each curve is vapor−liquid equilibrium (VLE), and the upper part is liquid−liquid equilibrium (LLE). Vapor− liquid−liquid equilibrium (VLLE) occurs obviously in our predictions at 623 and 633 K. Note that the experimental data on the top right is VLLE data.

4. CONCLUDING REMARKS In order to model the phase equilibrium of crude oil and supercritical water, we extend the crude characterization method (CM) to handle phase equilibrium calculations using the predictive Peng−Robinson 1978 (PPR78) model and develop a structural analysis approach to bridge the gap between crude CM and the PPR78 model. Thus, we propose a

temperature and pressure approaches the supercritical point of water, at which the association behavior of water is reduced. 3.5. Sensitivity Study of Structural Assumptions. Structural assumptions have been made in this study when no experimental measurements are available for the crude sample. We examine the sensitivity of the structural 8818

DOI: 10.1021/acs.iecr.5b02516 Ind. Eng. Chem. Res. 2015, 54, 8809−8820

Article

Industrial & Engineering Chemistry Research

(5) Choi, K.-H. Process to upgrade whole crude oil by hot pressurized water and recovery fluid. U.S. Patent US 12/633,232, 2009. (6) Ates, A.; Azimi, G.; Choi, K.-H.; Green, W. H.; Timko, M. T. The role of catalyst in supercritical water desulfurization. Appl. Catal., B 2014, 147, 144−155. (7) Timko, M. T.; Ghoniem, A. F.; Green, W. H. Upgrading and desulfurization of heavy oils by supercritical water. J. Supercrit. Fluids 2015, 96, 114−123. (8) Danesh, A.; Xu, D. H.; Todd, A. C. Comparative study of cubic equations of state for predicting phase behavior and volumetric properties of injection gas-reservoir oil systems. Fluid Phase Equilib. 1991, 63, 259−278. (9) Ruffier-Meray, V.; Ungerer, P.; Carpentier, B.; Courcy, J. P. Fractionation of hydrocarbons between oil and gas phases. Oil Gas Sci. Technol. 1998, 53, 379−390. (10) Turek, E. A.; Metcalfe, R. S.; Yarborough, L.; Robinson, R. L. Phase equilibria in CO2 multicomponent hydrocarbon systems Experimental data and an improved prediction technique. SPEJ, Soc. Pet. Eng. J. 1984, 24, 308−324. (11) Horton, S. R.; Hou, Z.; Moreno, B. M.; Bennett, C. A.; Klein, M. T. Molecule-based modeling of heavy oil. Sci. China: Chem. 2013, 56, 840−847. (12) Nigam, A.; Klein, M. T. A mechanism-oriented lumping strategy for heavy hydrocarbon pyrolysis - Imposition of quantitative structurereactivity relationships for pure components. Ind. Eng. Chem. Res. 1993, 32, 1297−1303. (13) Watson, K. M.; Nelson, E. F. Improved methods for approximating critical and thermal properties of petroleum fractions. Ind. Eng. Chem. 1933, 25, 880−887. (14) Lee, B. I.; Kesler, M. G. Generalized thermodynamic correlation based on 3-parameter corresponding states. AIChE J. 1975, 21, 510− 527. (15) Kesler, M. G.; Lee, B. I. Improve prediction of enthalpy of fractions. Hydrocarbon Process. 1976, 55, 153−158. (16) Riazi, M. R.; Daubert, T. E. Characterization parameters for petroleum fractions. Ind. Eng. Chem. Res. 1987, 26, 755−759. (17) Riazi, M. R.; Daubert, T. E. Prediction of molecular-type analysis of petroleum fractions and coal liquids. Ind. Eng. Chem. Process Des. Dev. 1986, 25, 1009−1015. (18) Peng, B. Molecular Modelling of Petroleum Processes; University of Manchester Institute of Science and Technology: 1999. (19) Gomez-Prado, J.; Zhang, N.; Theodoropoulos, C. Characterisation of heavy petroleum fractions using modified molecular-type homologous series (MTHS) representation. Energy 2008, 33, 974− 987. (20) Reiter, A. M.; Wallek, T.; Mair-Zelenka, P.; Siebenhofer, M.; Reinberger, P. Characterization of crude oil by real component surrogates. Energy Fuels 2014, 28, 5565−5571. (21) Gao, G. H.; Daridon, J. L.; Saintguirons, H.; Xans, P.; Montel, F. A simple correlation to evaluate binary interaction parameters of the Peng-Robinson equation of state - Binary light-hydrocarbon systems. Fluid Phase Equilib. 1992, 74, 85−93. (22) 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, 74−85. (23) Jaubert, J. N.; Mutelet, F. VLE predictions with the PengRobinson equation of state and temperature dependent k(ij) calculated through a group contribution method. Fluid Phase Equilib. 2004, 224, 285−304. (24) Jaubert, J. N.; Vitu, S.; Mutelet, F.; Corriou, J. P. Extension of the PPR78 model (predictive 1978, Peng-Robinson EOS with temperature dependent k(ij) calculated through a group contribution method) to systems containing aromatic compounds. Fluid Phase Equilib. 2005, 237, 193−211. (25) Vitu, S.; Jaubert, J. N.; Mutelet, F. Extension of the PPR78 model (Predictive 1978, Peng-Robinson EOS with temperature dependent kij calculated through a group contribution method) to systems containing naphtenic compounds. Fluid Phase Equilib. 2006, 243, 9−28.

Figure 17. Saturate CO2 composition in Athabasca bitumen at different temperatures and pressures for different structural assumptions: Case 1 (solid lines), Case 2 (dashed lines), and Case 3 (dashed-dotted lines).

generalized group contribution pseudocomponent methodology for phase behavior modeling of mixtures of petroleum fluids and pure components. We validate our modeling using mixing of bitumen with one of three representative pure components: propane, CO2, and water. Our model predicts accurate results for propane and CO2 and results with reasonable and decreasing errors (as we move toward the critical point) for water because of its self-association characteristics. The results of our analysis show that the PPR78 model with a simple correction for water can predict reliable BIPs among crude PCs and pure components over a wide range of temperatures and pressures. Note that our model is potentially capable of accurately predicting the phase behaviors of mixtures of crude oil and other pure components, such as N2, H2, and H2S (further validations are needed for these components). This model can be used to study phase behaviors of petroleum fluids in reservoirs, oil recovery using supercritical CO2, crude upgrading using supercritical water, etc.



AUTHOR INFORMATION

Corresponding Author

*(P.H.) E-mail: [email protected]. Phone: +1 617 715 4348. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS The authors thank Saudi Aramco for supporting this work (Contract Number 6600023444). The support and encouragement of Dr. Ki-Hyouk Choi, Dr. Randall Field, Prof. William H. Green, and Prof. Michael T. Timko are greatly appreciated.



REFERENCES

(1) Orr, F. M.; Taber, J. J. Use of carbon-dioxide in enhanced oilrecovery. Science 1984, 224, 563−569. (2) Jaubert, J. N.; Avaullee, L.; Pierre, C. Is it still necessary to measure the minimum miscibility pressure? Ind. Eng. Chem. Res. 2002, 41, 303−310. (3) Berkowitz, N.; Calderon, J. Extraction of oil sand bitumens with supercritical water. Fuel Process. Technol. 1990, 25, 33−44. (4) Sato, T.; Adschiri, T.; Arai, K.; Rempel, G. L.; Ng, F. T. T. Upgrading of asphalt with and without partial oxidation in supercritical water. Fuel 2003, 82, 1231−1239. 8819

DOI: 10.1021/acs.iecr.5b02516 Ind. Eng. Chem. Res. 2015, 54, 8809−8820

Article

Industrial & Engineering Chemistry Research

hydrocarbons and associating compounds. Fluid Phase Equilib. 2011, 306, 112−123. (44) Xu, X. C.; Jaubert, J. N.; Privat, R.; Duchet-Suchaux, P.; BranaMulero, F. Predicting binary-interaction parameters of cubic equations of state for petroleum fluids containing pseudo-components. Ind. Eng. Chem. Res. 2015, 54, 2816−2824. (45) Hasan, M. U.; Ali, M. F.; Bukhari, A. Structural characterization of Saudi Arabian heavy crude-oil by NMR-spectroscopy. Fuel 1983, 62, 518−523. (46) Cyr, N.; Mcintyre, D. D.; Toth, G.; Strausz, O. P. Hydrocarbon structural group analysis of Athabasca asphaltene and its g.p.c. fractions by C-13 n.m.r. Fuel 1987, 66, 1709−1714. (47) Peng, D.; Robinson, D. B. A new 2-constant equation of state. Ind. Eng. Chem. Fundam. 1976, 15, 59−64. (48) Elliott, J. R.; Lira, C. T. Introductory Chemical Engineering Thermodynamics, 2nd ed.; Prentice Hall: New York, 2012. (49) Robinson, D. B.; Peng, D. Y. The Characterization of the Heptanes and Heavier Fractions for the GPA Peng-Robinson Programs; Gas Processors Association: Tulsa, 1978. (50) ExxonMobil. Crude Oil Information; July 2, 2015. Available from http://www.exxonmobil.com/crudeoil/about_crudes_kearl.aspx. (51) Petrov, A. A. Petroleum Hydrocarbons; Springer-Verlag: Berlin, 1987. (52) Speight, J. G. The Chemistry and Technology of Petroleum; Taylor & Francis: West Conshohocken, 1999. (53) Sabbagh, O.; Akbarzadeh, K.; Badamchi-Zadeh, A.; Svrcek, W. Y.; Yarranton, H. W. Applying the PR-EoS to asphaltene precipitation from n-alkane diluted heavy oils and bitumens. Energy Fuels 2006, 20, 625−634. (54) Badamchi-Zadeh, A.; Yarranton, H. W.; Svrcek, W. Y.; Maini, B. B. Phase behaviour and physical property measurements for VAPEX solvents: Part I. Propane and Athabasca bitumen. J. Can. Pet. Technol. 2009, 48, 54−61. (55) Badamchi-Zadeh, A.; Yarranton, H. W.; Maini, B. B.; Satyro, M. A. Phase behaviour and physical property measurements for VAPEX Solvents: Part II. Propane, carbon dioxide and Athabasca bitumen. J. Can. Pet. Technol. 2009, 48, 57−65. (56) Mehrotra, A. K.; Svrcek, W. Y. Bitumen density and gas solubility predictions using the Peng-Robinson equation of state. AOSTRA J. Res. 1985, 1, 15−229. (57) Amani, M. J.; Gray, M. R.; Shaw, J. M. Volume of mixing and solubility of water in Athabasca bitumen at high temperature and pressure. Fluid Phase Equilib. 2013, 358, 203−211.

(26) Vitu, S.; Privat, R.; Jaubert, J. N.; Mutelet, F. Predicting the phase equilibria of CO2+ hydrocarbon systems with the PPR78 model (PR EOS and k(ij) calculated through a group contribution method). J. Supercrit. Fluids 2008, 45, 1−26. (27) Privat, R.; Jaubert, J. N.; Mutelet, F. Addition of the nitrogen group to the PPR78 model (Predictive 1978, Peng Robinson EOS with temperature-dependent k(ij) calculated through a group contribution method). Ind. Eng. Chem. Res. 2008, 47, 2033−2048. (28) Privat, R.; Jaubert, J. N.; Mutelet, F. Addition of the sulfhydryl group (-SH) to the PPR78 model (predictive 1978, Peng-Robinson EOS with temperature dependent k(ij) calculated through a group contribution method). J. Chem. Thermodyn. 2008, 40, 1331−1341. (29) Privat, R.; Jaubert, J. N.; Mutelet, F. Use of the PPR78 model to predict new equilibrium data of binary systems involving hydrocarbons and nitrogen. Comparison with other GCEOS. Ind. Eng. Chem. Res. 2008, 47, 7483−7489. (30) Privat, R.; Mutelet, F.; Jaubert, J. N. Addition of the hydrogen sulfide group to the PPR78 Model (Predictive 1978, Peng-Robinson equation of state with temperature dependent k(ij) calculated through a group contribution method). Ind. Eng. Chem. Res. 2008, 47, 10041− 10052. (31) Privat, R.; Jaubert, J. N. Addition of the sulfhydryl group (-SH) to the PPR78 model: Estimation of missing group-interaction parameters for systems containing mercaptans and carbon dioxide or nitrogen or methane, from newly published data. Fluid Phase Equilib. 2012, 334, 197−203. (32) Qian, J. W.; Jaubert, J. N.; Privat, R. Phase equilibria in hydrogen-containing binary systems modeled with the Peng-Robinson equation of state and temperature-dependent binary interaction parameters calculated through a group-contribution method. J. Supercrit. Fluids 2013, 75, 58−71. (33) Qian, J. W.; Jaubert, J. N.; Privat, R. Prediction of the phase behavior of alkene-containing binary systems with the PPR78 model. Fluid Phase Equilib. 2013, 354, 212−235. (34) Qian, J. W.; Privat, R.; Jaubert, J. N.; Duchet-Suchaux, P. Enthalpy and heat capacity changes on mixing: Fundamental aspects and prediction by means of the PPR78 cubic equation of state. Energy Fuels 2013, 27, 7150−7178. (35) Qian, J.-W.; Privat, R.; Jaubert, J.-N. Predicting the phase equilibria, critical phenomena, and mixing enthalpies of binary aqueous systems containing alkanes, cycloalkanes, aromatics, alkenes, and gases (N2, CO2, H2S, H2) with the PPR78 Equation of State. Ind. Eng. Chem. Res. 2013, 52, 16457−16490. (36) Plée, V.; Jaubert, J. N.; Privat, R.; Arpentinier, P. Extension of the E-PPR78 equation of state to predict fluid phase equilibria of natural gases containing carbon monoxide, helium-4 and argon. J. Pet. Sci. Eng. 2015, DOI: 10.1016/j.petrol.2015.03.025. (37) Abdoul, W.; Rauzy, E.; Péneloux, A. Group-contribution equation of state for correlating and predicting thermodynamic properties of weakly polar and nonassociating mixtures - Binary and multicomponent systems. Fluid Phase Equilib. 1991, 68, 47−102. (38) Jaubert, J. N.; Privat, R.; Mutelet, F. Predicting the phase equilibria of synthetic petroleum fluids with the PPR78 approach. AIChE J. 2010, 56, 3225−3235. (39) Jaubert, J. N.; Privat, R. Relationship between the binary interaction parameters (k(ij)) of the Peng-Robinson and those of the Soave-Redlich-Kwong equations of state: Application to the definition of the PR2SRK model. Fluid Phase Equilib. 2010, 295, 26−37. (40) Hajiw, M.; Chapoy, A.; Coquelet, C. Hydrocarbons - water phase equilibria using the CPA Equation of state with a group contribution method. Can. J. Chem. Eng. 2015, 93, 432−442. (41) Gros, H. P.; Bottini, S.; Brignole, E. A. A group contribution equation of state for associating mixtures. Fluid Phase Equilib. 1996, 116, 537−544. (42) Ferreira, O.; Brignole, E. A.; Macedo, E. A. Modelling of phase equilibria for associating mixtures using an equation of state. J. Chem. Thermodyn. 2004, 36, 1105−1117. (43) Sánchez, F. A.; Pereda, S.; Brignole, E. A. GCA-EoS: A SAFT group contribution model-Extension to mixtures containing aromatic 8820

DOI: 10.1021/acs.iecr.5b02516 Ind. Eng. Chem. Res. 2015, 54, 8809−8820