Derivation of Molecular Representations of Middle Distillates

2378

Energy & Fuels 2005, 19, 2378-2393

Derivation of Molecular Representations of Middle Distillates Zhanyao Ha,†,‡ Zbigniew Ring,‡ and Shijie Liu*,†,§,# Department of Chemical and Material Engineering, University of Alberta, Edmonton, Alberta T6G 2G6, Canada, National Center for Upgrading Technology, 1 Oil Patch Drive, Devon, Alberta T9G 1A8, Canada, Alberta Research Council, 250 Karl Clark Road, Edmonton, Alberta T6N 1E4, Canada, and Faculty of Paper Science and Engineering, SUNY College of Environmental Science and Forestry, 1 Forest Drive, Syracuse, New York 13210 Received April 4, 2005. Revised Manuscript Received July 21, 2005

The molecular representation of hydrocarbon mixtures is critical to the advanced kinetics modeling of refining conversion processes; however, the achievement of such a representation is considered a significant challenge. The isomeric lump in a homologous series sets the analytical limit in analytical characterization of middle and heavy distillates. This paper proposes a new procedure for de-lumping detailed analytical information generated using a gas chromatographyfield-ionization mass spectrometry (GC-FIMS) method into a molecular representation. As a result, concentration distributions of the various molecules in the sample of interest are calculated. This paper presents a deterministic computer-assisted procedure that automatically (i) generates hydrocarbon molecules according to literature-based set of rules, (ii) selects the hydrocarbon molecules that are most thermodynamically stable (likely to exist), (iii) optimizes the threedimensional geometry of those molecules and calculates their thermodynamic properties, and (iv) calculates the concentration distributions of those molecules. A separate validation calculation involves (i) prediction of the physical properties for pure hydrocarbons using quantitative structure-property relationship (QSPR) correlations; and (ii) prediction of bulk physical properties of this mixture using the calculated concentrations, properties of its pure hydrocarbon components, and suitable mixing rules. This procedure was applied to find molecular representations of five middle-distillate samples and the results were validated through the estimation and comparison of such simulated and measured physical properties as density, refractive index, and simulated distillation curves (the latter of which are used as a consistency check). Good agreement was observed between the predicted and measured properties for all five middle-distillate samples. This agreement validates the molecular characterization algorithm, at least for the purpose of bulk property prediction in middle distillates.

1. Introduction Accurate characterization of feedstocks has an important role in kinetics modeling and integrated process optimizations. In refining conversion processes such as hydrocracking and fluid catalytic cracking (FCC), the reactivity and product yield distributions are strongly dependent on the molecular composition of feedstocks. Feedstocks of different molecular compositions have different cracking patterns. Advances in the characterization of petroleum fractions are closely related to the development of new more-sensitive and less-expensive analytical techniques. The composition of light petroleum fractions (i.e., naphtha) can be measured at the molecular level using the Detailed Hydrocarbon Analysis (DHA) method (ASTM D6733), where unambiguous * Author to whom correspondence should be addressed. E-mail address: [email protected]. † Department of Chemical and Material Engineering, University of Alberta. ‡ National Center for Upgrading Technology. § Alberta Research Council. # Faculty of Paper Science and Engineering, SUNY College of Environmental Science and Forestry.

structural assignments can go as far as C9 isoparaffins, C7 cycloparaffins, and C10 aromatics. However, the compositions of middle distillates and vacuum gas oils can only be measured, in terms of hydrocarbon class distributions, by boiling point (BP)1 or carbon number (CN).2,3 One of the analytical techniques that is particularly useful for middle and heavy distillates is gas chromatography-field-ionization mass spectrometry (GC-FIMS). Because of the formidable chemical complexity of petroleum fractions, the state-of-the-art mechanistic kinetics models of refining conversion processes, such as the Structural Oriented Lumping (SOL) model,4 the Linear Free Energy Relationship (LEFR) model,5 and the Single-Event Kinetics model,6 have resorted to (1) Boduszynski, M. M. Energy Fuels 1988, 2 (5), 597-613. (2) Bouquet, M.; Brument, J. Fuel Sci. Technol. Int. 1990, 8 (9), 961986. (3) Chasey, K. L.; Aczel, T. Energy Fuels 1991, 5, 386-394. (4) Quann, R. J.; Jaffe, S. B. Ind. Eng. Chem. Res. 1992, 31, 24832497. (5) Broadbelt, L. J.; Stark, S. M.; Klein, M. T. Ind. Eng. Chem. Res. 1994, 33, 790-799. (6) Hillewaert, L. P.; Dierickx, J. L.; Froment, G. F. AIChE J. 1988, 34, 17-24.

10.1021/ef050097g CCC: $30.25 © 2005 American Chemical Society Published on Web 09/10/2005

Molecular Representations of Middle Distillates

various types of lumping to reduce the computational burden associated with solving those models. A strong assumption in lumping isomers for the purpose of product property prediction is that the physical properties of the hydrocarbon isomers within a lump are identical. This assumption is not true for most hydrocarbons. Many thermophysical properties of isomers vary widely with their molecular structures.7 The catalytic reaction mechanisms, reactivity, and selectivity can also be quite different among isomers (particularly when shape-selective catalysts are involved8). For example, the cracking of three C6 isomers (2-methylpentane, 3-methylpentane, and 2,3-dimethylbutane) on USHY shows that these three undertake very different reaction paths and yield different cracked-product distributions.9 Furthermore, in mechanistic kinetics models, the reactivity is often correlated with reaction indices that are dependent on the structure of compounds. Therefore, it is important to determine detailed distributions of isomers within an isomeric lump. Neurock et al.10 developed a Monte Carlo method to characterize complex mixtures. In their method, bulk elemental analysis and nuclear magnetic resonance (NMR) measurements were used to build the probability curves for various structural attributes of hydrocarbon molecules, such as side chains, aromatic rings, and naphthenic rings. The Monte Carlo method was then used to randomly sample the probability distributions and construct the molecules, thereby forming a set of molecules to fit the analytical bulk data. Molecules (on the order of 105) were sampled to represent the petroleum mixtures, including heavy oil, vacuum gas oil, and asphaltenes. Further efforts had been made to derive a smaller set (10-100) of molecules based on the Monte Carlo representations using the quadrature method.11 Random molecule generation suffers from an inefficiency that is associated with frequent generation of redundant molecules and occurs if due consideration is not given to their thermodynamic stability; it may lead to molecules that are unlikely to exist. However, to date, probably the most significant weakness of the work on deriving molecular representations is that a relatively small amount of analytical information is used as its basis. This results in multiple nonunique molecular representations and makes it difficult to test and independently verify the resulting molecular representations; as a result, this process introduces significant uncertainty. Molecular composition is only partially reflected by bulk properties such as gravity, elemental content, and NMR analysis. Kim et al.12 investigated the compositions of three FCC feedstocks that had very similar bulk properties. Through mass spectrometry, they found that the hydrocarbon class distributions in those feedstocks were significantly different. Using high(7) API Technical Data Book; American Petroleum Institute: Washington, DC, 1992. (8) Krambeck, F. J. Kinetics and Thermodynamic Lumping of Multicomponent Mixtures; Astarita, G., Sandier, S. I., Eds.; Elsevier Science Publishers BV: Amsterdam, 1991; pp 111-129. (9) Wojciechowski, B. W. Catal. Rev.sSci. Eng. 1998, 40 (3), 209328. (10) Neurock, M. N.; Nigam, A.; Trauth, D.; Klein, M. T. Chem. Eng. Sci. 1994, 49, 4153-4177. (11) Campbell, D. M.; Klein, M. T. Appl. Catal., A 1997, 160, 4154. (12) Kim, H. N.; Verstraete, J. P.; Virk, P. S.; Fafet, Oil Gas J. 1998, 96, 85-88.

Energy & Fuels, Vol. 19, No. 6, 2005 2379

Figure 1. Molecular de-lumping algorithm for molecular characterization.

ionizing-voltage mass spectrometry (ASTM D2786/ D3239), Ramaswamy et al.13 also found that the hydrocarbon class distributions of two vacuum gas oils were entirely different, although the two samples had almost the same contents of saturates, aromatics, and polars (ASTM D2549). In fairness, it is usually difficult, and sometimes impossible, to get analytical information that goes beyond bulk analysis and BP distribution, particularly for heavy materials such as residues. However, it is possible to get detailed characterization data for middle and heavy distillates (e.g., high-resolution mass spectrometry or gas chromatography with low- or medium-resolution spectrometry). In principle, accurate molecular representations should be based on as much information as possible to limit the number of possible solutions and/or assumptions made to derive molecular representations. Our work, which is described below, used distributions of hydrocarbon class versus carbon number, which were determined using GC-FIMS as the starting point, and developed a deterministic procedure that generated a set of molecules using molecular construction rules and a thermodynamic stability criterion, and then this procedure was used to determine the concentrations of all individual isomers. The procedure was validated using five compositionally dissimilar middle-distillate samples through the estimation of selected physical properties of the corresponding five molecular mixtures and comparing those properties with independently generated measurements. 2. Methodology The objective of this work was to develop a deterministic method for molecular characterization of petroleum middle-distillate diesel blending components. The resulting method is briefly described below for middle distillates as an example. Starting from a reconciled “characterization matrix” that was measured by GCFIMS, which determines 19 hydrocarbon homologous series distributed by carbon number (CN), a delumping algorithm is used to derive a molecular representation (see Figure 1). For each measured isomeric lump in the distillate sample, a group of isomers is generated, using a set of molecular generation rules that are derived from information (scattered throughout the open literature) on molecular substitution patterns in various hydro(13) Ramaswamy, V.; Singh, I. D.; Krishna, R. Indian J. Technol. 1989, 27, 85-88.

2380

Energy & Fuels, Vol. 19, No. 6, 2005

Ha et al.

Table 1. Summary of the Bulk Properties and Main Components for Five Middle-Distillate Samples Saturates sample

total

nP/IPa

S1 S2 S3 S4 S5

76.98 73.52 79.51 74.17 54.74

46.56 33.60 5.58 8.48 18.62

a

Aromatics CycloP

total

1As

2As

3As

Arom-S

IBP/FBP (°C)

refractive index, RI

density

MW

30.42 39.92 73.92 65.69 36.12

23.02 26.48 20.49 25.83 45.26

22.24 20.83 19.28 24.41 33.98

0.74 5.58 1.19 1.37 11.18

0.00 0.01 0.00 0.00 0.00

0.04 0.06 0.02 0.05 0.10

97/277.3 120.9/344.5 118/334.1 114.8/338 151.6/299.5

1.4436 1.4593 1.4584 1.4652 1.4755

0.7985 0.8291 0.8357 0.8476 0.8537

150 188 176 176 175

The terms nP and IP are defined as follows: nP ) n-paraffins and IP ) iso-paraffins.

carbons. The generated molecules are optimized geometrically on-line, using the MOPAC 2002 PM3 procedure, and their Gibbs free energies of formation also are calculated using MOPAC. A subset of isomers to represent the isomeric lump is selected, using a preset Gibbs-free-energy-of-formation threshold, which is used as an adjustable parameter of the method to control the number of molecules in the representation. To improve the computational efficiency, the heat of formation (the dominant part of the Gibbs free energy of formation for most hydrocarbons) is used as the threshold instead of the free energy in our procedure.5 With the list of representative isomers at hand, the distribution of isomer concentrations within the isomeric lump can be calculated by minimizing the Gibbs free energy of the lump, subject to the stoichiometric constraint and another constraint in the form of the BP distribution of the lump.14 The latter is a reflection of independently measured partial knowledge of the lump composition. At this point, a molecular representation of each isomeric lump is available that contains a number (adjustable) of specific, most thermodynamically favorable (i.e., most abundant in the material of interest) molecules of known concentrations. The molecular representation of the middle distillate is the collection of the molecular representations of all isomeric lumps put together. To assess how good this molecular representation is, in terms of representing its physical properties, the following procedure was used. For each molecule in the molecular representation, the specific gravity (SG) and refractive index (RI) were estimated using quantitative structure-property relationship (QSPR) correlations that have been developed by Ha et al.15 Using the calculated concentrations of individual molecules in the representation, their estimated properties, and simple mixing rules, bulk physical properties of the molecule mixtures were calculated and compared with the corresponding measured bulk properties. The procedure outlined above consists of several tasks, which include (i) analytical measurements and data reconciliation, (ii) molecule generation of isomers and on-line simulation, (iii) QSPR predictions for physical properties and thermodynamic calculations, (iv) isomer distribution within each measured isomeric lump, and (v) molecular refining to derive the minimum set of molecules. 2.1. Analytical Measurements and Data Reconciliation. Five compositionally diverse middle distillates were analyzed to validate our molecular-representation method. Samples S1-S3 were obtained from an

Edmonton refinery: Sample S1 was a heavy naphtha (boiling at 130-230 °C) with a high paraffinic content, sample S2 was a highly paraffinic light distillate (boiling at 200-300 °C), and sample S3 was a highly naphthenic light distillate. Sample S4 was a hydrotreated middle distillate derived from Canadian bitumen with a high naphthenic content. Sample S5 was a highly aromatic light distillate from an offshore crude. These samples were analyzed to determine the average molecular weight (MW), SG, RI, simulated distillation (SimDis), and hydrocarbon-versus-CN and hydrocarbon-versusBP distributions that were determined using GCFIMS. The bulk properties and the concentrations of the most important hydrocarbon components of all five samples are summarized in Table 1. All the analytical measurements were conducted at the National Centre for Upgrading Technology (NCUT). The density at 15.6 °C was measured via ASTM method D4052,16 using a DMA4 PAAR densitometer, and that value was used to calculate the corresponding SG value. The RI value at 25 °C was measured via ASTM method D1218,16 using an Abbe refractometer. The MW value was determined using the freezing point depression (FPD) method. SimDis results were obtained using the ASTM D288716 procedure and a Hewlett-Packard model HP 6890 gas chromatograph. The method for hydrocarbon class analysis by GC-FIMS was developed at NCUT.17 Tables 2 and 3 show the GC-FIMS reports (mass fractions of various lumps) for sample 1 as determined using the CN distribution and the BP distribution, respectively. Because tetraaromatics and benzonaphthothiophenes were not detected in the five middledistillate samples, only 16 hydrocarbon classes are reported in these tables. The two-dimensional (2D) mass fraction-versus-BP distributions of individual hydrocarbon types, as measured by GC-FIMS, had to be made consistent with the corresponding SimDis data (total mass distribution as determined by boiling point). The required consistency between the two distributions was achieved in terms of the both the retention time (RT) and mass scales. The RT-BP relationship of most hydrocarbons in petroleum middle distillates is very similar to that of n-paraffins (the experimental basis of the ASTM D2887 SimDis method for petroleum fractions16). The consistency of the RT scale was achieved mathematically by matching the retention times of the n-paraffinic calibration standard analyzed on the two respective instruments and converting the two RT scales into one BP scale. The SimDis test relies on a flame ionization detector (FID), where the response factors for all hydrocarbons are

(14) Ha, Z.; Ring, Z.; Liu, S. Estimation of isomeric distributions in petroleum fractions. Energy Fuels 2005, 19, 1660-1672. (15) Ha, Z.; Ring, Z.; Liu, S. Energy Fuels 2005, 19, 152-163.

(16) 2001 Annual Book of ASTM Standards; American Society for Testing and Materials, West Conshohocken, PA, 2001. (17) Briker, Y.; Ring, Z.; Iacchelli, A.; McLean, N.; Rahimi, M.; Fairbridge, C. Energy Fuels 2001, 15, 23-37.



Table 2. Gas Chromatography-Field Ionization Mass Spectroscopy (GC-FIMS) Results, Relative to Carbon Number (CN) Distribution: Sample 1 hydrocarbon class

C6

C7

C8

C9

C10

C11

C12

C13

C14

C15

sum

saturates total paraffins isoparaffins n-paraffins total cycloparaffins monocycloparaffins dicycloparaffins polycycloparaffins aromatics total monoaromatics alkylbenzenes benzocycloalkanes benzodicycloalkanes total diaromatics naphthalenes biphenyls naphthocycloalkanes fluorenes total triaromatics phenanthrenes phenanthrocyclolalkanes total aromatic sulfur benzothiophenes dibenzothiophenes

0.16 0.06 0.02 0.03 0.10 0.10 0.00 0.00 0.01 0.01 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

1.17 0.45 0.17 0.28 0.72 0.50 0.22 0.00 0.25 0.25 0.25 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

3.14 1.52 0.74 0.78 1.61 1.12 0.49 0.00 1.74 1.74 1.74 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

8.74 4.80 2.05 2.75 3.95 2.75 1.18 0.01 8.10 8.10 6.89 1.20 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

19.18 13.93 8.39 5.54 5.26 3.52 1.69 0.04 5.59 5.48 4.17 1.31 0.00 0.11 0.11 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

16.30 11.02 6.20 4.82 5.28 3.27 1.81 0.20 3.14 2.80 1.60 1.18 0.01 0.30 0.30 0.00 0.00 0.00 0.00 0.00 0.00 0.04 0.04 0.00

12.78 6.80 4.45 2.35 5.98 3.21 2.36 0.40 2.66 2.37 0.94 1.39 0.03 0.29 0.22 0.07 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

10.03 4.93 3.43 1.50 5.09 3.01 1.66 0.43 1.21 1.18 0.53 0.60 0.04 0.04 0.01 0.03 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

4.87 2.64 2.27 0.37 2.23 1.37 0.65 0.21 0.31 0.31 0.15 0.14 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

0.60 0.40 0.39 0.02 0.20 0.12 0.06 0.01 0.02 0.02 0.01 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

76.98 46.56 28.10 18.46 30.42 18.98 10.13 1.30 23.02 22.24 16.30 5.85 0.09 0.74 0.64 0.10 0.00 0.00 0.00 0.00 0.00 0.04 0.04 0.00

approximately the same. In contrast, the GC-FIMS analysis relies on a MS detector, where those factors vary significantly with the hydrocarbon classsand within a hydrocarbon classswith CN.17 Therefore, the two raw BP distributions are still not directly comparable until the respective response factors are applied. Such a treatment of raw GC-FIMS data results in a by-BP-distributed characterization matrix, such as that shown in Table 3. Summing each column in the reconciled characterization matrix results in a mass distribution equivalent to SimDis (FIMS-generated SimDis). The actual SimDis and FIMS-generated SimDis curves for sample 1, which are compared in Figure 2, show good agreement. These data reconciliation procedures were applied to all the samples, and full internal consistency was demonstrated (see Figures 2-6). 2.2. Molecular Generation and Simulation. 2.2.1. Molecular Generation Rules. In previous studies on the molecular representations of hydrocarbon mixtures, the molecules used for those representations were randomly generated from structural building blocks.10 A large set of those molecules then was optimized to match the derived distribution functions of structural attributes until a set of bulk physical properties of the mixture was sufficiently close to the measured values. In this work, we explored a deterministic way of generating molecules to improve the computational efficiency of the method. The way to build molecules was based on reported occurrences in petroleum and quantum-chemical simulations (for those molecules that have not been isolated and quantified). Each cell of the by-CN characterization matrix was assumed to be an isomeric lump of a specific hydrocarbon class at a given CN. The number of possible isomers within the lump increases rapidly with the number of C atoms contained. Even for light distillates (CN < 21), the number of all possible isomers for a specific CN can be in the millions.18 However, not all isomers are equally probable to be (18) Read, R. C. The enumeration of acyclic chemical compounds. In Chemical Application of Graph Theory; Balaban, A. T., Ed.; Academic Press: New York, 1976.

present in the isomeric lump in significant concentrations. The occurrences of individual isomers may be dependent on their thermodynamic stability,19 their origin, and the maturing history of the crude. For alkanes, the 2- and 3-methyl derivatives are the most abundant, and the 4-methyl derivative is present in small amounts. It is generally accepted that moderately branched paraffins dominate over highly branched paraffins.20 Tissot and Welte19 conclude that the most frequent configuration of paraffins is one tertiary C atom (2-methyl or 3-methyl). The configuration with two tertiary C atoms is less abundant and that with one quaternary C atom is usually very rare. Small alkanes are not inherently biogenic molecules. They are generated through thermal degradation and cracking of C-C bonds of either kerogen or larger bitumen molecules. The compositional investigations of 18 crude oils by Martin et al.21 also showed the prevailing occurrences of n-paraffins and 2- to 3-methyl iso-paraffins in naphthas of those crudes. Apparently, the nature favors simple paraffinic structures, at least in the middledistillate range. However, contrary to this observation, the thermodynamic data (both heat of formation and free energy of formation) of paraffin isomers in the naphtha range show that the 2,2-dimethyl branched paraffins are thermodynamically most stable, whereas n-paraffins and other 2-methyl and 3-methyl branched paraffins are less stable. Figure 7 illustrates the reported free energies of formations for C9 isoparaffins.7 A possible explanation for this apparent inconsistency is that the concentrations of less-stable alkanes were “frozen” at relatively high levels, because of the kinetics barriers, after they evolved due to thermal cracking of larger crude molecules. NMR analyses for crude oils and their distillates showed that the tertiary aliphatic (19) Tissot, B. P.; Welte, D. H. Petroleum Formation and Occurrence; Springer-Verlag: Berlin, 1984. (20) Speight, J. G. The Chemistry and Technology of Petroleum, 3rd Edition; Marcel Dekker: New York, 1999. (21) Martin, R. L.; Winters, J. C.; Williams, J. A. In Proceedings of the 6th World Petroleum Congress, Frankfurt, Germany, 1963, pp 231260.

2382


Ha et al.

Table 3. Gas Chromatography-Field Ionization Mass Spectroscopy (GC-FIMS) Results, Relative to BP Distribution: Sample 1a hydrocarbon class

IBP-100

100-110

110-120

120-130

130-140

140-150

150-160

160-170

170-180

180-190

saturates total paraffins isoparaffins n-paraffins total cycloparaffins MonocycloP DicycloP PolycycloP aromatics total monoaromatics alkylbenzenes BenzocycloP BenzodicycloP total diaromatics naphthalenes biphenyls NaphthocycloP fluorenes total triaromatics phenanthrenes phenanthrocycloP total aromatic sulfur benzothiophenes dibenzothiophenes

0.36 0.31 0.08 0.23 0.05 0.05 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

0.69 0.34 0.06 0.28 0.35 0.35 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

0.53 0.48 0.48 0.00 0.05 0.05 0.00 0.00 0.25 0.25 0.25 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

1.26 0.79 0.11 0.67 0.48 0.48 0.00 0.00 0.07 0.07 0.07 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

1.95 0.79 0.76 0.02 1.16 1.15 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

1.33 0.60 0.60 0.00 0.73 0.70 0.02 0.01 1.22 1.22 1.22 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

5.62 3.51 0.81 2.71 2.11 1.61 0.49 0.01 0.90 0.90 0.90 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

7.35 4.81 4.81 0.00 2.54 1.85 0.68 0.01 2.41 2.41 2.41 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

10.44 6.84 1.39 5.45 3.60 2.66 0.92 0.02 3.45 3.45 3.44 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

7.51 5.07 5.07 0.00 2.44 1.39 1.02 0.03 1.70 1.70 1.63 0.07 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

hydrocarbon class

190-200

200-210

210-220

220-230

230-240

240-250

250-260

260-270

270-280

sum

saturates total paraffins isoparaffins n-paraffins total cycloparaffins MonocycloP DicycloP PolycycloP aromatics total monoaromatics alkylbenzenes BenzocycloP BenzodicycloP total diaromatics naphthalenes biphenyls NaphthocycloP fluorenes total triaromatics phenanthrenes phenanthrocycloP total aromatic sulfur benzothiophenes dibenzothiophenes

9.33 5.43 1.85 3.58 3.90 2.67 1.11 0.12 2.00 1.99 1.65 0.34 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

6.21 3.34 3.34 0.01 2.86 1.19 1.57 0.11 2.04 2.04 1.75 0.29 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

7.87 4.78 1.86 2.92 3.09 1.43 1.47 0.18 2.47 2.36 1.40 0.96 0.00 0.11 0.11 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

5.31 2.56 2.54 0.01 2.75 1.42 1.12 0.21 1.36 1.36 0.66 0.70 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

7.16 4.84 2.70 2.14 2.32 1.15 0.91 0.26 2.10 1.86 0.44 1.41 0.01 0.20 0.19 0.01 0.00 0.00 0.00 0.00 0.00 0.04 0.04 0.00

2.32 1.11 1.11 0.00 1.22 0.52 0.51 0.19 1.60 1.48 0.28 1.18 0.03 0.12 0.11 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

1.50 0.92 0.50 0.42 0.58 0.25 0.22 0.11 0.97 0.79 0.15 0.61 0.03 0.18 0.10 0.08 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

0.20 0.03 0.03 0.01 0.17 0.05 0.07 0.05 0.46 0.33 0.05 0.26 0.03 0.13 0.12 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

0.02 0.01 0.00 0.01 0.02 0.00 0.01 0.01 0.03 0.02 0.00 0.02 0.00 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

76.98 46.56 28.10 18.46 30.42 18.98 10.13 1.30 23.02 22.24 16.30 5.85 0.09 0.74 0.64 0.10 0.00 0.00 0.00 0.00 0.00 0.04 0.04 0.00

a

Boiling point (BP) ranges are expressed in units of degrees Celsius (°C).

carbon accounts for ∼10% of the total carbon and the branch index (CH3/CH2) of saturates ranged from 0.1 to 0.3.13,22,23 We limited the isoparaffins in the middle distillate to have, at most, three branches with a maximum branch length of C4, and no quaternary carbon was allowed. Similar rules were made by Mizan and Klein24 in their study of catalytic hydroisomerization, to rationalize the product yield distributions. Based on the fact that only mono-cycloparaffins with five- and six-membered rings have been isolated from light petroleum fractions and that thermodynamic studies show that naphthenic rings with five and six C atoms are thermodynamically most stable,20 it can be assumed (22) Japanwala, S.; Chung, K. H.; Dettman, H. D.; Gray, M. R. Energy Fuels 2002, 16, 477-484. (23) Khan, H. U.; Sharma, R. L.; Nautiyal, S. P.; Agrawal, K. M.; Schmidt, P. Pet. Sci. Technol. 2000, 18, 889-899. (24) Mizan, T. I.; Klein, M. T. Catal. Today 1999, 50, 159-172.

that only these mono-cycloparaffins are present in appreciable concentrations. Methyl-cyclohexanes and methyl-cyclopentanes are frequently the most abundant in the light fractions (CN < 10). Monocycloalkanes and dicycloalkanes generally account for half of the total cycloparaffins for CN > 10. Tricycloalkanes are most likely to have a structural arrangement similar to perhydrophenanthrene. The abundance of tricycloalkanes decreases with CN beyond 20. The structures of tetracycloalkanes and pentacycloalkanes are directly related to the tetracyclic steroids and pentacyclic triterpenes.19 Generally, tetracycloalkanes and pentacycloalkanes are most abundant in immature crude oils. They are beyond the range of middle distillates and are reported in gas oils. The aromatics in crude oils usually include 1-5 condensed aromatic rings substituted with a small number of short paraffinic chains. Among them, alky-



Figure 2. “GC-FIMS-generated SimDis” and “Simulated SimDis”, compared to experimental SimDis for sample 1.


lated homologues of benzene, naphthalene, and phenanthrene are the most abundant.19 They often include 1-3 C atoms in substituent chains. For instance, the most abundant alkyl benzene in the gasoline fraction is frequently toluene or xylene, whereas benzene is usually less abundant. The same is true for naphthalene and phenanthrene, where ethyl- or propyl-substituted molecules are the most abundant. When several structural arrangements are possible, as for molecules with three or more aromatic rings, only a few of them are favored in crude oils. For example, alkyl-phenanthrenes are largely dominant triaromatics. Among the five isomeric configurations of tetra-aromatics, cata-condensed chrysene is likely to be the most prevalent form, based on the origin of ring structures in petroleum.19 Naphthenoaromatics are particularly abundant, compared to purely aromatic structures, in shallow immature crude oils. The aromatic types become dominant after a

significant thermal evolution. Bicyclic indane, tetralin, and their methyl derivatives are usually abundant. Tricyclic tetrahydrophenanthrene and derivatives are also quite common. Tetracyclic and pentacyclic molecules, which are mostly related to steroid and triterpenoid structures with 1-3 aromatic rings, are important in the high boiling fractions, such as vacuum gas oils.19 The distribution of alkyl groups may vary with the number of groups on a ring core, the carbon number, and the degree of branching of the substituents. With the NMR measurements on a physically isolated monoaromatic fraction of petroleum, Mair and Barnewall25 deduced that, on average, alkyl benzenes in the C13C15 range had one methyl substituent and one longer (25) Mair, B. J.; Barnewall, J. M. J. Chem. Eng. Data 1964, 9 (2), 282-292.

2384


Ha et al.



chain with a methyl branch. By comparing the MS fragmentation pattern of both synthetic and petroleum alkylbenzenes in the C20-C40 range, Hood et al.26 concluded that the alkyl substituent structures are not as diverse as suggested by the number of possible isomers. They observed that petroleum alkylbenzenes are composed almost entirely of those that have a single long substituent chain and 1-4 methyl substituents on the benzene ring. It is interesting that the distribution of substituents differs between aromatics and naphthenic rings within the same molecule. The short chains (methyl and ethyl) seem to be prevalent substituents of the aromatic portion of the ring cluster, whereas a limited number (1 or 2) of longer chains may be attached to the naphthenic rings. The total number of chains, (26) Hood, A.; Clere, R. J.; O’Neal, M. J. J. Inst. Pet. 1959, 45 (426), 168-173.

which is generally 4-6, as well as their lengths, increases with CN.26 A molecular simulation study of the thermodynamic stabilities of various aromatic isomers also showed a tendency for a higher degree of substitution in alkylated aromatics to be more stable. Figure 8 shows the thermodynamic stability ranking for C10 alkyl-benzenes. The figure indicates that the simulated values of the standard heats of formation differ from the measured values. However, these differences do not affect the stability ranking, which is the same in both cases. Therefore, the simulated results can be used for stability ranking. Figures 7 and 8 show that, with the same number of branches, the molecules with wellseparated substituents are more thermodynamically stable than those with closely located substituents. As shown in Figures 7 and 8, 2,6-dimethyl-heptane and 1,3-




Figure 7. Thermodynamic stability rankings, based on standard free energies of formation for C9 isoparaffins. (Notations: M, methyl; E, ethyl.)

dimethyl-5-ethyl-benzene (2,2-dimethyl-heptane excluded) are the most stable structures for dimethylheptanes and tribranched C10 benzenes, respectively. Based on literature data and the previously mentioned observations, the following molecular generation rules were applied in this work: (a) There are a maximum of three branches in isoparaffinic structures. (b) The maximum length of an isoparaffin branch is CN ) 4. (c) Quaternary carbon is not allowed in saturate structures.

(d) There are a maximum of four substituents on an aromatic/naphthenic ring core. (e) Only one nonmethyl branch is allowed in the case of multiple substitutions on a ring core. (f) Multiple substituents are located on the main chain or ring structure as far apart as possible. (g) Only five- or six-membered naphthtenic rings are allowed. (h) Only one nonmethyl branch is allowed on the naphthenic rings of naphthenoaromatics. 2.2.2. Molecular Fingerprint in MOPAC Input File Format. Following the molecular generation rules, hydrocarbon molecules are constructed in the MOPAC

2386


Ha et al.

Figure 8. Thermodynamic Stability Ranking Based on Standard Heats of Formation for C10 alkyl-benzenes. (Notations: M, methyl; E, ethyl; B, butyl.)

Figure 9. MOPAC input file for ethane.

input file format. Figure 9 shows an example of the MOPAC input file for ethane. The first line includes the functional key words, from which the calculation task is loaded. Lines 2 and 3 are comments. Line 4 and the following lines specify individual atoms in the molecule, the bond length, bond angle, and the dihedral angle among atoms, with the optimization flag (0 or 1) following behind. The MOPAC input file offers three advantages in the course of our molecular characterization algorithm. It (i) contains a detailed three-dimensional description of the molecular structure, (ii) enables viewing (using CS Chem3D, etc.) of individual molecules, and (iii) allows on-line molecular simulation. Broadbelt et al.27 proposed an algorithm for translating a 2D molecular graph into a three-dimensional (3D) internal coordinate representation of the MOPAC input. A similar algorithm was used in our work; however, molecules were directly assembled from the available main structures (chain or ring core) and substituents. The initial bond lengths were assigned based on the types of atoms and types of bonds between them. The bond angle is dependent on (27) Broadbelt, L. J.; Stark, S. M.; Klein, M. T. Chem. Eng. Sci. 1994, 49, 4991-5010.

the hybridization of the parent atom. The specification of the dihedral angle is related to the types of parent, angle, and dihedral atoms and the order of the atoms being visited.27 2.2.3. Molecular Simulation, Selection, and Thermodynamic Calculation. Semiempirical quantum-chemical calculations (AM1 and PM3 procedures) provide a reasonable balance between the computational efficiency and reliability of the estimates for a wide range of electronic and thermodynamic properties. Using the PM3 procedure, Stewart28 reported that the accuracy of the heats of formation for all organic compounds was +0.21 kcal/mol. Therefore, the PM3 procedure was used for geometry optimization in this work. By definition, ∆H°f is the ideal-gas-phase heat of formation at 298 K of one mole of a compound from its elements in their standard state. The number of molecules generated from the rules listed in the previous section is very large. For instance, hundreds of isomers can be built for C20 cyclic hydrocarbons. To limit the number of molecules to a manageable level, we developed a selection algorithm to screen out those less thermodynamically stable (as determined by the Gibbs free energy, approximated by the standard heat of formation, and calculated after the PM3 spatial optimization). The molecule with the lowest heat of formation was taken as a reference. The molecules, whose standard heat of formation differed from the reference molecule by less than a user-selected threshold value, were selected to the molecular representation of the isomeric lump at hand. The threshold, 2.0 kcal/mol in this work, is an adjustable parameter that (28) Stewart, J. J. P. Semiempirical molecular orbital methods. In Reviews in Computational Chemistry; Lipkowitz, D. B., Boyd, D. B., Eds.; VCH Publishers: New York, 1990; Chapter 2.



controls the total number of molecules in the representation. The molecules whose heats of formation differed from the reference molecule by more than 2.0 kcal/mol were discarded. Thermodynamically, a difference of 2.0 kcal/mol in the free energy or heat of formation between two isomers yields the 94:6 (concentration ratio) of one over the other.29 For all the molecules selected to represent the isomeric lump, thermodynamic properties such as enthalpy, heat capacity, and entropy were calculated using MOPAC 2002. From the vibrational, rotational, internal, and translational contributions, MOPAC calculated ideal gas heats of formation and entropies at a prescribed temperature. These, in turn, were used to calculate the standard Gibbs free energies of formation.30

[

∆G°f ) ∆H°f - T So - nCS°C -

( ) ] nH S° 2 H2

(1)

where nC and nH are the number of C and H atoms, respectively, S°C is the entropy of graphite at 298.15 K (S°C ) 1.361 cal K-1 mol-1), and S°H2 the entropy of hydrogen gas at 298.15 K (S°H2 ) 31.211 cal K-1 mol-1). 2.2.4. Generation and Manipulation of Hydrocarbon Molecules. A molecular generation algorithm was created to construct hydrocarbon molecules for all the hydrocarbon classes identified by GC-FIMS (see Table 2). The flow diagram of this algorithm, shown in Figure 10, indicates that the hydrocarbon types of interest ranged from paraffins (hydrogen deficiency Z ) 2) to tetra-aromatics (Z ) -24), with Z decreasing in evennumbered increments, corresponding to individual hydrocarbon types, over the boiling range of middle distillates (CN ) 5-21). As mentioned previously, including three aromatic sulfur types (benzothiophenes, dibenzophiophenes, and benzonaphthothiophenes), 19 type classes (HCType) were identified by the GC-FIMS method used. Because the generation of molecular representations of sulfur-containing compounds was beyond the scope of this work, only 16 purely hydrocarbon classes are considered below. The names of those classes are abbreviated here as follows: nP (normal paraffin), iP (iso-paraffin), 1N (mono-cycloalkane), 2N (dicycloalkane), 3N (tricycloalkane), 1A (alkylbenzene), 1A1N (benzo-cycloalkane), 1A2N (benzo-dicycloalkane), 2A (naphthalene), A_A (biphenyl), 2A1N (naphthocycloalkane), ANA (fluorene), 3A (phenanthrene), 3A1N (phenanthro-cycloalkane), 4Ap (pyrene), and 4Ac (chrysene). First, the isomeric lump is selected for the first hydrocarbon type and the minimum CN within this hydrocarbon type that matches the boiling range of the sample at hand. Following the molecular generation algorithm in Figure 10, a permutation subroutine is called to generate all the possible isomers for a selected Z-CN isomeric lump, based on the molecular building blocks and generation rules. The generated isomers are stored in an isomer identification matrix (IsomerID) from which the corresponding MOPAC input file is (29) Wade, L. G. Organic Chemistry; Prentice-Hall: Englewood Cliffs, NJ, 1987; p 132. (30) Stull, D. R.; Westrum, E. F., Jr.; Sinke, G. C. The Chemical Thermodynamics of Organic Compounds; Wiley: New York, 1969; p 205.

Figure 10. Programming diagram for molecular generation.

generated by calling subroutines (described below) corresponding to individual hydrocarbon types. For each constructed molecule, the MOPAC executable is called on-line to optimize its 3D geometry using the PM3 procedure. The original MOPAC input file is then replaced with the geometry-optimized file, to be used for thermodynamic calculations later. The PM3 optimization procedure also calculates the standard heat of formation (∆H°f) for the optimized molecule. The ∆H°f value for each molecule is compared with the isomer that has the smallest ∆H°f. If the difference is larger than the user-provided threshold, this molecule is discarded. Otherwise, it is accepted to represent the isomeric lump at hand. This prescreening procedure eliminates molecules that are thermodynamically unfavorable, thus improving the computational efficiency of the thermodynamic calculations to follow. The molecular generation proceeds until the last isomer in IsomerID is generated. The program then moves to the next isomeric lump (CN) until the final CN is reached. At this point, the molecular generation for a selected hydrocarbon type (Z series) is finished. The program then selects the next hydrocarbon type and repeats the aforementioned molecular generation and pre-screening procedures until the last Z series is processed. Finally, all the representative molecules are selected. These molecules are stored in a molecular bank for further manipulation.

2388


Ha et al.

Figure 11. Programming diagram for iso-paraffin generation. Table 4. Branch Types Attached to the Main Structure (Chain or Ring) CN in branch

type 1

1 2 3 4 m

methyl ethyl n-propyl n-butyl n-Cm-a

type 2

type 3

type 4

/

/

/

iso-propyl iso-butyl iso-Cm-a

sec-butyl

ter-butyl

a The terms “n-C -” and “iso-C -” respectively represent linear m m and isobranched-linear substituents with m C atoms.

2.2.4.1. Molecular Generations of Paraffins. Generating paraffins is different from building cyclic hydrocarbon molecules. Writing a MOPAC input file for an n-paraffin molecule (HCType ) 1) is relatively straightforward. When the molecular size (CN) is known, a single n-paraffin molecule can be easily constructed by adding C atoms one after another with the corresponding H atoms attached. The MOPAC input file for ethane is shown in Figure 9 as an example. A flow diagram of the algorithm for building isoparaffin molecules is shown in Figure 11. At a given CN, isoparaffin isomers (HCType ) 2) are assigned a main straight chain and several branches and then are assembled by adding branches on the straight chain. The MOPAC input file of the main chain is exactly the same as that of n-paraffins, except for the main chain carbon number (MCCN). Different types of branches are also allowed for the substitutes branch length (LBr) longer than C3. n-propyl and iso-propyl are allowed for LBr ) 3, whereas four different butyls (n-, iso-, sec-, and ter-) are used for LBr ) 4. The branch types used in the molecular construction are listed in Table 4. MOPAC input characteristics are generated separately for these

different types of branches, which will be used for the overall molecular construction. As shown in Figure 11, iso-paraffin isomers are permuted first for a given CN. For the iso-paraffins with a limited number of isomers (CN e 8), the complete set of isomers is determined and placed in the IsomerID (e.g., two, four, and eight isomers for C5-, C6-, and C7isoparaffins, respectively). IsomerID is an n × 12 matrix, where n represents the number of isomer permutations and the number 12 corresponds to the CN, MCCN, number of branches (#Br), lengths of 3 branches (3LBrs), types of 3 branches (3TBrs), and positions of 3 branches (3PBrs), accordingly. For the isoparaffins with CN g 9, the isomer permutations are observed to be automatically subject to the constraints from the molecular generation rules (max #Br ) 3, max LBr ) 4, no quaternary carbon, only one nonmethyl branch is allowed). The permutations proceed until all isomers are processed. The parts of the MOPAC input files corresponding to the main chain and branches are generated separately using the information contained in IsomerID. Those file fragments, which correspond to various branches, then are inserted into the main chain MOPAC input file one by one to complete the construction of the MOPAC input file for the isoparaffin isomer at hand. The constructed initial input file is used by the MOPAC software for geometry optimization and calculation of ∆H°f for stability screening. As mentioned previously, the mono-methyl-branched isoparaffins are abundant in crudes but are not thermodynamically favored (see Figure 7). To prevent these molecules from being screened out, they are excluded from the screening procedure and selected directly. The rest of the molecules are screened using the same ∆H°f threshold. As shown in Figure 11, this procedure repeats until the last isomer is constructed and then the program returns to the main program (Figure 10), increments CN, and repeats until all the isoparaffins have been constructed. 2.2.4.2. Molecular Generations of Cyclic Hydrocarbons. The principle for generating cyclic hydrocarbon molecules is similar to that for building iso-paraffins, except that the main structure here is the ring core. Some cyclic hydrocarbons could have more than one type of ring-core structure at a given CN (in this work, only the five- and six-membered rings are considered). For the cyclic hydrocarbons, for Z values from 0 to -24, in middle distillates, the ring structure arrangements are limited, despite the relative large number of rings (e.g., four ring-core structures corresponding to tetra-aromatics, Z ) -24). The likely-to-exist ring-core structures, corresponding to the hydrocarbon types, HCType ) 3-16 with Z ) 0 to -24, are depicted in Table 5. For the homologous series with complex ring structures (Z ) -4, -10, -14, -20, and -22), molecular simulations have been conducted to determine the five that are most thermodynamically favorable. Those five are selected as the ring cores to assemble the corresponding molecules. As shown in Figure 12, HCType and CN are the required inputs. For a given HCType and CN, the permutation procedure selects a ring core from the data presented in Table 5 and then all the possible isomers are permuted using the rules for the selected ring core.



Table 5. Ring-Core Structures of Cyclic Hydrocarbons Representing Molecules Identified by GC-FIMS

2390


Figure 12. Programming diagram for generation of cyclic hydrocarbons.

When the total length of substituents is less than C3, the permutation routine goes through the nonsymmetrically available positions on the ring core and directly fills the n × 14 IsomerID matrix, which containing n isomers. The 14 entries of each isomer are as follows: CN, #Br, four LBrs, four TBrs, four PBrs. If the total CN of all the substitutes is >3, the permutation routine selects #Br ) 1, 2, 3, or 4 and then permutes the branches over the available sites (as numbered in the ring-core data presented in Table 5). Following the molecular generation rules, the substitutes are distributed over the ring core as far apart as possible with sorted permutations of LBr, TBr, and PBr. The results are placed in matrix IsomerID. The procedure repeats until the last position, and then the last type, are reached. The MOPAC input files for a selected ring-core structure is generated first. Then, from the substitution attributes contained in IsomerID, cyclic hydrocarbons are generated by adding the input file fragments representing the selected branches to the MOPAC input file of the ring core. The same types of substituents as those used for isoparaffins were used to construct cyclic hydrocarbons, except that substituents longer than C5 (LBr g 5) were allowed. Based on the occurrences of substituted cyclic structures, only linear and isobranched-linear substitutes of LBr g 5 are utilized in the construction of cyclic molecules, as shown in Table 4. As done previously, the generated molecules are screened using the MOPAC-generated ∆H°f value. This

Ha et al.

procedure repeats until the last isomer is constructed. If there is more than one ring-core structure listed in Table 5 for the selected HCType, the program goes to the next ring core and repeats the permutation and assembling procedures described above until all the isomers for the last ring core are generated. The routine then returns back to the main program (see Figure 10) and loops through CN and HCType. All the spatially optimized and stability-screened molecules are stored in the molecule bank in the form of the MOPAC input files. These files are later used by the QSPR models to estimate their physical properties (BP, density, and RI). Thermodynamic calculations are conducted by MOPAC to calculate the entropy S, and then the free energy of formation ∆G°f. With these properties assigned, the selected molecules are ready to be manipulated to estimate the isomer concentration distribution and bulk property predictions, as described below. 2.3. Predictions of the Boiling Point, Specific Gravity, and Refractive Index. When the concentrations of the individual molecules are known, predictive models for BP, SG, and RI and suitable mixing rules can be used to calculate the bulk properties for the molecular representation. Unfortunately, most group contribution models that are available in the open literature do not offer the accuracy that we pursued. Therefore, we developed a set of QSPR models for BP, SG, and RI estimation15 for pure hydrocarbons. The QSPR approach relies on correlating the property in question with the molecular structure of the compound captured in the form of several descriptors. Those can be divided into several groups: constitution, topology, geometry, electrostatics, quantum chemistry, and thermodynamics. Using statistical methods and large data sets (including 186 saturated and 200 aromatic hydrocarbons), we developed a set of empirical multilinear models that were demonstrated to be significantly better than the group contribution models commonly used for property estimation. The quality of these models is reflected in Table 6. Except for the SG model for aromatics (R2 ) 0.9881), the correlation coefficients of the models were better than 0.99. The standard deviations of these models were less than (6.2 °C for BP, (0.008 for SG, and (0.005 for RI. The leave-oneout cross-validation and external test data sets were also used to validate the models, with excellent results. 2.4. Isomer Distribution within an Isomeric Lump. It is impractical and probably impossible to experimentally determine detailed molecular compositions of petroleum fractions that boil at >200 °C. Modern analytical techniques provide, at best, compositional information in terms of the concentrations of isomeric lumps. One possible way to “de-lump” this information is to distribute the concentration of the isomeric lump among individual isomers, assuming thermodynamic equilibrium. Based on the assumption of an ideal solution that is considered appropriate for petroleum mixtures, Smith and Missen31 presented a thermodynamically rigorous treatment of equilibrium among isomers in a closed system. However, compositional analyses of various crude oils revealed that the (31) Smith, W. R.; Missen, R. W. Chemical Reaction Equilibrium Analysis: Theory and Algorithm; Wiley: New York, 1982.



Table 6. Performance of Quantitative Structure-Property Relationship (QSPR) Models on Boiling Point (BP), Specific Gravity (SG), and Refractive Index (RI) Predictions for Hydrocarbons Correlation Coefficient model

Fisher number, F

R2

R2(CV)a

s

s (test)

r% (test)b

BP for saturates BP for aromatics

13332.1 6351.6

0.9979 0.996

0.9976 0.9954

(5.87 °C (6.2 °C

(6.1 °C (7.24 °C

0.68 0.81

SG for saturates SG for aromatics

1925.4 1820.9

0.991 0.9881

0.9894 0.9863

(0.007 (0.008

(0.01 (0.01

1.0 0.8

RI for saturates RI for aromatics

3054.3 2052.26

0.9921 0.9902

0.9902 0.9882

(0.004 (0.005

(0.006 (0.01

0.4 0.53

a

Cross-validation correlation coefficient. b Average relative error (expressed as a percentage).

reported concentration distributions of alkane isomers in many crude oils were quite far from thermodynamic equilibrium.19 As reported elsewhere,14 we have developed a procedure to estimate the concentration distributions of individual isomers within the isomeric lump. This can be achieved by the minimization of the free energy of the lump subject to the stiochiometric constraint and an additional constraint representing partial knowledge of the composition that is measured independently of the GC-FIMS matrix. In this work, we used the molal average boiling point (BPlump) of the lump as this additional constraint. As a result, the problem becomes that of a search for the equilibrium composition, subject to the extra average BP constraint. Mathematically, this corresponds to constrained minimization of the free energy (∆G°n) of the system

min[∆G°n] ) min subject to

{

[

N

∑ i)1

N

xi∆G°i + RT

xi ln xi ∑ i)1

is due to the averaging used. In principle, the entire measured distribution could be used in our de-lumping procedure. 2.5. The Minimum Molecular Representation. Using the molecular representation generated above, it is possible to calculate several bulk and distributed properties of the sample at hand. Sorting all the molecules by calculated BPs and summing the corresponding mass will generate an equivalent SimDis curve. After the physical properties of interest here are either predicted using the QSPR correlations or calculated from the available molecular information, the bulk properties (SG, RI, and MW) of the mixture can be predicted using the following simple mixing rules:

1

)

SG

]

1

xwi

∑i SG

(3)

i

)

MW

xwi

∑i MW

(4) i

and

∑xi ) 1 ∑xiBPi ) BPlump

RI )

∑i xwiRIi

(5)

(2)

It was proven analytically that this problem has a unique solution.14 This approach was applied to predict the isomeric distribution among hexane/heptane isomers using the standard free energy of formation and the normal boiling point of the isomers involved.14 The predicted results were compared with the reported data of 18 geologically different crude oils.21 The calculated concentration distributions matched the reported data satisfactorily, except for three younger crudes formed during the tertiary Cenozoic Era (South Houston, Wafra, and Wilminton). In particular, this new procedure can predict the key isomers with average relative prediction errors of 15% and 12% for hexane and heptane isomers, respectively. These key isomers, three out of five for hexanes and four out of nine for heptanes, account for more than 92%/90% of the total hexane/ heptane isomers reported.21 It should be mentioned here that as CN increases, BPlump may become insufficient to solve this problem. However, BPlump is calculated from a BP distribution of the lump and involves some loss of information that

where xwi is the mass fraction of molecule i. It is worth noting that this mixing can be performed over the entire molecular representation to obtain bulk SG, MW, and RI values or over narrow boiling-range intervals to obtain BP distributions of these properties rather than bulk values. This feature is important for the simulation of phase equilibrium, because commercial phase equilibrium software usually requires physical properties of the individual pseudo-components (defined as narrow boiling-point fractions) rather than bulk properties of the flashed mixture. In this work, we chose to use the bulk properties calculated for the representation to validate our procedure by comparing those properties with experimental data. If the stability screening threshold (2.0 kcal/mol) is adequate, it should be possible to find good agreement between the predicted and measured bulk properties. We also made an attempt to determine the minimum number of molecules that reproduces the bulk properties with the desired accuracy. For computational efficiency, for a molecular representation that has already been developed, the number of molecules was adjusted by setting the threshold value for the smallest concentration allowed rather than using the heat of formation threshold and re-derive the molecular representation.

2392


Ha et al.

Table 7. Simulated Bulk Properties, Compared with Experimental Results, for Five Samples Molecular Weight, Mw (amu)

Density, F @ 15.6°C (g/mL)

Refracive Index, RI @ 25°C

Deviation in SimDis (°C)

sample

initial set of molecules

measured

calculated

measured

calculated

measured

calculated

IBP

5%-95% recovery range

FBP

S1 S2 S3 S4 S5

801 1588 1531 1716 1416

150 188 176 176 175

148.4 187.8 171.2 175.3 179.9

0.7985 0.8291 0.8357 0.8476 0.8537

0.7986 0.8361 0.8434 0.8477 0.8533

1.4436 1.4593 1.4584 1.4652 1.4755

1.4449 1.4653 1.4619 1.4672 1.4767

9.2 4.8 18.2 11.8 2.7

2.0 3.0 2.4 3.8 2.6

7.8 14.2 2.2 7.2 2.6

The accuracy of bulk property prediction can be measured as the tolerance:

Tol ) |SGi - SG0| + |RIi - RI0| + |(Dev•SimDis)i (Dev•SimDis)0| (6) where

Dev•SimDis )

1

N

∑ (xw,pred - xw,exp)j

N j)1

(7)

The previously noted subscripts “0” and “i” refer to the initial and minimized molecular representations, respectively. N is the number of boiling-point intervals in the experimental SimDis curve chosen for comparison, and xw is the mass fraction of the interval j. The three terms in eq 6 have similar orders of magnitude; hence, we did not use any weighting factors. As a result, the smallest molecular representation can be found for a given tolerance. 3. Molecular Representation and Bulk Property Predictions As mentioned previously, we tested the validity of the de-lumping procedure by comparing the predicted and measured bulk physical properties of the representations and the corresponding samples, respectively. The selected samples covered a variety of compositions (i.e., highly paraffinic (S1), highly naphthenic (S3), and highly aromatic (S5)). With the initial screening threshold, ∆H°f ) 2.0 kcal/mol, the initial representations ranged from 801 molecules to 1716 molecules for samples S1-S5, as listed in Table 7. The density and RI for samples S1-S5 were calculated from their molecular representations, using eqs 3 and 5. In addition, the average MW of samples S1-S5 were directly estimated from their GC-FIMS characterization matrixes (see example in Table 2) using eq 4, where xwi is the measured mass concentration of the isomeric lump i. These calculated properties agree well with the corresponding measurements for the five middle-distillate samples, as compared in Table 7 again. The maximum absolute prediction errors were (0.008 g/mL for density and (0.006 for RI, equivalent to the standard deviation of the QSPR models used for the estimation of the properties of pure compounds. The average prediction error for both bulk density and RI was only (0.003 for all five samples. The agreement between the predicted and measured average MW for samples S1S5 was determined to be satisfactory, with the average absolute error of (2.44 amu, which is comparable with the experimental error of the freezing-point depression method (1%-2%). The maximum calculated deviation was only 4.9 amu.

Furthermore, as another means of verifying the analytical data and the characterization method used in this work, two comparisons of the calculated and measured SimDis curves are shown in Figures 2-6. The first one is the comparison between the GC-FIMSgenerated (“FIMS-gen SimDis”) and the measured SimDis data, as introduced in the previously presented data reconciliation section. The second comparison is between the SimDis predicted from representative molecules (“Simulated SimDis”) and the measured ones. Generally, good agreement was observed in all these comparisons. The deviations at the initial boiling point (IBP) (0.5 wt % cut point) and final boiling point (FBP) (99.5 wt % cut point) from simulations for each sample are also shown in Table 7, together with the average deviation over the 5%-95% recovery range. The average deviations for IBP and FBP over all five samples were 9.3 and 6.8 °C, respectively, whereas the average deviation in the 5%-95% recovery range was only 2.7 °C. For comparison, the repeatabilities of the ASTM D2887 SimDis measurements were 6.0 and 5.0 °C for IBP and FBP, respectively. The repeatability in the 5%95% recovery range is 2.5 °C. The corresponding reproducibilities were 23.0 and 13.5 °C for IBP and FBP, and 6.5 °C over the 5%-95% range. The good agreement between predicted and measured bulk properties indicates that a value of ∆H°f ) 2.0 kcal/mol was a safe choice for the molecule screening threshold. 4. Sensitivity of Estimated Properties to the Size of Molecular Representations Next, we explored the consequences of reducing the number of molecules in the representation for the estimated bulk properties. As indicated in Tables 7 and 8, the tolerance of 10-4 reduced the number of molecules by a factor between 12% (sample 1) and 37% (sample 3). Table 8 lists the number of molecules of each hydrocarbon type, given a tolerance of 10-4. The losses of the overall mass due to the removal of the minor molecules were

Derivation of Molecular Representations of Middle Distillates

Recommend Documents