© H. Aguilar-Cisneros et al., published by IFP Energies nouvelles, 2017

Notation

A : Parameter of the equation of state

A₁-A₃: Coefficients in Equation (5)

a₁-a₄: Coefficients in Equations (11)-(14)

B : Parameter of the equation of state

b₁-b₄: Coefficients in Equations (11)-(14)

C _k : Carbon number of defined hydrocarbons

c₁: Coefficient in Equation (11)

d₁: Coefficient in Equation (11)

FG : Number of functional groups types

F_Obj: Objective Function

H/C: Hydrogen carbon ratio

MW : Molecular Weight

N : Number of data points

P : Pressure

P_c: Critical Pressure

p_c: Contribution to P_c by functional groups

P_Sat: Saturation Pressure

R : Universal constant of gases, 8.31441 J/mol K

RAD%:

T : Temperature

T_c: Critical temperature

V : Molar volume, m³/mol

V_c: Critical volume, m³/mol

v_Ci: Contribution to V_c by functional groups

x : Mole fraction

Z : Compressibility factor

Greek letters

α : Coefficient in Equation (9)

β : Coefficient in Equation (9)

χ : Coefficient in Equation (9)

Δv : Group volume increment in Equations (4) and (5)

Φ : Fugacity coefficient

ρ : Mass density, kg/m³

ν : Number of functional groups

ω : Acentric factor

Superscripts

Cal: Calculated data

L: Liquid phase

Exp: Experimental data

Subscripts

c: Critical property

k : Carbon number k

i : Functional group i or pure component i

j : Undefined petroleum fractions j or pure component j

Introduction

Hydrogen solubility predictions are essential part in several processes in the petroleum and chemical industries; for instance, hydrocracking and hydrogenation processes are used to transform heavy oil into more valuable products. Hydrogen solubility data are required for designing and operating such processes and it is also needed in the corresponding kinetic models (Saajanlehto et al., 2014). Several authors have pointed out the importance of taking into account Vapor-Liquid Equilibria (VLE) in reaction modeling and its influence in conversion predictions (Chen et al., 2011; Hook, 1985; Pellegrini et al., 2008). A detailed review about VLE of hydrogen and petroleum fraction systems was presented by Chávez et al. (2014). They presented experimental data and modeling methods, most of which were conducted for systems composed of hydrogen and defined hydrocarbons, with few systems formed by hydrogen and undefined petroleum fractions. This may be attributable to the lack of experimental data, making a hard task to develop accurate prediction methods. Besides, binary systems containing hydrogen all exhibit a type III phase behavior in the classification scheme of Van Konynenburg and Scott (Privat and Jaubert, 2013) which cannot be predicted with cubic equations of state using classical mixing rules because of the quantum nature of hydrogen and the mixture size-asymmetry (Deiters, 2013). Despite these issues, most methods found in literature use cubic Equations Of State (EOS) because of their mathematical simplicity and their capability of modeling multicomponent mixtures over a wide range of temperatures and pressures. To overcome such limitations several authors have proposed to fit EOS parameters to VLE data (Qian et al., 2013); however, in general, their capability of prediction is reduced with such regressions. New approaches have been suggested in order to preserve the predictive features of cubic equations of state, such as that described by Jaubert et al. (2010). In such approach, binary interaction parameters are predicted by using a group contribution method, where classical cubic equations of state and excess free energy models (EOS/g^E) are combined, being applied successfully to binary mixtures of hydrogen with defined hydrocarbons (Qian et al., 2013).

A further challenge is the mathematical characterization of the undefined petroleum fraction; that is, the determination of its physical-chemical properties by mathematical means. In such characterization process, it is important to find a convenient number of pseudocomponents in which the undefined petroleum fraction will be divided. According to Lin et al. (1985), better results are obtained as the number of pseudocomponents increases when VLE of systems formed by hydrogen and undefined fractions are modeled. On the other hand, Chávez concluded that 30 pseudocomponents are enough to obtain good agreement between calculated and experimental data. The characterization procedure also involves estimation of the corresponding critical properties entering into cubic equations of state. These properties are frequently calculated from correlations (Firoozabadi, 1988; Lee and Kesler, 1975; Pedersen et al., 2004), but they do not always lead to the same results, because each one was developed for different sets of fluids at limited conditions. Thus, recently, group contribution methods have been recently extended to pseudocomponents, showing a good performance (Carreón-Calderón et al., 2012; Xu et al., 2015).

Carreón-Calderón et al. (2012) presented a new characterization procedure based on group contribution methods. They used selected functional groups to assign a hypothetical chemical structure to an undefined petroleum fraction by a minimization process of its corresponding Gibbs free energy. This methodology was applied to simulate VLE of heavy (Carreón-Calderón et al., 2014) and light hydrocarbons (Uribe-Vargas et al., 2016). Good results were achieved in both cases without making use of correlations to calculate critical properties or adjusted binary interaction parameters. In particular, they modeled vapor-liquid equilibria of heavy petroleum fluids, where no significant differences were shown regardless of the number of pseudocomponents used in modeling. These results show the feasibility of applying this characterization procedure to VLE of hydrogen with heavy undefined petroleum fraction. In this work, we add hydrogen functional group to the set of functional groups original proposed and suggest modifications to the characterization procedure in order to determine the hydrogen solubility properly. The entire methodology was tested with experimental data found in literature: Venezuelan Heavy Coking Gas Oil (HGCO) (Ji et al., 2013), Karamay Atmospheric Residue (KRAR) (Ji et al., 2013), Liaohe Atmospheric Residue (LHAR) (Ji et al., 2013), Venezuelan Atmospheric Residue (VNAR) (Ji et al., 2013), Canadian Light Virgin Gas Oil (CLVGO) (Cai et al., 2001), Canadian Heavy Virgin Gas Oil (CHVGO) (Cai et al., 2001) and Athabasca Bitumen (AB) (Lal et al., 1999).

1 Characterization Procedure

The approach presented by Carreón-Calderón et al. (2012) is used to determine the properties of the undefined petroleum fractions. In that methodology, a hypothetical structure is calculated with a set of selected functional groups, which depend on the group contribution method used. Carreón-Calderón et al. (2012, 2014) have suggested two group contribution methods; Joback-Reid (JR) (Joback and Reid, 1987) and Marrero-Gani (MG) (Marrero and Gani, 2001). JR method was suggested to simulate VLE of light petroleum fluids (Uribe-Vargas et al., 2016), while MG method showed better functionality in heavy petroleum fluids heavier than 14 API (Carreón-Calderón et al., 2014). In the latter, the C₇₊ fraction, which has unknown composition, was split in 3, 6 and 9 pseudocomponents, founding no significant difference between calculations and experiments regardless of the fraction splitting. Accordingly, in this work, the undefined petroleum fractions are treated as a single pseudocomponent in order to have a simpler characterization, where their corresponding properties are calculated through MG method.

Table 1 shows the functional groups of the MG method and the equivalent groups of the GCVOL (Elbro et al., 1991; Ihmels and Gmehling, 2003) method for determining liquid densities, as originally proposed by Carreón-Calderón et al. (2012). A quick review of these functional groups was made and it was found that the last functional groups are not equivalent each other (aC ≠ AC − C). Therefore, in this work the last functional group (aC) was changed to aC − C functional group from the same MG method.

An optimization process is used to estimate a hypothetical chemical structure of the undefined petroleum fraction, where its Gibbs free energy is the objective function to be minimized. Since the undefined fraction is treated like a pure component, minimize the liquid fugacity coefficient is equivalent to minimize the Gibbs energy. Thus, the objective function is:$$\underset{{\upsilon }_1,{\upsilon }_2,\cdots ,{\upsilon }_{\mathrm{FG}} \ge 0}{\min } {\varphi }_i^L\left[T_{\mathrm{C}i}\left({\upsilon }_1,{\upsilon }_2,\cdots ,{\upsilon }_{\mathrm{FG}}\right),P_{\mathrm{C}i}\left({\upsilon }_1,{\upsilon }_2,\cdots ,{\upsilon }_{\mathrm{FG}}\right),{\omega }_i\left({\upsilon }_1,{\upsilon }_2,\cdots ,{\upsilon }_{\mathrm{FG}}\right)\right]$$(1)where T_C, P_C and ω are the critical temperature, critical pressure and acentric factor, respectively, with FG being the number of functional groups types. The expression depends on the cubic EOS selected; in this work, the Peng-Robinson cubic Equation Of State PR-EOS (Peng and Robinson, 1976) is used to determine the fugacity coefficient as follows:$$\ln \left({\varphi }_i^L\right)=\left(Z^L-1\right) \ln \left(Z^L-B\right)-\frac{A}{2\sqrt{2} B}\ln \left(\frac{Z^L+\left(\sqrt{2}+1\right)B}{Z^L-\left(\sqrt{2}-1\right)B}\right)$$(2)where Z ^L  = PV ^L ⁄(RT) is the compressibility factor of the liquid phase, A and B are PR-EOS parameters. The properties of undefined petroleum fractions are usually given at standard conditions; hence, the pressure P and temperature T are set equal to 0.101325 MPa and 288.15 K, respectively; v ^L is the molar liquid volume and R is the gas constant.

This minimization problem is constrained by the molecular weight and density of the corresponding undefined petroleum fraction:$$\sum _{j=1}^{\mathrm{FG}}{\upsilon }_j {\mathrm{MW}}_j-{\mathrm{MW}}_i=0$$(3) $$\sum _{j=1}^{\mathrm{FG}}{\upsilon }_j \Delta v_j-\frac{{\mathrm{MW}}_i}{{\rho }_i^L}=0$$(4)

In the above expression, subscripts j and i indicate the functional group j of the undefined fraction i, MW is the molecular weight, υ_j, the number (coefficient) of a functional group, Δv_j the volume increment, the liquid density. The Δv_j is expressed as a function of temperature by means of$$\Delta v_j=A_1+A_{2}T+A_3{T}^2$$(5)where A₁, A₂ and A₃ are the parameters of the GCVOL approach (Ihmels and Gmehling, 2003).

Finally, to solve the minimization problem, Equation (2) requires critical properties (T_C, P_C and ω). These properties are calculated using Equations from (6) to (8), where, the temperature is in K, pressures in MPa and volume in m³/mol (Marrero and Gani, 2001).$$T_{\mathrm{C}}=231.239 \ln \left(\sum _{i=1}^{\mathrm{FG}}{\upsilon }_i t_{\mathrm{C}i}\right)$$(6) $$P_{\mathrm{C}}=0.1{\left(\sum _{i=1}^{\mathrm{FG}}{\upsilon }_i p_{\mathrm{C}i}+0.108998\right)}^{-2}+0.59827$$(7) $$V_{\mathrm{C}}=\left(7.95+\sum _{i=1}^{\mathrm{FG}}{\upsilon }_i v_{\mathrm{C}i}\right)\times {10}^{-6}$$(8)

In previous equations t_Ci, p_Ci and v_Ci represent contributions to the critical temperature, the critical pressure, and the critical molar volume because of the presence of a defined functional group i in a molecule.

One of the major contributions of the characterization procedure described above is the assignment of a hypothetical chemical structure to an undefined petroleum fraction. This allows us to handle the undefined fraction as hypothetical pure component and, in principle, to make use of any thermodynamic model exclusively developed for known components and their mixtures. Thus, in addition to PR-EOS, Huron-Vidal (Michelsen, 1990) mixing rule together with the UNIversal quasi-chemical Functional group Activity Coefficients (UNIFAC) approach, as activity coefficient model, are used also in this work for calculating VLE of hydrogen and heavy petroleum fluids, where the group-interaction parameters Ψ_nm in the UNIFAC approach are given by$${\Psi }_{\mathrm{nm}}=\mathrm{e}\mathrm{xp}-\frac{{\alpha }_{\mathrm{nm}}+{\beta }_{\mathrm{nm}}T+{\chi }_{\mathrm{nm}}{T}^2}{T}$$(9)

The parameters α_nm, β_nm and χ_nm are specific parameters for a couple of functional groups n and m, whose values are obtained from literature (Horstmann et al., 2005). Notice that α_nm ≠ α_mn, β_nm ≠ β_mn, χ_nm ≠ χ_mn. It is worth pointing out that hydrogen is handled as a single functional group.

2 Solubility of Hydrogen in Defined Hydrocarbon Components

Before proceeding with the calculation of VLE of systems composed of hydrogen and heavy petroleum fluids, the hydrogen solubility in defined hydrocarbons was estimated in order to evaluate the thermodynamic models and their corresponding binary parameters described above.

Table 2 shows the 13 binary systems involved in this study, which are formed by nine paraffinic and four aromatic compounds. The range of experimental temperatures and pressures are also included and the comparison of calculated with experimental values are given in Table 3.

As observed, hydrogen solubility estimation is smaller than the experimental values in all cases. The difference is higher as the carbon number of the component increases; the Relative Average Deviation (RAD) increases from 13% for pentane to 75% for hexatetracontane, whereas, in the case of aromatic compounds increases from 32% to 52% for benzene and pyrene, respectively. Because of these results, group-interaction parameters Ψ_nm were calculated again by fitting calculations to the available VLE experimental data in Table 2. This fitting process requires to identify the functional groups forming the defined components. In case of paraffinic compounds, they can be represented by -CH₃ and -CH₂ functional groups, and aromatic compounds by ACH and AC functional groups. Because -CH₃ and AC belongs to the main group 1 (-CH₂) and 3 (ACH) respectively, only interaction parameters between -CH₂ and H₂, ACH and H₂, and ACH and CH₂ are required.

Table 4 shows the original α_nm parameters between main functional groups.

In the fitting process, for simplicity, parameter α_nm is adjusted while β_nm and χ_nm parameters keep their original values. The fitted group-interaction parameters are obtained by a minimization process, where the objective function (F_Obj) is given by$$F_{\mathrm{Obj}}=\frac{\sqrt{\sum _{i=1}^N{\left[P_{{\mathrm{Sat}}_i^{\mathrm{Exp}}}-P_{{\mathrm{Sat}}_i^{\mathrm{Cal}}}({\alpha }_{\mathrm{nm}},{\alpha }_{\mathrm{mn}})\right]}^2}}{N}$$(10)

Here and are the experimental and calculated saturation pressure on the bubble point respectively and N represent the number of data points. The result of this fitting process is a group of new parameters (α_nm, α_mn) for each binary mixtures (H₂-hydrocarbon). Figures 1a and 1b show the trend of α_nm and α_mn parameters for paraffinic and aromatic compounds as function of the carbon number. The correlated parameters take the form of:$${\mathrm{\alpha }}_{{\mathrm{H}}_2-\mathrm{C}{\mathrm{H}}_2}=a_1\ln \left(b_1{\mathrm{C}}_k+c_1\right)+d_1$$(11) $${\mathrm{\alpha }}_{\mathrm{C}{\mathrm{H}}_2-{\mathrm{H}}_2}=a_2\ln \left({\mathrm{C}}_k\right)+b_2$$(12) $${\mathrm{\alpha }}_{{\mathrm{H}}_2-\mathrm{ACH}}=a_3\ln \left({\mathrm{C}}_k\right)+b_3$$(13) $${\mathrm{\alpha }}_{\mathrm{ACH}-{\mathrm{H}}_2}=a_4\ln \left({\mathrm{C}}_k\right)+b_4$$(14)where, C_k is the carbon number of the defined compounds. and are the group-interaction parameters between -CH₂ and H₂ functional groups, and are group-interaction parameters between ACH and H₂ functional groups. Finally, coefficients a_n, b_n, c_n, d_n are shown in Table 5.

Similar results have been obtained by Florusse et al. (2003); they used a Statistical Associating Fluid Theory (SAFT) model to predict hydrogen solubility in heavy alkanes, where the corresponding cross-binary interaction parameters were obtained as a function of the carbon number as well.

Figures 2a and 2b show experimental (♦) and calculated hydrogen solubility data for pyrene and hexatetracontane, respectively. The original group-interaction parameters underestimate hydrogen solubility (-----) whereas the predicted values are in better agreement with experimental data (—).

Table 3 summarizes the RAD of all defined components with the original group-interaction parameters and with those given by Equations from (11) to (14).

3 Solubility of Hydrogen in Undefined Petroleum Fractions

As first step, critical properties and molecular pseudostructures of the seven undefined petroleum fractions studied in this work were determined using the characterization procedure described in Section 1. Table 6 shows density, molecular weight and hydrogen-carbon ratio of the seven undefined fractions used in this work. As will be discussed later, the last two properties play an important role in modeling such systems with hydrogen. In all cases, the experimental error reported in the literature for hydrogen solubility was ±5% (Cai et al., 2001; Ji et al., 2013), except for AB system where such experimental error was not reported (Lal et al., 1999).

Once the molecular pseudostructure is formed with the functional groups from Table 1, it is possible to make a distinction between paraffinic, naphthenic and aromatic contributions through to the no-integer coefficients (υ_j) resulting from the minimization process. Notice that although nine functional groups are used in the characterization procedure (see Tab. 1), there are only two main functional groupS according to the UNIFAC method: first six functional belong to the main -CH₂ functional group, whereas the rest of them belong to the main ACH functional group. Accordingly, the calculated coefficients υ_j can be added in order to obtain the carbon number according to these two main functional groups. Therefore, Equations from (11) to (14) can be used to estimate group-interaction parameters between the H₂ and the functional groups calculated by the characterization procedure.

As illustration, Figure 3 shows the experimental and predicted hydrogen solubility in CHVGO (at 603.15 K) which has a molecular weight of 350 g/mol. The black diamonds (♦) represent the experimental values and the crossed bars represent the experimental reported error of ±5%. Calculations using the original group-interaction parameters (— —) underestimates the hydrogen solubility, as it is also observed for the defined hydrocarbons from Table 2; in this case the RAD is 51%. When our proposed correlations, Equations from (11) to (14), are used to estimate group-interaction parameters, the predicted values (⋅⋅⋅⋅) remains below experimental data, but the RAD is reduced up to 25%. To improve the accuracy of hydrogen solubility prediction, a modification of the original minimization process of characterization is suggested. This modification consists of including another constraint, in addition to Equations (3) and (4). We believed that this constraint should somehow account for the chemical nature of the undefined fractional. Hence, the Hydrogen-Carbon ratio (H/C) is proposed as a third constraint. According to experimental observations made by Ji et al. (2013), this parameter may be related to the aromaticity of undefined fraction.

Figure 4 shows the coefficients (υ_j) obtained from the original characterization procedure and the modified one with the hydrogen-carbon ratio as constraint. The latter gives smaller values for the aromatic functional groups and larger values for all other functional groups. According to Park et al. (1996) hydrogen solubility decreases as the number of aromatic rings increases. Hence, we believe that the inclusion of the H/C ratio may help to construct a pseudostructure in better agreement with its aromatic nature. Figure 3 shows the predicted hydrogen solubility (—) using the modified characterization procedure. As seen, predicted hydrogen solubilities are in better agreement with experimental data (RAD = 5%).

The modified characterization procedure is summarized in Figure 5, where the stages are described in a flowchart. Here (H⁄C)_j is the hydrogen-carbon ratio of a functional group j and (H⁄C)_Exp is the experimental hydrogen-carbon ratio of a i undefined petroleum fraction.

As observed in Table 6, three undefined fractions have molecular weights below 360 g/mol, whereas the rest of them have values above 500 g/mol. Although good results were achieved in the former, hydrogen solubility was overestimated in the latter. When the calculated pseudostructures are reviewed, it is found that they are formed mainly by paraffinic and cyclic functional groups and by a small portion of aromatic ones. This seems to be the reason why the hydrogen solubility is overestimated. Hence, a sensitivity analysis was made to know the effect of the H/C ratio in the characterization process from Figure 5.

In Figure 6a, three different pseudostructures of AB depending on the H/C ratio are shown. When H/C parameter is set equal to the experimental value, the pseudostructure is mainly formed by paraffinic and naphthenic functional groups and hydrogen solubility is overestimated (RAD > 50%). In contrast, hydrogen solubility is underestimated (RAD > 20%) when H/C parameter is set equal to 30% below the experimental value. In latter case, aromatic functional groups predominate over paraffinic and naphthenic ones. Figure 6b shows the trend of the RAD as function of the used H/C ratio. It can be seen that there is minimum value when H/C parameter is set around 80% of the experimental value; its corresponding pseudostructure is shown in Figure 6a. At this point, it is important to emphasize that the characterization procedure do not intend to make a characterization at molecular level, but look for an average structure that allow to reproduce bulk properties (Carreón-Calderón et al., 2012, 2014; Uribe-Vargas et al., 2016).

As with the A/B undefined petroleum fraction, the same sensitivity analysis was made for the VNAR, KRAR, and LAHR fractions. Figure 7 illustrates that minimum values of the RAD are achieved when H/C parameter lies between 80 and 85% of the experimental value. Therefore, we suggested to set the H/C ratio in this interval in order to obtained better results for heavy undefined petroleum fraction with molecular weights above 500 g/mol.

The results of hydrogen solubility prediction in Athabasca Bitumen (AB) at different temperatures and pressures are illustrated in Figures 8a and 8b. Here, the experimental value of the H/C ratio is 1.51, but, according to that previously discussed, it was set equal to 1.208, obtaining an overall RAD of 6.8%. In Table 7 all results are summarized, the overall RAD being equal to 14.12%.

Conclusions

The characterization procedure presented by Carreón-Calderón et al. (2012, 2014) was extended to predict hydrogen solubility in heavy undefined petroleum fractions including H₂ as another functional group, and the hydrogen-carbon ratio as third constraint in the characterization process. Group-interaction parameters between hydrogen and the main functional groups of the UNIFAC approach were calculated through a set of correlations based on experimental vapor-liquid equilibria data of binary mixtures of hydrogen and defined hydrocarbons. These modifications were enough for prediction purposes with molecular weights below 500 g/mol. For undefined petroleum fraction with molecular weights above 500 g/mol, the hydrogen-carbon ratio was modified to 80% of its experimental value. The inclusion of the hydrogen-carbon ratio was revealed as important property in order to construct the hypothetical structures in better accordance with its molecular nature.

Acknowledgments

The authors express their gratitude to the Mexican Institute of Petroleum and SENER-CONACyT for financial support through the Y.61006 research project.

References

All Tables

All Figures