A three-parameter markov model for sedimentation II. Simulation of transit times and comparison with experimental results

Powder Technology,

221

49 (1987) 227 - 240

A Three-Parameter Markov Model for Sedimentation II. Simulation of Transit Times and Comparison with Experimental Results D. K. PICKARD* Department

of Statistics, Harvard University, Cambridge, MA 02138 (U.S.A.)

E. M. TORY Department

of Mathematics and Computer Science, Mount Allison University, Sackville, N.B. EOA 3C0 (Canada)

and B. A. TUCKMAN** Department

of Statistics, Harvard University, Cambridge, MA 02138 (U.S.A.)

(Received June 19,1986)

SUMMARY

An adequate description of the global behavior of sedimenting suspensions requires a knowledge of not only the mean velocity, but also the variation and autocorrelation of the velocities of individual particles. Much of the published information on these two parameters is encoded in transit times. For example, Koglin measured the length of time for individual spheres to traverse a fixed distance x and found that sedimentation velocities based on these transit times satisfied a log-normal distribution. Computer simulations, which yield the distribution of transit times as a function of the three parameters of the parent distribution, are an effective method of decoding these data. Our three-pammeter Markov model yields first-passage times which are approximately log-normal, and the convergence to log-normality as x increases is remarkably fast. Thus, Koglin’s data provide strong support for the Markov model, experimentally verifying theoretical predictions based (via simulations) on it. INDIVIDUAL NOMENA

VARIATION

AND GLOBAL

PHE-

Theoretical and experimental studies of sedimentation generally assume that there is *Recent address: Department of Mathematics and Statistics, Queen’s University, Kingston, Ont. K7L 3N6 (Canda). Deceased July 28,1986. **Present address: Sloan School of Management, Massachusetts Institute of Technology, Cambridge, MA 02139 (U.S.A.) 0032-5910/87/$3.50

a unique velocity which is representative of all particles at that concentration. Flux theory is based on this assumption [l], and even those theoretical treatments which recognize the variability of particle velocities simply use it to compute a mean velocity. In reality, sedimentation is an enormously complicated system of interacting particle paths [2, 31. The remarkable spatial and temporal variability of particle velocities in approximately monodisperse suspensions of spheres [ 4 - 71 suggested to us that existing models were not rich enough to handle this complexity. In 1967, therefore, we began work on a stochastic approach which explicitly recognized and dealt with this diversity. This approach succeeds because it effectively models the deterministic chaos of the system of differential equations governing particle velocities [8] and provides a simple mechanism for treating individual trajectories [9]. The Markov model which resulted [lo, 111 was formalized [12 - 141 and used to explain experimental results [ 3, 151. It is precisely the behavior of individual particles which is at variance with standard theoretical and experimental results. Data on particle trajectories have been largely ignored, yet they contain far more information than the interface velocities on which attention has been lavished. Indeed, we have shown [14] that the interface velocity of a monodisperse suspension need not correspond to the mean velocity of particles at the initial concentration. This result, which is incomprehensible in the context of flux theory [3], arises because an adequate description of global behavior requires the incorporation of @ Elsevier Sequoia/Printed in The Netherlands

228

individual variation and autocorrelation. Other phenomena which arise from the significant variability of velocity include [ 151 the depletion of upper levels in sedimentation [ 161, the diffuse interface at low concentrations [17], and the effect of concentration gradients on mean settling velocities [ 181.

EXPERIMENTAL

RESULTS

In contrast to the general Markov model [12 - 141, which provides a framework within which phenomena can be understood [ 151, the three-parameter model [ 171 makes specific predictions about the distribution and autocorrelation of velocities. Quantitative tests of these predictions are difficult because the data base is woefully inadequate. To the best of our knowledge, only four major studies [4 - 7,16,19 - 211 have examined the behavior of individual particles. The first three measured the times to traverse a fixed distance. Independently of our work and prior to its publication, Koglin carried out an experimental study of transit times in dilute monodisperse suspensions [ 19,201. Solids concentrations ranged from 8.9 X 10m5to 2.6 X 10e2 by volume [19 (Fig. 4), 20 (Fig. l)]. Using

(1) and estimating its mean and variance in the usual way,

(2) he concluded (using the KolmogorovSmirnov statistic [ 221) that his values of Y were consistent with a normal distribution. He also found that the Y-values for particles traversing a sequence of distance sections were correlated [19]. The correlation decreased exponentially as the distance between intervals increased [19]. Furthermore, the distribution of Y evolved from an initial to a steady state. The transit times in these studies [4 - 7,19, 201 vividly illustrate the diversity of individual velocities. In Koglin’s work, the

existence of a distribution with a substantial variance, the correlation between mean velocities of a particle over successive distance intervals, and the existence of a steady state toward which initial velocity distributions evolve, all strongly support the Markov interpretation [lo - 15,171. However, the theoretical prediction of transit times proved to be an intractable problem, so it was not clear whether the log-normality of the distribution of Koglin’s transit times was evidence for or against the three-parameter model [17]. Computer simulations of transit experiments are an effective method of deciphering these data, testing the theory, and extracting the three parameters [ 171 of the parent velocity distribution.

THE MARKOV

MODEL

Brenner [2] has shown that configuration, rather than concentration, determines particle velocities at low Reynolds numbers. We contend that configuration acts primarily through concentration and view remaining effects as contributing velocity variation. Thus, we replace the apparent chaos of intractable complexity [8,9] with the formal chaos of stochastic models parametrized by concentration. This decisive step leads to a new simplicity in which particle velocities at a given concentration are treated as independent observations of the same stochastic process; i.e., conditioning on concentration decouples the spatial interaction. Velocities are described by Markov processes, the simplest structures incorporating velocity continuity and autocorrelation [23]. For particles at concentration 8, velocity transitions from u to u, during time increments t, are governed by (conditional) probability density functions (pdfs), he(u, t I u). Slurry evolution then proceeds incrementally as follows: particles consult their current concentrations and independently choose their new velocities using the relevant transition pdfs, generating new positions and hence new concentrations. When concentration remains relatively constant, the Markov processes converge to steady state [12]; i.e.,

Mu, tlu)

uniformly -

atI

(3)

229

Thus, the Markov processes are ergodic

[231.

Applications of the Markov model require the specification of the transition functions he for all initial velocities, time increments, and concentration parameters. Since direct estimation of these would require inconceivable amounts of data, parametric models are needed. Diffusion processes provide a rich and natural class of candidates.

MONTE CARLO SIMULATION

An attractive and powerful feature of the Markov model is the ease with which it may be simulated; in essence, the slurry evolution described above is a simulation algorithm. During small time increments r, particle accelerations from u to ZJmay be treated as constant at (LJ- u)/T, so the position and velocity of each particle at time t + T may be deduced from the corresponding values at time t, once the new velocities are specified. These velocities are governed by the transition pdfs he(u, ~1u) where the applicable concen-

tration 8 is determined by the configuration of all particles at time t. Monte Carlo simulation of entire slurries yielding complete trajectories for all particles is therefore conceptually straightforward. Practically, the limiting computational factor is the particle-by-particle concentration updating. By restricting attention to the vertical components of position and velocity vectors, this bottleneck can be eliminated. Concentration updating then reduces to convolution which submits to the fast Fourier transform [ 241. Consequently, systems involving several thousand particles can be efficiently simulated [9]. Figures 1 and 2 illustrate such simulations on an Appolo Domain (DN600) high-resolution graphics work station. Approximately one thousand particles are displayed as they settle toward the packed bed. Concentration, velocity, and flux profiles are shown alongside the sedimenting particles. The onedimensional trajectories are smeared across the horizontal dimension and allowed to overlap, thereby implying a third dimension. Visual realism is enhanced by the use of a separate zero-mean velocity process (which need not

b Fig. 1. Simulation of the sedimentation of a monodisperse suspension. Initial uniform concentration: 15% solids by volume. Elapsed time: 200 dimensionless units. The upper (lower) line indicates the starting (finishing) line for transits of the distance x.

230

PROFILES

Fig. 2. Simulation of the sedimentation of a dilute monodisperse suspension. Initial uniform concentration: 2% solids by volume. Elapsed time: 50 dimensionless units. Note the fuzzy interface and the depletion of the upper levels.

depend on concentration) to model horizontal motion. By automating frame-by-frame photography, simulated real-time motion pictures of sedimentation have been obtained. Such films will be useful in assessing complex slurry phenomena. Simulation of steady state is particularly simple and informative. Concentration and hence the transition pdfs remain fixed, so it suffices to characterize the trajectory of a single generic particle; i.e., system size is irrelevant. Our present concern is the simulation and analysis of transit phenomena in steady state. The positions X(t) and velocities V(s) of our generic particle then clearly satisfy X(t) = X(0) + IV(s)

ds

(4)

0

Since much of the data on individual trajectories takes the form of transit or first-passage times, let T(x) denote the time required by our generic particle to traverse the fixed distance x. Referring to Figs. 1 and 2, T(x) is the time required to traverse the distance

between the horizontal lines. When negative velocities are possible, the ambiguity is resolved by considering the first downward crossing at each line. Consequently, T(x) = inf[t > 0:X(t) =X(O) + x]

(5)

and hence X=

J 0

T(x) V(s) ds

(6)

In our stochastic model, the velocities V(s) are random variables, so some subtle technical issues [25] arise in considering X(t) and T(x) via eqns. (4) and (5). Discussion of individual trajectories necessitates a pathwise interpretation for the stochastic integral in eqn. (4). Consequently, we must require the Markov velocity process [ V(s):s 2 0] to have a version [25] with continuous sample paths [25]. This assumption also guarantees that X(t) and T(x), as defined above, are bona fide random variables (i.e., may be studied stochastically [25]). More generally, under this assumption, analytic tools may be applied within realizations, and stochastic

231

tools across them, without caveats on the interpretation and/or validity of the results. Notice that the velocity process [V(t)] is assumed to be Markov, so the resulting position process [X(t)] is not [25]. However, by eqn. (4), the vector process describing the joint values of position and velocity is Markov [25]. The pathwise equation t+7 X( t + T) = X(t) + .j- V(s) ds t

(7)

provides the basis for simulating particle trajectories. Assuming constant acceleration during the short time interval [t, t + T] yields X(s) at all interim times in terms of V(t), V(t + T), and X(t). Hence, the entire trajectory may be deduced from the velocities V(m), which may be easily generated from the transition pdfs he(u, ~1u). Treating acceleration as constant in eqn. (7) is equivalent to approximating the integral by the trapezoidal rule. One might hope to improve the accuracy by using more sophisticated numerical procedures. However, higher-order methods make neither analytic nor stochastic sense. Such methods depend for their enhanced accuracy on the smoothness of the integrand. In Markov diffusion processes, the infinitesimal variance is always 0(r), so

vtt +7) -

var

[

v(t)

1V(t) = u

7

1=

0(7-l)

(8)

and hence [ V(s): s > 0] is nowhere differentiable [25]. Furthermore, the joint position and velocity process is Markov, so predictorcorrector methods, which use more than the most recent observations, are stochastically inappropriate. It may also be possible to compute the joint Markov transition structure, X(t + T), V(t + 7)1X(t)

=x,

V(t) = u

(9)

In this case, the numerical integration of velocity to obtain position is unnecessary. Position-velocity skeletons [X(m), V(w): n > 0] can then be simulated directly.

THE THREE-PARAMETER

MODEL

The simplest velocity model which incorporates the basic features of sedimentation is a variant of the Ornstein-Uhlenbeck process [25]. Velocity is given by the stochastic differential equation dV=-_P(V--)dt+odW

(19)

where [W(t)] is a standardized Wiener process [25]. Basically, this states that velocities are drawn toward the value /.Iin proportion to their distances from p, and that random noise is superimposed. The resulting velocity process is an ergodic Markov diffusion [17, 251. Its transition structure

(11)

V(t + T)l V(t) = u

is Gaussian with mean m(v, 7) = p + (u - p)e+

(12)

and variance v(u, 7) =

o*( 1 - e-2p7)

20 The steady-state distribution is V-N i

p, $

(13)

(14)

1

and the autocorrelation of velocities is p(t) = epPt

(15)

This simple model was introduced in [ 171 which is considered to be I. By eqn. (4), the joint position-velocity process is both Markov and Gaussian. Indeed, the transition structure, eqn. (9), is given by the means (gX(t

+ 7)1X(t)

+ (1 - e+‘)(u

=x,

V(t) = u] =x + p7

- j.4) (16)

P e[v(t

+ 7)1X(t) =x, V(t) = u] = m(u, 7)

(17)

the variances Var[X(t

+ @IX(t)

=x,

V(t) = u]

= (207 - 3 + 4e+* - e-2pT)u2

2P3 var[v(t

+ 7)1X(t)

=x,

(18)

V(t) = u] = v(u, 7) (19)

232

structure of the model and clarifies the nature of its parameters. Let

and covariance

u u, = P

(20)

t, = pt

An important feature of this model is its ability to use experimental data; only the form of the governing pdfs is specified, leaving the three parameters to be determined experimentally. It has been shown that, in dilute suspensions, the particle/container diameter ratio (d/D) critically affects the level and variety of particle velocities [ 191. Consequently, each particle/container system generates a family of Markov processes parametrized by concentration; i.e., p, u, and fl are experimentally determined functions of the concentration parameter 0. Notice that it is not necessary to restrict this model to purely hydrodynamic particle interactions; other effects, such as electrostatic, could be incorporated directly via the parameters. In steady state, sequentially observed positions of a single particle at fixed times inherit multivariate normality from the velocity process. Moreover, particles act independently. Consequently, maximum likelihood estimation for p, a2, and /3is effective and straightforward [17]. In contrast, the distributional form for first-passage data is unknown and almost certainly complex; so this estimation problem is currently intractable. Analyses and comparisons of data are further complicated by the diversity of experiments, the difficulty in assessing steady state, and the use of suspensions which are not strictly monodisperse. Let u0 denote the Stokes velocity of the spheres. Taking

(21)

Px

x* = P

whence El* =

P* = 1

o/m

o*=-=

x

(23)

P

h, the coefficient of variation, is the only remaining parameter, the other two being absorbed in the scaling process. t, is measured in units of correlation intervals and 3t, in correlation lengths. This simplification shows that the coefficient of variation h is the fundamental parameter which determines the nature of sedimentgtion in steady state; the other two parameters (0 and p) simply determine the time and distance over which the phenomena occur. This is the key which enables us to decipher first-passage results.

FIRST-PASSAGE

TIMES

When velocities must be positive, firstpassage probabilities may be expressed in terms of the position process [X(t)] . Indeed, T(x) > t * X(t) < x

(24)

so PET(x) > t] = P[X(t)

< 3t]

(25)

i.e., the distribution functions of Z’(x) and X(t) sum to unity, the parameter in one being the variable in the other. In our three-parameter model, the velocity process [V(t)] (and hence the position process [X(t)]) is Gaussian, so in dimensionless form, X(t) is normally distributed with cx(t)

yields a dimensionless form [3] in which parametric values from different experiments can be compared. Although this dimensionless version is useful, a different standardization crystallizes the fundamental

I

(22)

= t

Var X(t) = 2h2(e+ - 1 + t) I

(26)

provided that X(0) = 0. When X, the coefficient of variation, is small, negative velocities are rare, so the distribution of T(x) is given approximately by

233

log t-ww

log l-wwo

Fig. 3. Theoretical approximation

P[T(x) > t] =

(

cp

of the pdf for log,[T(x)/x]

when the coefficient of variation h is small.

x-t

XJ2(e8

-

1 + t) 1

(27)

where Q is the cumulative distribution function for the standard normal. We fix h and CC and, from eqn. (27), calculate P[T(x)/x > t/x] for t > 0. From this, we obtain the pdf for log,[T(x)/x]. Figure 3 illustrates these pdfs for several values of x when h = 0.1 and X = 0.3. Notice that, for h = 0.1, the curves are very symmetrical and almost perfectly centered about log[T@)/3c] = 0. As X increases, the curves remain highly symmetrical, but the initial centering is much worse. This is due to slight changes in the symmetry in the tails. For larger values of h, negative velocities occur more frequently. Then, a particle can reach distance x and subsequently retreat, leaving X(t) < x but T(x) < t, thereby contradicting eqn. (24). The relationships between X(t) and T(x) are still important but less easy to analyse. In practice, the concentration range where X is large enough to generate a significant proportion of negative velocities is about 0.5% through 9% (see Fig. 4).

Fig. 4. Fraction of solids moving upward versus local solids concentration.

state velocity is unity, time is suitably scaled, and the fundamental parameter is the coefficient of variation X. Particles are started at the origin and ultimately pass milestones located at integral distance x (X = 1, 2, . . ., 16). Time increments were small enough that 20 were generally needed to traverse the distance between consecutive milestones. The simulation discretizes the Markov model and treats acceleration as constant during the small time increments 7. That is, AX is approximated by &

STEADY-STATE

SIMULATIONS

For simulations, we use the dimensionless form discussed above, so the mean steady-

=

[V(O) + V(T)lT 2

(28)

For this approximation, the global error in locating particles is 6(r2). Note, however, that

234

W&d

= C{& V-%)1VW13

&R(x) = =

which is exact _[26]. Furthermore, the variability of AX and its co-variability with V are remarkably close to the exact values (6(r3) and 6(~~), respectively [26]). When X[(n - l)r] x, a quadratic interpolation scheme should be used to estimate the exact crossing time (see the Appendix). When particles can have negative velocities, a crossing may occur between (n - 1)~ and nr even though X[ (n l)r], X(nr) < x. This is likely, in fact, when V(m) is negative and X(nr) is close to x. The quadratic interpolation also accommodates this situation. It cannot describe two (or more) changes in direction during a single time increment, but this event occurs rarely unless both the coefficient of variation and the time increment are large. To simulate transit times in steady state, it is necessary to start the particles off with the correct velocity distribution. Let S(x) denote the velocity of a particle at the (random) time T(x); i.e., S(x) is its speed as it first reaches the milestone x. Then S(0) is the appropriate initial distribution of velocities. Of course, in steady state, the distribution of S(x) is independent of x. Surprisingly, S is not the steady-state velocity distribution; conditioning on position, rather than time, yields a different result. Let R(x) denote the (positive) velocity of a random particle arriving at x. To an observer at x, fast and slow particles arrive in proportion to their speeds; i.e., rqdr)

Pd[R(x) = r] =

(30)

J we(u) du

where qe is the steady-state pdf of V. Then,

sr2qdr)dr O s dr

=l+A2

(31)

rqdr)

0

When there are no negative velocities, particles pass each milestone precisely once, and hence S(x) = R(x). In this case,

(32)

In our Gaussian model, negative velocities are always possible, so particles can pass the same milestone many times. If the trajectory crosses 0 and x only once, then S(x) = R(x); if it crosses either barrier more than once, S(x) uses the first values while R(x) weights all of them equally. Thus, S(x) and R(x) have different distributions because, as we shall see, first and subsequent arrivals have different distributions when the velocity process is Markov. The distribution of arrival velocities, as viewed by our observer at x, clearly depends only on the velocities actually seen, as indicated in eqn. (30). For the velocity on first arrival, the entire distribution plays a role. The distribution of S(x) may be simulated by starting the particles off with arbitrary velocities and considering the velocity distribution as the particles arrive at some distant point. This works because the velocity process is ergodic. However, this is not only expensive, but runs a very real risk of accumulating large numerical errors. We can achieve relatively instantaneous convergence by using R(x) for the initial velocities. Indeed, the distribution of arrival velocities showed remarkable stability for h < 1; there was some initial movement for h = 1, but convergence was instantaneous for h = 0.5 [26]. Essentially, for X < 1, there are relatively few multiple crossings, so R(x) approximates S(x) surprisingly well. In particular, eqn. (32) and VarR(x)

0

&R(x) =

&V

(29)

=7

&V3

= T

-

= P(1 - X2)

(33)

are good approximations for the mean and variance of the velocity at first crossing. The unexpected subtleties involving multiple crossings and arrival velocities were also studied via simulations. Table 1 summarizes the frequency of multiple crossings at steady state (based on 2000 particles per data point) for coefficients of variation between 0.1 and 1.0. Notice that double

235 TABLE 1 Velocities on arrival at milestones Crossing number

Coefficient of variation: h 0.1 0.2 0.3 0.4

’ Mean Standard deviation 2 Mean Frequency 3 Mean Frequency 4 Mean Frequency

1.01 0.10

1.04 0.19

1.09 0.29

* * * * * *

* * * *

* * * * * *

* *

0.5

0.6

0.7

0.8

0.9

1.0

1.15 0.36

1.24 0.43

1.34 0.50

1.45 0.57

1.56 0.63

1.70 0.70

1.81 0.76

0.30 < 0.1% * *

0.47 0.3% * *

0.55 1.0% 0.56 0.1% * *

0.73 2.0% 0.80 0.2% * *

0.87 3.4% 0.75 0.4% 0.71 < 0.1%

1.01 4.3% 0.94 0.6% 0.82 < 0.1%

1.09 5.7% 1.04 0.8% 0.90 < 0.1%

* *

* *

*Insufficient data for entry to be meaningful.

crossings are practically nonexistent for X < 0.4 and remain rare (<2%) until X reaches highest experimental values (h = 0.7). Table 1 also reports the mean velocities on first and subsequent arrivals. For X < 0.3, negative velocities are very rare, so the mean of S is given approximately by eqn. (32). Notice that second, third, etc., arrival velocities are indistinguishable from each other, but are substantially smaller than first crossings. This surprising result is easily understood in terms of the Markov model. If a particle is on its second (or third, etc.) crossing, then at some (probably recent) time in the past, it had a negative velocity. The autocorrelation of velocities then yields a relatively slow crossing. Only the most recently observed negative velocity matters, by the Markov property, so all subsequent crossings follow the same distribution. Calculation of the moments of the distribution of arrival velocities is straightforward, but tedious. Using eqns. (14) and (23) and performing the integrations indicated in eqn. (31), we get &R(X) = 1 + x2

x3 AB+X

(34)

where

0

A=@; and B = flmei/(*h’) Noting that

r3@(r) r tR*(x) = O rw&9 dr

dr (37)

0

we similarly obtain

var R(x)

= A*(1 - A*)

+ h3[(1 + 2h*)AB + X(1 + A*)]

(38)

(AB + ii)*

Both eqn. (34) and eqn. (38) have been algebraically simplified to permit comparison with eqns. (32) and (33) respectively. Owing to the magnitude of B, the additional terms are negligible for small h. Table 2 compares theoretical values of &Z(X) with those from the simulation. The latter are obtained from the means for first and subsequent arrivals, weighting each mean according to the frequency with which the crossing occurs. The agreement is excellent. Owing to the inclusion of the smaller velocities of second and subsequent crossings, the standard deviation R(x), also shown in Table 2, is slightly greater than that for first crossings (Table 1) .

(35) THE ASYMPTOTIC DISTRIBUTION

(36)

OF T(x)

Let T(x, y) denote the time required (in steady state) to travel between milestones x and y; i.e., the elapsed time between first

236 TABLE 2 Theoretical and experimental arrival velocities Arrival velocity

Coefficient of variation X 0.1 0.2 0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

Theoretical mean Experimental mean Theoretical standard

1.01 1.01 0.10

1.16 1.15 0.37

1.24 1.24 0.44

1.34 1.33 0.52

1.44 1.43 0.59

1.55 1.53 0.65

1.66 1.64 0.72

1.78 1.77 0.79

1.04 1.04 0.20

1.09 1.09 0.29

deviation

arrivals at x and Y. Travel from x to x + Y + z must go through 3c+ Y, so

Variance of First Passage Times

T(X, X + Y + 2) = iln(X,X + Y) + T(X + Y, X + Y + 2)

(39)

Taking expected values in eqn. (39) yields [T(Y + 2) = &T(Y) + &T(z)

(40)

since the distribution of T(x, Y) depends only on the journey’s length. Voyages from the origin to RX must pass each intermediate milestone, so T(nx) = ZT[jx,

0’ + 1)x].

(41)

These summands are iden ticdy distributed and dependent. Unless this dependence is very strong, Z’(m) will be approximately normally distributed by the Central Limit Theorem [23]. In fact, over large time-scales, the position process generated by eqn. (10) behaves like a Wiener process with drift /J [25]. This process has independent increments, so the dependence is sufficiently limited. Equation (40) may also be used to deduce exact mean transit times [26]

T(x)

&= [ X 1 =

-&CT(x) &s-l(x)

= p-1

(42)

The dependence in eqn. (41) prevents a similar analysis for Var T(x), but we can expect it to increase linearly in distance and standardized variance; i.e.,

kxa2

VarT(x)=

-

WP2

=

kxh2

(43)

,o 0

4

Dimensionless Distance

0

- Coefficient

I.2

lb

of Variation: xx’

Fig. 5. Approximate linearity of Var T(x).

Our simulations support rough linearity over a broad spectrum of x and h values. These results are summarized in Fig. 5, which is based on two thousand (2000) particles per data point. Actually, Fzincreases slowly (with xX2) from about 1.8 to its asymptotic value of 2.0 (as x + 00, T(x) is given approximately by the first-passage time in the Wiener process with drift p). Figure 6 compares normal and log-normal fits to simulated transit data (again based on 2000 particles per data point) at finite distances x. Under the null hypotheses, the Kolmogorov-Smirnov statistics will rarely (1% and 15%) appear above the horizontal lines. Consequently, the KolmogorovSmirnov statistics [22] reject normality for all distances x < 16 at virtually any level of significance. On the other hand, by x = 6, the log-normal fits are generating KS statistics of the size typically found under the null hypothesis (i.e., KS E [0.6, 0.71).

231 Kolmogorov Normal

Smlmov

and Lognormal

I

4

Distance

Estimated Density of log (l(x)/x)

Statistics Fits forT(x)

I.2 x

‘Fig. 6. Test of goodness-of-fit of the normal and lognormal distributions to the simulated first-passage data. The upper and lower horizontal lines are, respectively, the 1% and 15% upper percentage points of the standardized KS statistic.

We also obtained density estimates [26] (based on 2000 particles at X = 0.5) for log, [T(x)/x] to help assess the log-normal approximation. Initially, log[T(x)/x] is asymmetric, but the convergence in Fig. 7 is stunning. By x = 6, log[T(x)/x] is clearly normal, and at x = 16, the fit is truly remarkable. Table 3 summarizes the sample sizes necessary to detect discrepancies from the normal and log-normal fits. Since Koglin [20] based his analysis on 30 - 40 particles, his experimental results are in perfect harmony with our model.

Fig. 7. Convergence of the distribution of transit times to log-normality.

distances travelled in a fixed time, are normally distributed. The model provides a reasonable fit of our own limited data in I [17], which are of the latter form. Our simulations based on this model yield first-passage times which are asymptotically log-normal. Moreover, the convergence is so rapid that, even for x = 0.5, log-normality is

CONCLUSIONS

Our three-parameter Markov model predicts that instantaneous velocities, and hence TABLE 3 Sample sizes to detect non-normality and non-log-normality Distribution Distance x

Normal 1 4

6

8

12

16

Log-normal 0.5 1

2

3

4

5

Sample size ff = 0.15 cu=o.o5

25 30

150 200

180 240

250 330

415 550

90 120

210 280

415 550

835 1100

1480 2000

105 140

105 140

238

an effective fit for practical sample sizes. Koglin’s data support the three-parameter model and show that is describes the behavior of real slurries. Thus, the model predicts not only the general behavior of sedimenting slurries, but also the nature of the distribution of velocities. This guarantees the central role that the three-parameter model will play in future research on stochastic sedimentation.

A V 8 h /J V

P a 7 @

ACKNOWLEDGEMENTS

This paper is based in part on one presented at the 1st World Congress, Particle Technology, Nuremberg, April 16 - 18,1986. This work was supported by grants from the National Science Foundation (CPEM17673), the Natural Sciences and Engineering Research Council of Canada (A1268) and the Queen’s University Advisory Research Council (380-680).

discriminant time concentration parameter coefficient of variation of velocity mean of steady-state velocity variance of hs autocorrelation of velocities constant in eqn. (10) indicating strength of white noise time-increment standard normal cumulative distribution function

Script expectation operator order of magnitude notation Subscripts, superscripts and marks over symbols indicates first arrival f dummy index i * indicates a dimensionless variable 1 indicates an estimated variable _ indicates a sample mean

LIST OF SYMBOLS

A B i D fe g he i, j

k m

; P Qe Ik S s, t T UO u,

W

u

X

X Y Y, 2 P

function of X defined by eqn. (35) function of h defined by eqn. (36) critical value of n sphere diameter container diameter steady-state crossing velocity pdf distance from milestone velocity transition pdf dummy indices of summation constant mean of he number of summands normal distribution probability measure steady-state velocity pdf velocity random arrival velocity random velocity at first crossing times random transit times Stokes velocity velocities white noise position random particle position log-transformed transit times position autocorrelation parameter

REFERENCES 1 G. J. Kynch, Trans. Faraday Sot., 48 (1952) 166. 599. 2 H. Brenner, Chem. Eng. Sci., 19 (1964) 3 E. M. Tory and D. K. Pickard, Proc. 2nd World Congr. Chem. Eng., Montreal, Oct. 4 - 9, 1981, Vol. 5, 254. M.Sc. Thesis, London Univ. 4 R. P. Boardman, (1961). Proc. Symp. 5 B. H. Kaye and R. P. Boardman, Interaction between Fluids and Particles, Inst. Chem. Engr. (1962) 17. R. Johne, D. Zng. Dissertation, Universitiit Karlsruhe (1965). R. Johne, Chem. Eng. Tech.; 38 (1966) 428. E. M. Tory and D. K. Pickard, Proc. Atlantic Math Days 1983, University of New Brunswick, 1984, 84. 9 D. K. Pickard and E. M. Tory, in G. J. Umphrey and I. MacNeil (eds.), Advances in the Statistical Sciences, Vol. 4, Stochastic Hydrology, Reidel, Dordrecht, in press. at 10 D. K. Pickard and E. M. Tory, Paperpresented 22nd Can. Chem. Eng. Co&, Toronto, Sept. 20, 1972. at 11 D. K. Pickard and E. M. Tory, Paperpresented 24th Can. Chem. Eng. Conf., Ottawa, Oct. 23, 1974. 12 D. K. Pickard and E. M. Tory, J. Math. Anal. A&., 60 (1977) 349. 13 D. K. Pickard and E. M. Tory, J. Math. Anal. A&., 72 (1979) 150. 14 E. M. Tory and D. K. Pickard, J. Math. Anal. Appl., 86 (1982) 442.

239

15 E. M. Tory and D. K. Pickard, in B. M. Moudgil and P. Somasundaran (eds.), Flocculation, Sedimentation and Consolidation, A.I.Ch.E., 1986, p. 297. 16 B. Koglin, Dechema Monogr. 79l3, 235 (1976) 1589. 17 E. M. Tory and D. K. Pickard, Can. J. Chem. Eng., 55 (1977) 655. 18 W. Gosele and R. Wambsganss, Ger. Chem. Eng., 6 (1983) 351. 19 B. Koglin, D. Zng. Dissertation, Universitiit Karlsruhe (1971). 20 B. Koglin, Chem. Zng. Tech., 43 (1971) 761. 21 T. E. Pollock, Ph.D. Thesis, McMaster Univ. Hamilton, Ont. (1981). 22 M. A. Stephens, J. Amer. Stat. Assoc., 69 (1974) 730. 23 S. Karlin and H. M. Taylor, A First Course in Stochastic Processes, Academic Press, New York, 1975. 24 M. V. Mathews, &SC. Thesis, Harvard Univ. (1983). 25 S. Karlin and H. M. Taylor, A Second Course in Stochastic Processes, Academic Press, New York, 1981. 26 D. K. Pickard, E. M. Tory and B. A. Tuckmann, Technical Report SSP7, Harvard University, Cambridge MA (1985).

APPENDIX QUADRATIC TIMES

INTERPOLATION

OF TRANSIT

The problem is to determine, from the positions and velocities at the beginning and end of an increment, whether the particle crosses a given milestone. Consider the nth time increment. Suppose that initially (Le., at time (n - 1)~) the particle is x units from the milestone and has velocity u. Also, suppose that later (Le., at time YU) its velocity is u. Then the distances from the milestone during this increment are given by

(A3)

g(r)) = 0

Thus, interpolating transit times requires finding the roots and optima of g. Obviously, g has real roots if and only if the discriminant, A=U*+2x=

(A4)

7

is non-negative, in which case the roots are -u + @

(A5)

V = (u - u)/r When inequality (A2) holds, the choice -u+

$i

7)f = (u - u)/r 2x

W)

= u+-fl

is always non-negative, and is the smaller root when both are non-negative, and hence is our only candidate for first passage. Since g is quadratic, the minimal value of g for 0 d 77< T occurs either at one of the endpoints, or at the unique critical point c. The specific choice, of course, depends on the initial and final velocities, u and u. Indeed, i

mh

g(c)

if u > 0 > u

g(r)

if 0<-u
g(rl)=

or if u,u > 0

i

OQlj
Ig(O)

otherwise

At the end of the time-increment, tion is g(7) =x -

u+u 7

2

(A7)

the posi-

(A@

Obviously, g(r) < 0 implies the existence of positive roots, and hence A > 0. The position corresponding to the unique critical point is

g(q,=x-b/i[.+(y)f] dt

.

7.4‘7

(AlI The milestone is crossed during this increment if and only if

mb g(q)<0

(A2)

g(c) = x -

2(u - V)

W)

Then, g(c) < 0 is equivalent to A > 0. It follows that first passage occurs during our increment at time nf if and only if either

u+v

o
(i) 2

in which case, first passage occurs for the least non-negative r) satisfying

or

7> x

240

(ii) A > 0 with u > 0 > u. Implicitly, these conditions entail 0 < rlf G r. Notice that, in case (ii), the particle may cross the milestone at qf, but later retreat back across it before time T. Indeed, this will

occur frequently for slow-moving particles near the milestone. Consequently, quadratic interpolation is necessary, in simulations having negative velocities, to avoid overestimating first-passage times.