Simulation of sedimentation and fluidization of polydisperse suspensions via a Markov model

Chemical Engineering Science 61 (2006) 5575 – 5589 www.elsevier.com/locate/ces

Simulation of sedimentation and fluidization of polydisperse suspensions via a Markov model Monika Bargieł a,∗ , Elmer M. Tory b a Institute of Computer Science, AGH University of Science and Technology, al. Mickiewicza 30, 30-059 Kraków, Poland b Department of Mathematics and Computer Science, Mount Allison University, Sackville, NB, Canada E4L 1G6

Received 18 March 2006; received in revised form 26 April 2006; accepted 26 April 2006 Available online 9 May 2006

Abstract The particle-based approach to sedimentation is extended to include velocity fluctuations that result in hydrodynamic diffusion. The vector process describing the joint values of position and velocity is Markov. Thus, no integration of velocity is required. Height–velocity “skeletons” for each particle are generated from a bivariate-normal distribution with means, variances, and covariance that depend on three parameters. For each particle, there is a unique region in which the vector of species concentrations determines that particle’s parameters and hence its Markov process, but the concentrations in that region depend on the Markov processes of neighboring particles. Though only discrete values of height and velocity are generated, the model ensures that sample paths and particle velocities are continuous. Furthermore, steady-state velocities are normally distributed and velocity autocorrelations decay exponentially. Published experimental results indicate that both are excellent approximations. For polydisperse suspensions, the Markov model is much simpler than the standard hydrodynamic-diffusion model and represents the actual process much better. We simulate the sedimentation and fluidization of polydisperse suspensions and study the effects of two additional parameters: variance and autocorrelation decay rate of particle velocities. 䉷 2006 Elsevier Ltd. All rights reserved. Keywords: Simulation; Sedimentation; Fluidization; Hydrodynamic diffusion; Polydisperse suspensions

1. Introduction For many decades, theoretical treatments of sedimentation and fluidization have modeled suspensions as continua (Bustos et al., 1999; Drew and Passman, 1999). In the continuum model, the concentration of the kth species,
k , satisfies the system of first-order conservation laws (Berres et al., 2005)  j
j
k + q
k + fk (
) = 0, jt jz

k = 1, . . . , K,

(1)

where t is time, z is height, q is the applied fluidization velocity (zero for batch sedimentation),
= (
1 ,
2 , . . . ,
K )T is the vector of solid concentrations (volume fractions), and fk =
k vk is the flux density of the kth species in batch sedimentation. Let
=
1 +
2 + · · · +
K be the total solid concentration. Then, ∗ Corresponding author. Tel.: +48 12 617 3499; fax: +48 12 633 9406.

E-mail addresses: [email protected] (M. Bargieł), [email protected] (E.M. Tory). 0009-2509/$ - see front matter 䉷 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2006.04.036

sedimentation is the evolution of
(z, t), 0 < z
H, t 0, from
0 (z) to
max (z), where
0 (z) is the initial value of
, and
max is the value when
=
max (z), the maximum total solid concentration attained in the polydisperse system. This evolution is governed by the solid flux vector f = (f1 , f2 , . . . , fK )T . The two essential ingredients for predicting this evolution are a model equation vk (
), k = 1, . . . , K (a means of predicting settling velocities from solid concentrations) and a method of implementing the changes produced by the flux. (Upward velocities and displacements are positive.) Fluidization also involves transitions in
, but the emphasis is on steady-state values. Recently, sophisticated numerical methods have been developed to find
(z, t) and thereby simulate sedimentation and fluidization of polydisperse suspensions of rigid spheres (Berres et al., 2003, 2004, 2005; Bürger et al., 2001). Particle-based simulations are an alternative to these (continuum) methods. Bargieł et al. (2005) distributed particles uniformly over a height, H, and simulated the sedimentation of polydisperse suspensions. To extend their method to

5576

M. Bargieł, E.M. Tory / Chemical Engineering Science 61 (2006) 5575 – 5589

fluidization, we assume that vk (
i ) is the velocity of the ith particle with respect to the suspension and that this velocity is determined by the concentration of a thin layer (h = 0.25 mm) of suspension immediately below (above) this particle if vk (
i ) is negative (positive) (Bargieł et al., 2005; Bargieł and Tory, 2006). Then, the particle-based equivalent of Eq. (1) is   dZki (t) = q + vk (
i ) dt,

(2)

where
i is the value of
that applies to the ith particle. (To avoid confusion, species are denoted by subscripts and individual particles by superscripts.) More realistically, the particles were initially distributed randomly over H in other simulations (Bargieł et al., 2005; Bargieł and Tory, 2006). In their work, the only parameter was the mean velocity k (
i ). Each particle is assumed to move independently with velocity vk (
i ) = k (
i ), but
i is determined by the number of particles of each species in the rate-determining region of thickness h. (To ensure that the rise of the packed bed is essentially linear, the value of h should be as small as possible, consistent with an efficient simulation.) The initial randomness implies that
i and hence vk (
i ) are different for every particle. Subsequent movements of many particles determine the future values of vk (
i ) for each particle. Therefore, in a set of simulations with the same initial concentration,
0 , each simulation is different. These features label this simulation as a (degenerate) Pickard–Tory process (Pickard and Tory, 1987). The variability of velocities of individual particles in gas fluidized beds has long been recognized (Harrison et al., 1961). More recently, that of slowly sedimenting particles has been accounted for by incorporating a term for hydrodynamic diffusion (Davis, 1996; Ham and Homsy, 1988) by analogy to molecular diffusion (Ham and Homsy, 1988). Then, sedimentation and fluidization can be studied as advection–diffusion processes (Tory and Ford, 2004) defined by    j
k j j
j
+
k (q + k (
)) = Dk (
) k , jt jz jz jz

(3)

where Dk (
) is the diffusion coefficient for the kth species and  j j
j j Dk (
) Dk (
) = . jz j
j jz K

(4)

j =1

Ramirez and Galvin (2005) use a simplified version of Eq. (3). Their Eq. (9) contains a typographical error. (The plus sign should be a minus.) If Dk (
) is constant, Eq. (3) simplifies to their corrected equation. However, experiments by Nicolai et al. (1995) show that D(
) varies substantially in monodisperse suspensions. We would expect even more variability in a polydisperse suspension. Given the complexity of Eq. (3), we could transform it to its stochastic equivalent (Laso, 1994) and attempt to solve it as a Pickard–Tory process (Bargieł and Tory, 1996; Pickard and Tory, 1987). Generalizing the treatment of

Tory and Ford (2004), we write     j j
k j j Dk (
) = −
k Dk (
) jz jz jz jz +

j2 [Dk (
)
k ]. jz2

(5)

When Eq. (5) is substituted in Eq. (3), the latter has the form of a Fokker–Planck equation (Laso, 1994) that is equivalent to the stochastic differential equation (SDE) (Laso, 1994)   j i i i dZk (t) = q + k (
) + Dk (
) dt jz + (2Dk (
i ))1/2 dW (t),

(6)

where W (t) is a one-dimensional standardized Wiener process (Kloeden and Platen, 1992). Simulations using Eq. (6) represent Pickard–Tory processes in which the parameters are k and Dk , both determined by the random variable
i . Eq. (6) reduces to Eq. (2) in the absence of diffusion. Eq. (6) has two advantages over Eq. (3). Though it still contains the derivative shown in Eq. (4), it is simpler than Eq. (3). Also, the continuum approximation is not nearly as good for suspensions as it is for gases and liquids. Small-scale experiments typically involve 105 –107 particles. Distributing particles randomly in a settling column yields irregularities in
that are clearly visible in dilute suspensions (Ham and Homsy, 1988; Tory and Pickard, 1977). One-dimensional simulations in which one or two million particles are randomly distributed over a height, H, reflect this variability in initial concentration (Bargieł et al., 2005; Bargieł and Tory, 2006). Thus, these simulations more closely resemble the actual experiments. However, the complexity of the derivative is still a difficulty. j
j /jz can be approximated by a central difference (Bargieł and Tory, 1996; Tory, 2000), but it is unlikely that the derivatives jDk (
i )/j
j , k =1, . . . , K, are known accurately. The second problem is that hydrodynamic diffusion is unlike molecular diffusion. The sudden changes in the velocities of molecules contrast with the gradual changes for particles in suspensions. According to Eq. (6), position is continuous, but velocity is not. In fact, particle velocities are continuous (Tory et al., 1994) in creeping flow. Indeed, the coefficient of autocorrelation decays rather slowly (Tory, 2000). A much better method is available. The theoretical basis for a particle-based simulation incorporating both the variability and continuity of velocity for monodisperse suspensions (Pickard and Tory, 1977, 1979; Tory and Pickard, 1977) was laid a decade before the hydrodynamic-diffusion approach. It was subsequently extended to polydisperse suspensions (Tory and Pickard, 1982). 2. Theory Identical particles in a suspension settle with different velocities and the velocity of each particle changes with time. These features are incorporated in a Markov model (Pickard and Tory, 1977). The Markov property reflects the well-known result that

M. Bargieł, E.M. Tory / Chemical Engineering Science 61 (2006) 5575 – 5589

the velocity of a particle in creeping flow depends only on the current configuration of the system and not on any previous configuration (Happel and Brenner, 1986; Kim and Karrila, 1991). Since a given concentration represents many different configurations (and hence many different velocities), the deterministic effect of this variability (of configuration on particle velocity) is replaced by a random process. Kynch’s theory (Bustos et al., 1999; Kynch, 1952) contains only a single parameter, the mean velocity. The spatial and temporal variability of identical particles at the same solid concentration shows that (at least) two additional parameters are necessary. Based on their own experiments and those of Koglin (1971a,b, 1972, 1973, 1976), Tory and Pickard (1977) formulated a three-parameter Markov model for the one-dimensional sedimentation of a monodisperse suspension. The parameters are the mean (
i ) and variance 2 (
i ) of the steady-state distribution of velocities, and the “relaxation time” or its reciprocal, (
i ). This Markov model, extended to include fluidization, can be expressed as   dZ i (t) = q + V i (t) dt,

(7)

  dV i (t) = − (
i ) V i (t) − (
i ) dt 1/2  (
i ) dW (t). + 2(
i )

(8)

In this model, steady-state velocities of particles are normally distributed and the autocorrelation of a particle’s velocity decays exponentially. Theoretical and experimental evidence confirms that these are excellent approximations (Tory, 2000) from
= 0.02 (Tory and Pickard, 1977) to
= 0.40 (Nicolai et al., 1995). The model also predicts that the times for identical spheres to traverse a fixed distance will closely approximate a lognormal distribution (Pickard et al., 1987; Tory et al., 1994), as found experimentally by Koglin (1971a,b). Particle velocities are continuous, but nowhere differentiable (Karlin and Taylor, 1981; Pickard et al., 1987). The particles have continuous sample paths (Karlin and Taylor, 1981; Pickard et al., 1987). All three parameters can be expressed as dimensionless variables. Thus, ∗ = a/|u∞ | = 1/t ∗ where t ∗ (the dimensionless relaxation time) is expressed in terms of t0 , the time for a particle moving with its Stokes velocity (defined by Eq. (17)) to traverse a distance equal to its characteristic dimension, in this case the sphere radius a. Tory (2000) and Nicolai et al. (1995) estimated t ∗ and the dimensionless parameters ∗ = /u∞ and ∗ = /|u∞ | for monodisperse suspensions. The velocity process, given by Eq. (8), is Markov, so the resulting position process [Z i (t)] is not. However, the vector process describing the joint values of position and velocity is Markov (Karlin and Taylor, 1981; Pickard et al., 1987). Position–velocity “skeletons” can then be simulated directly (Pickard et al., 1987). This is the key feature of our simulation, because it makes numerical integration unnecessary. The joint Markov transition structure for a monodisperse suspension is given in the Appendix.

5577

3. Sedimentation and fluidization of polydisperse suspensions This Markov model is easily generalized to K species. Then   dZki (t) = q + Vki (t) dt,

(9)

  dVki (t) = − k (
i ) Vki (t) − k (
i ) dt  1/2 + 2k (
i ) k (
i ) dW (t),

(10)

where Vki is the velocity of the ith particle with respect to the suspension. The values of the parameters k (
i ), k (
i ), and k (
i ) are determined by the value of
in the thin layer (of thickness h) that applies to the ith particle. If the last value of Vki was negative (downward) or zero, this layer is directly below the particle; if its velocity was positive (upward), directly above. The particles are initially distributed randomly over a height H. Initial velocities are distributed according to normal distributions with means k (
0 ) and variances 2k (
0 ). This yields the initial values, Zki (0) = z0i , Vki (0) = v0i , for Eqs. (11)–(15). Subsequent values of Zki and Vki are generated from a bivariate-normal distribution with means   E Zki (t + )|Zki (t) = zi , Vki (t) = v i     = zi + q + k (
i )  + 1 − exp(−k (
i ))   × v i − k (
i ) /k (
i ), (11)   E Vki (t + )|Zki (t) = zi , Vki (t) = v i   = k (
i ) + v i − k (
i ) exp(−k (
i )),

(12)

variances   2Z = Var Zki (t + )|Zki (t) = zi , Vki (t) = v i  = 2k (
i ) − 3 + 4 exp(−k (
i ))  − exp(−2k (
i )) 2k (
i )/2k (
i ),

(13)

  2V = Var Vki (t + )|Zki (t) = zi , Vki (t) = v i   = 2k (
i ) 1 − exp(−2k (
i )) ,

(14)

and covariance   ZV = Cov Zki (t + ), Vki (t + )|Zki (t) = zi , Vki (t) = v i 2  = 1 − exp(−k (
i )) 2k (
i )/k (
i ).

(15)

5578

M. Bargieł, E.M. Tory / Chemical Engineering Science 61 (2006) 5575 – 5589

In this Pickard–Tory process, each particle has its unique environment (the region of thickness h = 0.25 mm immediately below or above it) and its own Markov process. Its three parameters are random variables that are functions of the random variable
i . The value of the latter depends on the number of particles of each species that are currently in the rate-determining region which, in turn, depends on the Markov processes of other particles. The values of the parameters for all particles are determined at t and each particle is given a new position and velocity at t + . These are determined by two random numbers from a bivariate-normal distribution with means, variances and covariance given by Eqs. (11)–(15). After all particles have been moved, the values of
i at t +  are calculated. In this way, height–velocity skeletons, Zki (m), Vki (m), m0, for each particle are generated in parallel. In our previous simulations (Bargieł et al., 2005; Bargieł and Tory, 2006), the time increment  was small enough that all particles traveled less than h. Since the distance increment is now a random variable, we set  = 0.25 s so that (k + 4k ) < h, k = 1, . . . , K, in all our simulations. Thus, there is a negligible probability that any particle travels more than h in any time step, but some particles may change direction. To avoid the possibility that a particle’s intermediate position would hit the packed bed, but its final position would lie above it (Pickard et al., 1987; Tory et al., 1994), we assume that k (
i ) = 0 for all spheres within h of the packed bed. Velocities are much less variable in a steep concentration gradient (Bürger and Tory, 2000; Luke, 2000; Mucha et al., 2004; Tory et al., 1992), so this is a good approximation. Though only position–velocity skeletons are generated, Eq. (10) ensures that particle paths and velocities are continuous, steady-state velocities are normally distributed, and velocity autocorrelations decay exponentially. As this is a joint Markov process, there is no numerical error of integration (Pickard et al., 1987). Experimentally, of course,
i and hence the values of the three parameters are not constant during . However, this is not a serious problem. Eq. (10) is an Ito equation (Pickard and Tory, 1987). The use of non-anticipating functions and left-hand endpoints in the formal definition of the Ito integral guarantee that errors incurred by treating
i as fixed do not accumulate (Karlin and Taylor, 1981; Pickard and Tory, 1987). Thus, Eq. (10) can be simulated. Though
i is fixed in the time interval [t, t + ), its value at t is based on the values of
in [zi − h, zi ] if the particle is descending and [zi , zi + h] if the particle is ascending. This scheme was validated in our earlier studies (Bargieł et al., 2005; Bargieł and Tory, 2006). In our simulations, as in experiments, particles move in and out of the uppermost region (Tory and Pickard, 1977), so the topmost particle keeps changing. Though individual trajectories are of some interest (Tory et al., 1994), we do not require them here.

4. Parametric values The Masliyah–Lockett–Bassoon (MLB) equation (Bürger et al., 2002; Lockett and Bassoon, 1979; Masliyah, 1979) for

sedimentation of polydisperse suspensions yields the means k (
i ) = u∞1 (1 −
i )n−2 

i i k (k − (
i )) − K j =1 j
j (j − (
)) × , (1 − f ) (16) where u∞1 = −gd 21 (1 − f )/(18f )

(17)

is the Stokes velocity of the largest species, k is the density and dk the diameter of spheres of the kth largest species, f and f are the density and dynamic viscosity of the fluid, k = dk2 /d12 ,
ij is the concentration of the jth species in the rate-controlling region for the ith particle,
=
1 +
2 +· · ·+
K is the total solid fraction, and (
) = f (1 −
) + 1
1 + 2
2 + · · · + K
K is the density of the suspension. Note that u∞1 < 0 when 1 > f . The last term in Eq. (16) arises from the flow of fluid caused by the movement (up or down) of particles. The exponent, n, is not a universal constant. It varies from 4.0 to 6.55 and can be chosen to fit the experimental data (Bargieł and Tory, 2006). Bürger et al. (2002) provide a detailed derivation of Eq. (16). The MLB equation is a generalization of the Richardson–Zaki equation (Richardson and Zaki, 1954),   n 
i = u∞ 1 −
i , (18) and reduces to it, via Eq. (20), when k = 1 and k = s , k = 1, . . . , K. Eq. (16) can be written (Bürger et al., 2002) as ⎡ k (
i ) = u∞1 (1 −
i )n−2 ⎣k (k −
·
i )

−

K 

⎤ j
ij (j −
·
i )⎦

1 ,

(19)

j =1

where k =k −f and
=(1 , 2 , . . . , K ) is the (row) vector of density differences. The MLB equation has several attractive features. It is relatively simple and is based on slip velocities (Bürger et al., 2002). If n is used as a best-fit parameter, the equation agrees very closely with several experimental results (Bargieł et al., 2005; Bargieł and Tory, 2006). It correctly predicts that polydisperse suspensions of spheres of equal density are stable (Berres et al., 2003). Aside from random fluctuations and wall effects, the projections of particles in stable suspensions are uniformly distributed in any horizontal plane in the settling tank. Unstable suspensions exhibit three-dimensional structures such as “blobs” and vertical “columns” and thus cannot be represented by a one-dimensional equation (Bürger et al., 2002). For a particular bidisperse suspension of species with different densities (Bargieł and Tory, 2006; Bürger et al., 2002; Law et al., 1987) the MLB equation correctly predicts stable and unstable values of
.

M. Bargieł, E.M. Tory / Chemical Engineering Science 61 (2006) 5575 – 5589

When 1 = 2 = · · · = K , Eq. (19) reduces to ⎛ ⎞ K  j
ij ⎠ . k (
i ) = u∞1 (1 −
i )n−1 ⎝k −

bivariate-normal distribution (Larsen and Marx, 2001). Using the formulae (20)

j =1

For a bidisperse suspension, Eq. (20) yields the two equations  n−1  1 (
i ) = u∞1 1 −
i 1 −
i1 − 2
i2 , (21)  n−1  2 (
i ) = u∞2 1 −
i 1 −
i1 /2 −
i2 ,

(22)

Unfortunately, little is known about the other two parameters. The variance of particle velocities in monodisperse suspensions has long been a source of puzzlement and controversy (Bergougnoux et al., 2003; Caflisch and Luke, 1985; Guazzelli, 2001; Koglin, 1971a,b, 1972, 1973, 1976; Mucha et al., 2004; Nicolai et al., 1995; Tory et al., 1992). Very little is known about the variances in polydisperse suspensions. It is reasonable to assume (Tory and Ford, 2004) that the dimensionless velocity standard deviations and dimensionless relaxation rates, respectively, satisfy ∗k = ∗1 and ∗k = ∗1 in monodisperse regions in sedimenting suspensions and fluidized beds. Naively, this would suggest that 2k = k 21 and k = 1 dk /d1 for a given value of
in the polydisperse case. However, in polydisperse regions, there is computational (Höfler, 1997) and experimental (Peysson and Guazzelli, 1999) evidence that the motion of large spheres increases the variance of smaller spheres beyond its value in monodisperse suspensions of the same concentration. In dilute suspensions, small spheres close to a large sphere are dragged downward by it (Tory and Kamel, 1997) while small spheres further away are affected by the resulting upward flow of fluid. Using physical scaling arguments, Peysson and Guazzelli (1999) suggest k = 1, 2,

Z = (−2 ln r)1/2 cos(2s),

(26)

V = (−2 ln r)1/2 sin(2s),

(27)

this algorithm generates Z and V , which are N (0, 1), from two random numbers, r and s, that are uniformly distributed over the unit interval. Since  1 cos(2s) sin(2s) ds = 0, (28) 0

since 1 = 1 by definition. For very small values of
, these have the same form as Batchelor’s (1982) equation, which in our terminology is  (23) k (
i ) = u∞k 1 + Sk1
i1 + Sk2
i2 , k = 1, 2.

k ≈
1/3 |u∞V |(dk /dV )1/2 ,

5579

(24)

where u∞V is the Stokes velocity of a virtual sphere of diameter  1/3 dV = (d13 + 1/2 d23 )/(1 + )1/2 , (25) = n2 /n1 , and nk is the number of spheres of species k in the region. Eq. (24) was in good agreement with their experimental values. Their
1/3 dependence contrasts with
1/2 derived theoretically (Caflisch and Luke, 1985; Tory et al., 1992). They did not measure values of . Ways of measuring  will be discussed later in this paper. 5. Algorithms We adapted the Box–Müller algorithm (Kloeden and Platen, 1992; Tory et al., 1994) to generate random numbers from a

Z and V are uncorrelated. Noting also that   1 2 1 1 2 cos (2s)2 ds = cos2 d = 1,  0  0 we see that  1 1 cos(2s) sin(2s + ZV )2 ds = sin ZV ,  0

(29)

(30)

where sin ZV = ZV /Z V = Cov( Z , V )

(31)

is the correlation coefficient (Larsen and Marx, 2001). Thus, Eq. (27) is replaced by

V = (−2 ln r)1/2 sin(2s + ZV ),

(32)

where ZV = Arcsin(ZV /Z V ).

(33)

Let EZ = E[Zki (t + )|Zki (t) = zi , Vki (t) = v i ]

(34)

and EV = E[Vki (t + )|Zki (t) = zi , Vki (t) = v i ].

(35)

Then zi (t + ) = EZ + Z Z

(36)

and v i (t + ) = EV + V V .

(37)

Using Eqs. (26) and (32), it is straightforward to verify that Z and V are N (0, 1) and that Eqs. (36)–(37) satisfy Eqs. (11)–(15). The algorithm for determining the solid concentration that applies to an individual particle is given by Bargieł and Tory (2006). 6. Simulations Random motion of particles occurs at all concentrations, but global effects are much greater when
is very small. Fortunately, it is then much easier to study the velocities and trajectories of individual particles. Thus, our examples utilize data from experiments at very low concentrations.

5580

M. Bargieł, E.M. Tory / Chemical Engineering Science 61 (2006) 5575 – 5589

6.1. Example 1. Sedimentation of a bidisperse suspension

500

Peysson and Guazzelli (1999) studied fluctuations in horizontal and vertical velocity components in a dilute bidisperse suspension of glass spheres. Only the vertical velocities are relevant to our simulations. In our notation, their mean velocities are height [mm]

2 (
) = u∞2 (1 + S21
1 + S22
2 ) = u∞2 (1 − 4
1 − 5.6
2 ).

(38)

400

300

200

(39)

(In their notation, the smaller species was designated as 1, whereas the opposite is true in our case.) Owing to the relatively small effect of
1 and
2 on these means, there is considerable uncertainty (Peysson and Guazzelli, 1999) about the coefficients in these equations. Peysson and Guazzelli found 2 /|2 | ≈ 0.5 for a monodisperse suspension with
= 0.01. This compares with 2 /|2 | ≈ 0.61 calculated by Bürger and Tory (2000) for the monodisperse suspension studied by Tory and Pickard (1977) in which
= 0.02. Setting = 0 in Eq. (25) and calculating 1 from Eq. (24), we see that their results are consistent with those of Bürger and Tory, in the sense that the different values of 2 /|2 | are in line with the
1/3 law in Eq. (24). Peysson and Guazzelli plotted 2 (
)/2 (
) and found that 2 /|2 | ≈ 1.02 at
1 =
2 = 0.01. As noted earlier, this higher value is expected. Though agreement in the values of  does not necessarily mean that the values of  will be the same, we use the value of  calculated by Bürger and Tory, yielding the dimensionless relaxation rate ∗ (0.02) ≈ 0.015 (Tory, 2000). In the concentration gradient at the top of their suspension, ∗ ≈ 0.035 (Tory, 2000). In creeping flow, the time for a single sphere to attain its Stokes velocity is much less than that for a similar change in the velocity of a sphere in a suspension. This implies that ∗ is large in very dilute suspensions. This is consistent with its larger value in the gradient (in which the concentration is lower). Choosing the form ∗ = C1 /(
2 + C2 ) and inserting the values ∗ (0.01) = 0.035 and ∗ (0.02) = 0.015, we obtain  ∗ = 7.875 × 10−6 /
2 + 1.25 × 10−4 . (40) Both the suspension–fluid interface at the top and the bi–monodisperse interface are diffuse, but the decreasing value of ∗ and the rapidly increasing value of ∗ with decreasing
ensure that discernible interfaces will form. The quantitative values are very uncertain, but the qualitative trends are correct. We simulated the sedimentation of the particles studied by Peysson and Guazzelli. As our simulation requires means for the entire range of concentrations, we used Eq. (20) with n=4.4, which closely approximates Eqs. (38) and (39) for very small values of
1 and
2 . In our notation, a1 =394
m, a2 =200
m, 1 = 2 = 2.50 g/cm3 , f = 1.09 g/cm3 , f = 13 Poise,
2 = 0.01, 0

1
0.02. Thus, u∞1 ≈ −0.0367 cm/s and u∞2 ≈ −0.00945 cm/s. Values of k were calculated from Eq. (24)

100

0

0

2000

4000

6000

time [s] Fig. 1. Height–time curves for the uppermost sphere of each species. The curves for the Markov model lie slightly above those calculated from the mean value. The solid lines near the bottom show the growth of the packed bed.

500

height [mm]

1 (
) = u∞1 (1 + S11
1 + S12
2 ) = u∞1 (1 − 5.6
1 − 4
2 ),

Example 1 Markov model mean velocity

490

Example 1 Markov model mean velocity 480

0

10

20

30

40

50

time [s] Fig. 2. Detail of curves for the first 50 s of sedimentation.

and values of k from Eq. (40) and k = ∗ u∞k /ak .

(41)

According to information provided in their paper, the suspension had an initial height H = 50 cm, and sedimented over an area A = 16 cm2 , giving an initial volume V = 800 cm3 . Let Vk = 4ak3 /3 be the volume of a sphere of the kth species. Then V1 ≈ 2.562 × 10−4 cm3 and V2 ≈ 3.351 × 10−5 cm3 . We used

M. Bargieł, E.M. Tory / Chemical Engineering Science 61 (2006) 5575 – 5589

500

500

500 100 s

100 s

490

490

480

480

480

470

460

height [mm]

490

height [mm]

height [mm]

100 s

450 0.0

470

460

0.01

0.02

5581

450 0.0

0.03

470

460

0.01

0.02

450 0.0

0.03

0.01





1

0.02

0.03




2

Fig. 3. Solid profiles at an early stage of sedimentation.

420

330

4000 s

320

140

400

310

130

390

370 0.0

height [mm]

410

380

(a)

150 2000 s

height [mm]

height [mm]

1000 s

300

290

0.01
2

280 0.0

0.02 (b)

120

110

0.01
2

100 0.0

0.02 (c)

0.01
2

0.02

Fig. 4. Solid profile of the smaller species at 1000, 2000, and 4000 s.

31,226 and 238,735 spheres of species 1 and 2, respectively, to obtain
1 =
2 = 0.01. Fig. 1 shows that fluctuations in velocity in the Markov model cause the uppermost small sphere (top line) to lie slightly above the interface calculated solely on the basis of the mean velocity (i.e., 1 = 2 = 0). Similarly, the lower dashed line shows that the position of the uppermost large sphere gradually diverges from the corresponding interface. For the small spheres, most of the divergence occurs early, as shown in Fig. 2, which is an enlargement of the top left corner of Fig. 1. The uppermost curve shows the “induction period”, which is well known in experimental work. This is the period, approximately 10 s in this simulation, in which the faster settling particles leave the upper region. Since the inflow of particles from below is insuf-

ficient to balance this outflow, a concentration gradient forms and gradually expands. In Fig. 2, the distance between the two lower curves suggests the formation of a gradient in the concentration of the larger species, as there is a difference between the positions of particles initially at the top in the deterministic case (where the concentration is still at roughly its initial value) and the topmost large particle (where
1 vanishes). Solid profiles confirm that gradients exist for both species. Profiles at 5, 10, 15, and 20 s (not shown) explain the unusual shape of the path of the topmost large sphere. At 5 s, there is already a gradient in the concentration of small spheres at the top of the suspension, so all the larger spheres quickly settle out of this region. By 10 s, all of the larger spheres are below this gradient, so they settle more slowly. As the gradient in the concentration

5582

M. Bargieł, E.M. Tory / Chemical Engineering Science 61 (2006) 5575 – 5589

190

360

500

500 s 180

480

340

170

470

450 0.0

height [mm]

350

330

0.01
1

310 0.0

0.02 (b)

160

150

320

460

(a)

1000 s

490

height [mm]

height [mm]

0s

0.01
1

140 0.0

0.02 (c)

0.01


0.02

1

Fig. 5. Solid profile for the larger species at 0, 500, and 1000 s.

of the larger spheres expands, the boundary represented by the uppermost large sphere descends more quickly (as would also be predicted by an advection–diffusion equation), producing the inflection point in the second dashed line. (As noted elsewhere, the sphere that is uppermost may change as settling proceeds.) In Figs. 3–5, which show the gradients for both species, each point represents the average concentration in a height interval of 0.25 mm (the value of h). Fig. 3 shows the solid profiles for both species at an early stage of sedimentation. Fig. 4 illustrates the evolution of the solid profile of the smaller species in the monodisperse region. The gradient gradually expands from about 9 mm at 100 s to 14 mm at 500 s (not shown), 17 mm at 1000 s, 22 mm at 2000 s, 28 mm at 4000 s, and 30 mm at 5000 s (not shown). The diffusion near the interface is counteracted to some extent by “self-sharpening” caused by the increase in |2 (
2 )| as
2 → 0. Although the trend is clear, the scatter means that there is considerable uncertainty about the quantitative values. Fig. 5 shows the solid profiles of the larger spheres at 500 and 1000 s. The scatter in the initial values, also shown in Fig. 5, emphasizes the difference between smallscale experiments with tens or hundreds of thousands of particles and the idealized smoothing of the continuum approach. The spreading of the interface between the bi- and monodisperse suspensions is evident from Figs. 3 and 5. The gradient expands from roughly 25 mm at 100 s to 40 mm at 500 s, and to 50 mm at 1000 s. Here also, the diffusion is counteracted to some extent by self-sharpening. The exact values are very uncertain, but the spread is faster for the large spheres because 1 (
) > 2 (
2 ). The shapes of the solid profiles are also different because 2 (
2 ) → 0 and ∗2 (
2 ) → 6.3 × 10−2 as
2 → 0, whereas 1 (
) remains finite and ∗1 (
) = ∗1 (
) increases more slowly as
→ (0,
2 )T . Though there have been many experimental studies of settling of polydisperse systems, most of them measured only the rates of fall of the various interfaces. However, experimental

studies of closely sized spheres show that velocity fluctuations play an important role in the formation of concentration gradients at the solid–fluid interface (Davis and Hassen, 1988; Lee et al., 1992). Lee et al. (1992) measured interface broadening of a suspension of rigid glass microspheres (average radius 67.9
m, standard deviation 4.0
m) in suspensions with 0.05


0.15. Their Figure 2 shows that the interface expands while moving downward as settling proceeds. As in our simulation, broadening is promoted by velocity fluctuations and opposed by self-sharpening. The slight polydispersivity also promotes broadening. 6.2. Example 2. Fluidization of a bidisperse suspension The second example is the fluidization of a hypothetical bidisperse suspension with a1 = 394
m (as before) and a2 = 385
m. Thus, 2 ≈ 0.9548. We used 489,174 particles of species 1 and 524,288 particles of species 2. As in Example 1, V1 ≈ 2.562 × 10−4 cm3 , so we set V = 12, 532.64 cm3 so that
1 =
2 = 0.01 initially. Also we set H = 100 cm, so A=125.3264 cm2 . (Since this is a one-dimensional simulation, the value of A is meaningless, but it is required to obtain the desired value of H.) As in Example 1, 1 = 2 = 2.50 g/cm3 , f = 1.09 g/cm3 , f = 13 Poise, and u∞1 ≈ −0.0367 cm/s. As these spheres are very nearly the same size, we used 1 = 2 =
1/3 |u∞1 |(1 + 2 )/2. We used Eq. (20) to calculate settling velocities, Eqs. (40) and (41) to obtain k and fluidized the suspension by setting q = 0.0336 cm/s. In fluidization (in contrast to sedimentation), no particles settle out. However, some particles may hit the bottom; these are reflected with the same value of |Vki |. It is supposed that, for this bidisperse system, the fluidized bed may eventually separate into two regions: an upper region consisting of smaller particles and a lower region of larger particles. Assuming that Eq. (18) holds for monodisperse regions, the theoretical

M. Bargieł, E.M. Tory / Chemical Engineering Science 61 (2006) 5575 – 5589

equilibrium concentrations are given by q + u∞1 (1 −
1 )n = 0.0336 − 0.0367(1 −
1 )4.4 = 0 ⇒
1 ≈ 0.0199

velocities of spheres of both species, so both of these positions rise quickly, but dilution from bed expansion causes both species to settle more rapidly, and the slopes of both curves decrease. The two curves diverge very slowly until about 40,000 s and then separate rapidly. However, this is somewhat misleading. Figs. 7–10 show the evolution of concentration gradients. In these figures only, each point represents the mean concentration in a height of 2 mm (8h). The solid profiles in Fig. 7 show that, at 10,000 s, the larger species is greatly depleted near the top of the suspension and for a considerable distance downward. As this depletion occurs, the smaller particles settle more rapidly because
and (especially)
1 have decreased. Particles that were left behind earlier get closer to those that were left behind later. Consequently, the concentration of the smaller particles increases. Bargieł et al. (2005) give a simple graphical proof of this “Smith effect” (Smith, 1966). In Fig. 7, the concentration of the smaller species has increased considerably beyond its initial value of 0.01. This effect is moderated by the expansion of the bed and the upward diffusion of particles near the top. In marked contrast to the gradient in the concentration of the larger species, that of the smaller extends over a much shorter distance at the top. Fig. 8 shows the profiles at 40,000 s. Particles of species 2 are largely absent from the lowest region, so the concentration of species 1 there has reached its steady-state value. The gradient for species 1 extends far into the upper levels, but the upper part contains relatively few particles of this species. Similarly, the gradient for species 2 extends downward into the lower region, but the number of small particles there is not large. By 80,000 s, the species have largely separated, as shown in Fig. 9. By 150,000 s, the system is clearly at steady state, as indicated by the leveling out of the heights of the uppermost particles in Fig. 6. Fig. 10 shows the final profiles and indicates that the theoretical values of
1 and
2 are obtained. The two gradients are sharper than those in Fig. 9,

(42)

and 0.0336 − 0.0367(0.9548)(1 −
2 )4.4 = 0 ⇒
2 ≈ 0.0095.

(43)

Fig. 6 traces the positions of the top sphere of each species. Initially, the fluidization velocity exceeds the settling

1600

height [mm]

1400

1200

1000

800

Example 2 Species 1 Species 2

600 0

50000

100000

5583

150000

time [s] Fig. 6. Height–time curves for the uppermost sphere of each species in a fluidized bed.

1000

1000

800

800

800

600

height [mm]

1000

height [mm]

height [mm]

10000 s

600

600

400

400

400

200

200

200

0 0.0

0.01

0.02
1

0.03

0 0.0

0.01

0.02
2

0.03

0 0.0

0.01

0.02

0.03




Fig. 7. Solid profiles at an early stage of fluidization. The initial concentrations were
1 =
2 = 0.01.

5584

M. Bargieł, E.M. Tory / Chemical Engineering Science 61 (2006) 5575 – 5589

1400

1200

1200

1200

1000

1000

1000

800 600

height [mm]

1400

height [mm]

height [mm]

40000 s 1400

800 600

800 600

400

400

400

200

200

200

0 0.0

0.01

0.02

0 0.0

0.03

0.01




0.02

0 0.0

0.03

0.01


2

1

0.02

0.03




Fig. 8. Solid profiles at an intermediate stage of fluidization.

1600

1400

1400

1400

1200

1200

1200

1000

1000

1000

800 600

height [mm]

1600

height [mm]

height [mm]

80000 s 1600

800 600

800 600

400

400

400

200

200

200

0 0.0

0.01

0.02

0.03

0 0.0


1

0.01

0.02

0.03

0 0.0


2

0.01

0.02

0.03




Fig. 9. Solid profiles at an advanced stage of fluidization.

but the variability of particle velocities prevents the complete separation of species with almost the same Stokes velocities. Ramirez and Galvin (2005) solved a related fluidization problem. Their particles had the same size and differed only slightly in density. Using a stiff differential-equation routine, they obtained solid profiles that clearly showed the effects of dispersion on the heavy–light boundary. Their bed height (corresponding to our uppermost sphere) showed a rise and fall somewhat similar to that of our topmost large sphere. Their initial and final concentrations were much greater than ours, so the two problems are not directly comparable. The experimental profiles for this problem are more similar to ours in that they clearly show the effects of dispersion at the solid–fluid interface. This is captured in earlier computed profiles, but not in the final.

They note that a common weakness of previous models has been the need to know a priori the final conditions. They stress that their approach does not require this knowledge and hence is truly predictive, assuming that the dispersion coefficient is known. Our simulation is also truly predictive; we impose a fluidization velocity and the system eventually finds its steady state. 6.3. Example 3. Buoyant monodisperse suspension below a clear fluid The third example is related to the situation in which a dilute suspension of particles lies above a clear fluid or a more dilute suspension. We simulate the equivalent case in which

M. Bargieł, E.M. Tory / Chemical Engineering Science 61 (2006) 5575 – 5589

5585

1600

1400

1400

1400

1200

1200

1200

1000

1000

1000

800

height [mm]

1600

height [mm]

height [mm]

150000 s 1600

800

800

600

600

600

400

400

400

200

200

200

0 0.0

0.01

0.02

0 0.0

0.03

0.01


1

0.02

0.03

0 0.0

0.01


2

0.02

0.03




Fig. 10. Solid profiles of fluidized bed at steady state.

1080

2400 monodisperse 500 s

monodisperse 10000 s

1060

2200 height [mm]

2300

height [mm]

1070

1050

2100

1040

2000

1030 0.0 (a)

0.01

0.02

1900 0.0

0.03




(b)

0.01

0.02

0.03




Fig. 11. Diffuse interface in a monodisperse suspension of buoyant particles ascending and diffusing upward: (a) early stage and (b) advanced stage.

a monodisperse suspension of buoyant particles lies below a clear fluid. Here, d = 1.171 mm, u∞ = 0.1276 mm/s (recall that upward velocities are positive), and
0 = 0.02. The simulation involves 570,584 spheres with a total volume of

479.92 cm3 in a very tall container with a cross-sectional area of 239.86 cm2 . Initially, the suspension is randomly distributed over the lowest 100 cm. The mean velocity of the spheres is given by Eq. (18) with n = 4.7. This value provides a good fit

5586

M. Bargieł, E.M. Tory / Chemical Engineering Science 61 (2006) 5575 – 5589

Table 1 Characteristics of polydisperse suspension Species

Number of spheres

Diameter (mm)

1 2 3 4 5 6 7 8 9 10

40,000 80,000 170,000 140,000 50,000 20,000 20,000 20,000 20,000 10,000

1.20 1.19 1.18 1.17 1.16 1.15 1.14 1.13 1.12 1.11

570,000


k × 103 1.509 2.943 6.097 4.895 1.704 0.664 0.647 0.630 0.613 0.299

k 1.0000 0.9834 0.9669 0.9506 0.9344 0.9184 0.9025 0.8867 0.8711 0.8566

20.001

for many monodisperse suspensions (Scott, 1984). The standard deviation of the velocity fluctuations is  =
1/3 u∞ , and  = 2∗ u∞ /d = 0.2179∗ /s, where ∗ is given by Eq. (40). Normally, in sedimenting suspensions or fluidized beds, hydrodynamic diffusion broadens the interface while the concentration dependence of the mean velocity (given by Eq. (18) in this case) tends to sharpen it. In this simulation, diffusion broadens the interface as usual, but Eq. (18) broadens it even more because the dilute region moves upward more rapidly. The solid profile in Fig. 11(a) shows that a well-defined gradient exists at 500 s. Fig. 11(b) shows that the gradient extends over almost 10 times the height by 10,000 s. 6.4. Example 4. Buoyant polydisperse suspension below a clear fluid We compare Example 3 with the upward sedimentation of a related polydisperse suspension with a very narrow size distribution. This suspension consists of 10 species with the properties described in Table 1. From columns 2 and 3, the arithmetic mean diameter is 1.171 mm. (More sophisticated means are available, Shannon et al., 1963; Scott and Mandersloot, 1979, but are unsuited to our purpose.) The Stokes velocity of a particle with this mean diameter is the same as in Example 3 (0.1276 mm/s). Velocities are given by Eq. (20) with u∞1 = (1.200/1.171)2 (0.1276 mm/s) ≈ 0.1340 mm/s. As in the previous example,
0 = 0.02. Since these spheres are very closely sized, for all of them we use  =
1/3 u∞ =
1/3 (0.1276) mm/s and  = 0.2179∗ /s with ∗ given by Eq. (40). (The size distribution is adapted from a random sample of 55 spheres measured by Verhoeven, 1963. We have reduced the size of the particles to use more of them, smoothed the distribution slightly by adding 20,000 spheres of species 7, and made the spheres buoyant.) Figs. 12(a) and (b), depicting the overall solid profiles of the polydisperse suspension, curve more strongly upward at the top (compared to the monodisperse case), reflecting the presence of the three largest species. Otherwise, the profiles are very similar to those in the monodisperse case. Fig. 13 shows the positions of the top spheres of species 1 (largest), 4

(average), and 10 (smallest) compared to the position of the top sphere in the monodisperse suspension. Note that the values for species 4 essentially coincide with those of the monodisperse suspension. As expected, there is some differential rising by size, even though the differences are small. The asymmetry of the distribution is reflected in the positions of the top spheres of the largest and smallest species compared to the average. 7. Discussion The Markov model for sedimentation was designed to simulate (Pickard and Tory, 1974) the behavior of suspensions via a Pickard–Tory process (Pickard and Tory, 1974, 1987), but early work necessarily focused on proofs (Pickard and Tory, 1977, 1979; Tory and Pickard, 1982) and comparisons of theoretical and experimental distributions of velocity and transit times (Pickard et al., 1987; Tory and Pickard, 1977; Tory et al., 1994; Tory and Hesse, 1996; Tory, 2000). Subsequently, Moore’s law (Voller and Porté-Agel, 2002) and efficient algorithms (Bargieł et al., 2005; Bargieł and Tory, 2006) made simulation a practical, flexible alternative to numerical methods applied to continuum models (Berres et al., 2004; Bürger et al., 2001). The parameters in the Markov model can be measured directly by following individual particles in the suspension experimentally (Nicolai et al., 1995) or computationally (Ladd, 1997). The diffusion coefficient in the advection–diffusion equation, Eq. (3), is derived from more fundamental variables (Nicolai et al., 1995; Tory, 2000) and represents a loss of information. In this equation, the complexity of the diffusion term increases rapidly with the number of species. In contrast, the Markov model remains relatively simple, depending only on the values of k (
), k (
), and k (
). Also, as noted earlier, the Markov model preserves more of the true nature of sedimentation by ensuring that particle velocities are continuous, the steady-state distribution of velocities is Gaussian, and the autocorrelation of velocity decays exponentially. (See Eq. (A.5) for the monodisperse case.) The simulation of fluidization (Example 2) shows that fluctuations in concentration do not grow with time. The best-fit values are very tightly determined throughout the experiment. Very few experimental values of k (
) and k (
) are known for polydisperse systems. However, it is possible to postulate a range of plausible values and then carry out simulations to determine whether these parameters have any significant effect on the global behavior of the suspension. Clearly, they have a great effect on individual particles, but the cumulative effect on the global behavior may be small (Pickard and Tory, 1979; Tory and Pickard, 1982). Thus, simulations can play a role in determining when it is appropriate to invest in experiments to determine the values of these additional parameters. Though we have used the MLB equation in this work and in our earlier papers (Bargieł et al., 2005; Bargieł and Tory, 2006), any k (
) can be used. The method can also be extended beyond creeping flow provided that the future velocity of a particle depends (at least primarily) on its current velocity (and not some previous velocity). We hope to explore such subjects as layer inversion (Moritomi et al., 1982),

M. Bargieł, E.M. Tory / Chemical Engineering Science 61 (2006) 5575 – 5589

5587

2400

1080 polydisperse 500 s

polydisperse 10000 s

1060

2200 height [mm]

2300

height [mm]

1070

1050

2100

1040

2000

1030 0.0

0.01

0.02


(a)

1900 0.0

0.03

0.01

0.02

0.03




(b)

Fig. 12. Diffuse interface in a polydisperse suspension of buoyant particles ascending and diffusing upward: (a) early stage and (b) advanced stage.

perse suspension with 1 , 2 , and 3 all different (Bargieł and Tory, 2006). Given that k (
) smoothes concentrations while k (
) smoothes changes in velocities, a stochastic simulation of sedimentation or fluidization of such suspensions should be feasible.

3000

height [mm]

Species 1 Species 4 Species 10 Monodisperse

Notation 2000

1000 0

5000

10000

15000

time [s] Fig. 13. Comparison of positions of uppermost spheres in monodisperse and polydisperse suspensions.

complete mixing (Moritomi et al., 1986), and multi-species segregation and dispersion in liquid-fluidized beds (Ramirez and Galvin, 2005) and bed classifiers (Chen et al., 2002). Using just k (
), we have already simulated the sedimentation of a tridis-

dk Dk f fk h H K q t u∞ u∞1 vk z

diameter of spheres of kth species diffusion coefficient for the kth species vector of solid flux densities flux density of kth species thickness of velocity-determining layer initial height of suspension number of species fluidization velocity time Stokes velocity Stokes velocity of the largest species (or densest species if all diameters are equal) velocity of particles of kth species height

Greek letters k k

reciprocal of the relaxation time for kth species dk2 /d12

5588

k f k k
2k 

k


M. Bargieł, E.M. Tory / Chemical Engineering Science 61 (2006) 5575 – 5589

mean velocity of kth species density of fluid density of kth species k −  f row vector of densities (1 , 2 , . . . , K ) variance of velocity of kth species time step total solid concentration of sedimenting particles concentration of kth species vector of solid concentrations (
1 ,
2 , . . . ,
K )T

Acknowledgments We thank the reviewers for many helpful comments. Appendix A Our treatment follows that of Pickard et al. (1987) with a slight change in notation that simplifies the equations. In the three-parameter model, velocity is given by the SDE dV = −(V − ) dt + (2)1/2  dW .

(A.1)

Velocities are drawn to the mean  in proportion to their distance from , and random noise is superimposed. The resulting velocity process is an ergodic Markov diffusion. Its transition structure (over a Markov time step  and starting with velocity v) V (t + )|(V (t) = v)

(A.2)

is Gaussian with mean m(v, ) =  + (v − ) exp(−),

(A.3)

and variance

(v, ) = 2 (1 − exp(−2)).

(A.4)

From Eqs. (A.3) and (A.4), the steady-state distribution of velocities is normal with mean  and variance 2 . The autocorrelation of velocities is r() = Cov[V (t + ), V (t)]/2 = exp(−).

(A.5)

The joint position–velocity process is both Markov and Gaussian. The transition structure Z(t + ), V (t + )|Z(t) = z, V (t) = v

(A.6)

is given by the means E[Z(t + )|Z(t) = z, V (t) = v] = z +  + [1 − exp(−)](v − )/,

(A.7)

E[V (t + )|Z(t) = z, V (t) = v] = m(v, ),

(A.8)

variances Var[Z(t + )|Z(t) = z, V (t) = v] = [2 − 3 + 4 exp(−) − exp(−2)]2 /2 , Var[V (t + )|Z(t) = z, V (t) = v] = (v, ),

(A.9) (A.10)

and covariance Cov[Z(t + ), V (t + )|Z(t) = z, V (t) = v] = [1 − exp(−)]2 2 /.

(A.11)

References Bargieł, M., Tory, E.M., 1996. Stochastic dynamic solution of nonlinear differential equations for transport phenomena. A.I.Ch.E. Journal 42, 889–891. Bargieł, M., Tory, E.M., 2006. An extension of the particle-based approach to simulating the sedimentation of polydisperse suspensions. International Journal of Mineral Processing, in press. Bargieł, M., Ford, R.A., Tory, E.M., 2005. Simulation of sedimentation of polydisperse suspensions: a particle-based approach. A.I.Ch.E. Journal 51, 2457–2468. Batchelor, G.K., 1982. Sedimentation in a dilute polydisperse system of interacting spheres. Part 1. General theory. Journal of Fluid Mechanics 119, 379–408. Bergougnoux, L., Ghicini, S., Guazzelli, É., Hinch, J., 2003. Spreading fronts and fluctuations in sedimentation. Physics of Fluids 15, 1875–1887. Berres, S., Bürger, R., Karlsen, K.H., Tory, E.M., 2003. Strongly degenerate parabolic–hyperbolic systems modeling polydisperse sedimentation with compression. SIAM Journal of Applied Mathematics 64, 41–80. Berres, S., Bürger, R., Karlsen, K.H., 2004. Central schemes and systems of conservation laws with discontinuous coefficients modeling gravity separation of polydisperse suspensions. Journal of Computational and Applied Mathematics 164–165, 53–80. Berres, S., Bürger, R., Tory, E.M., 2005. On mathematical models and numerical simulation of the fluidization of polydisperse suspensions. Applied Mathematical Modelling 29, 159–193. Bürger, R., Tory, E.M., 2000. On upper rarefaction waves in batch settling. Powder Technology 108, 74–87. Bürger, R., Fjelde, K.K., Höfler, K., Karlsen, K.H., 2001. Central difference solutions of the kinetic model of settling of polydisperse suspensions and three-dimensional particle-scale simulations. Journal of Engineering Mathematics 41, 167–187. Bürger, R., Karlsen, K.H., Tory, E.M., Wendland, W.L., 2002. Model equations and instability regions for the sedimentation of equal spheres. Zeitschrift für Angewandte Mathematik und Mechanik 82, 699–722. Bustos, M.C., Concha, F., Bürger, R., Tory, E.M., 1999. Sedimentation and Thickening. Phenomenological Foundation and Mathematical Theory. Kluwer Academic Publishers, Dordrecht. Caflisch, R.E., Luke, J.H.C., 1985. Variance in the sedimentation speed of a suspension. Physics of Fluids 28, 759–760. Chen, A., Grace, J.R., Epstein, N., Lim, C.J., 2002. Steady state dispersion of mono-size, binary and multi-size particles in a liquid fluidized bed classifier. Chemical Engineering Science 57, 991–1002. Davis, R.H., 1996. Hydrodynamic diffusion of suspended particles: a symposium. Journal of Fluid Mechanics 310, 325–335. Davis, R.H., Hassen, M.A., 1988. Spreading of the interface at the top of a slightly polydisperse sedimenting suspension. Journal of Fluid Mechanics 196, 107–134. Drew, D.A., Passman, S., 1999. Theory of Multicomponent Fluids. Springer, New York. Guazzelli, É., 2001. Evolution of particle–velocity correlations in sedimentation. Physics of Fluids 13, 1537–1540. Ham, J.M., Homsy, G.M., 1988. Hindered settling and hydrodynamic dispersion in quiescent sedimenting suspensions. International Journal of Multiphase Flow 14, 533–546.

M. Bargieł, E.M. Tory / Chemical Engineering Science 61 (2006) 5575 – 5589 Happel, J., Brenner, H., 1986. Low Reynolds Number Hydrodynamics. Martinus Nijhoff, Dordrecht. Harrison, D., Davidson, J.F., de Kock, J.W., 1961. On the nature of aggregative and particulate fluidization. Transactions of the Institution of Chemical Engineers 39, 202–211. Höfler, K., 1997. Simulation and modeling of mono and bidisperse suspensions. Ph.D. Dissertation, Universität Stuttgart. Karlin, S., Taylor, H.M., 1981. A Second Course in Stochastic Processes. Academic Press, San Diego. Kim, S., Karrila, S.J., 1991. Microhydrodynamics, Principles and Selected Applications. Butterworth-Heineman, Boston. Kloeden, P.E., Platen, E., 1992. Numerical Solution of Stochastic Differential Equations. Springer, Berlin. Koglin, B., 1971a. Untersuchungen zur Sedimentationsgeschwindigkeit in neidrig konzentrierten Suspensionen. Ph.D. Dissertation, Universität Karlsruhe. Koglin, B., 1971b. Statistische Verteilung der Sedimentationsgeschwindigkeit in neidrig konzentrierten Suspensionen. Chemie Ingenieur Technik 43, 761–764. Koglin, B., 1972. Settling rate of individual particles in suspension. In: Groves, M.J., Wyatt Sargent, J.L. (Eds.), Particle Size Analysis. Society for Analytical Chemistry, London, pp. 223–235. Koglin, B., 1973. Dynamic equilibrium of settling velocity distribution in dilute suspensions of spherical and irregularly shaped particles. Proceedings of the Conference on Particle Technology, Chicago, pp. 265–271. Koglin, B., 1976. Zum mechanismus der Sinkgeschwindigkeitserhöhung in niedrig konzentrierten Suspensionen. In: 27 Vorträge des 1. Europäischen Symposions “Partikelmesstechnik”, Nürnberg, 17.–19 September 1975, Dechema Monographien, vols. 1589–1615, 79B, pp. 235–250. Kynch, G.J., 1952. A theory of sedimentation. Transactions of Faraday Society 48, 166–176. Ladd, A.J.C., 1997. Sedimentation of homogeneous suspensions of nonBrownian spheres. Physics of Fluids 9, 491–499. Larsen, R.J., Marx, M.L., 2001. An Introduction to Mathematical Statistics and its Applications. third ed. Prentice-Hall, Upper Saddle River, NJ. Laso, M., 1994. Stochastic dynamic approach to transport phenomena. A.I.Ch.E. Journal 40, 1297–1311. Law, H.-S., Masliyah, J.H., MacTaggart, R.S., Nandakumar, K., 1987. Gravity separation of bidisperse suspensions, light and heavy particle species. Chemical Engineering Science 42, 1527–1538. Lee, S., Jang, Y., Choi, C., Lee, T., 1992. Combined effect of sedimentation velocity fluctuation and self-sharpening on interface broadening. Physics of Fluids A 4, 2601–2606. Lockett, M.J., Bassoon, K.S., 1979. Sedimentation of binary particle mixtures. Powder Technology 24, 1–7. Luke, J.H.C., 2000. Decay of velocity fluctuations in a stably stratified suspension. Physics of Fluids 12, 1619–1621. Masliyah, J.H., 1979. Hindered settling in a multi-species particle system. Chemical Engineering Science 34, 1166–1168. Moritomi, H., Iwase, T., Chiba, T., 1982. A comprehensive interpretation of solid layer inversion in liquid fluidized beds. Chemical Engineering Science 37, 1751–1757. Moritomi, H., Yamagishi, T., Chiba, T., 1986. Prediction of complete mixing of liquid-fluidized binary solid particles. Chemical Engineering Science 41, 297–305. Mucha, P.J., Tee, S.-Y., Weitz, D.A., Shraiman, B.I., Brenner, M.P., 2004. A model for velocity fluctuations in sedimentation. Journal of Fluid Mechanics 501, 71–104.

5589

Nicolai, H., Herzhaft, B., Hinch, E.J., Oger, L., Guazzelli, É., 1995. Particle velocity fluctuations and hydrodynamic self-diffusion of sedimenting nonBrownian spheres. Physics of Fluids 7, 2–23. Peysson, Y., Guazzelli, É., 1999. Velocity fluctuations in a bidisperse sedimenting suspension. Physics of Fluids 11, 1953–1955. Pickard, D.K., Tory, E.M., 1974. A Markov model for sedimentation. 24th Canadian Chemical Engineering Conference. Pickard, D.K., Tory, E.M., 1977. A Markov model for sedimentation. Journal of Mathematical Analysis and Applications 60, 349–369. Pickard, D.K., Tory, E.M., 1979. Experimental implications of a Markov model for sedimentation. Journal of Mathematical Analysis and Applications 72, 150–176. Pickard, D.K., Tory, E.M., 1987. A Markov model for sedimentation: fundamental issues and insights, In: MacNeill, I.B., Umphrey, G.J. (Eds.), Advances in the Statistical Sciences, Stochastic Hydrology, vol. IV, D. Reidel, Dordrecht, pp. 1–25. Pickard, D.K., Tory, E.M., Tuckman, B.A., 1987. A three-parameter Markov model for sedimentation II. Simulation of transit times and comparison with experimental results. Powder Technology 49, 227–240. Ramirez, W.F., Galvin, K.P., 2005. Dynamic model of multi-species segregation and dispersion in liquid fluidized beds. A.I.Ch.E. Journal 51, 2103–2108. Richardson, J.F., Zaki, W.N., 1954. Sedimentation and fluidization. Part I. Transactions of the Institution of Chemical Engineers (London) 32, 35–53. Scott, K.J., 1984. Hindered settling of a suspension of spheres. CSIR Report CENG 497, Chemical Engineering Research Group, CSIR, Pretoria, South Africa. Scott, K.J., Mandersloot, W.G.B., 1979. The mean particle size in hindered settling of multisized particles. Powder Technology 24, 99–101. Shannon, P.T., Stroupe, E., Tory, E.M., 1963. Batch and continuous thickening. Basic theory, Solids flux for rigid spheres. Industrial and Engineering Chemistry Fundamentals 2, 203–211. Smith, T.N., 1966. The sedimentation of particles having the dispersion of sizes. Transactions of the Institution of Chemical Engineers 44, T153– T157. Tory, E.M., 2000. Stochastic sedimentation and hydrodynamic diffusion. Chemical Engineering Journal 80, 81–89. Tory, E.M., Ford, R.A., 2004. Simulation of sedimentation of bidisperse suspensions. International Journal of Mineral Processing 73, 119–130. Tory, E.M., Hesse, C.H., 1996. Theoretical and experimental evidence for a Markov model for sedimentation. In: Tory, E.M. (Ed.), Sedimentation of Small Particles in a Viscous Fluid. Computational Mechanics Publications, Southampton, pp. 241–281. Tory, E.M., Kamel, M.T., 1997. Mean velocities in polydisperse suspensions. Powder Technology 93, 199–207. Tory, E.M., Pickard, D.K., 1977. A three-parameter Markov model for sedimentation. Canadian Journal of Chemical Engineering 55, 655–665. Tory, E.M., Pickard, D.K., 1982. Extensions and refinements of a Markov model for sedimentation. Journal of Mathematical Analysis and Applications 86, 442–470. Tory, E.M., Kamel, M.T., Chan Man Fong, C.F., 1992. Sedimentation is container-size dependent. Powder Technology 73, 219–238. Tory, E.M., Bargieł, M., Honeycutt, R.L., 1994. A three-parameter Markov model for sedimentation III. A stochastic Runge–Kutta method for computing first-passage times. Powder Technology 80, 133–146. Verhoeven, J., 1963. Prediction of batch settling behaviour, B. Chemical Engineering Thesis, McMaster University. Voller, V.R., Porté-Agel, F., 2002. Moore’s law and numerical modeling. Journal of Computational Physics 197, 698–703.