Stability and Aggregation of Asphaltenes in AsphalteneâResin

674

Energy & Fuels 2004, 18, 674-681

Stability and Aggregation of Asphaltenes in Asphaltene-Resin-Solvent Mixtures A. Ortega-Rodriguez,⊥ Y. Duda,†,‡ F. Guevara-Rodriguez,† and C. Lira-Galeana*,⊥ Branch of R&D on Deep Water E&P and Branch of Molecular Engineering Research, Mexican Institute of Petroleum, Apartado Postal 14-805, 07730 Mexico City, Mexico, and Institute of Condensed Matter Physics, National Academy of Sciences of Ukraine, Lviv 79011, Ukraine Received October 6, 2003

The Integral Equation Theory (IET) of fluids and molecular dynamics (MD) simulations have been used to study the stability and aggregation of asphaltenes in asphaltene-resin mixtures with different host solvents (n-heptane, toluene, and pyridine). The theory is based on an approximate interaction-potential model that represents the asphaltene-asphaltene (A-A), asphaltene-resin (A-R), and resin-resin (R-R) interactions in a solvent of dielectric constant . This interaction potential is used within the Ornstein-Zernike Hypernetted Chain (OZ-HNC) approximation for calculating radial distribution functions, structure factors, the limits of material stability, and the phase diagrams of asphaltene-containing systems, where the peptizing behavior of resins is studied as a function of the ratio of resin-to-asphaltene molecules in each host solvent. Comparisons of the theoretical OZ-HNC and MD simulations confirm the good accuracy of the proposed approach, which suggests its suitability as an alternative calculation method for predicting the phase behavior of asphaltene-containing oil fluids.

Introduction Asphaltenes are the heaviest fraction of petroleum whose detailed chemical constitution and thermodynamic behavior have been the subject of different investigations over the past 60 years.1-3 These compounds cause serious operational problems in oil recovery and transport, by plugging the reservoir rock or production pipeline, as a consequence of their molecular aggregation and later precipitation/deposition from the oils.4,5 As the oil industry moves toward deeper reservoirs, the probability of encountering asphaltene problems and the costs associated with their remediation will likely increase soon. Asphaltenes are defined as the heaviest components of carbonaceous sources that are insoluble in normal paraffins such as n-pentane or n-heptane but are soluble in an excess amount of aromatics such as benzene or toluene.6,7 They are a collection of polydisperse molecules consisting mostly of polynuclear aromatic units, with varying proportions of aliphatics and alicyclic moieties and small amounts of heteroatoms. Petroleum * Author to whom all correspondence should be addressed. Phone: +52-55-9175-6507. Fax: +52-55-9175-7225. E-mail: [email protected]. ⊥ Branch of R&D on Deep Water E&P, Mexican Institute of Petroleum. † Branch of Molecular Engineering Research, Mexican Institute of Petroleum. ‡ National Academy of Sciences of Ukraine. (1) Nellensteyn, F. J. The Science of Petroleum; Oxford University Press: Oxford, U.K., 1938; Vol. 4, p 2760. (2) Yen, T. F. Prepr.sAm. Chem. Soc., Div. Pet. Chem. 1972, 17 (4), 102. (3) Ortega-Rodriguez, A.; Cruz, S. A.; Gil-Villegas, A.; GuevaraRodriguez, F.; Lira-Galeana, C. Energy Fuels 2003, 17, 1100. (4) Leontaritis, K. J. Presented at the Production Operations Symposium of the Society of Petroleum Engineers, Oklahoma City, OK, March 13-14, 1989; Paper No. SPE 18892. (5) Shields, D. Offshore 2000, (Sept.), 84.

resins, on the other hand, are a fraction of the deasphalted component of crude oils that adsorbs in silica gel and is extracted with the action of polar solvents.8 Two main interpretations on the physical state of asphaltenes in crude oil have been proposed in the literature. One considers asphaltenes as a class of compounds that are dissolved in the surrounding medium (i.e., the rest of the oil), and one must consider that they may precipitate after the oil solubility falls below a certain threshold.9,10 The second hypothesis on the state of asphaltenes in oil fluids considers a specific stabilizing effect of resin molecules: the asphaltenes are considered to be insoluble colloidal solids that are peptized by adsorbed resin molecules in their surface.11-13 According to this approach, the partitioning of the resins between the surfaces of the asphaltene colloids and the supernatant solution is likely to govern the phenomenon of asphaltene separation. Flory-Huggins-type of solution theories and different equations of state are some thermodynamic models that have been proposed to represent the aforementioned two conceptions of asphaltene phase behavior.14-16 Validation of these models is often performed by tuning some (6) Speight, J. G. The Chemistry and Technology of Petroleum; Marcel Dekker: New York, 1991. (7) Standard Test Method for Separation of Asphalt into Four Fractions. ASTM Standard D 4124, American Society for Testing and Materials, Philadelphia, PA, 1988. (8) Koots, J. A.; Speight, J. G. Fuel 1975, 54, 179. (9) Hirschberg, A.; de Long, L. N.; Schipper, B. A.; Meyers, J. G. Soc. Pet. Eng. J. 1984, 24, 283. (10) Andersen, S. I. Fuel Sci. Technol. 1994, 12, 1551. (11) Christy, A. A.; Dahel, B.; Kvalheim, O. M. Fuel 1989, 68, 40. (12) Pfeiffer, J. P.; Saal, R. N. J. Presented at The Sixteenth Colloid Symposium at Stanford University, Stanford, CA, July 6-8, 1939; 139. (13) Murgich, J. Pet. Sci. Technol. 2002, 20, 1029.

10.1021/ef0301699 CCC: $27.50 © 2004 American Chemical Society Published on Web 03/25/2004

Asphaltenes in Asphaltene-Resin-Solvent Mixtures

of the model parameters to (bulk) precipitation data. This parameter estimation method may be limited to the pressure-temperature-composition (P-T-x) conditions considered in the data reduction step; thus, it may be limited for application to a broader set of conditions. As a useful alternative to the aforementioned procedures, one can resort to recent molecular simulation techniques and statistical thermodynamic theories, which provide a way to obtain “experimental” data on the thermodynamic behavior of complex fluid systems at the colloidal length scale.17-19 In contrast to the macroscopic models, molecular theories provide a description to the thermodynamic behavior of a system from microscopic considerations, from which the global phase behavior of a fluid can be “tracked from the roots,” and predictions to a wide variety of conditions are theoretically justified. Two useful tools of particular application within this category are the molecular dynamics (MD) simulations,31,36 and the Ornstein-Zernike (OZ) integral equations with different closures.20 In this work, we present a method to describe the aggregation behavior and phase stability of colloidal asphaltene-resin-solvent systems using a molecularbased approach. Our description begins from the A-A, A-R, and R-R interaction potentials. The parameters of such potentials are also adjustable; however, those adjustments are made only at the initial (molecular) stage of our study. After the basic (microscopical) condition of the system is defined, the macroscopic behavior of the system (i.e., the asphaltene aggregation states, the dispersing power of solvents, and the phase diagrams of asphaltene precipitation) is readily obtained by solving a set of integral equations from statistical thermodynamics, whose foundation is particularly useful in describing the properties of colloidal systems.20,21 The proposed method is shown to be able to reproduce all observed trends of the experimental (macroscopic) behavior of asphaltene precipitation in real fluids. The results of our method are compared to MD simulations. Good accuracy of the Integral Equation Theory (IET) method is observed, in comparison to the MD simulations. We first describe the method used to obtain the interaction potentials for the A-A, A-R, and R-R pairs that are embedded in a host medium (solvent) using molecular mechanics (MM) calculations. The calculation of the pair distribution functions and structure factors by solving the multicomponent Ornstein-Zernike Hypernetted Chain Approximation (OZ-HNC) integral equations is later described. A comparison of these results with those predicted from MD simulations is presented and discussed. The equations for calculating the locus of points of the limit of intrinsic stability from the calculated structure factors is outlined. Sample (14) Burke, N. E.; Hobbs, R. E.; Kashou, S. F. J. Pet. Technol. 1990, 42, 1440. (15) Nghiem, L. X.; Coombe, D. A. SPE J. 1997, 2, 170. (16) Buenrostro-Gonzalez, E.; Gil-Villegas, A.; Wu, J.; Bazua, E.; Lira-Galeana, C. AIChE J., in press, 2004. (17) Murgich, J.; Rodriguez, J.; Aray, Y. Energy Fuels 1996, 10, 68. (18) Rogel, E. Langmuir 2002, 18, 1928. (19) Rogel, E. Energy Fuels 1997, 11, 920. (20) Lee, L. L. Molecular Thermodynamics of Nonideal Fluids; Butherworth: Boston, 1988. (21) McQuarrie, D. A. Statistical Mechanics; Harper and Row: New York, 1976.

Energy & Fuels, Vol. 18, No. 3, 2004 675

calculations in a model oil suspension in polar and nonpolar solvents are finally presented. Theory Intermolecular Potentials. We have investigated the stability of a binary asphaltene-resin mixture interacting through the effective spherical potential Uij(r), given by the expression derived by Ortega-Rodriguez et al.:22

Z1Z2 -Rr C2e-γ/r -βr [e Uij(r) ) + (C1re )] r 2r6

4

(1)

where r is the distance between the centers of mass of each molecule, and Z1 and Z2 are the number of atoms in each asphaltene or resin molecule. The parameters in eq 1 (e.g., C1,C2, R, β, and γ) were obtained by OrtegaRodriguez et al. by fitting eq 1 to the corresponding MM numerical results as predicted from the INSIGHT II and DISCOVER 3.0 suite of programs23 for the model asphaltene and resin molecules shown in Figure 1. Details are given elsewhere.22 In the MM calculation, the multidimensional Newton equations of motion are solved for a set of different relative orientations and distances between each pair of molecules. The expression in eq 1, which distinguishes among Coulomb, dispersion, and London contributions, was obtained by fitting its predictions to the corresponding simulation results. The assumptions in eq 1 and its satisfactory application to the study of the dissolution ability of different solvents in model dispersions and to predict the precipitation phase diagram of a problematic asphaltenic live-oil fluid have been presented elsewhere.3,16 Figure 2 shows the functional representation of eq 1 when toluene is the host solvent. Solution to the OZ Integral Equation. Within the framework of the IET of fluids, two key quantities of a binary asphaltene-resin system are the corresponding total and direct correlation functions (hij(r) and cij(r), respectively,20 where the subscripts i and j indicate the component type; “A” will refer to asphaltenes and “R” will refer to resins). For these functions, one can write the following OZ integral equation:20 A,R

hij(r) ) cij(r) +

∑k Fk∫ drb′hik(|rb - br ′|)ckj(rb′)

(2)

where Fk is the number density of the k species (asphaltene or resin). The concentration of asphaltene is defined as

χ)

NA NA + NR

(3)

where Ni is the number of particles of the type i, and N ) NA + NR. Then, FA ) χF and FR ) (1 - χ)F, where F is the total number density of particles. To solve eq 2, and to calculate the pair distribution functions gij(r) and the partial structure factors Sij(k), we use the Hypernetted Chain (HNC) closure relation:20 (22) Ortega-Rodriguez, A.; Cruz, S. A.; Ruiz-Morales, Y.; LiraGaleana, C. Pet. Sci. Technol. 2001, 19, 245. (23) Molecular Simulations, Inc., San Diego, CA, 2000.

676

Energy & Fuels, Vol. 18, No. 3, 2004

Ortega-Rodriguez et al.

Figure 2. Intermolecular potentials in toluene ( ) 2.38) calculated from eq 1.

sion of ∼10-7. Most of the calculations were performed on a grid with a real-space step of ∆r ) 0.01σA and a total number of points of N0 ) 4096. We have performed some checks of the numerical stability of the results for those cases where convergence was more difficult to obtain, using N0 values of 2048 and 8192 and various ∆r. No significant variation was observed with these different grids. Calculation of Pair-Correlation Functions and Structure Factors. From the knowledge of hij(r) (and hˆ ij(k)), the radial distribution functions and the partial structure factors are calculated. The partial structure factors are given by the standard relation20

Sij(k) ) δij + xFiFjhˆ ij(k)

Figure 1. Molecular models for asphaltenes and resins used in this work: (a) three-dimensional view of an asphaltene molecule reported in ref 37, (b) transverse view of an asphaltene molecule reported in ref 37, and (c) a resin molecule proposed in ref 38.

gij(r) ) hij(r) + 1 ) exp(-βTUij(r) + θij(r))

(4)

where we introduce the well-known auxiliary function θ(r) ) h(r) - c(r); βT ) (kBT)-1. Equations 2 and 4 are solved together for a given density F, concentration χ, and organic solvent (dielectric constant in eq 1) through an iterative procedure that has been described in detail.20 Here, we only briefly discuss the solution algorithm: all r functions are discretized in N0 regularly spaced points ri ) i∆r (∆r is the space step). The functions are assumed to be zero above the cutoff distance N0∆r. Similarly, the k functions are known at the points kj ) j∆k, where ∆k ) 2π/(N0∆r). A cycle is started with an auxiliary function, θ0ij(r). This function has a faster converging rate, in particular, for extremely difficult cases that correspond to the region near the spinodal line. The integral equation gives cij(r) and cˆ ij(k) (through a Fourier transform). The OZ equation in k-space yields θij(k). Finally, an inverse Fourier transform gives the output function θ1ij(r). The process is iterated until the input “0” and output “1” functions converge with a prescribed preci-

(5)

where δij is the Kronecker delta. From the partial structure factors, the concentrationconcentration structure factor Sχχ(k) is calculated through the following relation:

Sχχ(k) ) χ(1 - χ)[(1 - χ)SAA(k) + χSRR(k) 2xχ(1 - χ)SAR(k)] (6) Note that the structure factors predicted by eqs 5 and 6 can, in principle, be measured by small-angle X-ray scattering (SAXS) or small-angle neutron scattering (SANS) techniques,24 from which further comparisons between the predicted and experimental results of Sχχ(k) at the colloidal length scale can be optionally performed. Calculation of the Stability Limit. As it has been shown,20 the long-wavelength limit of concentrationconcentration structure factor, Sχχ(k f 0), is proportional to the osmotic compressibility and is related to the Gibbs free energy of mixing through the following relation:

( ) ∂2GM ∂χ2

) T,P

NβT-1 Sχχ(0)

(7)

In the presence of a phase separation, large local fluctuations occur in the mixture, so that the osmotic compressibility and, consequently, Sχχ(0) will diverge, which implies that (24) Lindner, P.; Zemb, Th. Neutrons, X-Rays and Light; Elsevier Science Publishers: Amsterdam, 2002.


( ) ∂2 G M ∂χ2

)0


(8)

T,P

The locus of points on the F-χ plane in which eq 8 is satisfied is known as the spinodal curve of the asphaltene-resin-solvent mixture; i.e., whenever two minima are observed in the energy of mixing, the system is unstable.25 Results and Discussion To locate the concentration for demixing of an asphaltene-resin dispersion, we search for singularities of Sχχ(0) by solving the integral equation at a fixed total fluid density F for a sequence of physical states with increasing asphaltene concentration, beginning with χ ) 0. Unlike other works,26,27 where the analytical solution of the Percus-Yevick approximation has been used, in the case of the HNC approximation, one cannot approach the state where Sχχ(0) f ∞. We first considered an asphaltene-resin mixture that had been dissolved in toluene ( ) 2.38). By reporting the values of Sχχ(0) as a function of the asphaltene concentration at sufficiently low fluid densities, we obtained a bell-shaped curve (see Figure 3), which is obviously zero at the pure-resin component limit and, in the present case of toluene as a solvent, has a maximum at χ ≈ 0.04. As the density F increases, and the spinodal state is approached, Sχχ(0) increases for every value of χ, at a rate that is higher, with respect to the maximum of the curve. The rate of numerical convergence becomes slower and slower, until no solution of eq 2 can be found. For the fluid density F ) 0.0124, it is impossible to obtain the HNC solution of eq 2 in a range of concentrations located approximately at χ ≈ 0.03. This region becomes wider and wider as the density increases. However, we did not attempt to make any asymptotic analysis of these behaviors, which could had been attempted by assuming that Sχχ(0) diverges according to a power law and trying to determine the critical exponent in the HNC case. A careful study of these details goes beyond the scope of this paper, but perhaps it would require a large number of integral-equation results very near the critical points, i.e., those points that are stable on the limit of stability. The spinodal curves obtained as a limit of HNC solution are presented in Figure 4. As shown in this figure, asphaltenes are more soluble in toluene than in heptane, as observed in the experiments.28 In the case of pyridine, which is known to be a strong solvent ( ) 12.3), we were unable to find a region of asphaltene precipitation; e.g., total dissolution, according to the experimental evidence, is predicted. The inset figure in Figure 4 presents the spinodal curves in the P-χ plane. The asphaltene-resin osmotic pressure is calculated using the virial equation of state:20 (25) Sheu, E. Y. Energy Fuels 2002, 16, 74. (26) Duda, Y.; Lira-Galeana, C. J. Colloid Interface Sci., 2003, in review. (27) Duda, Y. J. Colloid Interface Sci. 1999, 213, 498. (28) Buenrostro-Gonzalez, E.; Andersen, S. I.; Garcia-Martinez, J. A.; Lira-Galeana, C. Energy Fuels 2002, 16, 732.

Figure 3. Infinite wavelength limit of the concentrationconcentration structure factor Sχχ(0), as a function of the asphaltene concentration, χ, in toluene at different mixture densities (from F ) 0.0126 to F ) 0.012).

Figure 4. Spinodal curves obtained as a limit of the OZ-HNC solution. Inset presents the values of osmotic pressure on the spinodal curves calculated from the virial eq 9. Dashed and solid curves correspond to heptane and toluene, respectively.

βT P F

)1-

βT

2

2

∑∑

6 i)1j)1

FiFj

∫

dUij(r) 4πr3 gij(r) dr (9) 0 dr r

where Uij(r) and gij(r) are given by eqs 1 and 4, respectively. Besides examining the behavior of Sχχ(0), which is a “global” measure of the deviation from the mean concentration value, it may be more interesting to gain insight into the “local” structure of a mixture near phase separation or near the boundary of material stability, by examining the shape of the radial distribution functions, gij(r). Figure 5 shows the radial distribution function of the asphaltene-asphaltene pair, gAA(r), in toluene. At a fixed mixture density of F ) 0.12, four different asphaltene concentrations have been considered. As one can see, the approaching of the asphaltene concentration to the critical value leads to an increase in oscillations, i.e., at χ ) 0.01588 and 0.05405, gAA(r) has long-range correlations above unity, which is a clear indication of asphaltene aggregation.27 A suitable method to investigate the local composition of a mixture is based on the analysis of the local mole

678



Figure 6. OZ-HNC results for the sum of the local mole fractions xs(r) as a function of asphaltene concentration χ at a mixture density of F ) 0.0122 in toluene. Inset shows the A-A coordination number (nAA(r)) under the same conditions. Dashed, dash-dotted, and solid lines correspond to concentrations of 0.1, 0.06, and 0.04625, respectively. Figure 5. Radial distribution functions of the A-A pair in toluene obtained from HNC solution at a mixture density of F ) 0.012 and different asphaltene concentrations χ. Panels a and b show the approach to the two-phase region from the right-hand and left-hand sides of the bell curve of Figure 4, respectively.

fractions. The local mole fraction xii(r0) of a binary system is defined by30

xij(r0) )

nij(r0) njj(r0) + nij(r0)

(for i ) A,R; j * i) (10a)

xjj(r0) ) 1 - xij(r0)

(10b)

where

nij(r0) ) 4πFj

∫0r

0

dr r2gij(r)

(11)

is a running coordination number, which gives the average number of type-i particles surrounding a type-j particle, out to a distance r0. Therefore, xij(r0) represents the fraction of atoms of species i, which are in a sphere of radius r0 around a type-j atom that is taken as the origin. The sum of the local mole fractions taken at an appropriate distance r0,

xs(r) ) xAA(r0) + xRR(r0)

(12)

is used to characterize the deviations of the mean number of particles of the two types corresponding to specific cross interactions between unlike molecules. For an ideal mixture, xs(r0) ) 1 for every value of r0, because the two types of particles are equally distributed in proportion to the composition of the mixture. For strongly nonideal systems, we expect a microscopic tendency toward mixing or demixing. In the case of (29) Duda, Y.; Vakarin, E.; Alejandre, J. J. Colloid Interface Sci. 2003, 258, 10. (30) Nakanishi, K.; Okazaki, S.; Ikari, K.; Higuchi, T.; Tanaka, H. J. Chem. Phys. 1982, 76, 629.

demixing phase separation, the local average number of pairs of particles of the same type predominates over the number of pairs of unlike particles; therefore, the local mole fractions have larger values than those expected for an ideal system and the resulting value for the sum xs(r0) is larger than unity. In this work, we present the local mole fraction as a function of distance. The result shown in Figure 6 shows a tendency toward demixing, with a separation of the A-A and R-R pairs when r g 10 Å. The drastic increase in xs(r) as the right stability boundary is approached clearly indicates that the mixture tends toward a microscopically inhomogeneous state. Interestingly, the asphaltene running coordination number behavior (see the inset of Figure 6) confirms the drastic agglomeration of asphaltene when r g10 Å. To determine the boundary curves of asphalteneresin mixtures directly, one must perform a Grand Canonical or Gibbs ensemble Monte Carlo simulation,31 which is well-known to be a very time-consuming task, taking into account the finite-size-effect problems. In this work, we have performed MD simulations, which, in an indirect way, confirms the results obtained from the OZ-HNC solution. Most features of the canonical MD simulation used here have previously been detailed elsewhere.3 Accordingly, we confine the description of the simulation to its barest essentials, except where necessary to detail a new aspect. All simulations were performed in a cubic box with periodic boundary conditions, at a temperature of T ) 273 K, with a time step of 1.0 fs. In Figure 7, a comparison of structure predicted by the OZ-HNC approach and that obtained from a MD simulation is presented for two different sets of mixture density and composition. As can be seen, in both cases considered, the integral equation approach reproduces the behavior of gij(r) well. Thus, one may expect that the boundary curves obtained within the OZ-HNC (31) Frenkel, D.; Smith, B. Understanding Molecular Simulation; Academic Press: London, 1996.



Figure 7. (a) Radial distribution functions obtained from the OZ-HNC results (lines) and MD simulations (symbols) for (a) a mixture density of F ) 0.009 and an asphaltene concentration of χ ) 0.19825 and (b) a mixture density of F ) 0.0105 and an asphaltene concentration of χ ) 0.12. For both conditions, ) 1.92.

approximation (Figure 4) are of the same level of accuracy. To verify this hypothesis, we have performed the MD simulation runs for two fixed asphaltene concentrations of χ ) 0.19825 and 0.12 in toluene (68 and 58 asphaltene particles in a cubic box, respectively), changing the mixture density. During the simulation, the radial distribution functions, the mean energy (βTE), the concentration of asphaltene monomer ) Nmon (Xmon A A /NA), the normalized frequency distribution of asphaltene cluster size (f(S)), and asphaltene mean cluster sizes (Z) have been monitored. The algorithm of Sevick et al.32 has been used to calculate f(S), and Z, as discussed elsewhere,33 can be represented by the relation NA

Z)

∑1 S2ns NA

∑1 SnS

(13)

where S denotes the cluster size and ns is the mean number of clusters of size S (S > 1). The term f(S) is defined as f(S) ) ΣW j nS(j)/W, where W is the total number of MD configurations. By definition, each asphaltene particle within a cluster is connected (either directly or indirectly) to every other particle in the cluster. The determination of whether particles are directly connected follows from the definition of bound pairs and is easily determined by sampling all distinct pairs of asphaltenes in the system. In the present work, we have considered two asphaltene particles to be directly connected if the distance between them is smaller than rbAA ) 6.5 Å. We have verified (by considering a larger “bound” distance, rbAA ) 7 Å) that our conclusions are not dependent on the choice of rbAA. Note that the range of density values (F ) 0.006-0.012) has been chosen in such a way to check the one-phase and two-phase regions, in accordance with the OZ-HNC predictions in Figure 4. The radial distribution functions for the pairs A-A, A-R, and R-R obtained from the simulation in the case of toluene as the host medium for χ ) 0.19825 and different fluid densities are shown in the left-hand part of Figure 8. As expected, the first peak in gAR(r) and gRR(r) grows as the density increases. However, the value of the first peak of gAA(r) has a minimum at F ) 0.009. This interesting feature will be related below to other measurements that indirectly indicate a transformation of the system from a homogeneous state to an inhomogeneous state at F ≈ 0.009 (χ ) 0.19825). The right-hand portion of Figure 8 presents the frequency of the asphaltene cluster size, f(S). One can observe that the frequency of small asphaltene clusters increases as the density decreases and has a maximum value at F ) 0.0105 (almost the same as that at F ) 0.009) and then decreases. Meanwhile, both the cluster size and frequency of big clusters increase. For the smallest density considered (F ) 0.006), the highest frequency of the biggest cluster of 68 particles is observed. Recently, qualitatively similar behavior of the cluster-size frequency of a simple model mixture34 and of real asphaltene dispersion35 has been observed. As an illustration of the extent of predicted asphaltene aggregation and asphaltene-resin separation at this density, Figure 9 shows a snapshot of one of the MD configurations. Note that the same “visual” analysis of MD configurations, at densities of F > 0.008 (for χ ) 0.19825) and F > 0.009 (for χ ) 0.12), did not permit us to find the asphaltene-resin separation, because the stability of the liquid-liquid interface close to the spinodal curve is affected by high concentration fluctuations, especially in the cubic simulation box.29 (32) Sevick, E. M.; Monson, P. A.; Ottino, J. M. J. Chem. Phys. 1988, 88, 1198. (33) Stauffer, D.; Aharony, A. M. Introduction to Percolation Theory; Taylor and Francis: London, 1994. (34) Zamudio, A.; Duda, Y.; Medina-Noyola, M. Phys. Lett. A 2002, 305, 258. (35) Carnahan, N. F.; Quintero, L.; Pfund, D. M.; Fulton, J. L.; Smith, R. D.; Capel, M.; Leontaritis, K. Langmuir 1993, 9, 2035. (36) Tuckerman, M.; Berne, B. J.; Martyna, G. J. J. Chem. Phys. 1990, 97, 3. (37) Zajac, G. W.; Sethi, N. K.; Joseph, J. T. Scanning Microsc. 1994, 8 (3), 463. (38) Murgich, J.; Abanero, J. A.; Strausz, O. P. Energy Fuels 1999, 13, 278.

680



Figure 8. Radial distribution functions (left-hand side) and frequency of asphaltene cluster sizes (right-hand side) obtained from MD simulation results at different mixture densities in toluene. Asphaltene concentration is χ ) 0.19825.

In Figure 10, we present the results obtained from MD simulations for mean energy, concentration of asphaltene monomer, and asphaltene mean cluster sizes. All the curves have peculiarities near a mixture density of F ) 0.009. From Figure 10, one can conclude that the formation of big clusters of asphaltenes is

inhibited by concentration effects, and that the mostremarkable changes in the cluster size correspond to concentrations of 0.009, the cluster size still decreases, but its size changes very little. Moreover, the concentration value of 0.009 corresponds to a minimum in the mean energy; therefore,


Figure 9. Snapshot of a MD-derived system configuration at F ) 0.006 and χ ) 0.19825.

Figure 10. Results of MD simulation for the mean energy (βTE), concentration of asphaltene monomer (Xmon A ), and asphaltene mean cluster sizes (Z), as a function of the mixture density, at an asphaltene concentration of χ ) 0.19825.

this value optimizes the formation of clusters. On the other hand, in Figure 11 one can see that the asphaltene diffusion decreases as concentration increases, mainly by the collective effect of the cluster to which it belongs. Interestingly, the behavior of the asphaltene mean cluster size resembles the behavior of the self-diffusion coefficient in Figure 11, which indicates the relation between asphaltene aggregation and its diffusion, i.e., an increase of Z with mixture density leads to an increase in the self-diffusion coefficient. This could mean that, after separation of the resin and asphaltene molecules, the unpeptized asphaltene particle feels more free to move in the pure asphaltene neighborhood. When toluene acts as the solvent, the densities of the asphaltene-resin dispersion instability are F ) 0.008 ( 0.001 (for χ ) 0.19825) and F ) 0.009 ( 0.001 (for χ ) 0.12). As one can see, these data are in very good agreement with the results of the OZ-HNC solution shown in Figure 4.


Figure 11. Self-diffusion coefficient of asphaltenes dissolved in toluene at χ ) 0.198251 (solid lines) and χ ) 0.12 (dashed lines), calculated by MD simulations, as a function of the mixture density F. Triangles and circles correspond to simulation times of 10 and 20 ps, respectively.

in asphaltene-resin-solvent mixtures using an Integral Equation Theory (IET) approach (the Ornstein-Zernike Hypernetted Chain (OZ-HNC) approximation) and molecular dynamics (MD) simulations. This method is particularly suited to model asphaltene aggregation phenomena in oil fluids or model suspensions at the colloidal-length scale, where structural data from smallangle X-ray scattering (SAXS) or small-angle neutron scattering (SANS) experiments could be available. The present method requires approximate characterization data of the asphaltene and resin molecular units in oil fluids as input to the model and provides a detailed and fully predictive monitoring of the aggregation behavior of asphaltene-containing fluids as a function of solvent properties, density, and composition. The reliability of our IET (OZ-HNC) predictive method for the onset of phase stability and aggregation in asphaltene-resinsolvent model mixtures is supported by a reasonably good agreement with both experimental evidence and MD simulation. Further results are available upon request.39 The method presented in this paper is simple, and it can be used as a useful alternative for estimating the precipitation/aggregation properties of asphaltenecontaining petroleum fluids in reservoir, production, and transport engineering. Acknowledgment. This work was supported by the Molecular Engineering Research Branch of the Mexican Institute of Petroleum (IMP), under Contract Nos. D.00337, R.00077, and D.00324. We thank Dr. Eric Sheu of Vanton Research Laboratory (Lafayette, CA) for stimulating discussions. A. O.-R. is indebted to the National Council of Science and Technology of Mexico (CONACyT) and to IMP for a fellowship. EF0301699

Concluding Remarks We have presented, for the first time, a thorough theoretical analysis of phase stability and aggregation

(39) Ortega-Rodriguez, A. Thermodynamics of Asphaltene Precipitation in Petroleum Mixtures. Ph.D. Thesis, National University of Mexico, Mexico City, Mexico, 2004.

Stability and Aggregation of Asphaltenes in AsphalteneâResin

Recommend Documents

Stability and Aggregation of Asphaltenes in AsphalteneâResin