Influence of the Potential Well on the Breakage ... - ACS Publications

Article pubs.acs.org/Langmuir

Influence of the Potential Well on the Breakage Rate of Colloidal Aggregates in Simple Shear and Uniaxial Extensional Flows Zhiqiang Ren,† Yogesh M. Harshe,‡ and Marco Lattuada*,† †

Adolphe Merkle Institute, University of Fribourg, Chemin des Verdiers 4, 1700-Fribourg, Switzerland Department of Chemical Engineering, Delft University of Technology, Julianalaan 136, 2628 BL Delft, The Netherlands

‡

S Supporting Information *

ABSTRACT: In this work we build on our previous paper (Harshe, Y. M.; Lattuada, M. Langmuir 2012, 28, 283−292) and compute the breakage rate of colloidal aggregates under the effect of shear forces by means of Stokesian dynamics simulations. A library of clusters made of identical spherical particles covering a broad range of masses and fractal dimension values (from 1.8 to 3.0) was generated by means of a combination of several Monte Carlo methods. DLVO theory has been used to describe the interparticle interactions, and contact forces have been introduced by means of the discrete element method. The aggregate breakage process was investigated by exposing them to well-defined shear forces, generated under both simple shear and uniaxial extensional flow conditions, and by recording the time required to reach the first breakage event. It has been found that the breakage rate of clusters was controlled by the potential well between particles as described by DLVO theory. A semiempirical Arrhenius-type exponential equation that relates the potential well to the breakage rate has been used to fit the simulation results. The dependence of the breakage process on the radius of gyration, on the external shear strength, and on the fractal dimension has been obtained, providing a very general relationship for the breakage rate of clusters. It was also found that the fragment mass distribution is insensitive to the presence of electrostatic repulsive interactions. We also clarify the physical reason for the large difference in the breakage rate of clusters between simple shear and the uniaxial extensional flow using a criterion based on the energy dissipation rate. Finally, in order to answer the question of the minimum cluster size that can break under simple shear conditions, a critical rotation number has been introduced, expressing the maximum number of rotations that a cluster exposed to simple shear could sustain before breakage.

1. INTRODUCTION Controlling the formation of clusters in colloidal suspensions is an aspect of utmost importance in many industrial processes. Guaranteeing the stability of colloids in suspensions can sometimes be a difficult task, and it represents a challenge from both an experimental and a theoretical point of view.1,2 Introducing surface charges to provide electrostatic repulsion among the suspended particles is a common and practical way to control the stability of colloidal suspensions.3 There are situations of practical relevance where the destabilization of a colloidal suspension is a required step in their processing. In the case of electrostatically stabilized particles, the range of electrostatic repulsions can be controlled through the addition of electrolytes, which compress the electric double layer surrounding them. When the range of repulsive interactions © 2015 American Chemical Society

decreases below a critical value, the attractive force, usually of van der Waals origin, may gain the upper hand and the particles will start to aggregate, generating colloidal clusters.4 In industrial processes in which the aggregation of colloids is performed, the process is carried out in the presence of shear forces. Shear accelerates aggregation5−8 but can also overcome the attraction force between particles and break clusters into fragments.8,9 Often a dynamic balance between the aggregation and the breakup is established.10−12 The investigation of such a balance is extremely challenging13−16 because both the growth17 and the breakage phenomena18 of clusters need to Received: December 24, 2014 Revised: April 24, 2015 Published: May 5, 2015 5712

DOI: 10.1021/la504966y Langmuir 2015, 31, 5712−5721

Article

Langmuir

interparticle interactions (modeled trough DLVO theory), and applied shear rate (simple shear γ̇s = 20 000−100 000 s−1; uniaxial extensional flow γ̇e = 5000−60 000 s−1) has been investigated. The effect of repulsive interactions on the fragment mass distribution was also investigated. To obtain a general mathematical expression to predict the breakage rate, a wide range of simulations have been conducted to obtain statistically significant results. As a result, the final expression in this paper can be used to investigate the clusters in two of the most common flow fields used, i.e., uniaxial extensional flow and simple shear. While the former has been suggested to represent what happens below the Kolmogorov length scale under turbulent conditions,34 the latter is the most commonly used in controlled experiments. A second objective of this work is to further clarify the differences in breakage rate between the two different types of flow. To achieve this objective, an approach based on the energy dissipation rate was used to compare different effects of simple shear and the uniaxial extensional flow on the breakage of a cluster, and their relationships were interpreted from an energy point of view. Finally, a novel approach to obtaining the critical mass of clusters exposed to simple shear forces has been proposed by introducing a critical rotation number, i.e., expressing the maximum number of rotations that a cluster can withstand before breaking.

be taken into account at the same time. While shear effects on aggregation rate are relatively well established,19 their role in the breakage process is less well understood. A rigorous analysis of the process from a rigorous standpoint is very challenging because of the strong intertwining of the structure of the clusters, the pattern of the external shear, and the interactions between particles, including hydrodynamic and van de Walls attraction and electrostatic repulsion. Detailed simulations are under these conditions of great help. Many papers have been published on the investigation of the breakage process through both theoretical modeling20−23 and controlled experiments. 24−26 Even though under some conditions good agreement has been obtained between modeling and experiments,25−28 the detailed relationship between the breakage rate and external conditions is still to be fully unraveled. Previous research mostly focused on the dynamics of restructuring and on the long-term breakage behavior of clusters to retrieve the final distributions of fragments. However, the precise evaluation of the breakage rate under well-defined conditions has seldom been discussed. Higashitani et al.29 investigated the breakage process with a 2-D model implementing the discrete element method to account for the force balance. They proposed a power-law dependence of the size of the fragment produced from an aggregate breakup on the magnitude of the external shear. Inamuro et al.30 investigated the cluster breakage under external shear flows using the lattice Boltzmann method and found that the ratio of the hydrodynamic drag force to the cohesive force evidently has an influence on the breakage behavior. Nishiyama et al.31 further investigated the breakage phenomena through the lattice Boltzmann method combined with Brownian motion and found that ratio of the fluid force to the maximum interparticle force is the key factor in determining the breakage. Harada et al.2 also investigated the cluster breakage through Stokesian dynamics. The restructuring phenomena of the clusters were simulated, and a power law model describing the evolution of the coordination number with time and applied shear rate was proposed. In our previous work,32 we used Stokesian dynamics to investigate the breakage rate of a cluster, obtained as the inverse of the time required for the first breakage event to occur, and found a power relationship among the breakage rate of clusters, the magnitude of applied shear, and the radius of gyration. But the interactions between particles were assumed to have only an attractive van der Waals component, while the electrostatic repulsion was not accounted for in the simulations. Recent theoretical calculations suggest that electrostatic repulsive interactions affect the breakage rate exponentially.33 Because of the intrinsic complexity of the underlying physics, a general mathematical expression which relates the breakage rate of clusters to their geometrical features, to the interparticle interactions, and to the properties of the surrounding fluid is lacking. This work aims at filling the existing gap in that the breakage rate of colloidal clusters with different fractal dimension values,25 whose particles (assumed to have a radius of 1 μm) are subject to both attractive and repulsive interactions according to DLVO theory, has been investigated under simple shear and uniaxial extensional flows. A semiempirical exponential equation has been proposed on the basis of the results from Stokesian dynamics simulations of the breakage process. The dependence of the breakage rate on the aggregate size (Nsphere = 30−215), fractal dimension (df = 1.8−3.0),

2. SIMULATION DETAILS 2.1. Hydrodynamic Interactions. In the simulations carried out in the present work, hydrodynamic interactions were modeled by means of Stokesian dynamics (SD).35 This method has been developed to model interactions among a large number of particles dispersed in a Newtonian fluid exposed to linear flow fields with a very low particle Reynolds number (Rep ≪ 1). Compared to common finite element as well as lattice Boltzmann methods,13 SD can obtain much faster and accurate results. The core of this method is the evaluation of the grand resistance matrix, an 11Nsphere × 11Nsphere matrix, which provides a linear relationship between the forces − torques − stresslets and the translational velocities − angular velocities − shear rate of (acting on) the particles: ⎛ Fn̅ ⎞ ⎛ Un̅ − U̅ ∞ ⎞ ⎜ ⎟ ⎜ ⎟ ⎜Tn̅ ⎟ = η R̿ ·⎜ Ω̅ − Ω̅ ∞⎟ f ⎜ ⎟ ⎜⎜ n ⎟⎟ ⎜ ⎟ ⎝ −E̿ ∞ ⎠ ⎝ Sn̿ ⎠

(1)

Here, R̿ is the grand resistance matrix, which is built on the basis of the multipole expansion of the integral representation of Stokes equations, and ηf is the viscosity of the fluid. The forces (F̅n), torques (T̅ n), and stresslets (Sn̿ ) acting on each particle are a linear combination of the relative translational (U̅ n − U̅ ∞) and the relative angular (Ω̅ n − Ω̅ ∞) velocity of all particles with respect to the undisturbed fluid velocity computed at their centers and the rate of strain tensor (E̿ ∞) of the undisturbed fluid. The resistance matrix R̿ and its inverse mobility matrix M̿ contain all of the information about hydrodynamic interactions among the particles in the system. Durlofsky et al.36 developed explicit expressions for all of the matrix elements. Stokesian dynamics is based on a finite multipole expansion of the integral representation of Stokes equations, and a truncation of the expansion was made in building the matrices. The truncated expressions cannot include both the short-range 5713

DOI: 10.1021/la504966y Langmuir 2015, 31, 5712−5721

Article

Langmuir

were also included in the simulation. The tangential force expression was reported by Pantina and Furst38 in their experimental work for two particles in close contact. A more detailed mathematical formulation accounting for the contact tangential forces (F̅j) and critical bending moment (T̅ j) is proposed by Becker et al.:39

lubrication interactions and the long-term many-body interactions at the same time. Durlofsky et al.36 therefore proposed the following approximated equation to include both long- and short-range interactions together in one matrix ∞

∞

R̿ = (M̿ )−1 + R̿ 2B − R̿ 2B

(2)

where M̿ ∞ is the far-field mobility matrix, R̿ 2B is the exact twobody resistance matrix, and R̿ ∞ 2B is the inverse of the two-body far-field mobility matrix including only the first-order terms, added pairwise. By doing this, both lubrication interactions and long-range many-body interactions were included in the R̿ matrix, which can then be used to accurately model particle suspensions. 2.2. Particle Interactions. When two charged particles in aqueous solution approach each other, their interaction potential can be described by means of the Derjaguin− Landau−Verwey−Overbeek (DLVO) theory. The potential energy of interactions (V) between two particles is split into three contributions: the van der Waals attraction, the electrostatic repulsion, and the short-range Born repulsion.3 Feke et al.37 gave a general equation which can be used to calculate the overall potential between two spherical particles,

Fj̅ = K t(ξij̅ − ξji̅ )Tj̅ = 2R pK t nij̅ ξij̅

Here Kt is the tangential rigidity, ξ̅ij is the spring elongation, which is the tangential displacement vector between connected particles i and j, and ni̅ j is the unit normal vector between the centers of particles i and j. The spring is considered to be broken when the spring elongation ξ̅ exceeds its critical value ξ̅max, estimated from the critical bending moment, in which case sliding occurs. 2.3. Clusters Library. The aggregate geometries used in the simulations were generated by a combination of different Monte Carlo methods, as explained in previous work.25 There are three important parameters defining each cluster: the number of particles in Nsphere, the radius of gyration Rg, and the fractal dimension df. The radius of gyration is a measure of the cluster size, and the fractal dimension is related to the density of the cluster. The fractal dimension is typically used to define statistically self-similar objects.40 In the case of clusters, its values range from 1.0 for a straight chain of particles to 3.0 for a Euclidean object. However, it should be noted that the clusters used in our simulations are too small to truly show selfsimilarity. Therefore, the meaning of the fractal dimension is that of a scaling exponent relating cluster size to its mass. Low fractal dimension values indicate open structures and are typical of clusters obtained under diffusion-limited or reaction-limited stagnant conditions, whereas high values mean dense structures and are usually encountered in shear-driven aggregation. Images of typical clusters used in this work are shown in Figure S9. The scaling relationship is shown below:

−1 ⎡ 2 V s2 − 4 ⎤ 2 = + 2 + ln ⎢ 2 ⎥ Ah 6 ⎣s − 4 s s2 ⎦ + NR

Fcψ0

( ) exp(−κ(s − 2))

tanh2

κ

4RT 2

− 2s 2 + 60 s 2 + 14s + 54 ⎤ 1 ⎡ s − 14s + 54 + N12 ⎢ + + ⎥ 7 7 s ⎣ (s − 2) s (s + 2)7 ⎦ 2

(3)

Here s is the spheres’ center-to-center distance made dimensionless by the particle radius, while Ah is the Hamaker constant, ψ0 is the potential on the surface, R is the ideal gas constant, Fc is the Faraday constant, and T is the absolute temperature in Kelvin. There are three main parameters in this equation. The first parameter NR measures the ratio of the electrostatic forces to the van der Waals forces, NR =

Nsphere

64πR p3RTC0 Ah

(5)

where σ12 is the separation distance at which the total potential is zero, which ranges from 0.208 to 0.5 nm in the simulations. The third parameter κ is the ratio of the particle radius to the diffuse double-layer thickness: 1 1 = κ Rp

εrε0RT 2Fc 2C0

⎛ R g ⎞d f = k f ⎜⎜ ⎟⎟ ⎝ Rp ⎠

(8)

Here kf is a packing coefficient, which is normally close to 1. More details can be found in the literature.25 In this work, clusters with a broad range of the number of particles and fractal dimensions were used from the clusters library for breakage simulations. In the clusters library, open clusters (df = 1.8 and 2.1) were generated by a cluster−cluster aggregation algorithm, relatively dense clusters (2.1 < df < = 2.5) have been generated with a tunable fractal dimension algorithm starting from open clusters, and very dense clusters (2.5 < df < = 3.0) have been produced by a densification algorithm starting from clusters with df = 2.5. 2.4. Flow Field. In our simulations, both a simple shear flow field and a uniaxial extensional flow field have been considered. The simple shear flow was created by imposing a velocity gradient in the x and y directions

(4)

where Rp is the particle radius and C0 is the concentration of electrolyte in the fluid. The second parameter N12 expresses the ratio of the Born repulsive forces to the van der Waals forces, 6 4! ⎛⎜ σ12 ⎞⎟ N12 = 4 10! ⎜⎝ R p ⎟⎠

(7)

⎡0 1 0⎤ ⎢ ⎥ ∇us = γṡ ⎢ 0 0 0 ⎥ ⎣0 0 0⎦

(6)

Here, ε0 is the vacuum permittivity, and εr is the relative static permittivity of water. Besides these physical−chemical interactions, mechanical contact forces, simulated using the discrete element method,

(9)

where γ̇s the simple shear rate and suffix s indicates simple shear. 5714

DOI: 10.1021/la504966y Langmuir 2015, 31, 5712−5721

Article

Langmuir The uniaxial extensional flow was imposed as a compressive deformation in the y and z directions and a uniaxial extensional deformation in the x direction ⎡2 0 0 ⎤ ⎢ ⎥ ∇ue = γė ⎢ 0 −1 0 ⎥ ⎣ 0 0 −1⎦

Here τi is the breakage time for cluster i, n is the total cluster number, and kb is the breakage rate averaged over all clusters. The values of the simulation parameters used are listed in Table 1. Table 1. Set of Parameters Used in the Simulations

(10)

Rp = 1 μma ηf = 1.0 × 10−3 Pa.s Kt = 6.65 × 10−4 N/m kBO = 1.381 × 10−23 J/K

where γ̇e the shear rate of uniaxial extensional flow and subscript e indicates uniaxial extensional flow. In the present work we have investigated the following range of strain rates for simple shear and uniaxial extensional flows: γ̇s = 20 000−100 000 s−1; γ̇e = 5000−60 000 s−1. The range of the strain rates was so chosen to observe cluster breakage in physically realistic timeframes. Furthermore, for each simulation, the origin of the coordinate system was selected to be at the center of mass of the simulated aggregate. One should also notice that at the high shear rate values chosen all effects of Brownian motion have been neglected in the simulations. In fact, for the corresponding Peclet values, defined as

Pe =

a

(11)

where kBO is the Boltzmann constant and σ is the shear stress, for simple shear given by σs = ηfγ̇s and for uniaxial extensional flows to σe = ηfγ̇e, the Pe value is between 2.2 × 104 and 4.5 × 105 in our simulations. 2.5. Simulation Parameters. The procedure followed to perform the simulations closely follows what has been reported in our previous work.32 Before starting each simulation, an equilibration process was implemented to lower the clusters’ energy. During this equilibration process, no flow was applied, and only minimal particle displacements were observed. After the equilibration step, the flow field was imposed on a cluster and the actual Stokesian dynamics simulation was started. The velocities of all particles composing the cluster were obtained. The displacements of the particles were calculated by integrating the velocities by means of a fourth-order Runge− Kutta method, and very small time steps (10−9−10−8 s)41 were used to obtain good accuracy. To avoid particle overlap, a minimum particle distance of δmin = 0.3 nm was introduced into the simulations.41 If during the simulation the separation distance decreases below the δmin value, it is set back to the δmin value. The distance between particles was monitored throughout the complete duration of each simulation. A bond was considered to be broken when the interparticle distance between any two particles initially connected was found to reach a chosen limiting distance of 2 nm.32 The breakage of a cluster was determined by looking at the number of fragments through the analysis of all bonds. As soon as two fragments were present, the time of breakup was determined. The time of breakup was recorded for each simulation, and its reciprocal value provides the breakage rate. For each set of cluster morphology (df) and cluster mass (Nsphere), 100−300 clusters were randomly chosen from the clusters library, and for each simulated cluster, the breakage rate and produced fragments were determined. Then the obtained breakage rate was averaged over all clusters for the same cluster mass and structure.

kb =

1 n

n

∑ i=1

1 τi

Ah = 1.33 × 10−20 δmin = 0.3 nm εr = 80

Unless stated otherwise.

3. RESULTS AND DISCUSSION 3.1. Potential Well between Particles. The goal of this work was to perform detailed simulations of the breakage rate of fractal clusters composed of uniformly sized spherical particles. The interaction potential among particles was modeled by DLVO theory. In our previous work,32 we considered a similar problem, but only the attractive van der Waals interactions were accounted for. Here instead we primarily focused on the role played by the repulsive interactions. As mentioned in the previous section, clusters with different Nsphere and Rg were generated using Monte Carlo simulations and stored in a cluster library. Because the Monte Carlo algorithms do not explicitly consider interparticle interactions, clusters chosen directly from this library are not in a local minimum of potential energy at the start of the simulation. Therefore, an equilibration process was implemented before each simulation started to reduce the potential energy of clusters. The equilibration process was found not to change the relative positions of the particles significantly. Figure 1 shows a typical DLVO potential curve between two charged spherical particles in an aqueous medium, where the

6πσR p3 kBOT

Fc = 9.648 C/mol R = 8.314 J/(mol·K) ξmax = 24 nm ε0 = 8.854 × 10−12 F/m

Figure 1. Plot of dimensionless energy (V/kBOT) as a function of the interparticle separation for Rp = 1 μm, ψ0 = 0.03 V, C0 = 10 mol/m3.

contributions of van der Waals attraction, of electrostatic repulsion due to the electric double layer, and of short-range Born repulsion were added together. From Figure 1 we can observe that, as two particles approach, the overall energy initially increases until an energy peak (Vmax) is reached, which is due to the influence of electrostatic repulsion. The peak is followed by an energy minimum (Vmin) at a very short distance, usually referred to the primary minimum, before the Born repulsion starts to influence the potential energy. After the equilibration process, the nearest-neighbor particles inside a cluster will stay at Vmin to reduce the system energy. When the particles stay at Vmin, corresponding to a typical distance of less than 0.5 nm, it is assumed that the particles are

(12) 5715

DOI: 10.1021/la504966y Langmuir 2015, 31, 5712−5721

Article

Langmuir touching each other. It is further assumed that 2 nm, as discussed in our previous paper, has been taken as the criterion to identify particle separation because at this distance the overall force between a pair of particles has decayed to 2% of its maximum value. When a cluster is exposed to an external shear, it will start to deform and break into fragments if the applied shear is strong enough to provide sufficient energy to break bonds. The energy needed to separate a particle pair is given by the potential well (Vb). The minimum dimensionless energy required to break one particle−particle connection in a cluster is this potential well, given by Vb = (Vmax − Vmin)/kBOT

(13) Figure 3. kb vs Vb under uniaxial extensional flow for Nsphere = 60, df = 2.5, σe = 10 Pa, and four different values of the surface potential ψ0 = 0.02, 0.03, 0.04, and 0.05 V respectively. For each surface potential value ψ0, simulations for electrolyte concentrations C0 = 0.01, 0.1, 1, and 10 mol/m3 have been performed.

kBOT is used to make the energy dimensionless. Thus, the potential well calculated through DLVO theory should determine the cluster breakage rate. Finding a relationship between the potential well (Vb) and cluster breakage rate (kb) under an external shear has been attempted in a recent publication.33 An Arrhenius-type relationship is a good candidate for describing this relation because it describes how molecules change from one energy state to another by going across a potential well. It has also been used to model the rate of aggregation of particles subjected to repulsive interactions as well as the rate of escape of particles from a potential well.42 Thus in the following section we develop a relationship between Vb and kb following Arrhenius’ empirical formalism. 3.2. Relationship between Potential Well and Cluster Breakage Rate. Figures 2 and 3 show the relationship

between kb and Vb could be described by using an Arrhenius empirical formalism: d(ln k b) = −B d(Vb)

(14)

Here B is a constant for given values of Rg/Rp, df, and σ. B is also independent of ψ0 and C0. This suggests that the breakage process requires at least two particles to escape from the minimum of the potential energy curve and overcome the maximum of the well. The meaning of B is that of an inverse temperature in the Arrhenius equation, though much more complex. This B is a function of three basic parameters related to clusters geometry and to the external forces: Rg/Rp, df and σ. The details of how to uniquely estimate B is investigated in the following. As an example of the dependence of parameter B on the different physical quantities, we considered the case of clusters with df = 1.8 exposed to simple shear conditions. All other results in this paper were obtained using the same procedure. In Figure S1 (Supporting Information), Rg/Rp and df of the clusters were fixed, and the influence of σs was investigated. From Figure S1 we can observe that there is a liner relationship between ln kb and Vb when σs was varied from 20 to 100 Pa. Clearly the slope of each line is different, which means that d(ln kb)/d(Vb) changes when varying the external shear. Figure S2 shows instead the relationship between ln σs and (ln kb)/d(Vb) for different Rg/Rp values. It is apparent from Figure S2 that there is a linear relationship between these two quantities. The influence of the radius of gyration Rg/Rp on the breakage rate kb was investigated in a fashion similar to that for the external shear. The calculated results are plotted in Figures S3 and S4. In Figure S3 the external shear and the fractal dimension were fixed and Rg/Rp was changed from 6.6 to 14.9. We can see from Figure S3 that ln kb is also linear with respect to Vb for different Rg/Rp values, but again with different slopes. In Figure S4 the relationship between ln(Rg/Rp) and d(ln kb)/ d(Vb) was plotted, showing once more a linear relation. Combining the results in Figure S1−S4, a linear regression method was used to find the relationships among (ln kb)/ d(Vb), ln(Rg/Rp), and ln σ. This relation can be written in the following linear form,

Figure 2. kb vs Vb under simple shear for Nsphere = 60, df = 1.8, σs = 100 Pa, and four different values of the surface potential ψ0 = 0.02, 0.03, 0.04, and 0.05 V, respectively. For each surface potential value ψ0, simulations for electrolyte concentrations C0 = 0.01, 0.1, 1, and 10 mol/m3 have been performed.

between the potential well (Vb) and cluster breakage rate (kb) under both simple shear (Figure 2) and uniaxial extensional flow (Figure 3). In these two figures the potential well was calculated for different combinations of surface potential ψ0 values and electrolyte concentration C0. Figure S10 in SI shows Figure 3 redrawn including error bars, computed from the standard deviation of the breakage rate constants obtained from the population of clusters. The error bars show the intrinsic and non-negligible variability of the breakage rate constant values. First, we observe how the breakage rate decreases as the potential well increases. From these two figures we can observe that d(ln kb)/d(Vb) is a constant in a very wide range of potential well values (Vb). This means that the relationship 5716

DOI: 10.1021/la504966y Langmuir 2015, 31, 5712−5721

Article

Langmuir ⎞ ⎛ Rg d(ln k b) = −⎜⎜α′ ln + β′ ln σ + m′⎟⎟ d(Vb) Rp ⎠ ⎝

Parameters A and B are geometry- and hydrodynamic-related parameters which represent the influence of external shear and the morphology of clusters on the breakage rate, while potential well Vb is a physical−chemical parameter which represents the influence of the solvent and the material properties of particles. All of the parameters in eq 19 can be found either in this paper, in ref 32, or in handbooks. This equation can be used directly in simulations or experiments to predict the cluster breakage rate under a broad variety of external conditions. It should be noted that we have selected only one of the many expressions of the DLVO potential proposed in the literature. In fact, there are many different functional forms of the repulsive electrostatic potential and also of the attractive van der Waals interaction potential, including retardation effects.43 In the Supporting Information we have provided convincing proof that the results just shown are independent of the exact form of the potential, but only with respect to the value of the potential well. In particular, we have demonstrated that the breakage rate is almost insensitive to the long-range part of the interaction potential and depends only on the short part, i.e., the one that determines the potential well. This allows us to claim that the results obtained here are of very general validity. 3.3. Influence of the Primary Particle Size. In this section the influence of the particle size on the derivative of the breakage rate d(ln kb)/d(Vb) is investigated. The influence of the particles diameter on the breakage rate is also complex. The main reasons are that, for the same cluster, changing Rp also changes the potential energy between primary particles because of the change in the surface area. Furthermore, the influence of external shear becomes stronger as the overall particle surfaces increase. But the exponential form of the kernel of the system will be maintained. Harada et al.2 have introduced parameters named cluster cohesive strength Ah/Rp3 and fragmentation number (Fa = σRp3/Ah) to quantitatively assess the influence of the particles radius in clusters exposed to fluid flow. Harshe et al.32 have also reported a dependence of the breakage rate on the particle size scaling as (Rp0/Rp)3, with Rp0 = 1 μm being the reference size of the particles. We have decided to address the dependence of the cluster breakage rate on the primary particle size. The quantity d(ln kb)/d(Vb) is plotted as a function of particle radius under constant (Rp0/Rp)3, df, and σ in Figure 4. From Figure 4 we can see that when (Rp0/Rp)3 increases, d(ln kb)/d(Vb) increases linearly. d(ln kb)/d(Vb) should therefore be rescaled by the

(15)

where α′, β′, and m′ are fitting parameters which can be obtained from the regression of the results in Figures S3 and S4. The fitted values of α′, β′, and m′ for obtained values of d(ln kb)/d(Vb) are summarized in Table 2. Table 2. Regression Results for Breakage Parameters df

α′

m′

1.8 2.5 3.0

1.33 × 3.37 × 3.99 ×

1.8 2.5 3.0

1.17 × 1.63 × 1.72 ×

Simple Shear Flow 10−2 −7.88 × 10−4 −2 10 −4.42 × 10−3 10−2 −1.11 × 10−2 Uniaxial Extensional Flow 10−2 −8.44 × 10−4 10−2 −2.04 × 10−3 10−2 −2.51 × 10−3

β′ −2.32 × 10−3 −5.54 × 10−3 −5.03 × 10−3 −2.59 × 10−3 −3.41 × 10−3 −3.51 × 10−3

To calculate the breakage rate, eq 15 can be rewritten in exponential form as eq 16. This exponential relationship relating the breakage rate to the external shear, the morphology of the clusters, and the potential well can be obtained from the integration and is shown below: k b = Ae

−(α′ ln

Rg + β′ ln σ + m′)Vb Rp

(16)

Here A is a pre-exponential constant which was comes from integration. It can be obtained by taking any of the simulated results as an initial condition. In this work we will not discuss how to obtain A but just give our previous work as a reference,32 where detailed investigations of the breakage rate of clusters under only attractive interactions were reported. An empirical expression relating Rg/Rp and σ to kb has been obtained as follows:

⎛ R g ⎞α k b = m⎜⎜ ⎟⎟ σ β ⎝ Rp ⎠

(17)

Here α, β, and m are constants, and their values can be found in ref 32. These values were obtained without the electrostatic interaction and can be taken as initial conditions (when ψ0 = 0) for the integration of eq 15. By applying this initial condition, the expression of A can be written as follows:

A=

⎛ R g ⎞α m⎜ R ⎟ σ β ⎝ p⎠ ′

Rg

′

′

e−(α ln R p + β ln σ + m )Vb,0

(18)

Here Vb,0 is the potential well when ψ0 = 0 and can be obtained from the DLVO potential theory. A final general expression that could describe the relationship between the cluster breakage rates kb and the potential well Vb is the following: k b = Ae−BVb

Here B is given by Rg B = α′ ln + β′ ln σ + m′ Rp

(19) Figure 4. d(ln kb)/d(Vb) as a function of (Rp0/Rp)3, Rp0 = 1 μm, Rp = 0.5−2.0 μm, Nsphere = 60, df = 1.8, σs = 60 Pa, ψ0 = 0.03 V, and C0 = 0.01−10 mol/m3.

(20) 5717

DOI: 10.1021/la504966y Langmuir 2015, 31, 5712−5721

Article

Langmuir

Figure 5. Normalized fragment mass distributions in simple shear: (a) Nsphere = 100, df = 1.8, σs = 80 Pa, ψ0 = 0.03 V, and C0 = 0.01−10 mol/m3. (b) Nsphere = 100, df = 2.7, σs = 60 Pa, ψ0 = 0.03 V, and C0 = 0.01−10 mol/m3. Symbol indicate data obtained from Stokesian dynamics, and dashed lines represent the fitting of the data using the Schultz−Zimm distribution.

inverse of the cube of the particle size (Rp0/Rp)3, while the other parameters should remain unchanged. This conclusion is in accordance with the results in refs 2 and 32. Therefore, constant B in eq 19 can be rewritten as eq 21. ⎛ R p0 ⎞3⎛ ⎞ R ⎟⎟ ⎜⎜α′ ln g + β′ ln σ + m′⎟⎟ B = ⎜⎜ Rp ⎝ Rp ⎠ ⎝ ⎠

better design of industrial processes. In this section we will investigate further the breakage constant dependence on the type of flow using an approach based on the energy dissipation rate. For a generic incompressible Newtonian fluid flow field, we consider its Jacobian matrix, i.e., the gradient of the velocity field vector

(21)

⎡ ∂1u1 ∂2u1 ∂3u1 ⎤ ⎢ ⎥ ∇u ̅ = ⎢ ∂1u 2 ∂2u 2 ∂3u 2 ⎥ ⎢ ⎥ ⎣ ∂1u3 ∂2u3 ∂3u3 ⎦

3.4. Fragment Mass Distribution. The fragment mass distribution generated by clusters exposed to shear forces under only attractive interactions has been investigated in our previous work32 and successfully fitted using a Schultz−Zimm distribution. The parameters in this distribution have been obtained by regression as a function of Rg/Rp, df, and σ in Harshe’s paper.32 In this work, the same distribution model was selected, and the influence of the potential well on the Schultz− Zimm distribution was investigated. In Figure 5 the relationship between the average mass fraction of the smaller fragment (Nsfr) and the number fraction (nfr) of the smaller fragments is plotted. From Figure 5 we can see that when Rg/Rp, df, and σ were kept constant and the potential well was changed, the fragment mass distribution remained the same as in the case of purely attractive interactions and can therefore still be described by the same Schultz−Zimm distribution. The physical reason for this behavior is that breakage will occur at a location corresponding to the weakest link inside a cluster. Different cluster morphologies will have different weakest links. External shear can change the morphologies of the clusters and further change the weakest links distribution, but a change in the potential well between particles will not have any obvious influence on the weakest link distribution inside a cluster. Thus, the curves shown in Figure 5 are statistically the same, and the same fitting parameters reported in ref 32 can also be used for the evaluation of the fragment mass distribution in the presence of repulsive interactions. 3.5. Relation between Simple Shear and Uniaxial Extensional Flow. From the regression results in Table 2 we can observe that the fitting parameters under both simple shear and uniaxial extensional flow show similar trends, but the actual values are significantly different. In fact, at a given shear rate, uniaxial extensional flow is much more effective at breaking clusters. Because simple shear and uniaxial extensional flow are the two most commonly investigated flow patterns, understanding their effect on the breakage of the clusters will be important both from a fundamental perspective and in the

(22)

where ui is the component of velocity vector u̅ parallel to axis i and ∂j denotes the partial derivative with respect to space coordinate j. This matrix can be decomposed into a symmetric part, which describes the strain, and an antisymmetric part, which describes rotation. The strain can deform the cluster and finally induce breakage. From an energy point of view, the strain can increase the potential energy of a cluster, causing it to break. At the same time, the rotational part makes a cluster rotate but does not generate strain between particles inside a cluster. From the energy point of view, it can provide angular kinetic energy, but it cannot increase the potential energy of a cluster, which means that it does not have a direct influence on the breakage of a cluster. For the simple shear flow used in this paper, the symmetric and antisymmetric parts of the Jacobian matrix can be written as eq 23, ⎡0 ⎢ ∇us̅ = γṡ ⎢ 0 ⎣0 ⎡0 γṡ ⎢ = ⎢1 2 ⎣0

1 0 0 1 0 0

= Es + ωs

0⎤ ⎥ 0⎥ 0⎦ ⎡ 0 1 0⎤ 0⎤ γṡ ⎢ ⎥ ⎥ 0 ⎥ + ⎢− 1 0 0 ⎥ 2 ⎣ 0 0 0⎦ 0⎦ (23)

where Es and ωs are the strain rate tensor and the rotational tensor, respectively. For the uniaxial extensional flow, the Jacobian matrix can be written as eq 24. This is a symmetric matrix, and it does not have an antisymmetric part, implying that it equals its shear rate tensor. 5718

DOI: 10.1021/la504966y Langmuir 2015, 31, 5712−5721

Article

Langmuir ⎡2 0 0 ⎤ ⎢ ⎥ ∇ue̅ = γė ⎢ 0 −1 0 ⎥ = Ee ⎣ 0 0 −1⎦

(24)

As discussed in section 3.1, the breakage of a cluster is dominated by the potential well between two particles calculated through DLVO theory. A cluster breaks when the external shear provides enough energy to let the particles overcome the potential well. And this energy has to be provided by the strain part of the external shear. Therefore, to properly compare the effects of simple shear and uniaxial extensional flow on the breakage of a cluster, it is necessary to compare the energy that these two flow patterns can provide to a cluster. The parameter that we decided to look at in order to evaluate how the external shear gives energy to a fluid system is the viscous energy dissipation rate, which expresses the rate at which the external flow provides energy to a fluid system. According to the definition of turbulence found in textbooks,44 the time-averaged bulk energy dissipation rate per unit mass in a fluid system is defined as follows: η ε = 2 f ⟨||∇E ||22 ⟩ (25) V

Figure 6. Ratio of the breakage rate constants under simple shear and uniaxial extensional flow as a function of the applied stress ratio. One can observe that the point where the breakage rate ratio approaches 1 corresponds to values of the stress ratio of about 4 (±10%), compared to the theoretical value of 3.46. See the text for a detailed discussion.

the cluster deformation direction will always deviate from the shear strain direction. The deviation between the cluster deformation direction and shear strain direction will lower the efficiency of the energy dissipation rate in the cluster. Additionally, there are deviations in the simulation results due to the statistical variations in cluster morphology because different clusters all have different weak points. These differences cause some deviations in the simulated results from the energy dissipation criterion. 3.6. Critical Rotation Numbers. Our previous study32 has indicated that there is a critical hydrodynamic stress for each aggregate mass, below which clusters exposed to simple shear conditions do not break but undergo restructuring. The relationship between the critical stress and the critical aggregate mass was empirically modeled as follows:

Here ε is the energy dissipation rate, V is the unit volume, vf is the kinematic viscosity of the fluid, which is defined as νf = ηf /ρf

(26)

and notation ∥···∥2 indicates the L2 norm of a matrix in Hilbert space, defined as 3

|| A ||2 =

3

∑ ∑ (aij)2 i=1 j=1

(27)

where aij are the elements of matrix A. A high-energy dissipation rate means that clusters can quickly receive enough energy to overcome the potential well, and this leads to fast breakage. According to this idea, different flow patterns with the same energy dissipation rate should lead to the same breakage rate. In agreement with eq 25, the same energy dissipation rate ε for cluster breakage implies that the norm for the strain rate tensor has to be the same. This gives us a criterion to evaluate the ratio between the shear stresses (or the shear rates) of simple shear and uniaxial extensional flow necessary to have the same energy dissipation rate and therefore the same breakage rate. From the equality of the norm, σs/σe = (12)1/2. This means that in order to obtain the same rate of breakage in the two flow patterns, the shear rate in simple shear should be about 3.46 times higher than in uniaxial extensional flow. We performed simulations to verify this criterion. Figure 6 shows the simulated results from Stokesian dynamics for σs/σe versus ks/ke. From Figure 6 we observe that Stokesian dynamics simulations indicate that the simple criterion provided by the energy dissipation rate performs remarkably well and that the true ratio between the shear stresses σs/σe required to have identical breakage rates (ks/ke = 1) is only around 10−20% higher (with some fluctuations) than the theoretical estimate of (12)1/2. This deviation can be explained in two ways. The first possibility is that in uniaxial extensional flow the strain direction and the cluster deformation direction are always the same, so that clusters will more easily break along the direction of strain. Instead, for the simple shear flow, because of the influence of the rotation,

Nc = ψ (σc)−p

(28)

Here Nc is the critical aggregate mass, σc is the critical hydrodynamic stress, and ψ and p are fitting parameters. Harshe et al.32 have plotted the relationship between the stable critical mass and the hydrodynamic stress and found a power law relation between them. The reason for the existence of a critical hydrodynamic stress is that simple shear consists of both deformation and rotation. The deformation part can break the clusters, but at the same time the rotation part can increase the fractal dimension of the clusters due to restructuring.25 For the clusters with the same particle number, denser clusters are more difficult to break. When the external shear is decreased, the breakage rate also decreases so that clusters are exposed to simple shear for a longer time and can rotate for a longer time. If the rotation time is long enough, the fractal dimension of the clusters can increase. In fact, as shown in the video in the Supporting Information, when exposed to simple shear, clusters are first stretched, then rotated, and then compressed. This means that if the external shear falls below a critical value, then clusters keep on rotating and become progressively denser. Under this condition, they are very unlikely to break because of the restructuring process. The argument just presented can be rendered more quantitative by looking at the characteristic time for rotation. For simple shear flow in the x−y plane, the rotation rate can be calculated from the following equation: 5719

DOI: 10.1021/la504966y Langmuir 2015, 31, 5712−5721

Article

Langmuir

1 Ω∞ = − γṡe3 (29) 2 Here e3 is the unit vector in the z direction. From eq 29 we know that the rotation rate for simple shear is −γ̇s/2 in the z direction, so the number of rotations that a cluster undergoes before breaking can be calculated by the following equation, γs,ċ Ncrot = 4πk b,c (30)

shear, using a library of clusters generated from Monte Carlo simulations. As in our previous work, the breakage rate was defined as the time required for a cluster to undergo its first breakage event. The simulations are based on Stokesian dynamics, which can accurately account for hydrodynamic interactions of the particles mediated by the fluid. The interactions among the particles forming the clusters were modeled according to DLVO theory, while contact forces were simulated using a discrete element approach. By performing a large number of simulations, the influence of the repulsive electrostatic potential well on the breakage rate of the clusters under external shear forces was investigated in detail. It was found that the breakage rate exponentially decreases as the depth of the potential well, which two particles need to overcome in order to break a bond, increases. The first major finding of this work is that an Arrhenius-type exponential relationship between the potential well depth and breakage rate was found to satisfactorily reproduce the simulation results. The influence of all other geometrical parameters of the clusters, such as the radius of gyration, the fractal dimension, the size of the primary particles, and the external shear rate, was also investigated, and a very general correlation has been proposed for the breakage rate, with a broad range of validity. It has been further verified that the results are insensitive to the long-range part of the interparticle interaction potential, whereas they solely depend on the short-range part, which determines the potential well. This also implies that the results are insensitive to the exact expression used to model the interaction potential, as long as the depth of the potential well remains the same. The fragment mass distribution produced from the breakage event was found to be independent of the potential well and the extent of the repulsive interactions and could be described by a Schultz−Zimm distribution. A second result obtained was the quantification of the difference between simple shear and uniaxial extensional flow on the breakage rate of clusters, which has been thoroughly analyzed using a criterion based on the energy dissipation rate. It was found that in order to have the same breakage rate under the two different flow fields, the shear rate applied under simple shear has to be about 4.1 times higher than under uniaxial extensional flow, as indicated by a simple theoretical calculation and confirmed by our Stokesian dynamics simulations. Finally, a critical rotation number was introduced, indicating the maximum number of rotations that a cluster with a certain fractal dimension can undergo before breaking. When the rotation number for a cluster is higher than the critical value, it is very unlikely that the cluster will break because it will undergo restructuring as a consequence of a sufficient number of rotations. The critical rotation number is found to have an exponential dependence on the fractal dimension. The findings of this work will contribute substantially to a better understanding of the intricate relationship among interparticle interactions, shear forces, and morphology in the breakage process of fractal clusters.

Here Nrot c is the critical rotation number, γ̇s,c is the critical shear rate, and kb,c is the critical breakage rate. The meaning of this equation is the following. If clusters with a given fractal dimension make under external shear a number of rotations higher than this critical rotation number, they would keep on rotating, and it would be very unlikely that they would undergo a breakage event at this shear strength. We hereby suggest using the critical rotation numbers instead of the critical breakage stress to describe the critical point for clusters because the critical rotation number is independent of both the cluster size and the shear strength. Figure 7 shows the critical rotation

Figure 7. Critical rotation numbers of clusters computed using eq 28 under simple shear conditions for different cluster masses and potential well values as a function of the cluster fractal dimension.

numbers for different clusters with a fractal dimension of 1.8− 2.7. From Figure 7 we can find that the critical rotation numbers are independent of the cluster size, the external shear strength, and also the potential well. And the relation between the critical rotation number and the fractal dimension can be fitted by the following equation, Ncrot = 0.0456e 2.08d f

(31)

With eq 31 we can calculate the critical rotation numbers for any cluster fractal dimension. And if we want to calculate the critical hydrodynamic stress or the critical breakage rate, we only need to simultaneous solve eqs 19, 30, and 31. Namely, the breakage rate of a cluster is computed by means of eq 19, while the critical rotation number is computed via eq 31. Once the critical rotation number is obtained, the critical hydrodynamic stress can be computed from its definition, which is given by eq 30. With these three equations, we can recover the critical hydrodynamic stress and the critical breakage rate for clusters of any size.

■

ASSOCIATED CONTENT

S Supporting Information *

A video showing a cluster deformation under simple shear conditions, exemplifying the progressive densification process. Four figures discussed in the main text showing fitted simulations data are reported. Some characteristic clusters used in simulations, depicted for clarity. Detailed discussion of the influence of the functional form of the interaction potential

4. CONCLUSIONS In the present work, we performed detailed simulations to determine the breakage rate of fractal clusters under external 5720

DOI: 10.1021/la504966y Langmuir 2015, 31, 5712−5721

Article

Langmuir

(19) Dickinson, E. Structure and rheology of colloidal particle gels: Insight from computer simulation. Adv. Colloid Interface Sci. 2013, 199, 114−127. (20) Keaveny, E. E. Fluctuating force-coupling method for simulations of colloidal suspensions. J. Comput. Phys. 2014, 269, 61−79. (21) Zinchenko, A. Z.; Davis, R. H. Growth of multiparticle aggregates in sedimenting suspensions. J. Fluid Mech. 2014, 742, 577− 617. (22) Singh, A. Dynamics of suspensions of spherical doublets in simple shear and pressure driven flow. Chem. Eng. Sci. 2013, 104, 17− 24. (23) Fellay, L. S.; Twist, C.; Vanni, M. Motion of rigid aggregates under different flow conditions. Acta Mech 2013, 224, 2225−2248. (24) Yim, H. J.; Kim, J. H.; Shah, S. P. Cement particle flocculation and breakage monitoring under Couette flow. Cem. Concr. Res. 2013, 53, 36−43. (25) Ehrl, L.; Soos, M.; Lattuada, M. Generation and Geometrical Analysis of Dense Clusters with Variable Fractal Dimension. J. Phys. Chem. B 2009, 113, 10587−10599. (26) Kusters, K. A.; Wijers, J. G.; Thoenes, D. Aggregation kinetics of small particles in agitated vessels. Chem. Eng. Sci. 1997, 52, 107−121. (27) Babler, M. U. A collision efficiency model for flow-induced coagulation of fractal aggregates. AIChE J. 2008, 54, 1748−1760. (28) Babler, M. U.; Moussa, A. S.; Soos, M.; Morbidelli, M. Structure and Kinetics of Shear Aggregation in Turbulent Flows. I. Early Stage of Aggregation. Langmuir 2010, 26, 13142−13152. (29) Higashitani, K.; Iimura, K. Two-dimensional simulation of the breakup process of aggregates in shear and elongational flows. J. Colloid Interface Sci. 1998, 204, 320−327. (30) Inamuro, T.; Ii, T. Lattice Boltzmann simulation of the dispersion of aggregated particles under shear flows. Math Comput. Simul. 2006, 72, 141−146. (31) Nishiyama, T.; Yasuda, S.; Inamuro, T. Lattice Boltzmann simulation of the dispersion of aggregated Brownian particles under shear flows. Eur. Phys. J.: Spec. Top. 2009, 171, 145−149. (32) Harshe, Y. M.; Lattuada, M. Breakage Rate of Colloidal Aggregates in Shear Flow through Stokesian Dynamics. Langmuir 2012, 28, 283−292. (33) Conchuir, B. O.; Zaccone, A. Mechanism of flow-induced biomolecular and colloidal aggregate breakup. Phys. Rev. E 2013, 87, 032310. (34) Batchelor, G. K. Mass-Transfer from Small Particles Suspended in Turbulent Fluid. J. Fluid Mech. 1980, 98, 609−623. (35) Brady, J. F.; Bossis, G. Stokesian Dynamics. Annu. Rev. Fluid Mech. 1988, 20, 111−157. (36) Durlofsky, L.; Brady, J. F.; Bossis, G. Dynamic Simulation of Hydrodynamically Interacting Particles. J. Fluid Mech. 1987, 180, 21− 49. (37) Feke, D. L.; Prabhu, N. D.; Mann, J. A.; Mann, J. A. A Formulation of the Short-Range Repulsion between Spherical Colloidal Particles. J. Phys. Chem. 1984, 88, 5735−5739. (38) Pantina, J. P.; Furst, E. M. Elasticity and critical bending moment of model colloidal aggregates. Phys. Rev. Lett. 2005, 94, 138301. (39) Becker, V.; Briesen, H. Tangential-force model for interactions between bonded colloidal particles. Phys. Rev. E 2008, 78, 061404. (40) Rogak, S. N.; Flagan, R. C. Stokes Drag on Self-Similar Clusters of Spheres. J. Colloid Interface Sci. 1990, 134, 206−218. (41) Ekiel-Jezewska, M. L.; Gubiec, T.; Szymczak, P. Stokesian dynamics of close particles. Phys. Fluids 2008, 20, 063102. (42) Markutsya, S.; Fox, R. O.; Subramaniam, S. Characterization of sheared colloidal aggregation using Langevin dynamics simulation. Phys. Rev. E 2014, 89, 062312. (43) Berg, J. C. An Introduction to Interfaces and Colloids: The Bridge to Nanoscience; World Scientific Publishing: Singapore, 2010; pp 525− 615. (44) Batchelor, G. K. An Introduction to Fluid Dynamics; Cambride University Press: Cambridge, U.K., 2000; pp 131−173.

on the cluster breakage rate. The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/la504966y.

■

AUTHOR INFORMATION

Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Financial support by the Swiss National Science Foundation (grant number PP00P2133597/1) and by the Adolphe Merkle Institute is gratefully acknowledged.

■

REFERENCES

(1) Garland, S.; Gauthier, G.; Martin, J.; Morris, J. F. Normal stress measurements in sheared non-Brownian suspensions. J. Rheol. 2013, 57, 71−88. (2) Harada, S.; Tanaka, R.; Nogami, H.; Sawada, M.; Asakura, K. Structural change in non-fractal particle clusters under fluid stress. Colloids Surf., A 2007, 302, 396−402. (3) Richardi, J. One-dimensional assemblies of charged nanoparticles in water: A simulation study. J. Chem. Phys. 2009, 130, 044701. (4) Moussa, A. S.; Lattuada, M.; Conchuir, B. O.; Zaccone, A.; Morbidelli, M.; Soos, M. Flow-Induced Aggregation and Breakup of Particle Clusters Controlled by Surface Nanoroughness. Langmuir 2013, 29, 14386−14395. (5) Zaccone, A.; Wu, H.; Gentili, D.; Morbidelli, M. Theory of activated-rate processes under shear with application to shear-induced aggregation of colloids. Phys. Rev. E 2009, 80, 051404. (6) Zaccone, A.; Gentili, D.; Wu, H.; Morbidelli, M. Shear-induced reaction-limited aggregation kinetics of Brownian particles at arbitrary concentrations. J. Chem. Phys. 2010, 132, 134903. (7) Fellay, L. S.; Vanni, M. The effect of flow configuration on hydrodynamic stresses and dispersion of low density rigid aggregates. J. Colloid Interface Sci. 2012, 388, 47−55. (8) Vanni, M.; Gastaldi, A. Hydrodynamic Forces and Critical Stresses in Low-Density Aggregates under Shear Flow. Langmuir 2011, 27, 12822−12833. (9) Zaccone, A.; Soos, M.; Lattuada, M.; Wu, H.; Babler, M. U.; Morbidelli, M. Breakup of dense colloidal aggregates under hydrodynamic stresses. Phys. Rev. E 2009, 79, 061401. (10) Soos, M.; Wang, L.; Fox, R. O.; Sefcik, J.; Morbidelli, M. Population balance modeling of aggregation and breakage in turbulent Taylor-Couette flow. J. Colloid Interface Sci. 2007, 307, 433−446. (11) Arosio, P.; Rima, S.; Lattuada, M.; Morbidelli, M. Population Balance Modeling of Antibodies Aggregation Kinetics. J. Phys. Chem. B 2012, 116, 7066−7075. (12) Prakash, A. V.; Chaudhury, A.; Barrasso, D.; Ramachandran, R. Simulation of population balance model-based particulate processes via parallel and distributed computing. Chem. Eng. Res. Des. 2013, 91, 1259−1271. (13) Schlauch, E.; Ernst, M.; Seto, R.; Briesen, H.; Sommerfeld, M.; Behr, M. Comparison of three simulation methods for colloidal aggregates in Stokes flow: Finite elements, lattice Boltzmann and Stokesian dynamics. Comput. Fluids 2013, 86, 199−209. (14) Wilson, H. J. Stokes flow past three spheres. J. Comput. Phys. 2013, 245, 302−316. (15) Soos, M.; Sefcik, J.; Morbidelli, M. Master curves for aggregation and gelation: Effects of cluster structure and polydispersity. Ind. Eng. Chem. Res. 2007, 46, 1709−1720. (16) Dargazany, R.; Itskov, M. Yield behavior of colloidal aggregates due to combined tensile-bending loads. Phys. Rev. E 2012, 85, 051406. (17) Conchuir, B. O.; Harshe, Y. M.; Lattuada, M.; Zaccone, A. Analytical Model of Fractal Aggregate Stability and Restructuring in Shear Flows. Ind. Eng. Chem. Res. 2014, 53, 9109−9119. (18) Soos, M.; Ehrl, L.; Babler, M. U.; Morbidelli, M. Aggregate Breakup in a Contracting Nozzle. Langmuir 2010, 26, 10−18. 5721

DOI: 10.1021/la504966y Langmuir 2015, 31, 5712−5721