CHEMICAL ENGINEERING TRANSACTIONS VOL. 57, 2017 A publication of The Italian Association of Chemical Engineering Online at www.aidic.it/cet Guest Editors: Sauro Pierucci, Jiří Jaromír Klemeš, Laura Piazza, Serafim Bakalis Copyright © 2017, AIDIC Servizi S.r.l. ISBN 978-88-95608- 48-8; ISSN 2283-9216 Molecular Representation of Molecular Distillation Cuts of Vacuum Residue by Spectrometry Ultra-High Resolution and Conventional Analytic Claudia X. Ramíreza*, Juan E. Torresc*, Diana Catalina Palacio Lozanob, Juan P. Arenas-Diazb, Enrique Mejia-Ospinob, Viatcheslav Kafarova, Alexander Guzmánc a Faculty of Physicochemistry, Universidad Industrial de Santander- School of Chemical Engineering, Carrera 27 9th street, Bucaramanga-Santander. A,A. 678, Colombia b Faculty of Sciences, Universidad Industrial de Santander- School of Chemical Engineering, Carrera 27 9th street, Bucaramanga-Santander. A,A. 678, Colombia c Ecopetrol S.A., Instituto Colombiano del Petróleo, kilometro 7 via Piedecuesta, Piedecuesta- Santander. A,A. 4185, Colombia. claudia.ramirez2@correo.uis.edu.co, juan.torres@ecopetrol.com.co A set of 200 representative molecules for a Vacuum residue (VR) fractionated with molecular distillation (MD) was obtained. A built in-house wiped film molecular distillation unit that operates in continuous mode at pressures below 1E-5 Pa was used to fractionate a Colombian heavy VR, four cuts were obtained changing the temperature of the system. First, each cut was characterized by conventional analytical, using conventional standardized methodologies (ASTM Methods). Second, each cut was characterized by ultra-high mass resolution spectrometry, Fourier transform ion cyclotron resonance mass spectrometry (FT-ICR MS), combination of positive electrospray ionization ((+)ESI), negative electrospray ionization ((−)ESI) and positive atmospheric pressure photoionization ((+)APPI). Information obtained was employed for reconstruction and generation of set molecules that represents the original vacuum residue by contribution methods. A multi-feed strategy was used; it was based on a same set of representative molecules for each cut of MD but changing its composition among them. Optimization of molecular composition was obtained with simulated annealing. The error in estimation of bulk properties was lower than five percent when the FT-ICR-MS data was incorporated as input information. 1. Introduction Nowadays, light crude oil reserves are declining and the future production scenarios show a higher contribution of unconventional heavy oils (Charry Sanchez & Betancourt Torcat, 2016). This type of feeds has a lack of molecular detail in its characterization due to the large number of molecules that make it up (Read, 1976). That lack of detail is the driving force to make research and to develop projects on characterization and upgrading technologies of heavy and extra-heavy crude oils. Standard characterization techniques applied to heavy oils cuts provide a set of physicochemical properties such as TBP ("true boiling point"), density, refractive index (IR), distribution of aromatic molecules, % carbon (%C), % Hydrogen (%H), % Conradson (CCR) Carbon and molecular weight -WT, SARA analysis (Saturates, Aromatics, Resins and asphaltenes), among others. However, these properties are not sufficient for proper characterization due to the complexity of the components in vacuum residua (Poveda, Molina, & Agreda, 2014). A modern technique applied in the study of heavy crude and residues is Fourier transform ion cyclotron resonance mass spectrometry (FT-ICR-MS) which is a complementary technique. This is the most powerful technique that can provide a more detailed and vast information on chemical composition of very heavy and complex hydrocarbons (Marshall & Rodgers, 2004), information as carbon number, double bound equivalent (DBE), chemical formula and distribution of components. DOI: 10.3303/CET1757179 Please cite this article as: Ramirez C.X., Torres J.E., Palacio Lozano D.C., Arenas-Diaz J.P., Mejia-Ospino E., Kafarov V., Guzman A., 2017, Molecular representation of molecular distillation cuts of vacuum residue by spectrometry ultra-high resolution and conventional analytic, Chemical Engineering Transactions, 57, 1069-1074 DOI: 10.3303/CET1757179 1069 mailto:claudia.ramirez2@correo.uis.edu.co Due heavy oils have molecules of high molecular mass and low volatility, new highly efficient methods for the separation of these type complex feeds has been investigated, developed and represented. The supercritical fluid extraction and fractionation (SFEF) method and molecular distillation (DM) ( Orrego-Ruiz, Mejía-Ospino & Guzmán, 2015) are techniques that allow separation of the heaviest components without reaction in the system. The representation of molecules present into of heavy oil has been performed by Neurock and others (1994), represented a residue heavy by structural attributes; number of ring aromatics, number of side chains, number of carbons cata-condensed and peri-condensed, etc. Each attribute can be represented by a probability distribution function (PDF), which is adjusted according to the error of bulk properties of characterized feed. The set of molecules generated was obtained using two algorithms. First algorithm is Stochastic reconstruction allows generate the set of equimolar structures. From structures, properties of each molecule are performed by group contribution and subsequently mixing properties are obtained. Simulated annealing was the second algorithm employed for adjustment of the molecules compositions. In more recent years Zang and others (2014) used information obtained by the fragmentation technique collision induced dissociation (CID) in FT-ICR MS (carbon number and number of aromatic rings) as restrictions into the model of representation structural developed by them. The configuration of cores obtained in this work was used as a comparative parameter but not as a restrictive parameter. In this work is proposed as methodology for obtaining a set representative of molecules, first the VR fractionation in four cuts by DM. Second, the analysis of each cut by ultra-high resolution mass spectrometry (APPI+/ ESI+/ESI- FT-ICR MS) and conventional analytics. Finally, with a characterization detailed, a set of representative structures was obtained through representation for attributes, Incorporating information of FT- ICR-MS as the number maximum of carbons, aromatic and naphthenic rings. Also, from this information were provide the PDF's for aromatics, resins, nitrogen compounds and sulfur compounds in molecular reconstruction-MR (reconstruction by contribution methods) as input parameters. All molecules generating will be used for development of process models. 2. Materials and Characterization A Colombian vacuum residue (354+ °C) was fractionated using a built in-house wiped film molecular distillation unit that operates in continuous mode at pressures below 1E-5 Pa. Four cuts were obtained by changing the temperature of the system: a non-distillable sample, or residue 687 +°C AET (atmospheric equivalent temperature), and three distillates cuts with boiling points from IBP-603, 603-645 and 645-687 °C AET, . Each cut obtained was characterized with conventional standardized methods and unconventional analytical techniques. Table 1. Conventional analytic for cuts obtained from molecular distillation of Colombian vacuum residue. Colombian Vacuum Residue 1 IBP-603 603-645 645-687 687+ Density, ρ (g/cm3) 1,0138 1,0096 1,0428 1,0758 Saturates (wt fraction) 0,16 0,13 0,06 0,02 Aromatics (wt fraction) 0,68 0,70 0,52 0,23 Resins (wt fraction) 0,16 0,16 0,32 0,40 Asphaltenes (wt fraction) 0,01 0,01 0,11 0,36 Carbon, (wt fraction) 0,8657 0,8531 0,8592 0,8214 Sulfur (wt fraction) 0,0373 0,0369 0,04 0,0404 Nitrogen (Fraction) 0,0012 0,0014 0,0019 0,0027 RCC (wt fraction) 5,53 9,04 22,9 38,9 Waxes (wt fraction) 0,003 0,005 0,007 0,025 MW (g/mol) 370 446 631 1207 TBP, 0.5 % wt (°C) 409 496 551 576 5 % wt 447 523 575 626 50 % wt 530 570 638 757 90 % wt 582 602 804 865 99 % wt 654 645 898 926 2.1. Conventional standardized methods Four cuts obtained from molecular distillation unit were characterized with bulk properties: density was measured by means of a digital densimeter (ASTM D4052); molecular weight by vapor pressure osmometry 1070 (VPO) (ASTM 2503); C %wt, H %wt, N %wt, S %wt by elemental analysis (ASTM D5291); saturates, aromatics, resins and asphaltenes weight fraction by SARA compositional analysis (ASTM D2007); TBP by Simdis (ASTM D1160/ ASTM D7169); waxes content by wax content determination (UOP 46); sulfur content (ASTM D1552); total nitrogen (ASTM D3228). Properties for these cuts are given in Table 1. According to information obtained by TBP, the initial point for all feeds is over 400 °C, this indicate that current structures could have more than 20 carbon atoms. The final points were calculated through probability extrapolation method. 2.2. Nonconventional analytical techniques The nonconventional analytic used for characterization was ultra-high mass resolution spectrometry, a combination of positive electrospray ionization ((+) ESI), negative electrospray ionization ((−) ESI) and positive atmospheric pressure photoionization ((+) APPI) Fourier transform ion cyclotron resonance mass spectrometry (FT-ICR MS). FT-ICR MS analyses were performed using a 15 T SolariX FT-ICR mass spectrometer of Bruker Daltonics (Billerica, MA). Nitrogen was used as the drying and nebulizing gas. All prepared samples were directly injected with a syringe pump. For APPI, samples were diluted to 0.1 mg/mL in toluene. The nebulizing temperature was 300°C at a constant pressure of 0.7 bar. The drying temperature was 200°C at a flow rate of 4.5 L/min. An ion accumulation time in the collision cell of 0.010 s was set, which was followed by a time-of-flight of 0.1 s to transfer the ions to the ICR cell. For ESI, samples were diluted in 40:60 (v/v) toluene/methanol to yield a 0.1 mg/mL solution. Acetic acid was added (10 µL to each 1 mL of sample solution) to facilitate protonation of basic compounds by positive ion mode, and for negative ion mode, ammonium hydroxide was added (10 µL to each 1 mL of sample solution) to promote deprotonation of acidic compounds. In all spectra, time-domain transient signals were accumulated and averaged (100 scans) to enhanced the signal-to-noise ratio of each spectrum in four mega data size. The data obtained by ultra-high mass spectrometry was processed and analyzed, getting a set of probability distribution parameters employed as input parameters in the molecular reconstruction. The aim of using multiple ionization methods was to provide the greater information of composition of heavy feeds. 0 10 20 30 40 0 10 20 30 40 0 10 20 30 40 10 20 30 40 50 60 70 80 90 100 0 10 20 30 40 10 20 30 40 50 60 70 80 90 100 10 20 30 40 50 60 70 80 90 100 10 20 30 40 50 60 70 80 90 100 DB E Class HC/ IBP-603 VR1 /APPI Class HC/ 603-645 VR1 /APPI Class HC/ 645-687 VR1 /APPI Class HC/ 687+ VR1 /APPI Class S1/ IBP-603 VR1 /APPI DB E Class S1/ 603-645 VR1 /APPI Class S1/ 645-687 VR1 /APPI Class S1/ 687+ VR1 /APPI Class N1/ IBP-603 VR1 /ESI+ DB E Class N1/ 603-645 VR1 /ESI+ Class N1/ 645-687 VR1 /ESI+ Class N1/ 687+ VR1 /ESI+ Class N1/ IBP-603 VR1 /ESI- DB E Number of Carbons Class N1/ 603-645 VR1 /ESI- Number of Carbons Class N1/ 645-687 VR1 /ESI- Number of Carbons Class N1/ 687+ VR1 /ESI- Number of Carbons Figure 1: Distributions obtained for APPI/ ESI(+)/ESI(-) FT-ICR MS Distributions of aromatics classes (HC) and sulfur compound (S1) are showed in Figure 1, these compounds are detected through ionization method (+) APPI; besides, distributions of nitrogen compounds (N1) are detected according to basicity: pyridines through ESI(+) while pyrroles with ESI(-). These classes were used for obtaining the distribution parameters and constraints used in the generation of molecules in VR molecular reconstruction. In Figure 1 also is showed the relation between aromaticity (DBE) and carbon number for each detected molecule. The carbon number interval for all cuts is between 10-80; as well as aromatic rings number can take values less than 11. Based on the information obtained, the distributions of lighter cuts have 1071 a left bias while heavier cut (687+) presents a right bias, which indicates structures with greater carbon number when the cut becomes heavier. Moreover, if DBE value is constant and carbon number increases the aliphatic component of molecule increase. Otherwise, if carbon number is constant and DBE increases, indicate structure with more aromaticity and a less aliphatic component. A set of parameters and restrictions were determined with FT-ICR MS information: the probability density functions for aromatic and resin compounds; sulfur distribution, and nitrogen distribution. The maximum number of carbons per molecule was 80 and the maximum number of aromatics rings were constrained to ten. 3. Structural Attributes Modeling The basis for obtaining a set of representative structures is the global representation through attributes. Any structure present in the crude could be represented by its structural attributes, such as number of ring aromatics and number of cata-condensed carbons, etc. Each attribute in mixture can be represented by a probability distribution function (PDF). The PDF provides a quantitative probability of finding a structural attribute. A set of PDF could represent structural attributes considered in the molecule generation and the values of attributes for each structure: besides, it can provide properly the identification of a molecule (chemical rules) and finally to set structures that depict the feed. Therefore, each density function employed for representing a structural attribute allows to move the occurrence probability of events of a random experience to a numeric characteristic. According to the literature, the gamma distribution is the best one that represents the structural attributes of heavy oil cuts, while a combination of gamma, exponential, and chi-square distributions represents better light cuts (Klein el al., 2005). For molecular reconstruction is required to know the PDF parameters but these are unknown, and then these are assumed at the first time and later are calculated again to achieve a minimum error in the system. Molecular reconstruction of feeds from the information obtained by the characterization techniques uses Monte Carlo simulation, which is a quantitative technique that uses statistics and computers to mimic the behavior of non-dynamic random real systems (Faulín & John, 2012). This technique is much more accurate when more number of random components are generated (molecules). Once the molecules are represented by appropriate structural attributes, the PDF corresponding to these molecular attributes need to be optimized to match the experimental data. In this case was obtained the first PDF by information ESI+/ESI-/APPI+ FT-ICR MS and was implemented as constraints in the objective function.When a representative set of molecules is achieved, the next step is to estimate its properties. This calculation was assessed by contributions methods. The structural attributes considered in this work are CH, CH2, CH3 in paraffinic chains; number of naphthenic rings; CH2 and CH naphthenic; number of aromatic rings; aromatic carbon type; cata-condensed and peri-condensed carbon; number of sheets in asphaltenes; -SH, aromatic S (thiophene); aromatic N (pyridine and pyrrole) and -NH2. The objective function for both optimization algorithms is the same (Equation 1); in the first algorithm, the error between the predicted and calculated properties is minimized when changing the structures that represent the system. The second algorithm minimizes the error of the predicted and experimental properties, but changing the composition of these in the sample subject to the structures obtained in the first algorithm. The approach of nonlinear optimization of molecular reconstruction is as follows: 𝑂𝑏𝑗𝑒𝑡𝑖𝑣𝑒 𝐹𝑢𝑛𝑐𝑡𝑖𝑜𝑛 (𝑃) = 𝑀𝑖𝑛𝑖𝑚𝑖𝑧𝑒 = ( ∑ (𝑇𝐵𝑃𝑖𝑒𝑥𝑝−𝑇𝐵𝑃𝑖𝑐𝑎𝑙𝑐) ≠𝑓𝑟𝑎𝑐𝑡𝑖𝑜𝑛𝑠 𝑖=1 𝐴 ) + ( (𝜌𝑒𝑥𝑝−𝜌𝑐𝑎𝑙𝑐) 𝐵 ) + ( (𝑀𝑊𝑒𝑥𝑝−𝑀𝑊𝑐𝑎𝑙𝑐) 𝐶 ) + ( (𝑊𝑎𝑥𝑒𝑥𝑝−𝑊𝑎𝑥𝑐𝑎𝑙𝑐) 𝐷 ) + ( ∑ (𝑆𝐴𝑅𝐴𝑖 𝑒𝑥𝑝−𝑆𝐴𝑅𝐴𝑖 𝑐𝑎𝑙𝑐) ≠𝐶𝑜𝑚𝑝𝑜𝑛𝑒𝑛𝑡𝑠 𝑖=1 𝐸 ) + ( (%𝑆𝑒𝑥𝑝−%𝑆𝑐𝑎𝑙𝑐) 𝐻 ) + ( (𝑁𝑒𝑥𝑝𝑏á𝑠𝑖𝑐 −𝑁𝑐𝑎𝑙𝑐𝑏á𝑠𝑖𝑐) 𝐼 ) + ( (𝑁𝑒𝑥𝑝𝑇𝑜𝑡𝑎𝑙 −𝑁𝑐𝑎𝑙𝑐𝑇𝑜𝑡𝑎𝑙) 𝐽 ) + ( (%𝐶𝑒𝑥𝑝−%𝐶𝑐𝑎𝑙𝑐) 𝐻 ) + ( (𝑅𝐶𝐶𝑒𝑥𝑝−𝑅𝐶𝐶𝐶𝑐𝑎𝑙𝑐) 𝐻 ) + ( ∑ (𝑃𝐷𝐹𝑖𝑒𝑥𝑝−𝑃𝐷𝐹𝑖𝑐𝑎𝑙𝑐) 𝑀 𝑃𝑎𝑟𝑎𝑚𝑒𝑡𝑒𝑟𝑠 𝑖=1 𝐾 ) (1) The constants A, B,C ,D ,E ,F ,G ,H ,I ,J ,H and K present in objective function are weigh parameters and Subscripts exp correspond to experimental values and calc to calculate values. 3. Molecular Reconstruction of Molecular Distillation Cuts In this work was obtained a set of 200 molecules that mimics the properties a Colombian VR. Parameters of Aromatic distribution were calculated from each class associate (HC, N1, S1) and distribution parameters of gamma PDF (α, β) were adjusted by FT-ICR MS information; in the Table 2 is presented the results for all evaluated cuts. Each pair of parameters for each distribution was obtained for a specific ionization method. The PDF for non-basic nitrogen compounds was obtained from ESI-/FT-ICR MS and the PDF for basic nitrogen compounds from ESI-/FT-ICR MS and the PDF for sulfur compounds, aromatic distributions and aromatic distributions in resins from APPI+/FT-ICR MS. 1072 In molecular reconstruction can be used up to 76 distribution parameters, of which the 18 first ones were used in generation of saturated and aromatic structures and resins. The following parameters were used for generation of asphaltenes, with the constraint of maximum seven sheets present in the archipelago asphaltenes, 10 aromatics rings present in any structure aromatic and a maximum of 80 carbons presents in each structure. Each continental or archipelago asphaltene have its own set of distributions and may be different from the other. Table 2. Parameter of PDF obtained by FT-ICR MS Colombian Vacuum Residue 1 Parameters IBP- 603 603- 645 645- 687 687+ Parameters IBP- 603 603- 645 645- 687 687+ α- Non Basic Nitrogen Compounds 6,1 5,7 5,7 5 β- Sulfur Compounds 0,27 0,29 0,3 0,36 β- Non Basic Nitrogen Compounds 0,15 0,17 0,17 0,25 α- Aromatic rings 4,2 4,3 4,5 4 α- Basic Nitrogen compounds 5,6 4,9 4,4 4,3 β- Aromatic rings 0,27 0,29 0,3 0,36 β- Basic Nitrogen compounds 0,17 0,2 0,24 0,28 α- Aromatic rings in Resins 6,7 5,5 4,5 3,85 α- Sulfur Compounds 4,2 4,3 4,5 4 β- Aromatic rings in Resins 0,47 0,72 0,95 1,11 4. Results and discussion Molecular reconstruction of VR feed was made from the information obtained by molecular distillation and conventional and non-conventional analytical methods. The multi-feed option was used for holding the same set of molecules in each molecular distillation cut. Between cuts have common structures themselves (it was observed with FT-ICR MS) but the difference among cuts is its composition (Figure 2). A set of molecules was obtained with adjusting of calculated bulk properties, the error between experimental and calculate properties is shown in Figure 3. These model molecules can characterize the VR feed due to its properties and structure, which are coherent with its conventional and non-conventional analytics. Distributions of molecules for classes (HC, S1, and N1 by molecular reconstruction are shown in Figure 4. 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95100 0,00 0,05 0,10 0,15 Co m po sit io n (F ra ct io n) Estructures Number IBP-603 603-645 645-687 687+ Figure 2. Mass fraction of model molecules and error of bulk properties Figure 3. Mass fraction of model molecules and error of bulk properties Distributions obtained with molecular reconstruction are compared with distributions obtained for APPI/ ESI(+)/ESI(-) FT-ICR MS, both were depict in function of DBE and carbon number but the first takes into account composition while the second takes into account the ionization intensity. As shown in Figure 4, a general behavior of the obtained structures is found, an increase in the number of carbons in the molecules as the DM cut becomes heavier. The DBE in nitrogenous classes have a higher degree of polydipersity, the same maner this class present an increase in the number of carbons towards the heaviest cut. Comparing the structures of classes HC, S1 and N1 with the data obtained by FT-ICR MS can be found that the molecules generated by molecular reconstruction are within the structures obtained with greater intensity by FT-ICR MS. The composition behavior by cuts with the same set of molecules tends to be similar among the lighter ones (IBP-603 and 603-645). 0 1 2 3 4 5 6 E rr o r (% ) IBP-603 603-645 645-687 687+ 1073 0 10 20 30 40 0 10 20 30 40 0 10 20 30 40 10 20 30 40 50 60 70 80 90 100 0 10 20 30 40 10 20 30 40 50 60 70 80 90 100 10 20 30 40 50 60 70 80 90 100 10 20 30 40 50 60 70 80 90 100 D B E Class HC IBP-603 MR Class HC 603-645 MR Class HC 645-687 MR Class HC 687+ MR D B E Class S1 IBP-603 MR Class S1 603-645 MR Class S1 645-687 MR Class S1 687+ MR D B E Class N1 Pyrrole IBP-603 MR Class N1 Pyrrole 603-645 MR Class N1 Pyrrole 645-687 MR Class N1 Pyrrole 687+ MR Number of Carbons Class N1 Basic IBP-603 MR D B E Class N1 Basic 603-645 MR Number of Carbons Class N1 Basic 645-687 MR Number of Carbons Number of Carbons Class N1 Basic 687 + MR Figure 4. Distributions obtained with molecular reconstruction 5. Conclusions As result, a set of 200 molecules that mimics the properties of molecular distillation cuts of Colombian VR with an average error between the calculated and experimental properties no greater than 5% was obtained. The set of molecules obtained will be used in the development of simulation models in commercial simulators. The implemented multi-feed strategy of molecular reconstruction could represent four molecular distillation cuts of Colombian VR with different characteristics by the same set of molecules, the difference between these were their composition. From the obtained information through ionization methods with FT-ICR MS applied in feeds the aromaticity constraints and limits in carbon number were included in molecular reconstruction model, which was very useful in generation of asphaltenes structures. The type of asphaltene obtained through molecular reconstruction for the VR feed was continental. Acknowledgments The development of this work was supported by the Research Agreement No 5211770 Universidad Industrial de Santander – ICP (ECOPETROL S.A.) and COLCIENCIAS under the project No. FP44842-039-2015. Juan Arenas-Diaz and Diana Palacio thank to COLCIENCIAS for the Graduate Fellowships. References American Society for Testing and Materials, ASTM, 2012 ,Annual Book of ASTM Standards 2012. Astm. Charry-Sanchez J., & Betancourt-Torcat A., 2016, Environmental and Economics Trade-Offs for the Optimal Design of a Bitumen Upgrading Plant. Industrial & Engineering Chemistry Research, 55, 11996−12013. Faulín J., & Juan A., 2012, Simulación de Monte Carlo con Excel. Mexico. Klein M. T., Bertolacini R. J., Broadbelt L. J., & Group F., 2005, Molecular Modeling in Heavy Hydrocarbon Conversions. Molecular Modeling in Heavy Hydrocarbon Conversions. CRC Press. Marshall A. G., & Rodgers R. P., 2004, Petroleomics: The Next Grand Challenge for Chemical Analysis. Accounts of Chemical Research, 37, 53-59. Neurock M., Nigam A., Trauth D., & Klein M., 1994, Molecular Representation of Complex Hydrocarbon Feedstocks Through Efficient Characterization and Stochastic Algorithms. Chemical Engineering Science, 49, 4153-4177. Orrego-Ruiz J. A., Molina D., Mejía-Ospino E., & Guzmán A, 2015, Understanding the Molecular Information Contained in the Infrared Spectra of Colombian Vacuum Residua by Principal Component Analysis. In Analytical Methods in Petroleum Upstream Applications, 275-300, CRC Press. Poveda J. C., Molina D. R., & Agreda E. P., 2014, 1H- and 13C-NMR structural characterization of asphaltenes from vacuum residua modified by thermal cracking. CT&F-Ciencia, Tecnología y Futuro, 5, 49-59. Read R.C., 1976, The enumeration of acyclic chemical compounds. Chemical Applications of Graph Theory, 25- 61. Zhang L., Hou Z., Horton S. R., Klein M. T., Shi Q., Zhao S., & Xu C., 2014, Molecular representation of petroleum vacuum resid, Energy and Fuels, 28(3), 1736–1749. 1074