Do We Understand Hydrodynamic Dispersion in Sedimenting Suspensions?

Do We Understand Hydrodynamic Dispersion in Sedimenting Suspensions?

Francisco Ricardo da Cunha

Otávio Silva Rosa

Aldo João de Sousa

Department of Mechanical Engineering, University of Brasília

Campus Universitário

70910-900 Brasília, DF. Brazil

[e-mail]

Microstructural changes, that is an important feature for the understanding of the velocity variance in sedimentation is investigated with numerical simulations. The simulations are used to describe velocity fluctuations and hydrodynamic dispersion in a suspension of interacting point-particles sedimenting in a rectangular box with periodic sides and impenetrable bottom and top. It is observed how the positions of the particles evolve in a finite container. The suspension that was initially random in the gravity direction only, tends to be fully randomized as a result of the relative arrangements of the particles and the hydrodynamic interactions between them. The computer simulations, based on statistics over a significant number of particle configurations, suggest velocity variances and diffusivities dependent on the size of the simulated system but with anisotropy in velocity fluctuations and diffusion coefficients nearly independent of the box size.

Keywords: Velocity fluctuations, sedimentation, microstructure, dispersion

Introduction

Suspended particles subject to sedimentation do not generally move relative to the fluid with a constant velocity, but instead experience diffusion-like fluctuations in velocity due to interactions with neighboring and the resulting variations in the configurations of the suspended particles. These velocity fluctuations are observed in non-Brownian suspension flows with very small particle Reynolds numbers. Such fluctuations have a long-time behavior characteristic of diffusion processes and their effect is now called hydrodynamic self-dispersion (Cunha 1995). A small particle (with a diameter less than about 1mm) usually diffuses due to random impacts of the fluid molecules in which it is suspended. A larger particle lacks this Brownian motion, but often exhibits a diffusive motion originated from its long-range hydrodynamic interactions with other large particles in the suspension. This dispersion phenomenon is important for understanding of mixing process which inhibit separation (Davis 1996). The related phenomenon of shear-induced hydrodynamic diffusion in sheared monodisperse suspensions of spheres have been investigated experimentally (e.g. Leighton and Acrivos 1987) and theoretically (e.g. Cunha & Hinch 1996a; Loewenberg & Hinch 1996, Almeida 1998).

As a suspension contains many particles dispersed in a fluid, the motion of one particle creates velocity and pressure _field that exert forces on the neighboring particles and affect their motion. At low Reynolds numbers such a disturbance flow is propagating via fluid by the mechanism of vorticity diffusion, and because of this we say that the particles interact hydrodynamically. Long-range multibody hydrodynamic interactions play a key role in the motion of an individual sphere settling in the midst of a suspension of non-Brownian spheres. Indeed, each sphere undergoes a random-walk motion due to the fluid velocity disturbances caused by the surrounding ones (Ham & Homsy 1988). The dimensional analysis leads to a number of simple and powerful results that explain this greatly enhanced mixing in applications such as sediment transport (Hinch 1988, Cunha 1997). After Batchelor (1972) the average velocity in sedimentation can be successfully predicted theoretically. In the present time, however, there is a discrepancy between experiments (Nicolai & Guazzelli 1995) and numerical simulations (Ladd 1993, Koch 1994, Cunha (1995) and Ladd (1997)) regarding the dependence of the sedimentation velocity variance on the container size.

Experiments performed by Nicolai & Guazzelli (1995) using well-stirred suspensions showed that the values of the velocity fluctuations did not vary significantly when the vessel width was increased by a factor of four. Cunha (1995) reported that computer simulations of velocity variances in sedimentation using random and independently initial distributions of particles in all directions of space yield results which did increase with increasing the container size. The numerical simulation results follow the theory and scaling O(/a) predicted by Caflish & Luke (1985), Hinch (1988). Here, U_s= 2Dr a² g/9m is the fall speed of an isolated sphere of radius a, f is the volume fraction, and is the size of the container, m is the fluid viscosity, Dr denotes the difference between the density of the solid particles and fluid, g is the acceleration due to gravity. The screening mechanism of Koch and Shaqfeh (1991) is the only presently available theory that could lead to velocity fluctuations independent of the vessel size. However, the experiments (Nicolai et. al. 1995) and the lattice-Boltzmann numerical simulations of large-scale systems by Ladd (1997) were unable to verify the hypothesis that there is a deficit of one particle surrounding any given particle of the suspension assumed by this screening theory. At very low concentrations in a thin box, Segre, Herbolzheimer and Chaikin (1997) found a f^1/3 dependence, and an independence of the wider of the horizontal dimensions if it exceeded a certain correlation length _c = 10af ^--1/3 . Curiously this observed correlation length in the velocity fluctuations is somewhat greater than the narrower of the horizontal dimensions.

Fluctuations are often sensitive to subtle changes in the underlying suspension microstructure, which are difficult to observe experimentally. There should be some suspension microstructure that leads to velocity fluctuations depending on the particle volume concentration only. Searching for other microstructures than the ones already simulated by Cunha & Hinch (1996b), a new kind of initial distribution of particles is proposed and tested in this work. Numerical simulation is then used to examine the microstructural changes directly, as well to predict their effect on fluctuations in particle velocity and particle hydrodynamic dispersion. The long-range nature of the hydrodynamic interactions in sedimentation requires care for the simulation of an infinite suspension. At the present, the most appropriate way to mathematically treat such kind of problem is by the use of periodic boundary conditions, which represents the sedimenting box as a spatially periodic array of identical cells. In this article, efficient simulations of sedimentation that incorporate a strong anisotropic suspension microstructure are described and the particle velocity fluctuations and dispersion are properly quantified, presented and analyzed. Part of the results presented in this article has been summarized in a conference paper (Silva Rosa & Cunha 1998).

Description and Formulation of the Problem

Scaling Argument and Particle Distribution

Consider a random monodisperse dilute suspension, where is natural to suppose that for the purpose of considering the influence of neighboring particles on a test particle, it suffices to replace all the other particles in the suspension by point singularities such as point forces. Assume that m = Dr a³ is the net particle mass, and a typical statistical fluctuation in density , of a region of size >>a of such suspension scales as (Hinch 1988)

where N is a typical fluctuation in the number of the particles. Hence when a box of volume ³ containing N particles is divided into two equal parts by a vertical plane, one half of the box will contain (N/2 – ) particles, whereas the other half will contain (N/2 + ) This unbalance drives convection currents during the sedimentation process.

Now, the density number fluctuation leads to a fluctuation in the weight of mg

. At the time scale it takes vorticity to diffuse over the length , t ~ ² /m, this buoyancy force is balanced by the viscous drag associated with the driven flow . Hence

Then, by substituting g ~ a^-2mU_s/Dr and N ~ f( /a)³ into the above equation, one obtains

If the particles velocity remain correlated by a time t_c ~ (/), the hydrodynamic self-diffusivity, , scales as

The above scaling arguments help to explain how velocity fluctuations and hydrodynamic self-diffusivity in a random dilute monodisperse sedimenting suspension may be dependent on the system size (for more details see Cunha 1997). The scalings identifies a fundamental mechanism of dispersion and mixing operating in sedimentation and it provides a first insight on the mechanism by which horizontal density fluctuations lead to large velocity fluctuations in sedimentation.

In this article we propose an initial non-homogeneous distribution of particles that is random in the gravity direction only. The main purpose of this rather unnatural kind of distribution is to eliminate the presence of convection-driven secondary flows in the horizontal directions, what in turns would influence the velocity fluctuations of the particles. For this end, the initial configurations used in all simulations here are generated placing the particles "regularly" in the horizontal directions, and random and independently in the vertical (parallel to the gravity) direction inside a rectangular box of dimensions x x . The distribution is made so that each particle center is separated from the other ones by a distance greater than the dimensionless particle diameter in the same configuration, avoiding overlaps. A typical initial configuration for the concentration f = 3% and the dimensionless radius a/ =0.05 is shown in Figure (1).

Governing Equations

The long-range nature of the hydrodynamic interactions in sedimentation requires care for the simulation of an infinite suspension. At the present, the most appropriate way to mathematically treat such kind of problem is by the use of periodic boundary conditions, which represents the sedimenting suspension as a spatially periodic array of identical cells.

Consider rigid spherical particles sedimenting in a incompressible Newtonian fluid of dynamic shear viscosity m and density r, since the inertial effects of the fluid are negligible and the time scale is large compared to the viscous relaxation time (a²/n), the appropriate equations of the fluid motion in the usual Eulerian description of an inertial frame are the pseudo-steady Stokes equations

and

where u and P = p-r g ^. x are the velocity and the modified pressure fields. Equations (3) and (4), valid within the fluid, are supplemented by a non-slip boundary condition on the surface of each particle a =1, 2, , N, u = U_a +W_a x (x-x_a) on | x-x_a |= a. U_a and W_a are respectively the velocity and the angular velocity of the sphere a with radius a and at center x_a (t). It should be noted that in a dilute suspension of sedimenting particles the condition for the Stokes equation to be valid is Re << a/ where is the characteristic length scale of the interaction range (i.e. the mean interparticle distance). With ~ n^-1/3 ~ af^-1/3the particle Reynolds number must satisfies

In addition, we consider the limit case where Péclet number is high which corresponds to a suspension of non-Brownian particles,

where D is the ordinary diffusivity of a very dilute dispersion of independent spheres of radius a as first derived by Einstein (1956), k is Boltzmann's constant and T is the absolute temperature.

Owing to the long range of the velocity disturbance in a dilute sedimenting suspension, we will be interested in particle interactions at separation that are large compared with particle radius. Therefore, with an error O(f) in the sedimentation velocity of the particles, the simplest level of point-force approximation may be used for describing the hydrodynamic with fluid velocity governed by (Saffman 1973)

Here d(x-x_a) denotes Dirac's delta distribution, and f is the hydrodynamic force exerted on the fluid by a particle. In the absence of particle inertia, f^a = (4/3) pa³ Dr g; which is just the net force of gravity with the buoyancy removed.

We consider the condition with no mean flow, i.e.áuñ = 0. Taking an average of equation (5) we find an expression for the mean pressure gradient ,

so

which is the global balance between the average pressure gradient and the average force, áfñ =(1/N) ñS_af^a the particle exerted on the fluid. In particular (7) suggests that in point-particle approximation the pressure P in (5) may be adjusted to reflect the macroscopic increase in the density of the fluid fDr = n êf ê/g. The pressure is then decomposed into a periodic part p'(x) and a linear part áPñ,with gradient given by (7), P(x) = áPñ + p'(x).

Substituting the value of the reduced pressure P(x) into (5) one obtains

Periodic Stokes Flow

We shall show first the solution of the equations (6) and (8) within a rectangular box of dimensions x x 2h, satisfying the following periodic boundary conditions:

where u; v; w; are the velocity components of the uid motion in directions x; y; z respectively. For this purpose, it is convenient to introduce three-dimensional finite Fourier transforms of the velocity u(x) and the pressure field , according to equations

Now if we take the Fourier transform of (6) and (8), ones obtain respectively that

and

Considering the suspension monodisperse, taking the dot product with k in the transformed equation of motion (13), and using the continuity equation (12) this gives

and substituting the above expression back into (13) gives an expression for , namely

Providing that the solution must satisfy the imposed periodic boundary conditions (9), the pressure and velocity fields are expanded as three dimensional Fourier series as follows,

Hence the solution of the equations (6) and (8) are respectively:

and

where the sums S(x ê k) are given by

I = d_ije_ie_j denotes the unit second-rank tensor and V is the volume of the periodic box V=²x2h. The notation S(x ê k; x_o) represents the summation in which the vector x is switched by (x-x_a). The wave number vector k is defined as k =(b₁/, b₂/, b₃/2h) where b₁; b₂; b₃ are positive or negative integers (0; ± 1, ± 2;...). The exclusion of the term k = 0 is a direct result of the average pressure balancing the average force the particles exert on the fluid such as indicated in equation (7). The sums (19) are performed over a three-dimensional periodic lattice in which each (cell numbered by the index g) may contain N particles (numbered by the index a). The lattice points are given by the vector x_g = (g₁; g₂; g₃2h) with g₁; g₂; g₃ = 0; ± 1, ± 2..., and a particle has position vector r_a =x_a +x_g. Here x_a is the position vector of the particle with respect to the origin of its cell.

The difficulty in evaluating the series into (17) and (18) is always the question of their convergence. A great deal of attention to this particular problem was given by Ewald (1921) when examining the slow convergence of lattice sums arising from Coulomb interactions (which also decrease slowly with distance). He proposed a powerful technique which splits up the original sums (only conditionally convergent) into two more rapidly converging forms, in the physical space and Fourier spaces. Ewald's method not only improve the accuracy as well it is essential for simulating systems with periodic boundary conditions. In the existing literature on calculations of transport properties in suspension several published papers have already applied Ewald's technique successfully (e.g. Hasimoto 1959, Sangani and Acrivos 1982 Brady et al. 1988 and Beenaker 1986).

The next step is then an improvement of the convergence of the series (19) by applying the so-called Ewald summation technique and the Poisson summation formula (see Hasimoto 1959, Cunha 1995, Van de Vorst 1996 for details) to rewrite the sums in (19) as follows

where G_n(y) is the incomplete gamma function which satisfies the recurrence relations

After substitution of equations (20) and (21) in the fundamental solution (17) and (18), and working out all the derivatives, we arrive at the fundamental periodic solution for Stokes flows due to a lattice of f point forces.

which has been made dimensionless using a, U_s and 6pm U_s for the reference scales of length, velocity and force, respectively. P^ps,P^rs , J^ps and J^rs are defined as being dimensionless Green's functions in the physical (ps) and reciprocal (rs) spaces which conform with the periodicity of the flow (Cunha 1995):

The solution for p(x) and u(x) given in terms of the series (23) and (24) converge exponentially fast, with the convergence rate controlled by the arbitrary geometrical parameter l > 0. For optimal convergence, l should be chosen neither too small nor too large. Beenakker (1986) recommends l =p^1/2V^-1/3 for simple cubic lattice, giving the equal rates of convergence for the two sums.

Impenetrable Container

The method of images which plays so important a part in the mathematical theory of elasticity, fluid mechanics and heat transfer, is peculiarly adapted to the solution of the problem of the sedimentation of many interacting particles. The image system used in this work consists of stokeslets equal in magnitude but opposite in sign to the initial stokeslets. We then construct the appropriate Green's functions for the velocity and pressure field associated with a stokeslet actuating in the gravity direction. Now, all components of the velocity field (u; v; w) are periodic in x and y with period , the horizontal components u; v periodic in z with period h, but the vertical component w satisfying an impenetrable bottom and top condition of vanishing vertical velocity. Hence the following boundary conditions are imposed

Clearly the solution (24) satisfies the governing equations, but not the boundary conditions specified in (28). To overcome this, the procedure is essentially to consider a linear combination u(x)=u(x;x^s_a)-u(x;xⁱ_a) which satisfies the governing equations and the new boundary conditions (28). These two terms can be thought as the solution resulting from an initial stokeslet located at x^s_a=(x_o, y_o, z_o) and a fictitious stokeslet of equal magnitude but opposite sign at the image point xⁱ_a= (x_o, y_o, -z_o) (see Figure 2). Using such an image system we are able to obtain the solution for the Stokes flow induced by a lattice of stokeslets with side periodicity and impenetrable top and bottom.

where the tensor , and the function Q are defined as follows

Mobility Problem

Now the problem of N spherical particles free of inertia settling within an impenetrable box with periodic sides of dimensions xxh is considered. Let x_a denote the position of the particle a. Suppose an external force f is exerted on particle a and let U^a be its translation velocity. Then the equation (29) can be rewritten in an appropriate formulation of hydrodynamic interaction which relates the velocity U^a and the forces f as follows

with particle trajectories being obtained by integration of the kinematics equation

Here f^m =(0,0,1). The prime on the first sum means that the term a =m in cell g =1 has been excluded.

Equations (32) and (33) will be applied to examine the dynamics of N point-particles sedimenting and interacting hydrodynamically within a container of impenetrable boundaries. This type of formulation represents a mobility problem with hydrodynamic interactions O(N² ), calculated by using pairwise additivity (i.e. superposition of velocity) in the mobility matrix. We emphasize that our purpose here is not to perform detailed calculations of particle interactions. Rather, we aim to explore the physical processes giving rise to velocity fluctuations and its consequence when particles sediment. It will be shown next that point particle interactions suffice for the production of random particle migration needed to reproduce the phenomenon of hydrodynamic dispersion.

Results and Discussion

In this section, a brief description of some simulation features and the most relevant results of this work are presented.

Simulations Conditions and Preliminary Tests

In order to verify the global accuracy and to optimize the efficiency of the computer simulations several tests were performed before starting with the main computations. The optimum value for the l parameter had to be investigated given the particular features of the suspension initial microstructure. The l_optimum was examined for the concentrations of 1%, 3% and 4%, and for 27, 125 and 343 periodic boxes. The container aspect ratio used in all simulations was h/=3. The best value found under the conditions set was l=p (Silva Rosa, 1998). A few simulations were made with the sedimentation of one single particle in the box in order to verify some characteristics of the numerical routine and obtain preliminary results necessary for performing the main simulations. At first, the tests consisted of one particle of dimensionless radius a/ = 0.05 being released in a certain initial position close to the top and its motion followed as far as it reached the impenetrable boundary at the bottom. The integration of dX/dt = U was solved via a fourth order Runge-Kutta scheme. The magnitude of the non-dimensional time step in the Runge Kutta procedure was set up as Dt £ Ka/ - which is about the dimensionless time spent by a single particle to fall across its own radius. K = 0.5 gave enough precision for the present simulations. The value chosen for the one-particle simulations was Dt = 0.025.

The three-dimensional fluid motion was described by the dimensionless governing equation (29). The incomplete Gamma functions G_n(x_arg ), associated with the mobilities in the series, were numerically estimated by using asymptotic expansions for the limit conditions 0< x_arg <<n + 2(x_arg®0) and x_arg >>n + 2 (x_arg® ¥), and a continued fraction for moderate value of the argument x_arg . A cut-off parameter r_min was defined so that interactions between particles separated by a distance smaller than 2r_min are excluded so as to avoid the singularity in calculating the incomplete Gamma functions. To accelerate the calculations of Ewald summations and thereby obtain a considerable savings in computer time, the incomplete gamma functions were then tabulated for the values of n which appear in the solution series (29), n = (-1/2; 0; 1/2; 1), and their calculation required for each particle pair interaction were determined by simple liner interpolation away these tables.

Three simulations were made with 27, 125, and 343 periodic cells, and in all the cases, the perpendicular components of the particle velocity were found to be null during the whole sedimentation evolution. This happened because the test particle suffers little influence from the other particle in the lattice. The parallel component (i.e. gravity direction) remained constant and greater than zero when the particle was not near the top or the bottom. The information about the zone of influence of the top and bottom is essential to define a bulk region of the suspension in which the statistical data analysis should produce meaningful results. It was assumed that, in the following simulations only the particles comprised in the region 0.30 £ h/ £ 2.70 would be considered during the calculations of the suspension transport properties. By doing this, the statistical data analysis is performed in a region of the box where the variations in the individual velocity of a reference particle are only due to the pairwise interactions with the other particles. The differences of accuracy verified in the tests with 27, 125 and 343 periodic boxes were of the order of 0.1%. Thus, the arrangement of 27 periodic cells (3x3x3 boxes) proved to be the best choice.

Statistical Analysis of the Suspension

The evolution equation (32) and (33) are the heart of the dynamic simulation here. Given an initial configuration of particles, the equation (33) is integrated in time to follow the evolution of the suspension microstructure. In order to use the simulations to determine macroscopic properties of the sedimenting suspension we derive below corresponding average expressions.

The determination of the structure factor of the suspension is the simpler way of obtaining the density number fluctuations. It is defined as being

We were principally interested in looking at Fourier component corresponding to the box wave number k = (1/, 1/; 1/h) rather than particle-particle separation. The sedimentation velocity is given by the average velocity of the particles in the suspension defined as follows

The fluctuations in velocity of an individual particle in the suspension

is a measure of the deviation of the particle velocity from the mean velocity calculated over all particles in the configuration 1/N(t) or at each time step in dynamical simulation. Hence, we define the variances of the sedimentation velocity as

The normalized velocity fluctuations auto-correlation functions C_j were calculated in our simulations through the following expression

and self-dispersion coefficients were calculated by integrating the velocity auto-correlation functions

with the corresponding correlation t^c times being estimated by the relation

Here the angle brackets denote a sum over all particles, and an average over all configurations or realizations (i.e. an average over time in dynamic simulation).

Static Structure Factor and Velocity Fluctuations at Time t=0

The simulations with many particles started with measurements of the density number fluctuation, called the structure factor of the suspension. The hydrodynamic interactions between the particles were considered pairwise additivity interactions. This procedure results in computation times that scale as O(N²) per configuration or time step. Figure (3) shows the results of twenty simulations of the vertical structure factor with volume concentrations 1% £ f £ 5% and number of particles 54 £ N £ 1024. The averages were obtained over 1000 initial configurations in each case. The fact that all the mean values are pretty close to the unit indicates that the numerical process of locating particles random and independently in the gravity direction produced the expected stratification. The small deviations appear due to the finite size of the system. It is seen in Figure (3) that the vertical structure factor is essentially independent of the number of particles. In particular, a typical fluctuation in the particle density number is found to be /³ in agreement with the scale assumed by the scaling argument described in § 2. The perpendicular structure factor was exactly zero in every configuration of all simulations. This means that the 20000 initial distributions were regular in the horizontal directions, being the particles equally distant from each other in all configurations.

After having verified that the numerical routine did make the expected initial distributions, the investigation of the velocity fluctuation began. Thirty different simulations were made to compute the statistics over 150 independent initial configurations, an O(N²) operations per configuration or time step. This took about 7 hours of computation time in an Alpha Digital work-station with 333 MHz and 200 Mbytes of RAM. In these tests, called "Zero-Time Simulations", no evolution of the sedimentation in the time (and space) was produced. Because it is a mobility problem, the velocities and fluctuations could be calculated based only in the relative positions of the particles in each configuration. The box aspect ratio h/ was kept constant and equal to 3. The volume particle concentrations were varied throughout the range 0.5% £ f £ 5%, and the values of particle radius adopted, 0.02 £ a/ 0.07. In each configuration, the mean components of the velocity (of all particles comprised in the bulk region 0.3 £ h/ £ 2.7 defined above) and their variances were calculated. At the end of every simulation, the averages over the respective 150 initial distributions were then obtained. The parallel and the perpendicular to the gravity direction components of the velocity variance were made dimensionless dividing them by the squared Stokes velocity, U²_s.

In Figure (4), velocity variances in the vertical direction of one typical simulation are presented versus the number of the configuration. It was used, for each independent initial configuration, 176 spherical particles with dimensionless radius equal to a/ =0.05 - what gave a volume concentration of f = 3%. The ultimate goal of these simulations is to verify whether the velocity variances of the particular distributions considered in this work are depended or not on the size of the simulated system. For this purpose, the ensemble averages áU^'’2_^ñ/U²_s and áU^{’ 2}_úúñ/U²_s were evaluated and plotted versus the box parameter f /a predicted by the scaling argument in §2. Figure (5) shows the variances of velocity in a perpendicular direction, and, Figure (6), the variances of velocity in both perpendicular and parallel directions as a function of the scaling parameter f /a.

The results suggest that the velocity variances of particle distributions with stratification only in the vertical direction are best described by linear functions starting at the origin. The following equations were then obtained by linear regressions: áU^'’2_^ñ/U²_s 0.014f /a and áU^{’ 2}_úúñ/U²_s 0.650f /ª. The correlation coefficients were 71.5% and 92.0%, respectively. The values found for the variances of velocity in the horizontal directions stay between one third and one half of the ones achieved in totally random suspensions as simulated by Cunha, Sousa & Hinch (2000). In the vertical direction, they are between three and four times the related values in the referenced publication. Different initial distributions result in diverse mean statistical properties of the suspension. Hinch (1988) expected that if the horizontal fluctuations in the number of density were eliminated, as it was done in this work, both the horizontal and the vertical variances of velocity would be smaller than the values observed in the simulations of Cunha (1995). This prediction, however, was verified only in the direction perpendicular to the gravity. The velocity fluctuation anisotropies, (áU’ ²_úúñ/áU’ ²_^ñ)^1/2, are seen in Figure (7) versus f /a The mean value (»7), is greater than the experimental one reported by Nicolai et al. (1995), approximately 1.9, and than the result of the simulations run by Cunha & Hinch (1996b), Cunha, Sousa & Hinch (2000) which was about 2.5.

Time Developing of the Suspension

First, it was investigated the effect of lubrication forces on the suspension macroscopic properties. For this end, we added into the hydrodynamic force on a particle a an artificial short-range lubrication force acting between all pairs of particles, namely

where çx_i - x_jçis the interparticle space, and C₁ and C₂ determines the strength and range of the repulsive force to prevent particles from overlapping when the suspension evolves in time. Introduction of this extra repulsive force to prevent overlaps is not unrealistic because forces acting between particles in nature and in laboratory practice are often repulsive. Furthermore the pairwise addition of near-field lubrication forces in Stokesian dynamics simulations (Brady and Bossis, 1988) requires times steps prohibitively small to prevent overlaps. After several tests of this model (Silva Rosa, 1998) we conclude that such short range forces only produce a little effect on the suspension average properties. Moreover, computing Ewald's summation for many body interactions with the inclusion of lubrication force would require O(N³) operations (instead of N²) and smaller time steps which limits the size of the system that can be simulated. Thus, we argue that there is no reason to model lubrication force in dilute sedimenting suspension. Readers more comfortable with a continuum description of the suspension physics should know that in sedimentation rather than shear flow the settling of the particles induces a compensating backow of fluid that flows up through narrower channels of the lattices and not through the small gap between particles.

Typical pictures of the sedimentation process obtained with the dynamical simulation are displayed in Figures (8). The picture was collected at different times throughout one arbitrary realization between those ones simulated. The initial configuration used is exactly the same as the one in Figure (1), for f = 3%, a/ = 0.05 and 176 particles. The times presented in Figures (8) are made non-dimensional by the reference scale a/U_s. Before the sedimentation began, the suspension microstructure was formed by 16 vertical parallel lines of 11 particles each, with stratification within the lines. When the sedimentation started, a mixing began to happen destroying the initial regular distribution in the horizontal directions. This mixing is due to the particle velocity fluctuations, that are caused here only by the hydrodynamic interactions. If the volume concentration was f = 0%, there would not be any hydrodynamic fluctuation - what does not happen to thermal fluctuations (Brownian motion), which would exist even in zero particle concentration. Since the Reynolds and Stokes numbers are very low, there were not turbulent structures in the flow, nor particle inertia effects. The dimensionless time step used in the evolution was chosen to be K = 1/5, so that the effects of the velocity fluctuations could be analyzed in detail. Most of the simulations showed that a particle moves mainly downward in the gravity direction, while simultaneously responding to the vertical and lateral fluctuations caused by the continual change in configurations of its neighboring and resulting hydrodynamic interactions. These results provide some evidence that fluctuating positions behaves as it should be for a random walk process, justifying the description of particles migration in sedimentation as a hydrodynamic self-dispersive process.

Microstructural changes, that is the variations in the relative arrangements of the particles, are among the most important and interesting features of the sedimentation process. The distribution of neighboring particles around a reference particle can be characterized by the structure factor of the suspension. Figure (9) shows results of the statistics calculations of F_^(k; 0) for f = 3% and a/ = 0.05.

It is seen that the initial horizontal density fluctuations go way in time, indicating that the suspension tends to be randomized by the effects of the hydrodynamic interactions between the particles.

The long-time behavior of the fluctuations are described by the velocity fluctuation autocorrelation functions and the hydrodynamic self-diffusivities both parallel and perpendicular to the gravity direction. Figures (10) and ⁽¹¹⁾ shows the time development of the normalized velocity fluctuation autocorrelation functions for f = 3% and a/ = 0.05. The simulations reveal that these functions decay in time as a single exponential toward zero, indicating that particle velocity becomes totally uncorrelated after long time. The time to fall through » 60 a. The figures show also the error bars of these results. The results shows that the aspect ratio h/ = 3 is sufficient to reach asymptotic and provide adequate data to determine the hydrodynamic self-diffusivities. This indicates that the particles has enough time to sample the horizontal cross section significantly before settling out so that the diffusivities are long time behavior ones. Note that this occurs in a time scale smaller than the one that many particles have already reach on the impenetrable bottom. The results also suggested a correlation time of the vertical fluctuations as being approximately the time to fall through 20 particles radius which is approximately twice bigger than the correlation of the horizontal fluctuations.

The result corresponding to the hydrodynamic diffusivities both parallel and perpendicular to the gravity direction as function of time is displayed in Figure (12). The numerical simulation here determined values of the D_çç » 3aU_s which are slightly smaller than those reported by experiments and closely to the numerical result, D_çç » 2aU_s, of our previous simulations with fully random suspension considering the same parameter: h/ = 3, a/ = 0.05 and f = 3% (Cunha 1995). The experiments by Ham & Homsy (1988) found this coefficient increasing from about 2aU_s at f = 2.5% to 6aU_s at 5%, and the more recent experiments by Nicolai et al. (1995) reported such self-diffusivity as being approximately 5aU_s at 5%. One possible explanation for this difference is that the comparison with experiments should be made with the understanding that the values of our diffusivities depend on system size (see Figure 13) since they were calculated as the integral of the velocity autocorrelation functions. The laboratory experiments are also at dilute limit (f = 5%), but they have a box about 100 times the particle diameter, which is obviously impossible to copy here due to the limited size of our numerical system for larger number of particles. It should be mentioned that the vertical diffusivity obtained in the present simulation is much smaller than that one predicted by hydrodynamic screening theory of Koch and Shaqfeh (1991) who found D_çç = 0.52f ^-11 (» 17 at f= 3%).

Figure (13)

Since completing the work reported here (Silva Rosa 1998), the lattice-Boltzmann simulations of Ladd (1997) has come to our attention. Ladd reported numerical results of fluctuations and hydrodynamic dispersion in sedimentation for a large homogeneous suspension using 32768 particles (f= 10%) at finite Reynolds number (R_eW = 0.45), based on width of the periodic cell. His results show an anisotropy in correlation time equal 2.5 and the ratio of diffusivities equal to 23 for h/ = 4 that agree well with our numerical results. Surprisingly, this large system simulations closely follow the present results, showing the same system size dependence in velocity fluctuations, O((/a)^1/2), and hydrodynamic diffusivities, O((/a)^3/2), that do not agree wit th experimental results of Nicolai et al. (1995). Ladd also leaves open the question about the existence of a hydrodynamic screening for a random dilute monodisperse suspension in the way of the mass deficit predicted by Koch & Shaqfeh's theory (1991).

Concluding Remarks

In this article direct numerical simulations of a monodisperse dilute suspension of point particles with excluded volume sedimenting at low Reynolds number in a rectangular box with periodic sides and impenetrable bottom and top have been used to describe the microstructure, velocity fluctuations and dispersion of such suspension.

The numerical results reveal that putting the particles on a horizontal lattice and random vertically leads to the (f/a)^1/2 scale for the fluctuating velocity predicted by our scaling argument and by the numerical results of our previous simulations with fully random suspension. After times closely to velocity fluctuation correlation time the more regular initial distribution considered here was strongly randomized as result of the hydrodynamic interactions between the particles. We think that a next step would be to study fuzzy lattices. If the particles are exactly on lattice sites, then the fluctuations are zero. If they are randomly displaced from the lattice sites, then there will be fluctuations. If the randomness is large, then the system would be fully random and the fluctuations would be (f/a)^1/2.

The order of magnitude of the diffusivity parallel to gravity agreed well with those predicted experimentally for dilute suspension, but it depends on the size of the numerical system, increasing like f ^1/2 (/a)^3/2such as found by scaling arguments, theory and large-scale lattice-Boltzmann simulations. The simulations showed degree of anisotropy in velocity fluctuations independent of the system size and in close agreement with experimental measurements. Hydrodynamic self-diffusivity anisotropy was also independent of the system size, but five times larger than the anisotropy found by experiments.

The unsolved problems in suspension are now rather technical. One however is simple to pose and reveals a serious theoretical deficiency and this concerns fluctuations in sedimentation. While the average velocity in sedimentation can be successfully predicted theoretically, we are still unable to renormalize the rms fluctuations. Certainly the problem of determining the full multiparticle hydrodynamic interactions effects in larger numerical systems representing a sedimenting suspension is important and challenging. It is now clear that in a monodisperse sedimenting suspension free of inertia and colloidal forces the fine details of hydrodynamic interactions needs to be considered for changing suspension structure and to provide a hydrodynamic screening for the velocity fluctuations in sedimenting monodisperse dilute suspensions.

Acknowledgments

This work was supported in part by CAPES-Brasília/Brazil and CNPq-Brasília/Brazil. We would like to acknowledge the use of computational facilities at the Department of Mathematics and at Faculty of Technology in University of Brasília. We are grateful to Prof. Homero Luiz Picollo (CIC/UnB) for his valuable help providing the image analysis techniques. We express our thanks to Prof. John Hinch (DAMTP-Cambridge) for his helpful comments and suggestions during this work.

Cunha, F.R., Sousa, A.J. & Hinch, E.J., 2000 Numerical simulation of velocity fluctuations and dispersion of sedimenting particles. Chemical Engineering Communications (accepted).

Article received: February, 2001. Technical Editor: Angela Ourívio Nieckele.

• Almeida, M.H.P. 1998. Numerical simulation of particle deformation in dilute emulsions, (in portuguese), M.Sc. dissertation, University of Brasília, Brasília, DF, Brazil.
• Batchelor, G.K. 1972 Sedimentation in a dilute suspension of spheres. J. Fluid Mech. 240, 651.
• Brady, J.F. & Bossis, G., 1988. Stokesian dynamics. Ann. Rev. Fluid Mech. 20, 111.
• Brady, J.F.,Phillips, R.J., Lester, J.C., & Bossis, G. 1988 Dynamic simulation of hydrodynamically interacting suspensions. J. Fluid Mech. 195, 257.
• Beenakker, C.W.J. 1986 Ewald sum of Rotne-Prager tensor. J. Chem. Phys. 85, 1581.
• Caisch, R.E. & Luke, H.C. 1985 Variance in the sedimentation speed of suspension. Phys. Fluids 28, 759.
• Cunha, F.R. 1995 Hydrodynamic dispersion in suspensions. PhD Thesis, Department of Applied Mathematics and Theoretical Physics, Cambridge University.
• Cunha, F.R., & Hinch, E.J., 1996a Shear-induced dispersion in a dilute suspension of rough spheres. J. Fluid Mech., 309, 211.
• Cunha, F.R., & Hinch, E.J., 1996b Hydrodynamic self-dispersion of sedimentating non-Brownian spheres. VI Brazilian Congress of Engineering and Thermal Sciences 3, 1375.
• Cunha, F.R. 1997 On the fluctuations in a random suspension of sedimenting particles. J. of the Braz. Soc. Mechanical Sciences, 19 4, 474.
• Davis, R.H., 1996 Hydrodynamic diffusion of suspended particles: a symposium. J. Fluid Mech, 310, 325.
• Davis, R.H. and M.A. Hassen, M.A. 1988 Spreading of the interface at the top of a slightly polydisperse sedimenting suspension. J. Fluid Mech. 196, 107, Corrigendum," 1989 J. Fluid Mech. 202, 598.
• Einstein, A. 1956 Investigations on the theory of the Brownian Movement. Dover Publications, INC..
• Ewald, P.P. 1921 Die berechnung optischer und elektrostatischer gitterpotentiale, Ann. Phys. Lpz. 64, 253.
• Ham, J.M. & Homsy, G.M. 1988 Hindered settling and hydrodynamic dispersion in quiescent sedimenting suspensions, Int. J. Multiphase Flow 14, 533.
• Hasimoto, H. 1959 On the periodic fundamental solutions of the Stokes equations and their application to viscous flow past a cubic array of spheres. J. Fluid Mech. 5, 317.
• Hinch, E.J. 1988 Sedimentation of small particles. In Disorder and Mixing edited by E. Guyon, J.P. Nadal, and Y. Pomeau, Kluwer Academic, Dordrecht, p.153.
• Koch, D.L., & Shaqfeh, E.G. 1991 Screening in sedimenting suspensions. J. Fluid Mech 224, 275.
• Koch, D.L. 1994 Hydrodynamic diffusion in a suspension of sedimenting point particles with periodic boundary conditions. Phys. Fluids 6, 2894.
• Ladd, A.J.C. 1993 Dynamical simulation of sedimenting spheres. Phys. Fluids 5, 299.
• Ladd, A.J.C. 1997 Sedimentation of homogeneous suspensions of non-Brownian spheres. Phys. Fluids 9, 491.
• Lee, S., Jang, Y., Choi C. & Lee,T. 1992 Combined effect of sedimentation velocity fluctuation and self-sharpening on interface broadening. Phys. Fluids 4, 2601.
• Leighton, D.T., & Acrivos, A., 1987 The shear-induced migration of particles in concentrated suspensions. J. Fluid Mech., 181, 415.
• Loewenberg, M., & Hinch, E.J., 1996 Collision of two deformable drops in shear ow. J. Fluid Mech., 338, 299.
• Nicolai, H., Herzhaft, B., Hinch, E.J., Oger, L. & Guazzelli, E. 1995 Particle velocity fluctuations and hydrodynamic self-diffusion of sedimenting non-Brownian spheres. Phys. Fluids 7, 12.
• Nicolai, H. & Guazzelli, E. 1995 Effect of the vessel size on the hydrodynamic diffusion of sedimenting spheres. Phys. Fluids 7, 3.
• Saffman, P.G. 1973 On the settling speed of free and fixed suspensions. Stud. Appl. Math. 52, 115.
• Sangani, A.S. & 1982 Acrivos, A. Slow flow through a periodic array of spheres, J. Multiphase Flow 8, 343.
• Ségre, P.N., Herbolzeimer, E. and Chaikin, P.M. 1997 Long-range calculations in sedimentation. Phys. Rev. Lett. 79, 2574.
• Silva Rosa, O.L. 1998 Fluctuations and dispersion in sedimentation with vertical stratification. (in portuguese), M.Sc. dissertation, University of Brasília, Brasília, DF, Brazil.
• Silva Rosa, O.L. & Cunha, F.R. 1998 Fluctuations and dispersion in sedimentation with vertical stratification. "General papers in Fluid mechanics" (poster session paper 129) in AIChE Annual Meeting, Miami, USA and in Proceedings, 7th Brazilian Congress of Engineering and Thermal Sciences, Rio de Janeiro, Brazil.
• Van de Vorst, G.A.L. 1996 Integral formulation to simulate the viscous sintering of a two-dimensional lattice of periodic unit cells. Journal of Engineering Mathematics 30, 97.
• Xue, J.-Z., Herboltzheimer, E., Rutgers, M.A., Russel, W.B. & Chaikin, P.M. 1992 Diffusion, dispersion, and settling of hard spheres. Phys. Rev. Lett. 69, 1715.

Publication Dates

• Publication in this collection
  24 Sept 2002
• Date of issue
  2001

History

• Received
  Feb 2001