Preliminary Results on Molecular Modeling of Asphaltenes Using

Energy & Fuels 1996, 10, 97-107

97

Preliminary Results on Molecular Modeling of Asphaltenes Using Structure Elucidation Programs in Conjunction with Molecular Simulation Programs I. Kowalewski,*,† M. Vandenbroucke, and A. Y. Huc Division Ge´ ologie-Ge´ ochimie, Institut Franc¸ ais du Pe´ trole, 1-4 Avenue de Bois Pre´ au, Rueil-Malmaison, 92506 Cedex, France

M. J. Taylor TACIT, Postfach 1250, D-61186 Rosbach v.d.H, Germany

J. L. Faulon Sandia National Laboratories, Albuquerque, New Mexico 87185 Received June 9, 1995. Revised Manuscript Received August 8, 1995X

Molecular modeling using structure elucidation programs in conjunction with molecular simulation programs has been performed on asphaltene molecules, the heaviest fraction of crude oil, in order to obtain a chemical model allowing us to tentatively study their physicochemical properties. We have analyzed Boscan asphaltenes (Venezuela) derived from a marine source rock. The different steps of this molecular modeling are described. First, a 3-D chemical representation of Boscan asphaltene is defined from an analytical data set. Second, the results of molecular dynamic simulations indicate that only a few stable conformations are possible due to the high reticulation of the model of the asphaltene unit obtained.

Introduction During this last decade, asphaltene fractions of petroleum have been the focus of much attention due to the potential production problems during exploitation of fields caused by the viscous and flocculating nature of these products. In the extreme limit, one implication of asphaltene precipitation can be the nonproduction of oil due to a global reduction of reservoir permeability. But, since asphaltenes are poorly characterized on a molecular basis,1 the causes of precipitation are not well understood. In order to better understand the physicochemical behavior of asphaltenes, it is of interest to know more about their chemical structure in terms of macromolecular groups (primary bonding structure) and intraand intermolecular bonds (secondary bonding structure). In addition, another essential parameter is the determination of the average molecular weight in their natural environment, the other constituents of oils, and in an isolated state. Asphaltenes consist of a macrostructure made up of an aggregation of asphaltene “units” comprising Π and hydrogen bonds. They can be viewed as geomacromolecules including heterogeneous units randomly connected by covalent bonds (primary bonding structure); because of their mobility their secondary bonding structure is probably in an energy-minimized configuration. These units are constituted of more or less condensed aromatic rings * To whom correspondence should be addressed. † E-mail: [email protected]. X Abstract published in Advance ACS Abstracts, October 15, 1995. (1) Speight, J. G. Prepr.sAm. Chem. Soc. Div. Pet. Chem. 1989, 321-328.

0887-0624/96/2510-0097$12.00/0

carrying alkyl or naphthenic substituents with heteroatoms (especially N, S, O) interspersed within the system.2-5 The resin fraction of petroleum, containing smaller and more polar molecules, could play an important role in the dispersion of the asphaltenes in the crude oil. Some authors consider that asphaltenes form colloids which are solubilized by resins.6-9 Different asphaltene models have been already proposed. A classical model developed by Yen et al.10,11 describes aggregation of asphaltene “units”. These authors have proposed a three-dimensional structure for asphaltene characterized by the stacking of more or less polycondensed aromatic sheets linked by Π-Π bonds, hydrogen bonds, or sulfur and oxygen bridges.12 Several of these units are associated and form micellar aggregates. This model has been defined using data collected from solid phase compounds. Other models, (2) Yen, T. F. Prepr.sAm. Chem. Soc. Div. Pet. Chem. 1972, 17(4), F.102-F.114. (3) Yen, T. F. Prepr.sAm. Chem. Soc. Div. Pet. Chem. 1979, 24(4), 901-909. (4) Tissot, B. Rev. Inst. Fr. Pet. 1981, 36(4), 429-446. (5) Tissot, B.; Welte D. H. Petroleum Formation and Occurrence; 2nd ed.; Springer: Verlag, Berlin, 1984; p 538. (6) Pfeiffer, J. P.; Saal, R. N. J. J. Phys. Chem. 1940, 44, 139-149. (7) Speight, J. G.; Moschopedis, S. E. Prepr.sAm. Chem. Soc. Div. Pet. Chem. 1979, 24, 910-923. (8) Moschopedis, S. E.; Fryer, J. F.; Speight, J. G. Fuel 1976, 55, 187-192. (9) Pelet, R.; Behar, F.; Monin, J. C. Org. Geochem. 1986, 10, 839846. (10) Yen, T. F.; Erdmann, J. G.; Pollack, S. S. Anal. Chem. 1961, 33(11), 1587-1594. (11) Sadeghi, M. A.; Chilingarian, G. V.; Yen, T. F. Energy Sources 1986, 8, 99-123. (12) Ignasiak, T. M.; Bimer, J.; Samman, N.; Montgomery, O. P.; Strausz, D. S. Prepr.sAm. Chem. Soc. Div. Pet. Chem. 1979, 24(4), 1001.

© 1996 American Chemical Society

98

Energy & Fuels, Vol. 10, No. 1, 1996

published in the literature13-15 and developed from asphaltenes in solution, lead to a bidimensional planar organization (disks). More recently, Ravey et al.16-18 have investigated the macrostructure of asphaltene particles in solution using small angle X-ray scattering (SAXS) and small angle neutron scattering (SANS) and seem to confirm a bidimensional model for asphaltenes when dispersed in an efficient solvent. Their results suggest the existence of a colloidal and micellar nature for the asphaltene dispersion involving reversible association of asphaltene units. With flocculating agents, the association between these large primary particles would lead to the formation of a three-dimensional network. To investigate the physicochemical behavior of asphaltenes, molecular modeling offers new insights into the chemical structure of asphaltenes by taking into account the contribution of molecular groups and the bonds which are involved. This method seems to be of interest for studying their chemical structure and behavior, and for following the physicochemical relationships that exist between them and the other constituents of the oil. The aim of this paper is to present the pertinence of molecular modeling for the understanding of asphaltene molecules by combining structure elucidation programs (XMOL and SIGNATURE software) and molecular simulation programs (INSIGHT II/DISCOVER software). XMOL (developed by IFP) and SIGNATURE (developed by J. L. Faulon) allow a 3-D graphic representation of the asphaltene “molecule” generated from a consistent quantitative analytical data set. Then, this molecular model was used as an input for INSIGHT II/ DISCOVER software provided by Biosym Technologies. Starting from the final structure obtained by SIGNATURE, INSIGHT II (a graphic molecular modeling program) allows the building and the modification of the molecular system. DISCOVER (molecular simulation program), used in conjunction with INSIGHT II, calculates the minimum energy of the most stable conformation(s). Experimental Section Sample. To perform this study, we have isolated and analyzed the asphaltene fraction of Boscan crude oil (type II, Venezuela). This crude oil contains a high amount of asphaltenes (13%) and does not present any asphaltene precipitation during production operation. Analytical Procedure. The asphaltene fraction was prepared from the topped crude oil (fraction boiling above 210 °C) using n-heptane precipitation. The asphaltenes were then characterized by two types of analyses: Atomic Analyses and Molecular Weight. These analyses (elemental analysis, density, 13C NMR, and gel permeation chromatography) define the chemical properties of the asphaltene molecule by quantifying atoms and functional groups and molecular size. The density was determined using helium pycnometry at room temperature with primary vacuum (10-3 (13) Reerink, H. Ind. Eng. Chem. Prod. Res. Dev. 1973, 12, 82. (14) Reerink, H. Prepr.sAm. Chem. Soc. Div. Pet. Chem. 1971, 16(1), D18. (15) Dickie, J. P.; Yen, T. F. J. Colloid Interface Sci. 1969, 29, 475. (16) Ravey, J. C.; Espinat, D.; Ducouret, G. Fuel 1988, 67, 15601567. (17) Ravey, J. C.; Espinat, D. Prog. Colloid Polym. Sci. 1990, 81, 127-130. (18) Espinat, D.; Ravey, J. C. SPE-Oilfield Chem., Proc. of Int. Symp. Soc. Pet. Eng., New Orleans 1993, 365-373.

Kowalewski et al. Table 1. Elemental Analysis, Liquid 13C NMR, and Density for BOSCAN Asphaltenes elemental analysis

wt % C H O N S ashes

79.66 7.76 1.79 1.73 6.83 2.23

SIGNATURE parameters (atomic ratios) (h 116.9)d (o 1.685)d (np 1.86)d (s 3.215)d

13C

NMR (C %)

faa Caro.-H subst. Caro. cond. Caro. degree of condensation degree of substitution dipolar dephasingb CH CH2 + CQc CH3

39.6 11.3 13.9 14.4 4.5

density 1.24

0.551 28.5 5.5 42.6 12.8

a f × 100: carbon aromaticity %. b Caro.-H/f . c CQ: quaternary a a carbon. d INSIGHT II/DISCOVER nomenclature (see Glossary).

Torr). The liquid 13C NMR data were obtained on a MSL 400 Bruker spectrometer using a solution of 100 mg of asphaltenes in 2 mL of CDCl3 with 15 mg of iron acetylacetonate as relaxation agent. To calculate aromaticity, “one-pulse” sequences were used with proton inversed gated decoupling. Pulse time was 5 µs corresponding to a 60° angle with a time delay (D0) of 7 s. Other averaged structural parameters were determined using spin echo and attached proton test sequences with echo time (τ) of 8 and 3.125 ms, respectively, for saturated and aromatic carbons. 5000 acquisitions were accumulated. The results of atomic analysis are summarized in Table 1.

Care must be taken when determining asphaltene molecular weights because the analysis is complicated by the tendency of asphaltene molecules to associate with each other and with other surrounding constituents. It is known that according to the techniques which are involved (viscosity, vapor pressure osmometry, ultracentrifugation, small angle neutron scattering, ...), the molecular weight varies within a large range.8,19-23 We have calculated the average molecular weight in THF solution using gel permeation chromatography (GPC). Our analyses have been carried out using a HPLC pump (P2000 Spectra Physics) operating at a flow rate of 1 mL/min on two consecutive PLgel columns (1000 and 100 Å pore size) following a PLgel guard column. THF was chosen first because it is a good solvent for asphaltenes and second because it has a relatively high dipolar moment and provides relatively low molecular masses due to low association phenomena; the GPC chromatograms obtained for asphaltenes show a better defined profile distribution when compared to those obtained with other solvents (CH2Cl2, benzene, ...).23 In addition, and in order to minimize as far as possible intermolecular interactions, the asphaltene solution was very diluted. Using these experimental conditions, the measured molecular masses are assimilated to those of individual covalent “units” of asphaltene. The detection was done in absorbance mode monitoring at 260 nm (UV2000 Spectra-Physics). This wavelength is taken as the most appropriate to characterize the polynuclear aromatic units. The mass calibration of the system was fixed using polystyrene standards. In principle, standards of the same molecular type as the substance analyzed should be (19) Monin, J. C.; Pelet, R. Adv. Org. Geochem. 1983, 839-846. (20) Monin, J. C.; Vignat, A. Rev. Inst. Fr. Pet. 1984, 39(6), 821827. (21) Speight, J. G.; Wernick, D. L.; Gould, K. A. Overfield R. E. et al. Rev. Inst. Fr. Pet. 1985, 40(1), 51-61. (22) Sheu, R. Y.; Storm, D. A.; De Tar, M. M. J. Non-Cryst. Solids 1991, 131-133, 341-347. (23) Ducouret, G. Fractionnement des asphalte`nes par chromatographie d’exclusion ste´rique- Caracte´risation des fractions par diffusion des neutrons aux petits angles. Doctoral Thesis, University of Paris, VI, December 1987.

Molecular Modeling of Asphaltenes

Energy & Fuels, Vol. 10, No. 1, 1996 99

Table 2. Determination of Average Molecular Weight by GPC and Small Angle X-ray Scattering

Scheme 1. Molecular Analyses on Saturated C7+ Fraction

polydispersed cylinder model units radius (Å) asphaltenes MMpa MMwb number BOSCAN a

10000

30800

3

21.3

σc

thickness (Å)

0.725

10.8

b

Calculated by GPC. Calculated by small angle X-ray scattering. c Standard deviation of the distribution.

required to convert GPC chromatograms into molecular weight distributions and to calculate average molecular weights. It is obvious that polystyrenes do not represent the best standards for asphaltenes due to the large molecular polydispersity of the latter. However, the equivalent polystyrene weight data are given here to allow for some standardized reference and calculate relative molecular weight distributions. The areas of the recorded peaks were integrated using a GPC software developed at IFP. However, important tailing toward the low masses revealed the occurrence of significant adsorption phenomena with the PLgel. The average molecular weight MMp used in our model, reported in Table 2, was calculated by integrating the area corresponding to the first “peak” (higher molecular mass) in the chromatogram. The average molecular masses obtained in these conditions are likely to represent the smallest unit of asphaltenes which is assembled by covalent bonds. A second technique, small angle X-ray scattering, was used for measuring the average size distribution of asphaltenes. This technique is very sensitive to aromatic environments.24 Asphaltenes were diluted in toluene (2 wt %) and X-ray measurements were done at room temperature. The MW obtained by this technique (Table 2) are much higher than with the previous technique. Therefore, several of these covalent units separated by THF must be associated to account for the result obtained by small angle X-ray scattering. Molecular Analyses. These analyses allow the identification and the quantification of the molecular groups which constitute the “asphaltene molecule”. They were performed on the hydrocarbon fragments of the degradation products of the macromolecule obtained using anhydrous closed system pyrolysis. This technique has been described elsewhere.25,26 About 50 mg of asphaltenes were pyrolyzed in a closed gold reactor without water added. A series of isothermal heating at 10 MPa under argon during different times were performed in order to determine the appropriate experimental conditions to generate the maximum amount of hydrocarbons without secondary cracking of heaviest fraction. The best conditions were obtained for heating at 350 °C during 9 h. The molecular quantitation of the gas phase (H2, CO, CO2, C1-C5) was carried out by GC-TCD after a prior recovery and fractionation by a Toepler pump.25,27,28 The liquid pyrolysate was dissolved in n-pentane, recovered, quantified by GC/FID, and finally fractionated by minicolumn liquid chromatography into saturates, aromatics, and NSO compounds. The molecular quantitation of saturates C7+ and aromatics was performed according to the procedures developed at IFP (Schemes 1 and 2). The total C7+ saturated fraction was analyzed by gas chromatography with helium (He) as carrier gas using a FID detector set at 300 °C (Varian 3500) equipped with an on-column injector (temperature programmed from 50 to 300 °C at 180 °C/min). The capillary column was a DB1 30m (i.d. 0.32 mm; f.t. 0.1 µm). The GC was temperature programmed (24) Bardon, C.; Barre´, L.; Espinat, D.; et al. Fuel Sci. Technol. Int., in press. (25) Behar, F.; Leblond, C.; Saint-Paul, C. Rev. Inst. Fr. Pet. 1989, 44(3), 387-411. (26) Behar, F.; Kressmann, S.; Rudkiewicz, J. L.; Vandenbroucke, M. Org. Geochem. 1992, 19, 173-189. (27) Behar, F.; Ungerer, P.; Kressmann, S.; Rudkiewicz, J. L. Rev. Inst. Fr. Pet. 1991, 46(2), 151-181. (28) Prinzhofer, A.; Huc, A. Y. AAPG. Special Vol., in press.

from 50 to 100 °C at 10 °C/min and from 100 to 300 °C at 3 °C/min. Saturates were then fractionated into n-alkanes and iso-/cycloalkanes by insertion into 5 Å molecular sieves. The iso-/cycloalkanes were quantified (in wt %) by mass spectrometry on a Fisons Autospec with solid probe as introduction mode, using the Hood and O’Neal method.29 The operating conditions were as follows: ion source 270 °C; electron energy 70 eV; ionization mode: EI; voltage acceleration 8 kV; resolution 1000; acquisition time 3 s/decade. The solid probe was heated from 40 to 500 °C at 30 °C/min. The sulfur content of total aromatics was determined by elemental analysis. The aromatic fraction was separated according to the number of condensed aromatic rings by preparative normal phase HPLC (Varian Star 9010 pump) on Petrospher B (Chrompack, 250 × 4.6 mm). The solvent used for isocratic elution was n-hexane, and the flow rate was maintained at 1 mL/min. The elution time of each aromatic family was determined using a mixture of reference compounds and the chromatographic column was back-flushed after elution of pyrene. The detection was done using monitoring absorption at 254 nm (Varian 9050 UV detector). Four fractions were collected corresponding respectively to monoaromatics + thiophenes (F1), diaromatics + benzothiophenes (F2), fluorenes + dibenzothiophenes (F3), and phenanthrenes + polyaromatics (F4). Each of these subfractions was then quantified by solid probe-MS (Ultima Fisons) using Fisher group type analysis.30,31 A 0.5-µL volume of the sample diluted at 1/1000 with n-hexane was added to 2 µL of PFK (perfluorokerosene) and injected. The operating conditions of the mass spectrometer were as follows: mass range 20-600 amu, electron energy 70 eV; ionization mode: EI; voltage acceleration 8 kV; resolution 10 000; acquisition time 3 s/decade. The solid probe was heated from 40 to 650 °C at 50 °C/min. All the results of molecular analysis are presented in Table 3, b and c.

Software After the complete analytical data base was gathered, two structure elucidation programs XMOL and SIGNATURE were used. XMOL software, developed at IFP,32,33 was originally written for modeling geological macromolecules such as kerogens and asphaltenes. XMOL allows the prediction and the construction of average chemical structures of geomacromolecules in a 3-D space according to their quantitative physicochemical analyses. The software is divided into three modules. The (29) Hood, A.; O’Neal, M. J. Annual Book of ASTM Standards; ASTM D2786; ASTM: Philadelphia, Part 24, pp 367-374. (30) Fisher, I. P.; Fisher, P. Talanta 1974, 21, 867. (31) Bouquet, M.; Brument, J. Fuel Sci. Technol. 1990, 8(9), 961986. (32) Faulon, J. L.; Vandenbroucke, M.; Drappier, J. M.; Behar, F.; Romero, M. Org. Geochem. 1990, 16(4-6), 981-993. (33) Faulon, J. L. The`se de l’Ecole des Mines de Paris, January 1991.

100


Kowalewski et al.

Scheme 2. Molecular Analyses on Aromatic C7+ Fraction

first one, a prediction program, is based on the fingerprint of the macromolecule, i.e. equations allowing on the one hand to quantitatively select, from a library, molecular groups (aromatics, naphthenoaromatics, naphthenics, aliphatic chains, ...) according to molecular analyses (such as GC or MS of free hydrocarbons in related oils or Py-GC analyses of asphaltenes) and on the other hand to find the bonds that link the molecular groups in order to verify atomic analyses (elemental analysis, 13C NMR, IR, or UV spectroscopy, functional groups, ...). When analytical data were not available, assumptions had to be made. The program calculates and lists in a construction file the number of each molecular group and intergroup bonds selected according to the analyses. These numbers are computed for a given molecular weight of the macromolecule. The second module is an editor that enables to construct and store atom configurations (mass, diameter), bonds (length, direction), cyclic molecular groups (including steric energy minimization), and the intergroup bonds. The molecular structures built by the editor are stored in a library. Due to the fact that the kerogen libraries were not specifically adapted to asphaltene modelling, we had to create for this study a specific library which was as simple as possible and consistent with the type of molecular groups determined by molecular analyses. The third module is a construction program which allows to construct a 3-D isomer of the macromolecule calculated by the prediction program. This construction program works by taking molecular groups and intergroup bonds randomly from the construction file and linking them inside a chosen construction volume based on the measured density and molecular weight of the macromolecule. For each addition of a new molecular group, the nonintersection with previous groups is verified (minimization of the steric energy is not included in the software). It is possible to construct several isomers by changing the random sorting procedure in the program. Recently, Faulon developed a new computer-assisted structure elucidation program called SIGNATURE which is based on graph theory.34 SIGNATURE, unlike XMOL, is compatible with other molecular modeling software. This program allows large molecules to be assembled and has been already successfully applied to coal,35 lignin and lignocellulose.36 In order to perform asphalt(34) Faulon, J. L. J. Chem. Inf. Computer Sci. 1994, 34, 1204-1218 .

enes modeling we need first to use XMOL designed to calculate the number of each molecular group to be entered in the SIGNATURE file. The first task performed by SIGNATURE program is to compute the number of each fragments and interfragment bonds that best match with the quantitative data. Moreover, the program is able to determine the number of nonidentical models than can be constructed from a set of analytical data. One or more macromolecules randomly chosen among the possible isomers is constructed by applying the prediction program (which links the fragments with interfragment bonds). But SIGNATURE, having no builder, must be used in conjunction with other molecular modeling programs to build the model in a 3-D space. This was done using the INSIGHT II software provided by Biosym Technologies. The graphic molecular modeling program INSIGHT II is compatible with the output file obtained from SIGNATURE and allows the user to build and modify molecular systems interactively. Because the 3-D BOSCAN structure generated by the SIGNATURE program was not initially at a minimum state of energy, the model was submitted to molecular mechanics (MM) and molecular dynamics (MD) simulations. We performed these simulations using in conjunction with INSIGHT II and DISCOVER provided by Biosym Technologies. These simulations are based upon the use of a general valence force field which requires a set of equations and associated constants designed to reproduce molecular geometry and selected properties of test structures. For our study, we used the CVFF (covalent valence force field). MM was used to find the energy minimum closest to the initial conformation. During the energy minimization, the spatial positions of atoms are perturbed in a stepwise manner, until a minimum-energy conformation is identified. To identify this minimum, the minimization calculations were carried out using a conjugate gradient algorithm. Complex structures such as asphaltenes can be expected to display a number of low-energy conformations separated by more or less high-energy barriers. To study the possible stable conformations of Boscan asphaltene, MD were employed by application of quenched annealed dynamics (QAD) technique. In MD calculations, atoms are allowed to move with velocities depending on the temperature of the surrounding (35) Faulon, J. L.; Hatcher, PG.; Carlson, G. A.; Wenzel, K. A. Fuel Process. Technol. 1993, 34, 277-293. (36) Faulon, J. L.; Carlson, G. A.; Hatcher, PG. Org. Geochem. 1994, 21(12), 1169-1179.



Table 3. Analytical Files Entered in XMOL for BOSCAN Molecule: (a) Atomic Analyses, (b) Molecular Analyses, Saturated Hydrocarbons, and (c) Molecular Analyses, Aromatic Hydrocarbons (a) Atomic Analyses functional analysis 79.66 oxygen functional analysis (mequiv/g) 7.76 total hydroxyl 1.79 total carbonyl group 1.73 acid 6.83 acid H quinone 60.35 (c )a methoxyl 0 nitrogen distribution (%) 0.05 N-aromatic 0 sulfur distribution (%) 39.35 (cp)a S-aromatic 0.20 cyclosulfur and disulfur 0.05 thiol 0 28.5 3

elemental analysis (wt %) carbon hydrogen oxygen nitrogen sulfur 13C NMR (%) 30-35 ppm 56-58 ppm 72-76 ppm 92-96 ppm 128-130 ppm 154 ppm 178 ppm 210 ppm dipolar dephasing average no. condensed arom rings % of functional group linked to

0.20 0.24 0.16 0.22 0.0 0.0 100 80 20 0

2 aro

1 aro and 1 ali

2 ali

0 0 0 0 50

100 100 100 100 50

0 0 0 0 0

ester ketone ether amide sulfur

(b) Molecular Analyses: Saturated Hydrocarbons saturated hydrocarbons (%) alkanes

naphthenes

range of C

wt %

normal

ramified

isoprenoid

1 cycle

2 cycles

3 cycles

4 cycles

5 cycles

6 cycles

C1-C5 C6-C13 C14-C40

0.66 27.84 71.5

99 25 15

1 15 19.62

30 1.78

30 19.81

10.95

12.62

9.65

11.36

1

wt %

wt %

wt %

wt %

C

normal

ramified

C

normal

ramified

C

normal

ramified

C

normal

ramified

1 2 3 4 5 6 7 8 9 10

46.9 31.9 13.1 4.0 4.0 2.15 4.30 8.60 17.26 23.0

0 0 0 50 50 2.15 4.30 8.60 17.26 23.0

11 12 13 14 15 16 17 18 19 20

10.55 11.55 9.75 8.79 6.61 6.52 5.88 5.47 5.04 4.48

4.94 5.41 4.57 4.12 3.10 3.05 2.75 2.56 2.36 2.10

21 22 23 24 25 26 27 28 29 30

10.55 11.55 9.75 8.79 6.61 6.52 5.88 5.47 5.04 4.48

4.94 5.41 4.57 4.12 3.10 3.05 2.75 2.56 2.36 2.10

31 32 33 34 35 36 37 38 39 40

3.13 2.89 2.72 2.33 1.9 1.2 0.8 0.6 0.4 0.2

1.47 1.35 1.27 1.09 0.89 0.56 0.37 0.28 0.19 0.09

wt % pristane

wt %

74.5

phytane

25.5

(c) Molecular Analysis: Aromatic Hydrocarbons fraction

wt %

fraction

wt %

C6-C13

2.92

C14+

97.08

wt % benzene methylbenzene ethylbenzene xylene

wt %

3.72 16.28 4.18 27.44

trimethylbenzene naphthalene methylnaphthalene

aromatic structures CnH2n-pb

6 8 10 12 14 16 18 a

aromatic structures CnH2n-pb

wt %

p

25.11 6.97 16.30

p

3.3 3.4 22.3

20 22 24 26 28 30

3.1 5.1 10.8

wt % 7.9 11.7 7.2 6.4 6.9 3.6

7.4

INSIGHT II/DISCOVER nomenclature (see Glossary). b CnH2n-p: CnH2n-6 to CnH2n-10: monoaro; CnH2n-12 to CnH2n-16: diaro; CnH2n-18 to CnH2n-30: tri- and polyaro.

102


system. The motions of atoms can lead the structure to overcome potential energy barriers and not to remain trapped in a local energy minimum. QAD allows the system to overcome small energy barriers during the temperature relaxation and thus to settle the molecule into a lower energy minimum. The interactions between atoms are described by CVFF force field and the motions of atoms are determined stepwise by integrating Newton’s equation of motion.37 All the MD calculations were performed with an initial temperature set to 1100 K and for a period of time ranging from 0 to 150 ps. The temperature was then lowered from 1100 to 1 K in steps of 100 K (with 1 ps between drops). A structure was saved at each change of temperature (total of 10 structures). Furthermore, during each QAD simulation, the energy of the studied structure was minimized using the MM technique as described above. To summarize, we have first examined some possible structures of the BOSCAN asphaltene using computeraided structure elucidation and then studied the spatial conformations of the structures with a force field (CVFF) suitable for all atom types present in this asphaltene molecule. All the computations for this study have been performed on an Apollo workstation (XMOL) and on a Silicon Graphics Indigo-XZ4000 workstation (SIGNATURE version 1.3, INSIGHT II/DISCOVER version 2.3.0). The abbreviations used are described in the Glossary. Results and Discussion The XMOL software has been applied to BOSCAN asphaltenes as follows. All the analytical data obtained have been entered in analytical files (Table 3) and the consistency of the data verified. The following assumptions have been introduced to supplement missing analyses: 1. In the case of asphaltenes, the quantitative determination of functional groups is very difficult because they are present in very low concentration. Thus, the functional groups have been estimated by relying on the functional group distribution of type II kerogen standardized to the oxygen content of BOSCAN asphaltenes, with the assumption that quinone and methoxy groups are absent in these asphaltenes. The chemical degradation techniques used for functional group quantitative determination in coals and kerogens are described elsewhere.38 2. According to the literature on oil heavy ends,39-41 we have admitted that nitrogen was only in aromatic structures for the studied asphaltene. 3. For the sulfur compounds, according to the X-ray adsorption near edge structure spectroscopy (XANES) analyses performed at IFP, the distribution seems to be divided between thiophenic moieties and sulfide moieties such as disulfur bridges in a proportion of roughly 80% and 20% respectively. These figures are (37) Verlet, L. Phys. Rev. 1967, 159, 98-103. (38) Blom, L.; Edelhausen, L.; Van Krevelen, D. W. Fuel 1957, 36, 135-153. (39) Wilhelms, A.; Patience, R. L.; Larter, S. R.; Jorgensen, S. Geochim. Cosmochim. Acta 1992, 56, 3745-3750. (40) Schmitter, J. M.; Ignatiadis, I.; Dorbon M.; Arpino, P., et al. Fuel 1993, 63, 557. (41) Dorbon, M.; Ignatiadis, I.; Schmitter, J. M.; Arpino, P., et al. Fuel 1993, 63, 565.

Kowalewski et al. Table 4. Specifications for Construction construction parameters density no. of groups ratio of height in volume of construction random basis

1.24 50 2 1

in agreement with the literature42 and were used for BOSCAN asphaltenes. 4. Concerning the linkage of functional groups to aliphatic and/or aromatic carbons, for disulfur bridge distribution, we have equally assigned the bridges: 50% between two aromatic structures and 50% between one aromatic and one aliphatic structures. 5. For oxygen distribution, we have a priori assumed that carboxylic acids are 100% involved in aromatic structures. 6. Moreover, the molecular groups contained in the library have been adapted to account for the chemical composition of hydrocarbons released by pyrolysis of BOSCAN asphaltenes. We have based the library on naphthenoaromatic groups and, therefore, the average carbon number of molecular groups is C16 (naphthenoaromatic with four rings). Using the average molecular weight (MW) determined by gel permeation chromatography, this results in choosing approximately 50 molecular groups in the library for the construction file (Table 4). Once the assumptions were defined, a parallelepipedic box was selected for the total volume of construction with a height over square side in a ratio of 1/2 according to the dimensions obtained by small angle X-ray scattering. After running XMOL, the best construction file was selected (Table 5) and the asphaltene 3-D representation achieved (Figure 1). The bulk chemical formula and the molecular weight given by XMOL for BOSCAN asphaltene are summarized in Table 6. We can notice that the constraint of 50 molecular groups using the molecular groups as defined in the library was pertinent: the calculated molecular weight is similar to the experimental molecular weight obtained by GPC and obviously, the same is true for the elemental atomic ratios. The consistency of aromaticity (Csp2/Csp3 + Csp2) is excellent. The graphic display shows a major contribution of cycloalkanes and alkyl chains compared to aromatic rings as observed with analytical results. However, this model describes only an individual covalently bonded “unit” of the asphaltene molecule at an isolated state, and not a macrostructure made of an aggregation of asphaltene units linked by Π or H bonds. As can be observed on the 3-D representation, the asphaltene molecule constructed with XMOL is characterized by a globular shape. This globular shape of asphaltenes should not be considered as the stable form in oil. It is directly related to the construction procedure involved in the program to generate the 3-D graphic display. By selecting another volume of construction (flatter parallelepipedic for instance), the resulting shape would have changed. As a matter of fact, the 3-D molecular model represented here does not take into account the electronic constraints imposed by the environment (solvent, resins, aromatics, saturated hydrocarbons) which is a determinant factor. The XMOL (42) Kasrai, M.; Michael Bancroft, G.; Brunner, R. W., et al. Geochim. Cosmochim. Acta 1994, 58(13), 2865-2872.



an average error less than 3.7% for 100 C) which are listed in Table 7. These solutions have been obtained for structures containing between 100 and 1000 carbon atoms. The biggest error is due to the dipolar dephasing (average error per 100 C around 20%). This error has to be considered as acceptable (within the range of experimental error). Indeed, the dipolar dephasing value was determined using liquid 13C NMR (Table 1) which is less sensitive than the measurements performed using solid 13C NMR for asphaltenes. Figure 2 shows the relative error % on atomic ratios obtained for these 10 solutions; the minimum error % is obtained for the solution number 2. Therefore, the structure matching at best the analytical data is structure number 2 (molecular formula C647H756O11N12S20) and comprises 83 molecular fragments and 118 intergroup bonds. Its chemical structure is detailed in Table 8. Globally, the SIGNATURE solution is very close to the one obtained using XMOL software (Table 9) but the random selection of molecular groups is different. This difference can be explained as follows. First, the two libraries are not exactly similar, on the one hand regarding the molecular groups content, and on the other hand regarding the possible linkage between these fragments. With XMOL software, only two bonds are allowed per molecular group. In the case of SIGNATURE software, up to eight bonds are allowed for several fragments (i.e., a4n0 has seven bonding sites). Second, heteroatoms-containing molecular groups selected using the prediction program of XMOL (i.e., nitrogen-containing aromatics, sulfurcontaining aromatics) have to be manually introduced after obtaining the 3-D representation of the molecule. Each aromatic carbon replaced by an heteroatom must then be reincluded into a new aromatic group. Therefore, the molecular groups defined in the construction file (Table 5) have to be modified, specifically in the aromatic range. In our case, 12 nitrogen- and 15 sulfurcontaining aromatics (12 aroN and 15 aroS, respectively) were redistributed in a total of 27 aromatic carbons. Then the prediction program was applied to determine the nature and place of bonds between the fragments for this solution number 2. The result is the 3-D representation drawn in Figure 3. One can notice in this figure that some bonds have unrealistic lengths because the program SIGNATURE does not check the length of the cross-link bonds. With high cross-linked structures, such as our model, direct minimization leads to a conformation that is highly distorted. Therefore, we used SIGNATURE in conjunction with INSIGHT II/ DISCOVER. Starting from the Boscan molecule represented in Figure 3, we rebuilt a new molecule using the 3-D fragment library module of INSIGHT II. We divided the initial macromolecule into the largest possible fragments. This ensured the compatibility with

Table 5. Detailed Structure for the Best Solution of XMOL for BOSCAN Asphaltenea parameter

model

molecular fragment sn1 sn2 sn3 sn4 sn5 a1n0 a1n1 a1n2 a2n0 a2n1 a2n2 a3n1 a3n2 a4n0 sl2 sl9 sl10 sl15 si9 si21 si22 phenol interfragment bond L0 L4 L7 aroN aroS PDS S a

12 6 4 2 1 2 1 2 2 1 1 1 2 3 3 2 1 1 1 1 1 2 22 4 6 12 15 2 1

For abbreviations, please refer to the Glossary.

software does not calculate any energy minimization and thus the stability of the macromolecule which has been constructed is not checked. In order to obtain minimum energy conformations (most stable conformation), application of molecular mechanics and dynamics simulations has to be performed but XMOL is not compatible with other molecular modeling software. Therefore, we used SIGNATURE which is compatible with other molecular modeling software, XMOL having been used first in order to quantify the molecular groups which constitute the asphaltene unit. The input file of SIGNATURE was similar to the construction file given by XMOL and only some features were improved. On the one hand, the sulfur-containing aromatic molecular groups (15 aro-S) have been exactly defined as benzothiophenes (bt), dibenzothiophenes (dbt), dithiophenes (s2a), and naphthenothiophenes (Sa3n2) according to the quantitative data obtained by Fisher group type analysis. On the other hand, we introduced as molecular groups for nitrogen- and oxygen-containing aromatics respectively pyridine (na1), benzopyridines (na2) and oxygen-containing diaromatics (oa2). After running the elucidation program of SIGNATURE, we found 10 solutions in agreement with the analytical data (with

Table 6. Comparison between Experimental and Theoretical Data formula (XMOL)

MMw (XMOL)

MMp (GPC)

C618H722O10N12S20

9106

10000

atomic ratios O/Ca

H/C

N/Ca

S/Ca

aromaticity %

expt

XMOL

expt

XMOL

expt

XMOL

expt

XMOL

expt

XMOL

1.17

1.17

1.69

1.62

1.86

1.94

3.21

3.23

39.6

39.6

a

Multiplied by 100.

104


Kowalewski et al.

Figure 1. 3D representation of BOSCAN asphaltene obtained using XMOL software. Table 7. Solutions Obtained by Signature Equation for BOSCAN Asphaltene Containing between 100 and 1000 Carbon Atoms (the Best Solution is Solution No. 2) solution no. 1 2 3 4 5 6 7 8 9 10

molecular formula C637H745O10N12S20 C647H756O11N12S20 C659H770O11N12S22 C661H773O11N12S22 C663H775O11N13S22 C671H784O11N12S22 C674H788O12N13S23 C692H809O12N12S22 C723H845O12N12S24 C753H880O12N13S24

error per 100 C atoms (atoms)a

dipolar dephasing error per 100 C

calcd Mw

3.04(0.19)b

20.0 20.5 19.0 20.0 17.0 19.0 19.0 17.0 18.0 19.0

9357 9504 9726 9753 9793 9884 9986 10177 10649 11058

3.07(0.18) 2.85(0.18) 3.03(0.21) 3.60(1.37) 3.20(0.51) 2.90(0.22) 3.00(0.65) 3.19(0.74) 3.07(0.35)

a The error is the average error between the solution and the experimental data entered in the Signature input file. The solution no. 2 is the lowest-error solution. b The value of the error is underestimated not taking into account the dipolar dephasing.

the CVFF force field used in INSIGHT II/DISCOVER. The fragments were then connected using distance constraints. The charges were fixed to 0 and the potential types defined (potential type “op” was defined manually). Each fragment was defined as a subset and identified by different colors as represented in Figure 4. Periodically throughout the construction, due to some very high bond lengths, the various subsets were subjected to rapid minimization using conjugate gradients . In fact, some bond lengths were adjusted using mild distance constraints to avoid the distortion of the fragments during minimizations; then the minimized fragments were reconnected step by step. The geometry of the structure was gradually optimized according this

process. Once the structure completed, it was relaxed, the bond constraints (bond lengths constrained to 5 Å) were converted to bonds, and the structure was reminimized. In order to improve the model (by decreasing the total energy), the structure was modified manually by moving atoms in 3-D space, constraining or fixing them, and allowing the rest of the molecule to move during energy minimization. In this manner, bonds through rings were avoided and the chain positions optimized. This complete procedure is summarized in Figure 4. Last, the model was subjected to a final energy minimization using conjugate gradients to locate the minima of this conformation. The final structure



Figure 2. Relative error % on atomic ratios obtained for the SIGNATURE solutions ((Xtheoretical - Xexptl/Xexptl) × 100). Table 8. Detailed Structure for the Best Solution Obtained Using SIGNATURE for BOSCAN Molecule (Solution No. 2 of Table 7) parameter Signature for 100 C h o np s cp c cp(cpcph . ) molecular fragment a1n2 a1n3 a2n1 sn1 sn2 sn3 sn4 bt dbt s2 s2a Sa3n2 na1 na2 oa2 ester aro-o aro-oh sl1 sl2 sl3 sl5 sl9 sl11 sl12 sl15 sl16 sl23 si5 si10 sp19 H interfragment bond ali-s ali-o ali-h aro-csp2 aro-ali aro-aro a

model

error per 100 C atoms (analytical data - model data)

116.9 1.68 1.86 3.21 38.55 60.40 28.5

0.05 0.44 0.00 0.12 -0.02 0.43 20.5

1 2 3 7 3 3 1 3 1 2 5 2 8 4 2 3 2 1 3 4 3 4 1 2 1 1 1 1 5 1 1 2

3.07a

4 5 0 3 105 1

Total average error for 100 C.

is shown in Figure 5. We found a conformation with a total energy of about 7100 kcal/mol or 4.9 kcal/(atom mol).

Figure 3. 3D representation of BOSCAN asphaltene obtained using SIGNATURE software. Table 9. Comparison between the Best Solutions Obtained from XMOL and SIGNATURE Software software

formula

calcd Mw

XMOL SIGNATURE

C618H722O10N12S20 C647H756O11N12S20

9106 9504

After this construction step, MD using QAD technique was employed to study the closest energy minimum, starting from the final conformation obtained previously (Figure 5). The curve represented in Figure 6 displays the QAD simulation. At the beginning of the simulation, the total energy at an initial temperature of 1100 K was increased to about 16 000 kcal/mol. Then by freezing the molecule to 1 K, the total energy for the most stable conformation was decreased to about 7000 kcal/mol or 4.8 kcal/(atom mol). It means that the energy remained relatively constant as the structure moved between more or less equivalent conformations. Superimposing the different structures saved during the decrease of temperature, we can study the motions of the molecular fragments during the dynamic simulation. As we could expect, the core of the molecule is relatively rigid compared to the surrounding free aliphatic chains due probably to the high reticulation of the system obtained (high cross-link density). The fact that there was little improvement in the final energy when using MD at 1100 K and then decreasing to 1 K tells us that we were already near the bottom of the current energy well when we started the run. Because only few changes were found in the conformation during this MD experiment, we assume that we probably did not provide enough energy to the molecule; the starting temperature being too low, the molecule could not escape from the current energy well, possibly due to energy barriers resulting from the structure of the molecule. Another explanation is that we did not allow the simulation to run for times long enough to have the opportunity to escape the current energy well. More tests of MD simulations with higher starting temperature and longer time would be necessary to check these assumptions. But this may be due to an unrealistic model with too high cross-linking and in this case MD simulations using higher temperature and longer time should yield same results. To decrease this cross-linking, work is in progress to create a new library allowing only two or three bonds per molecular group.

106


Kowalewski et al.

Figure 4. Rebuilding and minimization of BOSCAN using in conjunction INSIGHT II/DISCOVER. (Each number corresponds to a fragment of the total molecule which is separately minimized. For the total molecule, each fragment is identified by different colors.)

Figure 6. QAD simulation.

Figure 5. Minimized BOSCAN asphaltene “unit”.

Conclusion To the best of our knowledge, this paper is the first to use molecular modeling techniques to study asphaltene structures. This effort for setting up a protocol of asphaltene molecular modeling has been reasonably successful on the BOSCAN sample. The procedure involves several steps including complementary software. The first step is the gathering of an analytical data set on asphaltenes. The second step is the use of the XMOL software to verify the consistency of the data set and to have a first 3-D asphaltene representation. Then, the SIGNATURE software generates construction

files compatible with INSIGHT II/DISCOVER, which allow to find a stable conformation with relatively low energy and to study the most stable(s) conformation(s). Due certainly to the reticulation of the model found for this asphaltene, there are probably very few possible conformations and MD simulations do not improve the structure substantially. Therefore, an attempt to decrease this high reticulation (by decreasing the number of bonds allowed per molecular group) is in progress. However, to improve the realism and the rigor of our model, the perspective is to take into account a solvent environment. This should help for the understanding of the physicochemical behavior of asphaltenes in oil production conditions. Work is in progress to examine the structural fluctuations of a solvated asphaltene using periodic boundary conditions (PBC) which refers to the simulation of the interactions of a molecular system in a periodic lattice of identical subunits. Benzofuran has been chosen as the model for simulating resin


interactions. By application of such conditions, the influence of the environment can be investigated. Acknowledgment. We wish to express our gratitude to Biosym Technologies for its assistance and the free use of INSIGHT II/DISCOVER software. More specifically, we thank M. K. Seiti from Biosym. The contributions of IFP Geochemistry colleagues are gratefully acknowledged for the gathering of the analytical data set. We thank D. Espinat and J. C. Roussel for the results of SAXS-SANS measurements and 13C NMR, respectively. Glossary sn1 sn2 sn3 sn4 sn5 a1n0 a1n1 a1n2 a2n0 a2n1 a2n2 a3n1 a3n2 a4n0 sl1 sl2 sl3 sl5 sl9 sl10

naphthalene with 1 cycle naphthalene with 2 cycles naphthalene with 3 cycles naphthalene with 4 cycles naphthalene with 5 cycles benzene naphtheno (1)-aromatic (1) naphtheno (2)-aromatic (1) naphthene naphtheno (1)-aromatic (2) naphtheno (2)-aromatic (2) naphtheno (1)-aromatic (3) naphtheno (2)-aromatic (3) pyrene methane ethane n-propane (n-C3) n-pentane (n-C5) n-nonane (n-C9) n-decane (n-C10)

Energy & Fuels, Vol. 10, No. 1, 1996 107 sl11 sl12 sl15 sl16 sl23 si5 si9 si10 si21 si22 sp19 H phenol or aro-oh aro-o L0 L4 L7 aroN na1 na2 aroS bt dbt s2a Sa3n2 PDS or s2 S h o np s cp c cp(cpcph . )

n-C11 n-C12 n-C15 n-C16 n-C23 isopentane (i-C5) isononane (i-C9) i-C10 i-C21 i-C22 pristane hydrogen phenol ether-benzene molecular groups linked by 0C molecular groups linked by 4C molecular groups linked by 7C nitrogen-containing aromatic pyridine carbazole sulfur-containing aromatic benzothiophene dibenzothiophene disulfur-containing aromatic naphtheno (2)-sulfur-containing triaromatic disulfur bridge sulfur hydrogen bonded to C sp3 oxygen sp2 aromatic nitrogen sulfur sp2 aromatic carbon sp3 aliphatic carbon dipolar dephasing EF950106T

Preliminary Results on Molecular Modeling of Asphaltenes Using

Recommend Documents