Aggregation in Colloidal Suspensions: Evaluation of the Role of

Article pubs.acs.org/JPCB

Aggregation in Colloidal Suspensions: Evaluation of the Role of Hydrodynamic Interactions by Means of Numerical Simulations A. Tomilov,† A. Videcoq,*,† M. Cerbelaud,‡ M. A. Piechowiak,† T. Chartier,† T. Ala-Nissila,§,∥ D. Bochicchio,⊥ and R. Ferrando⊥ †

SPCTS, UMR 7315, ENSCI, CNRS; Centre Européen de la Céramique, 12 rue Atlantis, 87068 Limoges cedex, France Institut des Matériaux Jean Rouxel (IMN), Université de Nantes, CNRS, 2 rue de la Houssinière, BP32229, 44322 Nantes cedex 3, France § Department of Applied Physics and COMP CoE, Aalto University School of Science, P.O. Box 11000, FI-00076 Aalto, Espoo, Finland ∥ Department of Physics, Brown University, Providence, Rhode Island 02912-8143, United States ⊥ Dipartimento di Fisica and CNR-IMEM, Via Dodecaneso 33, Genova I-16146, Italy ‡

S Supporting Information *

ABSTRACT: Numerical simulations constitute a precious tool for understanding the role of key parameters influencing the colloidal arrangement in suspensions, which is crucial for many applications. The present paper investigates numerically the role of hydrodynamic interactions on the aggregation processes in colloidal suspensions. Three simulation techniques are used: Brownian dynamics without hydrodynamic interactions, Brownian dynamics including some of the hydrodynamic interactions, using the Yamakawa−Rotne− Prager tensor, and stochastic rotation dynamics coupled with molecular dynamics. A system of monodisperse colloids strongly interacting through a generalized Lennard-Jones potential is studied for a colloid volume fraction ranging from 2.5 to 20%. Interestingly, effects of the hydrodynamic interactions are shown in the details of the aggregation processes. It is observed that the hydrodynamic interactions slow down the aggregation kinetics in the initial nucleation stage, while they speed up the next cluster coalescence stage. It is shown that the latter is due to an enhanced cluster diffusion in the simulations including hydrodynamic interactions. The higher the colloid volume fraction, the more pronounced the effects on the aggregation kinetics. It is also observed that hydrodynamic interactions slow down the reorganization kinetics. It turns out that the Brownian dynamics technique using the Yamakawa−Rotne−Prager tensor tends to overestimate the effects on cluster diffusion and cluster reorganization, even if it can be a method of choice for very dilute suspensions.

1. INTRODUCTION

material, such as the mechanical, the electrical, and the optical ones. Given the high number of parameters that influence the suspension behavior, numerical simulations constitute a precious tool to understand the role of each parameter. Many numerical studies concerning aggregation phenomena in colloidal suspensions have been performed using the Brownian dynamics (BD) technique, neglecting the many-body hydrodynamic interactions (HIs) between the colloidal particles.1−3 In several cases, this simulation technique is able to predict the structural properties of dilute suspensions in agreement with experimental observations. An example is the heteroaggregation of dilute suspensions of alumina and silica ceramic colloids4,5 that has been studied both experimentally and by BD, obtaining

Understanding the behavior of colloidal suspensions is of great importance in various technological processes. An important example is the field of ceramics. In fact, several ceramic processes commonly use the colloidal route. This is the case of casting or casting−coagulation processes, to obtain ceramic parts with specific geometry (for example, hip prosthesis). The colloidal route is also relevant to tape-casting processes that are employed, for example, to produce fuel cell electrodes. Other processes include stereolithography and inkjet printing for complex geometries. In order to improve ceramic processes or to develop new and reliable processes, it is crucial to understand in detail the behavior of colloidal suspensions, with special attention to their structural and rheological properties. Here, we focus on the first aspect, i.e., the arrangement of particles in the solvent, which is produced by complex aggregation phenomena, and eventually determine the final microstructure and the properties of the © 2013 American Chemical Society

Received: July 22, 2013 Revised: September 24, 2013 Published: October 21, 2013 14509

dx.doi.org/10.1021/jp407247y | J. Phys. Chem. B 2013, 117, 14509−14517

The Journal of Physical Chemistry B

Article

types of simulations could be complementary. Depending on the conditions, it could be more appropriate to use either BDYRP or SRD-MD. For example, the study of dilute suspensions with simulations where the solvent is modeled by particles could be very cumbersome compared to the study performed with an implicit description of the solvent. In this paper, we will thus compare the aggregation phenomena observed in BD (without and with HIs) with those observed in SRD-MD. The system used is composed of only one kind of particles interacting through a generalized Lennard-Jones potential with a deep well, whose width is chosen to reproduce the typical width of the alumina−silica interactions of our previous studies.5,14,18 The particle diameter is 600 nm, and the particles have a mass density of silica, as in our previous studies. Simulations are performed within the range of the colloid volume fraction, ϕc, in between 2.5 and 20%. The paper is structured as follows. In section 2, the different simulation techniques are briefly described. Section 3 is divided into five subsections, dealing with the simulation of the whole aggregation process by the three computational techniques. Here, we consider separately the initial stage of nucleation that is analyzed by monitoring dimer formation and the following stage in which small clusters diffuse, coalesce, and reorganize. Section 4 contains the summary and conclusions.

a good agreement for the overall shape and local order of the aggregates. Nevertheless, a natural question that arises concerns the influence of the HIs on the aggregation phenomena, especially for suspensions that are not dilute, in which HIs are expected to play a more important role. The present paper addresses this question by means of numerical simulations, comparing the results of BD with those of more sophisticated techniques that include HIs in different ways. In fact, several techniques designed to correctly reproduce the hydrodynamics in colloidal suspensions have emerged over the past decades. Some of these techniques are based on an improvement of the BD technique, in which HIs are included through a diffusion tensor.6−8 In other techniques, the solvent is taken into account explicitly, by effective fluid particles, as in dissipative particle dynamics or in stochastic rotation dynamics-molecular dynamics (SRD-MD). Some studies of the aggregation including HIs have recently been reported in the literature. In order to investigate the role of many-body HIs, Furukawa and Tanaka have used the fluid particle dynamics method,9 Whitmer and Luijten have used the hybrid stochastic rotation dynamics-molecular dynamics (SRDMD) technique,10 and Cao et al. the accelerated Stokesian dynamics.11 It turns out that significant hydrodynamic effects on cluster formation and gelation are observed for weak colloidal attraction, i.e., when the depth of the potential well is less than approximately 5 kBT.9,10 For a potential well depth of 8 kBT, it is still found that the systems with HIs form a percolated network at lower colloid volume fraction, and HIs slow down the reorganization of trimers from linear configuration to the minimum energy, triangular configuration.11 In the present paper, we study in detail the role of HIs in the different stages of the aggregation process. In fact, in typical colloidal aggregation processes, such as those studied in refs 5 and 12−14, there is an initial nucleation stage in which several small clusters form. This is followed by a coalescence stage, in which these small clusters diffuse in the solvent, collide with each other, and join to form larger aggregates. Our results will show that HIs have quite distinct effects on nucleation and coalescence stages. In several systems of interest for ceramic processing, colloids attract each other quite strongly. In our previous works, we have studied alumina−silica colloidal systems by means of BD simulations (without HIs),5,12−14 and from a comparison between experimental and numerical results, we have evaluated the potential well depth for the attraction between alumina and silica particles to be 14 kBT at room temperature.4,15 A further scope of the present study is therefore to understand how HIs affect the aggregation process in such a case of strong colloidal attraction, in which we expect cluster reorganization to be more difficult than in the cases considered in refs 9−11. In the following, we use three different simulation techniques: BD without HIs, BD with HIs, using the Yamakawa−Rotne−Prager tensor (BD-YRP),16,17 and SRDMD. The comparison of these techniques is interesting for understanding the effects of HIs on aggregation as a function of the colloid volume fraction. Indeed, the BD-YRP technique partially includes HIs, since it is valid only for relatively dilute suspensions and lubrication forces at short range are not included in this technique. Moreover, the computational cost of these simulation techniques is very different. In the SRD-MD technique, the solvent is treated explicitly by fluid particles, while, in BD-YRP, it is taken into account implicitly. These two

2. SIMULATION METHODOLOGY SRD-MD simulations are performed within the same framework as in our previous study.18 The solvent is simulated by SRD, with the three independent parameters chosen as suggested by Padding and Louis:19 the dimensionless mean free path of fluid particles, λ = 0.1, the average number of fluid particles per collision cell, γ = 5, and the rotation angle used in the collision steps, α = 90. In between two SRD collision steps, the colloids and the fluid particles evolve by molecular dynamics with colloid−colloid and colloid−fluid interaction potentials, while the fluid particles interact together during the collision steps only. The colloid−fluid interaction is described by a repulsive, short-range potential: ⎧ ⎛ σ ⎞12 ⎪ εcf ⎜ cf ⎟ (r ≤ rc ≡ 2.5σcf ) Vcf (r ) = ⎨ ⎝ r ⎠ ⎪ (r > rc) ⎩0

(1)

where the colloid−fluid interaction parameter, σcf, is chosen smaller than the colloid radius, ac (here ac = 300 nm). This allows the fluid particles to penetrate in between two colloids, so that spurious depletion attraction between colloids is avoided.19 This also mimics lubrication forces between the colloids. A value of σcf = 0.8ac is used here, and εcf = 2.5kBT. The different scales used in the SRD-MD simulations are the same as in ref 18 and are as follows. The linear size of the collision cells, a0, which fixes the length scale, is chosen as large as possible but small enough to correctly describe the hydrodynamics at the scale of the colloids: a0 = ac/2. The mass of the fluid particles, mf, i.e., the mass scale, is fixed in order to reproduce the mass density of water (ρf = 1000 kg m−3): mf = a03ρf/γ = 6.75 × 10−19 kg. The mass of a colloid corresponds to the mass of a silica sphere (mass density of ρc = 2200 kg m−3): Mc = 4πac3ρc/3 = 2.49 × 10−16 kg. The time scale of the SRD-MD simulation is adjusted in order to correctly reproduce the diffusion time. This is the time needed 14510



Article

In the following, “BD” refers to simulations of Brownian dynamics without hydrodynamics and “BD-YRP” to simulations of Brownian dynamics where hydrodynamics is taken into account through the YRP tensor. In all the simulations, the colloid−colloid, attractive potential is a generalized Lennard-Jones potential:

for a colloid to diffuse over a distance equal to its radius: τD = ac2/D0, where D0 is the colloid diffusion coefficient in infinite dilution. It is related to the friction coefficient, ζ, by D0 = kBT/ ζ. In the real physical system, the friction coefficient (assuming stick boundary conditions at the colloid surface) is 6πηac so that, with T = 293 K and η = 10−3 Pa s (viscosity of water), we obtain τD = 0.126 s. In the SRD-MD simulations, the friction coefficient generated by the simulated fluid can be estimated as a function of the SRD parameters (λ, γ, α), a0, mf, and t0, which denotes the time scale. Thus, the diffusion time can be expressed as a function of t0, and by equating to the real physical value of τD, the time scale t0 can be deduced. The complete derivation of this value is detailed in ref 18, and it provides t0 = 7.37 × 10−4 s. It is to be noticed that the SRD time step is ΔtSRD = λt0 = 7.37 × 10−5 s. As the length, mass, and time scales are now defined, the energy scale, i.e., kBT, is determined as well: kBT = a02mf/t02. This corresponds to a temperature of 2.02 × 10−3 K in the SRD-MD simulation. In the BD simulations, the Langevin equation of motion is integrated for each colloid.20 The effect of the solvent on the colloids is represented by the hydrodynamic drag, with a constant friction coefficient of 6πηac (with η = 10−3 Pa s), and by random forces, which represent the Brownian motion. Here, the temperature is fixed at T = 293 K. As already mentioned in the Introduction, hydrodynamic interactions can be introduced in Brownian dynamics simulations by using a diffusion tensor, which represents the diffusive or frictional effects occurring in the system. To better understand the effect of this tensor on the aggregation, we also performed Brownian dynamics simulations where we introduce the Yamakawa−Rotne−Prager (YRP) diffusion tensor, which is defined for spheres of equal dimension (radius ac) as follows:16,17 • when the distance between particles i and j is higher than their diameter, rij > 2ac: Dij = δij

⎧ ⎡ ⎛ σ ⎞ 2n ⎛ σ ⎞ n ⎤ cc cc ⎪ ⎪ 4εcc⎢⎜ ⎟ − ⎜ ⎟ ⎥ (r ≤ rc ≡ 1.6σcc) ⎝ ⎠ ⎝ ⎠⎦ r r ⎣ ⎨ Vcc(r ) = ⎪ ⎪0 (r > rc) ⎩

with σcc = 2ac. If not mentioned otherwise, n = 18 and εcc = 14kBT.

3. RESULTS AND DISCUSSION 3.1. Dimer Dissociation Time. As explained in ref 19, SRD-MD is compacting the hierarchy of time scales. In order to compare the aggregation processes between BD and SRD-MD simulations, there is an important point that has to be considered first. This point consists of checking whether the two following time scales are of the same order of magnitude in both simulation types: the diffusion time and the dimer dissociation time, i.e., the average time needed for a dimer to dissociate. According to our choice of the time scale in the SRD-MD technique, we know that the diffusion time is the same in both simulations. Therefore, we have to check whether the dimer dissociation time is also of the same magnitude. We have performed BD and SRD-MD simulations with the colloid−colloid potential of eq 4, with n = 6. The potential well depth has been chosen at εcc = 7kBT, so that the dissociation time is not too long and rather quickly reachable by the SRDMD technique, which is computationally demanding. The simulations are performed with two particles, initially placed at a distance corresponding to the minimum energy, and the system is left evolving until the two particles dissociate, i.e., until they escape from the potential well. The time it takes is registered and the dimer dissociation time is deduced from an average over more than 100 simulations of this kind. The same value of 11 s is obtained in both SRD-MD and BD techniques. Dimer dissociation is therefore taking place on similar time scales. As both the diffusion time and the dimer dissociation time are similar in SRD-MD and BD, we think that it is relevant to compare the aggregation kinetics in the two techniques to understand if this kinetics is affected by HIs. 3.2. Aggregation Kinetics. Here the aim is to study the aggregation kinetics at different colloid volume fractions in the simulations with and without the HIs. The aggregation kinetics is typically followed by the number of aggregates as a function of time. In order to have the same statistics, we have performed the simulations with a constant number of colloids (Nc = 500), varying the size of the simulation box to get different colloid volume fractions. The linear size of the simulation box is 69a0, 55a0, 48a0, and 44a0 for colloid volume fractions of approximately 5, 10, 15, and 20%. Initial configurations are generated with the SRD-MD technique. The particles are randomly placed at the beginning, after which the SRD-MD technique needs some equilibration steps in order to stabilize the temperature at the desired value. These steps are performed with only the repulsive part of eq 4, which produces better dispersed configurations than the initial random ones. In order to compare the aggregation kinetics with different simulation techniques, it is important to start with initial configurations of

rijr Tij ⎞ kBT kT ⎛ kT I + B ⎜⎜I + 2 ⎟⎟ − B 3 (2ac 2) 6πηac 8πηrij ⎝ rij ⎠ 8πηrij

⎛ r rT ⎞ ⎜ ij ij − 1 I⎟ ⎜ r2 3 ⎟⎠ ⎝ ij

(2)

• when the distance between particles i and j is lower than their diameter, rij ≤ 2ac: Dij =

rij rijr Tij ⎞ 3rij ⎞ kBT ⎛⎜⎛ 4ac ⎟ I − + ⎜ ⎟ 8 ⎠ 8 rij 2 ⎟⎠ 8πηac 2 ⎜⎝⎝ 3

(4)

(3)

where δij is the Kronecker symbol and I is the identity tensor. One of the most important approximations used to derive this tensor is to consider that hydrodynamic interactions are pairwise additive. Therefore, this tensor is generally appropriate for systems with relatively low colloid volume fractions. In these simulations, even when periodic boundary conditions are applied, we do not use the minimum image convention to calculate the diffusion tensor. This is done in order to guarantee its positiveness, as shown by Jardat et al.17 Integration of the equation of motion is done by the Ermark et al. algorithm.6,21 The random part of the displacement is obtained by a Cholesky decomposition of the YRP tensor, which is performed by a Lapack subroutine. 14511



Article

when the number of aggregates increases, the dimer formation predominates. In stage ii, when the number of aggregates decreases, the cluster−cluster coalescence predominates. Obviously, the higher the colloid volume fraction, the faster the aggregation kinetics. At the end of all the simulations, all the colloids form a single large aggregate. Figure 2 compares the three simulation techniques for each colloid volume fraction. Here, an additional simulation result is shown for ϕc = 2.5%. It has been obtained with Nc = 250 and a simulation box size of 69a0 (the simulation with 500 colloids was too cumbersome within the SRD-MD technique for this volume fraction). We observe that, for each volume fraction, the time at which a single final cluster is formed is approximately the same in all the different simulation techniques. Nevertheless, if we compare the cases with HIs (BD-YRP and SRD-MD) and without (BD), it is worth noticing that the aggregation kinetics in stage i is somewhat faster without HIs than with HIs. More precisely, the kinetics of dimer formation orders as follows as a function of the simulation technique: BD > BD-YRP > SRDMD. On the contrary, the kinetics is faster with HIs than without HIs in stage ii. This means that HIs tend to slow down the dimer formation and tend to speed up the cluster−cluster coalescence. Very recently, similar tendencies have been reported in the simulation of lipid membrane formation by BD and BD-YRP.22 Even if these tendencies are present for all the colloid volume fractions tested, the difference between the simulation techniques is more pronounced at high volume fractions, confirming that hydrodynamics plays a more important role when the concentration is increased. The slowing down of dimer formation in the case with HIs is thus compensated by the speedup of coalescence, so that the final time for obtaining a single aggregate is approximately the same with and without HIs. For stage i, dimer formation, as the diffusion coefficient of isolated colloids is the same in all the simulation techniques, the retardation effects observed with HIs are likely to be due to repulsive HIs between the approaching isolated colloids. For stage ii, in which HIs seem to accelerate the cluster−cluster coalescence, one can deduce that cluster diffusion should be influenced by HIs. We address these two points in the next two sections. 3.3. Dimer Formation. Here we consider dimer formation from the point at which the two approaching colloids are sufficiently close to feel the attractive Lennard-Jones potential so that they are likely to enter into contact with high probability. The kinetics of contact formation is followed within the different simulation techniques. Simulations with only two colloids are performed. They are initially placed at a small enough (center-to-center) distance (750 nm). The distance between the two particles is monitored as a function of time. The result, averaged over five simulations, is shown in Figure 3. The distance between the two colloids decreases until it reaches approximately the interparticle distance at the potential well bottom (21/18σcc ≃ 624 nm). At that moment, the dimer is formed. The time needed for that is larger in SRD-MD than in the simulations of BD type. This difference can be explained by the lubrication effects. Indeed, these short-range HIs are not taken into account in BD-YRP, whereas the SRD-MD technique includes these interactions. HIs in SRD-MD simulations induce in fact correlated motions, which retard direct contact of colloids. Each time a colloid makes a step toward the other, it also pushes the fluid particles against the

an equivalent level of dispersion. That is the reason why the configurations produced by the SRD-MD equilibration steps are used as initial configurations in the three simulation techniques (BD, BD-YRP, and SRD-MD). Figure 1 shows the number of aggregates formed as a function of time. Here an aggregate is defined as an assembly of

Figure 1. Aggregation kinetics followed by the number of aggregates as a function of time, for the colloid volume fractions 5, 10, 15, and 20% and within the three simulation techniques (a) BD, (b) BD-YRP, and (c) SRD-MD.

at least two particles. The curves show the following two-stage evolution. At the beginning, the colloids are dispersed and the aggregation starts by the formation of dimers. Then, isolated particles aggregate on existing clusters and finally aggregate growth continues by the coalescence of aggregates. In stage i, 14512



Article

Figure 2. Comparison of the aggregation kinetics in between the different simulation techniques for (a) ϕc = 2.5%, (b) ϕc = 5%, (c) ϕc = 10%, (d) ϕc = 15%, and (e) ϕc = 20%.

3.4. Diffusion Coefficient of Small Clusters. In section 3.2, we have noticed that HIs accelerate the cluster−cluster coalescence. This may be due to a different diffusion behavior of the clusters when HIs are included. In order to compare the diffusion of small clusters in the three simulation techniques, we have computed the diffusion coefficient of clusters made of 2, 3, 4, 5, 13, and 55 colloids. The two last values correspond to magic numbers at which the clusters present an icosahedral structure, which is highly symmetrical and quasi-spherical. Indeed, the cluster of 13 particles is made of one shell of 12 particles surrounding the central one, and the cluster of 55 particles is a two-shell cluster. The diffusion coefficient is computed from the center-of-mass mean-square displacement.

other colloid, which slows down dimer formation. This is attested by the number of fluid particles that have to be removed in between the two colloids, when the dimer is forming. A movie showing the fluid particles in a sphere of the same size as the colloids and centered in between them is provided in the Supporting Information. About 30 particles are removed from that sphere during the dimer formation. According to this result, the slowing down of the kinetics in the dimer formation regime observed in the BD-YRP simulations is due to longer range repulsive HIs between approaching colloids. In SRD-MD, the retardation effect is amplified by the additional short-range lubrication effects. 14513



Article

The three different techniques exhibit quite different diffusion behaviors for the small clusters, and it is not obvious to determine which one is the most accurate, because of the lack of experimental data. Nevertheless, the diffusion coefficient of a spherical cluster might be roughly estimated from the diffusion coefficient of the smallest sphere containing the cluster. This sphere has a radius of about 3.1ac and 4.4ac for the clusters of 13 and 55 particles, respectively. The corresponding values of the normalized diffusion coefficients are shown in Figure 4. We may think that the actual diffusion coefficient of the clusters is probably smaller than this estimate because of the asperities of the raspberry aspect of the clusters, and because of the uncorrelated part of the motion of its constituents (which causes the string slowing down of diffusion in BD). Within the BD technique, the addition of the individual frictions of all particles belonging to the cluster, which is an artifact of this technique, certainly slows down the cluster diffusion too much compared to what sounds reasonable in reality. The use of the YRP tensor considerably increases the diffusion coefficient of the clusters. However, it can be seen that the diffusion of the clusters of 13 and 55 particles in BD-YRP is slightly larger than the one of the smallest spheres containing the clusters. This is particularly true for the cluster of 13 particles. This may indicate that the diffusion of clusters is somewhat overestimated within the BD-YRP technique. Indeed, BD-YRP assumes pairwise additivity of HIs, which is certainly not accurate in the condensed phases that constitute the clusters. The SRD-MD technique provides intermediate cluster diffusion coefficients between BD and BD-YRP. The values found for the spherical clusters of 13 and 55 particles are smaller than the diffusion coefficient of the smallest spheres containing these clusters, which seems reasonable. 3.5. Trimer Reorganization Kinetics. In previous works, we have shown that the structure of the aggregates formed in the simulations results from the competition between the kinetics of coalescence and the kinetics of reorganization of clusters into energetically favorable configurations.13,14 In the first part of this article, we have already seen that HIs accelerate the kinetics of coalescence. In this part, we are interested in understanding the role of HIs on the kinetics of cluster reorganization. To study this, we applied the method proposed by Cao et al., who has done a similar study using the accelerated Stokesian dynamics technique.11 This method will be quickly described in the following; however, its complete description can be found in ref 11. The configuration of threeparticle aggregates (trimers) is followed as a function of time. The structure of the trimers is characterized by the angle θ in between the vector connecting the centers of the first two particles and the one connecting the midpoint between these particles and the center of the third particle (see Figure 5). In this part, each simulation contains 64 trimers. At the beginning of the simulation (initial conditions), the particles in the trimers are aligned (θ = 0°) and the linear size of the simulation box is of 90a0, so that the surface-to-surface distance between neighboring trimers is at least of 5.09ac. Then, during the simulation, the trimers reorganize themselves progressively (θ < 90) to a triangular configuration (θ = 90°), which is the most energetically favorable configuration. The angle distribution as a function of time is obtained from an average over 30 simulations and is reported in Figure 6. Whatever the simulation technique, we can see that the probability to find trimers in a triangular configuration obviously increases with time. The mean number of trimers with an angle θ < 90°

Figure 3. Center-to-center separation distance between two approaching colloids during the final stage of dimer formation.

An infinite simulation box is used in the techniques of BD kind, while a finite size box is used in the SRD-MD technique because of the fluid particles. In the latter, a linear box size of 64a0 has been used. The cluster diffusion coefficient, normalized by the diffusion coefficient of an isolated colloid (D0), is shown in Figure 4.

Figure 4. Diffusion coefficients of clusters normalized by the diffusion coefficient of an isolated colloid, as a function of the number of colloids in the cluster, N. The lines are guides to the eye, corresponding to the following power laws: 1/N for BD, 1/N0.37 for BD-YRP, and 1/N0.59 for SRD-MD. The ∗ symbols indicate the diffusion coefficients of spheres of radii 3.1ac and 4.4ac (see the text for details).

In BD, the cluster diffusion coefficient is D0/N, if N denotes the number of particles in the cluster. In this technique, the effective friction coefficient of the cluster is the sum of the friction coefficients of all the particles: N × 6πηac. This has already been observed in a system with highly size-asymmetric particles, where the diffusion coefficient of a large particle (radius a1) surrounded by n small particles (radius a2) has been found equal to D0/(1 + na2/a1).15 In SRD-MD, the cluster diffusion coefficient decreases less dramatically as N increases, because the effective drag coefficient depends on the dynamically changing shape of the clusters. We find out that it behaves as D0/N0.59 at large N. The same tendency is observed in the BD-YRP simulation where the cluster diffusion coefficient behaves as D0/N0.37. This faster diffusion of the small clusters both in SRD-MD and BD-YRP compared to BD explains the faster aggregation kinetics in the second stage, in which the clusters grow by cluster−cluster coalescence. 14514



Article

somewhat slower in SRD-MD than in BD.23 Surprisingly, it appears that the reorganization kinetics is slightly slower in BDYRP than in SRD-MD, in which lubrication is taken into account. It is worth noting that, at t = 0.12 s, about 50% less triangular configurations are found in BD-YRP as compared to BD. For the same instant of time (t = τD) and almost the same colloid−colloid interaction potential, a similar result has been obtained by Cao et al. from the comparison between Stokesian dynamics, with lubrication forces included, and BD (Figure 12b of ref 11). Thus, the absence of the lubrication forces in our BD-YRP simulations does not explain the discrepancy between SRD-MD and BD-YRP. Instead, this could be due to the approximation used in the diffusion tensor (either in BD-YRP or in Stokesian dynamics), which supposes that HIs are pairwise additive. Such techniques do not take into account the collective effect of the trimer motion on the hydrodynamics. As for the case of aggregate diffusion considered in the previous section, it seems that the effects of HIs are overestimated by BD-YRP. Expressions of the diffusion tensor including threebody effects have already been developed; however, their use in simulations is too cumbersome.21

Figure 5. Definition of the angle θ characterizing the structure of a trimer: (a) θ < 90° (intermediate configuration) and (b) θ = 90° (most energetically favorable configuration).

decreases progressively in favor of the number of trimers with an angle θ = 90°. A direct comparison of the simulation techniques is shown in Figure 7. If we compare the peaks at θ = 90° at the same simulation time, it can be seen that it is higher for BD simulations than for simulations with HIs. In BD simulations, the trimers reorganize themselves more quickly in a triangular configuration than in simulations with HIs. It turns out that HIs slow down the reorganization of clusters. This conclusion is in agreement with the results obtained by Cao et al. using Stokesian dynamics simulations.11 It is also consistent with a recent observation concerning the kinetics of crystallization in a binary colloidal system, in which the transformation of a disordered aggregate into an ordered crystallite appeared

4. SUMMARY AND CONCLUSIONS The influence of HIs on the aggregation process has been studied in a system with one kind of colloids strongly attracting with each other. Three different simulation techniques have been employed: BD, which neglects the HIs; BD-YRP, which allows some of the HIs to be included but is limited to dilute

Figure 6. Angular distributions of the trimers at three different times and obtained from the three different simulation techniques: (a) BD, (b) BDYRP, and (c) SRD-MD. 14515



Article

Figure 7. Angular distributions of the trimers obtained from the three different simulation techniques (BD, BD-YRP, and SRD-MD) at the following times: (a) 0.03 s, (b) 0.06 s, and (c) 0.12 s.

suspensions because of the assumption of pairwise additivity of the interactions; and SRD-MD. In order to compare the aggregation kinetics in between these techniques, we have first checked that two relevant time scales for the aggregation phenomena are similar, i.e., the diffusion time and the dimer dissociation time. It is found that the inclusion of HIs does not affect dramatically the aggregation kinetics if we consider the whole process. For instance, the time at which a single cluster is obtained in the simulation is essentially the same in the different techniques. Nevertheless, the analyses performed exhibit some effects of the HIs in the details of the aggregation processes. It is worth noticing that these effects tend to increase with the increase of the colloid volume fraction. At the beginning of the aggregation, HIs slow down nucleation of clusters because they hinder dimer formation. This phenomenon is more pronounced within the SRD-MD technique. We attribute this difference between BD-YRP and SRD-MD to the lubrication forces that are not present in BD-YRP. In the second stage, when the aggregates grow by cluster−cluster coalescence, we find an opposite behavior, which causes a compensation effect on total aggregation times. In fact, HIs speed up the aggregation kinetics in that stage. This happens because the clusters diffuse faster in BD-YRP and SRD-MD compared to BD. The diffusion of a cluster in BD is artificially slowed down, with the diffusion coefficient decreasing as 1/N. Indeed, it does not decrease linearly as a function of its size, but it does as a function of the number of particles constituting the cluster, because of the addition of all the individual particle frictions. This can be problematic in the case of heteoag-

gregation between highly size-asymmetric particles, where a large particle can be surrounded even by hundreds of small oppositely charged particles. On the other hand, the BD-YRP technique exhibits larger cluster diffusion coefficients than SRD-MD. We think that the cluster diffusion coefficients obtained by BD-YRP are somewhat overestimated, since they even exceed the diffusion coefficient of the smallest sphere containing the cluster. This overestimate may be due to the approximations made in the YRP-BD technique, i.e., the pairwise additivity of interactions. The reorganization kinetics is also very important, because its competition with the kinetics of cluster−cluster coalescence determines the final shape of the aggregates.13 We have analyzed here the reorganization kinetics of trimers. It turns out that HIs slow down the formation of the energetically optimized triangular configurations, as it has already been observed by the Stokesian dynamics technique.11 Again, the phenomenon is more pronounced in BD-YRP than in SRDMD. It is worth noticing that the combined effects of a faster coalescence kinetics and a slower reorganization kinetics will lead to more elongated and more branched aggregates when HIs are taken into account. This is crucial for the determination of percolation thresholds and the formation of gel-like structures. To conclude, according to these results, BD-YRP and SRDMD simulations show similar effects of HIs on aggregation phenomena, even if their amplitude can differ. It seems that BD-YRP slightly overestimates these effects. Nevertheless, this technique can still be used for very dilute suspensions, where SRD-MD is computationally very demanding. 14516



■

Article

(11) Cao, X. J.; Cummins, H. Z.; Morris, J. F. Hydrodynamic and Interparticle Potential Effects on Aggregation of Colloidal Particles. J. Colloid Interface Sci. 2012, 368, 86−96. (12) Cerbelaud, M.; Videcoq, A.; Abélard, P.; Pagnoux, C.; Rossignol, F.; Ferrando, R. Self-assembly of Oppositely Charged Particles in Dilute Ceramic Suspensions: Predictive Role of Simulations. Soft Matter 2010, 6, 370−382. (13) Cerbelaud, M.; Ferrando, R.; Videcoq, A. Simulations of Heteroaggregation in a Suspension of Alumina and Silica Particles: Effect of Dilution. J. Chem. Phys. 2010, 132, 084701. (14) Piechowiak, M. A.; Videcoq, A.; Ferrando, R.; Bochicchio, D.; Pagnoux, C.; Rossignol, F. Aggregation Kinetics and Gel Formation in Modestly Concentrated Suspensions of Oppositely Charged Model Ceramic Colloids: A Numerical Study. Phys. Chem. Chem. Phys. 2012, 14, 1431−1439. (15) Cerbelaud, M.; Videcoq, A.; Abélard, P.; Ferrando, R. Simulation of the Heteroagglomeration between Highly Sizeasymmetric Ceramic Particles. J. Colloid Interface Sci. 2009, 332, 360−365. (16) Yamakawa, H. Transport of Polymer Chains in Dilute Solution: Hydrodynamic Interactions. J. Chem. Phys. 1970, 53, 436−443. (17) Jardat, M.; Bernard, O.; Turq, P.; Kneller, G. Transport Coefficients of Electrolyte Solutions from Smart Brownian Dynamics Simulations. J. Chem. Phys. 1999, 110, 7993−7999. (18) Tomilov, A.; Videcoq, A.; Chartier, T.; Ala-Nissila, T.; Vattulainen, I. Tracer Diffusion in Colloidal Suspensions under Dilute and Crowded Conditions with Hydrodynamic Interactions. J. Chem. Phys. 2012, 137, 014503. (19) Padding, J. T.; Louis, A. A. Hydrodynamic Interactions and Brownian Forces in Colloidal Suspensions: Coarse-graining over Time and Length Scales. Phys. Rev. E 2006, 74, 031402. (20) Ermak, D. L. A Computer Simulation of Charged Particles in Solution. I. Technique and Equilibrium Properties. J. Chem. Phys. 1975, 62, 4189−4196. (21) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, U.K., 1987. (22) Ando, T.; Skolnick, J. On the Importance of Hydrodynamic Interactions in Lipid Membrane Formation. Biophys. J. 2013, 104, 96− 105. (23) Bochicchio, D.; Videcoq, A.; Ferrando, R. Kinetically Driven Ordered Phase Formation in Binary Colloidal Crystals. Phys. Rev. E 2013, 87, 022304. (24) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38.

ASSOCIATED CONTENT

S Supporting Information *

A movie showing the fluid particles in between the two approaching colloids during dimer formation. Only the fluid particles placed in a sphere of the same size as the colloids and centered in between them are shown. A finite size is attributed to the fluid particles even if they are actually point-like particles. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We acknowledge funding from the French National Research Agency (ANR) under the Young Researchers Project No. ANR-09-JCJC-0086, HYSISCO. T.A.-N. has been in part supported by the Academy of Finland through its COMP CoE Grant No. 251748. D.B. acknowledges support from the Programme VINCI 2011 of the French-Italian University on the subject “Studio dell’eteroaggregazione di particelle colloidali ceramiche tramite tecniche di ottimizzazione globale”. A.V. thanks T.A.-N. and Aalto University School of Science for a Visiting Professorship. The movie in the Supporting Information has been obtained by VMD, a molecular graphics program originally designed for the interactive visualization and analyses of biological materials, developed by the Theoretical Biophysics Group in the Beckman Institute for Advanced Science and Technology at the University of Illinois at Urbana−Champaign.24

■

REFERENCES

(1) Gu, L.; Xu, S.; Sun, Z.; Wang, J. Brownian Dynamics Simulation of Crystallization Dynamics of Charged Colloidal Particles. J. Colloid Interface Sci. 2010, 350, 409−416. (2) Zhang, Z.; Glotzer, S. Self-assembly of Patchy Particles. Nano Lett. 2004, 4, 1407−1413. (3) Kim, A. Y.; Hauch, K. D.; Berg, J. C.; Martin, J. E.; Anderson, R. A. Linear Chains and Chain-like Fractals from Electrostatic Heteroaggregation. J. Colloid Interface Sci. 2003, 260, 149−159. (4) Cerbelaud, M.; Videcoq, A.; Abélard, P.; Pagnoux, C.; Rossignol, F.; Ferrando, R. Heteroaggregation between Al2O3 Submicrometer Particles and SiO2 Nanoparticles: Experiment and Simulation. Langmuir 2008, 24, 3001−3008. (5) Piechowiak, M. A.; Videcoq, A.; Rossignol, F.; Pagnoux, C.; Carrion, C.; Cerbelaud, M.; Ferrando, R. Oppositely Charged Model Ceramic Colloids: Numerical Predictions and Experimental Observations by Confocal Laser Scanning Microscopy. Langmuir 2010, 26, 12540−12547. (6) Ermark, D.; Cammon, J. M. Brownian Dynamics With Hydrodynamic Interactions. J. Chem. Phys. 1978, 69, 1352−1360. (7) Bossis, G.; Brady, J. Self-diffusion of Brownian Particles in Concentrated Suspensions under Shear. Phys. Rev. E 1987, 53, 5044− 5050. (8) Phung, T.; Brady, J.; Bossis, G. Stokesian Dynamics Simulation of Brownian Suspensions. J. Fluid Mech.. 1996, 313, 181−207. (9) Furukawa, A.; Tanaka, H. Key Role of Hydrodynamic Interactions in Colloidal Gelation. Phys. Rev. Lett. 2010, 104, 245702. (10) Whitmer, J. K.; Luijten, E. Influence of Hydrodynamics on Cluster Formation in Colloid-Polymer Mixtures. J. Phys. Chem. B 2011, 115, 7294−7300. 14517


Aggregation in Colloidal Suspensions: Evaluation of the Role of

Recommend Documents