- Email: [email protected]
Predicting phase behavior of nonpolar mixtures from refractive index at ambient conditions
Fuel 237 (2019) 637–647
Contents lists available at ScienceDirect
Fuel journal homepage: www.elsevier.com/locate/fuel
Full Length Article
Predicting phase behavior of nonpolar mixtures from refractive index at ambient conditions Caleb J. Siscoa,b, Mohammed I.L. Abutaqiyaa,b, Hassan Serrano Saadea,c, Francisco M. Vargasa,b,
T ⁎
a
Department of Chemical and Biomolecular Engineering, Rice University, Houston, TX 77005, USA ENNOVA LLC, Stafford, TX 77477, USA c Engineering School of Lorena, University of São Paulo, Estrada Municipal do Campinho, s/n - Pte. Nova, Lorena - SP, Brazil b
A R T I C LE I N FO
A B S T R A C T
Keywords: Phase behavior Equation of state Refractive index Density Optical properties Nonpolar mixtures Isothermal compressibility Thermal expansion
In a recent work, Evangelista and Vargas proposed correlations relating critical properties for nonpolar hydrocarbons to simple measurements of refractive index or density performed at ambient conditions (Evangelista and Vargas, Fluid Phase Equilibria, 2018). These correlations are used in this work to develop two approaches for improving phase behavior predictions with cubic equations of state. First, the correlations are applied to mixed solvents, treated as pseudo-pure fractions, whose refractive index and molecular weight were measured experimentally. Phase behavior predictions are performed on these systems to test the utility of a lumped solvent approach to be used when compositional information on a solvent is unavailable. Second, the correlations are applied within the default mixing rules of the cubic equation of state framework to calculate binary interaction parameters algebraically. This approach requires only simple laboratory measurements performed at ambient conditions instead of the more customary saturation pressure data that is more laborious to obtain and that might not be readily determined for certain substances.
1. Introduction
b = Ωb
The cubic equations of state (CEOS) are arguably the most popular class of thermodynamic models for describing the phase behavior of nonpolar hydrocarbon systems. Despite significant research efforts in the past few decades to improve the accuracy of thermodynamic models for hydrocarbon systems, CEOS are still widely used because of their unique combination of simplicity, sufficient accuracy for many applications, acceptable computational speed, and their broad availability across various industrial thermodynamics software packages. A generalized cubic equation of state is shown in Eq. (1):
P=
RT a − −b (V + δ1 b)(V + δ2 b ) V
(1)
a = Ωa
Pc
−
a=
α = [1 + κ (ω)[1− T Tc
(2)
]]2
∑ ∑ xi xj aij = ∑ ∑ xi xj i
−
α
(3)
(4)
The critical temperature (Tc ), critical pressure (Pc ) and acentric factor (ω ) are component-specific parameters while δ1, δ2 , Ωa , Ωb are constants specific to the EOS. The values of these constants, as well as κ (ω) , are shown in Table 1 for two CEOS: Soave-Redlich-Kwong (SRK), and Peng-Robinson (PR) equations of state. The SRK [2] and PR [3] equations of state are the most commonly used CEOS in the oil and gas industry. To extend the CEOS to model mixtures, mole-based averages are − − used to calculate a and b from component values of a and b . The van der Waals mixing rules are:
with the variables a , α and b calculated by:
R2Tc2
RTc Pc
b =
j
i
aiaj (1−kij )
j
∑ xi bi i
(5)
(6)
where x i is the mole fraction of component i and kij is the binary interaction parameter between components i and j . The binary interaction parameters are often tuned to match saturation pressure data.
⁎
Corresponding author at: Department of Chemical and Biomolecular Engineering, Rice University, Houston, TX 77005, USA. E-mail addresses: [email protected] (C.J. Sisco), [email protected] (M.I.L. Abutaqiya), [email protected] (H. Serrano Saade), [email protected] (F.M. Vargas). https://doi.org/10.1016/j.fuel.2018.09.112 Received 16 July 2018; Received in revised form 10 September 2018; Accepted 22 September 2018 Available online 13 October 2018 0016-2361/ © 2018 Elsevier Ltd. All rights reserved.
Fuel 237 (2019) 637–647
C.J. Sisco et al.
Table 1 Constants for selected cubic equations of state.
Table 2 Supplier, purity, and experimental data measured at 25 °C and 1 bar for studied hydrocarbons [5].
EOS
δ1
δ2
Ωa
Ωb
κ (ω)
SRK
1
0
0.42748
0.08664
0.480 + 1.574ω−0.176ω2
PR
1+
1− 2
0.45724
0.07780
0.37464 + 1.54226ω−0.26992ω2
2
With the critical properties (Tc , Pc , ω ) and binary interaction parameters (kij ) defined for all components present in the mixture, modeling with CEOS can be performed across a wide range of temperatures and pressures, and phase behavior predictions – including saturation pressure, density, thermal expansion, solubility parameter, etc. – can be generated. In a recent work, Evangelista and Vargas [1] proposed correlations relating critical properties for nonpolar hydrocarbons to simple laboratory measurements performed at ambient conditions. These correlations are given by Eqs. (7)–(10) and referred to as the EV correlations in this paper:
Component
Supplier
Purity mol%
nD –
ρ g/cm3
n-Hexane n-Heptane n-Decane n-Dodecane Benzene Toluene p-Xylene
Alfa Aesar Sigma-Aldrich Sigma-Aldrich Sigma-Aldrich OmniSolv Sigma-Aldrich Alfa Aesar
99 ≥ 99 ≥ 99 ≥ 99 ≥ 99.7 ≥ 99.9 99
1.3724 1.3852 1.4097 1.4195 1.4975 1.4938 1.4929
0.6554 0.6797 0.7266 0.7452 0.8736 0.8622 0.8567
Refractive Index measurement uncertainty: u(T) = 0.06 °C, u(nD) = 0.0002 with a 95% confidence interval. Density measurement uncertainty: u(T) = 0.02 °C, u(ρ) = 0.0004 g/cm3 with a 95% confidence interval.
study.
2 3 2 7 TC (K ) = 373.49FRI MW
(7)
2.2. Modeling methodology
8 5 −1 PC (bar) = 2.70 × 10 4FRI MW
(8)
−8 9 6 5 C (cm3 mol) = 0.47FRI V MW
(9)
The classical example for applying equations of state to pure components or mixtures is to use defined components with known compositions and critical properties reported in the literature [6]. However, many applications of EOS models are for complex systems – such as crude oils – for which a full, detailed compositional analysis is infeasible to obtain experimentally and the critical properties of some components are either not reported or might not be readily measurable. For these cases, researchers have developed characterization techniques with correlations to calculate critical properties for the otherwise undefined fractions of the system. One such example is the characterization technique of Pedersen et al. [7] that uses density and molecular weight correlations to specify the critical properties for the heavy carbon number fractions of a crude oil [8]. A similar modeling approach is proposed in this section, with critical properties for pseudo-pure fractions calculated from refractive index or density measurements. The critical properties for each pseudocomponent are calculated using the EV correlations with the refractive index and molecular weight measurements that are reported in Table 3. Critical properties for pure components used in this modeling study are shown in Table 4. To compare the approach proposed in this section (called the “lumped” approach because the multicomponent solvent is treated as a single pseudo-component) to the standard modeling approach, saturation curves are generated with PR EOS for both approaches. The systems modeled in this section are formed by adding a light key component (propane, hydrogen sulfide or carbon dioxide) in various proportions to the solvent fractions reported in Table 3 with 1-methylnaphthalene (1-MN) as the heavy key component. In the lumped approach, critical properties of the pseudo-pure solvents are calculated from the EV correlations, while the standard approach uses defined components with critical properties reported in the literature [6]. A graphical representation of the standard (Std) and lumped (EV) modeling approach is shown in Fig. 1.
ω = 8.62 ×
−1 10−4FRI MW
(10)
where MW has units of g/mol and FRI is calculated from refractive index or density measurements performed at 20 °C and 1 atm. The expression relating FRI to refractive index is given as:
FRI =
nD2 −1 nD2 + 2
(11)
where nD is the refractive index measured at the yellow sodium-D line.If refractive index is not reported, then FRI can be calculated from a mass density correlation proposed by Vargas and Chapman [4] given in Eq. (12).
FRI = 0.5054ρ−0.3951ρ2 + 0.2314ρ3 ≈
ρ 3
(12) 3
where ρ is mass density in units of g/cm . Applying the EV correlations to the CEOS framework offers two unique approaches to improving phase behavior predictions. First, critical properties for multicomponent nonpolar solvent mixtures can be calculated from simple volumetric or optical measurements without knowledge of the composition of the mixture. Second, binary interaction parameters can be calculated algebraically from density or refractive index measurements at ambient temperature rather than tuned to bubble pressure data and then subsequently used to improve phase behavior predictions for mixtures. These approaches are investigated in subsequent sections. 2. Modeling thermodynamic properties using critical property correlations applied to mixed solvents 2.1. Experimental methodology
2.3. Modeling results Experimental data for density and refractive index of nonpolar hydrocarbon mixtures were recently reported by Wang et al. [5] Density measurements were performed with the Anton Paar DMA 4500 digital vibrating U-tube densitometer and refractive index measurements were performed with the automatic Anton Paar WR refractometer under the sodium-D wavelength with a measuring range of 1.30–1.72. Table 2 provides a list of suppliers and chemical purity for the studied components. Wang et al. [5] also performed refractive index and density measurements on the mixtures reported in Table 3. These mixtures constitute the pseudo-pure fractions that are used in this modeling
Phase diagrams for the studied systems are shown in Fig. 2 with critical points represented by open squares. The molar amounts of these systems are given by the ratios associated with each curve. For example, 80:20:10 means that 80 mol of the mixture consists of the light key, 20 mol is the mixed solvent and 10 mol is the heavy key. For the lumped solvent approach, the molar amount of the solvent fraction is precisely the second number in the amount vector. In the standard de-lumped approach, that solvent amount is subdivided into defined components and molar amounts according to the compositions reported in Table 3. 638
Fuel 237 (2019) 637–647
C.J. Sisco et al.
Table 3 Properties of selected mixtures measured at 20 °C and 1 bar to be used as pseudo-components [5]. Mixtures
B1 B2 B3 T1 T2 T3
Amount (mol%)
MW
nD
ρ
Comp 1
Comp 2
Comp 3
z1
z2
z3
g/mol
–
g/cm3
n-Heptane n-Heptane Toluene n-Hexane n-Hexane Toluene
n-Dodecane Toluene p-Xylene n-Decane n-Dodecane Benzene
– – – n-Dodecane Toluene p-Xylene
63.0 47.9 53.6 33.5 33.3 33.3
37.0 52.1 46.4 33.3 33.4 33.7
– – – 33.2 33.3 33.0
126.18 96.00 98.64 132.83 116.24 92.04
1.4047 1.4345 1.4957 1.4080 1.4254 1.4968
0.716 0.764 0.864 0.722 0.751 0.867
Table 4 Simulation parameters for studied pure components.
α=
Component
MW g/mol
Tc K
Pc bar
ω –
C3 H2S CO2 n-C6 n-C7 n-C10 n-C12 Benzene Toluene p-Xylene 1-MN
44.096 34.081 44.010 86.175 100.202 142.285 170.339 78.114 92.141 106.168 142.200
369.85 373.10 304.13 507.82 540.13 617.70 658.20 562.16 591.79 616.20 772.00
42.48 90.00 73.77 30.18 27.27 21.20 18.24 48.98 41.04 35.11 36.00
0.1520 0.1005 0.2230 0.2979 0.3498 0.4890 0.5750 0.2080 0.2630 0.3180 0.3416
Std
1 ⎛ ∂V ⎞ 1 ∂V κ=− ⎛ ⎞ V ⎝ ∂T ⎠P V ⎝ ∂P ⎠T
(13)
The phase diagrams (Fig. 2) and parity plots (Figs. 3–5) show very good agreement between the proposed modeling approach (EV) and standard modeling approach (Std), demonstrating that the EV correlations provide representative solvent simulation parameters that describe nonpolar solvent phase behavior in an almost identical manner to the standard modeling approach. This is not, however, a comparison to experimental data. These plots convey only that the EV correlations applied within the CEOS framework describes phase behavior accurately with respect to the standard approach for obtaining compositions and simulation parameters for the CEOS. For mixtures that are not described accurately using the standard modeling approach (Std), the proposed modeling approach (EV) will likely show similar inaccuracies.
EV
3. Discussion Bar plots for cohesive energy and solubility parameter for the pure components studied in this work are shown in Figs. 6 and 7, respectively, showing that PR EOS reproduces the trend in cohesive energy much better than it does for solubility parameter. The reference data are obtained from a correlation for solubility parameter proposed by Wang [9].
δ = 52.042FRI + 2.904
(14)
With the definition of solubility parameter as the square root of the residual internal energy multiplied by density, Wang’s solubility parameter correlation can be rewritten in terms of residual internal energy.
−U res =
MW (52.042FRI + 2.904)2 ρ
(15)
Based on solubility parameter alone, one might expect that 1-MN would be always the appropriate choice as heavy key for the selected mixtures and that the solvent components are sufficiently different from 1-MN that the latter will be largely driving the phase-split at the dew point. However, the poor volumetric property predictions that are a well-known drawback of the CEOS cause the solubility parameter of 1MN to be under-predicted, as compared to Wang’s correlation, and the cohesive energy of n-C12 is greater than that of 1-MN at some conditions. Because the form of the equation of state and the mixing rules remain consistent, deviations between the lumped (EV) and de-lumped (Std) approaches can be explained by only two factors: the accuracy of the EV correlations to describe the solvent to which they are applied and the implicit assumption that all components in the lumped solvent partition equivalently. Fig. 8, which shows plots for cohesive energy and solubility parameter of the solvent mixtures, suggests that the EV correlations describe nonpolar solvent mixtures almost equivalently to the standard de-lumped modeling approach. However, regardless of the accuracy of the correlations, the lumped (EV) approach will always deviate from the de-lumped (Std) approach to some degree because of the limitation of the partitioning assumption. The magnitude of the deviations, though, might be negligible
Fig. 1. Standard (Std) and lumped (EV) modeling approaches. In the standard approach (Std), the solvent composition is obtained from Table 3 and simulation parameters are obtained from Table 4. In the lumped approach (EV), the solvent is treated as a single pseudo-component with simulation parameters calculated from the EV correlations given the experimental refractive index measurements in Table 3. The light key (LK) and heavy key (HK) are treated the same in each approach.
In addition to saturation pressure, various properties of the bulk solvent phase are of interest, and it is important that the proposed modeling approach does not introduce significant error to these predictions. Figs. 3–5 show parity plots comparing thermal expansion coefficient, isothermal compressibility, and density, respectively, of the saturated liquid phase calculated by the proposed approach (EV) to the standard modeling approach (Std). The diagonal line on the parity plots represents perfect agreement between the EV and standard modeling approaches. Points lying below the 45° line on the parity plots signify under-predictions by the EV approach with respect to the standard approach. The thermal expansion coefficient (α ) and isothermal compressibility (κ ) are calculated by: 639
Fuel 237 (2019) 637–647
C.J. Sisco et al.
P (bar)
250
250
C3 / B1 / 1-MN
200
200
150
150 80:20:10
100
100 250
P (bar)
200
300
500
700
80:20:10
P (bar)
500
700
80:20:10
100
300
500
700
C3 / B3 / 1-MN
150
50:50:10
100
80:20:10
100
50:50:10
50 20:80:10
0 100
P (bar)
20:80:10
0
200
50 300
500
20:80:10
0
700
100 250
C3 / T1 / 1-MN
200
150
300
500
700
CO2 / T2 / 1-MN 80:20:10
150
100
100
80:20:10 50:50:10
50 100
300
500
700
100 250 200
80:20:10
100
20:80:10
0
H2S / T3 / 1-MN
150
50:50:10
50
20:80:10
0
P (bar)
50:50:10
250
H2S / B3 / 1-MN
150
200
80:20:10
50
20:80:10
300
700
H2S / B2 / 1-MN
100
50:50:10
100
250
500
150
0
200
300
200
50
250
100 250
CO2 / B2 / 1-MN
100
200
20:80:10
0
150
250
50:50:10
50
20:80:10
0
80:20:10
100
50:50:10
50
H2S / B1 / 1-MN
300
500
CO2 / T3 / 1-MN
80:20:10
150
50:50:10
100
50:50:10
50
50 20:80:10
0 100
300
500
700
20:80:10
0
700
100
300
500
700
T(K)
T(K)
Fig. 2. Phase diagrams generated with PR EOS for the solvents listed in Table 3 with various light keys (LK) and 1-MN as the heavy key (HK). Solid lines: lumped solvent approach (EV). Dashed lines: standard de-lumped approach (Std). Molar amount vectors given by [LK:Solvent:HK].
cases, the equi-partitioning assumption causes very little deviation to property predictions. However, if n-butane were added to the solvent, the bubble pressure calculation for the lumped and de-lumped approaches would not agree as closely. Propane and n-butane have similar volatilities, and the solubility of n-butane in the solvent has a nonnegligible effect on bubble pressure that would be obfuscated by lumping n-butane with the heavy components. The same argument applies for the dew pressure calculation. The more dissimilar the heavy key is from the components that constitute the solvent, the lesser impact lumping the solvent will have on dew pressure predictions. Though all systems studied show very good saturation pressure
depending on the system and how the solvent is lumped. In general, saturation pressure predictions for the lumped and delumped approaches agree closely when the component(s) most responsible for driving the phase-split are well characterized and their partitioning with respect to the solvent are only negligibly affected by lumping. For example, the calculated bubble pressure of a mixture containing propane with toluene, p-xylene and 1-MN will be essentially the same with the lumped and de-lumped approaches because the solubility of propane in the solvent is only negligibly affected by lumping the heavy components and only a small amount of the heavy components will partition into the vapor phase at the bubble point. In these 640
Fuel 237 (2019) 637–647
C.J. Sisco et al. 6
6
Fig. 4. Isothermal compressibility (κ) parity plots comparing the lumped solvent approach (EV) on the x-axis to the standard de-lumped modeling approach (Std) on the y-axis with PR EOS. The diagonal line represents perfect agreement between the EV and Std approach.
4. Calculating binary interaction parameters from critical property correlations
HS
CΎ
αStd (10-3 K-1)
αStd (10-3 K-1)
Fig. 3. Thermal expansion coefficient (α) parity plots comparing the lumped solvent approach (EV) on the x-axis to the standard de-lumped modeling approach (Std) on the y-axis with PR EOS. The diagonal line represents perfect agreement between the EV and Std approach.
process accurately in the dew region, it is important that there is sufficient resolution in the heavy components that are driving the phasesplit. When one of the solvent components that is important to this process is lumped in with other solvent components, then its effect is diluted and the prediction suffers. When the difference in the cohesive energy between the heavy key and solvent components is sufficiently large, which is the case for the solvents that do not contain n-C12 (B2, B3, T3), then the deviations between the lumped and standard approaches are negligible. The modeling study performed in this work is a proof-of-concept for the EV correlations and the lumped solvent approach. With the accuracy of the approach established, it is important to understand how and when the proposed approach can be applied to practical systems. The approach is primarily of interest for defining the simulation parameters of unknown fractions of a complex mixture. If all components and their compositions in a mixture are known in advance, then the standard approach would be always preferable unless the mixture exceeds some limitation on the number of components allowed in a simulation. In these cases, the EV correlations offer a method to lump similar components into a solvent fraction and calculate the critical properties of the pseudo-component. However, the primary reason these correlations were developed was to describe fractions of complex crude oils for which a detailed compositional analysis is not readily available. The crude oil can be lumped into as many or few fractions as desired and then the critical properties of each fraction can be obtained from simple refractive index or density measurements performed on each fraction at ambient conditions. The number of component fractions and whether the light end and/or heavy end requires more resolution depends on which properties are being predicted and the required accuracy of the predictions. If only bubble pressures are desired, then a simple lumping approach with all components heavier than n-hexane lumped into a single fraction is often sufficient, as shown by Abutaqiya et al. [10]. If dew pressures or asphaltene onsets are desired, then the heavy end must be characterized in more detail, as shown in the more customary crude oil characterization methodologies [7,11–14]. It was expected that solubility parameter would offer an indication as to why certain solvent mixtures showed consistently worse predictions than others, but there was no apparent trend in solubility parameter that suggested when the lumped solvent approach predictions might suffer. Cohesive energy, on the other hand, did show such a trend. When the difference in cohesive energy between the light and/or heavy key and all solvent components is large, then the lumped solvent approach will introduce only negligible error to saturation pressure predictions. The Wang correlation for residual internal energy in Eq. (15) (which is equivalent to cohesive energy) offers a means to check on whether de-lumping the crude oil fractions might be appropriate. Eq. (15) requires only refractive index (or density) and molecular weight measurements, and PR EOS with EV correlations was shown to match the residual internal energy correlation quite well, suggesting that cohesive energy can be used as a heuristic for determining the number of fractions into which the crude oil should be separated in order to generate simulation parameters that are representative of the actual mixture.
CO
4
2 B1 0 2 4 αEV (10-3 K-1)
CΎ
4 2 B2
CO
HS
6
0 6
CΎ
αStd (10-3 K-1)
6
αStd (10-3 K-1)
HS
0 0
4
2 B3 0
2 4 αEV (10-3 K-1) CO
HS
6
CΎ
4
2 T1 0
6
2 4 αEV (10-3 K-1) CO
HS
6
0 6
CΎ
αStd (10-3 K-1)
0
αStd (10-3 K-1)
CO
4
2 T2 0
2 4 αEV (10-3 K-1) CO
HS
6
CΎ
4
2 T3 0
0
2 4 αEV (10-3 K-1)
6
0
2 4 αEV (10-3 K-1)
6
4.1. Modeling methodology predictions – due mostly to the accuracy of the EV correlations – some solvents (B1, T1, T2) showed much more deviation between the lumped and standard approaches, particularly in the dew region, than the other solvents. The primary reason this deviation occurs is because the heavy key (1-MN) and the heaviest component in the solvent for these cases (n-C12) have similar cohesive energy. To describe the phase-splitting
The standard procedure for obtaining binary interaction parameters is by tuning to saturation pressure data for a given component pair. Other methodologies have been developed that use critical volumes [15] or a group contribution approach [16–20] to calculate binary interaction parameters, thereby eliminating the need for tuning. The approach presented in this work uses the EV correlations and the CEOS 641
Fuel 237 (2019) 637–647
C.J. Sisco et al.
mixing rule directly to calculate binary interaction parameters. Only a single data point for the binary mixture at any composition is required to implement this approach, but more data at different compositions can be used to refine the calculations. The necessary data for this approach are very easy to obtain, as the EV correlations require only simple measurements of a solvent at ambient conditions to define its critical properties. − Writing the expression for a from the van der Waals mixing rule for a binary mixture yields: −
a = x12 a1 + x 22 a2 + 2x1 x2 a1a2 (1−k12)
(16)
This is the standard mixing rule for CEOS, but the approach proposed here can be extended to any mixing rule. The optimal value of k12 is specific to a given component pair and the mixing rule. For binary mixtures of similar components with the van der Waals mixing rules, the optimal value for k12 is usually near zero while stronger deviations from zero are observed for dissimilar components. Large deviations from a k12 value of zero is generally interpreted as a deficiency in the mixing rule to describe the interaction energy of components 1 and 2. − If the process is reversed and a is already known for a given composition, then k12 can be calculated explicitly by rearranging the mixing rule to obtain: −
k12 = 1−(a−x12 a1−x 22 a2) (2x1 x2 a1a2 )
In this section, a methodology is proposed for calculating a from refractive index measurements and the EV correlations and then solving explicitly for k12 . To illustrate, consider an equimolar mixture of components 1 and 2. FRI can be calculated for this mixture either from refractive index or density measurements at 20 °C and 1 atm, and the molecular weight of the mixture can be measured experimentally or calculated from the compositions and pure-component molecular weights. With FRI and MW known, the critical properties of the mixture − (Tc, mix , Pc, mix , ωmix ) are calculated from the EV correlations and a is calculated the same as in Eq. (2), treating the binary mixture as a − − pseudo-component: a = Ωa (R2Tc, mix 2αmix ) Pc, mix . Because a is calculated for a specific mixture composition ( x1 = 0.5 for this example) and a1 and a2 are known from the critical properties of the pure components (either obtained from the literature or from the EV correlations), all terms on the right-hand side of Eq. (17) are known and k12 can be calculated explicitly. Additionally – because of the temperature-dependence implicit in the α term – k12 can be calculated as a function of temperature without any need for further experiments. There are no restrictions on the compositions of the binary mixtures or on the number of compositions for which the experiment must be performed. In principle, the calculation of k12 requires only one experimental measurement for refractive index or density at a known mixture composition of components 1 and 2, as illustrated in this example, but a more representative k12 value would likely be obtained if measurements for multiple compositions of components 1 and 2 were considered. To evaluate fairly the proposed methodology for calculating k12 , the γ −φ formulation is preferred to the more standard φ−φ formulation because the latter approach can sometimes show poor predictions for pure-component saturation pressures. The binary interaction parameter cannot correct pure-component properties, so showing poor predictions for these properties might mislead the reader about the usefulness of the proposed methodology. The γ −φ formulation for the equi-fugacity criterion is written as:
Cohesive Energy, -Ures/RT
30 25 20 15 10 5 0
Fig. 6. Cohesive energy for selected pure components at 20 °C and 1 bar from Eq. (15) [blue] and PR EOS [red]. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 22 21 20 19 Solubility Parameter,
(17) −
Fig. 5. Density (ρ) parity plots comparing the lumped solvent approach (EV) on the x-axis to the standard de-lumped modeling approach (Std) on the y-axis with PR EOS. The diagonal line represents perfect agreement between the EV and Std approach.
18 17 16 15 14 13 12
φ sat yi φiV̂ P = x i φiL̂ ⎛⎜ i ⎟⎞ Pisat ⎝ φi ⎠
Fig. 7. Solubility parameter for selected pure components at 20 °C and 1 bar from Eq. (14) [blue] and PR EOS [red]. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
φiV̂
(18)
and φi ̂ are the fugacity coefficients of the vapor and liquid, where respectively, calculated from the vapor composition {yi } and liquid composition {x i} at the pressure P and temperature T of interest. φisat is the fugacity coefficient of pure component i evaluated at Pisat and φi is 642
L
Fuel 237 (2019) 637–647
C.J. Sisco et al. 22 21 25
Solubility Parameter,
Cohesive Energy, -Ures/RT
30
20 15 10 5
20
19 18 17 16 15 14 13
0
12 B1
B2
B3
T1
T2
T3
B1
B2
B3
T1
T2
T3
Fig. 8. Cohesive energy and solubility parameter plots for solvent mixtures at 20 °C and 1 bar from PR EOS with standard (Std) approach [blue] and lumped solvent (EV) approach [red]. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
The calculation of k12 was completed via two approaches. In the EV1 approach, both the pure component critical properties and mixture properties were calculated with the EV correlations, while in the EV2 approach, the pure component critical properties were taken from the literature and mixture properties were calculated with the EV correlations.
Table 5 Density or refractive index values at 20 °C and 1 bar for selected binary mixtures. Mixtures
Amount (mol%)
MW
nD
ρ
Comp 1
Comp 2
z1
z2
g/mol
-
g/cm3
Benzene Toluene Cyclohexane n-Heptane Cyclohexane Cyclohexane
Toluene n-Decane n-Decane Toluene n-Dodecane n-Heptane
50.0 50.0 50.0 47.9 50.0 50.0
50.0 50.0 50.0 52.1 50.0 50.0
85.13 117.21 113.22 96.17 127.25 92.18
– – – 1.4345 – –
0.8721 0.7780 0.7471 – 0.7567 0.7225
6. Discussion In the EV1 approach, the critical properties for both the binary mixture and constituent pure components are calculated from the EV correlations. When the EV correlations describe mixture properties poorly, it is likely that it also describes one or both constituent pure component properties poorly. Because k12 is quantifying a difference between the mixture and pure component properties, it is not necessary that the EV correlations describe the mixture and pure component properties well for the EV1 approach to be effective. Rather, it is only important that the difference between the mixture and pure component properties is captured accurately. The EV1 approach benefits from a cancellation of errors when both the pure components and mixture are described poorly. In the EV2 approach, the critical properties of the binary mixture are calculated from the EV correlations, whereas the constituent pure component critical properties are taken from the literature. Because the literature values will necessarily produce good predictions for saturation pressure for the pure components, the EV correlations must also
the fugacity coefficient of pure component i evaluated at P . All fugacity coefficients are calculated using the PR EOS, and saturation pressure of pure component i (Pisat ) is calculated with correlations provided in the DIPPR database [6] and is a function of temperature only.
5. Modeling results P-xy curves for binary mixtures with k12 calculated from the mixing rule and k12 = 0 are compared to reference VLE data generated from NIST REFPROP software [21]. The binary mixtures studied are listed in Table 5. The results are presented in Figs. 9–14 and the mean absolute percent deviation is presented in Table 6. 0.16
0.40
298.15 K
P (bar)
P (bar)
323.15 K
0.30
0.12 0.08
0.20 0.10
0.04
0.00
0.00
0
0.2
0.4
x1,y1
0.6
0.8
1
0
1.00
0.4
x1,y1
0.6
0.8
1
0.6
0.8
1
2.00 348.15 K
373.15 K
0.80
1.60
0.60
1.20
P (bar)
P (bar)
0.2
0.40 0.20
0.80 0.40
0.00
0.00 0
0.2
0.4
x1,y1
0.6
0.8
1
0
0.2
0.4
x1,y1
Fig. 9. P-xy diagram for Benzene (1)/Toluene (2) with PR EOS. Open circles: NIST REFPROP [21]. Black dashed lines: k12 = 0 . Blue lines: EV1. Red lines: EV2. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 643
Fuel 237 (2019) 637–647
C.J. Sisco et al. 0.04
0.16
298.15 K
0.02
0.01
0.08 0.04
0.00
0.00 0
0.40
0.2
0.4
x1,y1
0.6
0.8
1
0
0.2
0.4
x1,y1
0.6
0.8
1
0.6
0.8
1
0.80
348.15 K
373.15 K 0.60
P (bar)
0.30
P (bar)
323.15 K
0.12
P (bar)
P (bar)
0.03
0.20
0.40
0.20
0.10
0.00
0.00
0
0.2
0.4
x1,y1
0.6
0.8
0
1
0.2
0.4
x1,y1
Fig. 10. P-xy diagram for Toluene (1)/n-Decane (2) with PR EOS. Open circles: NIST REFPROP [21]. Black dashed lines: k12 = 0 . Blue lines: EV1. Red lines: EV2. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 0.16
0.40
298.15 K
323.15 K 0.30
P (bar)
P (bar)
0.12
0.08
0.10
0.04
0.00
0.00 0 1.00
0.2
0.4
x1,y1
0.6
0.8
0
1
0.2
0.4
x1,y1
0.6
0.8
1
0.6
0.8
1
2.00
348.15 K
373.15 K
0.80
1.60
0.60
1.20
P (bar)
P (bar)
0.20
0.40 0.20
0.80 0.40
0.00
0.00 0
0.2
0.4
x1,y1
0.6
0.8
1
0
0.2
0.4
x1,y1
Fig. 11. P-xy diagram for Cyclohexane (1)/n-Decane (2) with PR EOS. Open circles: NIST REFPROP [21]. Black dashed lines: k12 = 0 . Blue lines: EV1. Red lines: EV2. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
approach performs poorly, such as the cyclohexane / n-dodecane mixture, one component yields much better saturation pressure predictions with the EV correlations than the other such that the errorcanceling effect is only partially realized. In cases such as benzene/ toluene, in which the EV correlations under-predict saturation pressures for both components, the deviations are roughly equal in magnitude and direction so that the error-canceling effect is more fully realized. Fortunately, the EV correlations are amenable for re-tuning to specific families of components or even a specific binary pair. All EV correlations have the same standard form given by:
produce good predictions for the mixture saturation pressure if the EV2 approach is to generate suitable values for k12 . When the EV2 approach fails to generate satisfactory values for k12 , it is due to the EV correlations describing one or both constituent pure components poorly and the deviations occurring in the same direction such that the errors are additive instead of negating. For example, saturation pressures for benzene and toluene are both under-predicted by the EV correlations, so the mixture saturation pressure is also under-predicted. Because the correlations-only approach (EV1) cancels the error in mixture saturation pressure with correspondingly poor pure component saturation pressure predictions, the EV1 approach outperforms the EV2 approach in these cases. In cases such as the toluene/n-decane mixture, the EV correlations over-predict saturation pressure for n-decane by roughly the same order of magnitude that they under-predict the saturation pressure for toluene. This yields a fortunate cancellation of errors for the mixture saturation pressure and correspondingly satisfactory results for k12 with the EV2 approach. To maintain consistency, it is best to use the EV1 approach always because of the error-canceling effect. The cases for which the EV1
Y = aFRI bMW c
(19)
where Y is the property of interest and a , b , and c are tuning constants. The tuning constants of the EV correlations can be fit to reproduce the properties for a specific set of compounds to increase the accuracy of certain binary systems. The same EV1 approach proposed here but with new constants can then be used to improve the results for the k12 calculation.
644
Fuel 237 (2019) 637–647
C.J. Sisco et al. 0.20
0.07
323.15 K
298.15 K 0.18 0.16
P (bar)
P (bar)
0.06
0.05
0.14
0.04
0.12
0.03
0.10 0
0.2
0.4
0.6
0.8
0
1
0.2
0.4
x1,y1 0.50
1.10
0.8
1
0.6
0.8
1
373.15 K
348.15 K 0.46
1.00
0.42
P (bar)
P (bar)
0.6
x1,y1
0.38
0.90 0.80
0.34
0.70
0.30 0
0.2
0.4
0.6
0.8
0
1
0.2
0.4 x1,y1
x1,y1
Fig. 12. P-xy diagram for n-Heptane (1)/Toluene (2) with PR EOS. Open circles: NIST REFPROP [21]. Black dashed lines: k12 = 0 . Blue lines: EV1. Red lines: EV2. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 0.16
0.40
298.15 K
323.15 K
0.30
P (bar)
P (bar)
0.12 0.08
0.04
0.20 0.10
0.00
0.00 0
0.2
0.4
x1,y1
0.6
0.8
1
0
1.00
2.00
0.60
1.20
P (bar)
P (bar)
1.60
0.40 0.20
0.4
x1,y1
0.6
0.8
1
0.6
0.8
1
373.15 K
348.15 K 0.80
0.2
0.80 0.40
0.00
0.00
0
0.2
0.4
x1,y1
0.6
0.8
1
0
0.2
0.4
x1,y1
Fig. 13. P-xy diagram for Cyclohexane (1)/n-Dodecane (2) with PR EOS. Open circles: NIST REFPROP [21]. Black dashed lines: k12 = 0 . Blue lines: EV1. Red lines: EV2. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
7. Conclusions
the lumped solvent modeling approach produces very good predictions for saturation pressure in the dew region. When the heavy key and solvent components have similar cohesive energy, dew pressure predictions suffer if the solvent is lumped. The EV correlations were also applied to binary mixtures to calculate binary interaction parameters algebraically from refractive index or density measurements at ambient conditions instead of the standard method of tuning to saturation pressure data, which are more laborious to obtain. Because the calculation of k12 requires that components 1 and 2 are known in advance, the accuracy of the EV correlations in describing the constituent pure components can also be known. The utility of this approach for generating good k12 values for a given component pair can then be estimated in advance. However, even if the EV correlations with constants given in this work do not describe a component pair with sufficient precision, it is relatively straightforward to tune the correlations to match a specific set of components and then generate k12 values using the same approach proposed in this work. The EV correlations were originally tuned to many nonpolar hydrocarbons of different families (n-alkanes, cyclics, and aromatics), which means that
In this work, critical property correlations recently proposed by Evangelista and Vargas [1] were evaluated for their accuracy in describing nonpolar solvent mixtures. The correlations, along with the proposed lumped solvent approach, show very good predictive capability for saturation pressures and bulk liquid properties like density, thermal expansion, and isothermal compressibility. This suggests that the approach is suitable to be applied to nonpolar crude oil fractions for which the composition is unknown. The EV correlations offer a simple expression for relating easily measured volumetric or optical properties to critical properties, and the lumped solvent approach shows small deviations with respect to the standard modeling approach. Additionally, it was shown that cohesive energy calculated by PR EOS with EV correlations matches well with the Wang correlation for cohesive energy, and differences in cohesive energy between the heavy key and the components comprising the solvent was shown to correlate with accuracy of the lumped solvent approach in matching dew point predictions. When the heavy key and solvent are sufficiently different, 645
Fuel 237 (2019) 637–647
C.J. Sisco et al. 0.40
0.15
323.15 K
0.13
0.35
0.11
0.30
P (bar)
P (bar)
298.15 K
0.09
0.25 0.20
0.07
0.15
0.05 0
0.2
0.4
0.6
0.8
0
1
0.2
0.4
0.90
1.80
348.15 K
0.80
0.8
1
0.6
0.8
1
373.15 K
1.60
0.70
P (bar)
P (bar)
0.6 x1,y1
x1,y1
0.60
1.40
1.20
0.50 0.40
1.00 0
0.2
0.4
0.6
0.8
0
1
0.2
x1,y1
0.4 x1,y1
Fig. 14. P-xy diagram for Cyclohexane (1)/n-Heptane (2) with PR EOS. Open circles: NIST REFPROP [21]. Black dashed lines: k12 = 0 . Blue lines: EV1. Red lines: EV2. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Table 6 Binary interaction parameters (kij ) and error analysis for selected binary mixtures. T / °C
Benzene + Toluene
Toluene + n-Decane
Cyclohexane + n-Decane
n−Heptane + Toluene
Cyclohexane + n-Dodecane
Cyclohexane + n-Heptane
25 50 75 100 25 50 75 100 25 50 75 100 25 50 75 100 25 50 75 100 25 50 75 100
AAPD in BP
kij
AAPD in DP
EV1
EV2
EV1
EV2
kij = 0
EV1
EV2
kij = 0
0.0001 0.0005 0.0004 0.0004 0.0119 0.0111 0.0103 0.0095 0.0166 0.0156 0.0146 0.0136 −0.0038 −0.0039 −0.0041 −0.0043 0.0335 0.0320 0.0305 0.0289 0.0003 0.0003 0.0003 0.0003
−0.1037 −0.1032 −0.1027 −0.1023 −0.0117 −0.0118 −0.0118 −0.0118 −0.0335 −0.0301 −0.0267 −0.0233 −0.0364 −0.0351 −0.0337 −0.0324 −0.0181 −0.0155 −0.0130 −0.0104 −0.0622 −0.0561 −0.0500 −0.0438
0.577 0.361 0.252 0.187 11.442 10.336 9.786 9.404 18.361 15.732 13.735 12.153 4.632 4.012 3.381 2.772 46.148 40.729 36.078 31.979 1.181 1.116 1.300 1.566
25.178 22.552 20.241 18.221 3.186 1.791 1.043 1.765 10.747 7.192 4.106 1.547 13.412 11.473 9.738 8.200 3.225 6.350 8.913 10.795 15.120 12.014 9.214 6.759
0.413 0.237 0.155 0.118 4.003 4.545 5.270 5.886 7.284 7.034 6.920 6.827 3.476 2.958 2.413 1.877 15.826 16.122 16.142 15.892 1.092 1.037 1.229 1.501
0.504 0.212 0.157 0.110 0.884 1.362 1.983 2.677 0.771 1.152 1.588 2.033 3.623 3.329 2.913 2.449 0.179 0.826 1.390 1.986 0.723 0.784 1.009 1.313
25.273 23.031 20.931 19.002 1.676 1.531 1.042 0.657 1.481 1.563 1.425 1.055 12.572 10.873 9.305 7.885 0.060 0.541 0.911 1.286 15.001 11.932 9.168 6.747
0.367 0.128 0.142 0.160 0.659 0.721 1.034 1.612 0.542 0.710 0.940 1.215 2.601 2.381 2.030 1.623 0.106 0.651 1.076 1.493 0.663 0.722 0.948 1.256
some components will not be as well-described by the correlations as others. Because all EV correlations have the same standard form, a systematic approach for tuning the constants to capture more precisely the properties of certain component sets can be easily implemented.
[6] [7]
References
[8] [9]
[1] Evangelista RF, Vargas FM. Prediction of the temperature dependence of densities and vapor pressures of nonpolar hydrocarbons based on their molecular structure and refractive index data at 20 °C. Fluid Phase Equilib 2018;468:29–37. [2] Soave G. Equilibrium constants from a modified Redlich-Kwong equation of state. Chem Eng Sci 1972;27(6):1197–203. [3] Peng D-Y, Robinson DB. A new two-constant equation of state. Ind Eng Chem Fundam 1976;15(1):59–64. [4] Vargas FM, Chapman WG. Application of the One-Third rule in hydrocarbon and crude oil systems. Fluid Phase Equilib 2010;290(1–2):103–8. [5] Wang F, Evangelista RF, Threatt TJ, Tavakkoli M, Vargas FM. Determination of
[10]
[11]
[12]
646
volumetric properties using refractive index measurements for nonpolar hydrocarbons and crude oils. Ind Eng Chem Res 2017. Wilding WV, Rowley RL, Oscarson JL. DIPPR® Project 801 evaluated process design data. Fluid Phase Equilib 1998;150:413–20. Pedersen KS, Blilie AL, Meisingset KK. PVT calculations on petroleum reservoir fluids using measured and estimated compositional data for the plus fraction. Ind Eng Chem Res 1992;31(5):1378–84. Panuganti SR, Vargas FM, Chapman WG. Property scaling relations for nonpolar hydrocarbons. Ind Eng Chem Res 2013;52(23):8009–20. Wang J. Predicting Asphaltene Flocculation in Crude Oils. Socorro, New Mexico: New Mexico Institute of Mining & Technology; 2000. Abutaqiya MIL, Panuganti SR, Vargas FM. Efficient algorithm for the prediction of pressure–volume–temperature properties of crude oils using the perturbed-chain statistical associating fluid theory equation of state. Ind Eng Chem Res May 2017;56(20):6088–102. Panuganti SR, Vargas FM, Gonzalez DL, Kurup AS, Chapman WG. PC-SAFT characterization of crude oils and modeling of asphaltene phase behavior. Fuel 2012;93:658–69. Tavakkoli M, Chen A, Vargas FM. Rethinking the modeling approach for asphaltene
Fuel 237 (2019) 637–647
C.J. Sisco et al.
[13]
[14]
[15]
[16]
[17]
precipitation using the PC-SAFT Equation of State. Fluid Phase Equilib 2016;416:120–9. Arya A, von Solms N, Kontogeorgis GM. Determination of asphaltene onset conditions using the cubic plus association equation of state. Fluid Phase Equilib 2015;400:8–19. Arya A, von Solms N, Kontogeorgis GM. Investigation of the gas injection effect on asphaltene onset precipitation using the cubic-plus-association equation of state. Energy Fuels 2016;30(5):3560–74. Gao G, Daridon J-L, Saint-Guirons 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. Jaubert J-N, Mutelet F. VLE predictions with the Peng-Robinson equation of state and temperature dependent kij calculated through a group contribution method. Fluid Phase Equilib 2004;224(2):285–304. Jaubert J-N, Vitu S, Mutelet F, Corriou J-P. Extension of the PPR78 model
[18]
[19]
[20]
[21]
647
(predictive 1978, Peng-Robinson EOS with temperature dependent kij calculated through a group contribution method) to systems containing aromatic compounds. Fluid Phase Equilib 2005;237(1–2):193–211. Soave G, Gamba S, Pellegrini LA. SRK equation of state: Predicting binary interaction parameters of hydrocarbons and related compounds. Fluid Phase Equilib 2010;299(2):285–93. Abudour AM, Mohammad SA, Robinson Jr. RL, Gasem KAM. Generalized binary interaction parameters for the Peng-Robinson equation of state. Fluid Phase Equilib 2014;383:156–73. Nasrifar K, Rahmanian N. Equations of state with group contribution binary interaction parameters for calculation of two-phase envelopes for synthetic and real natural gas mixtures with heavy fractions. Oil Gas Sci. Technol. 2018;73:7. Lemmon EW, Huber ML, McLinden MO. “NIST Standard Reference Database 23: Reference Fluid Thermodynamic and Transport Properties-REFPROP, Version 9.1.” In: Natl Std Ref Data Ser. NIST NSRDS -, May 2013.
Contents lists available at ScienceDirect
Fuel journal homepage: www.elsevier.com/locate/fuel
Full Length Article
Predicting phase behavior of nonpolar mixtures from refractive index at ambient conditions Caleb J. Siscoa,b, Mohammed I.L. Abutaqiyaa,b, Hassan Serrano Saadea,c, Francisco M. Vargasa,b,
T ⁎
a
Department of Chemical and Biomolecular Engineering, Rice University, Houston, TX 77005, USA ENNOVA LLC, Stafford, TX 77477, USA c Engineering School of Lorena, University of São Paulo, Estrada Municipal do Campinho, s/n - Pte. Nova, Lorena - SP, Brazil b
A R T I C LE I N FO
A B S T R A C T
Keywords: Phase behavior Equation of state Refractive index Density Optical properties Nonpolar mixtures Isothermal compressibility Thermal expansion
In a recent work, Evangelista and Vargas proposed correlations relating critical properties for nonpolar hydrocarbons to simple measurements of refractive index or density performed at ambient conditions (Evangelista and Vargas, Fluid Phase Equilibria, 2018). These correlations are used in this work to develop two approaches for improving phase behavior predictions with cubic equations of state. First, the correlations are applied to mixed solvents, treated as pseudo-pure fractions, whose refractive index and molecular weight were measured experimentally. Phase behavior predictions are performed on these systems to test the utility of a lumped solvent approach to be used when compositional information on a solvent is unavailable. Second, the correlations are applied within the default mixing rules of the cubic equation of state framework to calculate binary interaction parameters algebraically. This approach requires only simple laboratory measurements performed at ambient conditions instead of the more customary saturation pressure data that is more laborious to obtain and that might not be readily determined for certain substances.
1. Introduction
b = Ωb
The cubic equations of state (CEOS) are arguably the most popular class of thermodynamic models for describing the phase behavior of nonpolar hydrocarbon systems. Despite significant research efforts in the past few decades to improve the accuracy of thermodynamic models for hydrocarbon systems, CEOS are still widely used because of their unique combination of simplicity, sufficient accuracy for many applications, acceptable computational speed, and their broad availability across various industrial thermodynamics software packages. A generalized cubic equation of state is shown in Eq. (1):
P=
RT a − −b (V + δ1 b)(V + δ2 b ) V
(1)
a = Ωa
Pc
−
a=
α = [1 + κ (ω)[1− T Tc
(2)
]]2
∑ ∑ xi xj aij = ∑ ∑ xi xj i
−
α
(3)
(4)
The critical temperature (Tc ), critical pressure (Pc ) and acentric factor (ω ) are component-specific parameters while δ1, δ2 , Ωa , Ωb are constants specific to the EOS. The values of these constants, as well as κ (ω) , are shown in Table 1 for two CEOS: Soave-Redlich-Kwong (SRK), and Peng-Robinson (PR) equations of state. The SRK [2] and PR [3] equations of state are the most commonly used CEOS in the oil and gas industry. To extend the CEOS to model mixtures, mole-based averages are − − used to calculate a and b from component values of a and b . The van der Waals mixing rules are:
with the variables a , α and b calculated by:
R2Tc2
RTc Pc
b =
j
i
aiaj (1−kij )
j
∑ xi bi i
(5)
(6)
where x i is the mole fraction of component i and kij is the binary interaction parameter between components i and j . The binary interaction parameters are often tuned to match saturation pressure data.
⁎
Corresponding author at: Department of Chemical and Biomolecular Engineering, Rice University, Houston, TX 77005, USA. E-mail addresses: [email protected] (C.J. Sisco), [email protected] (M.I.L. Abutaqiya), [email protected] (H. Serrano Saade), [email protected] (F.M. Vargas). https://doi.org/10.1016/j.fuel.2018.09.112 Received 16 July 2018; Received in revised form 10 September 2018; Accepted 22 September 2018 Available online 13 October 2018 0016-2361/ © 2018 Elsevier Ltd. All rights reserved.
Fuel 237 (2019) 637–647
C.J. Sisco et al.
Table 1 Constants for selected cubic equations of state.
Table 2 Supplier, purity, and experimental data measured at 25 °C and 1 bar for studied hydrocarbons [5].
EOS
δ1
δ2
Ωa
Ωb
κ (ω)
SRK
1
0
0.42748
0.08664
0.480 + 1.574ω−0.176ω2
PR
1+
1− 2
0.45724
0.07780
0.37464 + 1.54226ω−0.26992ω2
2
With the critical properties (Tc , Pc , ω ) and binary interaction parameters (kij ) defined for all components present in the mixture, modeling with CEOS can be performed across a wide range of temperatures and pressures, and phase behavior predictions – including saturation pressure, density, thermal expansion, solubility parameter, etc. – can be generated. In a recent work, Evangelista and Vargas [1] proposed correlations relating critical properties for nonpolar hydrocarbons to simple laboratory measurements performed at ambient conditions. These correlations are given by Eqs. (7)–(10) and referred to as the EV correlations in this paper:
Component
Supplier
Purity mol%
nD –
ρ g/cm3
n-Hexane n-Heptane n-Decane n-Dodecane Benzene Toluene p-Xylene
Alfa Aesar Sigma-Aldrich Sigma-Aldrich Sigma-Aldrich OmniSolv Sigma-Aldrich Alfa Aesar
99 ≥ 99 ≥ 99 ≥ 99 ≥ 99.7 ≥ 99.9 99
1.3724 1.3852 1.4097 1.4195 1.4975 1.4938 1.4929
0.6554 0.6797 0.7266 0.7452 0.8736 0.8622 0.8567
Refractive Index measurement uncertainty: u(T) = 0.06 °C, u(nD) = 0.0002 with a 95% confidence interval. Density measurement uncertainty: u(T) = 0.02 °C, u(ρ) = 0.0004 g/cm3 with a 95% confidence interval.
study.
2 3 2 7 TC (K ) = 373.49FRI MW
(7)
2.2. Modeling methodology
8 5 −1 PC (bar) = 2.70 × 10 4FRI MW
(8)
−8 9 6 5 C (cm3 mol) = 0.47FRI V MW
(9)
The classical example for applying equations of state to pure components or mixtures is to use defined components with known compositions and critical properties reported in the literature [6]. However, many applications of EOS models are for complex systems – such as crude oils – for which a full, detailed compositional analysis is infeasible to obtain experimentally and the critical properties of some components are either not reported or might not be readily measurable. For these cases, researchers have developed characterization techniques with correlations to calculate critical properties for the otherwise undefined fractions of the system. One such example is the characterization technique of Pedersen et al. [7] that uses density and molecular weight correlations to specify the critical properties for the heavy carbon number fractions of a crude oil [8]. A similar modeling approach is proposed in this section, with critical properties for pseudo-pure fractions calculated from refractive index or density measurements. The critical properties for each pseudocomponent are calculated using the EV correlations with the refractive index and molecular weight measurements that are reported in Table 3. Critical properties for pure components used in this modeling study are shown in Table 4. To compare the approach proposed in this section (called the “lumped” approach because the multicomponent solvent is treated as a single pseudo-component) to the standard modeling approach, saturation curves are generated with PR EOS for both approaches. The systems modeled in this section are formed by adding a light key component (propane, hydrogen sulfide or carbon dioxide) in various proportions to the solvent fractions reported in Table 3 with 1-methylnaphthalene (1-MN) as the heavy key component. In the lumped approach, critical properties of the pseudo-pure solvents are calculated from the EV correlations, while the standard approach uses defined components with critical properties reported in the literature [6]. A graphical representation of the standard (Std) and lumped (EV) modeling approach is shown in Fig. 1.
ω = 8.62 ×
−1 10−4FRI MW
(10)
where MW has units of g/mol and FRI is calculated from refractive index or density measurements performed at 20 °C and 1 atm. The expression relating FRI to refractive index is given as:
FRI =
nD2 −1 nD2 + 2
(11)
where nD is the refractive index measured at the yellow sodium-D line.If refractive index is not reported, then FRI can be calculated from a mass density correlation proposed by Vargas and Chapman [4] given in Eq. (12).
FRI = 0.5054ρ−0.3951ρ2 + 0.2314ρ3 ≈
ρ 3
(12) 3
where ρ is mass density in units of g/cm . Applying the EV correlations to the CEOS framework offers two unique approaches to improving phase behavior predictions. First, critical properties for multicomponent nonpolar solvent mixtures can be calculated from simple volumetric or optical measurements without knowledge of the composition of the mixture. Second, binary interaction parameters can be calculated algebraically from density or refractive index measurements at ambient temperature rather than tuned to bubble pressure data and then subsequently used to improve phase behavior predictions for mixtures. These approaches are investigated in subsequent sections. 2. Modeling thermodynamic properties using critical property correlations applied to mixed solvents 2.1. Experimental methodology
2.3. Modeling results Experimental data for density and refractive index of nonpolar hydrocarbon mixtures were recently reported by Wang et al. [5] Density measurements were performed with the Anton Paar DMA 4500 digital vibrating U-tube densitometer and refractive index measurements were performed with the automatic Anton Paar WR refractometer under the sodium-D wavelength with a measuring range of 1.30–1.72. Table 2 provides a list of suppliers and chemical purity for the studied components. Wang et al. [5] also performed refractive index and density measurements on the mixtures reported in Table 3. These mixtures constitute the pseudo-pure fractions that are used in this modeling
Phase diagrams for the studied systems are shown in Fig. 2 with critical points represented by open squares. The molar amounts of these systems are given by the ratios associated with each curve. For example, 80:20:10 means that 80 mol of the mixture consists of the light key, 20 mol is the mixed solvent and 10 mol is the heavy key. For the lumped solvent approach, the molar amount of the solvent fraction is precisely the second number in the amount vector. In the standard de-lumped approach, that solvent amount is subdivided into defined components and molar amounts according to the compositions reported in Table 3. 638
Fuel 237 (2019) 637–647
C.J. Sisco et al.
Table 3 Properties of selected mixtures measured at 20 °C and 1 bar to be used as pseudo-components [5]. Mixtures
B1 B2 B3 T1 T2 T3
Amount (mol%)
MW
nD
ρ
Comp 1
Comp 2
Comp 3
z1
z2
z3
g/mol
–
g/cm3
n-Heptane n-Heptane Toluene n-Hexane n-Hexane Toluene
n-Dodecane Toluene p-Xylene n-Decane n-Dodecane Benzene
– – – n-Dodecane Toluene p-Xylene
63.0 47.9 53.6 33.5 33.3 33.3
37.0 52.1 46.4 33.3 33.4 33.7
– – – 33.2 33.3 33.0
126.18 96.00 98.64 132.83 116.24 92.04
1.4047 1.4345 1.4957 1.4080 1.4254 1.4968
0.716 0.764 0.864 0.722 0.751 0.867
Table 4 Simulation parameters for studied pure components.
α=
Component
MW g/mol
Tc K
Pc bar
ω –
C3 H2S CO2 n-C6 n-C7 n-C10 n-C12 Benzene Toluene p-Xylene 1-MN
44.096 34.081 44.010 86.175 100.202 142.285 170.339 78.114 92.141 106.168 142.200
369.85 373.10 304.13 507.82 540.13 617.70 658.20 562.16 591.79 616.20 772.00
42.48 90.00 73.77 30.18 27.27 21.20 18.24 48.98 41.04 35.11 36.00
0.1520 0.1005 0.2230 0.2979 0.3498 0.4890 0.5750 0.2080 0.2630 0.3180 0.3416
Std
1 ⎛ ∂V ⎞ 1 ∂V κ=− ⎛ ⎞ V ⎝ ∂T ⎠P V ⎝ ∂P ⎠T
(13)
The phase diagrams (Fig. 2) and parity plots (Figs. 3–5) show very good agreement between the proposed modeling approach (EV) and standard modeling approach (Std), demonstrating that the EV correlations provide representative solvent simulation parameters that describe nonpolar solvent phase behavior in an almost identical manner to the standard modeling approach. This is not, however, a comparison to experimental data. These plots convey only that the EV correlations applied within the CEOS framework describes phase behavior accurately with respect to the standard approach for obtaining compositions and simulation parameters for the CEOS. For mixtures that are not described accurately using the standard modeling approach (Std), the proposed modeling approach (EV) will likely show similar inaccuracies.
EV
3. Discussion Bar plots for cohesive energy and solubility parameter for the pure components studied in this work are shown in Figs. 6 and 7, respectively, showing that PR EOS reproduces the trend in cohesive energy much better than it does for solubility parameter. The reference data are obtained from a correlation for solubility parameter proposed by Wang [9].
δ = 52.042FRI + 2.904
(14)
With the definition of solubility parameter as the square root of the residual internal energy multiplied by density, Wang’s solubility parameter correlation can be rewritten in terms of residual internal energy.
−U res =
MW (52.042FRI + 2.904)2 ρ
(15)
Based on solubility parameter alone, one might expect that 1-MN would be always the appropriate choice as heavy key for the selected mixtures and that the solvent components are sufficiently different from 1-MN that the latter will be largely driving the phase-split at the dew point. However, the poor volumetric property predictions that are a well-known drawback of the CEOS cause the solubility parameter of 1MN to be under-predicted, as compared to Wang’s correlation, and the cohesive energy of n-C12 is greater than that of 1-MN at some conditions. Because the form of the equation of state and the mixing rules remain consistent, deviations between the lumped (EV) and de-lumped (Std) approaches can be explained by only two factors: the accuracy of the EV correlations to describe the solvent to which they are applied and the implicit assumption that all components in the lumped solvent partition equivalently. Fig. 8, which shows plots for cohesive energy and solubility parameter of the solvent mixtures, suggests that the EV correlations describe nonpolar solvent mixtures almost equivalently to the standard de-lumped modeling approach. However, regardless of the accuracy of the correlations, the lumped (EV) approach will always deviate from the de-lumped (Std) approach to some degree because of the limitation of the partitioning assumption. The magnitude of the deviations, though, might be negligible
Fig. 1. Standard (Std) and lumped (EV) modeling approaches. In the standard approach (Std), the solvent composition is obtained from Table 3 and simulation parameters are obtained from Table 4. In the lumped approach (EV), the solvent is treated as a single pseudo-component with simulation parameters calculated from the EV correlations given the experimental refractive index measurements in Table 3. The light key (LK) and heavy key (HK) are treated the same in each approach.
In addition to saturation pressure, various properties of the bulk solvent phase are of interest, and it is important that the proposed modeling approach does not introduce significant error to these predictions. Figs. 3–5 show parity plots comparing thermal expansion coefficient, isothermal compressibility, and density, respectively, of the saturated liquid phase calculated by the proposed approach (EV) to the standard modeling approach (Std). The diagonal line on the parity plots represents perfect agreement between the EV and standard modeling approaches. Points lying below the 45° line on the parity plots signify under-predictions by the EV approach with respect to the standard approach. The thermal expansion coefficient (α ) and isothermal compressibility (κ ) are calculated by: 639
Fuel 237 (2019) 637–647
C.J. Sisco et al.
P (bar)
250
250
C3 / B1 / 1-MN
200
200
150
150 80:20:10
100
100 250
P (bar)
200
300
500
700
80:20:10
P (bar)
500
700
80:20:10
100
300
500
700
C3 / B3 / 1-MN
150
50:50:10
100
80:20:10
100
50:50:10
50 20:80:10
0 100
P (bar)
20:80:10
0
200
50 300
500
20:80:10
0
700
100 250
C3 / T1 / 1-MN
200
150
300
500
700
CO2 / T2 / 1-MN 80:20:10
150
100
100
80:20:10 50:50:10
50 100
300
500
700
100 250 200
80:20:10
100
20:80:10
0
H2S / T3 / 1-MN
150
50:50:10
50
20:80:10
0
P (bar)
50:50:10
250
H2S / B3 / 1-MN
150
200
80:20:10
50
20:80:10
300
700
H2S / B2 / 1-MN
100
50:50:10
100
250
500
150
0
200
300
200
50
250
100 250
CO2 / B2 / 1-MN
100
200
20:80:10
0
150
250
50:50:10
50
20:80:10
0
80:20:10
100
50:50:10
50
H2S / B1 / 1-MN
300
500
CO2 / T3 / 1-MN
80:20:10
150
50:50:10
100
50:50:10
50
50 20:80:10
0 100
300
500
700
20:80:10
0
700
100
300
500
700
T(K)
T(K)
Fig. 2. Phase diagrams generated with PR EOS for the solvents listed in Table 3 with various light keys (LK) and 1-MN as the heavy key (HK). Solid lines: lumped solvent approach (EV). Dashed lines: standard de-lumped approach (Std). Molar amount vectors given by [LK:Solvent:HK].
cases, the equi-partitioning assumption causes very little deviation to property predictions. However, if n-butane were added to the solvent, the bubble pressure calculation for the lumped and de-lumped approaches would not agree as closely. Propane and n-butane have similar volatilities, and the solubility of n-butane in the solvent has a nonnegligible effect on bubble pressure that would be obfuscated by lumping n-butane with the heavy components. The same argument applies for the dew pressure calculation. The more dissimilar the heavy key is from the components that constitute the solvent, the lesser impact lumping the solvent will have on dew pressure predictions. Though all systems studied show very good saturation pressure
depending on the system and how the solvent is lumped. In general, saturation pressure predictions for the lumped and delumped approaches agree closely when the component(s) most responsible for driving the phase-split are well characterized and their partitioning with respect to the solvent are only negligibly affected by lumping. For example, the calculated bubble pressure of a mixture containing propane with toluene, p-xylene and 1-MN will be essentially the same with the lumped and de-lumped approaches because the solubility of propane in the solvent is only negligibly affected by lumping the heavy components and only a small amount of the heavy components will partition into the vapor phase at the bubble point. In these 640
Fuel 237 (2019) 637–647
C.J. Sisco et al. 6
6
Fig. 4. Isothermal compressibility (κ) parity plots comparing the lumped solvent approach (EV) on the x-axis to the standard de-lumped modeling approach (Std) on the y-axis with PR EOS. The diagonal line represents perfect agreement between the EV and Std approach.
4. Calculating binary interaction parameters from critical property correlations
HS
CΎ
αStd (10-3 K-1)
αStd (10-3 K-1)
Fig. 3. Thermal expansion coefficient (α) parity plots comparing the lumped solvent approach (EV) on the x-axis to the standard de-lumped modeling approach (Std) on the y-axis with PR EOS. The diagonal line represents perfect agreement between the EV and Std approach.
process accurately in the dew region, it is important that there is sufficient resolution in the heavy components that are driving the phasesplit. When one of the solvent components that is important to this process is lumped in with other solvent components, then its effect is diluted and the prediction suffers. When the difference in the cohesive energy between the heavy key and solvent components is sufficiently large, which is the case for the solvents that do not contain n-C12 (B2, B3, T3), then the deviations between the lumped and standard approaches are negligible. The modeling study performed in this work is a proof-of-concept for the EV correlations and the lumped solvent approach. With the accuracy of the approach established, it is important to understand how and when the proposed approach can be applied to practical systems. The approach is primarily of interest for defining the simulation parameters of unknown fractions of a complex mixture. If all components and their compositions in a mixture are known in advance, then the standard approach would be always preferable unless the mixture exceeds some limitation on the number of components allowed in a simulation. In these cases, the EV correlations offer a method to lump similar components into a solvent fraction and calculate the critical properties of the pseudo-component. However, the primary reason these correlations were developed was to describe fractions of complex crude oils for which a detailed compositional analysis is not readily available. The crude oil can be lumped into as many or few fractions as desired and then the critical properties of each fraction can be obtained from simple refractive index or density measurements performed on each fraction at ambient conditions. The number of component fractions and whether the light end and/or heavy end requires more resolution depends on which properties are being predicted and the required accuracy of the predictions. If only bubble pressures are desired, then a simple lumping approach with all components heavier than n-hexane lumped into a single fraction is often sufficient, as shown by Abutaqiya et al. [10]. If dew pressures or asphaltene onsets are desired, then the heavy end must be characterized in more detail, as shown in the more customary crude oil characterization methodologies [7,11–14]. It was expected that solubility parameter would offer an indication as to why certain solvent mixtures showed consistently worse predictions than others, but there was no apparent trend in solubility parameter that suggested when the lumped solvent approach predictions might suffer. Cohesive energy, on the other hand, did show such a trend. When the difference in cohesive energy between the light and/or heavy key and all solvent components is large, then the lumped solvent approach will introduce only negligible error to saturation pressure predictions. The Wang correlation for residual internal energy in Eq. (15) (which is equivalent to cohesive energy) offers a means to check on whether de-lumping the crude oil fractions might be appropriate. Eq. (15) requires only refractive index (or density) and molecular weight measurements, and PR EOS with EV correlations was shown to match the residual internal energy correlation quite well, suggesting that cohesive energy can be used as a heuristic for determining the number of fractions into which the crude oil should be separated in order to generate simulation parameters that are representative of the actual mixture.
CO
4
2 B1 0 2 4 αEV (10-3 K-1)
CΎ
4 2 B2
CO
HS
6
0 6
CΎ
αStd (10-3 K-1)
6
αStd (10-3 K-1)
HS
0 0
4
2 B3 0
2 4 αEV (10-3 K-1) CO
HS
6
CΎ
4
2 T1 0
6
2 4 αEV (10-3 K-1) CO
HS
6
0 6
CΎ
αStd (10-3 K-1)
0
αStd (10-3 K-1)
CO
4
2 T2 0
2 4 αEV (10-3 K-1) CO
HS
6
CΎ
4
2 T3 0
0
2 4 αEV (10-3 K-1)
6
0
2 4 αEV (10-3 K-1)
6
4.1. Modeling methodology predictions – due mostly to the accuracy of the EV correlations – some solvents (B1, T1, T2) showed much more deviation between the lumped and standard approaches, particularly in the dew region, than the other solvents. The primary reason this deviation occurs is because the heavy key (1-MN) and the heaviest component in the solvent for these cases (n-C12) have similar cohesive energy. To describe the phase-splitting
The standard procedure for obtaining binary interaction parameters is by tuning to saturation pressure data for a given component pair. Other methodologies have been developed that use critical volumes [15] or a group contribution approach [16–20] to calculate binary interaction parameters, thereby eliminating the need for tuning. The approach presented in this work uses the EV correlations and the CEOS 641
Fuel 237 (2019) 637–647
C.J. Sisco et al.
mixing rule directly to calculate binary interaction parameters. Only a single data point for the binary mixture at any composition is required to implement this approach, but more data at different compositions can be used to refine the calculations. The necessary data for this approach are very easy to obtain, as the EV correlations require only simple measurements of a solvent at ambient conditions to define its critical properties. − Writing the expression for a from the van der Waals mixing rule for a binary mixture yields: −
a = x12 a1 + x 22 a2 + 2x1 x2 a1a2 (1−k12)
(16)
This is the standard mixing rule for CEOS, but the approach proposed here can be extended to any mixing rule. The optimal value of k12 is specific to a given component pair and the mixing rule. For binary mixtures of similar components with the van der Waals mixing rules, the optimal value for k12 is usually near zero while stronger deviations from zero are observed for dissimilar components. Large deviations from a k12 value of zero is generally interpreted as a deficiency in the mixing rule to describe the interaction energy of components 1 and 2. − If the process is reversed and a is already known for a given composition, then k12 can be calculated explicitly by rearranging the mixing rule to obtain: −
k12 = 1−(a−x12 a1−x 22 a2) (2x1 x2 a1a2 )
In this section, a methodology is proposed for calculating a from refractive index measurements and the EV correlations and then solving explicitly for k12 . To illustrate, consider an equimolar mixture of components 1 and 2. FRI can be calculated for this mixture either from refractive index or density measurements at 20 °C and 1 atm, and the molecular weight of the mixture can be measured experimentally or calculated from the compositions and pure-component molecular weights. With FRI and MW known, the critical properties of the mixture − (Tc, mix , Pc, mix , ωmix ) are calculated from the EV correlations and a is calculated the same as in Eq. (2), treating the binary mixture as a − − pseudo-component: a = Ωa (R2Tc, mix 2αmix ) Pc, mix . Because a is calculated for a specific mixture composition ( x1 = 0.5 for this example) and a1 and a2 are known from the critical properties of the pure components (either obtained from the literature or from the EV correlations), all terms on the right-hand side of Eq. (17) are known and k12 can be calculated explicitly. Additionally – because of the temperature-dependence implicit in the α term – k12 can be calculated as a function of temperature without any need for further experiments. There are no restrictions on the compositions of the binary mixtures or on the number of compositions for which the experiment must be performed. In principle, the calculation of k12 requires only one experimental measurement for refractive index or density at a known mixture composition of components 1 and 2, as illustrated in this example, but a more representative k12 value would likely be obtained if measurements for multiple compositions of components 1 and 2 were considered. To evaluate fairly the proposed methodology for calculating k12 , the γ −φ formulation is preferred to the more standard φ−φ formulation because the latter approach can sometimes show poor predictions for pure-component saturation pressures. The binary interaction parameter cannot correct pure-component properties, so showing poor predictions for these properties might mislead the reader about the usefulness of the proposed methodology. The γ −φ formulation for the equi-fugacity criterion is written as:
Cohesive Energy, -Ures/RT
30 25 20 15 10 5 0
Fig. 6. Cohesive energy for selected pure components at 20 °C and 1 bar from Eq. (15) [blue] and PR EOS [red]. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 22 21 20 19 Solubility Parameter,
(17) −
Fig. 5. Density (ρ) parity plots comparing the lumped solvent approach (EV) on the x-axis to the standard de-lumped modeling approach (Std) on the y-axis with PR EOS. The diagonal line represents perfect agreement between the EV and Std approach.
18 17 16 15 14 13 12
φ sat yi φiV̂ P = x i φiL̂ ⎛⎜ i ⎟⎞ Pisat ⎝ φi ⎠
Fig. 7. Solubility parameter for selected pure components at 20 °C and 1 bar from Eq. (14) [blue] and PR EOS [red]. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
φiV̂
(18)
and φi ̂ are the fugacity coefficients of the vapor and liquid, where respectively, calculated from the vapor composition {yi } and liquid composition {x i} at the pressure P and temperature T of interest. φisat is the fugacity coefficient of pure component i evaluated at Pisat and φi is 642
L
Fuel 237 (2019) 637–647
C.J. Sisco et al. 22 21 25
Solubility Parameter,
Cohesive Energy, -Ures/RT
30
20 15 10 5
20
19 18 17 16 15 14 13
0
12 B1
B2
B3
T1
T2
T3
B1
B2
B3
T1
T2
T3
Fig. 8. Cohesive energy and solubility parameter plots for solvent mixtures at 20 °C and 1 bar from PR EOS with standard (Std) approach [blue] and lumped solvent (EV) approach [red]. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
The calculation of k12 was completed via two approaches. In the EV1 approach, both the pure component critical properties and mixture properties were calculated with the EV correlations, while in the EV2 approach, the pure component critical properties were taken from the literature and mixture properties were calculated with the EV correlations.
Table 5 Density or refractive index values at 20 °C and 1 bar for selected binary mixtures. Mixtures
Amount (mol%)
MW
nD
ρ
Comp 1
Comp 2
z1
z2
g/mol
-
g/cm3
Benzene Toluene Cyclohexane n-Heptane Cyclohexane Cyclohexane
Toluene n-Decane n-Decane Toluene n-Dodecane n-Heptane
50.0 50.0 50.0 47.9 50.0 50.0
50.0 50.0 50.0 52.1 50.0 50.0
85.13 117.21 113.22 96.17 127.25 92.18
– – – 1.4345 – –
0.8721 0.7780 0.7471 – 0.7567 0.7225
6. Discussion In the EV1 approach, the critical properties for both the binary mixture and constituent pure components are calculated from the EV correlations. When the EV correlations describe mixture properties poorly, it is likely that it also describes one or both constituent pure component properties poorly. Because k12 is quantifying a difference between the mixture and pure component properties, it is not necessary that the EV correlations describe the mixture and pure component properties well for the EV1 approach to be effective. Rather, it is only important that the difference between the mixture and pure component properties is captured accurately. The EV1 approach benefits from a cancellation of errors when both the pure components and mixture are described poorly. In the EV2 approach, the critical properties of the binary mixture are calculated from the EV correlations, whereas the constituent pure component critical properties are taken from the literature. Because the literature values will necessarily produce good predictions for saturation pressure for the pure components, the EV correlations must also
the fugacity coefficient of pure component i evaluated at P . All fugacity coefficients are calculated using the PR EOS, and saturation pressure of pure component i (Pisat ) is calculated with correlations provided in the DIPPR database [6] and is a function of temperature only.
5. Modeling results P-xy curves for binary mixtures with k12 calculated from the mixing rule and k12 = 0 are compared to reference VLE data generated from NIST REFPROP software [21]. The binary mixtures studied are listed in Table 5. The results are presented in Figs. 9–14 and the mean absolute percent deviation is presented in Table 6. 0.16
0.40
298.15 K
P (bar)
P (bar)
323.15 K
0.30
0.12 0.08
0.20 0.10
0.04
0.00
0.00
0
0.2
0.4
x1,y1
0.6
0.8
1
0
1.00
0.4
x1,y1
0.6
0.8
1
0.6
0.8
1
2.00 348.15 K
373.15 K
0.80
1.60
0.60
1.20
P (bar)
P (bar)
0.2
0.40 0.20
0.80 0.40
0.00
0.00 0
0.2
0.4
x1,y1
0.6
0.8
1
0
0.2
0.4
x1,y1
Fig. 9. P-xy diagram for Benzene (1)/Toluene (2) with PR EOS. Open circles: NIST REFPROP [21]. Black dashed lines: k12 = 0 . Blue lines: EV1. Red lines: EV2. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 643
Fuel 237 (2019) 637–647
C.J. Sisco et al. 0.04
0.16
298.15 K
0.02
0.01
0.08 0.04
0.00
0.00 0
0.40
0.2
0.4
x1,y1
0.6
0.8
1
0
0.2
0.4
x1,y1
0.6
0.8
1
0.6
0.8
1
0.80
348.15 K
373.15 K 0.60
P (bar)
0.30
P (bar)
323.15 K
0.12
P (bar)
P (bar)
0.03
0.20
0.40
0.20
0.10
0.00
0.00
0
0.2
0.4
x1,y1
0.6
0.8
0
1
0.2
0.4
x1,y1
Fig. 10. P-xy diagram for Toluene (1)/n-Decane (2) with PR EOS. Open circles: NIST REFPROP [21]. Black dashed lines: k12 = 0 . Blue lines: EV1. Red lines: EV2. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 0.16
0.40
298.15 K
323.15 K 0.30
P (bar)
P (bar)
0.12
0.08
0.10
0.04
0.00
0.00 0 1.00
0.2
0.4
x1,y1
0.6
0.8
0
1
0.2
0.4
x1,y1
0.6
0.8
1
0.6
0.8
1
2.00
348.15 K
373.15 K
0.80
1.60
0.60
1.20
P (bar)
P (bar)
0.20
0.40 0.20
0.80 0.40
0.00
0.00 0
0.2
0.4
x1,y1
0.6
0.8
1
0
0.2
0.4
x1,y1
Fig. 11. P-xy diagram for Cyclohexane (1)/n-Decane (2) with PR EOS. Open circles: NIST REFPROP [21]. Black dashed lines: k12 = 0 . Blue lines: EV1. Red lines: EV2. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
approach performs poorly, such as the cyclohexane / n-dodecane mixture, one component yields much better saturation pressure predictions with the EV correlations than the other such that the errorcanceling effect is only partially realized. In cases such as benzene/ toluene, in which the EV correlations under-predict saturation pressures for both components, the deviations are roughly equal in magnitude and direction so that the error-canceling effect is more fully realized. Fortunately, the EV correlations are amenable for re-tuning to specific families of components or even a specific binary pair. All EV correlations have the same standard form given by:
produce good predictions for the mixture saturation pressure if the EV2 approach is to generate suitable values for k12 . When the EV2 approach fails to generate satisfactory values for k12 , it is due to the EV correlations describing one or both constituent pure components poorly and the deviations occurring in the same direction such that the errors are additive instead of negating. For example, saturation pressures for benzene and toluene are both under-predicted by the EV correlations, so the mixture saturation pressure is also under-predicted. Because the correlations-only approach (EV1) cancels the error in mixture saturation pressure with correspondingly poor pure component saturation pressure predictions, the EV1 approach outperforms the EV2 approach in these cases. In cases such as the toluene/n-decane mixture, the EV correlations over-predict saturation pressure for n-decane by roughly the same order of magnitude that they under-predict the saturation pressure for toluene. This yields a fortunate cancellation of errors for the mixture saturation pressure and correspondingly satisfactory results for k12 with the EV2 approach. To maintain consistency, it is best to use the EV1 approach always because of the error-canceling effect. The cases for which the EV1
Y = aFRI bMW c
(19)
where Y is the property of interest and a , b , and c are tuning constants. The tuning constants of the EV correlations can be fit to reproduce the properties for a specific set of compounds to increase the accuracy of certain binary systems. The same EV1 approach proposed here but with new constants can then be used to improve the results for the k12 calculation.
644
Fuel 237 (2019) 637–647
C.J. Sisco et al. 0.20
0.07
323.15 K
298.15 K 0.18 0.16
P (bar)
P (bar)
0.06
0.05
0.14
0.04
0.12
0.03
0.10 0
0.2
0.4
0.6
0.8
0
1
0.2
0.4
x1,y1 0.50
1.10
0.8
1
0.6
0.8
1
373.15 K
348.15 K 0.46
1.00
0.42
P (bar)
P (bar)
0.6
x1,y1
0.38
0.90 0.80
0.34
0.70
0.30 0
0.2
0.4
0.6
0.8
0
1
0.2
0.4 x1,y1
x1,y1
Fig. 12. P-xy diagram for n-Heptane (1)/Toluene (2) with PR EOS. Open circles: NIST REFPROP [21]. Black dashed lines: k12 = 0 . Blue lines: EV1. Red lines: EV2. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 0.16
0.40
298.15 K
323.15 K
0.30
P (bar)
P (bar)
0.12 0.08
0.04
0.20 0.10
0.00
0.00 0
0.2
0.4
x1,y1
0.6
0.8
1
0
1.00
2.00
0.60
1.20
P (bar)
P (bar)
1.60
0.40 0.20
0.4
x1,y1
0.6
0.8
1
0.6
0.8
1
373.15 K
348.15 K 0.80
0.2
0.80 0.40
0.00
0.00
0
0.2
0.4
x1,y1
0.6
0.8
1
0
0.2
0.4
x1,y1
Fig. 13. P-xy diagram for Cyclohexane (1)/n-Dodecane (2) with PR EOS. Open circles: NIST REFPROP [21]. Black dashed lines: k12 = 0 . Blue lines: EV1. Red lines: EV2. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
7. Conclusions
the lumped solvent modeling approach produces very good predictions for saturation pressure in the dew region. When the heavy key and solvent components have similar cohesive energy, dew pressure predictions suffer if the solvent is lumped. The EV correlations were also applied to binary mixtures to calculate binary interaction parameters algebraically from refractive index or density measurements at ambient conditions instead of the standard method of tuning to saturation pressure data, which are more laborious to obtain. Because the calculation of k12 requires that components 1 and 2 are known in advance, the accuracy of the EV correlations in describing the constituent pure components can also be known. The utility of this approach for generating good k12 values for a given component pair can then be estimated in advance. However, even if the EV correlations with constants given in this work do not describe a component pair with sufficient precision, it is relatively straightforward to tune the correlations to match a specific set of components and then generate k12 values using the same approach proposed in this work. The EV correlations were originally tuned to many nonpolar hydrocarbons of different families (n-alkanes, cyclics, and aromatics), which means that
In this work, critical property correlations recently proposed by Evangelista and Vargas [1] were evaluated for their accuracy in describing nonpolar solvent mixtures. The correlations, along with the proposed lumped solvent approach, show very good predictive capability for saturation pressures and bulk liquid properties like density, thermal expansion, and isothermal compressibility. This suggests that the approach is suitable to be applied to nonpolar crude oil fractions for which the composition is unknown. The EV correlations offer a simple expression for relating easily measured volumetric or optical properties to critical properties, and the lumped solvent approach shows small deviations with respect to the standard modeling approach. Additionally, it was shown that cohesive energy calculated by PR EOS with EV correlations matches well with the Wang correlation for cohesive energy, and differences in cohesive energy between the heavy key and the components comprising the solvent was shown to correlate with accuracy of the lumped solvent approach in matching dew point predictions. When the heavy key and solvent are sufficiently different, 645
Fuel 237 (2019) 637–647
C.J. Sisco et al. 0.40
0.15
323.15 K
0.13
0.35
0.11
0.30
P (bar)
P (bar)
298.15 K
0.09
0.25 0.20
0.07
0.15
0.05 0
0.2
0.4
0.6
0.8
0
1
0.2
0.4
0.90
1.80
348.15 K
0.80
0.8
1
0.6
0.8
1
373.15 K
1.60
0.70
P (bar)
P (bar)
0.6 x1,y1
x1,y1
0.60
1.40
1.20
0.50 0.40
1.00 0
0.2
0.4
0.6
0.8
0
1
0.2
x1,y1
0.4 x1,y1
Fig. 14. P-xy diagram for Cyclohexane (1)/n-Heptane (2) with PR EOS. Open circles: NIST REFPROP [21]. Black dashed lines: k12 = 0 . Blue lines: EV1. Red lines: EV2. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Table 6 Binary interaction parameters (kij ) and error analysis for selected binary mixtures. T / °C
Benzene + Toluene
Toluene + n-Decane
Cyclohexane + n-Decane
n−Heptane + Toluene
Cyclohexane + n-Dodecane
Cyclohexane + n-Heptane
25 50 75 100 25 50 75 100 25 50 75 100 25 50 75 100 25 50 75 100 25 50 75 100
AAPD in BP
kij
AAPD in DP
EV1
EV2
EV1
EV2
kij = 0
EV1
EV2
kij = 0
0.0001 0.0005 0.0004 0.0004 0.0119 0.0111 0.0103 0.0095 0.0166 0.0156 0.0146 0.0136 −0.0038 −0.0039 −0.0041 −0.0043 0.0335 0.0320 0.0305 0.0289 0.0003 0.0003 0.0003 0.0003
−0.1037 −0.1032 −0.1027 −0.1023 −0.0117 −0.0118 −0.0118 −0.0118 −0.0335 −0.0301 −0.0267 −0.0233 −0.0364 −0.0351 −0.0337 −0.0324 −0.0181 −0.0155 −0.0130 −0.0104 −0.0622 −0.0561 −0.0500 −0.0438
0.577 0.361 0.252 0.187 11.442 10.336 9.786 9.404 18.361 15.732 13.735 12.153 4.632 4.012 3.381 2.772 46.148 40.729 36.078 31.979 1.181 1.116 1.300 1.566
25.178 22.552 20.241 18.221 3.186 1.791 1.043 1.765 10.747 7.192 4.106 1.547 13.412 11.473 9.738 8.200 3.225 6.350 8.913 10.795 15.120 12.014 9.214 6.759
0.413 0.237 0.155 0.118 4.003 4.545 5.270 5.886 7.284 7.034 6.920 6.827 3.476 2.958 2.413 1.877 15.826 16.122 16.142 15.892 1.092 1.037 1.229 1.501
0.504 0.212 0.157 0.110 0.884 1.362 1.983 2.677 0.771 1.152 1.588 2.033 3.623 3.329 2.913 2.449 0.179 0.826 1.390 1.986 0.723 0.784 1.009 1.313
25.273 23.031 20.931 19.002 1.676 1.531 1.042 0.657 1.481 1.563 1.425 1.055 12.572 10.873 9.305 7.885 0.060 0.541 0.911 1.286 15.001 11.932 9.168 6.747
0.367 0.128 0.142 0.160 0.659 0.721 1.034 1.612 0.542 0.710 0.940 1.215 2.601 2.381 2.030 1.623 0.106 0.651 1.076 1.493 0.663 0.722 0.948 1.256
some components will not be as well-described by the correlations as others. Because all EV correlations have the same standard form, a systematic approach for tuning the constants to capture more precisely the properties of certain component sets can be easily implemented.
[6] [7]
References
[8] [9]
[1] Evangelista RF, Vargas FM. Prediction of the temperature dependence of densities and vapor pressures of nonpolar hydrocarbons based on their molecular structure and refractive index data at 20 °C. Fluid Phase Equilib 2018;468:29–37. [2] Soave G. Equilibrium constants from a modified Redlich-Kwong equation of state. Chem Eng Sci 1972;27(6):1197–203. [3] Peng D-Y, Robinson DB. A new two-constant equation of state. Ind Eng Chem Fundam 1976;15(1):59–64. [4] Vargas FM, Chapman WG. Application of the One-Third rule in hydrocarbon and crude oil systems. Fluid Phase Equilib 2010;290(1–2):103–8. [5] Wang F, Evangelista RF, Threatt TJ, Tavakkoli M, Vargas FM. Determination of
[10]
[11]
[12]
646
volumetric properties using refractive index measurements for nonpolar hydrocarbons and crude oils. Ind Eng Chem Res 2017. Wilding WV, Rowley RL, Oscarson JL. DIPPR® Project 801 evaluated process design data. Fluid Phase Equilib 1998;150:413–20. Pedersen KS, Blilie AL, Meisingset KK. PVT calculations on petroleum reservoir fluids using measured and estimated compositional data for the plus fraction. Ind Eng Chem Res 1992;31(5):1378–84. Panuganti SR, Vargas FM, Chapman WG. Property scaling relations for nonpolar hydrocarbons. Ind Eng Chem Res 2013;52(23):8009–20. Wang J. Predicting Asphaltene Flocculation in Crude Oils. Socorro, New Mexico: New Mexico Institute of Mining & Technology; 2000. Abutaqiya MIL, Panuganti SR, Vargas FM. Efficient algorithm for the prediction of pressure–volume–temperature properties of crude oils using the perturbed-chain statistical associating fluid theory equation of state. Ind Eng Chem Res May 2017;56(20):6088–102. Panuganti SR, Vargas FM, Gonzalez DL, Kurup AS, Chapman WG. PC-SAFT characterization of crude oils and modeling of asphaltene phase behavior. Fuel 2012;93:658–69. Tavakkoli M, Chen A, Vargas FM. Rethinking the modeling approach for asphaltene
Fuel 237 (2019) 637–647
C.J. Sisco et al.
[13]
[14]
[15]
[16]
[17]
precipitation using the PC-SAFT Equation of State. Fluid Phase Equilib 2016;416:120–9. Arya A, von Solms N, Kontogeorgis GM. Determination of asphaltene onset conditions using the cubic plus association equation of state. Fluid Phase Equilib 2015;400:8–19. Arya A, von Solms N, Kontogeorgis GM. Investigation of the gas injection effect on asphaltene onset precipitation using the cubic-plus-association equation of state. Energy Fuels 2016;30(5):3560–74. Gao G, Daridon J-L, Saint-Guirons 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. Jaubert J-N, Mutelet F. VLE predictions with the Peng-Robinson equation of state and temperature dependent kij calculated through a group contribution method. Fluid Phase Equilib 2004;224(2):285–304. Jaubert J-N, Vitu S, Mutelet F, Corriou J-P. Extension of the PPR78 model
[18]
[19]
[20]
[21]
647
(predictive 1978, Peng-Robinson EOS with temperature dependent kij calculated through a group contribution method) to systems containing aromatic compounds. Fluid Phase Equilib 2005;237(1–2):193–211. Soave G, Gamba S, Pellegrini LA. SRK equation of state: Predicting binary interaction parameters of hydrocarbons and related compounds. Fluid Phase Equilib 2010;299(2):285–93. Abudour AM, Mohammad SA, Robinson Jr. RL, Gasem KAM. Generalized binary interaction parameters for the Peng-Robinson equation of state. Fluid Phase Equilib 2014;383:156–73. Nasrifar K, Rahmanian N. Equations of state with group contribution binary interaction parameters for calculation of two-phase envelopes for synthetic and real natural gas mixtures with heavy fractions. Oil Gas Sci. Technol. 2018;73:7. Lemmon EW, Huber ML, McLinden MO. “NIST Standard Reference Database 23: Reference Fluid Thermodynamic and Transport Properties-REFPROP, Version 9.1.” In: Natl Std Ref Data Ser. NIST NSRDS -, May 2013.