A New Stochastic Approach for the Simulation of ... - ACS Publications

Article pubs.acs.org/Langmuir

A New Stochastic Approach for the Simulation of Agglomeration between Colloidal Particles Christophe Henry,*,† Jean-Pierre Minier,‡ Jacek Pozorski,† and Grégory Lefèvre§ †

Institute of Fluid-Flow Machinery, Polish Academy of Sciences ul. Fiszera 14, 80-231 Gdańsk, Poland EDF R&D, Fluid Dynamics, Power Generation and Environment 6 quai Watier, 78400 Chatou, France § LECIME, CNRS-Chimie ParisTech, UMR 7575 11 rue Pierre et Marie Curie, 75005 Paris, France ‡

ABSTRACT: This paper presents a stochastic approach for the simulation of particle agglomeration, which is addressed as a two-step process: first, particles are transported by the flow toward each other (collision step) and, second, short-ranged particle−particle interactions lead either to the formation of an agglomerate or prevent it (adhesion step). Particle collisions are treated in the framework of Lagrangian approaches where the motions of a large number of particles are explicitly tracked. The key idea to detect collisions is to account for the whole continuous relative trajectory of particle pairs within each time step and not only the initial and final relative distances between two possible colliding partners at the beginning and at the end of the time steps. The present paper is thus the continuation of a previous work (Mohaupt M., Minier, J.-P., Tanière, A. A new approach for the detection of particle interactions for large-inertia and colloidal particles in a turbulent flow, Int. J. Multiphase Flow, 2011, 37, 746−755) and is devoted to an extension of the approach to the treatment of particle agglomeration. For that purpose, the attachment step is modeled using the DLVO theory (Derjaguin and Landau, Verwey and Overbeek) which describes particle−particle interactions as the sum of van der Waals and electrostatic forces. The attachment step is coupled with the collision step using a common energy balance approach, where particles are assumed to agglomerate only if their relative kinetic energy is high enough to overcome the maximum repulsive interaction energy between particles. Numerical results obtained with this model are shown to compare well with available experimental data on agglomeration. These promising results assert the applicability of the present modeling approach over a whole range of particle sizes (even nanoscopic) and solution conditions (both attractive and repulsive cases).

1. INTRODUCTION 1.1. Agglomeration Process. A large variety of industrial and academic systems involve suspensions where the formation of larger bodies through collision and adhesion of two smaller bodies can occur. With dependence on the context, this phenomenon is referred to as agglomeration, coagulation, flocculation, or aggregation. As underlined in recent reviews,1,2 agglomeration phenomena involve bodies of a different nature (i.e., solid particles, colloidal particles, biological organisms) and have applications in various fields (among which filtration technologies, water treatment facilities, or papermaking industry). In a process opposite to agglomeration, larger bodies can also be split in smaller bodies: this phenomenon is called fragmentation.3 Fragmentation generally occurs together with agglomeration (the rate of formation of larger bodies is then given by the balance between agglomeration and fragmentation). The present study aims at understanding, modeling, and simulating the phenomena related to particle agglomeration and fragmentation with particular attention to colloidal particles. Since fragmentation cannot occur without previous agglomeration, the present study is focused more specifically on the agglomeration of colloidal particles. As for the formation of multilayers of particles on a surface,4 the agglomeration process for colloidal particles results from © 2013 American Chemical Society

the coupling between two main interactions: (1) Particle−fluid interactions, which play a role in the motion of particles within a flow and govern the number of particle−particle encounters, as well as their time and spatial locations. (Gravity is irrelevant for colloidal particles but should be added when large inertial particles are considered.) This first step toward agglomeration is called the collision step. (2) Particle−particle interactions, which control whether colliding particles will adhere (in the case of attractive or weakly repulsive interactions) or simply bounce (for highly repulsive interactions). This second step toward agglomeration is thus called the adhesion (or attachment) step. 1.2. Modeling Approaches for Particle Agglomeration. From a modeling point of view, the first step (i.e., the collision step) gathers issues related to particle transport by a flow coupled with a proper description for the detection of interparticle collisions. For the adhesion step, particle−particle interactions are often described using the DLVO theory (named after the work of Derjaguin and Landau5 and Verwey and Overbeek6), which defines interparticle forces as the sum of van der Waals and double-layer electrostatic contributions.7−10 Received: June 27, 2013 Published: October 10, 2013 13694

dx.doi.org/10.1021/la403615w | Langmuir 2013, 29, 13694−13707

Langmuir

Article

or elastic collision) is then detailed in Treatment of the Collision. The validation of the method is then presented in Numerical Results for Particle Agglomeration. For that purpose, numerical results for particle agglomeration under diffusive motions are compared to various theoretical and experimental data. First, numerical predictions of the agglomeration of initially monodispersed particles undergoing Brownian motions in the fast aggregation regime (i.e., with attractive particle−particle interactions) are compared to the theoretical expectations obtained by Smoluchowski11 in Numerical Results for Rapid Agglomeration. Second, numerical results are compared to three different experimental studies, each of which is representative of a particular issue encountered in colloidal agglomeration. This allows one to assess the ability of the present stochastic method to predict complex situations such as the impact of the morphology of agglomerates (see Effect of Agglomerate Geometry on the Stability Ratio), surface roughness (see Impact of Surface Roughness on the Stability Ratio), and the effect of solution conditions (in Evolution of Particle Agglomeration with Particle Concentration) on macroscopic properties (average particle radius or stability ratio).

The range of these physicochemical interfacial forces (DLVO theory) is limited to a few tens of nanometers, a distance much smaller than typical interparticle distances in fluids. Due to the different nature and scales involved in the collision and adhesion steps, they are often treated separately when establishing a modeling approach.7 Two main modeling approaches have been proposed in the literature to address particle agglomeration: (1) Approach based on population balance equation (PBE), which consists of solving population balance equations for a set of characteristic particle classes chosen from the outset.7,11 Particle−fluid interactions are included in the collision kernel β, which describes the frequency of particle collisions in given hydrodynamic conditions, while particle−particle interactions influence the stability ratio (W), which is the ratio of the total number of collisions occurring divided by the number of collisions leading to an agglomerate. These kernels must be assumed beforehand and are therefore inputs of the method. (2) Particle-tracking approach, where the motions of a large number of particles are explicitly tracked and their collisions are detected and recorded. The difficulty that arises in such approaches is to choose a proper method to detect the collision between particles. Various methods have been proposed in the literature, among which exact deterministic methods12−14 valid in the ballistic regime and, more recently, stochastic methods15,16 valid in the diffusive regime. In contrast with PBE approaches, the collision and agglomeration kernels are now calculated from the simulations and are outputs of the approach. 1.3. Purpose of the Paper. The aim of the present study is to develop and assess a method for the simulation of agglomeration between colloidal particles that remains valid and tractable in complex situations (such as nonhomogeneous flows, turbulent flows, etc.). Due to the lack of analytical formulas for the collision kernel β in such complex situations, the present method is developed in the framework of particletracking approaches. The developed method thus proposes to determine the collision kernel in such cases and numerical results can thus be used as an input in classical PBE approaches. As mentioned above, classical deterministic algorithms have been devised. However, these deterministic algorithms are valid only in the ballistic regime - i.e., for time steps (Δt) much smaller than the particle relaxation time τp (Δt ≪ τp) - and are not adapted to simulate the agglomeration of nanoparticles. Therefore, the present method is in continuation of a previous work on a stochastic method for the detection and treatment of particle collisions (more details in ref 15 and in a forthcoming paper16) and proposes an extension of this method to simulate the whole agglomeration phenomenon. More specifically, since the stochastic method for the collision step has been validated in the case of particles undergoing diffusive motions, the accuracy of the present method will be assessed by comparing numerical results to experimental data obtained for particles undergoing diffusive motions. 1.4. Organisation of the Paper. Within this scope, the present paper is organized as follows. The modeling approach is first presented in Modeling Approach for Particle Agglomeration. For the sake of clarity, the stochastic model for the detection of collisions is briefly introduced in Stochastic Approach for the Detection of Particle Collisions followed by a description of the method used for the attachment step and its coupling with the collision step in Agglomeration Criterion. The numerical treatment of the collision event (agglomeration

2. MODELING APPROACH FOR PARTICLE AGGLOMERATION In this section, the stochastic method for the simulation of particle agglomeration is described. The equations of particle motions as well as the detection and treatment of particle collisions have been described in detail in a previous paper15 (and also in a forthcoming paper16), and therefore only a brief overview of the collision step is given here. Special attention is thus put on describing the attachment step as well as its coupling with the collision step. 2.1. Stochastic Approach for the Detection of Particle Collisions. As described in other papers,15,16 the stochastic method developed for the simulation of particle collisions is based on a Lagrangian approach where the motions of N particles are explicitly tracked. The first step in the simulation of particle collisions is to properly detect collisions for pairs of particles during the calculation (see Figure 1) and to evaluate the collision location in time and space (see Figure 2). For that purpose, the present work follows a stochastic method valid at large time steps.15,16 Particle Motions. In the framework of particle-tracking methods, the equations of particle motions in the diffusive regime [i.e., for time steps Δt much larger than the particle relaxation time τp (Δt ≫ τp)] are given by

Figure 1. Sketch of the stochastic approach used to detect particle collisions depending on the (a) initial and (b) final distances between a pair of particles. 13695


Langmuir

Article

B2x,i + B2x,j.15,16 The collision probability for each pair of particles is then given by the probability Pba(Rij) that the distance reaches a value lower or equal to the sum of the two particle radii (Rij = Ri + Rj) during one time step (see Figure 1). This probability is calculated using a so-called diffusive bridge and is given by15,19 ⎛ 2R ij(a + b − R ij) ⎞ ⎟ − 1 exp⎜ ⎝ Bx2,ijΔt ⎠ b Pa (R ij) = ⎛ 2ab ⎞ exp⎜ 2 ⎟ − 1 ⎝ Bx ,ijΔt ⎠

where a and b represent the initial and final distances between this pair of particles, that is their distance at the beginning and at the end of the time step. This approach has been validated previously for particles undergoing diffusive motions using large time steps (diffusive regime).15 Collision Location. Knowing that a collision occurred between a pair of particles, the collision location in space and time (see Figure 2) is determined using appropriate stochastic models valid in nonhomogeneous situations (more details in a forthcoming paper16). The first step toward a proper description of the collision event is to determine when it occurred. For that purpose, the time of collision (tc) is evaluated using a probabilistic approach based on the notion of the first passage time: similar to the detection of particle collision, the probability density function of the collision time, p(tc), is determined respecting the initial and final distances between particles and their relative motions during the time step.16 Then, the instant of collision being known, the position of both particles upon collision can be determined. Since the motions of the barycenter of a pair of colliding particles are not affected by the collision event (due to the conservation of momentum), the spatial location of the collision is evaluated by considering its trajectory Xg(t). The barycenter follows also a diffusive motion conditioned by its initial and final positions that are already known and the location of the collision is taken as the location of the barycenter at time tc.16 2.2. Agglomeration Criterion. The question that arises once a collision between a pair of particles has been detected is to determine whether the collision is elastic (bouncing) or leads to the formation of a larger agglomerate (adhesion). For that purpose, particle−particle interactions need to be described properly in the model. DLVO interactions. The DLVO theory considers that interbody interactions are the sum of the van der Waals (VDW) interaction and of the electrostatic double-layer (EDL) interaction.4,7,9,10,20 In the present work, VDW and EDL interactions are calculated using the following formulas: (1) retarded VDW potentials UVDW (or force FVDW) between two spheres have been estimated using Gregory’s formula21 valid for small separation distances, h ≪ Ri:

Figure 2. Sketch of the stochastic approach used to determine the collision time and its spatial location [based on the motions of the barycenter Xg(t)].

dX p(t ) = Us dt + Bx d W(t )

(1)

where Xp is the particle location, Us is the fluid velocity seen by the particle {Us(t) = Uf[t,xp(t)] where Uf(t,x) is the instantaneous fluid velocity field}, dW is an increment of the Wiener process and Bx is the diffusion coefficient [Bx = (2kBTτp/mp)1/2] for Brownian motions, with kB the Boltzmann constant, T the temperature, and mp the particle mass). It should be noted that particle velocities Up still exist but that they are fluctuating rapidly and, in a discrete sense, they are represented by a series of independent Gaussian variables. A first-order numerical scheme used to calculate the motion of each particle is given by (see also in refs 16−18): ⎧ X (t + Δt )≃ X (t ) + U (t )Δt + B Δt ξ p p s x x ⎪ ⎪ ⎨ 1 ξu ⎪ Up(t + Δt )≃ Us(t ) + Bx ⎪ 2 τp ⎩

(3)

(2)

(ρpd2p/18ρfνf)

with τp = the particle relaxation time related to drag forces only (dp being the particle diameter, ρp the particle density, ρf the fluid density, and νf the fluid kinematic viscosity), ξx and ξu representing independent random variables sampled in the standard normal distribution, and Δt is the time step of the simulation. In the following, only particles undergoing diffusive motions are considered (more particularly Brownian motions) and, in that case, the fluid velocity is taken as zero (fluid at rest). Perspectives on an extended model to simulate particle agglomeration in complex flows (with a velocity Us given by a fluid calculation) are discussed in the Conclusion. Detection of Particle Collisions. Once particle motions during a time step have been calculated, the collision for each pair of particles is evaluated. The difficulty that arises when evaluating a collision between a pair of particles undergoing diffusive motions during large time steps is that particle trajectories do not follow straight lines. Therefore, classical approaches based on the knowledge of the initial (meaning at the beginning of the time step) and final (meaning at the end of the time step) positions of particles are not suitable here.12−14 A new stochastic approach has thus been previously developed in order to account for the whole continuous trajectories of particles within one time step.15 The key idea is to study the relative motions between each pair of particles, which are also described by diffusive motions with a diffusion coefficient B2x,ij = 13696


Langmuir

Article

⎧ R iR j A ⎪UVDW = − Ham 6h R i + R j ⎪ ⎪ ⎡ λ ⎞⎤ 5.32h ⎛ ⎟⎥ × ⎢1 − ln⎜1 + ⎪ ⎝ ⎣ λ 5.32h ⎠⎦ ⎪ ⎨ AHam R iR j ⎪ ⎪ FVDW = 6h2 R i + R j ⎪ −1 ⎡ ⎪ ⎛ λ ⎞ ⎤ ⎜ ⎟ ⎢ ⎥ × − + 1 1 ⎪ ⎝ 5.32h ⎠ ⎦ ⎣ ⎩

(4)

with AHam the Hamaker constant between materials and λ the characteristic retardation wavelength (that measures the distance above which retardation effects play a role). (2) EDL interaction potential UEDL (or force FEDL) between two spheres is described using the formula from Bell et al.7,10 (valid at all separation distances provided that h ≪ κ−1 and κR > 5, where κ−1 = [ε0εrkBT/(2e2I)]1/2 is the Debye length and I the ionic strength):

Figure 3. DLVO interaction energy between two particles of 1 μm in the case of repulsive interactions: VDW interaction energy (red dashed line), EDL interaction energy (blue dotted and dashed line), and DLVO interaction energy (black solid line).

previous studies on particle deposition and re-entrainment22,23): agglomeration can occur only if the relative kinetic energy of the pair of colliding particles Ekin,ij is high enough to overcome the energy barrier Ebarr,ij:

R iR j(r − R i)(r − R j) ⎛ k T ⎞2 UEDL = 2πε0εR ⎜ B ⎟ × ⎝ ze ⎠ r[(R i + R j)r − R i2 − R j2] ×[Ωi ln(1 + Γ) + Ωj ln(1 − Γ)]

⎧ if E kin, ij ≥ E barr, ij ⇒agglomeration ⎪ ⎨ ⎪ ⎩ if E kin, ij < E barr, ij ⇒elastic collision

(5)

with Ωi =

Φi2

+

Φ2j

Therefore, the coupling between the collision and adhesion steps is rather straightforward and very simple to implement numerically, while still accounting exactly for the key notion of the energy barrier. It is worth adding that, due to the Lagrangian point of view, the variations of the energy barrier (for example, its random nature due to particle roughness) is handled without any approximation by the present approach. This point will be illustrated in the validation results. 2.3. Treatment of the Collision. The collision between a pair of particles being detected and the agglomeration criterion being calculated, the collision event (either agglomeration or elastic collision) has to be treated properly (see also Figure 4). Treatment of Particle Agglomeration. When the collision gives rise to agglomeration, one of the partners is removed from the simulation while the particle radius of the other partner is increased (see also Figure 4). At this stage, two issues arise: (1) the first one can be stated as follows. How do we properly calculate the DLVO interaction between two agglomerates? In the present paper, agglomerate−agglomerate interactions have been calculated, considering the interaction between two elementary (initial) particles that form the agglomerate. This assumption will be addressed more specifically later when comparing numerical results to experimental data from Snoswell et al.24 in Effect of Agglomerate Geometry on the Stability Ratio. (2) The second issue concerns the size of the agglomerate formed by the collision. In the present paper, the volume of the agglomerate is assumed to be equal to the sum of the two particle volumes and the agglomerate radius is thus given by Ragg = (R3i + R3j )1/3. This amounts to considering a fractal dimension of the agglomerate formed equal to 3 (compact aggregates).2,25 Consequently, polydispersity in particle sizes naturally results from agglomeration in the present modeling approach and the motions of agglomerates are then simulated respecting their corresponding size. Yet, the choice of such an aggregation law for the prediction of the

+ ΛΦΦ i j

Ωj = Φi2 + Φ2j − ΛΦΦ i j Λ=

Γ=

R i(r − R i) + R j(r − R j) R iR j (r − R i)(r − R j)

(6)

R j(r − R j) R i(r − R i) × e κ(R i + R j − r)

where Φi, and Φj are the reduced potentials Φ = (zeφ/kBT) [ϕ is the electric potential, commonly assumed to be equal to the zeta potential (ζ)], r is the center-to-center distance, ε0 is the dielectric permittivity of vacuum (ε0 = 8.854 × 10−12 C V−1 m−1), εR the dielectric constant of water (εR = 78.5), kB the Boltzmann constant (kB = 1.38 × 10−23 J K−1), T the temperature (in K), e the electron charge (e = 1.6 × 10−19 C), and z the valence of ions present in the solution. A typical DLVO interaction energy is displayed in Figure 3 for two 1 μm particles interacting (zeta potential equal to 25 mV in a 10 mM electrolyte, Hamaker constant of 1.0 × 10−20 J with a retardation wavelength of 100 nm). It can be seen that the VDW attraction is predominant at small and large separation distances, while the EDL repulsion dominates the interaction energy at intermediate separations. This peak in the DLVO curve corresponds to the so-called energy barrier, Ebarr, which can prevent incoming particles from reaching a contact point. Coupling between Attachment and Collision Step. The present particle-tracking approach has been developed in the case of large time steps, where particles travel distances much larger than the range of DLVO forces during one time step. Therefore, the natural approach is to consider DLVO interactions as collision boundary conditions (similarly to 13697


Langmuir

Article

Figure 4. Sketch of the stochastic approach used for the treatment of particle collisions.

resulting diameter may be debated and corresponds to one possibility among several possible evolutions. For example, if we consider an agglomerate formed by a large number of primary particles, we can expect this aggregate to also contain some void, or some internal porosity, whereas the present estimate will not only assume a global spherical shape but also will neglect this void content. This question will be reconsidered in Evolution of Particle Agglomeration with Particle Concentration. Treatment of Elastic Collisions. In the case where the energy barrier is higher than the kinetic energy of particles, an elastic collision is assumed, and the positions of colliding particles at the end of the time step are updated to account for particle motions after collision (see also Figure 4a). For that purpose, the motion of one particle from its position at the collision Xi(t + tc) is estimated as a new random−walk displacement over the rest of the time step (Δt−tc), using eq 1. The position of the second particle at the end of the time step is then obtained, respecting the position of the barycenter of the pair of particles at the end of the time step Xg(t + Δt) (which is unchanged due to the momentum conservation during the collision). The stochastic approach for collision has only been validated by comparing the number of collisions obtained numerically to theoretical expectations which have been derived in simplified situations for both initially monodispersed and polydispersed suspensions (more details will be available in ref 16). It has been shown that it provides accurate results of particle collisions in the diffusive regime provided that the root-meansquare (r.m.s.) displacement of particles during one time step, σx = Bx(Δt)1/2, is larger than the particle diameter, dp (σx ≫ dp).

presented of the effects of surface roughness as well as that of agglomerate morphology. 3.1. Numerical Results for Rapid Agglomeration. To assess the validity and accuracy of the present method, numerical calculations have been first performed with particles undergoing Brownian motion in the case of rapid agglomeration (i.e., for attractive particle−particle interactions). In that case, the limiting phenomenon for particle agglomeration is the collision step, since each collision ends up in the formation of a larger agglomerate. Approximate Analytical Expressions. In the framework of population balance approaches, approximate analytical expressions have been derived to describe the evolution of particle concentrations with time for rapid agglomeration.7,11 In the case of particles undergoing Brownian motions, the collision kernel between classes i and j is given by βij = 2kBT(Ri + Rj)2/(3RiRjμ). Next, considering that all collision kernels are equal (assumption valid at the initial stages of agglomeration), the evolution of the total number density of particles nT = Σi ∞ = 1ni starting from n0 initially monodispersed particles is given by7,11 n0 nT = 1 + βn 0 t (7)

3. NUMERICAL RESULTS FOR PARTICLE AGGLOMERATION In the following section, the new approach for particle agglomeration presented in this work is validated. For that purpose, numerical results obtained for particle agglomeration are compared to theoretical expectations on rapid agglomeration (see Numerical Results for Rapid Agglomeration) as well as to experimental data on the agglomeration of colloidal particles (see Comparison to Experimental Data on Particle Agglomeration). Since numerical calculations are performed for particles undergoing Brownian motions, accurate results are expected for smooth particles and further analyses are

Numerical Calculation with the Present Approach. Two numerical calculations using the stochastic approach for particle agglomeration have been performed by tracking the motions of 480 (respectively 1000) particles of 1 μm undergoing Brownian motions at room temperature (T = 296.15K) in a periodic box of size L = 500 μm. The corresponding initial volume density is n0 = 2.01 × 10−6 m−3 (n0 = 4.19 × 10−6 m−3, respectively). The time step used for the calculation of particle motions is Δt = 106τp = 3.03 s (i.e., within the diffusive regime) and the corresponding r.m.s. displacement of particles within one time step is larger than the particle diameter dp (here equal to 12dp). Upon collision, the two particles have been assumed to be irreversibly bound together (no fragmentation), and the size of

Equation 7 can be rewritten in terms of the characteristic aggregation time (or half-life of aggregation) τa = 3 μ/(4kBTn0). Similarly, formulas describe the evolution of the number of singlets (n1), doublets (n2), and triplets (n3) (see also ref 7). The number density of particles of class i is then given by ni =

13698

n0(t /τa)i − 1 (1 + t /τa)i + 1

(8)


Langmuir

Article

Figure 5. Theoretical (lines) and numerical (points) values for the concentration of primary particles n1 (red △), doublets n2 (green ▽), triplets n3 (blue □), and total particles nT (black ○) as a function of the reduced time t/τa in the case of rapid agglomeration.

Figure 6. Measurements of titania particle properties from TEM images and dynamic light scattering. Reprinted with permission from ref 24. Copyright 2005 Elsevier.

3.2. Comparison to Experimental Data on Particle Agglomeration. To further validate the present approach, numerical calculations have also been performed in conditions similar to those encountered experimentally and numerical results have been compared to experimental data. In the following, comparisons between experimental data and numerical results are presented for three cases, each of which is representative of a particular issue encountered in the agglomeration of colloidal particles. This validation procedure is therefore interesting in order to test the prediction capacities of the new approach in different situations. First we assess whether the stochastic approach can account correctly for the effect of the morphology of agglomerates. Second, we investigate if the approach captures the influence of surface roughness on the agglomeration rate. Third, we determine if the effect of the solution conditions (especially the particle concentration, pH, and ionic strength) on macroscopic properties (here the average particle radius) is properly reproduced. 3.2.1. Effect of Agglomerate Geometry on the Stability Ratio. Experimental Data. Snoswell et al.24 have experimentally studied the effect of various ionic strengths (ranging from 0.1 to 100 mM) and pH (6.3, 6.7, and 8.4) on the agglomeration of synthetic titania particles immersed in a

the larger agglomerate created is considered to be equal to Ragg (see Treatment of the Collision). In Figure 5, numerical results obtained with the stochastic approach (symbols) are compared to the theoretical values (straight lines). It can be seen that the evolution of the concentration of singlets, doublets, triplets, and total particles as a function of the reduced time t/τa is well-reproduced numerically. Indeed, at the initial stages of agglomeration, the number of singlets is decreasing (fragmentation is not considered), while the number of doublets (and triplets also) is increasing due to the agglomeration between primary particles. After some time, the number of doublets (and later triplets) starts to decrease due to the agglomeration of such agglomerates with other particles to form larger ones. It should be noted also that the accuracy is improved when the number of initial particles is increased (a better statistical convergence is achieved). Summary. This first validation confirms the ability of the present approach to capture the agglomeration rate of particles using large time steps that allow tractable calculations even in complex situations. Therefore, compared to classical particletracking approaches that are restricted to small time steps (i.e., Δt ≪ τp), the present approach yields similar results for much lower computational costs. 13699


Langmuir

Article

Numerical Results. Numerical results are shown in Figure 7, which displays the evolution of the stability ratio W with the

KCl electrolyte. For that purpose, the authors have measured the stability ratio (W), which is the ratio between the agglomeration rate in attractive conditions (i.e., the fast aggregation rate) and the measured aggregation rate in repulsive conditions. By definition, the stability ratio W is thus equal to 1 for pH and solution conditions where particle− particle interactions are attractive and is greater than one for repulsive particle−particle interactions. In other words, the higher the repulsion between particles, the higher the stability ratio W. This stability ratio is therefore often used to extract the solution conditions (pH and ionic strength) where particles are stable (i.e., where the stability ratio is high) (no agglomeration occurs). The size of titania particles was measured using multiangle light-scattering techniques24 (see Figure 6b), showing particles with a diameter of 300 nm approximately (90% of particles have a diameter ranging between 270 and 330 nm). However, as displayed in Figure 6a, typical transmission electron microscope (TEM) images of titania particles revealed that they are composed of smaller crystallites with a radius on the order of 10−15 nm. In addition, the zeta potentials of titania particles were derived experimentally from electromobility measurements. An empirical formula for the evolution of particle zeta potentials with ionic strength was then extracted, assuming small zeta potentials in a monovalent electrolyte. In that case, a linear relation between the logarithm of the potential and the ionic strength can be written24 as log(ψ) ∝ κ. Numerical Parameters. Numerical calculations have been performed using the stochastic approach for the simulation of particle agglomeration for Brownian particles immersed in a fluid at rest. The motions of titania particles with a mean diameter of 300 nm and a standard deviation of 10 nm have been tracked. The time step (Δt = 2.1 × 105 τp ≃ 0.057s) has been chosen so that the r.m.s. displacement of particles respects the criterion σx ≫ dp (here σx ≃ 10dp), and the size of the periodic box has been chosen equal to Lbox = 50σx (yielding Lbox = 150 μm). The particle mass concentration being equal to 0.023 kg m−3,24 the number of particles (Npart) in the simulation box is equal to 1220. Since only the stability ratio is evaluated numerically, calculations are performed using a particular treatment: once two colliding particles have a relative kinetic energy higher than the energy barrier (and thus will agglomerate), the numerical collision efficiency is incremented, but we do not merge this pair of particles into one bigger agglomerate and keep them as separate particles for the rest of the calculation. This treatment allows us to keep a constant number of particles within the domain and, therefore, to have a significant number of collisions (here more than 10000 collisions) within a reasonable time. Since the stability ratio is related to the ratio between the number of collisions and the number of collisions leading to agglomeration, this simplified treatment does not affect the predicted stability ratio and allows fast numerical predictions of W with a good-enough statistical precision. Similar results can also be obtained with a larger number of initial particles and a smaller elapsed time. The Hamaker constant used in the calculation of the van der Waals interactions has been chosen equal to 1.0 × 10−20J24 (this value has been thoroughly discussed by the authors since particles have a complex geometry). The zeta potentials of particles have been calculated using the empirical formulas obtained by Snoswell et al.24

Figure 7. Experimental data (points) and numerical values (lines) for the stability ratio (W) of 300 nm titania particles at various ionic strengths (I) for different pH values: pH 6.3 (black Δ), pH 6.7 (red □), and pH 8.4 (green ○).

pH value: it can be seen that pH values at which the stability ratio starts to rise are correctly reproduced by the numerical simulations. However, the slope of log W/log I is not recovered numerically. As discussed by the authors,24 the slope is highly correlated with the specific geometry of the particles. In the present case, we are dealing with particles with a diameter of 300 nm but which are actually composed of smaller 10−15 nm crystallites. Therefore, if the transport of these particles must indeed be evaluated with the diameter of the large composite particles (the 300 nm particles), the interaction between two such complex particles should be described by the interaction between smaller crystallites (having a diameter of around 10− 15 nm) instead of using the whole particle radius. To further confirm this trend, a second calculation has been performed considering that, upon collision, particle−particle interactions are now given by the interaction between two smaller 15 nm titania crystallites. The results are plotted in Figure 8, showing much better agreement with measurements. More generally, these good results confirm the fact that, when evaluating the adhesion between two agglomerates, the equivalent radius of the two aggregates (or particles) cannot be

Figure 8. Experimental data (points) and numerical values (lines) for the stability ratio (W) of titania particles (composed of smaller 15 nm nanoparticles) at various ionic strengths (I) for different pH values: pH 6.3 (black Δ), pH 6.7 (red □), and pH 8.4 (green ○). 13700


Langmuir

Article

Figure 9. Measured (points) and modeled (lines) zeta potentials of latex particles at different ionic strength for various pH values from Behrens et al.:29 I = 1 (black Δ), 10 (red ●), 100 (green ◇), extrapolated at I = 300 mM (blue line). Adapted from ref 29. Copyright 2000 American Chemical Society.

for 155 nm particles and Δt = 1.8 × 104τp ≃ 0.59 ms for 52 nm particles. The volume fraction being equal to 0.04,29 the number of particles in the simulation box is equal to 1081. Since only the stability ratio is evaluated numerically, colliding particles have been assumed not to form a larger agglomerate and remained thus separated (even when the criterion Ekin > Ebarr is respected). Particle−particle interactions have been evaluated using a value of the Hamaker constant of 0.9 × 10−20 J for all particle sizes10,29 and with a characteristic retardation wavelength of 1 μm (retardation effects thus play a minor role). The particle zeta potentials were measured experimentally at ionic strengths equal to 1, 10, and 100 mM by Behrens et al.,29 (see Figure 9, panels a and b). Since the zeta potentials of 155 nm particles at an ionic strength of 2 mM were not available, they have been extrapolated from available experimental data. For that purpose, we have followed the approximation retained by Snoswell et al.24 who considered a linear relation between the logarithm of zeta potentials and the ionic strength (this approximation is valid provided that the potentials are small enough |ψ| ≤ 50 mV):

used to perform a calculation since it does not account for the microscopic properties linked with agglomerate morphology. Similar trends have been underlined recently by Babick and Schießl,26,27 who calculated the interaction between two fractal aggregates. They have shown that the approximation of primary particle interaction (which considers only the interaction between the closest pair of primary particles) is valid at small separation distances and provided that κR ≫ 1. Within this paper, as a first approximation, we have thus chosen to calculate the interaction between two elementary particles (which make up larger particles or agglomerates). Yet, in reality, the intricate morphology of agglomerates may lead to multiple contacts (and thus interactions) between smaller particles composing each agglomerate. Summary. This first validation with experimental data shows that the present approach provides accurate results for the effect of particle agglomeration on macroscopic parameters (here the stability ratio which is related to the initial agglomeration rate). The use of a particle-tracking approach coupled with a proper description of particle−particle interactions also allows one to capture the effect of local nanoscopic properties linked with the morphology of the agglomerate. Similar conclusions have been reached recently using a formulation for the stability ratio (W) that is typically used in PBE approaches.28 3.2.2. Impact of Surface Roughness on the Stability Ratio. Experimental Data. Numerical results are compared to the experimental data from Behrens et al.,29 who studied the agglomeration of carboxyl latex particles in water. The authors measured the agglomeration of two sizes of latex particles (155 and 52 nm in radius) immersed in an electrolyte solution with various ionic strengths (1, 2, 10, and 100 mM) over a wide range of pH values (from 3 to 10). The agglomeration of Brownian latex particles was recorded using light scattering techniques and displayed in terms of the stability ratio (W). Numerical Parameters. Numerical calculations of the agglomeration of latex particles have been performed using the present approach for the simulation of particle agglomeration, considering Brownian particles immersed in a fluid at rest at room temperature. The time step has been chosen so that the r.m.s. displacement respects σx ≃ 5dp and the size of the periodic box has been chosen equal to Lbox = 10σx (Lbox = 15 μm for 155 nm particles and Lbox = 5.5 μm for 55 nm particles, respectively). This corresponds to Δt = 54 × 104τp ≃ 0.016 s

ψ = ψ0e−κh

(9)

with ψ0 and h being constant parameters determined using experimental data available at other ionic strengths. Similarly, to obtain the zeta potentials of 155 nm particles at an ionic strength of 300 mM, the same formula has been used. However, eq 9 is valid only for small zeta potentials, whereas, here, zeta potentials at I = 10 mM reach values lower than −100 mV for pH values larger than 7. Therefore, the obtained values at I = 300 mM were used with great care. As depicted in Figure 9a, the zeta potentials of the 155 nm particles obtained at I = 300 mM show a similar trend to those observed for 10 and 100 mM: the zeta potential is decreasing with increasing pH and a plateau value (here around −30 mV) is reached for large pH values (here above 6.7). Numerical Results. The numerical results obtained are plotted in Figures 10 and 11. It can be seen that for low pH values, particle zeta potentials are too small to give rise to significant repulsive particle−particle interactions, and the stability ratio is thus equal to 1 both experimentally and numerically. Yet, upon increasing the pH value, particle zeta potentials are higher, leading to more and more repulsive interactions and thus higher values of the stability ratio W. 13701


Langmuir

Article

agglomerates, which are composed of smaller particles. The presence of surface roughness on particles has a similar impact on the stability ratio as particle morphology since the interaction between the closest pair of asperities will govern the overall interaction. The interaction energy being proportional to the body radius, smaller asperities on a surface will reduce the repulsive interaction energy between particles (this is also true for particles interacting with rough surfaces in deposition and reentrainment phenomena22,23). In the present case, given that surface roughness was not characterized experimentally, the impact of surface roughness is analyzed using roughness parameters that appear plausible with other experimental data (for instance, Popa et al.30 measured a roughness up to 1.3 nm on latex particles while a surface roughness as high as 10 nm could be seen on hematite nanoparticles31). Asperities with sizes ranging from a few nanometers up to a few tens of nanometers and covering a small amount of the particle surface are thus introduced in the simulation. The interaction between two rough particles is then evaluated using a modeling approach which provides statistical information on the impact of surface roughness on the interaction energy. As in the deposition and re-entrainment of colloidal particles,22,23 particle roughness has been modeled considering perfectly spherical particles covered by hemispherical asperities. The number of asperities between the two colliding particles Nasp is given by a Poisson distribution with a mean value respecting the surface covered by asperities Scov and the asperity size Rasp (more details can be found in ref 23):

Figure 10. Experimental data (points) and numerical values (lines) for the stability ratio (W) of smooth 155 nm particles at various pH values for different ionic strengths: 2 mM I = 100 (black Δ), 10 (red ●), 100 (green □), and 300 mM (blue ▽).

Nasp =

R (2R part + R asp) ⎞ Scov 2R partR asp ⎛ ⎜⎜2R part − asp ⎟⎟ 4R part R asp 4R part ⎝ ⎠ (10)

Figure 11. Experimental data (points) and numerical values (lines) for the stability ratio (W) of smooth 52 nm particles at various pH values for different ionic strengths: 1 (black Δ), 10 (red ●), and 100 (green ◇).

The energy barrier between two rough particles is calculated assuming interaction energy to be additive and is given by the weighted sum of particle−particle and particle−asperities interactions. A first calculation has been performed introducing 5 nm asperities on the surface of particles with a surface coverage of 0.015%. The effects of such surface roughness on the stability ratio are shown in Figure 12: it can be seen that the numerical predicitons are in better agreement with experimental data,

Numerical results yield relatively correct predictions of the pH value at which the stability ratio starts to increase. However, the slope of the stability ratio with pH values in the case of highly repulsive particle−particle interactions is not always satisfactory (especially at an ionic strength of I = 300 mM). The particular case of the 155 nm particles agglomerating at an ionic strength of I = 300 mM has been studied carefully since only approximate extrapolations of the particle potentials are available. Using a particle potential of −30 mV as the plateau value at large pH conditions [determined using the formula (9)], the stability ratio has been found to be equal to 1.0 over the whole range of pH values (not shown). However, when the value of the plateau of particle potentials at high pH conditions is decreased to −35 mV, particles were found to be stable at pH values larger than 6.5 (similar results have been obtained with a plateau value lower than −35 mV), as displayed in Figure 10. Therefore, in the following, the zeta potential of 155 nm carboxyl latex particles at 300 mM is considered to be equal to −35 mV. The slope of the stability ratio with pH values is now investigated numerically. For that purpose, the influence of surface roughness on particle agglomeration is assessed. Indeed, it has been seen in Effect of Agglomerate Geometry on the Stability Ratio that the slope of the stability ratio with respect to the solution conditions can be explained by the morphology of

Figure 12. Experimental data (points) and numerical values (lines) for the stability ratio (W) of rough 155 nm particles (5 nm asperities covering 0.015% of the surface) at various pH values for different ionic strengths: 2 (black Δ), 10 (red ●), 100 (green □), and 300 mM (blue ▽). 13702


Langmuir

Article

Figure 13. Evolution of the stability ratio W of 155 nm particles at I = 300 mM with asperities of 5 nm present on particles with various surface coverage.

especially at I = 300 mM, where a kind of plateau can be observed for pH values larger than 7. It should be noted that similar results are obtained for a plateau value from −40 mV up to −50 mV. To further explore the impact of surface roughness on the stability ratio, numerical simulations have been performed using 155 nm particles, an ionic strength of 300 mM, and a zeta potential of −35 mV with various roughness characteristics. First, to study the influence of the surface coverage on particle agglomeration, the radius of asperities has been considered to be constant and only the surface coverage has been changed (see Figure 13). It can be seen that a surface coverage of 5.0% yields the lowest stability ratio (i.e., the lowest average energy barrier) and, in that case, the average number of asperities present between the two colliding particles is roughly equal to 3. For lower surface coverage, the probability to have at least one particle−asperity contact is reduced and the average energy barrier is thus increased (due to the high energy barrier between smooth particles). For higher surface coverage, the energy barrier is increasing due to multiple contacts with asperities (the higher the surface coverage, the higher the number of asperties playing a role in the interaction). Second, to evaluate the impact of the asperity radius on the stability ratio, various asperity sizes have been used in the calculation while the average number of asperity−particle contacts has been kept constant (here a value of 3.1 has been chosen). As shown in Figure 14, the presence of an asperity between the two colliding particles generally reduces the energy barrier encountered. This is due to the fact that, when an asperity lies between particles, the separation distance between particles bulk bodies is higher, leading to smaller particle−particle interactions. The interaction energy also depends on asperity− particle interactions, which are proportional to the Derjaguin prefactor RpartRasp/(Rpart + Rasp). Therefore, the larger the asperity size, the higher the stability ratio (except for very small subnanometer asperities where the bulk particle−particle distance is not high enough to lead to significant reductions of particle−particle interactions). From this numerical evaluation of the impact of surface roughness on the stability of latex nanoparticles, it can be seen that small protrusions (between 1 and 5 nm) covering a small amount of the surface (less than 1%) can play a significant role in the stability of particles. Therefore, particles that may appear as nearly spherical on TEM images can still exhibit unexpected

Figure 14. Evolution of the stability ratio (W) of 155 nm particles at I = 300 mM with the size of asperities present on particles (3.1 particle− asperity contacts are considered on average).

stability behavior related to small scale surface roughness (here nanometer-scale asperities). To estimate the overall effect of particle roughness on agglomeration rates, another numerical calculation has been performed for 155 nm particles using asperities of 5 nm (polydispersion in asperity size has been introduced with a standard deviation of 5 nm) covering 0.05% of the surface of particles. As shown in Figure 15, numerical results are in very good agreement with experimental data. It should be noted that similar results can be obtained with slightly different roughness characteristics that are still in the range of surface roughness measured in other experiments (1 nm for instance), confirming the dependence of particle agglomeration on realistic surface roughness characteristics (and not only to a very narrow range of surface roughness that may be questionable). In addition, we have checked that results are independent of the time step used in the simulation. For that purpose, calculations have been performed using time steps 10 times lower, 10 times larger, and 100 times larger than the reference time step. The results are plotted in Figure 16, showing that numerical results are independent of the time step. Summary. These promising results highlight the fact that relatively accurate numerical predictions of particle stability can be obtained when the presence of surface roughness on particles is properly accounted for (the roughness characteristics chosen here appear plausible with other experimental 13703


Langmuir

Article

32, and 12 nm), respecting the particle concentrations (16 mg/ L for 12 nm particles, 21 mg/L for 32 nm particles, and 44 mg/ L for 65 nm particles). Particles were assumed to undergo Brownian motions (at room temperature T ∼ 296 K). The time step has been chosen so as to respect the condition σx ≫ dp (here, a value of σx = 5dp was used), and the size of the box was taken as L = 100σx. Particle−particle DLVO interactions have been calculated using a value of the Hamaker constant of 2.0 × 10−20 J (Hamaker constants between 1.0 × 10−20 and 5.0 × 10−20 J have been reported for this system31,32). The zeta potentials of hematite particles at various solution ionic strengths and pH conditions have been measured experimentally (using a laser Doppler velocimetry setup31), and the average particle potentials used in the calculations respect the measured values (average and standard deviation of zeta potentials). Particle zeta potentials at an ionic strength of 40 mM have been extrapolated from available experimental data using eq 9 (see ref 24). Numerical Results. In a first attempt to evaluate the stability ratio with various solution conditions (pH and ionic strength), numerical calculations have been performed with the simplified approach, where agglomerating particles are kept separated. The numerical predictions for the evolution of the stability ratio with the ionic strength (at a pH of 5.7) are displayed in Figure 17: it can be seen that the agreement with experimental data is

Figure 15. Experimental data (points) and numerical values (lines) for the stability ratio (W) of rough 155 nm particles (5 nm asperities covering 0.05% of the surface) at various pH for different ionic strengths: 2 (black Δ), 10 (red ●), 100 (green □), and 300 mM (blue ▽).

Figure 16. Numerical values for the stability ratio (W) of rough 155 nm particles (5 nm asperities covering 0.05% of the surface) at various pH for different ionic strengths at various time steps [Δt/10 (dashed line), 10Δt (dashed dotted line), 100Δt (dotted line)]: 2 (black Δ), 10 (red ●), 100 (green □), and 300 mM (blue ▽).

data). This case underlines another advantage of the present approach: the effect of small nanoscale heterogeneities can be included in the modeling approach while allowing tractable simulations to be performed with good precision. 3.2.3. Evolution of Particle Agglomeration with Particle Concentration. Experimental Data. He et al.32 studied the agglomeration of hematite particles at various ionic strengths (ranging from 1 to 100 mM), various pH conditions (ranging from 5.7 to 9) for three different particle radii (65, 32, and 12 nm) and various particle concentrations. The size of particles was recorded at various times using light scattering techniques, and the stability ratios were also extracted from the experimental data. This experiment is, to the authors’ knowledge, one of the first to provide insight into the effect of particle concentration, solution conditions (pH and ionic strength), and particle size on both the stability ratio and on the evolution of particle size with time. Thus, even though some inconsistencies have been found in the original paper, experimental results are worth studying numerically (especially the dependence on particle concentration and the evolution of the average particle size with time). Numerical Parameters. Numerical calculations have been performed for the three types of monodispersed particles (65,

Figure 17. Experimental data (points) and numerical values (lines with crosses) for the stability ratio (W) of hematite particles at various ionic strengths for different particle radius: 65 (black Δ), 32 (red □), and 12 nm (green ○).

satisfactory for all particles over the whole range of ionic strengths investigated. It should be noted that some of the calculations of particle−particle interactions have been performed using the formula from Verwey and Overbeek,9 valid for κR < 5 (the EDL formula used in the present paper is indeed valid provided that κR > 5). Following this first attempt, further numerical calculations have been performed, taking into account particle agglomeration (the size of particles is increased when the criterion Ekin > Ebarr is respected). For that purpose, the motions of 912 particles of 65 nm were simulated in a 5 μm box. A first calculation has been performed at a solution ionic strength of 100 mM at pH 5.7 and with a concentration of 44 mg/L. Numerical results are compared to experimental data from He et al.32 (see blue continuous curve in Figure 18): it can be seen that the initial evolution of the mean particle diameter is wellcaptured by the model, but after a while, numerical results are highly influenced by statistical noise. In fact, since only N0 13704


Langmuir

Article

In particular, since we start with a monodispersed set of particles, the initial evolution of the mean radius is less sensitive to the uncertainties mentioned above and, therefore, the fact that the numerical curve reproduces very correctly, the experimental evolution can be regarded as quite satisfactory. At longer times the numerical prediction of the mean radius is somewhat at the lower level of experimental values (which is in line with the previous explanations) but still within the experimental range, showing that the trend is well-retrieved. Further assessments are also provided by considering various solution conditions and particle concentrations. The key result obtained by He et al.32 is the study of the influence of particle concentration on particle agglomeration. They measured the evolution of the mean radius of 65 nm particles at an ionic strength of 50 mM and a pH of 5.7 for various particle concentrations. However, an inconsistency has been found in the original article: particle concentration was said to be equal to 132 mg/L in the figure caption, whereas a value of 44 mg/L was indicated elsewhere. In the following, we have chosen to perform numerical calculations at the particle concentration provided and discussed within the text (44 mg/ L). To investigate numerically the effect of particle concentration on agglomeration, a second calculation has been made for a concentration of 14.6 mg/L (which respects the ratio of particle concentrations proposed in the figure caption of the original article). As shown in Figure 19, numerical predictions

Figure 18. Experimental data at 100 (cyan ○) and 50 mM (violet △) compared to numerical values (straight line) for the stability ratio (W) of hematite particles at pH 5.7 for 65 nm particles: calculation with N0 = 912 initial particles at I = 100 mM (blue line), averaged over 20 simulations with continuations of the first calculation with a domain twice larger at I = 100 mM (red dashed line), with an average at I = 50 mM (orange dashed line).

particles have been injected initially, the maximum particle size that can be reached is the one corresponding to an agglomerate of N0 particles. To remove this limitation, a second calculation has been performed: as soon as a single particle remains in the simulation, another calculation is started, considering eight such boxes which form a new simulation domain with a length that is twice longer. This numerical treatment provides approximate information on the statistics of particle agglomeration at longer times. It indeed allows one to perform calculations of the formation of large agglomerates with a constant reasonable CPU cost, while preventing the use of a very high number of particles as in a single calculation. This procedure has been repeated 20 times using different random numbers for each simulation, and average results are then extracted so as to obtain a reliable statistical prediction. The results are plotted in Figure 18, where the dotted smoothed curve (in red) corresponds to the average over the 20 simulations, showing a better agreement with experimental data at large time steps. The second set of data (experiments in purple and dotted curve) corresponds to the same calculation but at an ionic strength of 50 mM, which also shows good agreement between experiments and predictions. At this stage, it is worth recalling that there are still some uncertainties about whether the experimental and the numerical mean radii can be directly collated. Indeed, the experimental particle radii are obtained through light-scattering techniques which should be used carefully with polydispersed spherical suspensions due to their sensitivity to polydispersity in sizes (larger particles have a much larger scattering intensity than smaller ones). In contrast, numerical results are sensitive to the choice of the law retained to estimate the diameter, or radius, of an agglomerate from the radius of the primary particles, which make up this agglomerate (see the discussion in Treatment of the Collision). For example, since no internal porosity has been accounted for in the formation of large aggregates, the numerical estimations based on a spherical assumption and on the additivity of elementary particle volumes can lead to somewhat lower predictions for the radius of actual large aggregates. The comparison shown in Figure 18 and in the following figures should then be assessed with some care. Nevertheless, these comparisons do already reveal interesting quantitative trends.

Figure 19. Experimental data (points) and numerical values (lines) for the mean radius of 65 nm particles at various particle concentrations at pH 5.7 and I = 50 mM: 44 (blue △) and 14.67 mg/L (green □).

reproduce the experimental trends: the slower increase of the mean particle radius with decreasing particle concentration is well-captured in the simulations. Summary. The stochastic approach for particle agglomeration not only provides accurate predictions of the stability ratio but also provides correct evaluations of the evolution of particle sizes with time for various solution conditions and particle properties (zeta potentials and concentration). Therefore, the present method allows accurate simulations of particle agglomeration to be performed at reasonable computational costs.

4. CONCLUSION In this paper, a new stochastic numerical approach for the agglomeration of colloidal particles has been developed in the framework of particle-tracking approaches in flows. For that purpose, particle agglomeration has been modeled as the 13705


Langmuir

Article

related to particle agglomeration will pave the way for future studies of agglomerate fragmentation.

coupling between the collision step (transport of particles by the flow) and the attachment step (adhesion of particles). Following existing works on particle collisions (more details in ref 15 and in a forthcoming paper16), the transport step has been modeled using a stochastic approach that is valid in the diffusive regime (i.e., for time steps much larger than the particle relaxation time). Within this approach, the detection and treatment of particle collisions is based on the notion of continuous particle trajectories within a time step. The attachment step has been modeled using the DLVO theory, which describes the interaction between two particles as the sum of van der Waals and electrostatic forces. Following previous works on multilayer deposition,4,22 the coupling between the collision step and the attachment step is based on an energetic approach: particles adhere together if their relative kinetic energy is larger than the energy barrier encountered. The new specific aspects of the present paper are, therefore, a complete stochastic description and detailed simulations of the full agglomeration process, where agglomeration kernels are obtained as results. Numerical results have been compared to various experimental and theoretical data for the agglomeration of particles undergoing Brownian motions. Simulations were shown to provide satisfactory predictions of the stability ratio for various solution conditions (both pH and ionic strength) as well as accurate evaluations of the evolution of particle sizes with time and with particle concentration. In addition, the influence of particle roughness as well as agglomerate morphology on the stability ratio has been assessed using simplified models accounting for such local nanoscale properties. The good results obtained underline the fact that the present stochastic approach allows one to perform tractable simulations of particle agglomeration, even in the presence of small-scale heterogeneities. The promising results obtained with the present approach open interesting possibilities for the simulation of particle agglomeration over the whole range of particle sizes (even for nanoscopic particles) and time steps (from a ballistic to a diffusive regime) with a numerical method that induces reasonable computational costs. It is indeed believed that the approach developed here is applicable to more complex flows (such as nonhomogeneous flows or turbulent flows) when the fluid velocity is properly accounted for in eq 2 (for instance using DNS calculations of the fluid flow). As a result, the resulting stochastic approach can find potential applications in several fields since numerical predictions of the collision kernels obtained with the approach can be used as an input in classical PBE equations. Nevertheless, several issues still remain to be investigated in order to obtain a complete numerical model to simulate the whole process of particle agglomeration. Some of these issues are listed in the following. (1) Morphology of agglomerates: in the present paper, the size of the agglomerate formed after collision was chosen so that the apparent volume of the spherical agglomerate is equal to the sum of the volumes of the two colliding particles. This amounts to considering a fractal dimension of 3. However, the morphology of agglomerates has been shown to depend both on the hydrodynamic conditions and on interparticle forces.33 Such effects should thus be properly modeled in later simulations. (2) Fragmentation: fragmentation of agglomerates has not been considered in the present study. Yet, with the promising results obtained with the present approach, the level of understanding of the issues

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Maximova, N.; Dahl, O. Environmental implications of aggregation phenomena: Current understanding. Curr. Opin. Colloid Interface Sci. 2006, 11, 246−266. (2) Taboada-Serrano, P.; Chin, C.; Yiacoumi, S.; Tsouris, C. Modelling aggregation of colloidal particles. Curr. Opin. Colloid Interface Sci. 2005, 10, 132−132. (3) Dukhin, S.; Zhu, C.; Dave, R.; Yu, Q. Hydrodynamic fragmentation of nanoparticle aggregates at orthokinetic coagulation. Adv. Colloid Interface Sci. 2005, 114−155, 119−131. (4) Henry, C.; Minier, J.-P.; Lefèvre, G. Towards a description of particulate fouling: from single particle deposition to clogging. Adv. Colloid Interface Sci. 2012, 185−186, 34−76. (5) Derjaguin, B.; Landau, L. Theory of the stability of strongly charged lyophobic sols and of the adhesion of strongly charged particles in solution of electrolytes. Acta Physicochim. URSS 1941, 14, 633−662. (6) Verwey, E.; Overbeek, J. T. G. Theory of Stability of Lyophobic Colloids; Elsevier: Amsterdam, 1948. (7) Elimelech, M.; Gregory, J.; Jia, X.; Williams, R. Particle Deposition and Aggregation: Measurement, Modelling and Simulation; Butterworth Heinemann: Oxford, 1995. (8) Parsegian, V. A. van der Waals Forces: A Handbook for Biologists, Chemists, Engineers, and Physicists; Cambridge University Press: Cambridge, 2005. (9) Hunter, R. J. Foundations of Colloid Science, 2 ed.; Oxford University Press: Oxford, 2001. (10) Israelachvili, J. N. Intermolecular & Surface forces, 3 ed.; Academic Press: San Diego, 2011. (11) von Smoluchowski, M. Experiments on a mathematical theory of kinetic coagulation of colloid solutions. Z. Phys. Chem. 1917, 92, 129−168. (12) Chen, M.; Kontomaris, K.; McLaughlin, J. Direct numerical simulation of droplet collisions in a turbulent channel flow. Part I: collision algorithm. Int. J. Multiphase Flow 1998, 19, 1079−1103. (13) Sundaram, S.; Collins, L. R. Numerical considerations in simulating a turbulent suspension of finite-volume particles. J. Comput. Phys. 1996, 124, 337−350. (14) Sigurgeirsson, H.; Stuart, A.; Wan, W. Algorithms for particlefield simulation with collisions. J. Comput. Phys. 2001, 172, 766−807. (15) Mohaupt, M.; Minier, J.-P.; Tanière, A. A new approach for the detection of particle interactions for large-inertia and colloidal particles in a turbulent flow. Int. J. Multiphase Flow 2011, 37, 746−755. (16) Henry, C.; Minier, J.-P.; Mohaupt, M.; Profeta, C.; Pozorski, J.; Tanière, A. A stochastic approach for the simulation of collisions between colloidal particles at large time steps. Submitted to Int. J. Multiphase Flow, 2013. (17) Minier, J.-P.; Peirano, E.; Chibbaro, S. Weak first and second order numerical schemes for stochastic differential equations appearing in Lagrangian two-phase flow modelling. Monte Carlo Methods and Applications 2003, 9, 93−133. (18) Peirano, E.; Chibbaro, S.; Pozorski, J.; Minier, J.-P. Mean-field/ PDF numerical approach for polydispersed turbulent two-phase flows. Prog. Energy Combust. Sci. 2006, 32, 315−371. (19) Borodin, A.; Salminen, P. Handbook of Brownian Motion: Facts and Formulae; Ed. Birkhauser: Basel, 2002. (20) Liang, Y.; Hilal, N.; Langston, P.; Starov, V. Interactions forces between colloidal particles in liquid: Theory and experiment. Adv. Colloid Interface Sci. 2007, 134−135, 151−166. 13706


Langmuir

Article

(21) Gregory, J. Approximate expressions for retarded van der Waals interaction. J. Colloid Interface Sci. 1981, 83 (22) Henry, C.; Minier, J.-P.; Lefèvre, G.; Hurisse, O. Numerical study on the deposition rate of hematite particles on polypropylene walls: Role of surface roughness. Langmuir 2011, 27, 4603−4612. (23) Henry, C.; Minier, J.-P.; Lefèvre, G. Numerical study on the adhesion and reentrainment of nondeformable particles on surfaces: the role of surface roughness and electrostatic forces. Langmuir 2012, 28, 438−452. (24) Snoswell, D.; Duan, J.; Fornasiero, D.; Ralston, J. Colloid stability of synthetic titania and the influence of surface roughness. J. Colloid Interface Sci. 2005, 286, 526−535. (25) Derksen, J. Direct numerical simulations of aggregation of monosized spherical particles in homogeneous isotropic turbulence. AIChE J. 2012, 58, 2589−2600. (26) Babick, F.; Schießl, K.; Stintz, M. van-der-Waals interaction between two fractal aggregates. Adv. Powder Technol. 2011, 22, 220− 225. (27) Schießl, K.; Babick, F.; Stintz, M. Calculation of double layer interaction between colloidal aggregate. Adv. Powder Technol. 2012, 23, 139−147. (28) Bouhaik, I. S.; Leroy, P.; Ollivier, P.; Azaroual, M.; Mercury, L. Influence of surface conductivity on the apparent TiO2 nanoparticles: Application to the modeling of their aggregation kinetics. J. Colloid Interface Sci. 2013, 406, 75−85. (29) Behrens, S.; Christl, D.; Emmerzael, R.; Schurtenberger, P.; Borkovec, M. Charging and aggregation properties of carboxyl latex particles: Experiments versus DLVO theory. Langmuir 2000, 16, 2566−2575. (30) Popa, I.; Sinha, P.; Finessi, M.; Maroni, P.; Papastavrou, G.; Borkovec, M. Importance of charge regulation in attractive doublelayer forces between dissimilar surfaces. Phys. Rev. Lett. 2010, 104, 228301. (31) Schudel, M.; Behrens, S.; Holthoff, H.; Kretzschmar, R.; Borkovec, M. Absolute aggregation rate constants of Hematite particles in aqueous suspensions: A comparison of two different surface morphologies. J. Colloid Interface Sci. 1997, 196, 241−253. (32) He, Y.; Wan, J.; Tokunaga, T. Kinetic stability of hematite nanoparticles: The effect of particle sizes. J. Nanopart. Res. 2008, 10, 321−332. (33) Kim, S.; Lee, K.-S.; Zachariah, M.; Lee, D. Three-dimensional off-lattice Monte Carlo simulations on a direct relation between experimental process parameters and fractal dimension of colloidal aggregates. J. Colloid Interface Sci. 2010, 344, 353−361.

13707


A New Stochastic Approach for the Simulation of ... - ACS Publications

Recommend Documents