- Email: [email protected]
The distribution of velocities of a concentric sphere in a dilute dispersion of spheres sedimenting in a spherical container
159
Powder Technology, 74 (1993) 159-169
The distribution of’ velocities of a concentric sphere in a dilute dispersion of spheres sedimenting in a spherical container E. M. Tory, M. Bargiel Department of Mathematics and Computer Science, Mount Allison University, Sackville, NB, EOA 3C0 (Canada)
and M. T. Kamel Department of Mathematics, Statistics and Computer Science, Universe
of New Brunswick, Saint John, NB, E2L 4L5 (Canada)
(Received July 20, 1992; in revised form October 19, 1992)
Abstract We consider the motion of identical spheres in an incompressible fluid which is bounded by a spherical container wall. The test sphere is momentarily concentric with the container and the centres of the other spheres are distributed in the accessible region according to a Poisson process. We derive the third and fourth polyadic central moments of velocity and examine their components. These moments, together with the mean and covariance which are known from earlier work, provide a good description of the theoretical initial distribution. Estimates of the hard-sphere equilibrium distribution are obtained from computer simulations using a fixed number of hard spheres. Simulated values agreed closely with the theoretical, indicating that the Poisson assumption, which greatly simplifies the analysis, has little effect on the results. In small containers, the vertical component of velocity is skewed downward and is platykurtic. In large containers, the distributions of all components approach normality.
Introduction It has been generally believed that identical spheres settle with nearly equal velocities aside from a few cases where neighbours are unusually close, but experimental measurements, summarized by Tory and Pickard [l] and Tory et al. [2], show enormous variability. This is supported by theoretical calculations of the variance of particle velocities in dilute dispersions [2-61. The equations of creeping flow imply that the configuration of the system of particles and boundaries determines the distribution of particle velocities [7] and this distribution varies both spatially and temporally [8]. Furthermore, many phenomena which cannot be explained on the basis of Kynch’s theory [9] are accounted for by considering the distribution of particle velocities [l]. Though distributions of instantaneous velocities have not been measured, related measurements are compatible with a distribution which is approximately normal [lo, 111. Following Beenakker and Mazur [12] and Tory and Kamel [6], we consider the motion of identical spheres of unit radius in an incompressible fluid which is bounded by a spherical container wall of radius j3. The test sphere is momentarily concentric with the container
0032-5910/93/$6.00
and the centres of the other spheres are distributed in the accessible region according to a Poisson process with intensity A= 34/(4~) where 4 is the solids fraction [6]. We allow these spheres to overlap each other but not the test sphere or the boundary. In these circumstances, the mobility tensor assumes a particularly simple form [12] and the dimensionless mean and variance are easily calculated [6]. As noted by Tory and Kamel [6], the spherical case .is a useful prototype for more practical studies. Here we extend our theoretical calculations to the third and fourth central moments of velocity, thereby determining the skewness and kurtosis of the distribution. We use computer simulations to approximate the equilibrium hard-sphere distribution and compare the computed moments with the theoretical values. The use of the Poisson distribution leads to a variable number of spheres in the container. This is consistent with the idea of sampling from a much larger volume and leads to unconditional expectations which are simpler than expectations conditioned on a given number, N, of spheres in the volume [2,6]. However, the practical effect is nil because the difference consists entirely of terms which are too small to merit calculation [2].
0 1993 - Elsevier Sequoia. All rights reserved
160
Mean and variance of velocity Integration of eqn. (19) of Tory and Kamel [6] and insertion of the result into their eqn. (18) yields the dimensionless mean velocity 8(u)=
-b,++[2-9p-l+a(~-“)P
(1)
which agrees with their eqn. (20) but, unlike it, contains the term in p-l. To the order of approximation considered by Tory and Kamel [6], this agrees with the result of Beenakker and Mazur [12]. Since b0 is just R plus wall corrections, we are more interested in the second term which represents the mean effect of the other spheres on the test sphere. The variances of the components of velocity are given by [6] Co+)
= ii Var(u,) +J Var(z+,) + kk Var(u,) =
4
Here we are using the convenient notation of Frankel and Brenner [15], but Tory’s definition [14], which differs slightly from theirs, is more suited to our needs. Then [5, 61
[(&B-gj+
(ii+j) G
Here b, represents the contribution of the ith sphere to the translational velocity of the test sphere while b. represents that of the test sphere to its own velocity. Noting [5, 61 that 8(bibj) = [&Y(b)]‘,i #j
(9
and counting the number of terms in each summation, we obtain
/3-‘)kk+a(8-2)].
(2) -AE(b3) -iV(NPolyadic moments of velocity
-N(N-
Since the particle velocity, u, is a vector, its central moments are polyadics [7]. For example, Cov(u) is a dyadic which provides information about the three components u,, u,, and u,. Its diagonality indicates that these components are uncorrelated and the values of the diagonal elements are the variances of the components. Defining the variance of velocity as the scalar quantity [6] Var(u) = Var(u,) + Var(%) + Var(u,)
(3)
implies that
1
(n times)
(5) operator
[14, 151
sum of the distinct permutations
1+3N+3N+3N(N-l)+N+3N(N-l) +N(N-
l)(N-- 2) = (N+ 1)’
(11)
as is obvious from the first line of that equation. It follows that no terms were lost in going from eqn. (8) to eqn. (10). Using eqn. (10) together with [6] $(z@)=
-b,-I@=(b)
(12)
(13) we obtain
~(4w1’lN = S(u’p) - pqup)~(u’pv)~ + 2[qup)]3
g#J -
=
-N{E’(b3) - [‘(b)a(b’)r+
2[8(b)13}
of b, b,. . .b,, . (6)
a[+- c)“l= Gwu
For example, [bIb;~=bIbj:+b2b,b2+b~b,.
The total number of terms in eqn. (8) is therefore
(14)
which is the third central moment given that there are N spheres in the accessible region. The unconditional third central moment is
lw2...ur =
(10)
8(~~~N)=b;+~b,,8(b)~+N~(b~)+N(N-l)[~(b)]~,
z/3z + z p-‘+@(p-“) . (4) =4 [ 112 We need a convenient notation to deal with higher moments because the expressions rapidly become unwieldy. Following Grad [13], we define the n-adic
and the symmetrization
2)[8(b)13
and
Var(u) = Tr[Cov(u)]
b”=bb...b
l)(N-
1)[P(b)8(b2)r
(7)
=-
- Gw
+ %wI - P131W
A@=‘(b3) - ~~(b)~(b’)~ + 2[8(b)13)
161
+ lw(4w’Pr
+ t~(dw - P13}
(15)
This result is obtained by expanding the expression in the first line, using eqn. (14) with 8(N) = AV, and noting that
~~[~-$(~lN)l[~(~lN)-P12~1N)=0
(16)
The fourth to seventh terms in the second part of eqn. (15) can be written as
AV)3][fF(b)]3= - hI$Y(b)]3.
In addition to the values given in the table in Tory and Kamel [6], we require a
2w
(i-k)3i3 sin 8 do de= & ji2k+j2klp + + k3.
(26)
(17)
The final integration (with respect to r) yields Ep(b3). Then eqns. (19) and (24) give (18)
(15), (17) and (18) yield
8[(u - /L)3]= -AK%(P)
(25)
0 0
and the last term (see Appendix) is
Equations
i=isin8cosw+jsinOsinw+kcosO.
G ss
- Plln
= - hVEO[GP(b)kqb2)~ - 3[8(b)13}
-$[(N--
but very tedious. As evident in eqn. (23), b has subunits R and (i.k)i where
3
- Lw4wwwv
S(lMu”ln?
is straightforward
g?(u - 4’1
(19)
which is much simpler than the conditional expectation given by eqn. (14). The techniques used to obtain eqn. (14) yield
(27)
SE+ - ~(#FVl
Writing
= N?(V) - MC?(b)8@3)P + N(N- l)~[~(zP)]2r + (2N-N2)U~(bZ)[~(b)]2~+
u-p=iu,+jU,+k(u,-p)
(3N2-- 6N)[8(b)p
(20) A straightforward but tedious application of the method used to obtain eqn. (19) leads to J%ZJ - d41= AJ@‘(b4) + (W2WV2)l?
(21)
in which the last term consists of three distinct permutations. This is also much simuler than the corresponding conditional moment. 1; deriving eqn. (21), we have used (see Appendix) 8[(N-- AV)‘] = 3(AVj2 + AV.
(22)
(W
and expanding (U-P)’ yields 27 terms. By symmetry, the expected values of twenty are zero. This leaves 8[(u - A31 = lF’WW(uI
- ~11
+ ~2W-~[~;(~z- cl)1+k3WUz - 1-4~1 (29) Equating terms in eqns. (27) and (29), we obtain g?~,z(k - PL)l
= ~b,2(~z- 41 = -6[gln(y)-
z
+ gfi-l]
(30)
and Calculation
Z:[(% - IQ31
of the third central moment
Given [6] that b= ;r-‘[k+(:.k)+ +
Since the third central moment is zero for a normal distribution [16], eqn. (31) is a measure of non-normality.
r-3[3(i4)i-R]
-$8’+;B’n (
+
;
1
~r’(P-‘-p-‘)[2k-(i.k)~,
Calculation (23)
Equation
of the fourth central moment
(25) yields
162 TT 2lT
3
G ss
(i.k)4i4
Equating coefficients in eqns. (33), (34) and (36), and using eqn. (2),
sin 8 dw d6’
0 0
Z@x) - 3[Var(u,)]’
+4+
&
[i2k2+j2kz~+
$
(i4+j4)+
-$-j
[i'j'lp
(32) Using the table in Tory and Kamel[6] together with eqns. (26) and (32), we can complete the first two integrations intheevaluationof8(b4).Thefinalintegrationissimplified by considering only those terms which are a( 1) to @‘(p- ‘). All of the former arise from the lower limit of integration. The final result is AV8(b4)
(( 8601
_#j _ -
10035200
+ SOB-'
= S($) - 3[Var(u,)]’ 8601 6561 P-1 10035200 + 985600
= -34
(37)
and 8[(k -
CL)~IW4412
(38) The expressions on the left-hand side of eqns. (37) and (38) are the fourth ‘cumulants of the distribution of the horizontal and vertical components of velocity [16]. Since these are zero for a normal distribution [ 161,these equations provide measures of non-normality.
(3i4+3j4+[izjq7 )
Covariances
(33) S[(u - p)“] can be evaluated from eqn. (6) of Tory and Kamel and eqns. (2), (21) and (33). The expansion of (u - p)4 contains 81 terms but, by symmetry, the expectations of 60 of these are zero. This leaves
and correlations
The magnitudes of the off-diagonal elements in eqns. (29) and (30) indicate the strengths of the covariances between different components. Since [17]
c0v(u,z&)= iY{[u*”- squx”)](u*- p)} = ~b,z(uz - CL)1
(39)
3
the values of Cov(u,2u,) and Cov($uJ eqn. (30). We also have
are given by
Cov(u,zu,‘) = @ux”4) - a(u,)s($) = +&$) From eqn. (6) of Tory and Kamel [6] and eqn. (2), AW(b2)
=
i2
Var(u,) +j’ Var($) + k2 Var(u,) .
E-
-Var(uJ 8601
4 10035200
Var(u,) +
$$$I-’
(35)
)
(40)
Thus
W?&W2)12~ = (AV)2{ii?(b2)8(b2) + 8[b$(b2)b]
+$(bb ‘bb ‘)}
-4
=
(_El+?!&-
‘I[%
-l*(Y)]}.
(41)
It follows from eqns. (37) and (40) that Var(u$ = Var($) = 8(u,4) - [Var(u,)]’
(36)
2: 2[Var(u,)]’ - 34 I~~~OO (
+ $$OB-l
1 (42)
163
The distribution of velocities
Thus we can calculate [17] ~rrP%
m(k - CL)1 cc)1= [Var@Z) x Var@z)]‘”
-
(43)
since Var(u, - p) = Var(u,). Also Corr(u$Q
=
Cov(u,‘u,“) Var(uz) *
(44
Similarly VarNu, -
Jm4 - Id31=IW(uJ13~
PYI
= 2fvar(uZ)12
-
&&j@, +fi-’
(45) and hence Corr[u,2(u, - p)*] =
cc+:@, - I-4’1 {Var(u,‘) Var[(u, - P)*]}~~
can be calculated.
For f12> 4-l,
CYorr(u$,‘) = =L
(46)
these become
~]u,(UZ - PII Jz Var(uX)[Var(uZ)]‘/2
= -28.17+-1nj?-3n
(51) 1
WY - A41 -3=5.54&-l/3-2. Wr(uz>12
In
(47)
q44) PW41*
-3=
(52)
components,
it is
-0.49C#+p-2.
(53)
Thus, the distributions of u, and u,. are very slightly leptokurtic [19]. Tory and Kamel [5, 61 characterized Var(u) by eqn. (4). While the choice for contraction of 8[(u -P)~] is less obvious, we take the dot product of the first two vectors in each triad. (The symmetry makes the choice immaterial.) This yields an overall measure of skewness
cov(uI”u,“)
16011 XY
2W-(4H2 -8.19x
.
(
For the horizontal
cOrr[u,Z(u, - p)] =
!$A
;pp-lqj-3qn
It follows from eqn. (38) that the distribution of u, is platykurtic [19]. The excess kurtosis [18] is
437337 8800
166133 +di
It follows from eqn. (23) and the symmetry of the expected spatial configuration that the distributions of u, and u,, are symmetrical for all values of 4 and p. Equation (31) indicates that u, is negatively skewed, i.e., toward downward velocities. Using Var(u,) from eqn. (2), we obtain the coefficient of skewness [18] for u,:
10-2+-1~-2
(48)
Sk(u) = -
3f.?
k
p
and Corr[u~(u, - p)“] =
~+a4 - Id21 2 Var(u,) Var(u,)
=0.556~$-‘p-~
56&’ 2=- 9
B-1 (
(49)
Equation (30) indicates that rapidly falling spheres have larger horizontal velocities than those which fall slowly or move upward. While this quantity grows in absolute value as p increases, it is asymptotically zero compared to the variances. As indicated by eqn. (40), the covariance of u,’ and u,” is almost zero. This appears to reflect two opposing effects. It is easily shown that Corr(sin’ w, cos* w) = - 1 .
(#+2p-3/21,
(50)
Hence large values of u,” are associated with small values of u,“. On the other hand, large values of u,’ and t$ are both associated with large values of (uZ- CL)’ as indicated by eqn. (41). It is clear from eqns. (47)-(49) that the correlations among all these velocities decrease rapidly toward zero with increasing p.
Tk
(54) )
This coefficient reflects the fact that all the features are related to the asymmetry in the vertical velocity. An alternative would be to weight the first term in eqn. (27) by the number of permutations, but this is more cumbersome. Of course, no contraction can retain all the information in eqn. (27). Since S[(u - p)*] consists of the pieces given by eqns. (33) and (36), there seems to be no simple formula for summarizing its value. However, we can summarize a measure of excess kurtosis. Equation (36) is closely related to the value subtracted in computing this quantity. This leaves the scalar invariant [20] of AI@(b4). The symmetry of all quantities in eqn. (33) means that we can use any consistent method of computing this invariant. Then [Pj2p contributes (i*i)(j.j) + (j j)(i*i) = 2 with similar contributions from [Pk’p and b”k’p. The result is
164
510151 250880 +p-1
Ku(u) =
other sphere is discarded and new random coordinates are chosen. The probability of placing these N spheres with no overlap is roughly
2 4
= 3.894-‘p -2
(55)
Equations (54) and (55) imply that the skewness and excess kurtosis are arbitrarily small for sufficiently large /3. It is also clear from eqn. (23) that pJn__P(b”)=d(l),
n=5, 6, . . . .
(56)
Consequently, we expect that the distribution of velocities will be approximately normal in containers which are very large compared to the size of the particles. For very dilute dispersions in small containers, the distributions of u, and u,, have slightly smaller tails and/or sharper peaks than the normal distribution [19]. In contrast, the distribution of u, is much more spread out and this persists to fairly large values of @. Simple calculations show that the conditional moments given by eqns. (14) and (20) are virtually identical to the unconditional moments expressed by eqns. (19) and (21) respectively. It is clear from eqn. (6) of Tory and Kamel [6] and eqn. (12) that eqn. (14) may be written as CY{[u --a(u]N)]3)N) = -NtP(b3). Using these same equations shows that
(57) together
Z{[u - S(u]N)]4]nr) -N’&(b4) +N(N-
with eqn. (19)
l)[[EP(b2)]2I”.
(58) Setting N=AV in eqns. (57) and (58) completes the comparison of eqns. (14) and (20) with (19) and (21) respectively.
Computer simulation of velocities The distribution of velocities in small containers can be estimated directly by simulation. Since a comparison between theoretical and empirical results for a Poisson process would serve only to test the accuracy of the random number generator, we replace that approximation by the hard-sphere model. The test sphere is placed at the centre of a spherical container of radius l/3
(59) and N spheres are randomly placed in the accessible region. Any attempted placement which overlaps an-
(l- $(l- $)...(l-y 4) =exp[ - (N;l)+l which corresponds to an approximate solution of the well-known birthday problem [22] with N people and a year of N/4 days. For += 10m3 and N= 1000, 2000, 5000 and 10 000, these probabilities are roughly 0.607, 0.368, 0.082 and 0.007 respectively. Thus, it is almost inevitable that some overlaps will occur in many of the simulations and hence that the overall hard-sphere distribution will differ slightly from the Poisson point approximation. However, these differences are important only insofar as they limit the number of spheres which are close to the central sphere. For $J= 10e3, the distribution of centres is very close to the equilibrium distribution of hard spheres. Thus, the distribution of velocities is correspondingly close. For this reason and because the simplifying assumptions which permit the derivation of theoretical values rapidly become untenable with increasing 4, we set += 10v3 for all simulations. (See Discussion.) Each data set consisted of 10 000 simulations. We obtained 15 sets for N= 1000 and 10 sets for N = 2000, 5000 or 10 000. Since theoretical values for the moments of u, and u,, are identical, and simulated values can differ only by chance fluctuations, we first compute the moments separately for each set, check that there are no significant differences, and then pool the values for u, and u,, On the graphs, they are all labelled as moments of u,. Figures l(a) and l(b) show the empirical distribution of u, and uy (combined) for N= 1000 and 10000, respectively for all data sets for these values of N. The symmetry is evident in both, while the greater spread is apparent in Fig. l(b). Figures 2(a) and 2(b) show the distribution of u, + b, for 1000 and 10 000 neighbours respectively. The skewness is clear in both and, once again, the greater spread is apparent. Values for N = 2000 and 5000 (not shown) fall nicely between those for 1000 and 10 000 in Figs. 1 and 2. For each configuration, we evaluated 1~+b, the effect of the other spheres on the central sphere. For each set of 10 000 configurations, we calculated estimates of @@)+b, and the second, third and fourth (conditional) central moments of the components of u. The sample means (i.e., the means for each set) are approximately normally distributed [19], so we used the t-distribution (with 14 or 9 degrees of freedom) to
165
1000 spheres
1000 spheres
c
0.0
(4
2,
+ bo
10000 spheres
10000 spheres
1
(“1
1
(b)
u,+
bo
Fig. 1. Frequency distribution of the horizontal components of velocity. Values for u, and u,. have been pooled. In all Figs. in this paper, += 10w3. (a) 1000 neighbours, (b) 10 000 neighbours.
Fig. 2. Frequency distribution of the vertical component of velocity. The constant, -b, representing the velocity of the central sphere in the absence of any neighbours, has been subtracted froin ui to highlight the effect of neighbours. (a) 1000 neighbours, (b) 10 000 neighbours.
calculate 95% confidence intervals for the overall means of each of these quantities [21]. Figure 3(a) shows the simulated values of $(u,) and 8(u,J, combined, while 3(b) shows those for a@,)+&. The central moments were calculated in the usual way [16, 181. Though we reduced the degrees of freedom by one to allow for the calculation of the mean from the data, this had no practical effect since each moment was based on 10 000 simulations. As noted above, the conditional moments differ negligibly from the unconditional. In each case, the simulated moments (taken to represent both the conditional and unconditional values) are compared to theoretical values for the unconditional moments. Figure 4(a) displays the values of Var(u,) and Var(u,), combined, while 4(b) and 4(c) show the values of Var(u,) and Var(u), respectively. Figure 5 shows the third moment. Figure 6(a) and (b) indicates the fourth cumulants of u, and uu (combined) and U, respectively.
In all cases but one, the confidence limits for a fixed number of hard spheres include or almost include the theoretical values for a Poisson process. In this one case, the confidence limits for the fourth cumulant of U, and u,, (combined) are very close together and the actual difference in values is extremely small. The small differences between theoretical and simulated values indicate that the simulation provides an empirical distribution of velocities which closely approximates that corresponding to a Poisson spatial distribution of centres.
Comparison
with other results
We have already noted [2] that our eqn. (4) is consistent with Caflisch and Luke’s [3] eqn. (3) which indicates that the variance of velocity is approximately proportional to
166
~,,,,,,,,,,,, -0.003
50
I ,,,,,, I ,,,,,,,,,,,,/,),,/,
0
100
150
200
0.00
250
B
(a)
0"
100
150
200
2
200
2 0
200
; ‘0
B
(a)
0.006
+
T--. N
F-Y N 3 z
‘~“~~“~1~~‘~‘(~“~“(~~‘~“I”~~~‘~‘~
50
ikO6
0.001
I)
0
-0.004
50
0
~"~""~I~"~~"~'I~'~"~~'~I~~~'~"" 100 150 200
@)
0.02
250
B
50
/“~~‘~“‘I~‘(~“~~‘I”~~‘~”
100
150 B
(b)
Fig. 3. Mean value of components of velocity (a) Horizontal. By symmetry, this value is exactly zero for both the hard-sphere and Poisson distribution of centres. Simulated values of u, and uy have been pooled. (b) Vertical. Tpe ordinate shows the effect of all neighbours on the velocity of the central sphere. By p = 225, the theoretical value is within 10v6 of the horizontal asymptote of 0.002.
NID$P’ = &#I .
(61)
Very recently, Ladd [23] calculated Cov(u) for spheres sedimenting in a square box with periodic boundary conditions. In our notation, his result for large N is I
Cov(u)=N*“#P[O.O66
2(ii +j) 4 0.823k.k] .
(62)
Since the numerical coefficients are different for different boundary conditions, we cannot expect quantitative agreement between eqns. (2) and (62). However, the coefficients are of comparable magnitude, viz., 0.072 3 versus 0.066 2 and 0.579 versus 0.823. Here we have used eqn. (61) and approximated our variance by the leading terms. Ladd’s Fig. 5 [23] shows Var(u,]N) 20.012 5 and Var(u,(N) = 0.153 for simulations with N = 108, 4 = 0.05, and a random distribution of sphere centres. (These
,
0.0 Cc)
50
~“~~“‘~I’~‘~“‘~~I”‘~J’~‘~I’~““”
100
150 P
Fig. 4. Variance of velocity. (a) Horizontal components. Values of u, and uYhave been pooled. (b) Vertical component. (c) Var(u) as defined by eqn. (3).
values are consistent with those shown in his Table 2 and Fig. 5 for steady-state values.) Insertion of this N and 4 in our eqns. (59) and (2) yields Var(u,)=0.0314 and Var(u,) = 0.175. (Though the equations for Var(u,) and Var(u,lN) differ, the correction factor [2] is neg-
167
2
0
z 0.00000
2 , -0.006 xN kc Q
I
-0.008
-O.O1 O5rn
250 B
Fig. 5. Third moment of velocity. The solid line and the small circles refer to the third central moment of u.. The dotted line gives the values of [fm(y)-
$$+
-$$$-t]?
which
is used in eqn. (54) to calculate the skewness of U.
ligiile even for N- 108.) Given the difference in boundary conditions and the fact that 4-0.05 is well beyond the range where our assumptions are valid, the orderof-magnitude agreement is the best that can be expected. There. are no experimental data which are directly comparable. with our theoretical results for Var(u,). However, an extensive analysis [2] of published data showed qualitative agreement between theory and experiment. Little work has been reported on the higher moments of velocity. Recently, however, Ladd [23] obtained values of the fourth cumulants of U, and U, at several concentrations. Dividing our eqns. (37) and (38) by 3 Var(u,) and 3 Var(u,) respectively yields the cumulants plotted in his Fig. 8. For r$= 0.05, our values, 0.068 and 0.669 respectively, are much larger than his because the Poisson approximation is then grossly inadequate, permitting too many close neighbours and thereby inflating the value of JY[(u-p)‘]. Nevertheless, Ladd’s results complement ours by extending the computed near-normality to dispersions with fewer particles and higher concentrations, viz., 108 spheres at 4=0.05: At moderate to large concentrations (0.25 Q 4 Q 0.50), nonGaussian effects were substantial, especially for the initial (equilibrium) distributions [23]. ,According to Caflisch and Luke [3], it can be shown that the distribution of NW% approaches a Gaussian as N? 03 with 4 held constant. Since Var(N- ‘%) = N - In Var(u)
(63)
the variance of N - “‘u, which is approximately proportional to &“, is finite in the limit of an infinite volume. It follows, of course, that the distribution of
0.000 @I
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,I,,:,,,,,,
50
100
150
200
250
B
Fig. 6. Fourth cumuhmts of components ofvelocity. (a) Horizontal. The theoretical curve, given by eqn. (37), has a horizontal asymptote of roughly - 2.57 X 10e6. (b) Vertical. The theoretical curve, given by eqn. (38), has a horizontal asymptote of roughly 1.85~10-~ which is approached very slowly on account of the logarithmic term.
u is approximately normal in all containers which are sufficiently large in every dimension.
Discussion Near-normality of the steady-state distribution of velocities implies that the distribution of distances travelled in an interval of time is approximately normal [lo] and the distribution of times to traverse a fixed distance is approximately log-normal [ll]. Thus, the simulated near-normality of steady-state distributions [23] is consistent with the experimental results of Tory and Pickard [lo] and Koglin [24]. The near-normality of both the initial and steady-state distributions is consistent with the three-parameter Markov model of Tory and Pickard [lo, 111 which postulates a onedimensional Langevin equation for u,.
168
An almost-normal distribution is well characterized by its mean and its second, third and fourth central moments. Thus, we made no effort to find higher moments. These would be prohibitively time-consuming for hand calculation, but might be accessible through symbolic manipulation routines. We used Maple [25] on a SUN SPARCstation 1+ to check the calculations for g’[(u -p>“], n =2, 3, 4. Taking i, j, k as scalar variables collapses the terms in [. B and simplifies the symbolic computations. Nevertheless, memory limitations made it necessary to compute the fourth moment in pieces. The symbolic computation gives all the terms for B(V), but the value of b expressed by eqn. (23) is only an approximation. Thus, we discarded all o@ -‘) terms. Higher moments could probably be calculated by pruning the expression to eliminate all terms which would be insignificant and by breaking down the integrals into several pieces. Despite the small values of p, velocities vary substantially. The computational effort for our simulations increases as /3’. Fortunately, the curves are already close to normal for N= 10 000 (/3=216) and the theoretical values of the moments indicate the approach to normality with increasing p. In our model, the effect of each neighbour is independent of that of every other neighbour. That is, we consider only two-body interactions (including the interaction of two bodies and the boundary). However, the central sphere interacts simultaneously with all N neighbours and all of these interactions are important
PI. Our model is an excellent approximation at 4 4 10m3. It should also be useful at higher concentrations, but two key simplifications break down fairly quickly as 4 increases. Particle positions are not independent if the placement of one particle interferes with the placement of another. Though this occurs at all values of 4, the effect is negligible at 4 Q 10T3. The approximation in eqn. (23) arises from neglecting higher-order terms. If these are included, the effect of one neighbour is no longer independent of that of another even if their positions are independent. These higher-order terms are negligible at 4~ 10w3. Any empirical distribution of initial positions will depend on the particular mixing process. For large values of 4, the simulation process which we have used would ditfer considerably [26] from the equilibrium distribution used by Ladd [23]. For r#~=lO-‘, the difference is inconsequential. At large values, the equilibrium distribution appears to be a better approximation to real processes. This equilibrium distribution can be simulated fairly easily by treating the spheres as a hard-sphere gas and using the central sphere and spherical boundary as hard surfaces. A randomly chosen configuration from
the equilibrium distribution then serves as the initial distribution of positions from which the moments can be computed numerically from equations containing higher-order terms [2, 231. Instead of applying our model at concentrations where it may still be a reasonable approximation (perhaps 4= 10m2), we have preferred to test it where it is still almost exact, i.e., at += 10B3. The close agreement between theory and simulation provides a secure starting point for simulations at higher concentrations. Our analysis utilizes the symmetry which exists when the test sphere is concentric with the container. This is a condition which applies only at this instant of time. Thus, our analysis cannot be used for studies of the approach to a steady-state distribution [lo, 23,271. This is best achieved by using periodic boundaries [23]. However, this work and earlier studies of the variance [2-6] indicate that the steady-state values obtained will depend critically on the size of the cell used.
Acknowledgements This research was supported in part by the Natural Sciences and Engineering Council of Canada through grants OPG’I268, OPG 5377, and EQP 92826. We thank C.‘ F. Chan Man Fong for his help in the’early stages of this study, and A. J. C. Ladd for permission to use his results prior to their publication.
List of symbols b
.
vector defined, by eqn. (23) ’ polyadic (tensor of rank, n) defined by eqn. (5) Corr(.> correlation coefficient between two components of velocity Cov(u) covariance dyadic whose scalar coefficients form the covariance matrix i; j orthogonal unit vectors in horizontal direction k unit vector in vertical direction excess kurtosis of velocity defined by eqn. (55) Vu(u)’ N number of spheres in test region p(y; Al/) Poisson distribution defined by eqn. (Al) r dimensionless distance from test sphere to another sphere r” ‘unit vector directed from test sphere toward another sphere skewness of velocity defined by eqn. (54) Sk(u) Tr trace (of matrix of coefficients of covariance dyadic) velocity of test sphere U horizontal components of velocity UX,:uy vertical component of velocity UZ b” ,’
169
V
Var(u)
dimensionless volume of test region variance of velocity of test sphere defined sum of variances of velocity components
Greek letters P
9 A Q c1 p w
dimensionless radius of spherical container polar angle of i number of spheres per unit volume volume fraction of solids expected value of u vertical component of ~1 longitude of i
script
&I-) @(*)
expectation conditional expectation terms of specified order
Subsctipts 0 indicating test sphere i indicating ith sphere
i k
12 C. W. J. Beenakker and P. Mazur, ffhys. Fluids, 28 (1985) 3203. 13 H. Grad, Commun. Pure Appl. Math., 2 (1949) 325. 14 E. M. Tory, Chem. Eng Sci., 24 (1969) 1637. 15 J. Frenkel and H. Brenner, J. Fluid Mech., 204 (1989) 97. 16 E. S. Keeping, Introduction to Stattrrical Znf&rence,Van Nostrand, Princeton, 1962. 17 M. H. De Groot, Probability and Statistics, Addison-Wesley, Reading, MA, 1975. 18 R. C. Weast (ed.), Handbook of Tables for prObabil&y and Statistics, Chemical Rubber Co., Cleveland, 1968, p. 7. 19 D. H. Sanders, A. F. Murph and R. J. Eng, Statistics. A Fresh Approach, McGraw-Hill, New York, 1977. 20 L. Brand, Vector Analysis Wiley, New York, 1957. 21 N. Gilbert, Statistics, Saunders, New York, 2nd edn., 1981. 22 W. Feller, An Introduction to Probability Theory and its Applications, Vol. 1, Wiley, New York, 2nd edn., 1957. 23 A. J. C. Ladd, Phys. Fluids A, in press. 24 B. Koglin, Chem. Zng. Tech., 43 (1971) 761. 25 B. W. Char, K. 0. Geddes, G. H. Gonnet, B. L. Leong, M. B. Monagan and S. M. Watt,Maple LanguageReference Manual, Waterloo Maple Publishing, Waterloo, 1991. 26 B. Widom, J. Chem. Phys., 44 (1966) 3888. 27 D. K. Pickard and E. M. Tory,J. Math. Anal. Appl., 60 (1977) 349.
indicating jth sphere indicating kth sphere Appendix
References 1 E. M. Tory and D. K. Pickard, in B. M. Moudgil and P. Somasundaran (eds.), Flocculation, Sedimentation and Consolidation, AIChE, New York, 1986, pp. 297-306. 2 E. M. Tory, M. T. Kamel and C. F. Chan Man Fong, Powder Technof., 73 (1992) 219. 3 R. E. Caflisch and J. H. C. Luke, Phys. Fluids, 28 (1985) 759. 4 E. M. Tory and D. K. Pickard, Powder TechnoL, 47 (1986) 39 and 48 (1986) 189. 5 E. M. Tory and M. T. Kamel, Powder TechnoZ., 55 (1988) 51. 6 E. M. Tory and M. T. Kamel, Powder Technol., 55 (1988) 187. 7 J. Happel and H. Brenner, Low Reynolds Number Hydrodynamics, Martinus Nijhoff, Dordrecht, Holland, 1983. 8 E. M. Tory and D. K. Pickard, Proc. 2nd World Congr. Chem. Eng., Montreal, 1981, Vol. 6, pp. 254-257. 9 G. J. Kynch, Trans. Faraday Sot., 48 (1952) 166. 10 E. M. Tory and D. K. Pickard, Can. 1. Chem. Eng., 55 (1977) 655. 11 D. K. Pickard, E. M. Tory and B. A. Tuckman, Powder Technol., 49 (1987) 227.
Central moments of the Poiwon distribution
In our notation, the Poisson distribution for the number of points (i.e. sphere centres) in the region (of volume v) accessible to them is given by [16]
(Al) The central moments can be calculated from
directly [16]
Alternatively, they may be obtained from the cumulants [16] which are all equal to AV. Then Var(N) = AV,
(A3
8[(iv-
(w
Avy] = hV
and $[(N-h~4]=AV+3(AV)2
W)
Powder Technology, 74 (1993) 159-169
The distribution of’ velocities of a concentric sphere in a dilute dispersion of spheres sedimenting in a spherical container E. M. Tory, M. Bargiel Department of Mathematics and Computer Science, Mount Allison University, Sackville, NB, EOA 3C0 (Canada)
and M. T. Kamel Department of Mathematics, Statistics and Computer Science, Universe
of New Brunswick, Saint John, NB, E2L 4L5 (Canada)
(Received July 20, 1992; in revised form October 19, 1992)
Abstract We consider the motion of identical spheres in an incompressible fluid which is bounded by a spherical container wall. The test sphere is momentarily concentric with the container and the centres of the other spheres are distributed in the accessible region according to a Poisson process. We derive the third and fourth polyadic central moments of velocity and examine their components. These moments, together with the mean and covariance which are known from earlier work, provide a good description of the theoretical initial distribution. Estimates of the hard-sphere equilibrium distribution are obtained from computer simulations using a fixed number of hard spheres. Simulated values agreed closely with the theoretical, indicating that the Poisson assumption, which greatly simplifies the analysis, has little effect on the results. In small containers, the vertical component of velocity is skewed downward and is platykurtic. In large containers, the distributions of all components approach normality.
Introduction It has been generally believed that identical spheres settle with nearly equal velocities aside from a few cases where neighbours are unusually close, but experimental measurements, summarized by Tory and Pickard [l] and Tory et al. [2], show enormous variability. This is supported by theoretical calculations of the variance of particle velocities in dilute dispersions [2-61. The equations of creeping flow imply that the configuration of the system of particles and boundaries determines the distribution of particle velocities [7] and this distribution varies both spatially and temporally [8]. Furthermore, many phenomena which cannot be explained on the basis of Kynch’s theory [9] are accounted for by considering the distribution of particle velocities [l]. Though distributions of instantaneous velocities have not been measured, related measurements are compatible with a distribution which is approximately normal [lo, 111. Following Beenakker and Mazur [12] and Tory and Kamel [6], we consider the motion of identical spheres of unit radius in an incompressible fluid which is bounded by a spherical container wall of radius j3. The test sphere is momentarily concentric with the container
0032-5910/93/$6.00
and the centres of the other spheres are distributed in the accessible region according to a Poisson process with intensity A= 34/(4~) where 4 is the solids fraction [6]. We allow these spheres to overlap each other but not the test sphere or the boundary. In these circumstances, the mobility tensor assumes a particularly simple form [12] and the dimensionless mean and variance are easily calculated [6]. As noted by Tory and Kamel [6], the spherical case .is a useful prototype for more practical studies. Here we extend our theoretical calculations to the third and fourth central moments of velocity, thereby determining the skewness and kurtosis of the distribution. We use computer simulations to approximate the equilibrium hard-sphere distribution and compare the computed moments with the theoretical values. The use of the Poisson distribution leads to a variable number of spheres in the container. This is consistent with the idea of sampling from a much larger volume and leads to unconditional expectations which are simpler than expectations conditioned on a given number, N, of spheres in the volume [2,6]. However, the practical effect is nil because the difference consists entirely of terms which are too small to merit calculation [2].
0 1993 - Elsevier Sequoia. All rights reserved
160
Mean and variance of velocity Integration of eqn. (19) of Tory and Kamel [6] and insertion of the result into their eqn. (18) yields the dimensionless mean velocity 8(u)=
-b,++[2-9p-l+a(~-“)P
(1)
which agrees with their eqn. (20) but, unlike it, contains the term in p-l. To the order of approximation considered by Tory and Kamel [6], this agrees with the result of Beenakker and Mazur [12]. Since b0 is just R plus wall corrections, we are more interested in the second term which represents the mean effect of the other spheres on the test sphere. The variances of the components of velocity are given by [6] Co+)
= ii Var(u,) +J Var(z+,) + kk Var(u,) =
4
Here we are using the convenient notation of Frankel and Brenner [15], but Tory’s definition [14], which differs slightly from theirs, is more suited to our needs. Then [5, 61
[(&B-gj+
(ii+j) G
Here b, represents the contribution of the ith sphere to the translational velocity of the test sphere while b. represents that of the test sphere to its own velocity. Noting [5, 61 that 8(bibj) = [&Y(b)]‘,i #j
(9
and counting the number of terms in each summation, we obtain
/3-‘)kk+a(8-2)].
(2) -AE(b3) -iV(NPolyadic moments of velocity
-N(N-
Since the particle velocity, u, is a vector, its central moments are polyadics [7]. For example, Cov(u) is a dyadic which provides information about the three components u,, u,, and u,. Its diagonality indicates that these components are uncorrelated and the values of the diagonal elements are the variances of the components. Defining the variance of velocity as the scalar quantity [6] Var(u) = Var(u,) + Var(%) + Var(u,)
(3)
implies that
1
(n times)
(5) operator
[14, 151
sum of the distinct permutations
1+3N+3N+3N(N-l)+N+3N(N-l) +N(N-
l)(N-- 2) = (N+ 1)’
(11)
as is obvious from the first line of that equation. It follows that no terms were lost in going from eqn. (8) to eqn. (10). Using eqn. (10) together with [6] $(z@)=
-b,-I@=(b)
(12)
(13) we obtain
~(4w1’lN = S(u’p) - pqup)~(u’pv)~ + 2[qup)]3
g#J -
=
-N{E’(b3) - [‘(b)a(b’)r+
2[8(b)13}
of b, b,. . .b,, . (6)
a[+- c)“l= Gwu
For example, [bIb;~=bIbj:+b2b,b2+b~b,.
The total number of terms in eqn. (8) is therefore
(14)
which is the third central moment given that there are N spheres in the accessible region. The unconditional third central moment is
lw2...ur =
(10)
8(~~~N)=b;+~b,,8(b)~+N~(b~)+N(N-l)[~(b)]~,
z/3z + z p-‘+@(p-“) . (4) =4 [ 112 We need a convenient notation to deal with higher moments because the expressions rapidly become unwieldy. Following Grad [13], we define the n-adic
and the symmetrization
2)[8(b)13
and
Var(u) = Tr[Cov(u)]
b”=bb...b
l)(N-
1)[P(b)8(b2)r
(7)
=-
- Gw
+ %wI - P131W
A@=‘(b3) - ~~(b)~(b’)~ + 2[8(b)13)
161
+ lw(4w’Pr
+ t~(dw - P13}
(15)
This result is obtained by expanding the expression in the first line, using eqn. (14) with 8(N) = AV, and noting that
~~[~-$(~lN)l[~(~lN)-P12~1N)=0
(16)
The fourth to seventh terms in the second part of eqn. (15) can be written as
AV)3][fF(b)]3= - hI$Y(b)]3.
In addition to the values given in the table in Tory and Kamel [6], we require a
2w
(i-k)3i3 sin 8 do de= & ji2k+j2klp + + k3.
(26)
(17)
The final integration (with respect to r) yields Ep(b3). Then eqns. (19) and (24) give (18)
(15), (17) and (18) yield
8[(u - /L)3]= -AK%(P)
(25)
0 0
and the last term (see Appendix) is
Equations
i=isin8cosw+jsinOsinw+kcosO.
G ss
- Plln
= - hVEO[GP(b)kqb2)~ - 3[8(b)13}
-$[(N--
but very tedious. As evident in eqn. (23), b has subunits R and (i.k)i where
3
- Lw4wwwv
S(lMu”ln?
is straightforward
g?(u - 4’1
(19)
which is much simpler than the conditional expectation given by eqn. (14). The techniques used to obtain eqn. (14) yield
(27)
SE+ - ~(#FVl
Writing
= N?(V) - MC?(b)8@3)P + N(N- l)~[~(zP)]2r + (2N-N2)U~(bZ)[~(b)]2~+
u-p=iu,+jU,+k(u,-p)
(3N2-- 6N)[8(b)p
(20) A straightforward but tedious application of the method used to obtain eqn. (19) leads to J%ZJ - d41= AJ@‘(b4) + (W2WV2)l?
(21)
in which the last term consists of three distinct permutations. This is also much simuler than the corresponding conditional moment. 1; deriving eqn. (21), we have used (see Appendix) 8[(N-- AV)‘] = 3(AVj2 + AV.
(22)
(W
and expanding (U-P)’ yields 27 terms. By symmetry, the expected values of twenty are zero. This leaves 8[(u - A31 = lF’WW(uI
- ~11
+ ~2W-~[~;(~z- cl)1+k3WUz - 1-4~1 (29) Equating terms in eqns. (27) and (29), we obtain g?~,z(k - PL)l
= ~b,2(~z- 41 = -6[gln(y)-
z
+ gfi-l]
(30)
and Calculation
Z:[(% - IQ31
of the third central moment
Given [6] that b= ;r-‘[k+(:.k)+ +
Since the third central moment is zero for a normal distribution [16], eqn. (31) is a measure of non-normality.
r-3[3(i4)i-R]
-$8’+;B’n (
+
;
1
~r’(P-‘-p-‘)[2k-(i.k)~,
Calculation (23)
Equation
of the fourth central moment
(25) yields
162 TT 2lT
3
G ss
(i.k)4i4
Equating coefficients in eqns. (33), (34) and (36), and using eqn. (2),
sin 8 dw d6’
0 0
Z@x) - 3[Var(u,)]’
+4+
&
[i2k2+j2kz~+
$
(i4+j4)+
-$-j
[i'j'lp
(32) Using the table in Tory and Kamel[6] together with eqns. (26) and (32), we can complete the first two integrations intheevaluationof8(b4).Thefinalintegrationissimplified by considering only those terms which are a( 1) to @‘(p- ‘). All of the former arise from the lower limit of integration. The final result is AV8(b4)
(( 8601
_#j _ -
10035200
+ SOB-'
= S($) - 3[Var(u,)]’ 8601 6561 P-1 10035200 + 985600
= -34
(37)
and 8[(k -
CL)~IW4412
(38) The expressions on the left-hand side of eqns. (37) and (38) are the fourth ‘cumulants of the distribution of the horizontal and vertical components of velocity [16]. Since these are zero for a normal distribution [ 161,these equations provide measures of non-normality.
(3i4+3j4+[izjq7 )
Covariances
(33) S[(u - p)“] can be evaluated from eqn. (6) of Tory and Kamel and eqns. (2), (21) and (33). The expansion of (u - p)4 contains 81 terms but, by symmetry, the expectations of 60 of these are zero. This leaves
and correlations
The magnitudes of the off-diagonal elements in eqns. (29) and (30) indicate the strengths of the covariances between different components. Since [17]
c0v(u,z&)= iY{[u*”- squx”)](u*- p)} = ~b,z(uz - CL)1
(39)
3
the values of Cov(u,2u,) and Cov($uJ eqn. (30). We also have
are given by
Cov(u,zu,‘) = @ux”4) - a(u,)s($) = +&$) From eqn. (6) of Tory and Kamel [6] and eqn. (2), AW(b2)
=
i2
Var(u,) +j’ Var($) + k2 Var(u,) .
E-
-Var(uJ 8601
4 10035200
Var(u,) +
$$$I-’
(35)
)
(40)
Thus
W?&W2)12~ = (AV)2{ii?(b2)8(b2) + 8[b$(b2)b]
+$(bb ‘bb ‘)}
-4
=
(_El+?!&-
‘I[%
-l*(Y)]}.
(41)
It follows from eqns. (37) and (40) that Var(u$ = Var($) = 8(u,4) - [Var(u,)]’
(36)
2: 2[Var(u,)]’ - 34 I~~~OO (
+ $$OB-l
1 (42)
163
The distribution of velocities
Thus we can calculate [17] ~rrP%
m(k - CL)1 cc)1= [Var@Z) x Var@z)]‘”
-
(43)
since Var(u, - p) = Var(u,). Also Corr(u$Q
=
Cov(u,‘u,“) Var(uz) *
(44
Similarly VarNu, -
Jm4 - Id31=IW(uJ13~
PYI
= 2fvar(uZ)12
-
&&j@, +fi-’
(45) and hence Corr[u,2(u, - p)*] =
cc+:@, - I-4’1 {Var(u,‘) Var[(u, - P)*]}~~
can be calculated.
For f12> 4-l,
CYorr(u$,‘) = =L
(46)
these become
~]u,(UZ - PII Jz Var(uX)[Var(uZ)]‘/2
= -28.17+-1nj?-3n
(51) 1
WY - A41 -3=5.54&-l/3-2. Wr(uz>12
In
(47)
q44) PW41*
-3=
(52)
components,
it is
-0.49C#+p-2.
(53)
Thus, the distributions of u, and u,. are very slightly leptokurtic [19]. Tory and Kamel [5, 61 characterized Var(u) by eqn. (4). While the choice for contraction of 8[(u -P)~] is less obvious, we take the dot product of the first two vectors in each triad. (The symmetry makes the choice immaterial.) This yields an overall measure of skewness
cov(uI”u,“)
16011 XY
2W-(4H2 -8.19x
.
(
For the horizontal
cOrr[u,Z(u, - p)] =
!$A
;pp-lqj-3qn
It follows from eqn. (38) that the distribution of u, is platykurtic [19]. The excess kurtosis [18] is
437337 8800
166133 +di
It follows from eqn. (23) and the symmetry of the expected spatial configuration that the distributions of u, and u,, are symmetrical for all values of 4 and p. Equation (31) indicates that u, is negatively skewed, i.e., toward downward velocities. Using Var(u,) from eqn. (2), we obtain the coefficient of skewness [18] for u,:
10-2+-1~-2
(48)
Sk(u) = -
3f.?
k
p
and Corr[u~(u, - p)“] =
~+a4 - Id21 2 Var(u,) Var(u,)
=0.556~$-‘p-~
56&’ 2=- 9
B-1 (
(49)
Equation (30) indicates that rapidly falling spheres have larger horizontal velocities than those which fall slowly or move upward. While this quantity grows in absolute value as p increases, it is asymptotically zero compared to the variances. As indicated by eqn. (40), the covariance of u,’ and u,” is almost zero. This appears to reflect two opposing effects. It is easily shown that Corr(sin’ w, cos* w) = - 1 .
(#+2p-3/21,
(50)
Hence large values of u,” are associated with small values of u,“. On the other hand, large values of u,’ and t$ are both associated with large values of (uZ- CL)’ as indicated by eqn. (41). It is clear from eqns. (47)-(49) that the correlations among all these velocities decrease rapidly toward zero with increasing p.
Tk
(54) )
This coefficient reflects the fact that all the features are related to the asymmetry in the vertical velocity. An alternative would be to weight the first term in eqn. (27) by the number of permutations, but this is more cumbersome. Of course, no contraction can retain all the information in eqn. (27). Since S[(u - p)*] consists of the pieces given by eqns. (33) and (36), there seems to be no simple formula for summarizing its value. However, we can summarize a measure of excess kurtosis. Equation (36) is closely related to the value subtracted in computing this quantity. This leaves the scalar invariant [20] of AI@(b4). The symmetry of all quantities in eqn. (33) means that we can use any consistent method of computing this invariant. Then [Pj2p contributes (i*i)(j.j) + (j j)(i*i) = 2 with similar contributions from [Pk’p and b”k’p. The result is
164
510151 250880 +p-1
Ku(u) =
other sphere is discarded and new random coordinates are chosen. The probability of placing these N spheres with no overlap is roughly
2 4
= 3.894-‘p -2
(55)
Equations (54) and (55) imply that the skewness and excess kurtosis are arbitrarily small for sufficiently large /3. It is also clear from eqn. (23) that pJn__P(b”)=d(l),
n=5, 6, . . . .
(56)
Consequently, we expect that the distribution of velocities will be approximately normal in containers which are very large compared to the size of the particles. For very dilute dispersions in small containers, the distributions of u, and u,, have slightly smaller tails and/or sharper peaks than the normal distribution [19]. In contrast, the distribution of u, is much more spread out and this persists to fairly large values of @. Simple calculations show that the conditional moments given by eqns. (14) and (20) are virtually identical to the unconditional moments expressed by eqns. (19) and (21) respectively. It is clear from eqn. (6) of Tory and Kamel [6] and eqn. (12) that eqn. (14) may be written as CY{[u --a(u]N)]3)N) = -NtP(b3). Using these same equations shows that
(57) together
Z{[u - S(u]N)]4]nr) -N’&(b4) +N(N-
with eqn. (19)
l)[[EP(b2)]2I”.
(58) Setting N=AV in eqns. (57) and (58) completes the comparison of eqns. (14) and (20) with (19) and (21) respectively.
Computer simulation of velocities The distribution of velocities in small containers can be estimated directly by simulation. Since a comparison between theoretical and empirical results for a Poisson process would serve only to test the accuracy of the random number generator, we replace that approximation by the hard-sphere model. The test sphere is placed at the centre of a spherical container of radius l/3
(59) and N spheres are randomly placed in the accessible region. Any attempted placement which overlaps an-
(l- $(l- $)...(l-y 4) =exp[ - (N;l)+l which corresponds to an approximate solution of the well-known birthday problem [22] with N people and a year of N/4 days. For += 10m3 and N= 1000, 2000, 5000 and 10 000, these probabilities are roughly 0.607, 0.368, 0.082 and 0.007 respectively. Thus, it is almost inevitable that some overlaps will occur in many of the simulations and hence that the overall hard-sphere distribution will differ slightly from the Poisson point approximation. However, these differences are important only insofar as they limit the number of spheres which are close to the central sphere. For $J= 10e3, the distribution of centres is very close to the equilibrium distribution of hard spheres. Thus, the distribution of velocities is correspondingly close. For this reason and because the simplifying assumptions which permit the derivation of theoretical values rapidly become untenable with increasing 4, we set += 10v3 for all simulations. (See Discussion.) Each data set consisted of 10 000 simulations. We obtained 15 sets for N= 1000 and 10 sets for N = 2000, 5000 or 10 000. Since theoretical values for the moments of u, and u,, are identical, and simulated values can differ only by chance fluctuations, we first compute the moments separately for each set, check that there are no significant differences, and then pool the values for u, and u,, On the graphs, they are all labelled as moments of u,. Figures l(a) and l(b) show the empirical distribution of u, and uy (combined) for N= 1000 and 10000, respectively for all data sets for these values of N. The symmetry is evident in both, while the greater spread is apparent in Fig. l(b). Figures 2(a) and 2(b) show the distribution of u, + b, for 1000 and 10 000 neighbours respectively. The skewness is clear in both and, once again, the greater spread is apparent. Values for N = 2000 and 5000 (not shown) fall nicely between those for 1000 and 10 000 in Figs. 1 and 2. For each configuration, we evaluated 1~+b, the effect of the other spheres on the central sphere. For each set of 10 000 configurations, we calculated estimates of @@)+b, and the second, third and fourth (conditional) central moments of the components of u. The sample means (i.e., the means for each set) are approximately normally distributed [19], so we used the t-distribution (with 14 or 9 degrees of freedom) to
165
1000 spheres
1000 spheres
c
0.0
(4
2,
+ bo
10000 spheres
10000 spheres
1
(“1
1
(b)
u,+
bo
Fig. 1. Frequency distribution of the horizontal components of velocity. Values for u, and u,. have been pooled. In all Figs. in this paper, += 10w3. (a) 1000 neighbours, (b) 10 000 neighbours.
Fig. 2. Frequency distribution of the vertical component of velocity. The constant, -b, representing the velocity of the central sphere in the absence of any neighbours, has been subtracted froin ui to highlight the effect of neighbours. (a) 1000 neighbours, (b) 10 000 neighbours.
calculate 95% confidence intervals for the overall means of each of these quantities [21]. Figure 3(a) shows the simulated values of $(u,) and 8(u,J, combined, while 3(b) shows those for a@,)+&. The central moments were calculated in the usual way [16, 181. Though we reduced the degrees of freedom by one to allow for the calculation of the mean from the data, this had no practical effect since each moment was based on 10 000 simulations. As noted above, the conditional moments differ negligibly from the unconditional. In each case, the simulated moments (taken to represent both the conditional and unconditional values) are compared to theoretical values for the unconditional moments. Figure 4(a) displays the values of Var(u,) and Var(u,), combined, while 4(b) and 4(c) show the values of Var(u,) and Var(u), respectively. Figure 5 shows the third moment. Figure 6(a) and (b) indicates the fourth cumulants of u, and uu (combined) and U, respectively.
In all cases but one, the confidence limits for a fixed number of hard spheres include or almost include the theoretical values for a Poisson process. In this one case, the confidence limits for the fourth cumulant of U, and u,, (combined) are very close together and the actual difference in values is extremely small. The small differences between theoretical and simulated values indicate that the simulation provides an empirical distribution of velocities which closely approximates that corresponding to a Poisson spatial distribution of centres.
Comparison
with other results
We have already noted [2] that our eqn. (4) is consistent with Caflisch and Luke’s [3] eqn. (3) which indicates that the variance of velocity is approximately proportional to
166
~,,,,,,,,,,,, -0.003
50
I ,,,,,, I ,,,,,,,,,,,,/,),,/,
0
100
150
200
0.00
250
B
(a)
0"
100
150
200
2
200
2 0
200
; ‘0
B
(a)
0.006
+
T--. N
F-Y N 3 z
‘~“~~“~1~~‘~‘(~“~“(~~‘~“I”~~~‘~‘~
50
ikO6
0.001
I)
0
-0.004
50
0
~"~""~I~"~~"~'I~'~"~~'~I~~~'~"" 100 150 200
@)
0.02
250
B
50
/“~~‘~“‘I~‘(~“~~‘I”~~‘~”
100
150 B
(b)
Fig. 3. Mean value of components of velocity (a) Horizontal. By symmetry, this value is exactly zero for both the hard-sphere and Poisson distribution of centres. Simulated values of u, and uy have been pooled. (b) Vertical. Tpe ordinate shows the effect of all neighbours on the velocity of the central sphere. By p = 225, the theoretical value is within 10v6 of the horizontal asymptote of 0.002.
NID$P’ = &#I .
(61)
Very recently, Ladd [23] calculated Cov(u) for spheres sedimenting in a square box with periodic boundary conditions. In our notation, his result for large N is I
Cov(u)=N*“#P[O.O66
2(ii +j) 4 0.823k.k] .
(62)
Since the numerical coefficients are different for different boundary conditions, we cannot expect quantitative agreement between eqns. (2) and (62). However, the coefficients are of comparable magnitude, viz., 0.072 3 versus 0.066 2 and 0.579 versus 0.823. Here we have used eqn. (61) and approximated our variance by the leading terms. Ladd’s Fig. 5 [23] shows Var(u,]N) 20.012 5 and Var(u,(N) = 0.153 for simulations with N = 108, 4 = 0.05, and a random distribution of sphere centres. (These
,
0.0 Cc)
50
~“~~“‘~I’~‘~“‘~~I”‘~J’~‘~I’~““”
100
150 P
Fig. 4. Variance of velocity. (a) Horizontal components. Values of u, and uYhave been pooled. (b) Vertical component. (c) Var(u) as defined by eqn. (3).
values are consistent with those shown in his Table 2 and Fig. 5 for steady-state values.) Insertion of this N and 4 in our eqns. (59) and (2) yields Var(u,)=0.0314 and Var(u,) = 0.175. (Though the equations for Var(u,) and Var(u,lN) differ, the correction factor [2] is neg-
167
2
0
z 0.00000
2 , -0.006 xN kc Q
I
-0.008
-O.O1 O5rn
250 B
Fig. 5. Third moment of velocity. The solid line and the small circles refer to the third central moment of u.. The dotted line gives the values of [fm(y)-
$$+
-$$$-t]?
which
is used in eqn. (54) to calculate the skewness of U.
ligiile even for N- 108.) Given the difference in boundary conditions and the fact that 4-0.05 is well beyond the range where our assumptions are valid, the orderof-magnitude agreement is the best that can be expected. There. are no experimental data which are directly comparable. with our theoretical results for Var(u,). However, an extensive analysis [2] of published data showed qualitative agreement between theory and experiment. Little work has been reported on the higher moments of velocity. Recently, however, Ladd [23] obtained values of the fourth cumulants of U, and U, at several concentrations. Dividing our eqns. (37) and (38) by 3 Var(u,) and 3 Var(u,) respectively yields the cumulants plotted in his Fig. 8. For r$= 0.05, our values, 0.068 and 0.669 respectively, are much larger than his because the Poisson approximation is then grossly inadequate, permitting too many close neighbours and thereby inflating the value of JY[(u-p)‘]. Nevertheless, Ladd’s results complement ours by extending the computed near-normality to dispersions with fewer particles and higher concentrations, viz., 108 spheres at 4=0.05: At moderate to large concentrations (0.25 Q 4 Q 0.50), nonGaussian effects were substantial, especially for the initial (equilibrium) distributions [23]. ,According to Caflisch and Luke [3], it can be shown that the distribution of NW% approaches a Gaussian as N? 03 with 4 held constant. Since Var(N- ‘%) = N - In Var(u)
(63)
the variance of N - “‘u, which is approximately proportional to &“, is finite in the limit of an infinite volume. It follows, of course, that the distribution of
0.000 @I
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,I,,:,,,,,,
50
100
150
200
250
B
Fig. 6. Fourth cumuhmts of components ofvelocity. (a) Horizontal. The theoretical curve, given by eqn. (37), has a horizontal asymptote of roughly - 2.57 X 10e6. (b) Vertical. The theoretical curve, given by eqn. (38), has a horizontal asymptote of roughly 1.85~10-~ which is approached very slowly on account of the logarithmic term.
u is approximately normal in all containers which are sufficiently large in every dimension.
Discussion Near-normality of the steady-state distribution of velocities implies that the distribution of distances travelled in an interval of time is approximately normal [lo] and the distribution of times to traverse a fixed distance is approximately log-normal [ll]. Thus, the simulated near-normality of steady-state distributions [23] is consistent with the experimental results of Tory and Pickard [lo] and Koglin [24]. The near-normality of both the initial and steady-state distributions is consistent with the three-parameter Markov model of Tory and Pickard [lo, 111 which postulates a onedimensional Langevin equation for u,.
168
An almost-normal distribution is well characterized by its mean and its second, third and fourth central moments. Thus, we made no effort to find higher moments. These would be prohibitively time-consuming for hand calculation, but might be accessible through symbolic manipulation routines. We used Maple [25] on a SUN SPARCstation 1+ to check the calculations for g’[(u -p>“], n =2, 3, 4. Taking i, j, k as scalar variables collapses the terms in [. B and simplifies the symbolic computations. Nevertheless, memory limitations made it necessary to compute the fourth moment in pieces. The symbolic computation gives all the terms for B(V), but the value of b expressed by eqn. (23) is only an approximation. Thus, we discarded all o@ -‘) terms. Higher moments could probably be calculated by pruning the expression to eliminate all terms which would be insignificant and by breaking down the integrals into several pieces. Despite the small values of p, velocities vary substantially. The computational effort for our simulations increases as /3’. Fortunately, the curves are already close to normal for N= 10 000 (/3=216) and the theoretical values of the moments indicate the approach to normality with increasing p. In our model, the effect of each neighbour is independent of that of every other neighbour. That is, we consider only two-body interactions (including the interaction of two bodies and the boundary). However, the central sphere interacts simultaneously with all N neighbours and all of these interactions are important
PI. Our model is an excellent approximation at 4 4 10m3. It should also be useful at higher concentrations, but two key simplifications break down fairly quickly as 4 increases. Particle positions are not independent if the placement of one particle interferes with the placement of another. Though this occurs at all values of 4, the effect is negligible at 4 Q 10T3. The approximation in eqn. (23) arises from neglecting higher-order terms. If these are included, the effect of one neighbour is no longer independent of that of another even if their positions are independent. These higher-order terms are negligible at 4~ 10w3. Any empirical distribution of initial positions will depend on the particular mixing process. For large values of 4, the simulation process which we have used would ditfer considerably [26] from the equilibrium distribution used by Ladd [23]. For r#~=lO-‘, the difference is inconsequential. At large values, the equilibrium distribution appears to be a better approximation to real processes. This equilibrium distribution can be simulated fairly easily by treating the spheres as a hard-sphere gas and using the central sphere and spherical boundary as hard surfaces. A randomly chosen configuration from
the equilibrium distribution then serves as the initial distribution of positions from which the moments can be computed numerically from equations containing higher-order terms [2, 231. Instead of applying our model at concentrations where it may still be a reasonable approximation (perhaps 4= 10m2), we have preferred to test it where it is still almost exact, i.e., at += 10B3. The close agreement between theory and simulation provides a secure starting point for simulations at higher concentrations. Our analysis utilizes the symmetry which exists when the test sphere is concentric with the container. This is a condition which applies only at this instant of time. Thus, our analysis cannot be used for studies of the approach to a steady-state distribution [lo, 23,271. This is best achieved by using periodic boundaries [23]. However, this work and earlier studies of the variance [2-6] indicate that the steady-state values obtained will depend critically on the size of the cell used.
Acknowledgements This research was supported in part by the Natural Sciences and Engineering Council of Canada through grants OPG’I268, OPG 5377, and EQP 92826. We thank C.‘ F. Chan Man Fong for his help in the’early stages of this study, and A. J. C. Ladd for permission to use his results prior to their publication.
List of symbols b
.
vector defined, by eqn. (23) ’ polyadic (tensor of rank, n) defined by eqn. (5) Corr(.> correlation coefficient between two components of velocity Cov(u) covariance dyadic whose scalar coefficients form the covariance matrix i; j orthogonal unit vectors in horizontal direction k unit vector in vertical direction excess kurtosis of velocity defined by eqn. (55) Vu(u)’ N number of spheres in test region p(y; Al/) Poisson distribution defined by eqn. (Al) r dimensionless distance from test sphere to another sphere r” ‘unit vector directed from test sphere toward another sphere skewness of velocity defined by eqn. (54) Sk(u) Tr trace (of matrix of coefficients of covariance dyadic) velocity of test sphere U horizontal components of velocity UX,:uy vertical component of velocity UZ b” ,’
169
V
Var(u)
dimensionless volume of test region variance of velocity of test sphere defined sum of variances of velocity components
Greek letters P
9 A Q c1 p w
dimensionless radius of spherical container polar angle of i number of spheres per unit volume volume fraction of solids expected value of u vertical component of ~1 longitude of i
script
&I-) @(*)
expectation conditional expectation terms of specified order
Subsctipts 0 indicating test sphere i indicating ith sphere
i k
12 C. W. J. Beenakker and P. Mazur, ffhys. Fluids, 28 (1985) 3203. 13 H. Grad, Commun. Pure Appl. Math., 2 (1949) 325. 14 E. M. Tory, Chem. Eng Sci., 24 (1969) 1637. 15 J. Frenkel and H. Brenner, J. Fluid Mech., 204 (1989) 97. 16 E. S. Keeping, Introduction to Stattrrical Znf&rence,Van Nostrand, Princeton, 1962. 17 M. H. De Groot, Probability and Statistics, Addison-Wesley, Reading, MA, 1975. 18 R. C. Weast (ed.), Handbook of Tables for prObabil&y and Statistics, Chemical Rubber Co., Cleveland, 1968, p. 7. 19 D. H. Sanders, A. F. Murph and R. J. Eng, Statistics. A Fresh Approach, McGraw-Hill, New York, 1977. 20 L. Brand, Vector Analysis Wiley, New York, 1957. 21 N. Gilbert, Statistics, Saunders, New York, 2nd edn., 1981. 22 W. Feller, An Introduction to Probability Theory and its Applications, Vol. 1, Wiley, New York, 2nd edn., 1957. 23 A. J. C. Ladd, Phys. Fluids A, in press. 24 B. Koglin, Chem. Zng. Tech., 43 (1971) 761. 25 B. W. Char, K. 0. Geddes, G. H. Gonnet, B. L. Leong, M. B. Monagan and S. M. Watt,Maple LanguageReference Manual, Waterloo Maple Publishing, Waterloo, 1991. 26 B. Widom, J. Chem. Phys., 44 (1966) 3888. 27 D. K. Pickard and E. M. Tory,J. Math. Anal. Appl., 60 (1977) 349.
indicating jth sphere indicating kth sphere Appendix
References 1 E. M. Tory and D. K. Pickard, in B. M. Moudgil and P. Somasundaran (eds.), Flocculation, Sedimentation and Consolidation, AIChE, New York, 1986, pp. 297-306. 2 E. M. Tory, M. T. Kamel and C. F. Chan Man Fong, Powder Technof., 73 (1992) 219. 3 R. E. Caflisch and J. H. C. Luke, Phys. Fluids, 28 (1985) 759. 4 E. M. Tory and D. K. Pickard, Powder TechnoL, 47 (1986) 39 and 48 (1986) 189. 5 E. M. Tory and M. T. Kamel, Powder TechnoZ., 55 (1988) 51. 6 E. M. Tory and M. T. Kamel, Powder Technol., 55 (1988) 187. 7 J. Happel and H. Brenner, Low Reynolds Number Hydrodynamics, Martinus Nijhoff, Dordrecht, Holland, 1983. 8 E. M. Tory and D. K. Pickard, Proc. 2nd World Congr. Chem. Eng., Montreal, 1981, Vol. 6, pp. 254-257. 9 G. J. Kynch, Trans. Faraday Sot., 48 (1952) 166. 10 E. M. Tory and D. K. Pickard, Can. 1. Chem. Eng., 55 (1977) 655. 11 D. K. Pickard, E. M. Tory and B. A. Tuckman, Powder Technol., 49 (1987) 227.
Central moments of the Poiwon distribution
In our notation, the Poisson distribution for the number of points (i.e. sphere centres) in the region (of volume v) accessible to them is given by [16]
(Al) The central moments can be calculated from
directly [16]
Alternatively, they may be obtained from the cumulants [16] which are all equal to AV. Then Var(N) = AV,
(A3
8[(iv-
(w
Avy] = hV
and $[(N-h~4]=AV+3(AV)2
W)