Molecular Dynamics Study of Model Molecules Resembling

Energy & Fuels 2008, 22, 2379–2389

2379

Molecular Dynamics Study of Model Molecules Resembling Asphaltene-Like Structures in Aqueous Organic Solvent Systems Tetyana Kuznicki,†,‡ Jacob H. Masliyah,‡ and Subir Bhattacharjee*,† Department of Mechanical Engineering, UniVersity of Alberta, Edmonton, Alberta, T6G 2G8, Canada, and Department of Chemical and Materials Engineering, UniVersity of Alberta, Edmonton, Alberta, T6G 2G6, Canada ReceiVed January 24, 2008. ReVised Manuscript ReceiVed April 8, 2008

Molecular dynamics (MD) simulations were used to study the aggregation of model molecules representing some structural features of asphaltenes in different solvents, namely, water, toluene, and heptane. Three different model molecules with similar molecular weights were developed as representatives of continental- and archipelago-type structures of asphaltenes. Two types of continental models were employed, one containing no charged entities and another containing one anionic terminal carboxylic group on an aliphatic chain. Analysis of the aggregate structures in single pure solvents indicate stacking of the polyaromatic rings as the dominant structural feature of the aggregates. Simulations were then conducted to study aggregation and partitioning of these model molecules in binary mixtures of water and toluene. A key observation from these simulations was that uncharged molecules did not congregate at the toluene-water interface, whereas charged terminal groups had a distinct affinity for the toluene-water interface. In the presence of charged entities, the resulting aggregate had a complex three-dimensional structure containing stacked polyaromatic rings tethered to the interface by the anionic terminal groups.

1. Introduction Asphaltenes, the heaviest fractions of petroleum, are described as a solubility class of molecules that precipitate in paraffinic solvents but dissolve in excess toluene.1 Generally, asphaltenes are thought of as a collection of diverse molecules consisting predominantly of polynuclear aromatic units, varying proportions of aliphatic chains, small amounts of heteroatoms (such as S, N, and O), and traces of elements, such as V and Ni. The hydrogen/carbon ratios between 1.0 and 1.2 are typical with the heteroatom content of a few percent.2 As a consequence of their extraordinary variety and complexity, asphaltenes have been considered as “enigmatic”,3 with several aspects of their physics and chemistry being highly debated. Asphaltenes cause serious operational problems in oil recovery and transport, ranging from plugging of reservoirs to fouling of production pipelines. Most of these problems can be traced to aggregation and precipitation of asphaltenes. Self-aggregation of asphaltenes has been described as the first step in the formation of small particles.4 It has been hypothesized that asphaltene molecular moieties associate to form micellar structures through hydrogen bonds.5 More recent studies state that charge transfer between molecules may be a mechanism * To whom correspondence should be addressed. Telephone: (780) 4926712. Fax: (780) 492-2200. E-mail: [email protected]. † Department of Mechanical Engineering. ‡ Department of Chemical and Materials Engineering. (1) Speight, J. G. The Chemistry and Technology of Petroleum (electronic resource), 4th ed.; CRC Press/Taylor and Francis: Boca Raton, FL, 2007. (2) Vicente, L.; Soto, C.; Pacheco-Sanchez, H.; Hernandez-Trujillo, J.; Martinez-Magadan, J. M. Fluid Phase Equilib. 2006, 239, 100–106. (3) Hortal, A. R.; Hurtado, P.; Martinez-Haya, B.; Mullins, O. C. Energy Fuels 2007, 21, 2863–2868. (4) Andersen, S. I.; Birdi, K. S. J. Colloid Interface Sci. 1991, 142, 497–502. (5) Ho, B.; Briggs, D. E. Colloids Surf. 1982, 4, 285–303.

leading to self-aggregation.6 Self-aggregation of asphaltenes in organic solvents has been observed by calorimetry,4 surface tension5 and interfacial tension measurements,7,8 vapor pressure osmometry,8 Langmuir trough experiments,9–13 small-angle neutron scattering,6,14–16 and transmission electron microscopy.17 The structure of asphaltene nanoaggregates in bulk pure solvents comprises a core region formed by stacking of a few asphaltene polyaromatic rings, surrounded by several pendant aliphatic chains.18,19 Such aggregates, owing to the presence of (6) Asphaltenes: Fundamentals and Applications; Sheu, E. Y., Mullins, O. C., Eds.; Plenum Press: New York, 1995. (7) Mohamed, R. S.; Ramos, A. C. S.; Loh, W. Energy Fuels 1999, 13, 323–327. (8) Yarranton, H. W.; Alboudwarej, H.; Jakher, R. Ind. Eng. Chem. Res. 2000, 39, 2916–2924. (9) Zhang, L. Y.; Lawrence, S.; Xu, Z. H.; Masliyah, J. H. J. Colloid Interface Sci. 2003, 264, 128–140. (10) Horvath-Szabo, G.; Masliyah, J. H.; Elliott, J. A. W.; Yarranton, H. W.; Czarnecki, J. J. Colloid Interface Sci. 2005, 283, 5–17. (11) Zhang, L. Y.; Lopetinsky, R.; Xu, Z. H.; Masliyah, J. H. Energy Fuels 2005, 19, 1330–1336. (12) Cadena-Nava, R. D.; Cosultchi, A.; Ruiz-Garcia, J. Energy Fuels 2007, 21, 2129–2137. (13) Zhang, L. Y.; Breen, P.; Xu, Z. H.; Masliyah, J. H. Energy Fuels 2007, 21, 274–285. (14) Gawrys, K. L.; Blankenship, G. A.; Kilpatrick, P. K. Langmuir 2006, 22, 4487–4497. (15) Jestin, J.; Simon, S.; Zupancic, L.; Barre, L. Langmuir 2007, 23, 10471–10478. (16) Verruto, V. J.; Kilpatrick, P. K. Energy Fuels 2007, 21, 1217– 1225. (17) Sharma, A.; Mullins, O. C. High-resolution transmission electron microscopy of asphaltenes. In Asphaltenes, HeaVy Oils and Petroleomics; Mullins, O. C., Sheu, E. Y., Hammami, A., Marshall, A. G., Eds.; Springer: New York, 2007; Chapter 10. (18) McLean, J. D.; Kilpatrick, P. K. J. Colloid Interface Sci. 1997, 196, 23–34. (19) Mullins, O. C.; Betancourt, S. S.; Cribbs, M. E.; Dubost, F. X.; Creek, J. L.; Andrews, A. B.; Venkataramanan, L. Energy Fuels 2007, 21, 2785–2794.

10.1021/ef800057n CCC: $40.75  2008 American Chemical Society Published on Web 06/19/2008

2380 Energy & Fuels, Vol. 22, No. 4, 2008

Kuznicki et al.

low-energy aliphatic chains at the periphery, should not show significant interfacial activity. However, asphaltenes are regarded as the key stabilizers of oil-water emulsions, acting as “natural surfactants” in petroleum.11,12,15,18 In this context, asphaltene aggregation at oil-water interfaces has received considerable attention. A commonly held notion is that asphaltene polyaromatic rings align themselves parallel to the interface,18 forming a two-dimensional monolayer structure,11,13,20 and over time, realignment of this layer provides a steric barrier toward rupture of the interface. More recently, three-dimensional (3D) structures forming at the interface are being postulated as reasons of the interfacial stability.12 A theoretical understanding of ashphaltene aggregation is quite difficult because these molecules cannot be treated as simple homogeneous entities that are amenable to common meso-scale aggregation theories based on models of colloidal interactions. A more detailed consideration of the molecular architecture of asphaltenes necessitates the use of molecular simulation techniques.21,22 In this context, the present study aims at testing the applicability of large-scale classical molecular dynamics (MD) simulations based on empirical force fields to describe aggregation and interfacial behaviors of asphaltenelike molecules. These simulations are still too computationally intensive to barely allow for the calculation of trajectories over small fractions of a microsecond. Hence, they only provide a glimpse of very short-range (high-frequency) phenomena in complex fluids. Application of such techniques to predict the behavior of asphaltene-like systems, which evolve through extremely slow dynamics, may be of limited value. However, it is of interest to explore how these simulations can predict some of the high-frequency processes, such as self-aggregation (micellization), partitioning at oil-water interfaces, or stacking of polyaromatic rings. No single molecular model is currently thought to represent the “average” asphaltene molecule. The number and manner in which aromatic rings are placed within a molecule can make a significant difference in properties of the resulting model.23 Despite debates with respect to a generalized shape for asphaltenes,24 it is common to structurally classify asphaltenes into “continental”- (pericondensed) and “archipelago”-type generic architectures. The continental model resembles a relatively flat disk-like molecule, with a center that has a dominant aromatic core and a periphery that is characterized by aliphatic chains3,19,22,23,25 The “archipelago” model consists of aliphatic chains interconnecting groups of small aromatic regions.21,22,29 The actual structure of these molecules is not well-defined, with even uncertainty regarding the number of aromatic rings present in the planar aromatic cores.26,27 Alkyl side chains are believed to have an average length of 5-6 carbons. Continental-type models have been used in molecular modeling studies to demonstrate the formation of aggregates by stacking of their aromatic regions.28 The present study is based on hypothetical molecular models that comprise some of the commonly postulated structural (20) Zhang, L. Y.; Xu, Z. H.; Masliyah, J. H. Ind. Eng. Chem. Res. 2005, 44, 1160–1174. (21) Murgich, J.; Abanero, J. A.; Strausz, O. P. Energy Fuels 1999, 13, 278–286. (22) Zhang, L. Q.; Greenfield, M. L. J. Chem. Phys. 2007, 127, 194502. (23) Murgich, J. Mol. Simul. 2003, 29, 451–461. (24) Strausz, O. P.; Safarik, I.; Lown, E. M.; Morales-Izquierdo, A. Energy Fuels 2008, 22, 1156–1166. (25) Schneider, M. H.; Andrews, A. B.; Mitra-Kirtley, S.; Mullins, O. C. Energy Fuels 2007, 21, 2875–2882. (26) Speight, J. G. J. Pet. Sci. Eng. 1999, 22, 3–15. (27) Sheremata, J. M.; Gray, M. R.; Dettman, H. D.; McCaffrey, W. C. Energy Fuels 2004, 18, 1377–1384.

Figure 1. Two-dimensional model structures employed in the simulations representing common structural and compositional aspects of petroleum asphaltenes: (a) continental model, (b) archipelago model, and (c) anionic continental model. The model (c) contains a COOterminal group in one of the aliphatic chains but is identical to the continental model (a) in all other aspects. The models explicitly show the heteroatoms N, O, and S.

features of asphaltenes. Three-dimensional models representing continental and archipelago structures, with very similar molecular weights, similar hydrogen/carbon ratios, and identical

Model Molecules of Asphaltene-Like Structures

Energy & Fuels, Vol. 22, No. 4, 2008 2381

Table 1. Different Types of Simulations Conducted in This Studya Type A: Single Solvent + 12 Continental + 12 Archipelago simulation

Nwater

A1 A2 A3

7987

Ntoluene

Nheptane

time (ns)

final size (nm)

306

10 10 10

6.5 × 6.5 × 6.5 5.4 × 5.4 × 5.4 5.3 × 5.3 × 5.3

779 Type B: Mixed Solvent + 12 Continental + 12 Archipelago

simulation

Nwater

Ntoluene

time (ns)

final size (nm)

initial packing

B1 B2 B3

7987 7987 13 965

885 948 714

10 10 10

15.2 × 5.4 × 5.4 15.9 × 5.3 × 5.3 8.4 × 8.4 × 8.4

aggregated in water dispersed in water dispersed in mixture

simulation

Nwater

Ntoluene

time (ns)

final size (nm)

initial packing

C1 C2

10 332 12 241

1453 1047

11.2 34

19.3 × 5.6 × 5.6 12.1 × 7.0 × 7.0

dispersed in water dispersed in mixture

Type C: Mixed Solvent + 12 Anionic Continental + 12 Archipelago

a The number of asphaltene molecules is constant and equal to 24 (12 molecules of each type) in every simulation. For simulations of types A and B, the continental model of the asphaltene was used. For simulations of type C, the anionic continental model was used.

anionic carboxyl group on an aliphatic side chain. Equimolecular mixtures of these model molecules (continental and archipelago) were added to pure solvents (water, toluene, and heptane), and their aggregation, diffusion, and structural changes were studied, employing MD simulations. After this, the aggregation and partitioning behavior of these model molecules in toluene-water systems were studied, focusing on the structure of the aggregates formed in these systems and their orientations and locations relative to the toluene-water interfaces. The key feature of these simulations is the consideration of sufficiently large numbers of molecules and sufficiently large systems. The simulations, although performed employing putative molecular models representing asphaltenes, provide considerable insight regarding their aggregation in pure solvents and structure at oil-water interfaces. The study primarily focuses on the structural features of the aggregates and, hence, presents the results in a qualitative manner based on simulation snapshots. 2. MD Simulation MD simulations were conducted using the simulation freeware GROMACS (version 3.3).30,31 For all MD simulations, the GROMOS96 force field32 was used. Because GROMACS is primarily geared for biomolecular simulations, adaptation of the molecular topologies, simulation box construction, and certain simulation details pertaining to petroleum-based systems requires further elaboration, which are briefly delineated below. 2.1. Molecular Models. Three solvents were employed in the simulations, namely, water, toluene, and heptane. Water was modeled as simple point charge (SPC) molecules. The hydrogen explicit topology for toluene was generated by supplying the coordinates of the molecule to PRODRG.33 The toluene topology explicitly considered the sp2 hybridized carbons for the aromatic ring. For heptane, the united atom model was adopted.34 In addition to the bonded interactions between neighboring sites, the nonbonded

Figure 2. Snapshots of aggregates formed in pure liquid phases after 10 ns: (a) water, (b) toluene, and (c) heptane. The continental model is shown in gray, and the archipelago model is shown as green. All 24 model asphaltene molecules seem to form a single aggregate in all solvents after 10 ns. The stacking of the polyaromatic rings are evident to different extents in all three solvents.

heteroatoms (S, N, and O) were created. A variation of the continental structure was also created, which had a terminal

(28) Rogel, E. Langmuir 2002, 18, 1928–1937. (29) Strausz, O. P.; Lown, E. M. The Chemistry of Alberta Oil Sands, Bitumens and HeaVy Oils; Alberta Energy Research Institute: Calgary, Alberta, Canada, 2003. (30) Berendsen, H. J. C.; Vanderspoel, D.; Vandrunen, R. Comput. Phys. Commun. 1995, 91, 43–56. (31) Lindahl, E.; Hess, B.; van der Spoel, D. J. Mol. Model. 2001, 7, 306–317. (32) Gunsteren, V. The GROMOS home page. GROMOS96 Manual, 2007. (33) Schuttelkopf, A. W.; van Aalten, D. M. F. Acta Crystallogr., Sect. D: Biol. Crystallogr. 2004, 60, 1355–1363, part 8. (34) Kuznicki, T.; Masliyah, J. H.; Bhattacharjee, S. Langmuir 2007, 23, 1792–1803.


Kuznicki et al.

Figure 3. Snapshots of a small region of the aggregate in water (Figure 2a) formed by three archipelago models and one continental model showing the (a) side and (b) top views. The second molecule from the bottom in (a) is the continental model. The axis along which the polyaromatic rings are stacked is depicted in (a) as a dashed line. Each molecule is represented using a different color. N, O, and S are shown explicitly as blue, red, and yellow beads, respectively. In (a), the self-stacking of an archipelago model is shown by A, while the stacking between a continental and an archipelago model is shown by B.

Figure 4. Radial distribution functions of the model asphaltene molecules in water, toluene, and heptane.

Figure 5. Mean separation between two consecutive polyaromatic rings along the axis of stacking in water, toluene, and heptane.

interactions between the sites were described by the Lennard-Jones 6-12 potential. All partial charges for these molecules were used in the topologies based on the GROMOS96 force field parameters.

Because actual asphaltene fractions may have variable amounts of both continental- and archipelago-type molecules, three model molecules representing some key structural and compositional features of asphaltenes of molecular weight around 800 were used in the simulations. These molecules are termed “continental” (C54H65NO2S; MW, 792), “archipelago” (C55H51NO2S; MW, 790), and “anionic continental” (C54H62NO4S; MW, 821) in the rest of this paper. The two-dimensional structures of these molecules are shown in Figure 1. Both the continental and anionic continental models have one large aromatic region per molecule plus several side chains. The only difference in these two models is that the anionic continental molecule contains a terminal COO- group on one of the side chains. The archipelago molecule consists of a long aliphatic ether chain connecting two aromatic ring structures. Three-dimensional structures of these model molecules were created and geometry-optimized using Accelrys Materials Studio molecular modeling software. The coordinates were then exported and submitted to PRODRG to generate the topology files necessary for MD simulations in GROMACS. The alkane chain of the molecules was modeled as united atom structures. All of the aromatic regions were modeled by considering the sp2 hybridized carbons. For the anionic asphaltenes, one sodium ion per molecule was included in this topology. In all simulations, the total number of model asphaltene molecules was kept constant at 24, obtained by mixing 12 molecules each of continental (or anionic continental) and archipelago molecules. The concentration (molarity) of the molecules differed from one simulation to another, depending upon the number of solvent molecules added to the simulation box. 2.2. Simulation Box and Initial Configuration. Three types of simulations were conducted in this work. First, the aggregation of a mixture of continental and archipelago model asphaltenes in pure solvents, namely, water, toluene, and heptane, was studied. These simulations are referred to as simulation type A. Second, aggregation and partitioning of these model asphaltenes in binary solvent systems comprising water and toluene was studied (type B). Finally, aggregation and partitioning of a mixture of archipelago and anionic continental model asphaltenes in a mixture of toluene



Figure 6. (a) Initial and (b) final t ) 10 ns snapshots of the molecular arrangements in simulation B1. The starting configuration in this simulation was an equilibrated aggregate in water based on simulation A1 (Figure 2a).

and water was studied (type C). A summary of the simulations conducted is given in Table 1. In type A simulations, the simulation boxes were constructed by first packing an ordered array of 12 molecules each of continental and archipelago asphaltene molecules (24 total) in a cubic box of edge 6.5 nm. After this, the box was filled randomly with the appropriate solvent. The simulation box dimensions reported in Table 1 represent the final equilibrated box. Simulations with heptane (A3, Table 1) were performed only to ensure that the asphaltenes aggregated in heptane to a larger extent than toluene. Accordingly, we did not conduct any binary solvent simulations involving heptane in this study. Simulations of type B were performed using water and toluene as the solvents, with 12 continental and 12 archipelago molecules added to this solvent mixture. For these simulations, three approaches were employed to construct the simulation box. In simulation B1 (Table 1), the equilibrated structure from simulation A1 (aggregated asphaltenes in water) was taken as the initial configuration, after which the box was expanded in the x direction and toluene molecules were added to the empty volume. Toluene molecules were initially added as a regular cubic array. We refer

to this type of box construction as “aggregated in water” in Table 1. In simulation B2, the asphaltene molecules were initially placed in a regular cubic lattice of edge 7 nm. After this, the box was filled with water. Then, the box was expanded in the x direction, and toluene molecules were added to the empty volume. This type of construction is referred to as “dispersed in water”. In simulation B3, the asphaltenes were placed in a cubic lattice of edge 8 nm, after which the box was filled with toluene. This system was energyoptimized for about 1 ns. Then, the entire box was filled up randomly with water. This box construction is referred to as “dispersed in mixture” because the asphaltene molecules in this geometry are initially in random contact with both water and toluene molecules. For type C simulations, the primary objective was to study the role of anionic terminal groups in one fraction of the asphaltenes. For these simulations, the 12 continental molecules were replaced by 12 anionic continental molecules (Figure 1c). No changes were made to the archipelago molecules. The initial configuration in simulation C1 was created employing a procedure similar to simulation B2. The initial configuration for simulation C2 was generated employing a procedure similar to simulation B3. In the


Kuznicki et al.

Figure 7. Snapshot of the model asphaltene molecules at t ) 10 ns in simulation B2. The molecules were initially dispersed in water. After 10 ns, the molecules migrate to the toluene phase. The toluene molecules are not shown to highlight the aggregate structure. The vertical dashed lines show the approximate locations of the water-toluene interfaces. The model asphaltene molecules remain completely within the toluene phase and are not at the water-toluene interfaces.

type C simulations, the box was rendered electroneutral by adding 12 Na+ ions (one for each COO-). Although these simulations were conducted for a longer duration (simulation C2 for 34 ns), we report results in this paper based on the first 10 ns of the run. The overall packing of molecules for the starting configuration in aqueous systems was performed by closely matching the mixture density to the experimental density of water (996 kg/m3). Energy profiles were examined to confirm the approach to equilibrium. During the model construction and subsequent simulations, 3D periodic boundary conditions were applied. Details of simulation box construction for the binary solvent systems employed in this study are available elsewhere.34 It should be noted that the initial box dimensions are not important in these simulations, because we employ temperature and pressure coupling, which modify the box volume to attain the appropriate densities and compressibilities of the solvents. The system energy and density were monitored during the simulations to ensure equilibration of the box. 2.3. Simulation Procedure. The simulations were performed in the isothermal-isobaric (NPT) ensemble at temperature T ) 298 K and pressure P ) 105 Pa (1 bar). Because the initial asphaltene molecular models were energy-optimized, static structure optimization was not performed. The MD simulations were conducted for a duration of at least 10 ns with a time step of 2 fs. Weak coupling of the molecules to a solvent bath of constant temperature was maintained using the Berendsen thermostat, with a coupling constant of τT ) 0.1 ps.35 The pressure was controlled using the Berendsen algorithm, with a pressure coupling constant of τp ) 1 ps, applying a compressibility of 4.5 × 10-5 kJ-1 mol nm3. Unless otherwise stated, the simulation parameters were generally maintained at the default values offered in GROMACS.36 For simulations involving (35) Berendsen, H. J. C.; Postma, J. P. M.; Vangunsteren, W. F.; Dinola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684–3690. (36) Spoel, D. V.; van Buuren, A. R.; Apol, E.; Meulenhoff, P. J.; Tieleman, D. P.; Sijbers, A.; Feenstra, K. The GROMACS Manual (version 3.3); Groningen University: Groningen, The Netherlands, 2004.

Coulomb interactions, particle mesh Ewald (PME) summation was employed to account for the long-range electrostatic interactions. The computations were quite intensive, owing to the large box sizes and the large number of atoms involved in simulations of types B and C. However, a large system was desirable, owing to the fact that we intended to minimize the influence of periodic boundary conditions and explore the aggregation behavior of a relatively large number of asphaltene molecules. A parallel version of the MD simulation program was employed on a computational cluster with four nodes. Each processor was allocated 8 GB of memory. To our knowledge, the type C simulations represent the largest box dimensions and possibly the largest times over which asphaltene-oil-water systems have been studied employing MD. Trajectories were used to compute the density distributions, system energy, radial distribution functions, mean square displacements of different molecules, and other associated parameters for the system. These were computed using appropriate postprocessing programs available in GROMACS,36 as well as using the visualization program VMD.37

3. Asphaltene Properties in Pure Solvents The aggregation behavior of the model asphaltenes in three different solvents, namely, water, toluene, and heptane is described in this section. For these simulations, we employed 12 molecules each of the continental and archipelago models (see parts a and b of Figure 1), all well-dispersed in the solvent as the initial configuration. Figure 2 depicts the overall aggregate structures formed after 10 ns in these three solvents. When dispersed in water, the model asphaltene molecules tend to form a droplet-like aggregate. The rate of aggregation (37) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graphics 1996, 14, 33.



Figure 8. Snapshots of the model asphaltene molecules in the toluene-water system for simulation B3. (a) Front view of the box at 10 ns, showing the aggregate inside the cylindrical droplet of toluene formed in water. All model asphaltene molecules are inside the toluene phase. The toluene molecules are not shown in this snapshot. (b) Side view of the simulation box at 10 ns, showing the aggregate formed in the toluene droplet. The water molecules surrounding the toluene droplet are not shown in this snapshot.

in water is quite high. It was found that the droplet-like structure depicted in Figure 2a is attained within 1.5 ns of simulation time. After this, the structure remains essentially stable during the remainder of the simulation (up to 10 ns). In toluene, the molecules appear to be slightly more dispersed than in water, although the presence of an aggregate is apparent from Figure 2b. In heptane, the model asphaltenes tend to coalesce to form larger aggregates, as shown in Figure 2c. In the case of toluene

or heptane, the aggregation rate is not as rapid as observed in water, with aggregation rates being the slowest in heptane. A salient feature of the aggregate structures in Figure 2 is the existence of partial ordering of the polyaromatic ring regions of the model asphaltenes into stacks. The stacking behavior is most pronounced in heptane followed by water and toluene. This stacking of aromatic rings can be attributed to π-π interactions between molecules containing aromatic moieties.


Kuznicki et al.

Figure 9. Partitioning and aggregation behavior of a mixture of anionic continental and archipelago molecules in toluene-water mixtures. (a) Simulation C1 and (b) simulation C2. The anionic continental molecules are shown in gray with the oxygen atoms shown as red dots. The archipelago molecules are rendered in green. The ratio of toluene/water is different in (a) and (b), yielding different thicknesses of the toluene films. The free sodium ions are shown as blue spheres and are only present in the aqueous phase.

It is well-known that asphaltenes aggregate in organic solvents even at very low concentrations. In aromatic solvents, such as toluene, aggregation at concentrations as low as 50 mg/L has been reported.38,39 Our simulations, which represent very high concentrations of model asphaltenes in the box, clearly indicate their aggregation in toluene, as well as stacking of the polyaromatic rings. Figure 3 depicts a snapshot of four model molecules in water from Figure 2a. It is apparent that the archipelago model partakes in different types of stacking, including self-stacking of the two aromatic rings of a single molecule (indicated by A in Figure 3a) and stacking with the aromatic ring of the continental molecule (shown as B in Figure 3a). Five aromatic rings are stacked in this configuration. It is observed that the rings do not exactly superimpose on each other (Figure 3b), and there are slight offsets in the arrangement. From the snapshot, it also appears that neither the ring nor side-chain heteroatoms play a significant role in the stacking process. In all our simulations, (38) Acevedo, S.; Ranaudo, M. A.; Pereira, J. C.; Castillo, J.; Fernandez, A.; Perez, P.; Caetano, M. Fuel 1999, 78, 997–1003. (39) Goncalves, S.; Castillo, J.; Fernandez, A.; Hung, J. Fuel 2004, 83, 1823–1828.

we observed a slightly lower tendency of the continental molecules toward stacking, which might be attributed to the steric interference of the aliphatic side chains surrounding the polyaromatic ring. Finally, in all cases, the solvent molecules are excluded from regions between the stacked asphaltenes. Figure 4 shows the radial distribution function, g(r), of the model molecules in water, toluene, and heptane. The distribution function was based on the equilibrated aggregate structure from the last 2 ns of simulations (from 8 to 10 ns). These plots exhibit two to three peaks. This indicates a strong short-range correlation between the molecules, supporting the aggregation process observed in Figure 2. The shape of the radial distribution functions is similar in all cases. However, the tendency to aggregate is much more pronounced in the case of water compared to heptane or toluene. Furthermore, the asphaltenes aggregate to a larger extent in heptane than in toluene. The appearance of the first peak of g(r) is due to the compact association of the polyaromatic rings through stacking within the aggregate. Interestingly, the position of the first peak in g(r) is virtually identical in all three solvents. This was further verified by explicitly calculating the mean distance between two stacked aromatic rings from the simulation trajectories. As



Figure 10. Detailed aggregate structure at the toluene-water interface in simulation C2. The anionic continental molecules (gray) assemble with their terminal charged atoms tethered to the interface. All other molecules aggregate around these tethered molecules, forming a complex 3D structure. The archipelago model molecules (green) form stacks that extend into the bulk toluene phase. The dashed line shows the approximate location of the water-toluene interface. The toluene molecules are not shown in the image.

shown in Figure 5, the distance between two stacked polyaromatic rings is indistinguishable, with a mean separation of about 0.4-0.45 nm in all three solvents. This implies that none of the solvents penetrate the stacked aromatic planes, owing to volume exclusion effects. Another interesting conclusion from these calculations is that aggregates comprising continental or archipelago models exhibit very similar structural features. The mean squared displacement (MSD) of various molecules (solvent and asphaltenes) in the simulations of type A were evaluated over a 10 ns interval. From the slopes of the MSD, the diffusivities of various species were evaluated. The diffusivities of water, toluene, and heptane were obtained as 4.05 ( 0.04 × 10-9, 0.89 ( 0.05 × 10-9, and 6.29 ( 0.39 × 10-9 m 2/s, respectively. The diffusivity of the solvents were comparable to values reported in the literature.40,41 The diffusivities of the continental and archipelago model molecules were found to be 4.0 ( 1.2 × 10-10 and 4.5 ( 1.6 × 10-10 m2/s in water. In toluene, the corresponding diffusivities were 1.2 ( 0.8 × 10-10 and 7.0 ( 2.0 × 10-11 m2/s. For these hypothetical molecules, it appears that the calculated diffusivities, although subject to considerable uncertainties, (40) Bhide, S. Y.; Berkowitz, M. L. J. Chem. Phys. 2005, 123, 224702. (41) Lee, S. H.; Lee, H.; Pak, H. Bull. Korean Chem. Soc. 1997, 18, 478–484.

owing to the small number of molecules (12) employed in accumulation of the averages, are qualitatively reasonable. First, the diffusivity of these molecules are about 1 order of magnitude lower than the solvent diffusivities. Second, the diffusivities of the two molecules are quite similar in water, which is commensurate with their similar molecular weights. Third, the diffusivity of the archipelago model is lower in toluene compared to the corresponding diffusivity of the continental molecule, a phenomenon that can be attributed to the interactions between the multiple polyaromatic rings of these molecules and toluene. Finally, one can compare the results with theoretical values of diffusivity based on a modified Stokes-Einstein equation of the form42 D)

kBT 6πµ(0.676〈R2 〉1⁄2)

(1)

where kB is the Boltzmann constant, T is the absolute temperature, µ is the solvent viscosity, and 〈R2〉1/2 is the root-meansquare radius of gyration of the molecule. Using the radius of gyration of the asphaltene molecules in our simulations, we obtain a diffusion coefficient of the continental molecules based (42) Cussler, E. L. Diffusion: Mass Transfer in Liquid Systems, 2nd ed.; Cambridge University Press: Cambridge, U.K., 1997.


on eq 1 to be 3.2 ( 0.6 × 10-10 m2/s in water. It should be noted that diffusivities of macromolecules are concentrationdependent and can either increase or decrease with concentration. In toluene, diffusivities of polymers, such as polystyrene, increase with concentration.42 Similarly, in water, diffusivities of proteins can increase with concentration.43 The model asphaltene concentrations in our simulation box are relatively high (about 0.34-1.2 g/L). These results indicate that the bulk transport characteristics of the molecules were adequately represented in the 10 ns simulations reported here. 4. Aggregation and Partitioning in Toluene-Water Mixtures We now consider the aggregation and partitioning of the model continental and archipelago molecules in binary solvent mixtures of toluene and water. The different simulations conducted for these studies are shown in Table 1 as type B simulations. Figure 6 depicts the initial and final structures of the aggregate in simulation B1. The final aggregate, obtained after 10 ns of simulation, did not change its structure appreciably. It remained suspended in the aqueous phase, never venturing into the toluene phase. The diameter of the roughly cylindrical aggregate was around 4 nm, and very small changes in the stacking of the polyaromatic rings, if at all, were observed during the intermediate stages of the run. To summarize, the aggregated molecules, when suspended in water, do not seem to migrate to the toluene phase within the time frame of the simulation. Furthermore, these aggregates do not show significant affinity for the toluene-water interface because the mean separation between the centroid of the aggregate and the nearest toluene-water interface increased during the simulation. The aggregate structure contains a core region formed by stacked polyaromatic rings surrounded by the aliphatic chains.19 It was observed that polyaromatic rings in the archipelago molecules can fold very efficiently and enhance the formation of such compact aggregates. Furthermore, it is evident that the surface of these aggregates, once formed, is relatively insensitive to the polarity of the solvent. Figure 7 shows a snapshot of the simulation box in simulation B2 at t ) 10 ns. During the simulation, the solvents partitioned to yield alternating layers of toluene and water phases. The toluene molecules are not shown in the figure to clearly depict the aggregated structure of the model molecules. The toluene-water interfaces are shown by the vertical dashed lines. In this case, we observe that every model molecule that was initially dispersed in water eventually migrated into the toluene phase. Within the toluene layer, the molecules aggregated through stacking of the polyaromatic rings. The aggregate remains suspended in the central region of the toluene layer and does not remain associated with the water-toluene interface. In fact, during the intermediate stages of the simulation, it was observed that the individual molecules migrate across the interface quite readily, without being trapped at the interfaces. When Figures 6 and 7 are compared, it becomes evident that individual asphaltenes can diffuse quite readily across water-toluene interfaces, in which case they prefer to migrate to the toluene phase, but the bulkier aggregates have a much smaller diffusivity and show a lower preference for the toluene phase. In parts a and b of Figure 8, snapshots of the final configuration of the molecules (at t ) 10 ns) during simulation B3 are shown. In this case, the molecules were initially dispersed (43) Bhattacharjee, S.; Kim, A. S.; Elimelech, M. J. Colloid Interface Sci. 1999, 212, 81–99.

Kuznicki et al.

in toluene and water was added randomly to the box, filling up the void regions of the box that could accommodate these molecules. Consequently, the initial structure represented a mixture of solvent (water and toluene) molecules. It is evident from the front view in Figure 8a that the two solvents partition readily, forming a distinct circular cylindrical droplet of toluene dispersed in water. The toluene molecules in the droplet are not shown to enable viewing of the model asphaltenes. The cylindrical droplet containing the model asphaltenes is more clearly shown in Figure 8b, which depicts the toluene molecules in the side view of the simulation box. When the snapshots of parts a and b of Figure 8 are compared, it is immediately apparent that the model molecules remain in the toluene phase, forming a cylindrical aggregate surrounded by toluene molecules on all sides. The overall structure resembles a rodlike droplet of toluene with a supporting backbone of molecules that have aggregated predominantly by stacking of their polyaromatic rings. In all simulations of type B, we encountered the somewhat surprising behavior of the chosen model molecules that they were never associated with the toluene-water interface. Instead, these molecules always preferred to remain in the bulk liquid, preferably toluene. When these model molecules are already present in an aggregated state in water, they do not seem to show any particular affinity to migrate to the toluene phase or the toluene-water interface. These observations might appear anomalous when considering the general notion in petroleum engineering that asphaltenes are surface-active molecules that impart structure and rigidity to the oil-water interfaces. Clearly, the continental and archipelago models used in the simulations do not exhibit this trend. This leads us to explore whether polar or charged moieties in some of these model molecules can render them affine to the water-oil interface. 5. Aggregation and Partitioning Assisted by Polar (Anionic) Moieties To study the influence of charged groups on the side chains of model asphaltene molecules, we conducted the type C simulations (Table 1), where the continental model molecules were modified to add a terminal carboxylic group on one of the side chains (Figure 1c). Each terminal carboxylic group was deprotonated to render it anionic. A total of 12 molecules of these anionic continental model were then mixed with 12 molecules of the archipelago model to study their aggregation behavior in water-toluene mixtures. For type C simulations, the initial configuration always contained individual molecules dispersed on a regular cubic lattice, after which the solvents were added to the box. In simulation C1, the molecules were initially solvated in water, followed by the addition of toluene. In simulation C2, the molecules were solvated by randomly adding both water and toluene to the simulation box. It should be noted that, in these simulations, the explicit anionic charges needed to be neutralized by adding 12 sodium ions (Na+) to the initial configuration. Also, PME summation was employed in the calculations to address the long-range nature of the nonbonded electrostatic interactions. Parts a and b of Figure 9 depict the simulation boxes at t ) 10 ns for simulations C1 and C2, respectively. These snapshots represent a film of toluene surrounded by two layers of water on both sides, similar to the equilibrated system shown in Figure 7. The sodium ions are shown in these images as blue spheres. Once again, the toluene molecules are not shown in these snapshots. The box dimensions are much larger in these simulations compared to simulations of type A and B. Further-


more, the thickness of the toluene film are different in simulations C1 and C2. In both simulations, the two solvents form distinct phases, with a large number of model asphaltenes present in the toluene film. In simulation C1 (Figure 9a), few of the model asphaltene molecules also form a nanoaggregate in the aqueous phase. Notably, in this simulation, all of the model molecules were dispersed only in the aqueous phase in the starting configuration. In simulation C2 (Figure 9b), where the molecules were initially randomly solvated by toluene and water, they all migrated into the toluene phase. The much thinner toluene film contains a closely packed aggregate of model asphaltenes. A remarkable feature of the aggregation and partitioning process evident in the snapshots in Figure 9 is the presence of the anionic terminal groups of the continental molecules at the toluene-water interfaces. In both simulations C1 and C2, we observe that these anionic continental molecules preferably associate with the interface, tethered through the deprotonated carboxyl groups. After the initial capture of these groups at the interface, stacking of polyaromatic rings of different model asphaltenes result in the formation of a complex 3D aggregate structure. A key observation in this case is that the anionic continental molecules are also instrumental in bringing the archipelago molecules to the interface. It is discernible that, in real systems, such an associated structure can easily provide considerable rigidity to the oil-water interface. A more detailed view of the aggregate is shown in Figure 10, which depicts a small region of the toluene-water interface corresponding to simulation C2. The toluene molecules are not shown explicitly. All of the model asphaltene molecules are present in the toluene phase. The oxygen atoms in the anionic continental molecules are depicted as red spheres. Only the terminal carboxyl groups (represented by two closely placed terminal oxygen atoms) of these molecules protrude into the aqueous phase through the interface (dashed line). The carboxylic aliphatic chains act as tethers, holding the continental molecules near the interface. The polyaromatic rings of these molecules can self-assemble through stacking, and as evident in Figure 10, the stacking of the rings can extend several layers deep into the toluene phase. Both uncharged and anionic molecules can participate in the stacking aggregation. It is therefore evident that charged groups, preferably located on the aliphatic chains of the molecules, play a key role in ascribing surface activity to them. The stacking of polyaromatic rings of these molecules does not occur at the interface but within a particular fluid phase. The stacking can lead to complex 3D aggregate architectures that can impart considerable rigidity to the oil-water interface. However, such a mechanism of interface


stabilization cannot be observed in the absence of charged groups on the model molecules. 6. Conclusions The MD simulations conducted in this work throw considerable light into the mechanism of aggregation of model molecules with asphaltene-like structures in pure solvents, as well as their aggregation and partitioning in binary solvent systems. The results presented are mostly qualitative, based on simulation snapshots, but provide graphic details of the behavior of these polyaromatic molecules in oil-water systems. Our simulations show that these model molecules tend to form aggregates in pure solvents, which resemble asphaltene nanoaggregates reported in several earlier studies. These aggregates are also less surface-active, particularly when no charged terminal groups are involved. The polyaromatic rings of both the continental and archipelago models prefer to stack in the toluene phase in binary solvent systems. The presence of ionic terminal groups on the aliphatic chains of specific asphaltene fractions can dramatically enhance the interfacial activity of the entire asphaltene population. These results indicate that the interfacial activity of these molecules are not ascribable to their polyaromatic ring regions but more likely to their terminal group chemistries. However, one should bear in mind that these simulations can only capture high-frequency (very fast) dynamics of such systems, and these results may not be indicative of their long-term behavior. Asphaltenes represent a very complex solubility class, and molecules with quite different chemical structures may be part of it. It is likely that only small fractions of these structures will behave as surface-active molecules. If this is the case, studies employing a single “average” model molecule based on bulk properties of asphaltenes will reflect mainly the characteristics of the larger fraction in a mixture, and the surface active properties will not be present in the resulting average structure. To address interfacial activity of asphaltenes, it is therefore much more pertinent to develop molecular models obtained from asphaltenes isolated from the films formed at the interfaces. Acknowledgment. The authors gratefully acknowledge financial support from the Canada Research Chairs (CRC) program, Natural Sciences and Engineering Research Council of Canada (NSERC), and Alberta Energy Research Institute (AERI). We are also grateful to the referees for their valuable comments and suggestions, which were immensely insightful and helpful in modifying our conclusions from the work. EF800057N

Molecular Dynamics Study of Model Molecules Resembling

Recommend Documents