Aggregate Structure in Heavy Crude Oil: Using a Dissipative Particle

Jul 21, 2010 - The alkyl chains are represented by DPD particles and strings connecting ..... One thing we should bear in mind there is no specificati...
0 downloads 0 Views 7MB Size
Download PDF
Energy Fuels 2010, 24, 4312–4326 Published on Web 07/21/2010

: DOI:10.1021/ef1003446

Aggregate Structure in Heavy Crude Oil: Using a Dissipative Particle Dynamics Based Mesoscale Platform Sheng-Fei Zhang,†,‡ Li-Li Sun,†,‡ Jun-Bo Xu,† Hao Wu,† and Hao Wen*,† †

State Key Laboratory of Multi-Phase Complex System, Institute of Process Engineering, Chinese Academy of Sciences, Beijing 100190, P. R. China, and ‡Graduate University of Chinese Academy of Sciences, Beijing 100049, P. R. China Received March 22, 2010. Revised Manuscript Received June 28, 2010

Heavy crude oil consists of thousands of compounds and much of them have a fairly large relative molar mass and complex structure. It is hard to learn the dynamic behavior of this fluid system at all atom models. The present study aims at constructing a mesoscale platform to explore aggregate behavior of asphaltenes in heavy crude oil. The aggregate structure in heavy crude oils was investigated by introducing rigid body fragments, which represents the significant presence of structures of fused aromatic rings in fractions such as asphaltenes and resins into dissipative particle dynamics (DPD). Another pressing task about how to determine the structure of the average model molecules and conservative force parameters was discussed in detail. With some regularity concerning the number of rings, the distribution of side chains and heteroatoms in average model molecules are revealed. Finally, we integrated the modified DPD program, model molecules, and the parameters selected for the preliminarily simulation of the heavy crude oil and emulsion system. The interlayer distance and the number of layers of the well-ordered structure in heavy crude oil are similar to some molecular dynamics works and supported by X-ray and transmission electron microscopy (TEM) experimental data. The relationship between the stability and the mass ratio among components of heavy crude oil is explored, and the result of our simulations fits the regularity Shell once published. In the emulsion system, the surfactant-like feature of asphaltenes and resins are observed. The preliminary simulation results demonstrate the validity of the rotational algorithm and parameters employed and encourage us to extend this platform to study the rheological and colloidal characteristics of heavy crude oils in the future. It is known that the main characteristics of heavy crude oil are directly related to the presence of aromatic compounds, such as asphaltenes and resins which have a relative molar mass as large as 1000.7 These macromolecules are a big challenge to the simulation technologies. A large amount of disputes exist on the actual contents and average molecular structures, since it is hard to get a full image from simply assembling the recognized structural fragments of the heavy fractions in heavy crude oils. However, a number of works have been done on the phenomena of self-aggregation for different crude oil systems at the molecular level. Rogel8 manifested that the stabilization energies of asphaltene and resin associates were mainly attributed to the van der Waals interactions among asphaltene and resin molecules by MD simulations. Three geometric piling structures of asphaltene self-association were observed by Pacheco-S anchez et al.9 using MD simulations in a monodispersed system of 96 hypothetical asphaltene molecules with nonperiodic boundary conditions. Takanohashi et al.10 found that temperature is important to the status of the hydrogen bond between asphaltene molecules and aromatic-aromatic stacking interactions in their heat-induced relaxation study of asphaltene aggregates. Considering the aggregation of model molecules of asphaltene, water, toluene, and heptane, Kuznicki et al.11

1. Introduction The research interest on heavy crude oil is increasing as the reservoir of conventional crude oil has decreased.1 The recovery and refining of heavy crude oil encounters various physicochemical phenomena concerning the aggregates formed by associated asphaltene molecules. The aggregate structure of asphaltene molecules in heavy crude oil influences the rheology, dispersity, stabilization of water-in-oil emulsions,2 poisoning of catalysts,3 fouling of hot metal surfaces, and extent of wettability alteration,4 and is a large concern of industry. Experimental investigation of the microstructure of heavy crude oil is proven to be difficult5 because of the strong dependence on the experimental conditions and also the absence of effective approaches.6 Fortunately, simulation technologies, such as molecular mechanics (MM), molecular dynamics (MD), Monte Carlo simulation (MC), Brown dynamics, dissipative particle dynamics (DPD), and MesoDyn, may provide details at the microscale level in advancing the knowledge in this field. *To whom correspondence should be addressed. (1) Marshall, A. G.; Rodgers, R. P. Proc. Natl. Acad. Sci. U.S.A. 2008, 105 (47), 18090–18095. (2) McLean, J. D.; Kilpatrick, P. K. J. Colloid Interface Sci. 1997, 196, 23–34. (3) Silverman, L. D.; Winkler, S.; Tiethof, J. A.; Witoshkin, A. Matrix Effects in Catalytic Cracking. In NPRA Annual Meeting, Los Angeles, CA, 1986. (4) Buckley, J. S.; Wang, J.; Creek, J. L. Solubility of the LeastSoluble Asphaltenes. In Asphaltenes, Heavy Oils, and Petroleomics; Springer: New York, 2007. (5) Merdrignac, I.; Espinat, D. Oil Gas Sci. Technol. 2007, 62, 7–32. (6) Kharrat, A. M.; Zacharia, J.; Cherian, V. J.; Anyatonwu, A. Energy Fuels 2007, 21, 3618–3621. r 2010 American Chemical Society

(7) Mullins, O. C. SPE J. 2008, 13, 48–57. (8) Rogel, E. Energy Fuels 2000, 14, 566–574. (9) Pacheco-Sanchez, J. H.; Zaragoza, I. P.; Martı´ nez-Magadan, J. M. Energy Fuels 2003, 17 (5), 1346–1355. (10) Takanohashi, T.; Sato, S.; Saito, I.; Tanaka, R. Energy Fuels 2003, 17, 135–139. (11) Kuznicki, T.; Masliyah, J. H.; Bhattacharjee, S. Energy Fuels 2008, 22, 2379–2389.

4312

pubs.acs.org/EF

Energy Fuels 2010, 24, 4312–4326

: DOI:10.1021/ef1003446

Zhang et al.

observed that uncharged molecules did not aggregate at the toluene-water interface, whereas charged terminal groups had an opposite tendency. He also concluded that the bottlenecks for MD resulted from the lack of unique suitable model molecules and force fields as well as a small length and time scale.12 To overcome the obstacle of adopting proper average model molecules representing compounds in heavy crude oil, Boek et al.13 developed an algorithm to generate quantitative molecular representations (QMRs) of asphaltenes based on the experimental data. This algorithm was used to study the nanoaggregates of one resin and two asphaltene structures by Headen et al.14 The formation of asphaltene multimers was observed in both toluene and heptane by calculating the average angle between polyaromatic planes as a function of separation, while some aggregates of resins could only be found in heptane. Apart from MD and MC simulations, Duda et al.15 developed an analytical solution of the associative OrnsteinZernike integral equation in the Percus-Yevick approximation to study the structural properties and phase separation of asphaltene in model oil. The relationship between instability and pair potential parameters, sizes, density, and composition of the system was revealed. Qin et al.16 indicated that the consideration of asphaltene precipitation led to a different prediction of oil recovery from their reservoir simulation results using solid models proposed by Nghiem et al.17 Many of the works on the aggregate structures, however, are limited to the modeling systems composed of a few asphaltene and resin molecules, which is difficult to apply to the whole scene of heavy crude oil. Adopting coarse-grained model molecules at the mesoscale level will keep the simulation from suffering these troubles. Aguilera-Mercado et al.18 observed the complex multimodal aggregation pattern of asphaltenes in their canonical MC simulation by treating asphaltenes as discotic seven-center Lennard-Jonesium molecules and resins as single spheres. However, their oversimplification of the structural features of asphaltenes and resins hinders its application to obtain further information about the microstructures in heavy crude oil. The dissipative particle dynamics (DPD), in which the model molecules could be decided according to the coarse-grained level, is a widely used technology at the mesoscale and can process a much larger spatial and temporal scale system than MD can. Applications have abounded in simulations of lipids,19 block copolymers,20 vesicle formation,21 and surfactants.22 It also exhibits the possibility to predict the rheology of the heavy crude oil in further work, benefiting from its dissipative force term.

The DPD simulations are primarily determined by the coarse-graining of model molecules and interaction parameters between DPD particles.21 One remarkable structural feature of heavy compounds in heavy crude oil is the portion of fused aromatic rings, which exhibits a characteristic rigidlike body. Many MD simulations9,11,13 have shown that the carbon atoms in fused aromatic rings approximately maintain their relative positions with respect to each other in one plane when running dynamics. Thus, it could be quite acceptable taking these fused aromatic rings as rigid bodies (we refer them as a rigid sheet in this paper) at the mesoscale. As a matter of fact, plenty of applications involving a rigid body (comprising of particles in a solid status) exist in DPD simulations. The colloids in nodal or square shape are treated as rigid bodies in the DPD module supported by Culgi,23 and the reflection boundary formed by a rigid body24,25 has been accepted as a useful technology in DPD simulations. Apparently, modifications on the DPD algorithm are necessary in order to hold the relative positions of DPD particles in rigid bodies, since all the molecules are represented by different types of freely moving particles in the standard DPD algorithm. In this paper, we integrated a rotational algorithm into the DPD program to evolve the motion of the rigid sheets. The procedure of determining the structure of model molecules was discussed, and some criteria were provided for judging the constructed models. The heavy crude oil and emulsion system are preliminarily explored to inspect the performance of the DPD platform. 2. Method Description 2.1. DPD Method. DPD is an approach proposed to describe the primarily hydrodynamic behavior of complex fluids.26-28 The evolution of the positions rij and impulses vij of all interacting particles is obtained by solving the Newtonian equation of motion dri =dt ¼ vi ð1Þ X D R S ðFC ð2Þ fi ¼ ij þ Fij þ Fij þ Fij Þ ¼ mi dvi =dt i6¼j

where mi is the mass of particle i, fi is the interactions exerted on particle i by other particles and can be simply written in terms of conservative, dissipative, random, and spring forces, repreD R S sented by FC ij , Fij , Fij , and Fij , respectively. ( aij ð1 - rij Þ^rij rij < rc FC ð3Þ ij ¼ 0 rij g rc ( - γωD ðrij Þð^rij 3 vij Þ^rij rij < rc D Fij ¼ ð4Þ 0 rij g rc ( σωR ðrij Þθij ^rij rij < rc R Fij ¼ ð5Þ 0 rij g rc

(12) Kuznicki, T.; Masliyah, J. H.; Bhattacharjee, S. Energy Fuels 2009, 23, 5027–5035. (13) Boek, E. S.; Yakovlev, D. S.; Headen, T. F. Energy Fuels 2009, 23, 1209–1219. (14) Headen, T. F.; Boek, E. S.; Skipper, N. T. Energy Fuels 2009, 23, 1220–1229. (15) Duda, Y.; Lira-Galeana, C. Fluid Phase Equilib. 2006, 241, 257– 267. (16) Qin, X.; Wang, P.; Sepehrnoori, K.; Pope, G. A. Ind. Eng. Chem. Res. 2000, 39, 2644–2654. (17) Nghiem, L. X.; Hassam, M. S.; Nutakki, R. Efficient Modeling of Asphaltene Precipitation. In SPE Annual Technical Conference and Exhibition, Houston, TX, October 3-6, 1993. (18) Aguilera-Mercado, B.; Herdes, C.; Murgich, J.; M€ uller, E. A. Energy Fuels 2006, 20, 327–338. (19) Groot, R. D.; Rabone, K. L. Biophys. J. 2001, 81, 725–736. (20) Qian, H.-J.; Lu, Z.-Y.; Chen, L.-J.; Li, Z.-S.; Sun, C.-C. Macromolecules 2005, 38, 1395–1401. (21) Yamamoto, S.; Maruyama, Y.; Hyodo, S.-a. J. Chem. Phys. 2002, 116, 5842–5849. (22) Rekvig, L.; Kranenburg, M.; Vreede, J.; Hafskjold, B.; Smit, B. Langmuir 2003, 19, 8195–8205.

FSij ¼ - cðrs - rij Þ^rij

ð6Þ

where rij is the distance between particle i and j, rc is the cutoff distance, rs is the length of the spring involved, aij is the (23) Culgi BV. Culgi, version 4.0.0; 2009. (24) Kim, J. M.; Phillips, R. J. Chem. Eng. Sci. 2004, 59, 4155–4168. (25) Revenga, M.; Zuniga, I.; Espanol, P. Comput. Phys. Commun. 1999, 122, 309–311. (26) Espanol, P.; Warren, P. Europhys. Lett. 1995, 30, 191–196. (27) Groot, R. D.; Warren, P. B. J. Chem. Phys. 1997, 107, 4423–4435. (28) Hoogerbrugge, P. J.; Koelman, J. M. V. A. Europhys. Lett. 1992, 19, 155–160.

4313

Energy Fuels 2010, 24, 4312–4326

: DOI:10.1021/ef1003446

Zhang et al.

conservative force parameter, γ is the friction factor, σ defines the fluctuation amplitude, c is the spring constant, and ωD(rij) and ωR(rij) are weight functions. To satisfy the fluctuation-dissipation theorem, the relationship between dissipative and random force parameters should be ensured as σ2 ¼ 2γkB T

ð7Þ

ωD ðrij Þ ¼ ½ωR ðrij Þ2

ð8Þ

and the relationship between weighting functions as ( rij 1rij < rc 2 D R ω ðrij Þ ¼ ½ω ðrij Þ ¼ rc 0 rij g rc

Table 1. Average Structure Parameters of Four Compounds of Most Crude Oil in Chinaa components

NH/NC

MW

RA

LP

saturate aromatic resin asphalteneb

1.9-2.0 1.5-1.7 1.4-1.5 1.1-1.3

600-900 700-1100 700-1200 800-1200

0 2-4 5-7 7-14

4-7 4-5 3-4

a NH/NC, number ratio of hydrogen and carbon atoms; MW, average relative molar weight; RA, number of aromatic rings; LP, average length of alkyl side-chain. b Separated by n-heptane.

ð9Þ

chains of asphaltene molecules were believed to have an average length of 5-6 carbons in asphaltene molecules and much more in other heavy compounds.39 Resins and the aromatics fraction, also containing aromatic rings, were like asphaltene of the continental type. In contrast, the smaller ring systems and longer alkyl side chains and the less number of heteroatoms and metals in functional groups exist in resins and aromatic molecules. Some structural characterization data of typical Chinese crude oil were listed in Table 1. Obtaining the coarse-grained model molecules is essential for the study based on the DPD method. Liew40 derived coarsegrained models from microscopic based on the pair-potential model. In this work, however, we construct model molecules at a rather low coarse-graining level by grouping three water molecules into one DPD particle, as some correlative works done at the mesoscale.19,21 The physical units could be derived from the promissory coarse-graining level. The mass scale can be defined as m ¼ Nm mwater ð11Þ

θij in eq 5 is a randomly fluctuating variable with Gaussian statistics and takes the following form in this work, with ξij as a random number drawn from a uniform distribution with zero mean: ξij θij ¼ pffiffiffiffiffiffi ð10Þ Δt The rotational algorithm for the motion of rigid body was adopted in our DPD program and has been shown in the Appendix in detail. 2.2. Model Molecules and Their Coarse-Graining Method. Generally, there are two strategies when we come to study the heavy crude oil through simulation works. Using average model molecule for each fraction of heavy crude oil is a common one as done in many works by MD9,29,30 since we still have no answer to how many compounds are in heavy crude oil or what are their structures exactly like. Another way is to construct a statistical representation of the oil sample which seemed a promising approach.13,31-34 At present, the performance of the mesoscale platform but not a specialized oil sample is what we are concerned with mostly. Thus, the previous strategy is taken and four models for heavy crude oil in terms of SARA fractionation were constructed. However, it is quite convenient to adopt the latter strategy in the future when the characterization data of the oil sample is available. Despite debates with respect to the structures of those heavy fractions like asphaltene and resin, it was common to take the idea that asphaltene molecules are shaped like a hand, the aromatic rings as the palm and the alkyl side chains as the fingers.7 Apart from the continental architectures (the hand model), another widely accepted model is archipelago architectures, which were thought to be composed of several smaller fused aromatic rings connected by bridge chains.5,35,36 Ruiz-Morales37 concluded the ring system roughly consisted of 4-10 fused aromatic rings to account for their color and other optical properties.7,37 What is more, there used to be 3-6 naphthenic rings closely coalescing with the fused aromatic core.38,39 Alkyl side

where Nm denotes the number of water molecules in a DPD particle and mwater the mass of the water molecule. According to Groot,41 the length scale Rc in angstroms and the time scale τ in picoseconds can be evaluated, respectively, as Rc ¼ 3:107ðFNm Þ1=3

ð12Þ

τ ¼ ð14:1 ( 0:1ÞNm 5=3

ð13Þ

and

Actually, the mass, length and time scales in physical units are m = 8.97  10-20 kg, Rc = 6.46  10-10 m, and τ = 8.8  10-11 s when we get the choice of Nm = 3 and F = 3. It is easy to know that the spatial and temporal range of the DPD simulation could be in the range of nanometers and micrometers. The fused aromatic rings, however, could not be exactly processed by grouping adjacent carbons as alkyl chains for its complex topology and proved to be the main challenge for the coarse-graining process. To our knowledge, little work has been done by DPD involving this kind of structure. Some particular concepts were adopted as follows to avoid splitting the fused aromatic rings directly. First, peri-condensation was taken as the basic pattern to construct the moiety of fused aromatic rings since it was more likely to yield a stable structure.42 To maintain its topologic features, the fused hexa-particle rings were chosen. All the efforts at preserving their topologic features could reproduce the offset π-stacked geometry phenomenon easily (see Figure 10), which had been observed in some MD simulations.9,30 Meanwhile, the rigidity of these aromatic rings was taken into account. Finally, the coarse-graining level should also be consistent with what we declared before. Figure 1 shows a coarse-graining model for the fused aromatic ring, which was constructed by the rigid sheet of the hexa-particle ring. This

(29) Vicente, L.; Soto, C.; Pacheco-Sanchez, H.; Hernandez-Trujillo, J.; Martinez-Magad an, J. M. Fluid Phase Equilib. 2006, 239, 100–106.  (30) Pacheco-S anchez, J. H.; Alvarez-Ramı´ rez, F.; Martı´ nezMagad an, J. M. Energy Fuels 2004, 18, 1676–1686. (31) Trauth, D. M.; Stark, S. M.; Petti, T. F.; Matthew Neurock, T.; Klein, M. T. Energy Fuels 1994, 8, 576–580. (32) Campbell, D. M.; Klein, M. T. Appl. Catal., A 1997, 160, 41–54. (33) Hou, Z.; Bennett, C. A.; Klein, M. T.; Virk, P. S. Energy Fuels 2010, 24, 58–67. (34) Sheremata, J. M.; Gray, M. R.; Dettman, H. D.; McCaffrey, W. C. Energy Fuels 2004, 18, 1377–1384. (35) Murgich, J. Mol. Simul. 2003, 29, 451–461. (36) Spiecker, P. M.; Gawrys, K. L.; Trail, C. B.; Kilpatrick, P. K. Colloids Surf., A 2003, 220, 9–27. (37) Ruiz-Morales, Y. Molecular Orbital Calculations and Optical Transitions of PAHs and Asphaltenes; Springer: New York, 2007; Vol. 4, pp 95-137. (38) Zijun, W.; Guohe, Q.; Wenjie, L.; Jialin, Q. Acta Petrolei Sinica 1999, 15 (6), 39–46. (39) Wenjie, L.; Guohe, Q.; Yuezhu, C. Acta Petrolei Sinica 1991, 7 (4), 1–11.

(40) Liew, C. C.; Mikami, M. Chem. Phys. Lett. 2003, 368, 346–351. (41) Groot, R. D. Langmuir 2000, 16, 7493–7502. (42) Andreatta, G.; Goncalves, C. C.; Buffin, G.; Bostrom, N.; Quintella, C. M.; Arteaga-Larios, F.; Perez, E.; Mullins, O. C.; Sheu, E. Y.; Hammami, A.; Marshall, A. G. Energy Fuels 2005, 19, 1282–1289.

4314

Energy Fuels 2010, 24, 4312–4326

: DOI:10.1021/ef1003446

Zhang et al. Table 2. Parameters of Conservative Force

A B C D

A

B

C

D

20 30 32 37

30 26 45 50

32 45 15 32

37 50 32 25

2.3. Parameters in DPD Simulation. Three parameters to be determined in DPD simulations, conservative force parameter aij between DPD particles i and j, noise parameter σ, and dissipative parameter γ, will clearly influence the simulation results.21 The conservative force parameter aij is the most important since it includes all the chemical details of the objective to be modeled. The noise and dissipative parameters related to the system temperature and fluid viscosity, respectively, also govern the behavior of some particular systems. Groot and Warren27 proposed the relationship between the conservative force parameter aij and the Flory-Huggins χparameter in 1997. The DPD conservative force parameter was investigated as a function of the particle size by Maiti,43 and as function of cohesive energy density by Travis.44 However, it is still a challenge to determine the conservative force parameters for the particles representing the extremely complex fused aromatic rings in heavy crude oil. Most of the conservative force parameters are determined in this work by the Blends45 method with slight adjustments to ensure their rationality of reproducing the dispersed structure of the colloidal system. Table 2 presents the conservative force parameters obtained by the Blends calculations at 298 K. Attention should be paid to the conservative force parameter between particles of type A because strong constraints are exerted on the particles of type A, making them act like rigid bodies. According to Groot’s calculations,27 the conservative force parameter has an inverse relation with the particle density represented as

Figure 1. Coarse-graining model for the fused aromatic rings.

aAA ¼ 75kB T=F

ð14Þ

The conservative force parameters between particles of type A should be aAA = 25 as F = 3. We roughly determined aAA = 20 since the particles of type A are placed in the rigid sheets and the particle density of type A must be larger than the even value of the system. In this work, we set the noise parameter σ = 3, the dissipative parameter γ = 4.5, and the spring constant C = 6 to construct rather soft side chains in model molecules. The dimensionless temperature T is initialized as 1.0 and will be monitored at each step in our DPD simulations. The default choice of the periodical box is 20  20  20 when no compressibility effect appears apparently. At this choice, the dimension of the box is roughly equivalent to 12.92 nm in physical units. In some simulation schemes, however, a much larger box is employed to make sure enough target molecules exist. At the beginning of a DPD simulation, most particles are distributed arbitrarily but evenly in the box. All the aromatic rings treated like rigid sheets are placed at the same orientation initially with the difference that the centroids of them were distritributed randomly in the periodical box. 2.4. Simulation Schemes. Unless otherwise specified, the parameters mentioned in section 2.3 were used as default values. There were four sorts of schemes in our study in terms of different compounds of heavy crude oil. Schemes I and II simply contained two fractions and were arranged to illustrate the structure of model molecules of asphaltene and resin from three aspects, rings system (I-1-I-8), distribution of side alkyl chains

Figure 2. Coarse-grained model molecules for DPD simulations (hexa-particle rings representing fused aromatic rings in molecules of heavy fractions of asphaltene, resin, and aromatic considered as rigid in this work. Particles in these hexa-particle rings are of the same type, and three different colors are used to mark them for clarity.)

coarse-grained method, specialized for treating fused aromatic rings, will be examined in section 3.1 by comparing with some relative experimental observations. The alkyl chains are represented by DPD particles and strings connecting them as previous work did.21 We set the average length of the side chains of the asphaltenes to 1 or 2 and 2-3 for the resin, 3-4 for the aromatic, and 4-8 for the saturate individually. Three issues about the model molecular structure, the size of the rings, the distribution of the connecting side chains, and the amount of heteroatoms will be discussed in section 3.1 too. Four types of DPD particles, A, B, C, and D, were employed to construct model molecules. A DPD particle of type A represented the moiety of aromatic rings, type B for the alkyl chain, type C for the functional group containing heteroatoms, and type D for water molecules. All types of DPD particles were visualized by different colors, as shown in Figure 2. We especially set three different colors for type A particles, red for the rigid moiety in asphaltenes, green for the resins, and blue for the aromatics, which could significantly benefit the detection of the specialized fraction from the simulation results.

(43) Maitia, A.; McGrother, S. J. Chem. Phys. 2004, 120 (3), 1595– 1602. (44) Travisa, K. P.; Bankhead, M.; Good, K.; Owens, S. L. J. Chem. Phys. 2007, 127, 1–12. (45) Accelrys. Materials Studio, Blends, version 4.1; San Diego, CA, 2008.

4315

Energy Fuels 2010, 24, 4312–4326

: DOI:10.1021/ef1003446

Zhang et al.

Figure 3. Continental model molecules of asphaltenes and their aggregate structures in heptane (side chains concealed).

Figure 3b-d. The second set consisted of Scheme III-19-III-21, differing from the first set with that the model molecules employed had no side chains. The last set included Scheme III-20, III-22, and III-23. They shared the same model molecule but had a different conservative force parameter aAB. We set aAB = 30 in Scheme III31, aAB = 45 in Scheme III-22, and aAB = 60 in Scheme III-23. Scheme IV (IV-302-IV-308) is to observe the emulsification and reveal the surfactant-like function of asphaltene.

(I-9-I-11 and II-14), and heteroatoms (I-11-I-13). Scheme III, consisting of all five fractions of heavy crude oil, was to study the structure of asphaltene nanoaggregates (III-15-III-23) and stability of heavy crude oil (III-24-III-301). To make it clear, we conducted three sets of modeling to investigate if structures of the model molecules and the conservative force parameter influenced the interlayer distance in a well-ordered structure. The first set included Scheme III-16-III-18 with model molecules shown in 4316

Energy Fuels 2010, 24, 4312–4326

: DOI:10.1021/ef1003446

Zhang et al. Table 3. Simulation Schemes

Scheme I: Asphaltene þ Heptane; masph:mhep = 1:19 no.

simulation system

I-1 I-2 I-3 I-4 I-5 I-6 I-7 I-8 I-9 I-10 I-11 I-12 I-13

asphaltene, Figure 3a; heptane, Figure 2 e asphaltene, Figure 3b; heptane, Figure 2e asphaltene, Figure 3c; heptane, Figure 2e asphaltene, Figure 3d; heptane, Figure 2e asphaltene, Figure 5a; heptane, Figure 2e asphaltene, Figure 5b; heptane, Figure 2e asphaltene, Figure 5c; heptane, Figure 2e asphaltene, Figure 5d; heptane, Figure 2e asphaltene, Figure 6a; heptane, Figure 2e asphaltene, Figure 6b; heptane, Figure 2e asphaltene, Figure 7a; heptane, Figure 2e asphaltene, Figure 7b; heptane, Figure 2e asphaltene, Figure 7c; heptane, Figure 2e Scheme II: Resin þ Toluene; mres:mtol = 1:19

no.

simulation system

II-14

resin, Figure 6c; toluene, Figure 2g Scheme III: Heavy Crude Oil; masph:mres:maro:msat:mlig = 1:3:3:2:17 no.

simulation system

III-15 III-16 III-17 III-18 III-19 III-20 III-21 III-22 III-23 III-24-III-301

no parameter is special asphaltene, Figure 3b; resin, Figure 2b; aromatic, Figure 2c; saturate, Figure 2d; light fraction, Figure 2e asphaltene, Figure 3c; resin, Figure 2b; aromatic, Figure 2c; saturate, Figure 2d; light fraction, Figure 2e asphaltene, Figure 3d; resin, Figure 2b; aromatic, Figure 2c; saturate, Figure 2d; light fraction, Figure 2e asphaltene, 3-ring without side chains; resin, Figure 2b; aromatic, Figure 2c; saturate, Figure 2d; light fraction, Figure 2e asphaltene, 4-ring without side chains; resin, Figure 2b; aromatic, Figure 2c; saturate, Figure 2d; light fraction, Figure 2e asphaltene, 5-ring without side chains; resin, Figure 2b; aromatic, Figure 2c; saturate, Figure 2d; light fraction, Figure 2e asphaltene, 4-ring without side chains; resin:, Figure 2b; aromatic, Figure 2c; saturate, Figure 2d; light fraction, Figure 2e; aAB = 45 asphaltene, 4-ring without side chains; resin, Figure 2b; aromatic, Figure 2c; saturate, Figure 2d; light fraction, Figure 2e; aAB = 60 masph:maro = 1:3 in III-24-III-253, masph:maro = 1:10 in III-254-III-301

Scheme IV: Heavy Crude Oil þ Water; Asphaltene, Figure 2a; Resin, Figure 2b; Aromatic, Figure 2c; Saturate, Figure 2d; Light Fraction, Figure 2e; Water, Figure 2f masph:mres:maro:msat:mlig:mwat

no. IV-302 IV-303 IV-304 IV-305 IV-306 IV-307 IV-308 IV-309

1:3:3:2:9:4 1:3:3:2:9:10 1:3:3:2:9:18 1:3:3:2:9:36 1:3:3:2:9:70 1:3:3:2:9:1 3:3:3:2:7:1 6:3:3:2:4:1

3.1.1. Aromatic Rings System. Both the continental and archipelago models were considered reasonable46 and hence compared in this work. In light of the finding by both MD simulations35 and experimental observations by X-ray47 that the large fused aromatic rings in asphaltene molecules tended to form ordered structure (stacking phenomena) in aggregates, the size of the rings system in the model molecules could be decided with the assistance of the appearance of this typical structure in the simulation. 3.1.1.1. Continental Models. Four sets of asphaltene model molecules are considered with the continental model in simulation Schemes I-1-I-14, containing two rings as shown in Figure 3a, three rings in Figure 3b, four rings in Figure 3c, and five rings in Figure 3d in peri-condensation individually. The corresponding snapshots of aggregates of asphaltene model molecules, after 300 000 step simulations, are presented in Figure 3e,f,g,h, respectively. In the case of the three-ring model molecule, most of the rigid sheets tend to form a big hoop in the periodical box and

The model molecules for each fraction adopted in simulation was just as Figure 2 showed unless expressly stated in Table 3. The mass ratio of binary system of asphaltene þ heptane was masph/mhep = 1:19 and that of the heavy crude oil system was masph/mres/maro/msat/mlig = 1:3:3:2:17. More details are listed in Table 3.

3. Results and Discussion 3.1. Average Model Molecules. To our best knowledge, there was no discussion regarding the average model molecules of heavy crude oil at the mesoscale. Here we stated the procedures and regulations to decide on the model molecules from three aspects of the structural feature, the aromatic rings system, side chains, and heteroatoms. Our attention was focused on the structural feature of asphaltene, since asphaltene was considered to be the typical fraction in heavy crude oil and its average model molecule was much more complex than those of other fractions. (46) Mullins, O. C.; Betancourt, S. S.; Cribbs, M. E.; Dubost, F. X.; Creek, J. L.; Andrews, A. B.; Venkataramanan, L. Energy Fuels 2007, 21 (5), 2785–2794.

(47) Sadeghi, M.-A.; Chilingarian, G. V.; Yen, T. F. Energy Sources 1986, 8, 99–123.

4317

Energy Fuels 2010, 24, 4312–4326

: DOI:10.1021/ef1003446

Zhang et al.

the others to form another small hoop closing to the big hoop. The aggregated asphaltene model molecules exhibit a face-to-face or offset π-stacked geometry ordered structure. The largest number of layers in the well-ordered structures is 3. Obviously, the clusters formed by four-ring model molecules are much more ordered than those by the three-ring model molecules. The average number of layers in the wellordered structures is about six in this case, which agrees well with the investigations by X-ray.47,48 When the number of rings in the model molecules comes to five, most of the rigid sheets turn to the same orientation approximately. There are as many as 30 molecules in line to form one huge ordered cluster. These results reveal that (1) A minimum number of hexa-particle rings exist to form a well-ordered structure. There should be three or more rings in an asphaltene model molecule if the hexa-particle ring is used to represent the aromatic core of the asphaltene molecule. (2) The highly ordered aggregating as shown in Figure 3h is not supported by experimental observations.49,50 Hence, we could estimate the maximum number of hexa-particle rings is no more than five in an average model molecule. (3) The relative molar mass of asphaltenes should be much less than 750 when the number of hexa-particle rings is smaller than three. Similarly, it should be much larger than 1200 when the number of hexaparticle rings is more than 5. Thus, the estimation of the minimum and maximum number of hexa-particle rings agrees with the relative molar mass of asphaltenes, suggesting the correctness of this coarse-graining method. Figure 3e-h indicates that the typical structures of asphaltene aggregates, tending to stack in a face-to-face or offset π-stacked geometry, are influenced greatly by the hexa-particle rings in the asphaltene model molecule. The census of the included angles among the rigid sheets of the asphaltene model molecules in a designed distance could be used to describe the characteristics of the ordered structure in aggregates. Two rigid sheets are considered to form a wellordered structure when the included angle is less than 0.1 radians (∼5°), while they form a T-shape and have poorlyordered characteristics when the angle approaches 1.57 radians (∼90°). Figure 4 presents the included angle distributions obtained from Schemes I-2-I-4, corresponding to Figure 3f-h. In the case of the three-ring model molecule, the distribution peaks that appeared at ∼0.8 and ∼1.1 rad indicate that some irregular structures or T-shaped geometry are formed in the asphaltene aggregates, as shown as the black circle marked in Figure 3f,g. For the five-ring model molecule, the included angles appear mainly in range of 0-0.2 radians, reflecting the well-ordered aggregates as shown in Figure 3h. The aggregate structures of the fourring model molecules reveal a better ordered aggregate than the three-ring model molecules do and worse than the fivering. Figure 4 presents also a hint that the peak value at the special included angle (∼0.1 radians) may be used to judge the ordered status of the fused aromatic rings in the asphaltene aggregates. 3.1.1.2. Archipelago Models. Nowadays, it remains dubious regarding the archipelago models. Mullins7 insisted that their large relative molar mass and colorless character

disagree with observations from their experiments. Nevertheless, Liao et al. supported this model.51 The archipelago model molecules of asphaltene were constructed deliberately as Figures 5a-d shows in simulation Schemes I-5-I-8 for comparison with the continental models. The morphologies of the asphaltene aggregates presented in Figure 5e-h reveal that (1) well-ordered structures are observed in all archipelago models. Apart from the intermolecular aromatic rings, obviously, intramolecular aromatic rings also tend to form ordered structures. (2) The number of ordered layers in aggregate structures is quite large, as many as 6 layers in Figure 5f, 12 in Figure 5g, and 12 in Figure 5h. The higher number of layers demonstrates that the archipelago models are more liable to form an ordered structure than the continental models. (3) Differing from the continental model, the ordered aggregate structure is no longer a change apparent with the number of rings in the archipelago model molecules. The reason could be considered as three rings in the archipelago model are large enough and the bridge chains also make their contribution to form a large stable ordered structure, as shown in Figure 5b. 3.1.2. Side Chains. The length and distribution of the side chains around the fused aromatic rings could play an important part in holding asphaltenes from aggregating by its steric hinderance effect as in polymer systems.46,52 The length issue had been discussed in our previous report.53 Here we focus on the amount of side chains. We set the average length of the side chains of asphaltenes to 1 or 2 and 2-3 for the resin, 3-4 for the aromatic, and 4-8 for the saturate individually, which was in line with the structural characterization data mentioned in section 2.2. It is also necessary to consider the way the side chains connected to the ring systems carefully. Symmetrical and unsymmetrical distribution, linear chains, fork chains, and bridge chains are all to be involved in model molecules. The discussions here focused on the amount of side chains connecting to the fused aromatic rings in the model molecules

(48) Shirokoff, J. W.; Siddiqui, M. N.; Ali, M. F. Energy Fuels 1997, 11, 561–565. (49) Trejo, F.; Ancheyta, J.; Rana, M. S. Energy Fuels 2009, 23, 429– 439. (50) Sharma, A.; Groenzin, H.; Tomita, A.; Mullins, O. C. Energy Fuels 2002, 16, 490–496.

(51) Liao, Z.; Zhao, J.; Creux, P.; Yang, C. Energy Fuels 2009, 23, 6272–6274. (52) Ran, Q.; Somasundaran, P.; Miao, C.; Liu, J.; Wu, S.; Shen, J. J. Colloid Interface Sci. 2009, 336, 624–633. (53) Zhang, S. F.; Sun, L. L. ; Xu, J. B.; Zhou, H.; Hao, W. Acta Phys.Chim. Sin. 2010, 26 (1), 57–65.

Figure 4. The included angle distributions of the three-ring, fourring, and five-ring model molecules of the asphaltenes.

4318

Energy Fuels 2010, 24, 4312–4326

: DOI:10.1021/ef1003446

Zhang et al.

Figure 7. Asphaltene model molecules containing different amounts of heteroatoms and the morphology of the aggregate.

which totally conflicts with the fact that asphaltenes are poorly miscible with heptane and should flocculate.54 Figure 6e reveals that all asphaltene molecules aggregate in a fluffy-sphere shaped cluster. The severe aggregate phenomena occures as a result of less side chains exist in the model molecule in Scheme I-10, as shown in Figure 6b. We may also understand the aggregating behavior at the stereohindrance effect point.55 Fused aromatic rings are eager to associate with each other for their poor miscibility with the alkyl environment, while the side chains could cover the surface of the fused aromatic core and hinder other cores to get too close.56,57 In the case of Scheme I-10, there are not enough side chains to balance the association tendency. Resin in heavy crude oils plays an important role in peptization of asphaltenes,56 and it also presents good miscibility with toluene.58 However, the hypothetical resin molecules we constructed with only two side chains, Figure 6c, show a stronge tendency to form aggregates in heptene as shown in Figure 6f. The reason could be same as to the one of Scheme I-9. The above discussions suggest that the proper amount of side chains in a model molecule of a specialized fraction (asphaltene, resin, and aromatic) could be estimated by observing its miscibility or dispersity among other fractions in heavy crude oils, which is attainable by experiments.

Figure 5. Archipelago model molecules of asphaltenes and their aggregate structures in heptane (side chains displayed).

Figure 6. Hypothetical unreasonalbe model molecules and their equilibrium morphology.

(54) Hammami, A.; Ratulowski, J. Precipitation and Deposition of Asphaltenes in Production Systems:A Flow Assurance Overview; Springer: New York, 2007; Vol. 23, pp 617-660. (55) Buenrostro-Gonzalez, E.; Groenzin, H.; Lira-Galeana, C.; Mullins, O. C. Energy Fuels 2001, 15 (4), 972–978. (56) Pereira, J. C.; L opez, I.; Salas, R.; Silva, F.; Fernandez, C.; Urbina, C.; L opez|, J. C. Energy Fuels 2007, 21, 1317–1321. (57) Porte, G.; Zhou, Ho.; Lazzeri, V. Langmuir 2003, 19, 40–47. (58) Ekholm, P.; Blomberg, E.; Claesson, P.; Auflem, I. H.; Sjoblom, J.; Kornfeldt, A. J. Colloid Interface Sci. 2002, 247 (2), 342–350.

of asphaltene. It was inevitable to draw the conclusion that asphaltene molecules containing more alkyl chains will have better dispersity, since alkyl chains are miscible with light fractions. The simulation result of Scheme I-9 presents, as shown in Figure 6d, a good dispersed system when we used the model molecule of asphaltene in Figure 6a. There are rare clusters containing more than three asphaltene molecules, 4319

Energy Fuels 2010, 24, 4312–4326

: DOI:10.1021/ef1003446

Zhang et al.

3.1.3. Heteroatom. The dynamic behavior of a molecule could be influenced greatly for its part containing heteroatoms which are supposed to be strongly related to the number of hydrogen bonds.59 There are many kinds of heteroatoms, and they can appear in the fused rings system or on the side chains. To understand the primary function of heteroatoms, we employ one sort of particle to represent all sorts of functional groups containing a heteroatom. No heteroatoms were arranged in model molecules of fractions in heavy crude oil except asphaltenes, in order to study the effect of heteroatoms on the aggregate structure. Three schemes were considered, without heteroatoms in asphaltene model molecules as shown in Figure 7a (Scheme I-11), containing some heteroatoms as in Figure 7b (Scheme I-12), and many as in Figure 7c (Scheme I-13). The corresponding morphologies of aggregates are shown in parts d, e, and f of Figure 7, respectively. The influence on the average aggregate number of the asphaltene model molecules and included angle distribution by the amount of heteroatoms can be observed by analyzing the trajectory files of the simulation. The average aggregate number in Schemes I-11, I-12, and I-13 is 3, 6, and 8, respectively, manifesting the average aggregate number will increase with the amount of heteroatoms in asphaltene model molecules. This tendency could also be observed visibly in Figure 7d-f. The reason could be due to the strong interactions between the heteroatomcontaining groups in the asphaltene molecules.59 Figure 8 illustrates that the rigid sheets in the aggregates will turn to a poorly-ordered structure as the amount of the heteroatom rises by comparing the included angle distribution peaks at ∼0.1 radians. The distribution peaks the appeared at ∼0.9 radians and ∼1.2 radians state that some asphaltenes molecules have roughly formed a T-shaped geometry. The phenomena could be seen vividly in Figure 7e,f. Thus, we consider that a relationship should exist between the formation of the T-shaped structure in asphaltene aggregates and the appearance of heteroatoms. A similar investigation had been made by J. Correa-Basurto’s work60 where they found a dihedral angle of 136.968° by ab initio calculations. 3.2. Heavy Crude Oil System. In our effort to certify the capability of the mesoscale platform in exploring the heavy crude oil system, several aspects related to this mixture were studied and comparison with experimental data had been conducted. 3.2.1. Nanoaggregate Structure. As Figure 9a (Scheme III15) shows, the fused aromatic rings (in red) of asphaltene molecules attract each other and form the cores of the nanoaggregates. The side chains (in black) outside the asphaltene molecules act as the repulsive part by asteric hinderance effect and repel additional asphaltenes molecules.46 For this reason, resin fractions, which seem to surround the aggregates, have been repulsed and actually do not participate in forming the nanoaggregates.46,61 On a larger scale, there are about three nanoaggregates (in the middle of Figure 9a) combining to form a cluster of nanoaggregate, which is a bit smaller than the one in Mullins’s modified Yen model.61 The deficiency of asphaltene molecules in this simulation system

could partly explain the difference. In this small cluster of nanoaggregates, we can find that, apart from being the external part, the side chains also act internally. Figure 9b, c shows an asphaltene nanoaggregate from Figure 9a, which contains eight asphaltene molecules and has as many as five molecules in the ordered structure. The main structural feature is fairly similar to the schematic structure proposed in Figure 9d and being supported by Mullins’s work.62 The nanoaggregate structure proposed in Figure 9d is not spacefilling but, instead, rather fluffy. Locally, several fused aromatic rings tend to collaborate orderly but they are in a rather disordered arrangement in total in Figure 9d. Still there are as many as three aromatic rings in the well-ordered status. As a matter of fact, there used to be three to seven layers in this kind of structure during our simulations, agreeing well with some experimental results.47,48 Figure 10 exhibits the distribution of asphaltene molecules and three typical aggregate structures are observed: (a) π-π repulsion results in a face-to-face geometry, (b) π-σ attraction results in a T-shaped geometry, and (c) σ-σ attraction results in an offset π-stacked geometry as in the rules summarized by Hunter and Saunders.63 The reproduction of these three typical geometries certifies the capacity of the mesoscale platform at integrating the interaction mechanism among the asphaltene molecules. 3.2.2. Interlayer Distance. There were three sets of modeling (III-15-III-23) conducted to investigate if the structures of the model molecules and the conservative force parameter influenced the interlayer distance in the well-ordered structure. The arrangement of modeling has been stated in section 2.4, and the results of these simulations are plotted in Figure S1 in the Supporting Information. In Scheme III-16-III-18, the maximum and minimum values of the interlayer distance appear in the cases of fourring and five-ring model molecules of asphaltene, respectively, as shown in Figure S1 in the Supporting Information. Figure S1 in the Supporting Information demonstrates that the interlayer distance has a somewhat weak relationship with the size of the fused aromatic rings. This may be due to the fact that the

(59) Mansoori, G. A.; Dynora, V.; Shariaty-Niassar, M. J. Pet. Sci. Eng. 2007, 58, 375–390. (60) Basurto, J. C.; Aburto, J.; Ferrara, J. T.; Torres, E. Mol. Simul. 2007, 33 (8), 649–654. (61) Mullins, O. C. Energy Fuels 2010, 24, 2179–2207.

(62) Andreatta, G.; Bostrom, N.; Mullins, O. C. Langmuir 2005, 21, 2728–2736. (63) Hunter, C. A.; Sanders, J. K. M. J. Am. Chem. Soc. 1990, 112 (14), 5525–5534.

Figure 8. The included angle distributions of asphaltene model molecules containing different amounts of heteroatoms.

4320

Energy Fuels 2010, 24, 4312–4326

: DOI:10.1021/ef1003446

Zhang et al.

Figure 9. Nanoaggregate structure in Scheme III-15 of heavy crude oil.75

hibit little smaller interlayer distances when compared with those in Scheme III-16-III-18. This difference should mainly be stemmed from the fact that long alkyl side chains connected with aromatic cores could cause the structures to be disordered and amorphous.49 Trejo et al.49 found that the poor-ordered structures determined by TEM have an interlayer distance of 0.355 nm corresponding to amorphous asphaltenes; however, the shortage of long chains may create a well-ordered structure which resembles graphite.49 The perfectly ordered layers used to have an interlayer separation of ∼0.335 nm. It is easy to find from Figure S1 in the Supporting Information that the largest interlayer distance occurs in the simulation Scheme III-22 with aAB = 45, neither in the Scheme III-20 with aAB = 30 nor the Scheme III-23 with aAB = 60. Apart from this phenomenon, the offset geometry is quite seldom observed during the simulations of Schemes III-20, III-22, and III-23. These results could reveal that the interlayer distance is independent of the conservative force parameters. The calculated value of the interlayer distance is distributed in the range of 0.375-0.405 nm and slightly larger than the experimental data. Two factors should be responsible for the difference. On the one hand, our simulation systems

Figure 10. Morphology of asphaltenes in heavy crude oil (other fractions have been removed for clarity).

interlayer distance rarely changes much even if the heavy crude oil samples could be from various parts of the world.47,48 The aggregate structures formed by the asphaltene model molecules (without side chains) in Scheme III-19-III-21 ex4321

Energy Fuels 2010, 24, 4312–4326

: DOI:10.1021/ef1003446

Zhang et al.

Figure 11. Shell’s SARA cross plot stability screen.54,64

Figure 12. A cross-plot of heavy crude oil from Simulations.

study asphaltenes in solvents like toluene and heptane. The rigid sheets in simulation tend to adopt offset π-stacked rather than face-to-face geometry in most situations and also form some significant irregular geometry. In other words, our simulation systems consist mostly of poorly ordered structures. By contrast, as Fernando Trejo49 concluded with the assistance of SEM/TEM, the rigid sheets of the asphaltenes in the vicinity tend to form well-ordered structures during their analysis of precipitated asphaltenes. The other one lies in the method we used to determine the interlayer distance. As described in Figure S2 in the Supporting Information, (d1 þ d2)/2 is adopted in this work to represent the interlayer distance and it should be different form what they adopted in experiments. 3.2.3. Concentration and Stability of Heavy Crude Oil. Asphaltene stability of crude oil is assessed by various techniques such as SARA screens and titration tests using n-heptane. Shell’s approach64 underlined two ratios, saturate/aromatic vs asphaltene/resin based on the SARA screen. Stankiewicz et al.64 and Hammami et al.54,64 presented their cross-plot of these ratios according to the exploration work on more than 230 sample oils worldwide, as shown in Figure 11. Most of the sample oils are, in Figure 11, in stable status when the ratio of asphaltene/resin is less than 0.2. As the ratio rises, however, rare sample oils remain stable unless an extremely low ratio of saturate/ aromatic is satisfied. One thing we should bear in mind there is no specification for masph/maro or masph/msat in Shell’s scheme. Two sets of simulation systems are involved in this work to check the performance of our DPD simulation on the stability prediction for crude oil. One set of simulation systems, Scheme III-24-III-253, has the mass fraction of asphaltene in whole oil maintained as 1:26 and the mass ratio of asphaltene/aromatic masph/maro = 1:3. Another set of simulation systems, Scheme III-254-III-301 has the mass ratio of asphaltene/aromatic masph/maro = 1:10. Here we choose light or medium crude oil for our interest in how the two ratios influence the stability of crude oil as well as the clarity of visualization.

The cluster size is taken as the parameter in this work for judging the status of a crude oil. The crude oil will be considered to be marginal when the maximum cluster size approaches eight asphaltene molecules and be stable or unstable when less or more than eight asphaltene molecules. The stable area of the first set of simulation systems (Scheme III-24-III-253) can be divided into three zones, marked as A, B, and C, respectively, in Figure 12. Zone A represents such a status that the crude oil will remain stable among the whole range of ratios of saturate/aromatic, demonstrating the stability of crude oil and the independence of the amount of saturate or aromatic because of enough resin in the oil. This phenomenon is pretty similar to the area marked “resin controlled” in Figure 11. In fact, resin is a key factor in determining stabilization of asphaltene clusters against aggregation.65 Resin molecules are able to produce a protective “wrapping”65 around the asphaltene clusters and prevent their aggregation by a steric repulsion mechanism. The peptization happens when asphaltene clusters are dispersed by plenty of resin, as frequently reported in the literature.66,67 In zone B, however, there is not enough resin to disperse asphaltene clusters and thus the ratio of saturate/ aromatic should be limited to low values to maintain the oil stability. Zone C is considered as a part of the stable phase since the simulation systems tend to be stable when the ratio of saturate/aromatic is set to an exceedingly low value, though there is no stable sample oils appearing in Shell’s cross-plot. This could be understood when we realize that asphaltene molecules could be in a good dispersed status in aromatic solvents, like toluene.68 Two curves separating stable and unstable oils have been traced in Figure 12, according to two set of simulations. The first one, as mentioned above, was with the ratio of asphaltene/aromatic = 1:3. The other set (Scheme III-254-III-301) was with the ratio of asphaltene/aromatic = 1:10. As can be seen from Figure 12, the second curve moved toward the upper right relative to the first one, demonstrating the expanded stable area. This could be explained by the decreased amount of the heaviest fraction and increased number of aromatic

(64) Stankiewicz, A. B.; Flannery, M. D.; Fuex, N. A.; Broze, G.; Couch, J. L.; Dubey, S. T.; Iyer, S. D.; Ratulowski, J.; Westerich, J. T. Prediction of asphaltene deposition risk in E&P operations. In Proceeding of 3rd International Symposium on Mechanisms and Mitigation of Fouling in Petroleum and Natural Gas Production, AIChE 2002 Spring National Meeting, New Orleans, LA, 2002; Vol. 47C, pp 410-416.

(65) Pereira, J. C.; L opez, I.; Salas, R.; Silva, F.; Fernandez, C.; Urbina, C.; L opez|, J. C. Energy Fuels 2007, 21, 1317–1321. (66) Carnahan, N. F.; Salager, J. L.; Ant on, R.; Davila, A. Energy Fuels 1999, 13, 309–314. (67) Porte, G.; Zhou, H.; Lazzeri, V. Langmuir 2003, 19, 40–47. (68) Anglea, C. W.; Yicheng, L.; Hassan, H.; Lueb, L. Fuel 2006, 85, 492–506.

4322

Energy Fuels 2010, 24, 4312–4326

: DOI:10.1021/ef1003446

Zhang et al.

Figure 13. The equilibrated morphologies of emulsion systems (Scheme IV-302-IV-306).

and saturate molecules in the system. Additionally, the comparison among these curves manifests that the stability

of oil cannot be decided directly by two ratios as in Figure 12. Detail of the composition of the oil is required. 4323

Energy Fuels 2010, 24, 4312–4326

: DOI:10.1021/ef1003446

Zhang et al.

Figure 15. The equilibrated morphology in Scheme IV-308.

Figure 14. The equilibrated morphology in Scheme IV-307.

3.3. Emulsion. More complex fluid system as water-oil emulsion was explored by mixing water model molecules and heavy crude oil system. In Schemes IV-302-IV-306, which were designed to show the change of morphology of emulsions as the ratio of oil to water altered, we maintain the mass ratio of oil as 1:3:3:2:9 and adjust the amount of water. From the snapshots (Figure 13), it was found that the morphology of emulsions experience stages of w/o (Figure 13a,f), lamellae (Figure 13b,g), bulk (Figure 13c,h), rod (Figure 13d,i), and o/ w (Figure 13e,j), which exactly resemble those of surfactants in the water-oil system.69 Basically, the heavy fractions stay in the oil phase as all five cases display. The asphaltene fraction gets the most close to the aqueous phase and forms the interfacial layer, followed by the resin fraction. Virtually, few aromatic molecules distribute among the interfacial layer. All the morphology of emulsions confirm the ability of asphaltene to act as an emulsifier and also exhibit the resin’s function as a much weaker emulsification.70 There is no stacking or ordered structures appearing in Figure 13a, and the interfacial film consists of only one layer of asphaltene molecules. The reason could relate to the lack of asphaltenes which can also cause part of the surface of the water sphere to be uncovered. The interfacial film turns to a plane in Figure 13b and is formed by one layer of asphaltene molecules too. The even lower amount of asphaltene could explain that. The equilibrated morphology of lamellae shows that most asphaltene molecules are orientated parallel to the interface. It also can be seen in this figure that some resin molecules appear near the interfacial film, reflecting its feature of a weak surfactant agent. As we increase the amount of asphaltene molecules and decrease those of light fractions correspondingly and keep the ratio of other compounds in emulsion at the same time in Scheme IV-307-IV-309, the number of water molecules dispersed in the oil phase raise slightly as seen in Figures 14a, 15a, and 16a. This may relate to the solubilization of asphaltenes. Additionally, the surplus asphaltene molecules could split the original one water sphere into two small ones

Figure 16. The equilibrated morphology in Scheme IV-309.

as Figures 15 and 16 display. The phenomena reflect that asphaltenes tend to disperse around the interfacial film rather than in the oil phase. Thus, they increase the surface of the interfacial film by forming more interfaces. The stacking phenomena appear as the ratio of asphaltene comes to a very high level as in Scheme IV-309. From the enlarged chart of one water sphere in Figure 16c, we can see many ordered structures containing two layers that appear near the interface, indicating a much thicker interfacial film is formed. It is common to take the idea that ahigher concentration of asphaltene in the oil phase will increase the amount of asphaltene adsorbed at the oil-water interface. Meanwhile, well ordered structures developed among those asphaltene aggregates in the interfacial film means enhanced interfacial film strength and stability of the emulsion.71 As to the dehydration and desalting operation, it could be feasible to develop additions which could prohibit surplus asphaltene adsorption at the oil-water interface efficiently and facilitate further coalescence among drops. 4. Conclusions A mesoscale platform based on the DPD method was constructed and applied to model the aggregating behavior of asphaltenes in heavy crude oils for the first time. To reflect the distinct characteristic of the fused aromatic rings in

(69) Alexandridis, P.; Olsson, U.; Lindman, B. Langmuir 1998, 14, 2627–2638. (70) Ali, M. F.; Alqam, M. H. Fuel 2000, 79, 1309–1316.

(71) Eley, D. D.; Hey, M. J.; Lee, M. A. Colloids Surf. 1987, 24 (2-3), 173–182.

4324

Energy Fuels 2010, 24, 4312–4326

: DOI:10.1021/ef1003446

Zhang et al.

The time evolution of a quaternion Q obeys the following relation, : 1 ðA4Þ Q ¼ Bω b ðtÞ 2

hydrocarbons, a rigid sheet was introduced into model molecules. Model molecules for each fraction in heavy crude oil were determined according to structural data in the literature and our simulation experience. Three types of associated geometry were observed as some MD work did. The interlayer distance of ordered structures obtained from modeling was a bit larger than 0.36 nm reported by pure asphaltene observation. This value remained no matter what kind of structure of model molecules and repulsive parameter were chosen. These results agreed with experimental data well, illustrating the correctness of the model molecules to some extent. In addition, the coarse-graining process saved time greatly, and the results from the DPD calculation were similar to those from the MD work, making it quite clear that the developed platform was a kind of quick MD. The advantage of this platform also was embodied by the flexibility of extending statistics as we did in analyzing the included angle in ordered structures and computing the size of clusters, benefiting from our original program. At last, we took the platform to study the heavy crude oil system and emulsion system preliminarily. The relationship between concentration and the stability of crude oil is consistent with Shell’s work showing the potential capacity of this platform to investigate a large and complex fluid system like heavy crude oil. As a matter of fact, the emulsified mechanism was studied recently and the structural characteristic of the emulsifier and demulsifier was identified. Most of all, the rheology of heavy crude oils and emulsion will be investigated in further work, which could shed some useful light on reducing the viscosity of these complex fluids.

where B and ωb(t) are defined as follows: 0 1 q0 - q1 - q2 - q3 B C 1B q1 q0 - q3 q2 C B C B ¼ B C 2 @ q2 q3 q0 - q1 A q3 - q2 q1 q0 1 0 B ωb ðtÞ C x C ω b ðtÞ ¼ B @ ωby ðtÞ A ωbz ðtÞ

0

0 1 - q2 2 þ q1 2 - q3 2 þ q0 2 2ðq3 q0 - q1 q2 Þ 2ðq1 q3 þ q2 q0 Þ B C B - 2ðq q - q q Þ q2 2 - q1 2 - q3 2 þ q0 2 2ðq1 q0 - q2 q3 Þ C @ A 2 1 3 0 - 2ðq2 q3 þ q1 q0 Þ - q22 þ q21 þ q23 þ q20 2ðq1 q3 - q2 q0 Þ

ðA7Þ There were several rotational algorithms of the existing integrating quaternion, such as Fincham’s implicit leapfrog algorithm72 and Ruocco and Sampoli’s integral leapfrog algorithm.73 We stated the main steps of an explicit leapfrog algorithm, which saved time and proved to be accurate enough. Another version of the midstep leapfrog was also implemented. As we all know, the standard leapfrog algorithm could be expressed like this74 : vðt þ Δt=2Þ ¼ vðt - Δt=2Þ þ f ðtÞΔt=m ðA8Þ rðt þ ΔtÞ ¼ rðtÞ þ vðt þ Δt=2ÞΔt

where LS(t) is the angular momentum and TS represents the sum of torques. Since the rotational velocity at each step is required to compute the dissipative force (which is quite different from the calculation executed in conservative force fields), we rewrite eq A10 as ðA12Þ Ls ðt þ Δt=2Þ ¼ Ls ðtÞ þ Ts Δt=2

Appendix: Quaternion Method The unit quaternion, depicting the orientation, is defined as follows ðA1Þ Q ¼ q0 þ q1 i þ q2 j þ q3 k

Ls ðt þ ΔtÞ ¼ Ls ðtÞ þ Ts Δt

ðA13Þ

Q(t), ωb(t), and LS(t) are known at the initial stage of iteration, _ thus we could use eq A4 to gain the Q(t). Then with substitution ~ þ of this into eq A14, Q(t þ Δt/2) and the orientation matrix A(t Δt/2) correspondingly could be estimated. : Qðt þ Δt=2Þ ¼ QðtÞ þ QðtÞΔt=2 ðA14Þ

ðA2Þ

and qi 2 ¼ 1

ðA9Þ

where t represents the current time and Δt is the time step. We can write a leapfrog algorithm for rotational motion in similar form as follows: ðA10Þ Ls ðt þ Δt=2Þ ¼ Ls ðt - Δt=2Þ þ Ts Δt : Qðt þ ΔtÞ ¼ QðtÞ þ Qðt þ Δt=2ÞΔt ðA11Þ

Supporting Information Available: Relationships between the interlayer distance and the model structure and the repulsive parameter and a sketch to illustrate the way of computing the interlayer distance in the simulation. This material is available free of charge via the Internet at http://pubs.acs.org.

3 X

ðA6Þ

A rotation matrix could be expressed by a set of quaternion, ~ ¼ A

Acknowledgment. We acknowledge the financial support from the National Natural Science Foundation of China (20776141), Foundation for Innovative Research Groups of the National Natural Science Foundation of China (20821092) and State Key Laboratory of Multiphase Complex System. We thank Prof. Johannes GEM Fraaije of Leiden university for fruitful discussions and RuiXing Wang from Culgi for valuable suggestions. We are grateful to Prof. Han Zhou from Research Institute of Petroleum Processing of SINOPEC for the use of Material Studio.

8 2 i ¼ j2 ¼ k2 ¼ ijk ¼ - 1 > > > < ij ¼ - ji ¼ k > jk ¼ - kj ¼ i > > : ki ¼ - ik ¼ j

ðA5Þ

(72) Fincham, D. Mol. Simul. 1992, 8, 165–171. (73) Ruocco, G.; Sampoli, M. Mol. Phys. 1994, 82, 875. (74) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, U.K., 1987. (75) Andreatta, G.; Goncalves, C. C.; Buffin, G.; Bostrom, N.; Quintella, C. M.; Larios, F. A.; Perez, E.; Mullins, O. C.; Sheu, E. Y.; Hammami, A.; Marshall, A. G. Energy Fuels 2005, 19, 1282–1289.

ðA3Þ

i¼0

Q means a rotation around the vector q1i þ q2j þ q3k with a rotation angle of 2arccos q0 physically. 4325

Energy Fuels 2010, 24, 4312–4326

: DOI:10.1021/ef1003446

Zhang et al.

Combining with the expression of L (t þ Δt/2) by eq A12, the angular velocity at midstep in the body coordinate could be obtained. ~ þ Δt=2Þ Ls ðt þ Δt=2Þ ðA15Þ ω b ðt þ Δt=2Þ ¼ ðIb Þ - 1 3 Aðt 3 S

Displacement of particles in the rigid body can be updated as ~ þ ΔtÞ - 1 Rb Rsi ðt þ ΔtÞ ¼ ½Aðt 3 i

So far as now, the procedure listed above is referred to as the explicit leapfrog algorithm commonly. If one requires that : Qðt þ Δt=2Þ ¼ QðtÞ þ Qðt þ Δt=2ÞΔt=2 ðA19Þ

_ þ Δt/2) could be calculated according to eq A16. Then Q(t : 1 ðA16Þ Qðt þ Δt=2Þ ¼ Bω b ðt þ Δt=2Þ 2 According to eq A11, the quaternion Q(t þ Δt) and orientation ~ þ Δt) could correspondingly be updated. Then the matrix A(t rotational velocity in space coordinate could be updated as ~ þ ΔtÞ - 1 ω s ðt þ ΔtÞ ¼ ½Aðt -1 ~ b s 3 ½I ðt þ ΔtÞ 3 Aðt þ ΔtÞ 3 L ðt þ ΔtÞ

ðA18Þ

Equation A19 can be solved by taking eq A14 as the beginning value, repeating eqs A15, A16, and A19 until the value of Q(t þ Δt/2) converges. The procedure with the special iteration is the so-called midstep implicit leapfrog algorithm.

ðA17Þ

4326