European Journal of Mechanics B/Fluids 50 (2015) 19–26

Contents lists available at ScienceDirect

European Journal of Mechanics B/Fluids

journal homepage: www.elsevier.com/locate/ejmflu

On the settling behaviour of polydisperse particle clouds in

viscous fluids
F. Bülow a,b,∗ , H. Nirschl a , W. Dörfler b
a
Institute of Mechanical Process Engineering and Mechanics, Karlsruhe Institute of Technology, 76131 Karlsruhe, Germany
b
Institute for Applied and Numerical Mathematics, Karlsruhe Institute of Technology, 76131 Karlsruhe, Germany

article info abstract

Article history: We investigate the behaviour of polydisperse particle clouds settling in a viscous fluid by means of
Received 31 October 2013 numerical simulations. The sedimentation process is described on examples of bidisperse systems as well
Received in revised form as fully polydisperse particle systems, pointing out parallels and differences to the monodisperse case.
20 October 2014
Based on our results we propose a simple formula which gives the maximum velocity of a polydisperse
Accepted 3 November 2014
Available online 10 November 2014
particle cloud settling in a viscous fluid. The proposed formula agrees with our numerical results even in
case of systems with high polydispersity.
Keywords:
© 2014 Elsevier Masson SAS. All rights reserved.
Fluid mechanics
Particle
Polydisperse
Sedimentation
Simulation
Stokesian dynamics

1. Introduction sedimenting cluster of particles have been identified [2]. These de-
pend on the particle Reynolds number, the cloud Reynolds number,
Sedimentation is part of many industrial processes in the fields the particle volume fraction inside the cluster and the size of the
of chemical, biological and environmental engineering. Waste wa- cluster. Research has mostly been concentrated on two regimes.
ter treatment and mineral processing are only two examples where On the one hand a regime where the cloud Reynolds number is
the treatment of highly polydisperse suspensions plays a major small but cannot be assumed to vanish and thus the particle evolu-
role. As opposed to the assumption of many models, particles are tion can be described by Oseen interactions. On the other hand the
normally not spatially uniformly distributed in such suspensions. Stokes regime, where the cloud Reynolds number is small enough
Due to hydrodynamic interactions and the polydisperse nature of to be assumed to be zero. In this regime, some particles can make
the particulate system, clusters of different sizes and lifetimes are an outward crossing of the settling cloud and leave it at its rear
formed. Such clusters, or particle clouds, can undergo a remark- end in a vertical tail. The rate of leakage, which describes the rate
able motion. In the absence of turbulence, spherical clouds behave
of leaked particles per time, has been studied extensively [3,4]. In
like drops of fluid settling in another fluid with lower mass den-
case of a higher but still small cloud Reynolds number the leakage
sity [1]. During sedimentation such a cloud evolves into a torus. The
decreases to an insignificant number of particles, while the evo-
inhomogeneous distribution of particles in the torus causes oscilla-
lution of the shape of the cloud remains the same [5]. Apart from
tions which eventually lead to a break-up into two or more smaller
qualitative investigations and descriptions of the shape evolution
clouds, which then undergo the same evolution. This process con-
of particle clouds, most works deal with the sedimentation veloc-
tinues until the resulting clouds are so small that they disintegrate.
In recent years such particle clouds have drawn a considerable ity of particle clouds. Adachi et al. [6] were among the first who
amount of attention and have been investigated both experimen- examined the cloud evolution experimentally and gave a theoreti-
tally as well as numerically. Different regimes of evolution for a cal prediction for the sedimentation velocity of a cloud comprised
of a large number of particles. With the increase in computational
power, numerical simulations of settling clouds of many particles
∗ Corresponding author at: Institute of Mechanical Process Engineering and have become feasible. Nitsche and Batchelor [4] proposed a sim-
Mechanics, Karlsruhe Institute of Technology, 76131 Karlsruhe, Germany. Tel.: +49
ple formula for the sedimentation velocity of a spherical particle
72160842428. cloud, which was in good agreement with results from their nu-
E-mail address: [email protected] (F. Bülow). merical simulations of monodisperse particle clouds comprised of
http://dx.doi.org/10.1016/j.euromechflu.2014.11.003
0997-7546/© 2014 Elsevier Masson SAS. All rights reserved.
20 F. Bülow et al. / European Journal of Mechanics B/Fluids 50 (2015) 19–26

up to 320 particles. Bosse et al. [7] have simulated particle clouds where FH , TH , UP , ωP ∈ R3Np are the same variables as in Eqs. (1)
of several thousand particles, yielding results comparable to those and (2), with the difference that the vectors in Eq. (3) contain the
of Nitsche and Bachelor and their proposed formula in case of a components of all Np particles. The mobility matrix M ∈ R6Np ×6Np
low cloud Reynolds number. Bosse et al. [7] also found that for depends only on the particle configuration and the particle radii
an increasing cloud Reynolds number, the number of emerging explicitly. The fluid motion does not have to be computed, which is
secondary clouds rises. Apart from this the overall process stays a major advantage over grid-based methods as for instance a cou-
the same—a repeating cascade starting from an initially spherical pling of a finite volume method and the discrete element method.
cloud, which flattens, evolves into a torus eventually breaking up A nondimensionalization of the equations leads to the definition of
into smaller clouds, which then evolve in the same manner. A num- the Stokes number, which in our case is given by
ber of other works deal with experiments, numerical simulations
or theoretical considerations on particle clouds in viscous fluids 2 ρp
St = Re. (4)
[8,9,3]. They all discuss the topics of cloud behaviour, initial cloud 9 ρf
velocity and particle leakage. But they confine themselves to inves- ρp and ρf are the particle mass density and the fluid mass den-
tigations on monodisperse particle clouds. Abade and Cunha [10] sity, respectively. Re denotes the Reynolds number of the flow. If
have investigated the behaviour of slightly polydisperse particle we can assume Stokes flow with vanishing Reynolds number, we
clouds. They found that polydisperse clouds comprised of a small can assume
number of particles are less stable than comparable monodisperse
  the Stokes number to vanish as well, provided that
ρp ≈ O ρf . With this assumption the left hand sides of the dimen-
clouds. Apart from clouds of spherical particles, some works con- sionless form of Eqs. (1) and (2) become zero, resulting in Fh = −Fg
sider clouds of fibres [11] and fibre–particle mixtures [12]. Nev- and Th = 0. Substitution into Eq. (3) now yields a linear rela-
ertheless, we are not aware of any works explicitly investigating tion between gravity on the one hand, and particle velocities on
the settling behaviour of polydisperse clouds comprised of more the other. The resulting algorithm reads as follows: (i) compute M
than 1000 particles, neither numerically nor experimentally. The from particle positions and their radii, (ii) obtain the particle veloc-
present work deals with the behaviour of such clouds at vanishing ities from the substituted form of Eq. (3), (iii) integrate UP , which
Reynolds number and shall serve as a basis for further research in contains all translational velocities, to yield the new particle po-
this area. Basis for our findings is our parallel Stokesian Dynam- sitions. The problem we solve – obtaining particle velocities from
ics code, which is based on the open source software RYUON [13]. acting forces – is called mobility problem. For the construction of
Stokesian Dynamics is a numerical method developed and espe- the mobility matrix M in (3) higher moments in the expansion of
cially well-suited for the simulation of particle suspensions [14,15]. the force density in the expression for the fluid velocity field in
Our work is structured as follows. First we give a short introduc- Stokes flow are neglected [14]. We only keep the zeroth moment,
tion to the simulation method and we derive a formula yielding the the force FH , and the antisymmetric part of the first moment, the
maximum velocity of a polydisperse particle cloud settling at van- torque TH . This leads to an error when particles are in close prox-
ishing Reynolds number. This formula agrees well with our numer- imity, as it usually is the case in dense particle suspensions. Closely
ical simulations, as we show in the results section. Furthermore we related to the mobility problem is the so-called resistance prob-
present results of our investigations on the motion of polydisperse lem. It has to be solved when particle velocities are given and the
clouds on examples of bidisperse and fully polydisperse clouds. forces are unknown. This is usually the case when dense suspen-
sions are considered and lubrication effects cannot be neglected.
2. Methods There are several ways to include near-field lubrication effects for
this type of problem [18–20]. We chose to solve the mobility prob-
In this section we describe our simulation method and we de- lem without lubrication correction rather than solving the resis-
rive a formula which gives the settling velocity of a polydisperse tance problem with lubrication correction for the simple reason
spherical particle cloud. We confine ourselves to a brief discus- that the resulting algorithm is much faster than a solution of the
sion of the exact form of the numerical scheme we use, since the resistance problem and yet exact enough to capture the dynamics
method and extensions to it are well known. For readers not famil- of a particle cloud. It is also possible to include near field effects
iar with the Stokesian Dynamics method, we refer to [14,16,17]. by an introduction of a short-range repulsive force, as in [4,10,21].
This method of including hydrodynamic near-field interactions has
2.1. Numerical scheme the advantage that no inversion of the mobility matrix is neces-
sary while near-field effects are still taken into account. However,
We consider particles settling in a viscous fluid, thus the Stokes our numerical tests of two particles interacting hydrodynamically
equations suffice to describe the fluid motion. When interparticle as in [22,23] have shown that, at least with the proposed choice
forces are negligible, the equations of motion for suspended non- of parameters, this treatment produces unrealistic results. In the
Brownian particles solely under the influence of gravity are following we will call an algorithm solving the mobility problem
mobility scheme and a method for the solution of the resistance
d
mp Up = Fh + Fg , (1) problem resistance scheme, respectively.
dt
and 2.2. Cloud settling velocity of a polydisperse cloud
d
Ip ωp = Th . (2) Nitsche and Batchelor [4] have proposed a formula for the set-
dt
tling velocity Uc of a monodisperse spherical particle cloud at low
Here mp is the mass of a particle, Up its translational velocity, Fh
Reynolds number,
the hydrodynamic force acting on the particle, and Fg gravity. Ip
denotes the moment of inertia tensor of a particle, ωp its angular
 
6 ap
velocity, and Th the hydrodynamic torque acting on the particle.
Uc ≈ USt Np +1 , (5)
5 ac
Based on the integral expression for the fluid velocity field in Stokes
where USt denotes the Stokes velocity of a single particle with ra-
flow [15], a linear relation between hydrodynamic force, hydrody-
dius ap and ac is the radius of the particle cloud formed by the Np
namic torque and the translational and angular velocities of parti-
particles. Eq. (5) has been derived from the Hadamard–Rybczynski
cles suspended in a fluid can be found. This relation is given by
    equation (H–R equation). The H–R equation describes the termi-
FH UP nal velocity of a spherical drop in a quiescent fluid under the ac-
M = , (3)
TH ωP tion of gravity [24]. It is insignificant for the derivation whether
 F. Bülow et al. / European Journal of Mechanics B/Fluids 50 (2015) 19–26 21

Fig. 1. Particle configuration and relative deviation of our results using the mobility scheme from the results of Goldman et al. [22] against the reciprocal dimensionless
inter-particle spacing r.

the matter forming the spherical drop has a higher or lower mass 3.1.1. Error made in one time step
density than the suspending fluid. Thus, based on the observation To determine the accuracy of our simulation, we chose the
that particle clouds behave like settling drops of fluid [1], it seems example of two equal spheres settling in a given configuration in
natural to utilize the Hadamard–Rybczynski equation in order to a quiescent fluid ([22]; Fig. 1). This choice arises from the fact that
derive a formula which gives the settling velocity of spherical par- this example covers different particle configurations relative to the
ticle clouds. In order to derive a formula yielding the maximum settling direction, as well as different interparticle distances. The
settling velocity of a polydisperse spherical particle cloud, we start particle settling velocity obtained from the solution of the mobility
from the Hadamard–Rybczynski equation as well. The cloud set- problem without lubrication correction starts to deviate from the
tling velocity Uc , as we derive it, is given by results of Goldman et al. [22] by more than 1% at a distance of
Nf
one particle diameter between the particle surfaces, where the
6 ap,i deviation increases with shrinking gap size. Note that the deviation
Uc = Np,i USt ,i . (6)
5 i=1 ac is smaller for larger angles ϕ , i.e. particles settling side by side
(Fig. 1). The deviation is below 5% for all configurations, down to
USt ,i is the Stokes velocity of a spherical particle with radius ap,i , inter-particle gaps of 1% of a particle diameter.
Nf is the number of particle size fractions in the system forming
a cloud of radius ac . Fraction i contains Np,i particles. A detailed
3.1.2. Accuracy over a longer time span
derivation of formula (6) is given in the Appendix. The only as-
sumption for this formula is that the particle cloud is dilute. This In the previous example where we examined the error made
can be described by φ ≪ 1, where φ is the particle volume fraction in one time step, particles do not move relative to each other. To
in the spherical cloud. Because of this requirement, Eq. (6) seems investigate the deviation of the solution of the mobility problem
to be valid only for cases when φ is very small. But as we will show without lubrication correction from the solution of the resistance
in Section 3, it gives a good approximation for a broad spectrum problem with lubrication correction after a longer time span, we
of particle volume fractions. As a comparison of Eqs. (5) and (6) have run numerical experiments as described in [15]. Durlofsky
shows, assuming a monodisperse particle system, Eq. (5) contains et al. [15] conducted numerical simulations of the settling of four
an additional term equal to the settling velocity USt of an isolated and of eight spheres initially located at the corners of a square
particle. Using a different approach, Ekiel-Jezewska et al. [9] have and a cube, respectively. Under the action of gravity and at van-
derived a formula for the settling velocity of monodisperse parti- ishing Stokes number both configurations undergo a periodic mo-
cle clouds as well. Their formula also contains the additional term, tion. This motion has first been described by Hocking [25] for the
which they attribute to a slip between particles and ambient fluid. four-particle configuration. The spheres initially located at the top
Apart from this term, Eq. (6) coincides with Eq. (5) in the monodis- corners settle faster than the particles initially located at the bot-
perse case. We discuss this topic in Section 3.4. If all fractions of the tom corners. After the faster settling particles pass the slower par-
particulate system are known, e.g. in numerical experiments, the ticles all particles return to the corners of the initial square or
velocity resulting from Eq. (6) can be used as characteristic value cube. Only this time the locations of top and bottom particles are
for nondimensionalization. This yields a more stable simulation, switched. Now the particles initially located at the bottom corners
which saves time and resources. settle faster and return to their original positions. As in [15] we de-
termined the vertical position when all particles form a horizontal
3. Results and discussion line in case of four particles, and when they lie in a plane perpen-
dicular to the settling direction in case of eight particles initially lo-
3.1. Validation cated at the corners of a cube. The considered dimensionless edge
lengths are rinit = 3, 4, 5, the used level of approximation is the
We have validated our code on several examples, obtaining described FT-version. Our results using the resistance scheme with
good agreement. In this section we give some examples which are lubrication correction as in [15] agree very well with their results.
significant for this work. For all numerical experiments we assume When we use the proposed mobility scheme without lubrication
that the particle density is the same for all particles. The driving correction, our results deviate. This can be seen in Fig. 2. In case
force in all simulations is gravity, which scales as the third power of eight settling particles (Fig. 2, right), the deviation is larger for
of the particle radius. particles with shorter initial distance. This is because lubrication
22 F. Bülow et al. / European Journal of Mechanics B/Fluids 50 (2015) 19–26

Fig. 2. Left: Particle positions of four settling equal particles initially located at the corners of a square with initial edge length rinit = 3, 4, 5 when all particles form a line
perpendicular to the settling direction. Right: Particle positions of eight settling equal particles initially located at the corners of a cube with initial edge length rinit = 3, 4, 5
when all particles lie in a plane perpendicular to the settling direction. Only four of the eight particles are shown because of the reflectional symmetry in the (x, z)- and
(y, z)-planes. All lengths are given in dimensionless units.

a b with the results of Metzger et al. [3] shows good agreement, which
supports our choice to solve the mobility problem rather than the
time-consuming resistance problem. The characteristic stages of
the settling process are shown in Fig. 3, where (a) shows the ini-
tially spherical cloud, (b) the formed torus (including top view),
(c) the particle system after the torus broke into smaller clouds, and
c d finally (d) even smaller clouds, formed after the secondary clouds
broke up. To nondimensionalize the simulated system we use a
characteristic length scale L0 and a characteristic velocity U0 . For
the simulation of bidisperse systems we chose the particle radius
of the large size fraction as L0 . For simulations of fully polydisperse
systems with a given radius distribution, L0 was chosen such that it
was slightly larger than the mean of the distribution. As U0 served
Fig. 3. Characteristic stages during the sedimentation of a monodisperse particle the Stokes velocity corresponding to the chosen L0 . To perform in-
cloud with initial particle volume fraction of 5%. Particles in (c) and (d) have been
tegration in time we applied an embedded Runge–Kutta–Fehlberg
enlarged by a factor of 3 for better visibility. The vertical tail of leaked particles is
not shown. method with adaptive step-size control.

effects play a bigger role for the configuration with shorter initial 3.2. Evolution of a bidisperse cloud
distance. In the graph on the left-hand side of Fig. 2, where the ver-
tical positions of four particles are displayed when they first form a In order to investigate the nature of the cloud evolution for
horizontal line, one can see that the deviation is very small for the a bidisperse particle system we have investigated two different
square with shortest edge length, whereas the square with longest kinds of initial conditions—unmixed bidisperse clouds, where the
edge length, where lubrication effects should play a smaller role, is fraction of large particles is only distributed in the upper hemi-
larger. This results from the fact that the particles initially located sphere of the cloud, and well-mixed clouds, where no restriction
at the corners of a square with edge length 5 settle over a longer to the location of the particles inside the cloud has been imposed.
distance until they form a line perpendicular to the settling direc-
tion. Thus the accumulated error is larger, even though the error 3.2.1. Unmixed bidisperse cloud
in one time step is smaller than for a square with edge length 3. Fig. 4 shows the characteristic first stages of a bidisperse parti-
Solving the resistance problem with lubrication correction rather cle cloud consisting of 1000 large and 1000 small particles, where
than the mobility problem without lubrication correction might the particles are placed according to a uniform random distribution
lead to more precise results, but the computational time needed and a given volume fraction φ inside the initially spherical cloud.
to achieve them is several orders of magnitude higher, especially The restriction in this study is that large particles are only placed in
for polydisperse systems. At the current state of the art, this re- the upper hemisphere and small particles only in the lower hemi-
stricts simulations using a lubrication correction scheme based on sphere. We chose this initial configuration because it promotes
the resistance problem to the simulation of several hundreds to mixing of the two fractions during sedimentation more than other
some thousands of particles within days or months. If we are in- comparable configurations. As depicted in Fig. 4 the fraction of
terested in macroscopic properties such as the overall motion of a large particles forms a cloud, dives through the portion of small
particle cloud, the solution of the mobility problem without lubri- particles, and takes a number of small particles with it. Fig. 5 shows
cation correction suffices (see e.g. [3] where a point force model that the higher the radius ratio λ of small to large particles, the
was used). Provided that the code is parallelized, the method we more small particles get carried along by the fraction of large par-
chose allows for the simulation of the dynamics of several thou- ticles. For values of λ lower than 0.5 the difference between the
sand up to a million particles within days [26]. For such particle resulting curves is not as distinct as it is for higher radius ratios.
systems different effects can be observed compared to simulations The smaller a particle is, the less impact it has hydrodynamically
of relatively small particle systems (c.f. Section 3.3). on a larger particle. This effect is reflected by the entries of the mo-
bility matrix M in Eq. (3). Furthermore, as gravity scales with the
3.1.3. Settling of a monodisperse particle cloud third power of the particle radius, the driving force for the settling
A qualitative comparison of our results from a simulation of the of particles is lower for a smaller particle size. These two factors are
settling process of a monodisperse particle cloud of 3000 particles the main reasons leading to the described behaviour. Even though
 F. Bülow et al. / European Journal of Mechanics B/Fluids 50 (2015) 19–26 23

a b

c d
Fig. 6. Number of large particles Np,large inside an initially unmixed settling cloud
of Np = 2000 particles for an initial volume fraction φ = 0.075.

clouds form a torus and then break up into smaller clouds for all
values of ν . The time of the break-up varies, depending on the
exact initial particle configuration. Moreover, the break-up of the
particle cloud itself depends on the exact initial positions of the
particles. Fig. 7(a) shows the initial configuration of a well-mixed
particle cloud with radius ratio λ = 0.5 and an initial volume
fraction φ = 0.05. Fig. 7(b)–(d) shows snapshots of the break-
Fig. 4. Initial particle configuration in an unmixed cloud and subsequent stages.
up for three different initial particle configurations. The only dif-
The radius ratio of small to large particles is λ = 0.5, the initial volume fraction of
particles inside the spherical cloud is φ = 0.05 and the total number of particles is ferences between the clouds shown in Fig. 7(b)–(d) are the exact
Np = 2000. initial positions of the particles inside the spherical cloud domain,
which are determined by a Monte Carlo method. Differing from the
monodisperse case, a polydisperse particle cloud not always breaks
up into two smaller approximately evenly sized clouds. Instead it
can break up into one large and another small cloud, into two sim-
ilar clouds, or even into three clouds. The break-up into more than
two clouds has thus far only been observed for break-up events
at higher Reynolds numbers [7] and in experiments as well as in
numerical studies using a point-particle model [27]. In [27] poly-
dispersity was not analysed explicitly. In reality, particle systems
are never truly monodisperse. Our numerical results show that the
polydispersity seems to have an effect on the settling process. The
fraction of large particles determines the motion of the cloud up to
ν ≈ 75%. At this value the dominating fraction was indeterminable
in our numerical experiments. For percentages higher than 75%,
the fraction of small particles seems to determine the cloud be-
Fig. 5. Number of small particles Np,small inside an initially unmixed settling cloud
haviour. By the expression ‘determine the cloud behaviour’ we
of Np = 2000 particles for an initial volume fraction φ = 0.075. mean that the respective fraction forms a cloud and drags the par-
ticles of the other fraction with it. The evolution of the number of
smaller particles have less hydrodynamical impact on large parti- small particles inside the settling cloud, Np,small , for the test runs for
cles, the many small particles left behind for e.g. a radius ratio of which the first break-up event is shown in Fig. 7(b)–(d) is displayed
λ = 0.1 can trap large particles, which then cannot settle with the on the right hand side of the figure. Compared with unmixed clouds
cloud anymore. This effect can be seen in Fig. 6, where the evolution (cf. Fig. 5) the loss of small particles in the initial phase is not as pro-
of the number of large particles Np,large is shown. With increasing nounced. More small particles remain in the cloud until the first
radius ratio more small particles stay in the cloud such that e.g. at break-up event.
λ = 0.9 no trapping occurs. Figs. 5 and 6 show the number of par-
ticles inside the cloud until the first break-up event. At the point 3.3. Evolution of a fully polydisperse cloud
of break-up most small particles still trapped inside the particle
cloud are freed and start to settle slower. In order to further investigate the behaviour of polydisperse
particle clouds, we have conducted numerical experiments on fully
3.2.2. Mixed bidisperse cloud polydisperse systems with a particle radius distribution accord-
We have conducted numerical experiments on the probability ing to a given log-normal distribution. The considered clouds con-
of a break-up of a bidisperse particle cloud. This process is shown sisted of 1000, 2500 and 5000 particles, respectively, with particle
for the monodisperse case in Fig. 3 and has been reported in the and fluid properties such that the resulting Reynolds number and
mentioned literature [3,4]. Differing from the numerical experi- Stokes number were sufficiently small, i.e. the assumptions Re ≈ 0
ments described in the last subsection, particles inside the initial and St ≈ 0 hold. For a particle volume fraction of φ = 0.01 in a
cloud are now well mixed. We simulated the settling process of spherical cloud of 5000 particles, the characteristic stages of the
particle clouds of 1000 particles with a radius ratio of λ = 0.5 and cloud evolution are shown in Fig. 8. The depicted particle system
a volume fraction of φ = 0.05 in the initially spherical cloud, while has a radius distribution with resulting dimensionless expected
increasing the percentage ν of small particles in the initial cloud in value µ ≈ 0.72 and standard deviation σ ≈ 0.15 (correspond-
steps of 5% from 0% to 95%. In ten out of ten runs the described ing to the lognormal distribution LN (−0.35, 0.04)). Abade and
24 F. Bülow et al. / European Journal of Mechanics B/Fluids 50 (2015) 19–26

Fig. 7. Left: Initial particle configuration for φ = 0.05 in a mixed cloud (a) and snapshots of the break-up event for different initial configurations (b)–(d). Right: Number
of small particles Np,small inside the clouds shown in (b)–(d) until the first break-up. Run 1 corresponds to (b), run 2 to (c) and run 3 to (d). The total number of particles is
Np = 1000.

a b

c d

Fig. 8. Characteristic stages during the sedimentation of a polydisperse particle

cloud of Np = 5000 particles. Particles in (c) and (d) have been enlarged by a factor
of 3 for better visibility. Fig. 9. Comparison of values for the settling velocity of a bidisperse cloud obtained
from our numerical simulations, Uc ,n , and the prediction Uc ,a from formula (6). Here
φ = 0.05 and Np = 1000. The settling velocities have been nondimensionalized by
Cunha [10] state that polydisperse particle clouds are less stable the Stokes velocity U0 of a particle of the fraction of large particles.
than comparable monodisperse clouds. They disintegrate before a
torus can be formed. The clouds investigated in their work con- systems. Tested values for φ are 0.01, 0.025, 0.03, 0.04, 0.05, 0.06,
sisted of 300 particles. According to Metzger et al. [3] this number 0.07, 0.08, 0.09, 0.1, 0.15 and 0.2. The radius distribution was cho-
of particles is also in the monodisperse case too low for the cloud sen for all simulations such that the resulting dimensionless ex-
to undergo the characteristic evolution, whereas the likelihood for pected value was µ ≈ 0.72. Together with a standard deviation
clouds comprised of 1000 and more particles is about 100%. Our of σ ≈ 0.15 this yields a reasonably polydisperse radius spectrum
results show that also polydisperse clouds undergo this evolution, (the corresponding distribution is LN (−0.35, 0.04)). A compar-
if they consist of enough particles. ison of the cloud settling velocity of a bidisperse cloud, obtained
from our numerical simulations, and the corresponding value from
3.4. Cloud settling velocity Eq. (6) is depicted in Fig. 9. The volume fraction of particles in the
cloud was φ = 0.05 with a total number of Np = 1000 particles.
For a validation of formula (6) we have computed the mean The error bars showing the deviation due to different initial config-
of the particle velocities in our numerical experiments at a stage urations in runs with the same set of parameters show that there is
where the particle cloud still forms a sphere. At this time no par- a statistical component. Yet, the agreement of prediction and nu-
ticles have left the cloud and the cloud attains its maximum set- merical experiment is very good for all radius ratios λ and all values
tling velocity. The error bars along the curves in Figs. 9, 10 and 12 of ν . A comparison of our numerical results and the proposed pre-
result from the uniform random placement of particles inside the diction (6) for the case of a fully polydisperse cloud is depicted in
initial cloud, which influences the cloud’s settling velocity. How- Fig. 10. For these results we have varied the number of particles per
ever, this influence is small compared to the total settling velocity cloud, Np , and the solid volume fraction φ in the clouds. The devi-
of the cloud, as one can see by the width of the error bars. Their ation of the cloud settling velocity from the computed mean is ex-
width goes to zero as the number of small particles in the cloud is tremely low for all tested Np and φ . Furthermore, Fig. 10 shows that
increased, which indicates that the exact placement of particles in prediction and numerical results agree very well in case of low val-
the cloud loses importance with a growing number of small parti- ues of φ , what corresponds to a dilute cloud. The deviation between
cles. For every set of parameters we have carried out four simula- the corresponding curves increases for higher volume fractions.
tions to obtain the error bars. This is due to the assumption of a dilute cloud, which has been
made in the derivation of formula (6). Nevertheless, this formula
3.4.1. Velocity of polydisperse clouds still gives a good prediction, even for solid volume fractions of 20%.
To cover a wide spectrum of possible particle systems, we have
run extensive numerical tests on the already described well-mixed 3.4.2. Velocity of a monodisperse cloud
bi- and polydisperse particle systems. While we varied the per- In addition to the comparisons of analytical prediction and nu-
centage of small particles in bidisperse particle systems with fixed merical results presented in Figs. 9 and 10 we have run simulations
volume fraction, we varied φ in case of fully polydisperse particle of monodisperse clouds with a varying number of particles Np and
 F. Bülow et al. / European Journal of Mechanics B/Fluids 50 (2015) 19–26 25

Fig. 12. Comparison of different estimates for the settling velocity of a polydisperse
Fig. 10. Comparison of the particle cloud settling velocity for varying φ and differ-
cloud for Np = 5000 and varying φ . Uc ,num denotes our numerical results, Uc ,ana,poly
ent Np in a fully polydisperse cloud. Uc ,n denotes values from the numerical simula-
the use of formula (6) considering polydispersity, and Uc ,ana,mono the use of formula
tion, Uc ,a the prediction from formula (6). U0 has been defined in Section 3.1.3.
(6) assuming a monodisperse particle system. U0 has been defined in Section 3.1.3.

4. Conclusions

This work gives more insight into the sedimentation process

of polydisperse particle clouds in viscous fluids. Until now re-
search has been concentrated on the investigation and description
of monodisperse particle clouds, drops or blobs, as they are some-
times called. Due to the fact that we have parallelized our algo-
rithm, we are able to simulate particle systems of several thousand
particles with any radius distribution in a short time. This enables
us to investigate large polydisperse particle clouds, for which the
evolution differs from that of relatively small particle clouds com-
prised of several hundred up to one thousand particles. One re-
sult is that, given that the considered particle cloud consists of
Fig. 11. Comparison of the particle cloud settling velocity for varying φ and differ-
a large number of particles, a polydisperse particle cloud under-
ent Np in a monodisperse cloud. Uc ,n denotes values from the numerical simulation, goes the same remarkable motion as a sufficiently large monodis-
Uc ,a the prediction from formula (6). Particle velocities have been nondimensional- perse cloud does. An initially spherical cloud evolves into a torus
ized by the Stokes velocity U0 of a single particle. that starts to oscillate until it breaks up into smaller clouds, which
then undergo the same motion. For a smaller number of particles,
a varying initial solid volume fraction φ in the clouds. The tested namely about 1000 particles and less, polydisperse clouds are not
number of particles and initial volume fractions are the same as in as stable as their monodisperse counterpart. In an extreme case of
the fully polydisperse case (c.f. Section 3.4.1, Fig. 10). As in the poly- polydispersity, a bidisperse system, the evolution is dominated by
disperse case we obtain a good agreement of the prediction from the fraction with the highest number of particles in the system. Due
formula (6) and our numerical results, as one can see in Fig. 11. Yet
to a loss of particles during sedimentation, initially spherical par-
contrary to the polydisperse case, the additional USt -term in Eq. (5)
ticle clouds attain their maximum settling velocity while they still
would increase the agreement of the analytical prediction by for-
have their initial shape. Starting from the Hadamard–Rybczynski
mula (6) and the numerically computed cloud settling velocity.
equation we have derived a formula which yields the settling ve-
locity of such a particle cloud. By a comparison with numerical ex-
3.4.3. On formulas for monodisperse systems applied to polydisperse
periments we have shown that the proposed formula gives a very
systems
good estimate for the settling velocity of a polydisperse spherical
The use of e.g. the mean particle radius together with a for-
particle cloud.
mula proposed for monodisperse clouds leads to an error, which
can be significant when polydisperse clouds are considered. Fig. 12
shows our numerical results compared to values obtained from for- Acknowledgement
mula (6), once considering the given polydisperse radius distribu-
tion and another time assuming a monodisperse particle system. F. Bülow acknowledges partial support by the Deutsche
For the latter we use the mean particle radius of the polydisperse Forschungsgemeinschaft (DFG) under grant Ni 414/14-1.
particle system described in Section 3.4.1 together with formula (6)
which in this case differs from Eq. (5) only in the additional term Appendix. Derivation of the cloud settling velocity
equal to the settling velocity of a single particle (cf. Section 2.2).
When polydispersity is not considered, the equation gives for a We start from the Hadamard–Rybczynski equation (H–R equa-
polydisperse system of 5000 particles a resulting cloud velocity tion),
with a relative deviation from the numerical value of 10% and more
for all volume fractions φ , while Eq. (6) considering polydisper- 2 ρf − ρc a2c g λ + 1
 
sity agrees very well with the value obtained from the simulation. Uc = , (A.1)
3 µf 3λ + 2
Therefore, an equation derived for monodisperse particle systems,
like Eq. (5), should only be used for truly monodisperse particle with the radius of the spherical cloud, ac , the dynamic fluid viscos-
clouds. If the particle size fractions of the system are known or ity µf and λ = µc /µf , the ratio of dynamic viscosity of the cloud
can be approximated, formula (6) which considers polydispersity to dynamic fluid viscosity. ρf and ρc are the mass density of the
yields good estimates. fluid and the mass density of the cloud, respectively. g denotes the
26 F. Bülow et al. / European Journal of Mechanics B/Fluids 50 (2015) 19–26

gravitational constant. Utilizing (ρf − ρc ) = φ(ρf − ρp ), where φ References

is the volume fraction of the particles with mass density ρp in the
spherical cloud, the H–R equation is equivalent to [1] G. Machu, W. Meile, L. G. Nitsche, U. Schaflinger, J. Fluid Mech. 447 (2001)
299–336.
[2] G. Subramanian, D. Koch, J. Fluid Mech. 603 (2008) 63–100.
2 φ ρf − ρp a2c g λ + 1
 
Uc = . (A.2) [3] B. Metzger, M. Nicolas, E. Guazzelli, J. Fluid Mech. 580 (2007) 283–301.
3 µf 3λ + 2 [4] J.M. Nitsche, G.K. Batchelor, J. Fluid Mech. 340 (1997) 161–175.
[5] F. Pignatel, M. Nicolas, E. Guazzelli, J. Fluid Mech. 671 (2011) 34–51.
If φ and λ are known, this is the equation we are looking for. Usu- [6] K. Adachi, S. Kiriyama, N. Yoshioka, Chem. Eng. Sci. 33 (1978) 115–121.
[7] T. Bosse, L. Kleiser, C. Härtel, Phys. Fluids 17 (2005) 1–17.
ally λ is not known a priori. Making the assumption φ ≪ 1 yields [8] S. Alabrudzinski, M.L. Ekiel-Jezewska, D. Chehata-Gomez, T.A. Kowalewski,
λ ≈ 1. Together with Phys. Fluids 21 (2009) 073302.
[9] M.L. Ekiel-Jezewska, B. Metzger, E. Guazzelli, Phys. Fluids 18 (2006) 038104.
Nf  3 [10] G.C. Abade, F.R. Cunha, Comput. Methods Appl. Mech. Engrg. 196 (2007)
 ap,i
φ= Np,i , (A.3) 4597–4612.
[11] J. Park, B. Metzger, E. Guazzelli, J. Butler, J. Fluid Mech. 648 (2010) 351–362.
ac
i =1 [12] M. Feist, F. Keller, H. Nirschl, W. Dörfler, Can. J. Chem. Eng. 89 (2011) 682–690.
[13] K. Ichiki, RYUON (2013) http://www.ryuon.sourceforge.net.
this yields [14] J.F. Brady, G. Bossis, Annu. Rev. Fluid Mech. 20 (1988) 111–157.
[15] L. Durlofsky, J.F. Brady, G. Bossis, J. Fluid Mech. 180 (1987) 21–49.
Nf [16] K. Ichiki, J. Fluid Mech. 452 (2002) 231–262.
6 ap,i
Uc = Np,i USt ,i . (A.4) [17] A. Sierou, J.F. Brady, J. Fluid Mech. 448 (2001) 115–146.
5 i =1 ac [18] B. Cichocki, M.L. Ekiel-Jezewska, E. Wajnryb, J. Chem. Phys. 111 (1999)
3265–3273.
[19] D.J. Jeffrey, Phys. Fluids A 4 (1) (1992) 16–29.
Here, [20] D.J. Jeffrey, Y. Onishi, J. Fluid Mech. 139 (1984) 261–290.
[21] F.R. Cunha, G.C. Abade, A.J. Sousa, E.J. Hinch, J. Fluids Eng. 124 (2002) 957–968.
2 (ρf − ρp )ap,i g
2
[22] A.J. Goldman, R.G. Cox, H. Brenner, Chem. Eng. Sci. 21 (1966) 1151–1170.
USt ,i = (A.5) [23] C. Pozrikidis, Acta Mech. 194 (2007) 213–231.
9 µf [24] E.E. Michaelides, Particles, Bubbles & Drops, World Scientific, Singapore, 2006.
[25] L. Hocking, J. Fluid Mech. 20 (1964) 129–139.
is the Stokes velocity of a spherical particle with radius ap,i . Np,i is [26] F. Bülow, P. Hamberger, H. Nirschl, W. Dörfler, Comput. & Fluids (2014).
the number of particles in one of the Nf particle size fractions of [27] A. Myłyk, W. Meile, G. Brenn, M.L. Ekiel-Jeżewska, Phys. Fluids (1994-present)
the particle system. 23 (2011) 063302.