PHYSICAL REVIEW E 86, 051403 (2012)

Effect of hydrodynamic interactions on rapid Brownian coagulation of colloidal dispersions
Yuki Matsuoka,1,* Tomonori Fukasawa,2 Ko Higashitani,2 and Ryoichi Yamamoto2,†
1

Production Engineering Research Laboratory, Sumitomo Bakelite Co., Ltd., Shizuoka 426-0041, Japan2
Department of Chemical Engineering, Kyoto University, Kyoto 615-8510, Japan(Received 19 July 2012; published 14 November 2012)

The rate of rapid Brownian coagulation is investigated for dispersions of spherical particles with particle volume
fractions ranging from φp = 0.003 to 0.1 by the direct numerical simulation method. This method explicitly
considers hydrodynamic interactions (HIs) between particles by simultaneously solving for the motions of thedispersed particles and the host fluid. In the dilute limit, the rate of rapid Brownian coagulation decreases toapproximately 0.3–0.5 times the theoretical Smoluchowski rate. We compare this result with results of previously
reported experiments and theoretical predictions and find a strong correlation between them. This demonstrates
that HIs between particles significantly reduce the coagulation rate. Moreover, the volume fraction dependenceof the coagulation rate indicates that the coagulation rate increases with increasing volume fraction. At highparticle volume fractions, the initial coagulation stage is affected by heterogeneous coagulation process before
the steady state is reached.
DOI: 10.1103/PhysRevE.86.051403

PACS number(s): 82.70.−y, 82.20.Wt, 83.10.Pp, 61.43.Hv

I. INTRODUCTION

Since Smoluchowski developed his theory of the Brownian
coagulation process in colloidal dispersions based on population balance equations of coagulates and coagulation kinetics[1], many theoretical, experimental, and simulation studies
have investigated the validity of this theory. A major deficiencyof the Smoluchowski theory is that it predicts a coagulation ratethat is approximately twice as large as the actual coagulation
rate of colloidal dispersions [2–6]. This discrepancy has been
mainly attributed to the theory’s failure to consider interparticle interactions. Thus, a correction factor for the coagulation
rate that considers interparticle interactions (hydrodynamic
lubrication, Derjaguin-Landau-Verwey-Overbeek interaction[7,8], etc.) has been proposed [9–11]. Many experimental
examinations of this correction have been conducted [2–6,12].
It has been suggested that the lower observed aggregation rateis due mostly to hydrodynamic lubrication forces between
particles that hinder two particles approaching each other.
However, these experimental examinations depend on the
proposed correction factor itself. Consequently, evaluation of
the hydrodynamic term in the correction factor is replaced byevaluation of other interparticle interactions, such as Van der
Waals attractions. Thus, additional direct evaluations of thesehydrodynamic effects would be valuable.
In the Smoluchowski formalism, particles are assumed to be
spherical and the coagulation process is limited to simple twobody collisions between coagulates of the same size. Theseassumptions are only valid in the early stages of coagulation in
dilute dispersions. In addition to the coagulation process itself,
the formation of fractal structures in the coagulate clusters
significantly affects the coagulation rate in the later stages ofcoagulation. An increase in the coagulation rate in the laterstages of coagulation is observed experimentally. A theoryof coagulation kinetics that is applicable to the later stages*†

ymatsuoka@sumibe.co.jp
ryoichi@cheme.kyoto-u.ac.jp1539-3755/2012/86(5)/051403(9)of coagulation has been proposed. This model considers thesteric effect of coagulation on the coagulation rate [13,14].
The coagulation rate is known to be enhanced in the very
early stages of coagulation in which the spatial distributionof particles has not yet reached a steady state [1,15]. As
the volume fraction of dispersed particles increases, thecharacteristic coagulation time (i.e., the time for the initial
particle number to be halved) becomes shorter. At a certain
volume fraction, the coagulation time becomes comparable
to the transition period in which the effects of the unsteadystate are significant. Therefore, coagulation times at higher
volume fractions are greatly affected by the effects of theunsteady state. An increase in the coagulation rate at highparticle volume fractions has been experimentally observed[3]. However, a detailed mechanism for this enhancement in
the coagulation rate has yet to be elucidated.
Computer simulations are potentially superior to physicalexperiments for investigating these problems since simulations avoid difficulties such as preparing an initial dispersed
particle distribution and observing aggregate forms at highconcentrations. They can thus provide essential insights thatcan be used to evaluate proposed theories. Several simulationsof colloidal coagulation have been reported, including somethat employ Brownian dynamics (BD) and Langevin dynamics(LD) [14,16–19]. In these simulations, the coagulation rate
has been observed to increase with increasing particle volume
fraction [14,16–18]. Lattuada recently proposed a correction to
the Smoluchowski equation for high volume fractions by usinga result derived from trapping theory for the diffusion-limited
reaction of chemical species. His theory and results appear to
be consistent when assessed by a BD simulation [19].
The BD and LD simulation methods are effective forcalculations of large systems because they have low computational costs. However, their handling of hydrodynamic
interactions (HIs), such as lubrication forces, that act between
particles is inadequate and may give rise to significant
discrepancies between calculation and experimental results.
This is especially true for colloidal particles dispersed inviscous liquids, such as water, solvents, or resins. In particular,
simulations performed without considering HIs (such as the051403-1©2012 American Physical Society

MATSUOKA, FUKASAWA, HIGASHITANI, AND YAMAMOTO

PHYSICAL REVIEW E 86, 051403 (2012)

BD method) yield coagulation rates in the dilute limit that arevery similar to the theoretical value given by the Smoluchowski
theory, which does not consider HIs. Although BD simulations
with average hydrodynamic corrections have recently been
conducted, a reduction in the coagulation rate such as thatfound in experiments at low volume fractions was notobserved [18].
Brady and Bossis developed the Stokesian dynamics (SD)simulation method that correctly incorporates hydrodynamic
effects in the low Re limit [20]. The SD method represents
short-range lubrication forces and long-range hydrodynamic
forces by the sum of two-body analytic interactions. Thismethod has been adapted to the rheology of colloidal
dispersions and sedimentation phenomena in which fluidhydrodynamics significantly affect particle dynamics [20,21].
Although the SD method has provided much valuable information regarding these phenomena, the assumption that HIs canbe treated as the sum of two-body interactions fails to consider
many-body HIs, which may be important in coagulationof high volume fractions [22,23]. Consequently, to gain acomplete understanding of these coagulation processes, it isdesirable to consider many-body HIs.
In recent years, several numerical methods have beendeveloped that accurately simulate dispersion in a variety ofsituations, including those described above [22,24–27]. These
dispersion modeling techniques are all based on the same approach that involves determining fluid motion simultaneouslywith particle motion. We refer to this approach as the direct
numerical simulation (DNS) approach. This approach enables
us to accurately consider full time-dependent many-body HIs.
In this study, we apply a direct numerical scheme based on
the smoothed profile (SP) method [24] to a monodisperse
dispersion of neutrally buoyant attractive Brownian particles ina stationary fluid. In SP method, the Navier-Stokes equation for
the fluid motion is discretized on a regular grid and Newton’s
equations for particle motion are solved simultaneously withfluid motion.
In the present study, we examine the rate of rapid Brownian
coagulation for a low particle volume fraction φp = 0.003 andfor a relatively high volume fraction φp = 0.1 using the DNS
method, which accurately accounts for HIs. Specifically, we
focus on whether our simulation gives a lower coagulation ratefrom that predicted by the Smoluchowski theory. Because oursimulation explicitly considers full many-body HIs, successful
reproduction of the reduced coagulation rate would imply thatthe reduction is due to hydrodynamic effects. In addition,
the volume fraction dependence of the coagulation rate isinvestigated. In these analyses, the focus is mainly on earlystage of coagulations. However, at high particle volume
fractions, our results show deviations from the early time
behavior. Some interpretations to this effect will be given as
well.

where nk is the concentration of particles with size k attime t and β(i,j ) is the collision frequency factor between
particles with sizes i and j . Since the collision frequency factor
depends on the individual properties of a particular motion
and on particle interactions, its functional form is generally
unknown except for some ideal cases. Once the collision
frequency factors and initial conditions have been determined,
it is possible to obtain information about the coagulation rates
and cluster size distributions during coagulation by solving
the simultaneous equations Eq. (1) for each cluster size.
Smoluchowski considered particle collisions in terms of therelative diffusion of one particle with respect to another andshowed that a simple form of the collision frequency factor
in the early stages of coagulation of initially monodisperse
suspensions could be determined by solving the diffusion
equation in spherical coordinates [1]. The particle spatial
distribution of i-fold clusters around j -fold clusters ni (r,t)
and the collision frequency factor β(i,j ) are expressed by √r−Rij
4Dij t
2Rij
Rij
−ξ 2
+√
ni (r,t) = ni (t) 1 −e dξ , (2)
r
πr 0


∂ni
β(i,j )ni (t) = 4π Rij2 Dij(3)∂r r=RijRij= 4π Rij Dij ni (t) 1 + ,(4)
π Dij tII. KINETIC THEORY OF RAPID COAGULATION

In general, the kinetic equation of Brownian coagulation is
given by
i=k−1
∞

dnk
1 
β(i,j )ni nj −β(k,i)nk ni ,=dt
2 i=1,i+j =k
i=1

(1)

where r is the distance between the two clusters and ni (t) =ni (∞,t),Rij = Rc,i + Rc,j ,Dij = Di + Dj , where Rc,i , Di
are the collision radius and the translational diffusion constantof the cluster with size i, respectively. The last term in Eq. (4)represents the unsteady effect that increases the coagulation
rate during the initial stages of coagulation. This effectreflects the unsteady state of the particle spatial distribution
from the beginning of coagulation to the establishment ofa steady-state particle distribution. The time dependence ofthe particle spatial distribution expressed by Eq. (2) is shownin Fig. 1. Equation (2) assumes a spatially uniform particle
distribution as the initial condition. As coagulation starts andproceeds, i-fold clusters diffuse toward j -fold clusters, leading
to collisions and a subsequent reduction in the density nearj -fold clusters. Accordingly, ∼Rij2 /Dij , the gradient of particle
distribution at r/Rij = 1 is often steep and the particle flux
becomes large. At later times, the area with reduced density
extends to larger interparticle distances. Ultimately, little further change in the particle distribution occurs and a steady-stateparticle distribution, ni (r,t) = ni (t)(1 − Rij /r), is established.
In our case, the initial particle distribution is not uniform,
particularly at high volume fractions. In addition, interparticle
interactions due to HIs and direct interparticle interactions
may alter the time dependence of the particle distribution.
The unsteady effects on the coagulation rate are expected tobe more complicated than was proposed by Smoluchowski.
These complex initial conditions or boundary conditions in
the diffusion equation have been studied, especially in thefield of the kinetics of chemical reactions in solution [28].
Assuming that the Stokes-Einstein relation, Di =kB T /6π μRH,i , and the condition in the early stages ofcoagulation, Rij /2 = Rc,i = Rc,j = RH,i = RH,j , both holdand a steady-state particle distribution is established (i.e.,

051403-2

EFFECT OF HYDRODYNAMIC INTERACTIONS ON RAPID . . .

coagulation process for a system of monodisperse sphericalparticles in which the particles interact with one anotherthrough short-range attractive interactions.
We consider a monodisperse dispersion of Np sphericalparticles with diameter σ = 2a in a Newtonian host fluid,
where a is the radius of a particle. The particles interact via amodified 200:100 Lennard-Jones potential as

1.2
1
0.8
ni(r,t) / ni(t)PHYSICAL REVIEW E 86, 051403 (2012)

VLJ (r) = 4[(σ/r)200 − (σ/r)100 ],0.6*
t/t =0t/t*=0.001*
t/t =0.01*
t/t =0.1*
t/t =1*0.4
0.2

t/t =10ni(t)(1-Rij /r)

0
0

1

2

3

4

5

r / Rij
FIG. 1. The time dependence of particle spatial distribution isexpressed by Eq. (2). Vertical and horizontal axes are normalizedby ni (t) and Rij , respectively. Particle distributions at t/t ∗ =0,0.001,0.01,0.1,1,10 are shown, where t ∗ = Rij2 /Dij . The steadystate particle distribution at t/t ∗  1 is also plotted.

t  Rij2 /Dij ), Eq. (4) can be simplified as β(i,j ) = 2K 0 =8kB T /3μ. Here, kB is the Boltzmann constant and T is thetemperature. Furthermore, μ is the viscosity of the fluid, RH,iis the hydrodynamic radius of a cluster with size i, and K 0 =4kB T /3μ is the so-called Smoluchowski rate constant for rapid
Brownian coagulation. Smoluchowski solved Eq. (1) for the
early stages of the coagulation of an initially monodisperse,monomer suspensions. The total particle concentration, Nt ,and the concentration of particles of each cluster size, ni , arerespectively given by1
1t−= Kt =,Nt (t) N0
N0 tc
 i−1 

t
t i+1
1+
ni (t) = N0
,
tc
tc

(5)
(6)

where N0 is the initial total particle concentration and thecoefficient K on the right-hand side of Eq. (5) represents thecoagulation rate in terms of the total number concentration.
The index i  1 in Eq. (6) represents the cluster size (numberof particles in the cluster). For the ideal conditions proposedby Smoluchowski, K = K 0 . The coagulation time tc is definedas the time for the total particle concentration to decrease tohalf the initial particle concentration, Nt (tc ) = N0 /2. The time
constant tc represents a useful time scale for identifying theearly stages of coagulation.
III. SIMULATION METHOD

We use a direct numerical scheme based on the SP
method to consider the hydrodynamic many-body interaction
between particles [24,29]. This method has been adapted torheological studies and studies of the gelation of colloidal
particle dispersions; it has been confirmed reasonably well
[23,30]. In the present study, we consider the rapid Brownian(7)

where r is the distance between two particles and  isthe interaction strength. The large power index in Eq. (7)produces a short interaction range, corresponding to the rapid
coagulation of sticky particles. In this study, the value of  isa fixed at 7.03kB T . This attractive potential strength is higher
than the thermal fluctuation strength kB T and is sufficiently
high to prevent coagulated particles from reversibly dispersing
through thermal fluctuations.
The time evolution of particle i with mass Mi = π σ 3 ρp /6
and moment of inertia Ii is governed by Newton’s equationsof motion as
CV
Mi V̇ i = F H
i + Fi + Gi ,Ṙi = V i ,I i · ˙ i = N H
i + Gi ,(8)(9)

where ρp is the density of the particle, Ri is the position of
the particle i, V i is the translational velocity of the particle,
Hand i is the rotational velocity of the particle. F H
i ,N i are thehydrodynamic forces and torques exerted by the fluid on eachparticle. F Ci is the interparticle force arising from the potentialof Eq. (7). G V
i ,G i are respectively random forces and torques
due to thermal fluctuations. These random fluctuations are
assumed to be Markovian and determine the temperature T .
The procedure for determining the temperature is describedin [31]. In the SP method, the velocity and pressure fields,v(x,t) and p(x,t), are defined according to three-dimensional
Cartesian grids that consist of fluid and particle domains.
Each domain in the grids is distinguished by a smoothedfunction φ(x,t), which is 1 in the particle domains and 0
in the fluid domains. These domains are separated by a thininterfacial domain of thickness ξ . The time evolution of thevelocity field is governed by the Navier-Stokes equation with
the incompressibility condition ∇ · v = 0 given byρf (∂t v + v · ∇v) = −∇p + μ∇ 2 v + ρf φ f p ,(10)where ρf is the density of the fluid and φ f p is the body force
that ensures the rigidity of the particles and the appropriate
nonslip boundary conditions at the fluid/particle interface. Thisphenomenon is explained in more detail elsewhere [24,29].
The unit of length in the lattice spacing is ; the otherfundamental units are μ and ρf . Similarly, the unit of timeis τ = ρf 2 /μ and the unit of energy is 0 = μ2 /ρf .
Unless otherwise stated, we set = 1, τ = 1, μ = 1, ρf =1, ρp = 1, a = 4, σ = 8, ξ = 2, and the Lennard-Jones time
unit τM = [(Mi σ 2 )/]1/2 ≈ 131. Assuming dispersions ofneutrally buoyant particles of radius 1 μm in water at room
temperature, our unit lengthand time τ correspondto 0.25 μm and 0.0625 μs. The second-order Runge-Kutta
scheme is used to integrate Newton’s equations. The NavierStokes equation is discretized by a Fourier spectral schemein space and by a second-order Runge-Kutta scheme in time.

051403-3

MATSUOKA, FUKASAWA, HIGASHITANI, AND YAMAMOTO

PHYSICAL REVIEW E 86, 051403 (2012)

The discretized time step is h = 0.074 81. Although this value
is chosen from the stability condition of the Navier-Stokes
equation, it is suitable for the particle equations because his smaller than the period of the particle’s vibration in thepotential well τV = 2π [Mi /V̈LJ (21/100 σ )]1/2 ≈ 11.7. Thus, asufficiently short time step relative to the typical time scale forparticle’s motions is adopted.

30

10 (1/Nt-1/N0)

25

15
2
0
10

-2

VLJ / kBT

-4
-6

5

-8
0.95

1

1.05

1.5

2

2.5

1.1

0

-5

-10
0

0.5

1

3

r/σ
FIG. 2. Interparticle interaction potential. The particles interact
via a modified 200:100 Lennard–Jones potential expressed by Eq. (7).
The dashed line indicates an that is energy 1.17kB T higher than thepotential minimum. The inset shows an enlargement about r/σ = 1.

20

15

-3

IV. RESULTS AND DISCUSSION

The simulations were performed in a three-dimensional
cubic box whose side length is L with periodic
boundary conditions. All simulations assumed a constant volume system L = 256, where the number ofparticles was Np = 187,281,375,625,938,1877,3754, and6258. These conditions respectively correspond to φp =0.003,0.0045,0.006,0.01,0.015,0.03,0.06, and 0.1 where
φp = π σ 3 Np /6L3 . For spherical particles, a = 4 and ξ = 2.
The thermal fluctuation forces on the particles are controlled
for the temperature kB T /0 = 5.97. The particle spatial
distribution is initially random. To increase the statisticalsignificance, five simulations were conducted for different
random particle distributions. The total number of calculation
steps is set to 25 000, which corresponds to 1874τ per
calculation.
We employed the following criterion to determine aggregation of two particles. When the center-to-center distance ofthe particles is less than the threshold distance Rth = 1.012σ ,the particle pair is included in the same cluster. As shownin Fig. 2, this threshold value corresponds to an interparticle
distance where the potential energy is only 1.17kB T higher
than the minimum of the interparticle interaction potential VLJ .
The aggregates obtained consist of primary particles and form
fractal structures at later stages of coagulation. One particle in

φp=0.003
φp=0.0045
φp=0.006
φp=0.01
φp=0.015
φp=0.03
φp=0.06
φp=0.1
Smoluchowski

10

5

0
0

400

800

1200

1600

t
FIG. 3. Relation between 1/Nt (t) − 1/N0 and t. Dots represent
the averages obtained for five different initial configurations. The
dashed line shows the Smoluchowski theory given by Eq. (5), where
K = K 0 = 4kB T /3μ.

each doublet can rotate freely around the other, allowing forrearrangements of the cluster structure.
Figure 3 shows the temporal change of the reciprocal of
the total particle concentration for different initial particle
volume fractions. The slope represents the rate of coagulation.
Smoluchowski’s theoretical values are also plotted (dottedline) in this figure. The rapid coagulation rate given in Fig. 3 isless than the theoretical value predicted by Smoluchowski’s
theory for the dilute condition φp = 0.003. Furthermore,increasing the volume fraction causes a corresponding increasein the slope of the plot in the early stages of aggregation.
At higher volume fractions, φp  0.03, an additional, drastic
increase in the slope is observed in the later stages ofaggregation.
For similar experimental estimations of the coagulation
rate at the early stages of coagulation [3], data in the timeinterval from t = 0 to tc are used for least-squares fitting of thecoagulation rate K by Eq. (5). For φp  0.015, the total particle
number was not halved in the computed times and the fitting
was performed over the whole data from the first time to thelast. For φp = 0.03,0.06,0.1, tc = 1745,661,288, respectively.
Figure 4 shows the volume fraction dependence of the ratio of
the rate in the early stages of coagulation to Smoluchowski’s
coagulation rate K/K 0 . The ratio asymptotically approaches
a certain value as it moves closer to the dilute limit volume
fraction. The coagulation rate obtained is estimated to beapproximately 0.3 to 0.5 times lower than that predicted by
Smoluchowski for the dilute limit. This agrees relatively wellwith previous experimental results (e.g., the open point in
Fig. 4) for which K/K 0 is reported to lie between 0.4 and0.6 [2–6]. Nevertheless, our values are slightly lower thanthose obtained experimentally.
The deviation of the observed coagulation rate from Smoluchowski’s theory can be mainly accounted for by the fact that
Smoluchowski’s theory does not consider interactions between

051403-4

EFFECT OF HYDRODYNAMIC INTERACTIONS ON RAPID . . .

PHYSICAL REVIEW E 86, 051403 (2012)

100

1.4
1.2

Smoluchowski

10-1

ni / N0K/K

01

φp=0.003
φp=0.0045
φp=0.006
φp=0.01
φp=0.015
φp=0.03
φp=0.06
φp=0.1
Eq. (6)This work
Higashitani et al. (1979)
0Eq. (13) (α =0.326 )

0.8
0.6

10

-2

10

-3

Cluster Size
(Fold)

1
2

3

Eq. (11)
0.4
0.2
10-4

10-3

10-2

10-1

φp
FIG. 4. Volume fraction dependence on the ratio of the coagulation rate K/K 0 . The coagulation rate K is obtained by fitting the data
between t = 0 ∼ tc with Eq. (5). The filled points represent calculated
values in this study and the open point represents an experimental
result of Higashitani et al. [4]. The arrows indicate Smoluchowski’s
result and αBr = 0.5 estimated by Eq. (11). The solid line represents
Eq. (13).

particles. This theory ignores interparticle interactions,whereas attractive forces due to long-range Van der WaalsLondon potentials and repulsive forces due to HIs are known
to greatly affect actual coagulation processes. The presence of
these interactions necessitates corrections to Smoluchowski’s
theory. The effect on the interparticle potential V (s) withrespect to the coagulation rate was introduced as a stability
factor by Fuchs, where s = r/a [9]. Spielman then added
the hydrodynamic lubrication effect between particles to thestability factor as a correction factor for the diffusion constantC(s) [10]. According to their corrections, the ratio of thetheoretical coagulation rate (including these interparticle interactions) to Smoluchowski’s rate of rapid Brownian coagulation
can be written as
  ∞
−1
C(s) V (s)/kB T
0αBr = K/K = 2
e
ds
.
(11)
s2
2
To calculate Eq. (11), Honig et al. used the Van der Waals
attractive potential generated by Hamaker [32] and the approximate expressions for HIs in Eq. (12) as V (s) and C(s) in
Eq. (11), respectively [11],C(s) = (6s 2 − 11s)/(6s 2 − 20s + 16).

(12)

These authors reported that αBr = 0.3–0.6, although this value
varies depending on the strength of Van der Waals attractions.
Their values are close to those obtained experimentally and thecorresponding Hamaker constant agrees well with theoreticalestimates [11]. Similar verifications have been conducted fora variety of particles and a variety of measurement techniques
to verify the validity of the theoretical correction, according to
Eq. (11) [4,6].
In our case, we estimated αBr using Eqs. (7), (11), and(12), where the integral range in Eq. (11) is limited from

10-3

10-2

10-1
t / t′c100101

FIG. 5. Temporal concentration changes for particles of different
sizes. Dots represent the results for up to a threefold cluster at φp =0.01,0.03,0.1. The horizontal axis is normalized by tc = 1/KN0 ,as is expected from coagulation rates K. The dashed line indicates
Eq. (6).

the distance of the Lennard-Jones potential minimum (i.e.,s = 21/100 σ/a = 2.014) to a sufficiently long distance. As
indicated by the arrow in Fig. 4, the reduction in the
coagulation rate due to interparticle interactions is estimatedto be approximately 50% of Smoluchowski’s value, which
does not include interparticle interactions. The obtained αBr
corresponds reasonably well with the present numerical resultsfor K/K 0 at lower volume fractions. The reduction in thisestimated value relative to Smoluchowski’s value is mainly
due to the hydrodynamic lubrication effect, as represented by
Eq. (12).
It should be emphasized that our simulation results withrespect to the coagulation rate reduction from Smoluchowski’s
theoretical value are obtained by calculating HIs between particles without employing any approximations. Conventional
studies that employ BD simulations often have discrepancies
between experimental results and simulation-based results andthe ratio of the coagulation rates K/K 0 is approximately unity
due to HIs not being considered in the simulations. In contrast,
our simulation results are consistent with current experimental
and theoretical results, clearly indicating that HIs between
particles play an important role in reducing the coagulation
rate.
Furthermore, as shown in Fig. 4, increases in the volume
fraction lead to gradual increases in the coagulation rate.
This result is qualitatively consistent with experimental andsimulation-based results [3,16,18,33]. The solid line shownin Fig. 4 is given by a formula based on the regression
equation for the ratio of the coagulation rates K/K 0 in thehigh concentration system proposed by Heine et al. [16] as

051403-5


K/K = α 1 +
0

0


2.5
−2.7
.
(− log φp )1 − φp(13)

MATSUOKA, FUKASAWA, HIGASHITANI, AND YAMAMOTO

10

PHYSICAL REVIEW E 86, 051403 (2012)04(a) φp=0.03
t/t′c=0.0210
t/t′c=0.210
t/t′c=0.419
t/t′c=0.629
t/t′c=0.839
t/t′c=1.05

10

-3

10

-4

φp=0.1
φp=0.06
φp=0.03
1/2
0.326[1+σ/(2πD1t) ]2.50ni / N0

-2

(a)

3

K′ / K

10-1

10

3.5

2
1.5
1
0.5

10-5
100

101

102

0
1
10

103

10

2

t

Cluster size i
10

4

0

(b) φp=0.1

-3

10

-4

3

(b)

φp=0.015
φp=0.01
φp=0.006
φp=0.0045
φp=0.003
0.326[1+σ/(2πD1t)1/2]

0

2.5
K′ / Kni / N0

-2

10

3.5

t/t′c=0.127t/t′c=1.27
t/t′c=2.53
t/t′c=3.80
t/t′c=5.07
t/t′c=6.33

10-1

10

3

10

2
1.5
1
0.5

10-5
100

101

102

0
1
10

103

10

2

3

10
t

Cluster size i
FIG. 6. Temporal variation of particle size distributions at(a) φp = 0.03, (b) φp = 0.1. Particle size distributions at t/τ =37.48,374.8,749.6,1124,1499,1874 are shown. Time is normalizedby tc . Lines indicate the lines of Eq. (6) calculated using the obtained
coagulation rates K.

Because Heine et al. obtained this formula from the calculation
results given by LD simulations in which HIs between particles
are not considered, the prefactor α 0 = 1 appears in their
original formula. Their equation reduces to Smoluchowskifs
result K/K 0 = 1 in the dilute limit (− log φp )−2.7
1. We
thus treat α 0 as the scaling factor. By least-squares fitting
α 0 to our result, we obtained α 0 = 0.326. This result agrees
well with the form of Eq. (13) despite differences in thecalculation methods for considering HIs. It thus seems thatHIs directly affect the magnitude of the coagulation rate but
they do not affect the dependence of the coagulation rate on theinitial particle concentration for these relatively low volume
fractions.

FIG. 7. Temporal variation of coagulation rate for (a) φp  0.03
and (b) φp < 0.03. Temporal coagulation rate K is obtained from thedifferential of 1/Nt at each time in Fig. 3. In (b), K are averaged
in time intervals of 100τ . The solid line indicates Smoluchowski’s
form for the coagulation
rate, including unsteady-state effects
[i.e., α 0 (1 + Rij / π Dij t)]. Here, the relationship α 0 = 0.326,Rij =σ,Dij = 2D1 = 2kB T /6π μa is assumed.

Figure 5 shows the relationship between cluster concentration and the size and time at each volume fraction. Time in
the figure is normalized by the estimated coagulation time tc =1/KN0 , where K is the coagulation rate obtained by the fitting.
This figure also shows a theoretical line given by Eq. (6) inwhich tc has been replaced with tc . As coagulation progresses,
the number of primary particles decreases and the populationsof doublet and triplet clusters increase. Next, the concentrationof doublet and triplet clusters decreases after it reaches amaximum at around t = tc . In the early stages of coagulation(i.e., t/tc < 1), the concentration of each cluster size is inreasonable agreement with Eq. (6). However, in later stages

051403-6

EFFECT OF HYDRODYNAMIC INTERACTIONS ON RAPID . . .

PHYSICAL REVIEW E 86, 051403 (2012)

104

2〈Rg〉 / a

1.6

1.4(a) t/τ ≤ 1220

φp=0.003
φp=0.0045
φp=0.006
φp=0.01
φp=0.015
φp=0.03
φp=0.06
φp=0.1

3

10

Nj

1.8

2(Rg / a)1.70
2(Rg / a)1.80

2

10

1.2
101
1

0.8
101

102

100
100

103
t

101
Rg,j / a102

104

FIG. 8. Temporal variation of number average radius of gyration.
Vertical axis is normalized by the primary particle’s radius a. The
averaging in Eq. (14) includes the contribution of primary particles
assumed to be clusters with Rg = a.

3

10

Nj

of coagulation (i.e., t/tc > 1), especially for higher volume
fractions φp  0.03, the concentration of doublet and triplet
clusters falls significantly below the theoretical line. This result
can be largely attributed to cluster-cluster or heterogeneous
coagulation, as shown below. Figure 6 shows the temporalchanges in the particle size distribution at φp = 0.03,0.1. Thelines in the figure indicate Eq. (6) calculated using the obtained
coagulation rates K. As seen in Fig. 5, data in the early stagesof coagulation match Eq. (6) for a broad range of cluster
sizes. However, in later stages of coagulation, the proportionof medium-sized clusters from twofold to several dozen folddecreases and larger clusters with sizes of several hundred fold
are formed. These behaviors are thought to correspond to thedrastic increase of the slope in the later stages of coagulationin Fig. 3. In this coagulation stage, coagulations between
clusters with nonspherical bulky shapes or between clusters
with heterogeneous collision diameters (e.g., between a large
cluster and the small particles around it) [13,16] are observedin our simulations. At the end of the simulation of φp = 0.1, asingle large cluster expands to include the entire system, likein gelation (not shown).
As seen in our simulation results, the coagulation process at
low particle volume fractions is still within the framework of
Smoluchowski’s assumptions until the end of the simulation
runs. However, at higher particle volume fractions, our results
show notable deviations from the early stage behavior. Next,
we discuss briefly on the late stage behavior.
The time at which the coagulation rate is measured affects
the coagulation rate and its volume fraction dependence.
Figure 3 indicates that the gradient of 1/Nt temporally increases in the later stages of coagulation for high volume fractions. Figure 7 shows the temporal change in the coagulation
rate, where the temporal coagulation rate K is obtained fromthe difference between 1/Nt with respect to time, at a particulartime in Fig. 3. At low volume fractions, short time changes of(b) t/τ > 1220
2(Rg / a)1.70
2(Rg / a)1.80

2

10

101

100
100

101
Rg,j / a102

FIG. 9. The relationships between the cluster size Ni and theradius of gyration of cluster Rg,j at φp = 0.1 for (a) t/τ  1220 and(b) t/τ > 1220. Equation (16) is shown by a solid line (df = 1.7)
and a dashed line (df = 1.8).

1/Nt are discrete due to the small numbers of particles, causing
numerical differential of 1/Nt to vibrate largely. K is thus averaged in time intervals of 100τ for φp < 0.03 [Fig. 7(b)]. For
φp  0.03 [Fig. 7(a)], the differentials are not averaged. The
solid line in the figure represents Smoluchowski’s expression
for the coagulation
rate, which includes unsteady-state effects,
α 0 (1 + Rij / π Dij t), like Eq. (4). Here, α 0 represents the ratioof the coagulation rate to that at steady state t  Rij2 /Dij . We
chose α 0 as the same value obtained from the fitting of Eq. (13)in Fig. 4. The relationship Rij = σ,Dij = 2D1 = 2kB T /6π μa
is also assumed. In Fig. 7(a), the initial coagulation rate
has a large value due to the unsteady-state particle spatial
distribution, as shown in Fig. 1. As time progresses, thecoagulation rate gradually decreases in a similar manner
to Smoluchowski’s form for the unsteady coagulation rate.
The coagulation rate subsequently has a local minimum

051403-7

MATSUOKA, FUKASAWA, HIGASHITANI, AND YAMAMOTO

PHYSICAL REVIEW E 86, 051403 (2012)

at approximately t = tc and then increases abruptly due tothe cluster-cluster or heterogeneous coagulations mentioned
above. These minimum values in the temporal coagulation
rate in Fig. 7 approximately correspond to the values of K/K 0
at φp  0.03 in Fig. 4. At φp < 0.03, the coagulation rates
tend to converge to constant values with time. In our case,the characteristic time constant for the unsteady-state regime
is roughly estimated at Rij2 /Dij ∼ 12π μa 3 /kB T = 404τ . As
the particle volume fraction increases, tc becomes shorter andbecomes comparable to this time constant. It is thus expectedthat at high particle volume fractions, the unsteady state in
the initial coagulation stage will be immediately followed by
the later cluster-cluster or heterogeneous coagulation process
before the unsteady state converges.
Figure 8 shows the time evolution of the number averageradius of gyration of clusters Rg :cluster growth for self-similar aggregates obey the following
scaling relationship:
Rg,j df
,
(16)
Nj ≈
a

c
1 
Rg,j ,
Nc j =1

N

Rg =

2
=
Rg,j

Nj



(14)

where df is the fractal dimension which characterizes thecluster geometric structures [34–36]. Figure 9 indicates therelationship between cluster sizes and radiuses of gyration ofclusters at φp = 0.1. Figures 9(a) and 9(b) show data for t/τ 1220 and t/τ > 1220, respectively. In the figure, lines of Eq.
(16), in which a prefactor is fixed two from the relation thatRg /σ = 2 for doublets, are shown for comparison. In Fig. 9(a),plotted data almost follow the scaling relationship. However,in Fig. 9(b), points with large cluster sizes, Nj  1000,
deviate largely from the scaling relationship. This suggests
the coexistence of very large clusters comparable to the systemsize with small ones. This was confirmed also by visualizationsof cluster configurations.
V. CONCLUSIONS2
jRi − RCM,j /Nj ,(15)

where Nc is the total number of clusters, Nj is the numberjof particles composing the j th cluster, Ri is the positionof the ith particle composing the j th cluster, and Rg,j and
RCM,j are the radius of gyration and the center of mass of thej th cluster, respectively. The calculation of Eq. (14) includes
primary particles assumed to be clusters with Rg = a. For
φp  0.03, radiuses of gyration of clusters increase abruptly
as time progresses. The time range for the growth in the radiusof gyration reasonably corresponds to those for increases ofcoagulation rates in Fig. 7(a). This clearly suggests that spatial
growth of clusters lead to increases of coagulation rate in thelater stages of coagulation. At φp = 0.1, the radius of gyration
saturates and has a peak at about t/τ = 1220. This reflects
the generation of very large clusters expanding the whole ofsystem as mentioned so far.
It is well known that the internal structure of colloidal
aggregates and coagulation types can be described by conceptsof fractal geometry. Several models and experimental results of

The rate of rapid Brownian coagulation for fractional
particle volume at dilute conditions φp = 0.003 to thatat high volume conditions φp = 0.1 is investigated using
the DNS method, which explicitly considers HIs between
particles by simultaneously solving the motions of both theparticles and the fluid. In the dilute cases, the rate of rapid
Brownian coagulation is reduced to approximately 0.3–0.5
times the theoretical value predicted by Smoluchowski. This
result is consistent with the results of previous experimental
and theoretical studies that took account of hydrodynamic
effects. Because exact HIs are explicitly incorporated in our
simulation, the reproduction of this reduction in coagulation
rate by our calculations clearly indicates that this observed
reduction in the coagulation rate is caused by hydrodynamic
effects. Moreover, the volume fraction dependence of thecoagulation rate indicates that the coagulation rate increases
with increasing volume fraction. At high particle volume
fractions, the unsteady state in the initial coagulation stage
is expected to be followed by a later, highly heterogeneous
coagulation process before the unsteady state converges.

[1] M.
Smoluchowski,Z.
Phys.
Chem.,
Stochiom.
Verwandtschaftsl 92, 129 (1917).
[2] J. W. T. Lichtenbelt, C. Pathmamanoharan, and P. H. Wiersema,J. Colloid Interface Sci. 49, 281 (1974).
[3] W. Hatton, P. Mcfadyen, and A. L. Smith, J. Chem. Soc., Faraday
Trans. 1 70, 655 (1974).
[4] K. Higashitani and Y. Matsuno, J. Chem. Eng. Jpn. 12, 460(1979).
[5] H. Sonntag and K. Strenge, Coagulation Kinetics and Structure
Formation (Plenum, New York, 1987).
[6] H. Holthoff, S. U. Egelhaaf, M. Borkovec, P.
Schurtenberger, and H. Sticher, Langmuir 12, 5541
(1996).

[7] B. V. Deryaguin and L. Landau, Acta Phys. Chem. 14, 633(1941).
[8] E. J. W. Verwey and J. T. G. Overbeek, Theory of
the Stability of Lyophobic Colloids (Elsevier, Amsterdam,
1948).
[9] N. Fuchs, Z. Phys. 87, 736 (1934).
[10] L. A. Spielman, J. Colloid Interface Sci. 33, 562 (1970).
[11] E. P. Honig, G. J. Roebersen, and P. H. Wiersema, J. Colloid
Interface Sci. 36, 97 (1971).
[12] K. Higashitani, T. Tanaka, and Y. Matsuno, J. Colloid Interface
Sci. 63, 551 (1978).
[13] T. Fukasawa and Y. Adachi, J. Colloid Interface Sci. 304, 115(2006).

i=1051403-8

EFFECT OF HYDRODYNAMIC INTERACTIONS ON RAPID . . .

PHYSICAL REVIEW E 86, 051403 (2012)[14] Y. Kusaka, T. Fukasawa, and Y. Adachi, J. Colloid Interface Sci.
363, 34 (2011).
[15] G. J. Roebersen and P. H. Wiersema, J. Colloid Interface Sci.
49, 98 (1974).
[16] M. C. Heine and S. E. Pratsinis, Langmuir 23, 9882 (2007).
[17] B. Buesser, M. C. Heine, and S. E. Pratsinis, Aerosol Sci. 40, 89(2009).
[18] G. Urbina-Villalba, J. Toro-Mendoza, A. Lozsán, and M. Garcı́aSucre, J. Phys. Chem. B 108, 5416 (2004).
[19] M. Lattuada, J. Phys. Chem. B 116, 120 (2012).
[20] J. F. Brady and G. Bossis, Annu. Rev. Fluid Mech. 20, 111(1988).
[21] D. R. Foss and J. F. Brady, J. Fluid Mech. 407, 167 (2000).
[22] H. Tanaka and T. Araki, Phys. Rev. Lett. 85, 1338 (2000).
[23] R. Yamamoto, K. Kim, Y. Nakayama, K. Miyazaki, and D. R.
Reichman, J. Phys. Soc. Jpn. 77, 084804 (2008).
[24] Y. Nakayama and R. Yamamoto, Phys. Rev. E 71, 036707(2005).
[25] J. T. Padding and A. A. Louis, Phys. Rev. E 74, 031402(2006).
[26] M. Fujita and Y. Yamaguchi, Phys. Rev. E 77, 026706 (2008).

[27] I. Pagonabarraga, B. Rotenberg, and D. Frenkel,
Phys.Chem.Chem.Phys. 12, 9566 (2010).
[28] S. A. Rice, in Diffusion-Limited Reactions, Compr. Chem.
Kinetics 25, edited by C. H. Bamford, C. F. H. Tipper, andR. G. Compton (Elsevier, Amsterdam, 1985).
[29] Y. Nakayama, K. Kim, and R. Yamamoto, Eur. Phys. J. E 26,
361 (2008).
[30] T. Iwashita and R. Yamamoto, Phys. Rev. E 80, 061402(2009).
[31] T. Iwashita, Y. Nakayama, and R. Yamamoto, J. Phys. Soc. Jpn.
77, 074007 (2008).
[32] H. C. Hamaker, Physica 4, 1058 (1937).
[33] T. M. Trzeciak, A. Podgóski, and J. C. M Marijnissen, inProceedings of the Fifth World Congress on Particle Technology(Orlando, FL, USA, 2006), p. 202d.
[34] R. Jullien, M. Kolb, and R. Botet, J. Phys. (Paris), Lett. 45, L977
(1984).
[35] D. A. Weitz, J. S. Huang, M. Y. Lin, and J. Sung, Phys. Rev.
Lett. 54, 1416 (1985).
[36] M. Y. Lin, H. M. Lindsay, D. A. Weitz, R. Klein, R. C. Ball, and
P. Meakin, J. Phys.: Condens. Matter 2, 3093 (1990).

051403-9