Simple Simulation Model for Exploring the Effects of Solvent and

Article Cite This: J. Phys. Chem. B 2019, 123, 6111−6122

pubs.acs.org/JPCB

Simple Simulation Model for Exploring the Effects of Solvent and Structure on Asphaltene Aggregation Nicholas J. H. Dunn,† Besha Gutama, and W. G. Noid* Department of Chemistry, The Pennsylvania State University, University Park, Pennsylvania 16802, United States

Downloaded via BUFFALO STATE on July 25, 2019 at 17:21:08 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: Asphaltenes are operationally defined as the fraction of crude oil that is soluble in toluene but insoluble in n-heptane. According to the Yen−Mullins model, typical asphaltenes are relatively small molecules consisting of a single aromatic core flanked by aliphatic chains. The Yen−Mullins model posits that asphaltene aggregation proceeds via a hierarchical mechanism involving small nanoaggregates with stacked aromatic cores surrounded by a corona of aliphatic tails. In this work, we introduce a coarse-grained (CG) model for investigating the physical picture underlying the Yen− Mullins model and, more generally, the effects of the solvent character and molecular structure upon asphaltene selfassembly. By representing proposed asphaltenes in united atom detail, this CG model accurately describes their shape and conformational properties. Conversely, the CG model mimics varying solvent conditions by modulating the effective attraction between aliphatic and aromatic groups. Given the simplicity of this model, we performed long, replicate simulations of 147 different asphaltene solutions. As proposed by the Yen−Mullins model, island-type molecules readily form stacked aggregates under conditions that promote aromatic interactions. Interestingly, the onset of nanoaggregation appears to be insensitive to the aliphatic tails, although these tails may sterically stunt further growth of nanoaggregates. Consequently, nanoaggregates form more readily and grow larger under conditions that promote both aliphatic and aromatic interactions. In contrast, archipelago-type molecules also form large aggregates, but they do not demonstrate significant stacking interactions. Thus, the CG model reasonably describes the physical intuition of the Yen− Mullins picture and may prove to be useful for exploring later stages of asphaltene aggregation.

■

INTRODUCTION Asphaltenes are a solubility class of molecules that have been termed the “cholesterol of petroleum” because of their propensity for aggregating in petrochemical mixtures.1,2 It is generally believed that asphaltenes exist as a colloidal suspension in crude oil and that they flocculate into viscoelastic masses when this suspension is destabilized by changes in processing conditions.3−6 The resulting floccs foul reservoirs, impede heat exchange equipment, and clog pipelines.5−7 It is estimated that asphaltene flocculation costs the petroleum industry billions of dollars each year.2,8,9 The precise molecular properties of asphaltenes have been extensively debated.10−17 Asphaltenes are operationally defined as the fraction of crude oil that is soluble in toluene and insoluble in n-heptane.2,7,11,18 Consequently, a given asphaltene sample may contain many diverse species, which considerably complicates the chemical analysis of their molecular properties.19−23 Furthermore, because of their high propensity for aggregation, it is challenging to distinguish between the properties of single molecules and those of asphaltene aggregates.24 As a result, two competing models have developed for characterizing asphaltene molecules.14,25,26 The continental (or island) model proposes that typical © 2019 American Chemical Society

asphaltene molecules are characterized by a single aromatic core surrounded by alkyl tails.17,27 Conversely, the archipelago model proposes that typical asphaltene molecules are characterized by several aromatic cores linked by alkyl chains.28,29 Recently, the continental model has gained considerable popularity,30−36 although the asphaltene fraction may well include both archipelago and continental structures.37,38 The Yen−Mullins (or modified Yen) model has emerged as a leading theory for describing representative asphaltene molecules and their aggregation.23,39−41 The Yen−Mullins model posits that typical asphaltenes are relatively small molecules (approximately 750 g/mol) with a single aromatic core surrounded by aliphatic tails.10 The Yen−Mullins model then describes asphaltene aggregation via a hierarchical mechanism that is driven by regular solution considerations.41−45 Under very dilute conditions of 100 mg/L asphaltenes in oil,46−51 asphaltenes self-assemble to form short, rod-like nanoaggregates of 7−10 molecules that stack to Received: May 6, 2019 Revised: June 21, 2019 Published: June 24, 2019 6111

DOI: 10.1021/acs.jpcb.9b04275 J. Phys. Chem. B 2019, 123, 6111−6122

Article

The Journal of Physical Chemistry B sequester their aromatic cores from a poor aliphatic solvent, while preserving their solubility via a corona of exposed aliphatic tails.41 At slightly higher concentrations of 2−5 g/ L,52−54 these nanoaggregates then cluster together in a fractal arrangement.55,56 This model is supported by experimental studies employing varied techniques such as two-step laser desorption ionization,30 NMR,31,32,49 alternating and direct current conductivity,46−48 and centrifugation.50,51 Given the difficulties in experimentally characterizing asphaltenes, molecular simulations provide a powerful tool for studying their behavior. In particular, all-atom (AA) molecular dynamics simulations have examined the nanoaggregation of model asphaltenes in various solvents, at interfaces, and under various thermodynamic conditions.20,22,57−67 Unfortunately, the computational expense of atomically detailed simulations necessarily limits the duration, size, and scope of these studies, although recent computing advances with graphics processing units have somewhat extended these limits.68 Moreover, in the dilute conditions relevant for nanoaggregation, the vast majority of computational effort is expended in the simulating solvent. Consequently, to observe significant aggregation on the length- and timescales that are computationally accessible to AA simulations, some studies have employed relatively high concentrations of asphaltenes and relatively small systems which may artificially bias the simulated aggregation mechanism or equilibrium morphology. The computational expense of AA models has motivated considerable interest in coarse-grained (CG) models for soft materials.69,70 By representing systems in reduced detail, CG models provide the necessary efficiency not only for simulating larger systems and longer timescales but also for obtaining improved statistics and for exploring more diverse conditions. Indeed, prior simulation studies have employed CG models to examine various aspects of asphaltene phase behavior and selfassembly.71−79 In this work, we introduce and investigate a very simple and highly efficient “toy” CG model for examining generic aspects of asphaltene phase behavior and self-assembly.80−82 This toy model is particularly well suited for exploring the regular solution considerations underlying the Yen−Mullins model, that is, that the hierarchical aggregation of asphaltenes is driven by a microphase separation of aromatic cores from a poor aliphatic solvent. Moreover, this model preserves many details of asphaltene structure, interactions, and fluctuations that are difficult to describe with analytic mean field or solution theory approaches.83 As illustrated in Figure 1, this model adopts a united atom description for each asphaltene molecule, while distinguishing between aliphatic and aromatic atoms. Because it accurately describes the shape and flexibility of asphaltene molecules, the model likely provides a reasonable description of their packing in aggregates. Conversely, the model adopts an implicit, mean field treatment of the surrounding solvent. Specifically, the model describes various solvent environments by varying the effective attraction between aliphatic and aromatic asphaltene atoms. Thus, the model is distinguished from prior CG models by combining a relatively detailed description of asphaltene molecules with an implicit, mean field description of the solvent. The simplicity and efficiency of the model enabled two long simulations for each of 147 different model asphaltene solutions. These simulations provide statistical sampling that may be conservatively estimated as equivalent to many milliseconds of conventional

Figure 1. Structures of model asphaltene compounds. Green and cyan sites indicate aliphatic tail, L, and aromatic core, R, sites, respectively. Panel (a) presents the ovalene-8 model, which features an ovalene aromatic core and is representative of island, or continental, asphaltenes. Panel (b) presents the bipyrene-8 model, which features a bipyrene aromatic core and is representative of archipelago asphaltenes. Both molecules have 32 R sites and 4 aliphatic tails each containing l = 8 L sites.

MD simulations with atomically detailed models. As a result, the simulations effectively characterize the influence of the molecular architecture and solvent conditions upon the stability and morphology of asphaltene aggregates. The remainder of the paper is organized as follows. The Methods section describes the CG model and simulations, as well as the metrics employed for analyzing the simulations. The Results section presents the results of these simulations. The Discussion section analyzes these results in the context of the Yen−Mullins model and other recent studies. Finally, the Conclusions section summarizes our main findings and suggests possible future directions.

■

METHODS CG Model. The CG model adopts a generic united atom representation to describe representative asphaltene molecules with two types of CG “sites.” The model represents each aliphatic carbon with an equivalent L site and each aromatic carbon with an equivalent R site. The L and R sites loosely correspond to aliphatic sp3 carbons and aromatic sp2 carbons, respectively, in the OPLS-AA potential.84 Figure 1 illustrates this CG representation for two representative asphaltene molecules, while representing L and R sites in green and teal, respectively. The model employs OPLS-AA potentials to describe the conformational flexibility of asphaltene molecules. The model rigidly constrains the distance between each pair of bonded sites to be 0.14 nm for R−R bonds, 0.151 nm for R−L bonds, and 0.1529 nm for L−L bonds, which correspond to the sp2− sp2, sp2−sp3, and sp3−sp3 bonds in OPLS-AA, respectively. Bond bending and torsional degrees of freedom are modeled with the OPLS-AA angle and dihedral potentials for the corresponding atom types. Consequently, the model preserves the rigidity and planarity of aromatic cores, as well as the flexibility and conformational tendencies of the aliphatic tails. Importantly, the CG model does not explicitly describe the solvent. Rather, the influence of the solvent is approximated in a simple mean field manner by modulating the effective attraction between pairs of asphaltene sites. We derive the site−site potentials from the standard Lennard-Jones (LJ) 12-6 potential ÅÄÅ ÑÉ Åi σ y12 i σ y6ÑÑ ULJ(r ) = 4ε0ÅÅÅÅjjj zzz − jjj zzz ÑÑÑÑ ÅÅk r { k r { ÑÑÑÖ (1) ÅÇ 6112


Article

The Journal of Physical Chemistry B where σ = 0.3525 nm and ε0 = 0.28451 kJ/mol, which correspond to the mean of the corresponding parameters for sp2 and sp3 carbons in the OPLS-AA potential. As illustrated in Figure 2, we employ the Weeks−Chandler−Anderson

asphaltene-solvent solution that is obtained as a function of the asphaltene coordinates after formally tracing over the solvent coordinates.86,87 In particular, the simulated pair potentials do not simply quantify the energy of direct interactions between asphaltene molecules but also incorporate temperature-dependent indirect and entropic contributions resulting from the solvent molecules that are not explicitly described by the model. CG Simulations. We performed all simulations with GROMACS 4.5.3 compiled with double floating point precision and an integration timestep of 2 fs.88 All simulations sampled the constant NVT ensemble with three-dimensional periodic boundary conditions. We maintained a constant temperature T = 303 K with the stochastic dynamics thermostat, while employing a coupling constant of 0.5 ps.89 We employed the LINCS algorithm to rigidly constrain all bonds.90 We implemented the nonbonded potentials in eq 5 by employing GROMACS table files with Ur and Ua specified by separate columns, such that Ua could be independently scaled to mimic the varying solvent quality. These table files ensured that the simulated forces smoothly vanished at a cutoff of 1.2 nm. Configurations were sampled from these simulations after every 1000 timesteps. We performed simulations of model asphaltenes with ovalene and bipyrene cores as representative examples of continental- and archipelago-type architectures, respectively. As illustrated in Figure 1, we decorated these cores with four symmetrically arranged aliphatic tails. We defined six model asphaltene molecules by considering tails with lengths, l, of 5, 8, or 11 carbons. We simulated each of the six model asphaltene molecules in a range of solvent conditions by varying the parameters (εR, εL). For each model molecule and each solvent condition considered, we simulated 25 identical asphaltene molecules in a volume that corresponded to a fixed 1.0 wt % (i.e., 6.8 g/L), assuming that the implicit solvent had the same density as heptane. For each system, we performed two independent replicate simulations that started from different initial configurations. The starting configuration for one replicate was generated by placing the molecules on a regular grid within the simulation box. The starting configuration for the second replicate was obtained after randomizing the initial lattice configuration during a short simulation with purely repulsive potentials. Unless otherwise specified, all reported results are obtained from averaging over these two replicate simulations. We simulated each system for a minimum of 60 ns and maximum of 400 ns. We terminated 150 of the 294 simulations prior to the 400 ns limit because they had either achieved a “stable state” or become “kinetically trapped” in a disordered aggregate. We deemed that a simulation had achieved a stable state when (1) the aggregate distributions for both replicates remained approximately constant during the course of the preceding 50 ns (i.e., the extent of aggregation, G, changed by less than 2.5 during the 50 ns simulation window) and (2) both replicates demonstrated qualitatively similar aggregate distributions. Conversely, we deemed that a simulation had become kinetically trapped when (1) the aggregate distributions remained approximately constant during the course of the preceding 50 ns and (2) the two replicates demonstrated qualitatively different aggregate distributions that indicated the formation of large irreversible aggregates. The Supporting Information compares simulated replicates for two systems and

Figure 2. Weeks−Chandler−Anderson decomposition of the LJ nonbonded pair potential, ULJ, into repulsive, Ur, and attractive, ε0Ua, components. The solid black, dashed red, and dashed-dotted blue curves present ULJ, ε0Ua, and Ur, respectively.

decomposition85 to split ULJ into purely attractive, Ua, and repulsive, Ur, components ULJ(r ) = Ur(r ) + ε0Ua(r )

(2)

l o oULJ(r ) + ε0 if r < r0 Ur(r ) = o m o o o if r ≥ r0 n0

(3)

where

l− 1 if r < r0 o o o o o Ä É ÅÅ Ñ Ua(r ) = o m ÅÅji σ zy12 ji σ zy6ÑÑÑ o o ÅÅjj zz − jj zz ÑÑ if r ≥ r0 4 o o ÅÅ r Ñ o o k r { ÑÑÑÖ n ÅÅÇk {

(4)

and r0 = 2 σ is the minimum of the original LJ potential. We then employ the potential 1/6

U2(r ;ε) = Ur(r ) + εUa(r )

(5)

to model the interaction between each “nonbonded” pair of sites, that is, between each pair of sites that are in distinct molecules and between each intramolecular pair that are separated by three or more bonds. We model a particular solvent environment by defining a pair of parameters (εR, εL), which specify the influence of the solvent upon the effective attraction between R and L sites. Each nonbonded pair of R sites interacts according to the potential URR(r) = U2(r;εR), whereas each nonbonded pair of L sites interacts according to the potential ULL(r) = U2(r;εL). “Cross-interactions” between nonbonded pairs of R and L sites are modeled by the potential URL(r) = U2(r;εRL) with the geometric mixing rule: εRL = εR εL . Importantly, because Ur is independent of ε, the magnitude of the attractive interactions does not alter the excluded volume interaction between the pair. Thus, each site remains in the same “size” and each molecule remains in the same “shape”, irrespective of the solvent environment. Finally, in closing this subsection, we emphasize again that this model does not explicitly describe solvent molecules but only the asphaltene solutes. The simulated potential is a very crude approximation to the free energy for the entire 6113


Article

The Journal of Physical Chemistry B

core attraction and those that are stabilized by tail−tail or core−tail attraction. Maximin Core Distance. Given two ovalene molecules, A and B, with aromatic cores, ( and ) , and aromatic sites that are indexed, a ∈ ( and b ∈ ) , respectively, we define the maximin core distance, rmaximin

illustrates the criteria for terminating simulations short of 400 ns. We initially selected solvent conditions to qualitatively explore the solvent space and range of aggregation behavior. After our initial simulations, we selected additional state points to delineate regions of solvent space demonstrating distinctive behavior. We performed a total of 294 simulations for 147 systems with a total production duration of 14.7 μs. Throughout the paper, we report times in terms of simulated nanoseconds based on the fixed integration timestep, dt = 2 fs, and the number of simulated timesteps. However, the elimination of solvent dramatically enhanced the sampling efficiency of the present CG simulations. Previous simulation studies indicated that asphaltene molecules diffuse through heptane or toluene with a diffusion constant on the order of 10−5 cm2/s.76 However, similar asphaltene molecules diffused through the present implicit solvent with observed selfdiffusion constants on the order of 1 cm2/s. Thus, it seems reasonable to conclude that the reported CG simulations provide comparable sampling to many milliseconds of equilibrium simulations with atomically detailed models. Aggregate Analysis. Aggregate Definition. Following Wang and Ferguson,76 we employ three intermolecular distances to characterize the interaction between molecules and the organization of aggregates. Figure 3 illustrates these

rmaximin = max[(max(min rab)), (max(min rab))] a∈( b∈)

b∈) a∈(

(6)

Operationally, we compute rmaximin as follows: (1) for each a ∈ ( , define ra→B as the distance from a to the nearest b ∈ ) , (2) define rA→B as the maximum of the resulting set {ra → B}a ∈ ( , (3) for each b ∈ ) , define rb→A as the distance from b to the nearest a ∈ ( , (4) define rB→A as the maximum of the resulting set {rb → A}b ∈ ) , and finally, (5) define rmaximin = max(rA→B,rB→A). Given two bipyrene-type molecules, A and B, we define an analogous distance for each of the four intermolecular core pairs and define rmaximin as the minimum of these four distances. When rmaximin is small, aromatic cores from two molecules are roughly parallel to each other and aligned face-to-face. This metric identifies aggregates with contiguous, aligned stacks of aromatic cores, as expected for the nanoaggregates proposed in the Yen−Mullins model. Conversely, aggregates with staggered or skew cores will not be identified by this metric. Extent of Aggregation. We quantify the extent of aggregation with the time average, G, of the mass-averaged aggregation number g291 G = ⟨g2(t )⟩ =

∑i ni(t )i 2 ∑i ni(t )i

(7)

Here, t indicates the simulated time, ni(t) is the number of aggregates of size i at time t, and the angular brackets denote an average over the sampled configurations from the two simulated replicates. Because each simulation contains 25 molecules, g2 can range from 1, in which case the molecules are all dispersed, to 25, in which case the molecules form a single aggregate. As noted above, different distance metrics identify different aggregates. Accordingly, Gmin, Gcore, and Gmaximin indicate the extent of aggregation according to the minimum distance, minimum core distance, and maximin aggregation criteria, respectively. Intra-Aggregate Alignment. We characterize the structure of aggregates by the second Legendre polynomial of the cosine of the angle formed between asphaltene cores

Figure 3. Illustration of the three distances used for identifying and classifying simulated aggregates. Given the indicated configuration for a pair of asphaltene molecules, the three labeled arrows indicate (a) rmin; (b) rcore; and (c) rmaximin.

distances. In each case, we deem a molecule to be part of an aggregate of a certain class if the corresponding distance, r, from the molecule to any molecule in the aggregate is less than rcut = 1.25r0 = 0.49 nm. Minimum Distance. Given two molecules, A and B, with sites indexed, a ∈ A and b ∈ B, respectively, we define the minimum distance, rmin, as the minimum distance from any a site to any b site. The minimum distance metric identifies disordered aggregates, as well as weakly bound clusters that form via chance encounters between molecules. Minimum Core Distance. Given two ovalene molecules, A and B, with aromatic cores, ( and ) , and aromatic sites that are indexed, a ∈ ( and b ∈ ) , respectively, we define the minimum core distance, rcore, as the minimum distance from any a site to any b site. Given two bipyrene-type molecules, A and B, we define an analogous distance for each of the four intermolecular core pairs and define rcore as the minimum of these four distances. When used in conjunction with rmin, rcore distinguishes between aggregates that are stabilized by core−

⟨P2(cos θ )⟩ =

1 (3 cos2 θ − 1) 2

contact

(8)

Here, θ is the angle formed by the aromatic cores (i.e., the angle formed by the vectors normal to two cores) of two contacting molecules. In the case of ovalene molecules, each pair of contacting molecules corresponds to a single pair of aromatic cores and, thus, a single angle is considered for each pair. In the case of bipyrene molecules, each pair of contacting molecules corresponds to six distinct pairs (i.e., four intermolecular pairs and two intramolecular pairs) of aromatic cores and, thus, six distinct angles are considered for each pair. The subscripted angular brackets in eq 8 indicate an average that weights equally each angle for each pair of contacting molecules in the sampled configurations. The second Legendre polynomial, P2(cos θ), ranges from −0.5, in which case the cores are perpendicular to one another, to 1, in which case the 6114


Article


aliphatic tails, each containing l = 8 L sites. For each circle in Figure 4, we performed two independent simulations of 25 ovalene-8 molecules under the indicated solvent conditions, (εL̅ , εR̅ ), in a volume corresponding to 1% weight fraction asphaltene. As described in the Methods section and Supporting Information, we simulated each replicate for 400 ns or until the simulated aggregates had reached a state that seemed unlikely to change during the remainder of the 400 ns. Thus, in conditions with significant aggregation, the simulations do not necessarily provide converged averages but rather reflect qualitative tendencies for asphaltene aggregation. The color of each circle indicates the extent of aggregation observed in the associated simulations, as quantified by the mass-averaged aggregation number, Gmin, for clusters identified by the minimum distance, rmin, criterion. Dark circles indicate minimal aggregation, brighter circles indicate increasing extent of aggregation, and the red “X” indicates solvent conditions that promote strong, irreversible aggregation. The solid black bars indicate the onset of significant aggregation, that is, Gmin ≥ 2. Asphaltenes remain dispersed under near ideal conditions with εL̅ , εR̅ ≈ 0 and undergo increasing aggregation as the “solvent quality” decreases farther from the origin. However, asphaltenes respond quite differently to variations in εL̅ and εR̅ . In aliphatic solvents (εL̅ = 0), asphaltenes aggregate for εR̅ ≥ 0.07, whereas in aromatic solvents (εR̅ = 0), asphaltenes do not aggregate until εL̅ ≥ 0.21. Thus, relatively weaker aromatic interactions are required to initiate aggregation for the ovalene8 model. This is consistent with the physical intuition of the Yen−Mullins model because the density and rigidity of the aromatic sites in the planar asphaltene core allow for many R− R interactions to simultaneously form with minimal entropic penalty. Conversely, because the aliphatic sites are more flexible and distributed more widely across the molecule, a stronger L−L attraction is necessary to offset the entropic penalty required for simultaneously bringing multiple intermolecular aliphatic pairs into contact. It is interesting that, for the simulated conditions, ovalene-8 molecules form larger aggregates under aromatic conditions than under aliphatic conditions. Presumably, the purely repulsive aliphatic tails retard the formation of larger aggregates in aliphatic conditions.92 Slightly increasing the aliphatic attraction from εL̅ = 0 to 0.05 significantly promotes aggregation even for εR̅ = 0.06 because now attractive interactions between the aliphatic tails also facilitate aggregation. Conversely, starting from aromatic solvents, slightly increasing εR̅ from 0 to 0.01 reduces the threshold for aggregation from εL̅ = 0.21 to 0.18 but does not lead to larger aggregates. Figure 5 presents representative configurations of simulated aggregates. The top, bottom, and middle rows correspond to simulations with aliphatic (εR̅ = 0.11, εL̅ = 0.0), mixed (εR̅ = εL̅ = 0.06), and aromatic (εR̅ = 0, εL̅ = 0.21) (implicit) solvents. The left column presents aggregates formed by the ovalene-8 asphaltene model described in Figures 1a and 4. In aliphatic solvents, ovalene-8 asphaltenes form columns with tightly stacked aromatic cores that closely resemble the prototypical Yen−Mullins nanoaggregates. Similar aggregates also form in mixed solvents, although in this case the aromatic cores appear to stack slightly less tightly. Conversely, in aromatic solvents, ovalene-8 forms disordered aggregates that are stabilized by contacts between aliphatic chains. Although aromatic interactions are purely repulsive in the aromatic solvent, aromatic

cores stack in a parallel fashion. In the case that cores show little tendency to align, ⟨P2(cos θ)⟩ ≈ 0.

■

RESULTS Here, we report the results of simulations with an implicit solvent CG model for asphaltene self-assembly. While the Methods section describes the model in detail, here we briefly summarize its essential features. The model explicitly represents each aromatic carbon with an R site and each aliphatic carbon with an L site. By capturing the rigidity of aromatic cores and flexibility of aliphatic tails, the model accurately describes the shape, conformational flexibility, and packing of representative asphaltene molecules. Conversely, the model adopts a simple mean field treatment for the surrounding petrochemical solvent. We mimic different solvent conditions by varying the dimensionless parameters εL̅ = εL/ kBT and εR̅ = εR/kBT, which quantify the effective attraction between aliphatic and aromatic sites, respectively, without changing the excluded volume of the sites. (Cross interactions between aliphatic and aromatic sites are governed by the geometric mixing rule εRL ̅ = εL̅ εR̅ .) In particular, we mimic aliphatic solvents by defining εL̅ = 0 and εR̅ > 0, such that there is an effective attraction between aromatic, but not aliphatic, sites. Similarly, we mimic aromatic solvents by defining εR̅ = 0 and εL̅ > 0, such that there is an effective attraction between aliphatic, but not aromatic, sites. We define “mixed” solvents by εR̅ , εL̅ > 0, such that both aromatic and aliphatic sites attract. Finally, we define the “ideal” solvent by εR̅ = εL̅ = 0. Given this very simple and efficient model, we explored the impact of solvent conditions upon the self-assembly of model asphaltenes. Figure 4 presents the resulting “solvent phase diagram” for the ovalene-8 molecule, which is representative of the island, or continental, asphaltene architecture. As indicated by Figure 1a, this molecule consists of a single ovalene aromatic core with 32 R sites that are surrounded by four

Figure 4. Solvent phase diagram quantifying the impact of solvent quality upon the self-assembly of the ovalene-8 asphaltene model. Each circle indicates the results of two independent simulations with 25 identical ovalene-8 molecules under the specified solvent conditions, i.e., εL̅ = εL/kBT and εR̅ = εR/kBT. The circle colors indicate the simulated mass-averaged aggregation number, Gmin, which employs the min criterion for identifying aggregates. Red X’s indicate solvent conditions that resulted in irreversible aggregation. The black bars indicate the onset of aggregation (i.e., Gmin ≥ 2) in aliphatic solvents along the y-axis and in aromatic conditions along the x-axis. 6115


Article


Figure 6. (a−f) Probability, PN, for observing aggregates of size N under the simulated conditions of Figure 5. The left and right columns present results for ovalene-8 and bipyrene-8 molecules, respectively. The top, middle, and bottom rows correspond to simulations in aliphatic solvents (εR̅ = 0.11, εL̅ = 0), mixed solvents (εR̅ = εL̅ = 0.06), and aromatic solvents (εR̅ = 0, εL̅ = 0.21), respectively. The solid black, dashed red, and dotted blue curves employ the min, core, and maximin criteria, respectively, for identifying aggregates.

Figure 5. (a−f) Representative aggregates observed under selected solvent conditions. The left and right columns present simulated aggregates for ovalene-8 and bipyrene-8 molecules, respectively. The top, middle, and bottom rows correspond to simulations in aliphatic solvents (εR̅ = 0.11, εL̅ = 0), mixed solvents (εR̅ = εL̅ = 0.06), and aromatic solvents (εR̅ = 0, εL̅ = 0.21), respectively.

cores occasionally contact and even align as a consequence of packing considerations. The right column presents representative configurations of aggregates formed in simulations of the bipyrene-8 asphaltene model under corresponding solvent conditions. As illustrated in Figure 1b, bipyrene-8 illustrates the archipelago architecture with 32 aromatic R sites equally divided between two smaller cores that are connected by a short aliphatic linker and surrounded by 4 aliphatic chains, each containing l = 8 L sites. In comparison to ovalene-8, bipyrene-8 forms considerably smaller aggregates in these aliphatic and mixed solvent conditions. Moreover, the bipyrene-8 aggregates are more disordered and demonstrate fewer stacking interactions. Apparently, the flexibility provided by a short aliphatic linker introduces a significant entropic penalty that dramatically reduces stacking interactions even in the presence of strong R−R attraction. In aromatic solvents, bipyrene-8 asphaltenes form large disordered aggregates that are stabilized by L−L attractions among densely packed aliphatic tails with relatively few R−R contacts. Interestingly, under these conditions, the asphaltenes appear to expel the aromatic cores to the surface of the aggregate as if beginning to form a small micelle. Figure 6 quantitatively analyzes the distribution of aggregates formed in these simulations. As in Figure 5, the left and right panels present results for ovalene-8 and bipyrene8 asphaltenes, respectively, while the top, middle, and bottom rows present results for aliphatic (εR̅ = 0.11, εL̅ = 0.0), mixed (εR̅ = εL̅ = 0.06), and aromatic (εR̅ = 0, εL̅ = 0.21) solvent conditions, respectively. Each panel presents the probability, PN, for an asphaltene to be found in an aggregate of size N. 1 n Specifically, we define PN = n n ∑t =t 1 nN , t , where the sum is

reported distributions are averaged over two independent simulations that demonstrate qualitative agreement, as illustrated by Figures S1 and S2 of the Supporting Information. As detailed in the Methods section, different distances can be computed between molecules, and these different distances lead to different criteria for identifying clusters. The solid black lines correspond to the least restrictive “min” distance criterion, which associates two molecules in a cluster whenever any two sites contact. The dashed red lines correspond to the “core” distance criterion, which associates two molecules in a cluster whenever their aromatic R sites contact. The dotted blue lines correspond to the most restrictive “maximin” distance criterion, which associates two molecules in a cluster whenever their entire aromatic cores contact. In the case of bipyrene asphaltenes, the core and maximin distances are applied to the closest intermolecular pair of cores, that is, two bipyrene molecules are deemed to be in contact when any of their cores contact. Clusters that are identified by the min criterion, but not by the core criterion, are stabilized by aliphatic contacts. Clusters that are identified by the min and core criteria, but not by the maximin criterion, are stabilized by a combination of aliphatic and aromatic contacts but do not demonstrate aromatic stacking. Clusters that are identified by all three criteria demonstrate tight core stacking as in Figure 5a. Figure 6a,c demonstrates that ovalene-8 asphaltenes demonstrate similar aggregation in aliphatic and mixed solvents. In both solvent conditions, ovalene-8 molecules aggregate via tight stacking interactions and monomers are only rarely observed. The attraction between aliphatic sites in mixed solvents facilitates the formation of larger clusters. However, the jagged distributions suggest that both simulations may be sampling slowly evolving distributions that might possibly evolve to a single aggregate containing all

t mol

evaluated over the nt configurations sampled from the two replicate simulations, nmol = 25 is the number of simulated molecules, and nN,t is the number of molecules in aggregates of size N in sampled configuration t. We emphasize that the 6116


Article


ovalene-8 asphaltenes and is identical to Figure 4, while Figure 7d presents results for the bipyrene-8 asphaltenes that are described by the right columns of Figures 5 and 6. We first consider self-assembly due to aromatic interactions in aliphatic solvents (εL̅ = 0, εR̅ > 0). Interestingly, the onset of aggregation appears relatively insensitive to chain length. For tails of length l = 5, 8, or 11, the critical threshold for aggregation in aliphatic solvents, ε*R̅ 0, is 0.07 for ovalene-type asphaltenes and 0.12−0.13 for bipyrene-type asphaltenes. As expected, because L sites are purely repulsive under these conditions, the average size of aggregates decreases with increasing tail length. However, it is somewhat surprising that, in aliphatic solvents with εL̅ = 0, ovalene- and bipyrene-type asphaltenes respond differently to increasing εR̅ . In particular, for the considered range, εR*̅ 0 ≤ εR̅ ≤ 0.16, ovalene aggregates do not grow noticeably larger with increasing εR̅ . In contrast, over the similar range ε*R̅ 0 ≤ εR̅ ≤ 0.17, bipyrene aggregates grow increasingly large with increasing εR̅ . Thus, under aliphatic conditions, dimers and small aggregates appear to form more easily for island asphaltenes than for archipelago asphaltenes, but archipelago asphaltenes aggregate to a greater extent with increasing εR̅ . Under these conditions, the aliphatic tails appear to limit the aggregation of island asphaltenes by sterically blocking “molecular recognition sites,”92 while the greater flexibility of archipelago cores may enable the formation of larger aggregates under conditions that promote sufficiently strong aromatic interactions. Next, we consider asphaltene self-assembly driven by aliphatic interactions in aromatic solvents (εR̅ = 0, εL̅ > 0). Asphaltenes with aliphatic tails of length l = 5 do not associate over the range of εL̅ considered. For asphaltenes with aliphatic tails of length l = 8 and 11, the critical threshold for aggregation in aromatic solvents ε*L̅ 0 decreases to 0.19 and 0.15, respectively. Thus, as expected, in aromatic solvents that stabilize aliphatic contacts, aggregation occurs more readily with increasing l. Interestingly, εL*̅ 0 is the same for both ovalene and bipyrene asphaltenes with the same l, which suggests that aliphatic aggregation is rather insensitive to the architecture of aromatic cores. Moreover, in striking contrast to aggregation in aliphatic solvents, very large aggregates form at the critical threshold εL̅ ≈ εL*̅ 0 in aromatic solvents. Even more strikingly, for asphaltenes with relatively long tails l = 11, the observed aggregates shrink as εL̅ increases further. This nonmonotonic behavior perhaps reflects the competition between intra- and intermolecular L−L interactions because the formation of intramolecular contacts would reduce the aliphatic surface available for forming intermolecular networks. However, it is also possible that this nonmonotonic trend reflects kinetic trapping that arises for smaller aggregates as the aliphatic attraction increases. Finally, Figure 7 demonstrates that asphaltenes generally respond to mixed solvents with εR̅ , εL̅ > 0 in a manner similar to that observed for ovalene-8 in Figure 4. As intuitively expected, decreasing solvent quality generally promotes aggregation. In particular, starting from aliphatic conditions with εL̅ = 0 and increasing εL̅ to 0.05 increases the stability and size of asphaltene aggregates. Similarly, starting from aromatic conditions with εR̅ = 0 and increasing εR̅ to 0.02 also increases the stability and size of asphaltene aggregates. Interestingly, these changes in solvent quality appear to have a greater stabilizing effect on bipyrene aggregates than ovalene aggregates.

25 simulated asphaltenes in the limit of infinite simulation time. Figure 6b,d demonstrates that bipyrene-8 asphaltenes also behave similarly in aliphatic and mixed solvents. Under these solvent conditions, bipyrene-8 asphaltenes remain almost entirely dispersed. The small clusters that are observed do form aromatic contacts between R sites but rarely demonstrate core stacking. Finally, Figure 6e,f demonstrates that ovalene-8 and bipyrene-8 asphaltenes behave similarly in aromatic solvents. Under these solvent conditions, neither asphaltene model forms stacked aggregates. Instead, they tend to form a single very large aggregate that is stabilized by aliphatic interactions with occasional chance contacts between aromatic sites, as suggested in Figure 5e,f. While Figures 4−6 analyze the importance of solvent quality and core architecture, Figure 7 investigates the importance of

Figure 7. (a−f) Impact of core architecture and tail length upon asphaltene self-assembly. Each panel presents a simulated solvent phase diagram for model asphaltenes as a function of εL̅ and εR̅ . The left and right columns present results for asphaltenes with ovalene and bipyrene cores, respectively. The top, middle, and bottom rows present results for asphaltenes with tails of length l = 5, 8, and 11, respectively. The circle colors indicate the simulated mass-averaged aggregation number, Gmin; red X’s indicate irreversible aggregation; and black bars indicate the onset of aggregation (i.e., Gmin ≥ 2).

aliphatic tails for asphaltene aggregation. The panels of Figure 7 present solvent phase diagrams analogous to Figure 4 for model asphaltene molecules with aliphatic tails of varying lengths. As in Figure 4, the circles indicate simulations with the specified solvent conditions, the colors indicate the massaveraged aggregation number Gmin averaged over two independent replicate simulations, the solid black lines indicate the onset of aggregation (i.e., Gmin ≥ 2), and the red X’s indicate the formation of large, kinetically trapped aggregates. The left and right columns present results for model asphaltenes with ovalene and bipyrene cores, respectively. The top, middle, and bottom rows present results for model asphaltenes with aliphatic tails of length l = 5, 8, and 11 L sites, respectively. In particular, Figure 7c presents results for 6117


Article

The Journal of Physical Chemistry B Figure 8 characterizes the structure of the aggregates described by Figure 7. As in Figure 7, the left and right

Figure 9. (a−f) Formation of Yen−Mullins nanoaggregates. Each panel presents a simulated solvent phase diagram for model asphaltenes as a function of εL̅ and εR̅ . The left and right columns present results for asphaltenes with ovalene and bipyrene cores, respectively. The top, middle, and bottom rows present results for asphaltenes with tails of length l = 5, 8, and 11, respectively. The circle colors indicate the simulated mass-averaged aggregation number, Gmaximin; red X’s indicate irreversible aggregation; and black bars indicate the onset of aggregation (i.e., Gmin ≥ 2).

Figure 8. (a−f) Alignment of aromatic cores within simulated aggregates. As in Figure 7, the left and right columns present results for asphaltenes with ovalene and bipyrene cores, respectively, while the top, middle, and bottom rows present results for asphaltenes with tails of length l = 5, 8, and 11, respectively. The circle colors indicate ⟨P2(cos θ)⟩core, where θ indicates the angle formed between pairs of aromatic cores within a single aggregate, and the average is performed over all simulated aggregates that are identified by the “core” metric for the specified system. The black bars indicate the onset of aggregation as indicated by Gmin ≥ 2.

However, as noted above, under these conditions, the nanoaggregates quickly shrink with increasing chain length and are very small for l = 11. Conversely, these aggregates are slightly stabilized and considerably larger when εL̅ increases to 0.05. Indeed, we only observe stacked asphaltene nanoaggregates that are larger than 7−10 molecules in such mixed solvents, which suggests that attractive aliphatic interactions may play an important role in the formation of larger asphaltene aggregates. We note, though, that the small aggregates may also be a consequence of the relatively small number of simulated asphaltenes and the finite, though relatively long, simulations.

columns present results for ovalene and bipyrene cores, respectively, while the top, middle, and bottom rows present results for aliphatic tails of length l = 5, 8, and 11. In Figure 8, the circle colors indicate the metric ⟨P2(cos θ)⟩core, which quantifies the alignment of the aromatic cores within the simulated aggregates identified by the core distance metric, that is, the aggregates identified by networks of aromatic R−R contacts. Clearly, ovalene aggregates form highly ordered structures under conditions that favor aromatic contacts. In fact, ovalene aggregates demonstrate strongly aligned aromatic cores whenever εR̅ > 0.05, even under conditions that do not promote significant aggregation. Conversely, bipyrene molecules demonstrate relatively little tendency to align their aromatic cores, although this alignment slightly increases with increasing εR̅ . Finally, Figure 9 identifies the conditions that favor the formation of the nanoaggregates proposed by the Yen−Mullins model. As in Figures 7 and 8, the left and right columns present results from simulations of ovalene- and bipyrene-type asphaltenes, respectively, while the top, middle, and bottom rows present results for asphaltenes with aliphatic tails of length l = 5, 8, and 11. In Figure 9, the circle colors indicate the mass-averaged size of aggregates defined by the most restrictive maximin criterion, which only counts aggregates formed by tight overlap between aromatic cores. Bipyrene-type asphaltenes do not form significant Yen−Mullins nanoaggregates under any of the simulated conditions. Ovalenetype asphaltenes form prototypical aggregates in aliphatic solvents promoting aromatic attractions (εL̅ = 0, εR̅ ≥ 0.07).

■

DISCUSSION In this work, we have introduced and simulated a very simple CG model for asphaltene self-assembly in petrochemical mixtures. The model represents each asphaltene molecule in united atom detail, while distinguishing between aliphatic and aromatic atoms. Conversely, the model implicitly represents the surrounding petrochemical environment and mimics different solvent conditions by varying the effective attraction between aliphatic and aromatic sites. Specifically, in mixed solvents both aromatic and aliphatic sites attract, in aliphatic solvents only aromatic sites attract, and in aromatic solvents only aliphatic sites attract. We anticipate that this model may prove particularly useful for studying generic aspects of asphaltene aggregation. The model reasonably describes the regular solution considerations that may drive asphaltene aggregation, that is, the microphase separation of aromatic groups from a poor aliphatic solvent.41−45 The model also realistically describes the shape, flexibility, conformational fluctuations, and packing of 6118


Article


We also explored the self-assembly of asphaltenes under aromatic conditions in which only aliphatic sites attract. For the asphaltenes considered in this study, a relatively strong aliphatic attraction is necessary for the simulated asphaltenes to aggregate under these conditions. The threshold for aggregation in aromatic solvents appears insensitive to the architecture of the aromatic cores and, as expected, systematically decreases with increasing aliphatic chain length. Under these conditions, both ovalene and bipyrene asphaltenes form large disordered aggregates that are stabilized by a network of aliphatic contacts. While aromatic cores are present in the center of the ovalene aggregates, the bipyrene aggregates appear to expel their cores to the surface, as if forming a micelle. It should be noted, of course, that asphaltenes are defined by their solubility in toluene. Consequently, the aggregates observed in simulated aromatic conditions either require a solvent that is “more aromatic” than toluene or they correspond to “pseudo-asphaltenes” that have longer aliphatic tails than the molecules that are defined by the asphaltene fraction of crude oil. Finally, it is important to emphasize the significant limitations of the present model. Most importantly, the model is fundamentally limited by its simple mean field description of the solvent, which represents a very crude approximation to the many-body, state-point-dependent potential of mean force that exactly describes the effective interaction between solutes by tracing over solvent coordinates.86,87 Consequently, the model is unlikely to accurately describe aggregates that are stabilized by direct interactions with solvent molecules. Similarly, this model provides limited insight into the role of solvent for driving asphaltene aggregation or for the temperature-dependence of nanoaggregation.49 Additionally, the elimination of the solvent precludes a realistic description of fast dynamical processes that involve solvent molecules on short timescales, although the model may provide a reasonable description of dynamics occurring on sufficiently long timescales.70 Moreover, this implicit solvent approximation hinders, if not completely precludes, mapping the simulated conditions specified by (εR̅ , εL̅ ) onto any particular real solvent. In particular, the simulated “purely” aromatic (εR̅ = 0) and aliphatic (εL̅ = 0) conditions correspond to rather “extreme” solvents in which the effective interactions between aromatic and aliphatic groups, respectively, are purely repulsive. Quite generally, the model should not be expected to provide quantitatively accurate predictions but rather only qualitative insight into generic aspects of asphaltene self-assembly.80−82 In addition to these fundamental limitations of the implicit solvent model, the present simulation study is obviously quite limited in scope. In the spirit of this mean field approximation, we considered model molecules with particularly simple and symmetric aromatic cores and aliphatic tails. In particular, the symmetry of the simulated aliphatic tails likely precludes the possibility of t-shaped interfaces between the cores and discourages off-center core stacking that might lead to more complex aggregate morphologies. Consequently, while prior simulation studies have reported aggregates with branched,63 curved,96 helical,75 and t-shaped77,97 morphologies, we observed only rather simple aggregates. Furthermore, our simulated model lacked strong, directional attractions that might result from the presence of heteroatoms.20,26 Moreover, in contrast to the polydiversity present in crude oil, our simulations considered only monodisperse asphaltene solu-

asphaltene molecules, which are not easily described by analytic mean field theories. Moreover, the implicit solvent model provides sufficient efficiency for rapidly exploring a wide range of solvent conditions, asphaltene geometries, and solution mixtures. This efficiency is particularly critical for modeling the dilute conditions in which nanoaggregates first form. Given this efficiency, we examined the phase behavior of 147 different asphaltene solutions over a cumulative simulation time of 14.7 μs. Furthermore, given the accelerated sampling provided by the smooth energy surface of the CG model, it is reasonable to estimate that these simulations provide similar sampling to many milliseconds of simulation time with an explicit solvent model. We simulated solutes with an ovalene core surrounded by four aliphatic tails as representative examples of the island model for asphaltenes. In aliphatic solvents, that is, under conditions in which only aromatic sites attract, the ovalene asphaltenes form prototypical nanoaggregates with stacked aromatic cores in accord with the Yen−Mullins model.23,39−41 Interestingly, the onset of aggregation appears to be independent of the length of the aliphatic tails that surround the aromatic core. This appears consistent with recent mass spectrometry experiments that observed similar aggregation numbers in very diverse asphaltene samples.93 However, the simulated nanoaggregates do not appear to grow larger with increasing attraction between aromatic sites. For the simulated aliphatic conditions, these nanoaggregates contain only 7−10 molecules. Moreover, the nanoaggregates shrink as the aliphatic tails become longer. Thus, in aliphatic solvents, the corona of surrounding aliphatic tails appears to stunt the growth of nanoaggregates via steric interactions, as suggested by prior simulation studies.92,94 Interestingly, stacked nanoaggregates grow largest and most readily in mixed solvents that promote both aromatic and aliphatic interactions, as may be expected in crude oil mixtures. Indeed, prior simulations have observed Yen−Mullins nanoaggregates in mixed heptane− toluene solvents.58,94 We also simulated solutes with a bipyrene core surrounded by four aliphatic tails, as representative examples of the archipelago model for asphaltenes. These asphaltenes also aggregate under aliphatic conditions, although bipyrene aggregation requires considerably stronger aromatic interactions due to the flexibility of the aliphatic linker connecting their smaller aromatic cores. As observed for ovalene asphaltenes, the threshold for bipyrene aggregation in aliphatic solvents appears insensitive to the length of their aliphatic tails. Interestingly, in contrast to ovalene aggregates, bipyrene aggregates become increasingly large with increasing aromatic attraction. However, bipyrene asphaltenes do not form prototypical nanoaggregates under any of the simulated conditions.34,36 These results complement recent simulations of archipelago-type asphaltenes by Wang et al., who observed the formation of large aggregates that were stabilized by contacts between the aromatic cores.95 Importantly, however, Wang et al. simulated more concentrated asphaltene solutions. Moreover, they simulated much larger archipelago-type asphaltenes with larger aromatic cores connected by longer aliphatic linkers.95 It would be interesting to examine whether our smaller bipyrene asphaltene models also form extended networks at these higher concentrations and to investigate the sensitivity of such networks to the size of the aromatic cores and aliphatic linkers. 6119


The Journal of Physical Chemistry B tions at a single asphaltene concentration. While acknowledging its significant limitations, we anticipate that the present CG model may be usefully adapted to efficiently examine many of these more complex considerations.

■

CONCLUSIONS

ACKNOWLEDGMENTS

■

REFERENCES

(1) Kokal, S. L.; Sayegh, S. G. Asphaltenes: the cholesterol of petroleum. Middle East Oil Show, SPE Paper No. 29787, 1995; pp 169−181. (2) Creek, J. L. Freedom of action in the state of asphaltenes: escape from conventional wisdom. Energy Fuels 2005, 19, 1212−1224. (3) Sheu, E. Y.; Mullins, O. C. Asphaltenes: Fundamentals and Applications; Plenum Press, 1995. (4) Spiecker, P. M.; Gawrys, K. L.; Kilpatrick, P. K. Aggregation and solubility behavior of asphaltenes and their subfractions. J. Colloid Interface Sci. 2003, 267, 178−193. (5) Pina, A.; Mougin, P.; Behar, E. Characterisation of asphaltenes and modelling of flocculation - state of the art. Oil Gas Sci. Technol. 2006, 61, 319−343. (6) Nabzar, L.; Aguiléra, M. E. The colloidal approach. a promising route for asphaltene deposition modelling. Oil Gas Sci. Technol. 2008, 63, 21−35. (7) Sheu, E. Y. Petroleum asphaltenes - properties, characterization, and issues. Energy Fuels 2002, 16, 74−82. (8) Rogel, E.; Ovalles, C.; Moir, M. Asphaltene stability in crude oils and petroleum materials by solubility profile analysis. Energy Fuels 2010, 24, 4369−4374. (9) Buenrostro-Gonzalez, E.; Lira-Galeana, C.; Gil-Villegas, A.; Wu, J. Asphaltene precipitation in crude oils: theory and experiments. AIChE J. 2004, 50, 2552−2570. (10) Groenzin, H.; Mullins, O. C. Molecular size and structure of asphaltenes from various sources. Energy Fuels 2000, 14, 677−684. (11) Wargadalam, V.; Norinaga, K.; Iino, M. Size and shape of a coal asphaltene studied by viscosity and diffusion coefficient measurements. Fuel 2002, 81, 1403−1407. (12) Strausz, O. P.; Peng, P.; Murgich, J. About the colloidal nature of asphaltenes and the mw of covalent monomeric units. Energy Fuels 2002, 16, 809−822. (13) Behrouzi, M.; Luckham, P. F. Limitations of size-exclusion chromatography in analyzing petroleum asphaltenes: a proof by atomic force microscopy. Energy Fuels 2008, 22, 1792−1798. (14) Herod, A. A.; Kandiyoti, R. Comment on ”Limitations of sizeexclusion chromatography in analyzing petroleum asphaltenes: a proof by atomic force microscopy” by M. Behrouzi and P. F. Luckham, Energy & Fuels, 2008, 22 (3), 1792-17989 doi: 10.1021/ ef800064q. Energy Fuels 2008, 22, 4307−4309. (15) Strausz, O. P.; Safarik, I.; Lown, E. M.; Morales-Izquierdo, A. A critique of asphaltene fluorescence decay and depolarization-based claims about molecular weight and molecular architecture. Energy Fuels 2008, 22, 1156−1166. (16) Mullins, O. C. Rebuttal to Strausz et al. regarding time-resolved fluorescence depolarization of asphaltenes. Energy Fuels 2009, 23, 2845−2854. (17) Mullins, O. C.; Martínez-Haya, B.; Marshall, A. G. Contrasting perspective on asphaltene molecular weight. this comment vs the

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.9b04275. Additional analysis of the difference between replicate runs of the solvent conditions shown in Figure 6 and phase diagram of the extent of aggregation by Gcore for all of the model asphaltenes examined here for comparison with Figures 7−9 (PDF)

■

■

Acknowledgment is made to the Donors of the American Chemical Society Petroleum Research Fund for the support (or the partial support) of this research. W.G.N. gratefully acknowledges the financial support of ACS PRF award 52100ND6 and REU grant CHE-1263053. Figures 1 and 5 employed VMD.98 VMD is developed with the NIH support by the Theoretical and Computational Biophysics group at the Beckman Institute, University of Illinois at Urbana-Champaign. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562.99 The computational resources required for this study were also provided in part by the Penn State Institute for Cyberscience and the Minnesota Supercomputing Institute.

We have developed a united atom, implicit solvent CG model for investigating the generic consequences of regular solution considerations for the self-assembly of asphaltenes in petrochemical mixtures. As expected from the Yen−Mullins model, our simulations indicate that island-type asphaltenes form nanoaggregates with stacked aromatic cores under conditions that promote aromatic attraction. Interestingly, the threshold for nanoaggregation appears insensitive to the length of the aliphatic tails that surround the aromatic core. However, these tails stunt the further growth of nanoaggregates under “ideal aliphatic” conditions for which there is no effective attraction between aliphatic sites. Consequently, the simulated nanoaggregates are most stable and grow most readily in mixed solvents that promote an effective attraction between both aliphatic and aromatic sites. Conversely, archipelago asphaltenes do not form prototypical nanoaggregates under any simulated conditions, although they do form large aggregates under conditions that promote strong attraction between either aromatic or aliphatic sites. Our simulation results generally agree with the expectations of the Yen−Mullins model, as well as prior simulation studies of nanoaggregate formation. The present study also motivates many potential future directions. Notwithstanding its limitations, we anticipate that the present model may prove useful for understanding, among other considerations, molecular asymmetry within asphaltenes, the effects of heteroatoms within asphaltene cores and tails, and polydiversity within asphaltene solutions. Moreover, we anticipate that the efficiency of this model may enable systematic simulation studies with larger systems that probe the later stages of asphaltene self-assembly under dilute conditions.

■

Article

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

W. G. Noid: 0000-0001-9675-8489 Present Address †

Minnesota Supercomputing Institute, University of Minnesota, Minneapolis, MN 55455. Notes

The authors declare no competing financial interest. 6120


Article

The Journal of Physical Chemistry B overview of A. A. Herod, K. D. Bartle, and R. Kandiyoti. Energy Fuels 2008, 22, 1765−1773. (18) Speight, J. G. The Chemistry and Technology of Petroleum, 4th ed.; CRC Press, 2007. (19) Hoepfner, M. P.; Limsakoune, V.; Chuenmeechao, V.; Maqbool, T.; Fogler, H. S. A fundamental study of asphaltene deposition. Energy Fuels 2013, 27, 725−735. (20) Li, D. D.; Greenfield, M. L. Chemical compositions of improved model asphalt systems for molecular simulations. Fuel 2014, 115, 347−356. (21) Wiehe, I.; Liang, K. Asphaltenes, resins, and other petroleum macromolecules. Fluid Phase Equilib. 1996, 117, 201−210 , Proceedings of the Seventh International Conference on Fluid Properties and Phase Equilibria for Chemical Process Design . (22) Hansen, J. S.; Lemarchand, C. A.; Nielsen, E.; Dyre, J. C.; Schroder, T. Four-component united-atom model of bitumen. J. Chem. Phys. 2013, 138, 094508. (23) Mullins, O. C.; Sabbah, H.; Eyssautier, J.; Pomerantz, A. E.; Barré, L.; Andrews, A. B.; Ruiz-Morales, Y.; Mostowfi, F.; McFarlane, R.; Goual, L.; et al. Advances in asphaltene science and the YenMullins model. Energy Fuels 2012, 26, 3986−4003. (24) Behrouzi, M.; Luckham, P. F. Limitations of size-exclusion chromatography in analyzing petroleum asphaltenes: a proof by atomic force microscopy. Energy Fuels 2008, 22, 1792−1798. (25) Murgich, J. Molecular simulation and the aggregation of the heavy fractions in crude oils. Mol. Simul. 2003, 29, 451−461. (26) Boek, E. S.; Yakovlev, D. S.; Headen, T. F. Quantitative molecular representation of asphaltenes and molecular dynamics simulation of their aggregation. Energy Fuels 2009, 23, 1209−1219. (27) Barré, L.; Jestin, J.; Morisset, A.; Palermo, T.; Simon, S. Relation between nanoscale structure of asphaltene aggregates and their macroscopic solution properties. Oil Gas Sci. Technol. 2009, 64, 617−628. (28) Herod, A. A.; Bartle, K. D.; Kandiyoti, R. Characterization of heavy hydrocarbons by chromatographic and mass spectrometric methods: an overview. Energy Fuels 2007, 21, 2176−2203. (29) Liao, Z.; Zhao, J.; Creux, P.; Yang, C. Discussion on the structural features of asphaltene molecules. Energy Fuels 2009, 23, 6272−6274. (30) Pomerantz, A. E.; Hammond, M. R.; Morrow, A. L.; Mullins, O. C.; Zare, R. N. Asphaltene molecular-mass distribution determined by two-step laser mass spectrometry. Energy Fuels 2009, 23, 1162−1168. (31) Klee, T.; Masterson, T.; Miller, B.; Barrasso, E.; Bell, J.; Lepkowicz, R.; West, J.; Haley, J. E.; Schmitt, D. L.; Flikkema, J. L.; et al. Triplet electronic spin states of crude oils and asphaltenes. Energy Fuels 2011, 25, 2065−2075. (32) Andrews, A. B.; Edwards, J. C.; Pomerantz, A. E.; Mullins, O. C.; Nordlund, D.; Norinaga, K. Comparison of coal-derived and petroleum asphaltenes by 13C nuclear magnetic resonance, DEPT, and XRS. Energy Fuels 2011, 25, 3068−3076. (33) Sabbah, H.; Morrow, A. L.; Pomerantz, A. E.; Zare, R. N. Evidence for island structures as the dominant architecture of asphaltenes. Energy Fuels 2011, 25, 1597−1604. (34) Schuler, B.; Meyer, G.; Peña, D.; Mullins, O. C.; Gross, L. Unraveling the molecular structures of asphaltenes by atomic force microscopy. J. Am. Chem. Soc. 2015, 137, 9870−9876. (35) Schuler, B.; Zhang, Y.; Collazos, S.; Fatayer, S.; Meyer, G.; Pérez, D.; Guitián, E.; Harper, M. R.; Kushnerick, J. D.; Peña, D.; et al. Characterizing aliphatic moieties in hydrocarbons with atomic force microscopy. Chem. Sci. 2017, 8, 2315−2320. (36) Schuler, B.; Fatayer, S.; Meyer, G.; Rogel, E.; Moir, M.; Zhang, Y.; Harper, M. R.; Pomerantz, A. E.; Bake, K. D.; Witt, M.; et al. Heavy oil based mixtures of different origins and treatments studied by atomic force microscopy. Energy Fuels 2017, 31, 6856−6861. (37) Chacón-Patiño, M. L.; Rowland, S. M.; Rodgers, R. P. Advances in asphaltene petroleomics. part 1: asphaltenes are composed of abundant island and archipelago structural motifs. Energy Fuels 2017, 31, 13509−13518.

(38) Chacón-Patiño, M. L.; Rowland, S. M.; Rodgers, R. P. Advances in asphaltene petroleomics. part 2: selective separation method that reveals fractions enriched in island and archipelago structural motifs by mass spectrometry. Energy Fuels 2018, 32, 314−328. (39) Yen, T. F.; Erdman, J. G.; Pollack, S. S. Investigation of structure of petroleum asphaltenes by X-ray diffraction. Anal. Chem. 1961, 33, 1587−1594. (40) Dickie, J. P.; Yen, T. F. Macrostructures of the asphaltic fractions by various instrumental methods. Anal. Chem. 1967, 39, 1847−1852. (41) Mullins, O. C. The modified Yen model. Energy Fuels 2010, 24, 2179−2207. (42) Rogel, E. Asphaltene aggregation: a molecular thermodynamic approach. Langmuir 2002, 18, 1928−1937. (43) Rogel, E. Thermodynamic modeling of asphaltene aggregation. Langmuir 2004, 20, 1003−1012. (44) Sirota, E. B. Physical structure of asphaltenes. Energy Fuels 2005, 19, 1290−1296. (45) Sirota, E. B.; Lin, M. Y. Physical behavior of asphaltenes. Energy Fuels 2007, 21, 2809−2815. (46) Sheu, E.; Long, Y.; Hamza, H. In Asphaltenes, Heavy Oils, and Petroleomics; Mullins, O. C., Sheu, E. Y., Hammami, A., Marshall, A. G., Eds.; Springer New York: New York, NY, 2007; pp 259−277. (47) Zeng, H.; Song, Y.-Q.; Johnson, D. L.; Mullins, O. C. Critical nanoaggregate concentration of asphaltenes by direct-current (DC) electrical conductivity. Energy Fuels 2009, 23, 1201−1208. (48) Goual, L. Impedance spectroscopy of petroleum fluids at low frequency. Energy Fuels 2009, 23, 2090−2094. (49) Lisitza, N. V.; Freed, D. E.; Sen, P. N.; Song, Y.-Q. Study of asphaltene nanoaggregation by nuclear magnetic resonance (NMR). Energy Fuels 2009, 23, 1189−1193. (50) Indo, K.; Ratulowski, J.; Dindoruk, B.; Gao, J.; Zuo, J.; Mullins, O. C. Asphaltene nanoaggregates measured in a live crude oil by centrifugation. Energy Fuels 2009, 23, 4460−4469. (51) Mostowfi, F.; Indo, K.; Mullins, O. C.; McFarlane, R. Asphaltene nanoaggregates studied by centrifugation. Energy Fuels 2009, 23, 1194−1200. (52) Goual, L.; Sedghi, M.; Zeng, H.; Mostowfi, F.; McFarlane, R.; Mullins, O. C. On the formation and properties of asphaltene nanoaggregates and clusters by DC-conductivity and centrifugation. Fuel 2011, 90, 2480−2490. (53) Anisimov, M. A.; Yudin, I. K.; Nikitin, V.; Nikolaenko, G.; Chernoutsan, A.; Toulhoat, H.; Frot, D.; Briolant, Y. Asphaltene aggregation in hydrocarbon solutions studied by photon correlation spectroscopy. J. Phys. Chem. 1995, 99, 9576−9580. (54) Yudin, I. K.; Anisimov, M. A. In Asphaltenes, Heavy Oils, and Petroleomics; Mullins, O. C., Sheu, E. Y., Hammami, A., Marshall, A. G., Eds.; Springer New York: New York, NY, 2007; pp 439−468. (55) Eyssautier, J.; Frot, D.; Barré, L. Structure and dynamic properties of colloidal asphaltene aggregates. Langmuir 2012, 28, 11997−12004. (56) Tanaka, R.; Sato, E.; Hunt, J. E.; Winans, R. E.; Sato, S.; Takanohashi, T. Characterization of asphaltene aggregates using X-ray diffraction and small-angle X-ray scattering. Energy Fuels 2004, 18, 1118−1125. (57) Ortega-Rodríguez, A.; Cruz, S. A.; Gil-Villegas, A.; GuevaraRodríguez, F.; Lira-Galeana, C. Molecular view of the asphaltene aggregation behavior in asphaltene-resin mixtures. Energy Fuels 2003, 17, 1100−1108. (58) Headen, T. F.; Boek, E. S.; Skipper, N. T. Evidence for asphaltene nanoaggregation in toluene and heptane from molecular dynamics simulations. Energy Fuels 2009, 23, 1220−1229. (59) Kuznicki, T.; Masliyah, J. H.; Bhattacharjee, S. Aggregation and partitioning of model asphaltenes at toluene-water interfaces: molecular dynamics simulations. Energy Fuels 2009, 23, 5027−5035. (60) Teklebrhan, R. B.; Ge, L.; Bhattacharjee, S.; Xu, Z.; Sjöblom, J. Initial Partition and Aggregation of Uncharged Polyaromatic Molecules at the Oil−Water Interface: A Molecular Dynamics Simulation Study. J. Phys. Chem. B 2014, 118, 1040−1051. 6121


Article

The Journal of Physical Chemistry B (61) Teklebrhan, R. B.; Ge, L.; Bhattacharjee, S.; Xu, Z.; Sjöblom, J. Probing structure-nanoaggregation relations of polyaromatic surfactants: a molecular dynamics simulation and dynamic light scattering study. J. Phys. Chem. B 2012, 116, 5907−5918. (62) Liu, J.; Zhao, Y.; Ren, S. Molecular dynamics simulation of selfaggregation of asphaltenes at an oil/water interface: formation and destruction of the asphaltene protective film. Energy Fuels 2015, 29, 1233−1242. (63) Lemarchand, C. A.; Hansen, J. S. Simple statistical model for branched aggregates: application to cooee bitumen. J. Phys. Chem. B 2015, 119, 14323−14331. (64) Lemarchand, C. A.; Greenfield, M. L.; Hansen, J. S. Dynamics and Structure of Bitumen-Water Mixtures. J. Phys. Chem. B 2016, 120, 5470−5480. (65) Teklebrhan, R. B.; Jian, C.; Choi, P.; Xu, Z.; Sjöblom, J. Competitive adsorption of naphthenic acids and polyaromatic molecules at a toluene-water interface. J. Phys. Chem. B 2016, 120, 12901−12910. (66) Headen, T. F.; Boek, E. S.; Jackson, G.; Totton, T. S.; Müller, E. A. Simulation of asphaltene aggregation through molecular dynamics: insights and limitations. Energy Fuels 2017, 31, 1108−1125. (67) Lv, G.; Gao, F.; Liu, G.; Yuan, S. The properties of asphaltene at the oil-water interface: a molecular dynamics simulation. Colloids Surf., A 2017, 515, 34−40. (68) Wang, S.; Xu, J.; Wen, H. The aggregation and diffusion of asphaltenes studied by gpu-accelerated dissipative particle dynamics. Comput. Phys. Commun. 2014, 185, 3069−3078. (69) Voth, G. A. Coarse-Graining of Condensed Phase and Biomolecular Systems; CRC press, 2008. (70) Peter, C.; Kremer, K. Multiscale simulation of soft matter systems. Faraday Discuss. 2010, 144, 9−24. (71) Aguilera-Mercado, B.; Herdes, C.; Murgich, J.; Müller, E. A. Mesoscopic simulation of aggregation of asphaltene and resin molecules in crude oils. Energy Fuels 2006, 20, 327−338. (72) Boek, E. S.; Headen, T. F.; Padding, J. T. Multi-scale simulation of asphaltene aggregation and deposition in capillary flow. Faraday Discuss. 2010, 144, 271−284. (73) Zhang, S.-F.; Sun, L. L.; Xu, J.-B.; Wu, H.; Wen, H. Aggregate structure in heavy crude oil: using a dissipative particle dynamics based mesoscale platform. Energy Fuels 2010, 24, 4312−4326. (74) Ruiz-Morales, Y.; Mullins, O. C. Coarse-grained molecular simulations to investigate asphaltenes at the oil-water interface. Energy Fuels 2015, 29, 1597−1609. (75) Hernández-Rojas, J.; Calvo, F.; Wales, D. J. Coarse-graining the structure of polycyclic aromatic hydrocarbons clusters. Phys. Chem. Chem. Phys. 2016, 18, 13736−13740. (76) Wang, J.; Ferguson, A. L. Mesoscale simulation of asphaltene aggregation. J. Phys. Chem. B 2016, 120, 8016−8035. (77) Wang, J.; Gayatri, M. A.; Ferguson, A. L. Mesoscale simulation and machine learning of asphaltene aggregation phase behavior and molecular assembly landscapes. J. Phys. Chem. B 2017, 121, 4923− 4944. (78) Desgranges, C.; Delhommelle, J. Prediction of the phase equilibria for island-type asphaltenes via HMC-WL simulations. J. Chem. Phys. 2018, 149, 072307. (79) Jiménez-Serratos, G.; Totton, T. S.; Jackson, G.; Müller, E. A. Aggregation behavior of model asphaltenes revealed from large-scale coarse-grained molecular simulations. J. Phys. Chem. B 2019, 123, 2380−2396. (80) Müller, M.; Katsov, K.; Shick, M. Biological and synthetic membranes: what can be learned from a coarse-grained description? Phys. Rep. 2006, 434, 113−176. (81) Deserno, M. Mesoscopic membrane physics: concepts, simulations, and selected applications. Macromol. Rapid Commun. 2009, 30, 752−771. (82) Schmid, F. Toy amphiphiles on the computer: what can we learn from generic models? Macromol. Rapid Commun. 2009, 30, 741−751.

(83) Subramanian, S.; Simon, S.; Sjö b lom, J. Asphaltene precipitation models: a review. J. Dispersion Sci. Technol. 2016, 37, 1027−1049. (84) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and testing of the opls all-atom force field on conformational energetics and properties of organic liquids. J. Am. Chem. Soc. 1996, 118, 11225−11236. (85) Weeks, J. D.; Chandler, D.; Andersen, H. C. Role of repulsive forces in determining the equilibrium structure of simple liquids. J. Chem. Phys. 1971, 54, 5237−5247. (86) Likos, C. N. Effective interactions in soft condensed matter physics. Phys. Rep. 2001, 348, 267−439. (87) Noid, W. G. Perspective: coarse-grained models for biomolecular systems. J. Chem. Phys. 2013, 139, 090901. (88) Pronk, S.; Páll, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D.; et al. GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit. Bioinformatics 2013, 29, 845−854. (89) Van Gunsteren, W. F.; Berendsen, H. J. C. A leap-frog algorithm for stochastic dynamics. Mol. Simul. 1988, 1, 173−185. (90) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: a linear constraint solver for molecular simulations. J. Comput. Chem. 1997, 18, 1463−1472. (91) Sedghi, M.; Goual, L.; Welch, W.; Kubelka, J. Effect of asphaltene structure on association and aggregation using molecular dynamics. J. Phys. Chem. B 2013, 117, 5765−5776. (92) Murgich, J.; Rodríguez, J.; Aray, Y. Molecular recognition and molecular mechanics of micelles of some model asphaltenes and resins. Energy Fuels 1996, 10, 68−76. (93) Wang, W.; Taylor, C.; Hu, H.; Humphries, K. L.; Jaini, A.; Kitimet, M.; Scott, T.; Stewart, Z.; Ulep, K. J.; Houck, S.; et al. Nanoaggregates of diverse asphaltenes by mass spectrometry and molecular dynamics. Energy Fuels 2017, 31, 9140−9151. (94) Aray, Y.; Hernández-Bravo, R.; Parra, J. G.; Rodríguez, J.; Coll, D. S. Exploring the structure-solubility relationship of asphaltene models in toluene, heptane, and amphiphiles using a molecular dynamic atomistic methodology. J. Phys. Chem. A 2011, 115, 11495− 11507. (95) Wang, J.; Gayatri, M.; Ferguson, A. L. Coarse-grained molecular simulation and nonlinear manifold learning of archipelago asphaltene aggregation and folding. J. Phys. Chem. B 2018, 122, 6627−6647. (96) Jian, C.; Tang, T. One-dimensional self-assembly of polyaromatic compounds revealed by molecular dynamics simulations. J. Phys. Chem. B 2014, 118, 12772−12780. (97) Pacheco-Sánchez, J. H.; Á lvarez-Ramírez, F.; MartínezMagadán, J. M. Morphology of aggregated asphaltene structural models. Energy Fuels 2004, 18, 1676−1686. (98) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual molecular dynamics. J. Mol. Graphics 1996, 14, 33−38. (99) Towns, J.; Cockerill, T.; Dahan, M.; Foster, I.; Gaither, K.; Grimshaw, A.; Hazlewood, V.; Lathrop, S.; Lifka, D.; Peterson, G. D.; et al. XSEDE: accelerating scientific discovery. Comput. Sci. Eng. 2014, 16, 62−74.

6122


Simple Simulation Model for Exploring the Effects of Solvent and

Recommend Documents