Energy Fuels 2010, 24, 4952–4960 Published on Web 08/30/2010
: DOI:10.1021/ef100123x
A Study of the Dynamic Evolution of Asphaltene Aggregate Size Distribution Using Monte Carlo Simulation Marjan Faraji and Ali R. Solaimany Nazar* Chemical Engineering Department, University of Isfahan, Isfahan, IRAN Received February 1, 2010. Revised Manuscript Received August 13, 2010
Monte Carlo (MC) simulation is used to determine the dynamic evolution of the asphaltene aggregates size distribution (ASD) in shear-induced organic solvents. This MC method is a time-driven simulation with periodic regulation of the number of particles. Asphaltene molecules are referred to as fractal-like aggregates that undergo simultaneous aggregation and breakage processes in shear-induced organic solvents. The effects of shear rate, initial particle size, and initial asphaltene concentration on the evolution of number average diameter are studied. The simulations are validated through comparison with the experimental data of the asphaltene ASD and numerical results of deterministic solution of population balance. These comparisons revealed that the MC simulation is capable of predicting the dynamic evolution of the ASD with reasonable accuracy.
The development of prediction models for asphaltene aggregation under flow conditions is important to the petroleum industry as it affects a variety of processing issues from oilfield production to the refinery. Evaluation of the inhibitor compounds performance and design of the proper operating conditions for avoiding asphaltene precipitation and deposition severely depend on the asphaltene aggregates size distribution (ASD). Therefore, a valid means of predicting asphaltene ASD is necessary and highly desirable.11 The dynamic evolution of the particle size distribution (PSD) in colloidal systems could be predicted from the solution of the population balance equation (PBE).12,13 PBE is a typical partially integro-differential equation. In the presence of aggregation and breakage, this equation is as follows:14
Introduction A major problem in the petroleum industry is the precipitation of asphaltene compounds. Asphaltenes are found in colloidal suspensions in oil mixtures and are stabilized by resin molecules. Under certain thermodynamic conditions, asphaltenes may be precipitated in the form of large aggregates. These aggregates can result in complex operating problems in production, transportation, and storage processes of crude oil. Asphaltene aggregates are not rigid spheres. As aggregation proceeds, particles grow as porous objects with highly irregular and open structures. They are usually referred to as fractal aggregates, characterized by their mass fractal dimensions.1 The fractal dimension is a measure of how the constituent particles of the aggregate fill the space. It was demonstrated that asphaltene aggregates generated in the flocculation/coagulation process possess fractal features.2 They have fractal dimensions in the range of 1.06-2.5 based on the aggregation mechanism.2-8 Asphaltene aggregate fractal dimensions have been calculated from various smallangle neutron scattering (SANS) studies to vary from 1.7 to 3 depending on the solvent, concentration, crude source, temperature, etc.9,10
Here, v is the particle volume, taken to represent size; n(v) is the size distribution such that n(v) dv is the number concentration of particles in the size range from v to (v þ dv); β(v, u) is the aggregation kernel constant between sizes v and u; S(u) is the breakage kernel of size u; and Γ(v|u) is the breakup distribution function of the particle of volume v resulting from the breakup of a particle of volume u. The first part in the ‘‘aggregation’’ term on the right-hand side of eq 1 represents
*To whom correspondence should be addressed. Phone: 0098-3117934027. Fax: 00983117934031. E-mail:
[email protected]. (1) Somasundaran, P.; Runkana, V. Int. J. Miner. Process. 2003, 72, 33–55. (2) Rahmani, N. H. G.; Dabros, T.; Masliyah, J. H. J. Colloid Interface Sci. 2005, 285, 599–608. (3) Sahimi, M.; Rassamdana, H.; Dabir, B. SPE J. 1997, 2, 157–169. (4) Rassamdana, H.; Sahimi, M AIChE J. 1996, 42, 3318–3332. (5) Sahimi, H.; Rassamdana, H. Physica A 1997, 240, 419–431. (6) Dabir, B.; Nematy, M.; Mehrabi, A. R.; Rassamdana, H.; Sahimi, M. Fuel 1996, 75, 1633–1645. (7) Rastegari, K.; Svrcek, W. Y.; Yarranton, H. W. Ind. Eng. Chem. Res. 2004, 43, 6861–6870. (8) Rahmani, N. H. G.; Dabros, T.; Masliyah, J. H. Energy Fuels 2005, 19, 1099–1108. (9) Tanaka, R.; Hunt, J. E.; Winans, R. E.; Thiyagarajan, P.; Sato, S.; Takanohashi, T. Energy Fuels 2003, 17, 127–134. (10) Gawrys, K. L.; Blankenship, G. A.; Kilpatrick, P. K. Langmuir 2006, 22, 4487–4497. r 2010 American Chemical Society
(11) Carlos, L.; Junior, R.; Ferreira, M. S.; Ramos, A. C. S. J. Pet. Sci. Eng. 2006, 51, 26–36. (12) Bove, S.; Solberg, T.; Hjertager, B. H. Chem. Eng. Sci. 2005, 60, 1449–1464. (13) Meimaroglou, D.; Roussos, A. I.; Kiparissides, C. Chem. Eng. Sci. 2006, 61, 5620–5635. (14) Blatz, P. J.; Tobolsky, A. V. J. Phys. Chem. 1945, 49, 77–79.
4952
pubs.acs.org/EF
Energy Fuels 2010, 24, 4952–4960
: DOI:10.1021/ef100123x
Faraji and Solaimany Nazar
the generation of the particle of volume v by the aggregation of smaller particles. The second part represents consumption of particle of volumex v by aggregation with other particles. The second term describes the gain and loss in number concentration respectively due to breakage. The first part in the ‘‘breakage’’ term describes the breakage of the particle of volume u as being bigger than the particle of volume v to yield particle of volume v, and the second part represents the breakage of particle of volume v with a statistical probability of S(v). In general, the available methods for solving these equations can be divided into deterministic and stochastic methods. The most popular deterministic methods, such as the method of moments,15-26 sectional method,27 discrete method,28 and discrete-sectional method29,30 are used to calculate the integrals of the population balance equation either through an appropriate discretization scheme or by quadrature. The method of moments is based on calculating some selected moments of the distribution instead of the entire distribution. The application of the moment definition to the overall population balance equation to derive results one ordinary differential equation representing each moment (say, 0, 1, 2, 3, etc.). Although the solution of the resultant system of ODEs is simple, the application of the technique requires the number of moments considered to be limited. Sectional method approximates the continuous PSD by a finite number of sections within which the PSD function is assumed to be constant. These types of methods demand the use of some ad hoc rules in order to satisfy the conservation of certain properties of the distribution, such as mass conservation, and at the same time to reduce the number of points required to represent the distribution.31 The most popular techniques for the solution of PBEs in general, and population balance models in particular, are the discretization techniques. In this case, the PBE is sought to be identically satisfied at a certain number of points only (at certain values of v within the domain spanning vmin and vmax). The simplest among the discretization techniques is the finite difference methods. In this technique, the derivatives (usually with respect to the internal coordinate v alone, but sometimes with respect to both t and v), are approximated by finite differences of a certain order (typically 1;Taylor’s expansion, or 2). This approximation results in either a system of
ODEs or algebraic equations. The discrete-sectional method is a combination of the discrete method and the sectional method. The method utilizes discrete method with finite number of bins to describe the lower end of PSD (discrete size region) and sectional method with proper number of sections to account for the other part of PSD (sectional size regime).The dynamic evolution of the PSD in a particulate process can also be obtained via stochastic Monte Carlo (MC) simulations. The stochastic methods do not require discretization of these integrals. This significantly simplifies the programming effort compared to normal numerical methods, even considering multiple mechanisms such as aggregation and breakage. In addition, their application to higher dimensional problems can be easily realized.32-34 Monte Carlo simulation methods have also been used to help elucidate the molecular structure. The disadvantage of the Monte Carlo method is that it is comparatively timeconsuming. The basis of MC simulation is the generation of a random number for an unknown variable in the target system. Dynamic evolution of the particle size distribution was studied by the properties of a finite sample of particle population.35-45 Applications of MC technique in agglomeration,46 restructuring,33,47 multicomponent aerosols,48 chemical reaction,49 and crystallization32 have been presented in recent years. Experimental measurements and theoretical predictions of the aggregation-breakage behavior of asphaltene aggregates in shear induced organic solvents have been the focus of a few recent works.2,50-54 The major motivations of the method development in the previous studies were based on different issues related to the effect of the fractal nature of the asphaltene structure on the asphaltene ASD. Asphaltene aggregates are assumed to be porous particles with constant porosity.50 (32) Gooch, J. R.; Hounslow, M. J. AIChE J. 1996, 42, 1864–1874. (33) Tandon, P.; Rosner, D. E. J. Colloid Interface Sci. 1999, 213, 273–286. (34) Rosner, D. E.; McGraw, R.; Tandon, P. Ind. Eng. Chem. Res. 2003, 42, 2699–2711. (35) Gillespie, D. T. J. Atmos. Sci. 1972, 29, 1496–1510. (36) Shah, B. H.; Ramkrishna, D.; Borwanker, J. D. AIChE J. 1977, 23, 897–904. (37) Goodson, M.; Kraft, M. J. Comput. Phys. 2002, 183, 210–232. (38) Vikhansky, A.; Kraft, M. J. Comput. Phys. 2004, 200, 50–59. (39) Kruis, F. E.; Maisels, A.; Fissan, H. AIChE J. 2000, 46, 1735– 1742. (40) Maisels, A.; Einar Kruis, F.; Fissan, H. Chem. Eng. Sci. 2004, 59, 2231–2239. (41) Zhao, H.; Zheng, C. Atmos. Environ. 2006, 40, 1510–1525. (42) Zhao, H.; Zheng, C.; Xu, M. Powder Technol. 2005, 154, 164– 178. (43) Zhao, H.; Zheng, C.; Xu, M. J. Colloid Interface Sci. 2005, 286, 195–208. (44) Zhao, H.; Zheng, C.; Xu, M. Appl. Math. Mech. 2005, 26, 953– 962. (45) Lin, Y.; Lee, K.; Matsoukas, T. Chem. Eng. Sci. 2002, 57, 2241– 2252. (46) Matsoukas, T.; Friedlander, S. K. J. Colloid Interface Sci. 1991, 146, 495–506. (47) Rosner, D. E.; Yu AIChE J. 2001, 47, 545–561. (48) Efendiev, Y.; Zachariah, M. J. Colloid Interface Sci. 2002, 249, 30–43. (49) Gillespie, D. T. J. Comput. Phys. 1976, 22, 403–434. (50) Rahmani, N. H. G.; Dabros, T.; Masliyah, J. H. Chem. Eng. Sci. 2004, 59, 685–697. (51) Rahmani, N. H. G.; Dabros, T.; Masliyah, J. H. Ind. Eng. Chem. Res. 2005, 44, 75–84. (52) Solaimany Nazar, A. R.; Rahimi, H. Energy Fuels 2008, 22, 3435–3442. (53) Solaimany Nazar, A. R.; Rahimi, H. Energy Fuels 2009, 23, 967– 974. (54) Rahimi, H.; Solaimany Nazar, A. R. Energy Fuels 2010, 24, 1088–1093.
(15) Sandu, A. Atmos. Environ. 2002, 36, 583–589. (16) Fan, R.; Marchisio, D. L.; Fox, R. O. Powder Technol. 2004, 139, 7–20. (17) Wright, D. L.; McGraw, R.; Rosner, D. E. J. Colloid Interface Sci. 2001, 236, 242–251. (18) McGraw, R. Aerosol Sci. Technol. 1997, 27, 255–265. (19) Zucca, A.; Marchisio, D. L.; Barresi, A. A.; Fox, R. O. Chem. Eng. Sci. 2006, 61, 87–95. (20) Marchisio, D. L.; Vigil, D. R.; Fox, R. O. Chem. Eng. Sci. 2003, 58, 3337–3351. (21) Dorao, C. A.; Jakobsen, H. A. J. Comput. Appl. Math. 2006, 196, 619–633. (22) Alopaeus, V.; Laakkonen, M.; Aittamaa, J. Chem. Eng. Sci. 2006, 61, 4919–4929. (23) Marchisio, D. L.; Vigil, R. D.; Fox, R. O. J. Colloid Interface Sci. 2003, 258, 322–334. (24) Marchisio, D. L.; Fox, R. O. J. Aerosol. Sci. 2005, 36, 43–73. (25) Diemer, R. B.; Olson, J. H. Chem. Eng. Sci. 2002, 57, 2211–2228. (26) Diemer, R. B.; Olson, J. H. Chem. Eng. Sci. 2002, 57, 2193–2209. (27) Gelbard, F.; Seinfeld, J. H. J. Colloid Interface Sci. 1980, 78, 485–501. (28) Landgrebe, J.; Pratsinis, S. Ind. Eng. Chem. Res. 1989, 28, 1474– 1481. (29) Wu, J.; Flagan, R. J. Colloid Interface Sci. 1988, 123, 339–352. (30) Landgrebe, J. D.; Pratsinis, S. E. J. Colloid Interface Sci. 1990, 139, 63–86. (31) Ramkrishna, D. Population Balances; Academic Press: San Diego, 2000.
4953
Energy Fuels 2010, 24, 4952–4960
: DOI:10.1021/ef100123x
Faraji and Solaimany Nazar
However, in general, from a fractal point of view, aggregates have various porosities, and the higher aggregate diameter results in higher porosity.52,53 Recently, an asphaltene aggregates fractal restructuring model was proposed to include the restructuring effects of shear induction into the asphaltene aggregates. It was assumed that aggregates fractal dimension changes during shear induction and that its variations could be described by a time-dependent function.54 Other authors made contributions to the simulation of interactions of asphaltene aggregates.55,56 Aguilera-Mercado, et al.55 modeled asphaltene molecules as “discotic seven-centered LennardJonesium molecules” and resins as spheres and then performed canonical Monte Carlo simulations to probe the effects of temperature and concentration on the aggregate cluster size. However, the dynamic prediction of asphaltene aggregation and breakage under shear conditions using MC simulation technique has not been reported before in the open literature. In the present study, the general PBE is solved to predict the dynamic evolution of asphaltene ASD using a time-driven MC algorithm. The predicted results are compared with the experimental data for a shear-induced asphaltene solution undergoing particle aggregation and breakage for the first time. Asphaltene molecules are considered as fractal aggregates.
where, β(Vi, Vj) in the above two equations represents the aggregation kernel. If the value of a generated random 0 , then particles number R1 is smaller than the value of Pi,j i and j can be aggregated together. Aggregation is characterized by the aggregation kernel, βi,j, which expresses the aggregation rate between i and j particles of volumes Vi and Vj. There are three main mechanisms presented in literature for aggregation of particles. These mechanisms are perikinetic aggregation, orthokinetic aggregation, and aggregation due to differential settling.1,50,59 The first one is due to Brownian motion (thermal motion). The orthokinetic one occurs due to fluid shear. Here, only particles of at least 1 μm are considered where the influence of thermal motion is negligible in comparison with shear effect.59 The aggregation due to differential settling happens when particles settle with different settling velocities. It is important in systems only in which particles are allowed to settle. In this work, orthokinetic aggregation is considered for aggregation of particles due to shear-induced condition presented in the model. The aggregation kernel from the rectilinear model for orthokinetic aggregation can be written in the form of eq 5:60
Model Description
where, G is the shear rate. Equation 5 only applies to rigid spheres. Although like asphaltene aggregates the aggregated particles are not rigid spheres, the predicted results from this relation are not in general in good agreement with the experimental results. For a fractal aggregate with a mass equivalent volume of Vi, the collision diameter dc,i, is related to the number of primary particles xi of diameter dP in it and is defined by:46
1=3
βi, j ¼ 0:31GðVi
The time-driven MC simulation is performed to predict of the dynamic evolution of asphaltene ASD. This method was developed by Liffman57 for coagulation and extended by Zhao58 to include the breakage mechanism. Once all the particles in the sample population are assigned to selected volumes according to the defined initial size distribution, the MC algorithm is initiated, and the effects of particle aggregation and breakage mechanisms on the dynamic evolution of the particle population are stochastically simulated. The simulation has followed the principal steps of the following procedures. Aggregation Event. To evaluate the occurrence possibility of the aggregation event of particle i, during the discrete time interval (simulation time step), Δt, a random number, R1, should be generated from a uniform distribution between 0 and 1. Particle i is under aggregation if:58 " # -Δt ð2Þ R1 e 1 - exp 2tagg, i
VS βðVi , Vj Þ
1=df
1=3
þ xj Þ3 ð5Þ
ð6Þ
where df denotes the fractal dimension of the aggregates. By considering this fact that asphaltene aggregates have a fractal-like structure, aggregation kernel can be rewritten as:59 1=df
βi, j ¼ 0:31GVP ðxi
1=df 3
þ xj
Þ
ð7Þ
When the aggregation event occurs between two particles of volumes Vi and Vj, an aggregate of volume Vk is produced and: ð8Þ Vk ¼ Vi þ Vj Breakage Event. A breakage event occurs if the following relation is satisfied:42,43,58 R2 eSi Δt
ð3Þ
ð9Þ
where R2 and Si are a generated random number from the uniform distribution between 0 and 1 and the breakage kernel, respectively. The breakage kernel, Si, is the rate of particle breakage and is related to the shear rate and aggregate size by:61
j
where, VS is the volume of simulation box. The aggregation partner i, particle j, is chosen with a probability:58 βðVi , Vj Þ P0 i, j ¼ P βðVi , Vj Þ
1=3
dc, i ¼ dP xi
In eq 2, tagg,i is the time required for aggregation event, which can be estimated as:58 tagg, i ¼ P
1=3
þ Vj Þ3 ¼ 0:31GVP ðxi
ð4Þ
Si ¼ A0 Gq Vi
1=3
ð10Þ
j
where q is a constant that is inversely proportional to the aggregate strength, and A0 is a proportionality constant that
(55) Aguilera-Mercado, B.; Herdes, C.; Murgich, J.; Muller, E. A. Energy Fuels 2006, 20, 327–338. (56) Barcenas, M.; Duda, Y. Phys. Lett. A 2007, 454–457. (57) Liffman, K. J. Comput. Phys. 1992, 100, 116–127. (58) Zhao, H.; Maisels, A.; Matsoukas, T.; Zheng, C. Powder Technol. 2007, 173, 38–50.
(59) Barthelmes, G.; Pratsinis, S. E.; Buggisch, H. Chem. Eng. Sci. 2003, 58, 2893–2902. (60) Saffman, P.; Turner, J. J. Fluid Mech. 1956, 1, 16–30. (61) Pandya, J. D.; Spielman, L. A. Chem. Eng. Sci. 1983, 38, 1983–1992.
4954
Energy Fuels 2010, 24, 4952–4960
: DOI:10.1021/ef100123x
Faraji and Solaimany Nazar
is determined experimentally. The breakage parameters A0 and q for each set of experiments are constant. These constants are adjusted to match the measured results. The constant A0 affects the relative strength of breakage relative to aggregation. The proportionality of the breakage kernel to the aggregate size have also been used by other authors62,63 and express the fact that aggregates with a high diameter are easier to break in sheared flow. By replacing di with the collision diameter, dc,i, in eq 10, the breakage kernel can be expressed for a fractal nature of asphaltene aggregates as:59 1=df 0 q 1=3 dc, i Si ¼ A G VP ð11Þ dP
Table 1. Toluene to Heptane Volume Ratio, Shear Rate, and Asphaltene Solid Volume Fraction of Samples in Each Run52 shear rate range, G (s-1)
asphaltene solid volume fraction (jS)
1:3 1:4 1:9
3.75-9.38 3.75-9.38 3.75-9.38
6.38 10-5 7.28 10-5 8.09 10-5
time step is set to an interval smaller than the fastest process. The sample refreshment technique allows the simulation to proceed without losing the statistical accuracy of the initial calculation. It is noted that in the case of aggregation, the doubling of the sample creates a new sample with identical statistics, thus avoiding the noise introduced by the random selection of the particles to fill the box, as is the case for the constant number MC.58 The above claims are tested thoroughly in a series of MC runs from different starting points during algorithm development. These tests revealed a significant repeatability of the distribution results, all with reasonable accuracy as the outcome of a single simulation run. Simulated Case Study. There are a few experimental data available dealing with the asphaltene ASD. Asphaltene solution in shear induced suspensions is considered in the simulated case study. The experimental data of the evolution of the asphaltene ASD recently presented by Solaimany Nazar and Rahimi 52 is used to validate the prediction of the presented model. In their work, asphaltenes were extracted from a crude sample of a well located in Bangestan petroleum field of Iran. The image-processing technique was used to study the evolution of asphaltene ASD in shear induced suspensions.52 Diluted solutions of asphaltene in solvent mixtures composed of toluene solvent and n-heptane precipitant were prepared to avoid the application limitation of the image processing technique, which is unsuitable for optically opaque concentrated asphaltene solutions. nHeptane solvent is a precipitant agent, causing asphaltene precipitation and aggregation in solution. Asphaltene powder was dissolved in toluene solvent to achieve a stock solution with concentration of 1 g/L. Then the samples were provided adding distinct volumes of n-heptane to toluene-asphaltene solution. The samples were prepared in different values of toluene to n-heptane volume ratio (T/H). The wait time between preparation of diluted asphaltene solutions in heptane/toluene and the beginning of the microscopy experiment was ∼12 h. The solutions aged sufficiently to reach a steady aggregate size. The shear rate was induced on the samples previously poured in an annular cell of an experimental apparatus while the inner wall of the cell rotated mechanically. The values of T/H volume ratio, shear rate range, and asphaltene solid particles volume fraction according to T/H volume ratio used in performing the experiments are presented in Table 1. The number of asphaltene particles used in the simulations is set to 100. The size of the simulation box is related to the asphaltene solid particles volume fraction, js, and the known initial distribution and number of particles. The relation between the volume of the simulation box and the asphaltene solid particles volume fraction is expressed as:
In this work, for convenience and simplicity of calculation, it is assumed that the binary breakage event assumes the breaking of an aggregate of volume Vi to be into two equal fragments of volume Vi/2.62 The binary breakage is the most widely used model for understanding the fragmentation of colloidal suspensions.52,59 For occurrence of breakage event, the time can be computed through:42,43,58 1 tbrk, i ¼ ð12Þ Si MC Methodology. In time-driven MC, the time increment, Δt, is fixed, and every simulation particle is tracked to determine whether it participates in an event within that time interval. MC simulation starts with a collection of asphaltene molecules, which are randomly positioned in the simulation box of volume VS. In general, the initial distribution of particles may be randomly generated or set if it has a specified distribution. In each time step of the simulation, each particle is examined to determine whether it participates in an either aggregation or breakage event or in neither aggregation nor breakage. The number of particles in the simulation box needs to be restored to its initial value with periodic regulation of the number of particles.58 If the number of particles during simulation increases, it is halved once it reaches the double value of the initial number of particles, and if the number in the simulation box decreases, then it is doubled when it reaches half of the initial value. Consequently, particle distribution does not change during the particle refreshing routine and ensures that the statistical accuracy of the simulation is not gradually lost. The simulation time step, Δt, should be less than or equal to both the minimum aggregation time and minimum breakage time:42,43,58 Δt e minfðtagg, i Þ, ðtbrk, i Þg
volume ratio of toluene to heptane (T/H)
ð13Þ
The computational precision of MC simulation depends on the number of MC cycles, so the time interval may be written as: ð14Þ Δt ¼ Rminfðtagg, i Þ, ðtbrk, i Þg where R, is a small constant, typically R e 0.01. By setting the time step on smaller than the fastest process in the simulation, ensures that at most one event occurs within time, Δt, for each particle in the simulation box. The stochastic computational precision depends on the number of MC loops and the total number of simulation particles. To increase the number of Monte Carlo loops, the
n Π P di3 i¼1 6 js ¼ VS
ð15Þ
where n is the total number of particles in the simulation box, and di is the mass equivalent rigid sphere aggregate diameter.
(62) Spicer, P. T.; Pratsinis, S. E. AlChE J. 1996, 42, 1612–1620. (63) Williams, M. M. R. Aerosol Sci. Technol. 1990, 12, 538–550.
4955
Energy Fuels 2010, 24, 4952–4960
: DOI:10.1021/ef100123x
Faraji and Solaimany Nazar
Figure 2. Image of the initial asphaltene particles for G = 5.63 (1/s) and T/H = 1:4.64
Figure 1. Asphaltene number average diameter evolution under different shear rates for T/H=1:4.
For an aggregate with collision diameter of dc,i, the mass equivalent rigid sphere aggregate diameter, di, may be computed from:52 ð1=3 - 1=df Þ
di ¼ dc, i xi
ð16Þ
Thus, the mass equivalent rigid sphere aggregate diameter is given by: ðdf =3Þ - 1 dc, i ð17Þ di ¼ dc, i dP The number average diameter of aggregates is defined by: P ni dc, i ð18Þ dav ¼ P ni
Figure 3. Initial asphaltene size distribution obtained from the image presented in Figure 2.
gates diameter due to the balance of the aggregation and breakage rates. Figure 2, shows an image of the initial asphaltene particles captured under G = 5.63 (1/s) and T/H = 1:4.64 Figure 3, presents the determined initial particle size distribution with an initial number average diameter of 110 μm, obtained from image analysis using Sigma Scan Pro 5 software on the image shown in Figure 2. The relative distribution of aggregates defined in this figure is the ratio of the number of aggregates with a certain diameter to the total number of aggregates. A simple relation between the particle number and the particle number density can be written as: Z vi þ 1 n dv ð19Þ Ni ¼
The initial distribution of asphaltene particles are determined by the first image captured at the initial time of measurement. The initial number average diameter was obtained from this image.52 The fractal structure dimension of the asphaltene sample was reported 1.6.52 The value of q= 2 was used in all calculations corresponding to asphaltene aggregates.52 The value of A0 = 0.00388 (m-1 s) is used for the simulated case study.52 We did not attempt to adjust the new values of parameters q and A0 to the model. Results and Discussions Effect of Shear Rate on ASD. The effect of shear rate on the predicted results of the asphaltene aggregates number average diameter is depicted in Figure 1. The initial asphaltene aggregate number average diameter is considered 130 μm. As expected, the model is very sensitive to the shear rate value due to the dependence of the aggregation and breakage kernels on this parameter. By inducing shear rate, aggregation and breakage processes are started simultaneously. At the earliest time, the aggregation phenomenon is predominant and the aggregates number average diameter grows with time and reaches a maximum value. At this time, the aggregates are large enough for their breakage probability to be higher than the aggregation probability. Therefore, the breakage phenomenon predominates. The aggregate number average diameter starts to decline beyond the maximum point. After some time, due to the increase in the population of smaller aggregates, aggregation mechanism is predominant again. Finally, the size of aggregates approaches a steady state size and no changes occur in the size of aggre-
vi
Figure 4 is a comparison of the calculated and experimental data52 for T/H = 1:4, under two shear rates of 3.75 and 7.51 (1/s). As shown, the variations of the number average diameter are well predicted by the model. In addition, as depicted in this figure, the variations of the results have qualitatively similar trends of the results in Figure 1. In Figure 5 the MC results for the relative asphaltene ASD is compared with those previously obtained by a deterministic solution of PBE and experimental data52 to study its accuracy at various aggregation time stages with G = 5.63 (1/s) and T/H=1:9. From the results of Figure 5, one can (64) Rahimi, H. Application of Population Balance Technique in Determination of Asphaltene Aggregate Size Distribution; M.S. Thesis, University of Isfahan: Isfahan, Iran, 2008.
4956
Energy Fuels 2010, 24, 4952–4960
: DOI:10.1021/ef100123x
Faraji and Solaimany Nazar
conclude that both approaches can be employed to calculate the ASD of asphaltene and qualitatively, the agreement is good on the whole. Table 2 compares the average percent deviation values of the predicted average asphaltene aggregate diameter from the experimental results at various aggregation time stages with G = 5.63 (1/s) and T/H = 1:9, obtained using both computational approaches results. Clearly, the MC simulation is in much better agreement with experimental data than is the deterministic solution of PB results. Figure 6, compares the cumulative asphaltene ASD results of the MC method under the deterministic solutions of PBE at G = 5.63 (1/s) and T/H = 1:9 conditions. In this figure, the vertical axis presents the ratio of the number of aggregates as smaller or equal to a certain diameter to the total number of aggregates. As shown, the MC simulation can provide relatively accurate predictions compared with the experimental Table 2. Average Percent Deviation Values of Monte Carlo and Deterministic Approaches for Prediction of ASD at Various Stages of Aggregation from Experimental Results at G = 5.63 (1/s) and T/H = 1:9
Figure 4. Comparison of the predicted asphaltene number average diameter evolution and experimental data for T/H = 1:4 under different shear rates; (a) for G = 3.75 (1/s), (b) for G = 7.514 (1/s).
time (min)
average percent deviation of Monte Carlo simulation from experimental results
average percent deviation of deterministic predictions from experimental results
4 15 25 40
58.80 28.17 17.65 41.97
97.07 35.09 114.67 57.25
Figure 5. Asphaltene ASD at various stages of aggregation for G = 5.63 (1/s) and T/H = 1:9.
4957
Energy Fuels 2010, 24, 4952–4960
: DOI:10.1021/ef100123x
Faraji and Solaimany Nazar
Figure 6. Cumulative asphaltene ASD at various stages of aggregation for G = 5.63 (1/s) and T/H = 1:9.
Figure 8. Predicted steady state asphaltene aggregate number average diameter vs shear rate for T/H = 1:4.
Figure 7. Predicted required time for aggregates to reach a maximum number average diameter and predicted values of maximum number average diameter vs shear rate for T/H = 1:4 and dav0 = 130 μm.
presented for T/H = 1:4 condition. The initial average diameter of aggregates is taken to equal 130 μm. Clearly, increase in the shear rate results in a smaller maximum average diameter at shorter times. Further results concerning the values of the steady state number average diameter with respect to shear rate for T/H = 1:4 condition are indicated in Figure 8. In fact, the steady average aggregate size decreases during progression of the breakage process with an increase in shear rate. Effect of Toluene to n-Heptane Volume Ratio. The T/H volume ratio dictates the asphaltene solid volume fraction as shown in Table 1. Due to the good dissolution of asphaltene in toluene solvent, the asphaltene solubility is proportional
data. Although, the MC model prediction is superior to the PBE method, there are significant deviations in the MC prediction from the experimental data in Figures 4-6, especially in a short time. Therefore, it seems the model still has room for improvement. A more realistic extension of the model would be to include the restructuring phenomenon of aggregates and the variation of their fractal dimension. Figure 7 illustrates the predicted required time and the values of maximum number average diameter for the aggregates to reach a maximum number average diameter as a function of shear rate respectively. The predicted results are 4958
Energy Fuels 2010, 24, 4952–4960
: DOI:10.1021/ef100123x
Faraji and Solaimany Nazar
Figure 11. Asphaltene number average diameter evolution for G = 9.38 (1/s) and different asphaltene solid volume fraction.
certain time, then it begins to decrease with time, and the diameter approaches a nearly constant value. However, for the bigger diameter the breakage mechanism is predominant at the earliest time. Hence, the initial number average diameter is a key parameter in the MC model for predicting the dynamic evolution of the aggregate size distribution, but this parameter is less important for predicting the steady state value. Effect of Initial Asphaltene Concentration. Figure 11 shows the time evolution of the asphaltene number average diameter at different asphaltene solid volume fractions for G = 9.38 (1/s). Here, in all the simulation runs the initial number average diameter of aggregate size is 100 μm to show the effect of initial asphaltene concentrations separately. As expected, the higher solid volume fractions, the faster growth of average diameter causes, and therefore a higher maximum average diameter is achieved. Moreover, the higher js leads to further reduction of the required time to achieve the maximum size.
Figure 9. Comparison of the predicted asphaltene number average diameter evolution and experimental data under G = 5.63 (1/s) for (a) T/H = 1:3 and (b) T/H = 1:4.
Conclusions For the first time, the time evolution of asphaltene ASD in simultaneous aggregation and breakage processes in shearinduced petroleum suspensions was studied using stochastic MC simulation. Asphaltene aggregates are considered as a fractal structure, and the effects of shear rate, initial asphaltene concentration, and initial number average diameter on the evolution of asphaltene ASD are studied. It was demonstrated that a balance between aggregation and breakage phenomena dictated the evolution trend during the early stages. The existing maximum point in the evolution trend depends on the magnitude of the fragmentation kernel, the shear rate, and the initial particle size distribution. However, the average cluster size declines to a steady value when the breakage is predominated. An increase in the solid asphaltene concentration results in a faster growth of aggregates; however, the steady state size of aggregates is insignificantly affected. The numerical results obtained from using MC simulations are relatively consistent with the experimental data available.
Figure 10. Evolution of number average diameter for G = 5.63 (1/s) and T/H = 1:4 with different initial average aggregate size.
to the T/H volume ratio. Therefore, as a function of this ratio, the out-of-solution particles are lowered. Figure 9 shows the effect of T/H volume ratio on the asphaltene number average of diameter evolution corresponding to G = 5.63 (1/s) for T/H = 1:3 and T/H=1:4. The higher T/H ratio yields a higher number average diameter. Since with a change in T/H ratio from 1:3 to 1:4, the concentration of solid asphaltene changes only from 64 to 73 ppm, the number average diameter will vary slightly. The model is efficiently capable of predicting the experimental results. Effect of Initial Number Average Diameter. Figure 10 presents the effect of the initial number average diameter of asphaltene aggregates on the dynamic evolution of ASD. At the beginning of the simulation corresponding to the smaller initial diameter of aggregates, the aggregation event is predominant. Finally, the number average diameter reaches a maximum at
Abbreviations A = Breakage parameter (m-1 s) d = Mass equivalent rigid sphere aggregate diameter (m) dav = Number average diameter of aggregate (m) dc = Collision diameter of aggregate (m) df = Fractal dimension (-) 0
4959
Energy Fuels 2010, 24, 4952–4960
: DOI:10.1021/ef100123x
Faraji and Solaimany Nazar
-l
G = Shear rate (s ) N = Number of aggregates (-) n = Total number of particles (-) p = Probability function (-) q = Breakage parameter (-) R = Random number (-) S = Breakage kernel (m3 s-l) t = Time s VP = Particle volume (m3) x = Number of primary particles (-)
Subscripts agg = Aggregation brk = Breakage i,j,k = Index 0 = Initial condition Acronyms ASD = Aggregate Size Distribution PBE = Population Balance Equation MC = Monte Carlo T/H =Toluene to n-heptane volume ratio
Greek Letters R = Constant β = Aggregation kernel (m3 s-l) Γ = Breakup distribution function Δ = Discrete interval jS = Solid volume fraction
Acknowledgment. The authors would like to acknowledge Isfahan Oil Refining Company for the financial support of the present study.
4960
