- Email: [email protected]
On the divergence problem in calculating particle velocities in dilute dispersions of identical spheres II. Effect of a plane wall
Powder
Technology,
55 (1988)
51 - 59
51
On the Divergence Problem in Calculating Particle Velocities in Dilute Dispersions of Identical Spheres II. Effect of a Plane Wall E. M. TORY Department of Mathematics EOA 3C0 (Canada)
and Computer
Science,
Mount Allison
University,
Sackville,
New Brunswick
and M. T. KAMEL* Division of Mathematics, Engineering New Brunswick E2L 4L5 (Canada)
and Computer
Science,
University
of New Brunswick,
Saint John,
This paper is dedicated to the memory of David K. Pickard (Received August 10,1987;
in revised form December 7,1987)
SUMMARY
Consider a very dilute random dispersion of identical spheres sedimenting slowly in a viscous fluid toward a horizontal plane wall. This wall imposes the return flow and yields the usual result for the effect of long-range interactions on mean velocity. Since the usual method of imposing return flow leads to the divergence of the variance of velocity, we used this method to calculate the covariance matrix. Very lengthy, but straightforward, analytical integration and algebraic manipulation yield a diagonal matrix whose entries are finite when the test sphere is a finite distance z from the plane. The trace of this matrix is the variance of velocity, given by Var(u)=$fJ
I
$z-
119 y-&
+o
1 iZ
I
This vanishes as the solids fraction 4 approaches zero, confirming a key postulate of the Pickard-Tory model for sedimentation. However, its linear dependence on z seems unphysical, and additional terms should be derived. Additional boundaries may also be necessary. INTRODUCTION
Treatments of physical processes often assume a smooth, uniform state. Yet many natural phenomena exploit the diversity of *On leave at The American University in Cairo, Cairo (Egypt). 0032-5910/88/$3.50
the environment: a crack becomes a chasm; a trickle becomes a torrent. When channelling occurs in slurries, they subside much faster than when the fluid is expelled normally [ 11. In sedimenting dispersions, random variations in the local concentration of particles lead to a distribution of particle velocities [2]. Small clusters in dilute dispersions are augmented by a ‘demixing’ process [3,4] and the initial distribution evolves toward a steady state [ 2, 31. Experimental measurements indicate that this distribution is approximately normal [ 5 - 71 with a mean which greatly exceeds the Stokes velocity [8,9]. It should be possible to derive these results from the creeping-motion equations [lo]. It is well known that the long-range interaction term leads to the divergence of the calculated sedimentation velocity of an unbounded dispersion of spheres [ 111. It is also well known that this problem can be overcome by a renormalization incorporating the return flow [ 111. More recently, it has been shown that a divergence problem also exists in the calculation of the variance of velocity [ 12, 131 and that this problem cannot be overcome by the same renormalization [ 131. It has been suggested that the problem might be resolved by imposing a constraint on the nature of the return flow [ 131. In practice, of course, dispersions are bounded and the boundaries impose the return flow [ 141. However, velocities of spheres near boundaries are atypical [lo], so the goal is an explicit evaluation of the influ@ Elsevier Sequoia/Printed in The Netherlands
52
ence of container wall on spheres fur from these walls. Recently, this goal was achieved, at least in part, for spheres sedimenting towards a plane wall, in an otherwise unbounded fluid [15]. Remarkably, the explicit wall effect on the mean velocity vanishes in the limit, yet the calculation agrees with Batchelor’s [ll] for the influence of long-range interactions [15]. This indicates that the implicit effect of the plane wall is the imposition of return flow. Furthermore, it appears that the wall also imposes the correction for the conditionally convergent term in rp3. The success of this method of treating longrange interactions in the calculation of the mean velocity of sedimentation suggests that it be used to compute the variance. In carrying out these computations, we make extensive use of the work of Tory and Pickard [13] and Beenakker and Mazur [ 151. MATHEMATICAL
FORMULATION
We examine the sedimentation of a semiinfinite random dispersion of identical spheres toward a horizontal plane wall located at the X-Y plane. Consider a test region (in this dispersion) which is a cylinder of height H and radius P + (I. Suppose that a test sphere (of radius a) is at height 1where 3a < 1 < H - 3a with H > 6~. The axis of the cylinder passes through the centre of the test sphere. A vertical cross-section of this cylindrical region is shown in Fig. 1. Using only spheres which lie entirely within the test region (i,e., ignoring the effect of all spheres outside it), we could compute the mean and variance of the velocity of the test sphere. Since the test region is finite, these moments would be finite. Taking larger and larger cylinders produces a sequence of values which may or may not converge. When the integrals required to evaluate the mean velocity are taken in a particular order, the sequence does converge. That is, the improper [k + (li -
+ (Zi -
1)R,‘r^i]
+ ~
Z)Sc’$i + 21~S~2k -
+ li)Si’s^i -
i k
1
I
---I
t-p
0
/-: 0
\
T
‘\I ,
7fa
I
L
;
-1 H
-1 20
I-
I .-.--A L____________!-__J I
j----
P+a
----I
Fig. 1. Vertical cross-section of cylindrical test volume. Sphere centres are distributed uniformly throughout the volume whose cross-section is the connected region enclosed by the two sets of dashed lines.
multiple integral exists. Furthermore, the limit as 1+ 00 yields the accepted value of the mean velocity [ll] to the appropriate order. The centres of spheres (other than the test sphere) which lie in the test region are contained in a cylinder (with a spherical hole in it) bounded by horizontal planes at a and H - a and the lateral surface at a distance P from the axis which passes through the test sphere. The hole arises because the spheres are impenetrable, so none of their centres is within 2a of the centre of the test sphere. The mobility tensors governing the test sphere are given by Beenakker et al. [ 161. Since no torque is applied, we need only the translational mobility tensors BOand Bi which are combined as pTT in [16]. A straightforward calculation, beginning with eqn. (6.1) of Beenakker et al. [ 161, shows that the test sphere contributes
to its own velocity, and each of the other spheres contributes
3~k - 3(t'i - Z)R,'~~] 61li(l+ Zi)S,3s^i]
+ i a3Sc5(l + li)‘[-k
-
5(1+ li)Si’$i]
(2)
53
to the velocity of the test sphere. Here, K is the force exerted by the fluid on each sphere, Ri is the vector from the test sphere to the ith, Ri = IRJ, ii = Ri/Ri, u. is the Stokes velocity, k is a unit vector in the vertical direction, and Ii is the height of the ith sphere. All positions are defined in terms of the centres of the spheres. Also, ii =
Ri - 21ik
(3)
si
where Si = (Rt
(4)
+ 4Z1i)1’2
Equation (2) contains terms of order (a/Ri)m(a/Si)” with m + n < 3 [16]. We now introduce dimensionless variables h = H/a - 1, p = P/a, ri = RJa, Si = Si/a and z = Z/a, and a dimensionless system of cylindrical co-ordinates (p$,c) centred at the origin of the test sphere. Then, as illustrated in Fig. 2, ri = ipi
COS 8i
+ ipi
Sbl8i
+
(5)
kci
s^i= {ri - 2(ci + z)k}/si Si =
[pi2
+ (si
+
(6)
22)2]“2
(7)
If there are, by chance, exactly N other spheres in the test region, the dimensionless velocity u of the test sphere is given by
3 - si-’ [k + sisi-‘Bi + 2(ri + z)2siA2k - 6z(s, + 2Z)SiP3s^i] 4
-
3 + :si-3 2
(ci + 2z)siP’$i -
f k
1
+ $ si-“(ci + 2z)2 [m-k - 5(ci + 2z)Si’ii] Y
I !
OUTLINE OF THE COMPUTATIONAL
(8) PROBLEM
Equation (8) can be written as
T L-+2 -2(C+z)k
I
Fig. 2. Vector diagram for st. The vector ri goes from the centre of the test sphere to the centre of the ith. The downward vector from the ith sphere reaches its mirror image below the plane.
u = -be
-
gbi i=1
(9)
Thus,
is a symmetric dyadic (second-rank tensor) whose coefficients form a 3 X 3 matrix of products of velocity components. Equation (10 can be written as
54
UU=bObO +bo
and
fb,+
Cov[t<~W,l = t(b)&(b)V@‘O
i=l
+
flb,i+i=l5
From eqn. (17),
C bib]
(11)
j#i
If bi and b, are independent and identically distributed (i.i.d.), then C(b@j) = c(b,)&(bj) = c(b)e(b)
i#j
(12)
It follows that &(uuW) = b,b, +
+ N[b,c(b)
(22) For a Poisson process, c (N) = Var(N) = XV, so eqns. (20 and (19) become
&(u) = -b,
v=7T
COV(UW) = e(UUlN) - [(UlN)C(UlN)
(14)
From eqns. (9) and (12), -NL(b)
(15)
Hence + NW,&(b)
+
t(bP,l
t(b)&(b)1
(17)
Note that the dyadic in the brackets is exactly analogous to the usual scalar expression for the variance. To ensure that the distribution of sphere centres is i.i.d., we ignore the spatial interference of the N spheres with each other, i.e., we allow them to overlap each other, but not the test sphere. Hence, within the accessible region, their centres are distributed according to a Poisson process with intensity X = 3$/( 47~) [l3]. Since eqn. (8) contains only two-sphere interactions, it follows that the his are also i.i.d. The unconditional values are obtained by iterating expectations. Thus [ 131,
C(u) = e[e @uvl
(19)
From eqn. (15), &(u) = -6,
-
&(N)&(b)
p2(h-l)-
p
1
(26)
p+-= h+m
where
G(w, h) =
JJ&(r)dV
(27)
V
The expected value for any g dependent upon the position of any one sphere in the accessible region is &v(g)
= ;
G(w, h)
(28)
It follows that 6(u) = -b.
-
z
G,(z)
(29)
and Cov(u) = g
G,(z)
(30)
(18)
and Cov(u) = Cov[ c (UIN)] + & [Cov(u IN)]
[
G(z) = lim G(z,p,h)
and -
(24)
The spheres within this region are located independently, each with uniform probability density l/V. Let r = (p,e,c) denote generically, in cylindrical co-ordinates, the random position of any one of these spheres. Define
(16) Cov(uW) =N[&(bb)
(23)
and
(13)
= b,b,
- AVc(b)
The dimensionless volume of the accessible region is
Now
e(u[N)L(ulN)
(N)[ L(bb) - &(b)C(b)]
t [Cov(ulN)] =
Cov(u) = AVL (66)
+ C(b)&
t(bb)l +W’J- l)t(b)&@)
C(ulN) = -b.
(21)
(20)
wheregi = b andg, = bb. Following Beenakker and Mazur [ 151, we divide the accessible region into three pieces: lower, middle, and upper with respect to height. Then, referring to Fig. 1 and the definitions of our dimensionless variables, we see that
55
h
P
277
+
gV)p
sss 2 0
de
dp
3
dt
(31)
0
The variable lower limit,Jm, for the middle region corresponds to the position of the centre of a sphere which touches the test sphere. Spheres in the lower and upper regions can take any position from the axis of the cylinder (p = 0) outwards to p. Then, G(z) = i2
i
jng(r)p
l-r0
dt’ dp dc + j -2
0
r
jng(r)p
JFF
0
de dp dS + j= p &(r)p 2
0
de dp dl
0
3
(32)
where the improper integrals are expected to exist in the sense of limits of the simplified expressions resulting from the finite case.
RESULTS
Computation of the moments is straightforward, but tedious. The starting point for both is b
=
4
r-lk
- ;(<
+
i
y3k
+ 2z)2s-5k
+
i
fr3r
+ ;({
-
i
{rp5r
-
i
s-‘k
-
i
+ 22)[1 + 32({ + z)]s_5s -
[I
+
T({
3(c
+
z)2]s-3k
+ 2#s7s
-
f sSP3S (33)
which is easily obtained from the term following Cr= I in eqn. (8). Calculation of the mean velocity uses g(r) = b in eqn. (32). The first integration shows that the terms in i and j are zero. The second integration yields terms which cancel exactly in the upper and lower regions. The second and third integrations for the central region are fairly lengthy, but algebraic manipulation greatly simplifies the results. The final result is simply
!I
-5+;-s+$
(34)
k
However, the terms in zP3 and zP5 would be changed by including higher-order terms in the expression for b. Disregarding these non-significant terms, our result agrees with that of Beenakker and Mazur [15]. Cov(u) is computed by substituting g(r) = bb in eqn. (32). Though the calculations are exceedingly lengthy (involving hundreds of pages), each individual integration can be done analytically and algebraic manipulations simplify the results considerably at each step. The first integration, which is relatively simple, shows that the covariance matrix is diagonal, i.e., the velocity components u,, uy , and u, are uncorrelated. For each section, the second and third integrations produce many terms of 0(z) or higher which cancel exactly. After these cancellations, only a few terms of O(1) or higher remain. These are shown in the Table. Many of these cancel exactly with terms from other sections or combine to yield terms of O(1). For large z, the final result is Cov(u)
= $I
gz-
g)(ii+jj)+
(Ez-
g)kk+O($)]
(35)
The trace of the covariance matrix is Tr[Cov(u)]
117 = I#I =.z
-
69 81 2 -160 (i-i + j-j) + ( 64
g)k.k
+ a(;)]
(36)
56 TABLE Summary
of results
Type
Section
Termsa
kk
Upper
27 --log(z+2)+ 16
kk
Lower
-
27 ;log2+ 81 -log(z-1)+ 32
81
-logz+ 16
kk
Middle
27 slog(z*+22)-
ii + jj
Upper
-
ii + jj
Lower
6392 A-1024
ii +jj
Middle
--
2972 1024
10240
aAll logarithms
249 512 1047 512
13159 5120
543
27
256
8
+
27s 32
-
27.~ Glog2+
+ ““slog(Z)+ 4096
e2(uollN) a=x,y,z
(37)
Since the covariance matrix is diagonal, we can define Var(UIN) = c Var(U&v) (Y
(38)
Then,
= t(u’ulN)
- &(ulN).~(uW)
(39)
is just the scalar invariant of eqn. (14). Following through the scalar invariants of subsequent equations, we see that Var(u) = Tr(Cov(u)]. Thus
[
27 plog2-
243.~ ---+-log2 512
;zlog(=2)
to the base e.
Var(u,IN) = E( z&W) -
Va.r(u)=cp
27 -log(z+l)+ 32
+ -log2-
which is just the scalar invariant of the covariante dyadic. Now,
Var(uIN)
543 256
-
27 -$0&*-l)-
+ -
19101
4052 512
gz-g +o (t 11 .z
(40)
DISCUSSION
The major result of this work is the calculated value for the variance of a dilute dispersion of identical spheres. Before discussing this result, we wish to make a few points regarding the mean velocity. Though set forth in a somewhat different way, our calculations for the mean confirm
the results of Beenakker and Mazur [ 151 by essentially the same method. They supply few details, but in so far as they are available, intermediate results from the two calculations agree. The null contribution of the upper and lower regions to eqn. (34) indicates that the mean velocity of a test sphere is determined entirely by the spheres cut by horizontal planes tangent to that sphere, i.e., those spheres whose centres lie between z - 2 and z + 2 when a test sphere of unit radius is centred at z. Since m
2n bp
ss
0 0
d6’ dp = 0
151> 2
(41)
this result is not affected by allowing $J(C)to vary for ItI > 2. That is, the concentration above and below this slice can be different from that within it without affecting the mean velocity. This justifies the assumption (Al in [2]) that configuration acts primarily through concentration with remaining effects incorporated in a stochastic structure. In particular, it supports the contention that it is often sufficient to use the concentration in thin horizontal slices because the concentra;ion and configuration there determine the nature of the return flow [2]. Thus, this result verifies Assumption A5 of [2] for very dilute dispersions of equal spheres.
57
Equation (40) indicates that the variance is finite, even in a semi-infinite dispersion, for all spheres which are at a finite distance from the horizontal plane wall. In particular, lim Var(u) = 0
(42)
4’0 confirming Assumption A8 of [17]. There are many intrinsic checks on the computations, including the cancellation of potentially divergent terms, the form of the algebraic simplifications, and the requirement that each component of velocity in each section have a positive variance. (Note that only the last entry of the Table is bounded as z -+ 00.) Also, each calculation was done at least twice by hand. As a final check, a locally enhanced version of muMATH (a symbolic manipulation package for microcomputers [IS]) was used to verify the calculations. Though wall and end effects are certainly important, it seems unlikely that the variance does indeed increase linearly with z. Laboratory experiments are normally carried out in vessels with diameters much smaller than their heights. Thus, wall effects dominate end effects, at least for single spheres [19]. Though wall effects are small for individual spheres, they profoundly affect the sedimentation of dilute dispersions [8, 91. There is also an end effect, but it is difficult to disentangle this from other influences [5]. The variability of particle velocities in dilute dispersions has been attributed to the closeness of the nearest neighbours, since the probability of the next-nearest neighbour also being close is negligible [ll]. When wall and end effects are absent, the velocity of a test sphere can also be enhanced by many spheres being slightly closer [4,20]. Since the calculated mean velocity remains constant with z to 6( l/z), a linear dependence of variance on z would imply that the velocity distribution becomes more and more spread out as we consider test spheres further and further from the horizontal plane. If this does not happen experimentally, we must seek an explanation in either the formulation of the problem or its solution. In the present formulation, the fluid is unbounded except for the plane wall. Thus, there is no free surface at the top and no bounding walls at the sides. It is possible that
the effect of these boundaries at infinity is not negligible. However, we have recently shown [21] that the variance in a spherical container [22,23] is unbounded as its radius becomes infinite. Thus, at least part of the difficulty appears to lie elsewhere. Interaction of spheres with a wall to the order given by Beenakker and Mazur [ 15,231 yields the appropriate mean correction, but it appears that higher-order interactions are necessary to yield a variance which is bounded for all spheres in a infinite dispersion. Specifically, we require interactions of two spheres with the wall. Three-sphere interactions affect the term in $2 [ll, 141 and it is the term in $J which requires correction. Beenakker et al. [ 16,221 provide a general method for deriving mobility tensors. We hope to use their technique to obtain specific expressions for higher-order interactions. Should these fail to yield a finite variance, we would have to reexamine the other possibilities.
ACKNOWLEDGEMENTS
This paper is based, in part, on a poster presentation at the First International Conference on Industrial and Applied Mathematics, Paris, June 29 - July 3, 1987. We thank M. Beattie and K. Sears for enhancing muMATH and using it to check our calculations. This research was supported by NSERC Grant A1268.
LIST OF SYMBOLS
a 6
b Cov(u)
G
Gl G2
g
radius of spheres translational mobility dyadic (second-rank tensor) vector defined by eqn. (8) covariance dyadic whose scalar coefficients form the covariance matrix function obtained by integrating g over test volume vector function obtained by integrating b over test volume dyadic function obtained by integrating bb wver test volume dummy function (stand-in for b = gl orbb
=&I
58
H h i, i Li k K 1 4
G I: P R R r
height of test cylinder dimensionless height = H/a - 1 dummy indices of summation orthogonal unit vectors in horizontal direction unit vector in vertical direction force exerted by the fluid on a sphere distance from test sphere to horizontal plane distance from ith sphere to horizontal plane exponent of a/R i number of spheres in test region exponent of CIlSi radius of accessible region dimensionless radius of accessible region distance from test sphere to another sphere vector from test sphere to another sphere dimensionless distance from test sphere to another sphere dimensionless vector from test sphere to another sphere unit vector in direction of r distance defined by eqn. (4) dimensionless distance defined by eqn. (7) dimensionless vector shown in Fig. 2 unit vector defined by eqn. (6) velocity of test sphere component of velocity (a, = x, y, orz) horizontal components of velocity vertical component of velocity Stokes velocity dyadic whose scalar coefficients form the 3 X 3 symmetric matrix of products of velocity components of the sphere dimensionless volume of test region variance of the velocity of the test sphere defined as the sum of the variances of the velocity components dimensionless distance from test sphere to horizontal plane vertical component of dimensionless distance from test sphere to another sphere angle in system of cylindrical coordinates centred at test sphere
h TT I-1 P
Script CY tv
number of sphere centres in a unit volume mobility tensor derived by Beenakker et al. [16] horizontal component of dimensionless distance from test sphere to another sphere volume fraction of solids expectation expectation computed over test volume conditional expectation terms of specified order terms which are small compared with those of specified order
Subscripts indicating test sphere 0 i indicating ith sphere
j
indicating jth sphere
REFERENCES 1 V. Mallareddy, M. Eng. Thesis, McMaster Univ. (1963). 2 D. K. Pickard and E. M. Tory, J. Math. Anal. A&., 60 (1977) 349. 3 B. Koglin, Proc. Conf. Particle Technology, Chicago (1973) 265. 4 M. T. Kamel, E. M. Tory and W. S. Jodrey, Powder Technol., 24 (1979) 19. 5 E. M. Tory and D. K. Pickard, Can. J. Chem. Eng., 55 (1977) 655. 6 D. K. Pickard, E. M. Tory and B. A. Tuckman, Powder Technol., 49 (1987) 227. 76. 7 B. Koglin, Chem. Zng. Tech., 43 (1971) 8 B. Koglin, D. Zng. Dissertation, Universitlt Karlsruhe (1971). Monogr., 79Z3, 235 (1976) 9 B. Koglin, Dechema 1589. Number 10 J. Happel and H. Brenner, Low Reynolds Hydrodynamics, Prentice-Hall, Englewood Cliffs, NJ, 1965. J. Fluid Mech., 52 (1972) 245. 11 G. K. Batchelor, 12 R. E. Caflisch and J. H. C. Luke, Physics of Fluids, 28 (1985) 759. 13 E. M. Tory and D. K. Pickard, Powder Technol., 47 (1986) 39 and 48 (1986) 189. AZChE H., 26 14 C. C. Reed and J. L. Anderson, (1980) 816. and P. Mazur, Physics of 15 C. W. J. Beenakker Fluids, 28 (1985) 767. 16 C. W. J. Beenakker, W. van Saarloos and P. Mazur, Physica, 127A (1984) 451. 17 E. M. Tory and D. K. Pickard, J. Math. Anal. Appl., 86 (1982) 442.
59
18 H. S. Wilf, Amer. Math. Monthly, 89 (1982) 4. 19 R. M. Sonshine, R. G. Cox and H. Brenner, Appl. Sci. Res., 16 (1966) 273, 325. 20 M. T. Kamel and E. M. Tory, Paper presented at 7th Can. Symp. Fluid Dynamics, Mount Allison Univ., June 2 - 4, 1986.
21 E. M. Tory and M. T. Kamel, submitted to Powder Technol. 22 C. W. J. Beenakker and P. Mazur, Physica, 131A (1985) 311. 23 C. W. J. Beenakker and P. Mazur, Physics of Fluids, 28 (1985) 3203.
Technology,
55 (1988)
51 - 59
51
On the Divergence Problem in Calculating Particle Velocities in Dilute Dispersions of Identical Spheres II. Effect of a Plane Wall E. M. TORY Department of Mathematics EOA 3C0 (Canada)
and Computer
Science,
Mount Allison
University,
Sackville,
New Brunswick
and M. T. KAMEL* Division of Mathematics, Engineering New Brunswick E2L 4L5 (Canada)
and Computer
Science,
University
of New Brunswick,
Saint John,
This paper is dedicated to the memory of David K. Pickard (Received August 10,1987;
in revised form December 7,1987)
SUMMARY
Consider a very dilute random dispersion of identical spheres sedimenting slowly in a viscous fluid toward a horizontal plane wall. This wall imposes the return flow and yields the usual result for the effect of long-range interactions on mean velocity. Since the usual method of imposing return flow leads to the divergence of the variance of velocity, we used this method to calculate the covariance matrix. Very lengthy, but straightforward, analytical integration and algebraic manipulation yield a diagonal matrix whose entries are finite when the test sphere is a finite distance z from the plane. The trace of this matrix is the variance of velocity, given by Var(u)=$fJ
I
$z-
119 y-&
+o
1 iZ
I
This vanishes as the solids fraction 4 approaches zero, confirming a key postulate of the Pickard-Tory model for sedimentation. However, its linear dependence on z seems unphysical, and additional terms should be derived. Additional boundaries may also be necessary. INTRODUCTION
Treatments of physical processes often assume a smooth, uniform state. Yet many natural phenomena exploit the diversity of *On leave at The American University in Cairo, Cairo (Egypt). 0032-5910/88/$3.50
the environment: a crack becomes a chasm; a trickle becomes a torrent. When channelling occurs in slurries, they subside much faster than when the fluid is expelled normally [ 11. In sedimenting dispersions, random variations in the local concentration of particles lead to a distribution of particle velocities [2]. Small clusters in dilute dispersions are augmented by a ‘demixing’ process [3,4] and the initial distribution evolves toward a steady state [ 2, 31. Experimental measurements indicate that this distribution is approximately normal [ 5 - 71 with a mean which greatly exceeds the Stokes velocity [8,9]. It should be possible to derive these results from the creeping-motion equations [lo]. It is well known that the long-range interaction term leads to the divergence of the calculated sedimentation velocity of an unbounded dispersion of spheres [ 111. It is also well known that this problem can be overcome by a renormalization incorporating the return flow [ 111. More recently, it has been shown that a divergence problem also exists in the calculation of the variance of velocity [ 12, 131 and that this problem cannot be overcome by the same renormalization [ 131. It has been suggested that the problem might be resolved by imposing a constraint on the nature of the return flow [ 131. In practice, of course, dispersions are bounded and the boundaries impose the return flow [ 141. However, velocities of spheres near boundaries are atypical [lo], so the goal is an explicit evaluation of the influ@ Elsevier Sequoia/Printed in The Netherlands
52
ence of container wall on spheres fur from these walls. Recently, this goal was achieved, at least in part, for spheres sedimenting towards a plane wall, in an otherwise unbounded fluid [15]. Remarkably, the explicit wall effect on the mean velocity vanishes in the limit, yet the calculation agrees with Batchelor’s [ll] for the influence of long-range interactions [15]. This indicates that the implicit effect of the plane wall is the imposition of return flow. Furthermore, it appears that the wall also imposes the correction for the conditionally convergent term in rp3. The success of this method of treating longrange interactions in the calculation of the mean velocity of sedimentation suggests that it be used to compute the variance. In carrying out these computations, we make extensive use of the work of Tory and Pickard [13] and Beenakker and Mazur [ 151. MATHEMATICAL
FORMULATION
We examine the sedimentation of a semiinfinite random dispersion of identical spheres toward a horizontal plane wall located at the X-Y plane. Consider a test region (in this dispersion) which is a cylinder of height H and radius P + (I. Suppose that a test sphere (of radius a) is at height 1where 3a < 1 < H - 3a with H > 6~. The axis of the cylinder passes through the centre of the test sphere. A vertical cross-section of this cylindrical region is shown in Fig. 1. Using only spheres which lie entirely within the test region (i,e., ignoring the effect of all spheres outside it), we could compute the mean and variance of the velocity of the test sphere. Since the test region is finite, these moments would be finite. Taking larger and larger cylinders produces a sequence of values which may or may not converge. When the integrals required to evaluate the mean velocity are taken in a particular order, the sequence does converge. That is, the improper [k + (li -
+ (Zi -
1)R,‘r^i]
+ ~
Z)Sc’$i + 21~S~2k -
+ li)Si’s^i -
i k
1
I
---I
t-p
0
/-: 0
\
T
‘\I ,
7fa
I
L
;
-1 H
-1 20
I-
I .-.--A L____________!-__J I
j----
P+a
----I
Fig. 1. Vertical cross-section of cylindrical test volume. Sphere centres are distributed uniformly throughout the volume whose cross-section is the connected region enclosed by the two sets of dashed lines.
multiple integral exists. Furthermore, the limit as 1+ 00 yields the accepted value of the mean velocity [ll] to the appropriate order. The centres of spheres (other than the test sphere) which lie in the test region are contained in a cylinder (with a spherical hole in it) bounded by horizontal planes at a and H - a and the lateral surface at a distance P from the axis which passes through the test sphere. The hole arises because the spheres are impenetrable, so none of their centres is within 2a of the centre of the test sphere. The mobility tensors governing the test sphere are given by Beenakker et al. [ 161. Since no torque is applied, we need only the translational mobility tensors BOand Bi which are combined as pTT in [16]. A straightforward calculation, beginning with eqn. (6.1) of Beenakker et al. [ 161, shows that the test sphere contributes
to its own velocity, and each of the other spheres contributes
3~k - 3(t'i - Z)R,'~~] 61li(l+ Zi)S,3s^i]
+ i a3Sc5(l + li)‘[-k
-
5(1+ li)Si’$i]
(2)
53
to the velocity of the test sphere. Here, K is the force exerted by the fluid on each sphere, Ri is the vector from the test sphere to the ith, Ri = IRJ, ii = Ri/Ri, u. is the Stokes velocity, k is a unit vector in the vertical direction, and Ii is the height of the ith sphere. All positions are defined in terms of the centres of the spheres. Also, ii =
Ri - 21ik
(3)
si
where Si = (Rt
(4)
+ 4Z1i)1’2
Equation (2) contains terms of order (a/Ri)m(a/Si)” with m + n < 3 [16]. We now introduce dimensionless variables h = H/a - 1, p = P/a, ri = RJa, Si = Si/a and z = Z/a, and a dimensionless system of cylindrical co-ordinates (p$,c) centred at the origin of the test sphere. Then, as illustrated in Fig. 2, ri = ipi
COS 8i
+ ipi
Sbl8i
+
(5)
kci
s^i= {ri - 2(ci + z)k}/si Si =
[pi2
+ (si
+
(6)
22)2]“2
(7)
If there are, by chance, exactly N other spheres in the test region, the dimensionless velocity u of the test sphere is given by
3 - si-’ [k + sisi-‘Bi + 2(ri + z)2siA2k - 6z(s, + 2Z)SiP3s^i] 4
-
3 + :si-3 2
(ci + 2z)siP’$i -
f k
1
+ $ si-“(ci + 2z)2 [m-k - 5(ci + 2z)Si’ii] Y
I !
OUTLINE OF THE COMPUTATIONAL
(8) PROBLEM
Equation (8) can be written as
T L-+2 -2(C+z)k
I
Fig. 2. Vector diagram for st. The vector ri goes from the centre of the test sphere to the centre of the ith. The downward vector from the ith sphere reaches its mirror image below the plane.
u = -be
-
gbi i=1
(9)
Thus,
is a symmetric dyadic (second-rank tensor) whose coefficients form a 3 X 3 matrix of products of velocity components. Equation (10 can be written as
54
UU=bObO +bo
and
fb,+
Cov[t<~W,l = t(b)&(b)V@‘O
i=l
+
flb,i+i=l5
From eqn. (17),
C bib]
(11)
j#i
If bi and b, are independent and identically distributed (i.i.d.), then C(b@j) = c(b,)&(bj) = c(b)e(b)
i#j
(12)
It follows that &(uuW) = b,b, +
+ N[b,c(b)
(22) For a Poisson process, c (N) = Var(N) = XV, so eqns. (20 and (19) become
&(u) = -b,
v=7T
COV(UW) = e(UUlN) - [(UlN)C(UlN)
(14)
From eqns. (9) and (12), -NL(b)
(15)
Hence + NW,&(b)
+
t(bP,l
t(b)&(b)1
(17)
Note that the dyadic in the brackets is exactly analogous to the usual scalar expression for the variance. To ensure that the distribution of sphere centres is i.i.d., we ignore the spatial interference of the N spheres with each other, i.e., we allow them to overlap each other, but not the test sphere. Hence, within the accessible region, their centres are distributed according to a Poisson process with intensity X = 3$/( 47~) [l3]. Since eqn. (8) contains only two-sphere interactions, it follows that the his are also i.i.d. The unconditional values are obtained by iterating expectations. Thus [ 131,
C(u) = e[e @uvl
(19)
From eqn. (15), &(u) = -6,
-
&(N)&(b)
p2(h-l)-
p
1
(26)
p+-= h+m
where
G(w, h) =
JJ&(r)dV
(27)
V
The expected value for any g dependent upon the position of any one sphere in the accessible region is &v(g)
= ;
G(w, h)
(28)
It follows that 6(u) = -b.
-
z
G,(z)
(29)
and Cov(u) = g
G,(z)
(30)
(18)
and Cov(u) = Cov[ c (UIN)] + & [Cov(u IN)]
[
G(z) = lim G(z,p,h)
and -
(24)
The spheres within this region are located independently, each with uniform probability density l/V. Let r = (p,e,c) denote generically, in cylindrical co-ordinates, the random position of any one of these spheres. Define
(16) Cov(uW) =N[&(bb)
(23)
and
(13)
= b,b,
- AVc(b)
The dimensionless volume of the accessible region is
Now
e(u[N)L(ulN)
(N)[ L(bb) - &(b)C(b)]
t [Cov(ulN)] =
Cov(u) = AVL (66)
+ C(b)&
t(bb)l +W’J- l)t(b)&@)
C(ulN) = -b.
(21)
(20)
wheregi = b andg, = bb. Following Beenakker and Mazur [ 151, we divide the accessible region into three pieces: lower, middle, and upper with respect to height. Then, referring to Fig. 1 and the definitions of our dimensionless variables, we see that
55
h
P
277
+
gV)p
sss 2 0
de
dp
3
dt
(31)
0
The variable lower limit,Jm, for the middle region corresponds to the position of the centre of a sphere which touches the test sphere. Spheres in the lower and upper regions can take any position from the axis of the cylinder (p = 0) outwards to p. Then, G(z) = i2
i
jng(r)p
l-r0
dt’ dp dc + j -2
0
r
jng(r)p
JFF
0
de dp dS + j= p &(r)p 2
0
de dp dl
0
3
(32)
where the improper integrals are expected to exist in the sense of limits of the simplified expressions resulting from the finite case.
RESULTS
Computation of the moments is straightforward, but tedious. The starting point for both is b
=
4
r-lk
- ;(<
+
i
y3k
+ 2z)2s-5k
+
i
fr3r
+ ;({
-
i
{rp5r
-
i
s-‘k
-
i
+ 22)[1 + 32({ + z)]s_5s -
[I
+
T({
3(c
+
z)2]s-3k
+ 2#s7s
-
f sSP3S (33)
which is easily obtained from the term following Cr= I in eqn. (8). Calculation of the mean velocity uses g(r) = b in eqn. (32). The first integration shows that the terms in i and j are zero. The second integration yields terms which cancel exactly in the upper and lower regions. The second and third integrations for the central region are fairly lengthy, but algebraic manipulation greatly simplifies the results. The final result is simply
!I
-5+;-s+$
(34)
k
However, the terms in zP3 and zP5 would be changed by including higher-order terms in the expression for b. Disregarding these non-significant terms, our result agrees with that of Beenakker and Mazur [15]. Cov(u) is computed by substituting g(r) = bb in eqn. (32). Though the calculations are exceedingly lengthy (involving hundreds of pages), each individual integration can be done analytically and algebraic manipulations simplify the results considerably at each step. The first integration, which is relatively simple, shows that the covariance matrix is diagonal, i.e., the velocity components u,, uy , and u, are uncorrelated. For each section, the second and third integrations produce many terms of 0(z) or higher which cancel exactly. After these cancellations, only a few terms of O(1) or higher remain. These are shown in the Table. Many of these cancel exactly with terms from other sections or combine to yield terms of O(1). For large z, the final result is Cov(u)
= $I
gz-
g)(ii+jj)+
(Ez-
g)kk+O($)]
(35)
The trace of the covariance matrix is Tr[Cov(u)]
117 = I#I =.z
-
69 81 2 -160 (i-i + j-j) + ( 64
g)k.k
+ a(;)]
(36)
56 TABLE Summary
of results
Type
Section
Termsa
kk
Upper
27 --log(z+2)+ 16
kk
Lower
-
27 ;log2+ 81 -log(z-1)+ 32
81
-logz+ 16
kk
Middle
27 slog(z*+22)-
ii + jj
Upper
-
ii + jj
Lower
6392 A-1024
ii +jj
Middle
--
2972 1024
10240
aAll logarithms
249 512 1047 512
13159 5120
543
27
256
8
+
27s 32
-
27.~ Glog2+
+ ““slog(Z)+ 4096
e2(uollN) a=x,y,z
(37)
Since the covariance matrix is diagonal, we can define Var(UIN) = c Var(U&v) (Y
(38)
Then,
= t(u’ulN)
- &(ulN).~(uW)
(39)
is just the scalar invariant of eqn. (14). Following through the scalar invariants of subsequent equations, we see that Var(u) = Tr(Cov(u)]. Thus
[
27 plog2-
243.~ ---+-log2 512
;zlog(=2)
to the base e.
Var(u,IN) = E( z&W) -
Va.r(u)=cp
27 -log(z+l)+ 32
+ -log2-
which is just the scalar invariant of the covariante dyadic. Now,
Var(uIN)
543 256
-
27 -$0&*-l)-
+ -
19101
4052 512
gz-g +o (t 11 .z
(40)
DISCUSSION
The major result of this work is the calculated value for the variance of a dilute dispersion of identical spheres. Before discussing this result, we wish to make a few points regarding the mean velocity. Though set forth in a somewhat different way, our calculations for the mean confirm
the results of Beenakker and Mazur [ 151 by essentially the same method. They supply few details, but in so far as they are available, intermediate results from the two calculations agree. The null contribution of the upper and lower regions to eqn. (34) indicates that the mean velocity of a test sphere is determined entirely by the spheres cut by horizontal planes tangent to that sphere, i.e., those spheres whose centres lie between z - 2 and z + 2 when a test sphere of unit radius is centred at z. Since m
2n bp
ss
0 0
d6’ dp = 0
151> 2
(41)
this result is not affected by allowing $J(C)to vary for ItI > 2. That is, the concentration above and below this slice can be different from that within it without affecting the mean velocity. This justifies the assumption (Al in [2]) that configuration acts primarily through concentration with remaining effects incorporated in a stochastic structure. In particular, it supports the contention that it is often sufficient to use the concentration in thin horizontal slices because the concentra;ion and configuration there determine the nature of the return flow [2]. Thus, this result verifies Assumption A5 of [2] for very dilute dispersions of equal spheres.
57
Equation (40) indicates that the variance is finite, even in a semi-infinite dispersion, for all spheres which are at a finite distance from the horizontal plane wall. In particular, lim Var(u) = 0
(42)
4’0 confirming Assumption A8 of [17]. There are many intrinsic checks on the computations, including the cancellation of potentially divergent terms, the form of the algebraic simplifications, and the requirement that each component of velocity in each section have a positive variance. (Note that only the last entry of the Table is bounded as z -+ 00.) Also, each calculation was done at least twice by hand. As a final check, a locally enhanced version of muMATH (a symbolic manipulation package for microcomputers [IS]) was used to verify the calculations. Though wall and end effects are certainly important, it seems unlikely that the variance does indeed increase linearly with z. Laboratory experiments are normally carried out in vessels with diameters much smaller than their heights. Thus, wall effects dominate end effects, at least for single spheres [19]. Though wall effects are small for individual spheres, they profoundly affect the sedimentation of dilute dispersions [8, 91. There is also an end effect, but it is difficult to disentangle this from other influences [5]. The variability of particle velocities in dilute dispersions has been attributed to the closeness of the nearest neighbours, since the probability of the next-nearest neighbour also being close is negligible [ll]. When wall and end effects are absent, the velocity of a test sphere can also be enhanced by many spheres being slightly closer [4,20]. Since the calculated mean velocity remains constant with z to 6( l/z), a linear dependence of variance on z would imply that the velocity distribution becomes more and more spread out as we consider test spheres further and further from the horizontal plane. If this does not happen experimentally, we must seek an explanation in either the formulation of the problem or its solution. In the present formulation, the fluid is unbounded except for the plane wall. Thus, there is no free surface at the top and no bounding walls at the sides. It is possible that
the effect of these boundaries at infinity is not negligible. However, we have recently shown [21] that the variance in a spherical container [22,23] is unbounded as its radius becomes infinite. Thus, at least part of the difficulty appears to lie elsewhere. Interaction of spheres with a wall to the order given by Beenakker and Mazur [ 15,231 yields the appropriate mean correction, but it appears that higher-order interactions are necessary to yield a variance which is bounded for all spheres in a infinite dispersion. Specifically, we require interactions of two spheres with the wall. Three-sphere interactions affect the term in $2 [ll, 141 and it is the term in $J which requires correction. Beenakker et al. [ 16,221 provide a general method for deriving mobility tensors. We hope to use their technique to obtain specific expressions for higher-order interactions. Should these fail to yield a finite variance, we would have to reexamine the other possibilities.
ACKNOWLEDGEMENTS
This paper is based, in part, on a poster presentation at the First International Conference on Industrial and Applied Mathematics, Paris, June 29 - July 3, 1987. We thank M. Beattie and K. Sears for enhancing muMATH and using it to check our calculations. This research was supported by NSERC Grant A1268.
LIST OF SYMBOLS
a 6
b Cov(u)
G
Gl G2
g
radius of spheres translational mobility dyadic (second-rank tensor) vector defined by eqn. (8) covariance dyadic whose scalar coefficients form the covariance matrix function obtained by integrating g over test volume vector function obtained by integrating b over test volume dyadic function obtained by integrating bb wver test volume dummy function (stand-in for b = gl orbb
=&I
58
H h i, i Li k K 1 4
G I: P R R r
height of test cylinder dimensionless height = H/a - 1 dummy indices of summation orthogonal unit vectors in horizontal direction unit vector in vertical direction force exerted by the fluid on a sphere distance from test sphere to horizontal plane distance from ith sphere to horizontal plane exponent of a/R i number of spheres in test region exponent of CIlSi radius of accessible region dimensionless radius of accessible region distance from test sphere to another sphere vector from test sphere to another sphere dimensionless distance from test sphere to another sphere dimensionless vector from test sphere to another sphere unit vector in direction of r distance defined by eqn. (4) dimensionless distance defined by eqn. (7) dimensionless vector shown in Fig. 2 unit vector defined by eqn. (6) velocity of test sphere component of velocity (a, = x, y, orz) horizontal components of velocity vertical component of velocity Stokes velocity dyadic whose scalar coefficients form the 3 X 3 symmetric matrix of products of velocity components of the sphere dimensionless volume of test region variance of the velocity of the test sphere defined as the sum of the variances of the velocity components dimensionless distance from test sphere to horizontal plane vertical component of dimensionless distance from test sphere to another sphere angle in system of cylindrical coordinates centred at test sphere
h TT I-1 P
Script CY tv
number of sphere centres in a unit volume mobility tensor derived by Beenakker et al. [16] horizontal component of dimensionless distance from test sphere to another sphere volume fraction of solids expectation expectation computed over test volume conditional expectation terms of specified order terms which are small compared with those of specified order
Subscripts indicating test sphere 0 i indicating ith sphere
j
indicating jth sphere
REFERENCES 1 V. Mallareddy, M. Eng. Thesis, McMaster Univ. (1963). 2 D. K. Pickard and E. M. Tory, J. Math. Anal. A&., 60 (1977) 349. 3 B. Koglin, Proc. Conf. Particle Technology, Chicago (1973) 265. 4 M. T. Kamel, E. M. Tory and W. S. Jodrey, Powder Technol., 24 (1979) 19. 5 E. M. Tory and D. K. Pickard, Can. J. Chem. Eng., 55 (1977) 655. 6 D. K. Pickard, E. M. Tory and B. A. Tuckman, Powder Technol., 49 (1987) 227. 76. 7 B. Koglin, Chem. Zng. Tech., 43 (1971) 8 B. Koglin, D. Zng. Dissertation, Universitlt Karlsruhe (1971). Monogr., 79Z3, 235 (1976) 9 B. Koglin, Dechema 1589. Number 10 J. Happel and H. Brenner, Low Reynolds Hydrodynamics, Prentice-Hall, Englewood Cliffs, NJ, 1965. J. Fluid Mech., 52 (1972) 245. 11 G. K. Batchelor, 12 R. E. Caflisch and J. H. C. Luke, Physics of Fluids, 28 (1985) 759. 13 E. M. Tory and D. K. Pickard, Powder Technol., 47 (1986) 39 and 48 (1986) 189. AZChE H., 26 14 C. C. Reed and J. L. Anderson, (1980) 816. and P. Mazur, Physics of 15 C. W. J. Beenakker Fluids, 28 (1985) 767. 16 C. W. J. Beenakker, W. van Saarloos and P. Mazur, Physica, 127A (1984) 451. 17 E. M. Tory and D. K. Pickard, J. Math. Anal. Appl., 86 (1982) 442.
59
18 H. S. Wilf, Amer. Math. Monthly, 89 (1982) 4. 19 R. M. Sonshine, R. G. Cox and H. Brenner, Appl. Sci. Res., 16 (1966) 273, 325. 20 M. T. Kamel and E. M. Tory, Paper presented at 7th Can. Symp. Fluid Dynamics, Mount Allison Univ., June 2 - 4, 1986.
21 E. M. Tory and M. T. Kamel, submitted to Powder Technol. 22 C. W. J. Beenakker and P. Mazur, Physica, 131A (1985) 311. 23 C. W. J. Beenakker and P. Mazur, Physics of Fluids, 28 (1985) 3203.