UFDC Home  myUFDC Home  Help 



Full Text  
SEDIMENTATION OF HARD SPHERE SUSPENSIONS AT LOW REYNOLDS NUMBER USING A LATTICEBOLTZMANN METHOD By NHANQUYEN NGUYEN A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2004 Copyright 2004 by NhanQuyen Nguyen This dissertation is for my parents, Peter and Hoa Thi Nguyen, for all their support. ACKNOWLEDGMENTS I would like to thank my advisor, Tony Ladd, for his guidance and understanding in giving me the opportunity to work on the sedimentation problem. The work was long and tedious, but his professionalism and hard work were an inspiration that kept me going. I would also like to thank my fellow colleagues in the C('i, ii,, Engineering department. With such good people around, it made the long hours of doing research bearable. I thank my lab members Jonghoon Lee, for ahlv , being helpful; Byoungjin Chun, for his congeniality; Berk Usta, for his generosity; and Murali Rangarv' i, for being a good guy. Thanks also go to Piotr Szymczak and Rolf Verberg, two postdocs whose infinite knowledge and wisdom I was fortunate to learn from. Finally, I thank my Mom, Dad, brothers, and sister for being a part of my life. TABLE OF CONTENTS Page ACKNOW LEDGMENTS ................................ iv LIST OF TABLES . . . . . . . . . vii LIST OF FIGURES ....... ................ ........ .... viii A B ST R A CT . . . . . . . . . xi CHAPTER 1 INTRODUCTION ................................. 1 1.1 Properties of a Sedimenting Suspension ................. 2 1.2 Fluid Dynamic Equations ......................... 3 1.3 Hydrodynamics at Low Reynolds Number ............... 5 1.3.1 LongRange Interactions and the Divergence Problem . 8 1.3.2 MultiParticle Interactions . . . . . 10 1.3.3 Particles Near Contact . . . . . 11 2 DESCRIPTION OF THE LATTICEBOLTZMANN MODEL . ... .13 2.1 Molecular Theory for Fluid Flow . . . . . 13 2.1.1 Boltzmann Equation . . . . . . 14 2.1.2 Connection between Microscopic and Macroscopic Dynamics 15 2.2 LatticeBoltzmann Method ........ . . . 17 2.2.1 Gas Kinetics to Fluid Dynamics: The C!i inp ,iEnskog Ex pansion . . . . . ... 22 2.2.2 SolidFluid Boundary Conditions ..... . . 30 2.2.3 Hydrodynamic Radius of a Particle .... . . 33 2.2.4 Particle M otion . . . . . . 37 3 INCORPORATING LUBRICATION INTO THE LATTICEBOLTZMANN METHOD . . . . . 40 3.1 Introduction . . . . . . . . 40 3.1.1 Surfaces Near Contact . . . . . . 42 3.1.2 Lubrication Forces . . . . . ... 43 3.1.3 Particle Wall Lubrication . . . . . 44 3.1.4 Cluster Implicit Method . . . . . 48 3.2 Conclusions . . . . . . . . 52 4 SEDIMENTATION OF HARDSPHERE SUSPENSIONS AT LOW REYNOLDS NUMBER . . . . . .... 53 4.1 Introduction . . . . . . . .... 53 4.1.1 Sedimentation Simulations . . . . . 56 4.1.2 Settling of Pairs of Particles . . . . . 57 4.2 Velocity Fluctuations . . . . . . . 59 4.2.1 TimeDependent Velocity Fluctuations . . ... 59 4.2.2 SteadyState Settling . . . . . . 63 4.2.3 Numerical Errors . . . . . . 66 4.3 Microstructure . . . . . . . . 68 4.3.1 Particle Concentration . . . . . 69 4.3.2 Structure Factor . . . . . . 71 4.3.3 Mean Settling Speed . . . . . . 75 4.4 Conclusions . . . . . . . . 76 5 SEDIMENTATION OF A POLYDISPERSE SUSPENSION . . ... 77 5.1 Introduction . . . . . . . . 77 5.1.1 Stratification . . . . . . . 77 5.1.2 Velocity fluctuations . . . . . . 79 5.1.3 Structure Factor . . . . . . 81 5.1.4 Stratification in Laboratory Experiments . . .. 82 5.2 Conclusions . . . . . . . . 84 6 SUM M ARY . . . . . . . . . 85 REFEREN CES . . . . . . . . . 88 BIOGRAPHICAL SKETCH . . . . . . . 94 LIST OF TABLES Page 21 Variance in the computed drag force pF = V< F2 > < F >2/ < F > for a particle of radius a moving along a random orientation with respect to the grid . . . . . ..... . .... . 35 22 Hydrodynamic radius ahy (in units of Ax) for various fluid viscosities; a is the input particle radius . . . . . . . 35 31 Lubrication ranges for various kinematic viscosities, determined for a sphere of radius a = 8.2Ax. The ranges are determined separately for the nor mal, hN, tangential, hT, and rotational, hR, motions. . . ..... 46 41 Steady state particle velocity fluctuations in a monodisperse suspension as a function of system width . . . . . . 64 51 Steadystate particle velocity fluctuations in a 10W'. polydisperse suspen sion as a function of system width. . . . . . 79 LIST OF FIGURES Page 21 The 18 velocity vectors (ci) in the latticeBoltzmann model. . .... 17 22 Location of boundary nodes for a curved surface (a) and two surfaces in near contact (b). . . . . . . . .... 29 23 Moving boundary node update. The solid circles are fluid nodes, squares are boundary nodes with the associated velocity vector, and open circles are interior nodes . . . . . . . . 31 24 Drag force F as a function of time, normalized by the drag force on an isolated sphere, Fo 6= 67aU. Time is measured in units of the Stokes tim e, ts = /U . . . . . . . . 34 25 Actual (solid lines) and hydrodynamic (dashed lines) surfaces for a parti cle settling onto a wall . . . . . . . 36 31 Normal force on a particle of input radius a settling onto a horizontal pla nar surface at zero Reynolds number, with F0 67,/.,, U. Solid symbols indicate: v* = 1/6 (circles), v* = 1/100 (triangles), and v* = 1/1200 (squares). Results including the lubrication correction are shown by the open sym bols . . . . . . . . 45 32 Tangential force on a particle settling next to a vertical planar surface at zero Reynolds number. Simulation data for drag force divided by Fo=0 C,/., U, is indicated by solid symbols. Results including the lubrication correction are shown by the open symbols. . . . . . . 46 33 Torque on a particle settling next to a vertical planar surface at zero Reynolds number. Simulation data for torque divided by To87,/.,? U, is indicated by solid symbols. Results including the lubrication correction are shown by the open symbols . . . . . . . 47 34 Torque on a particle rotating next to a vertical planar surface at zero Reynolds number. Simulation data for torque divided by To87/'~., Q, is indicated by solid symbols. Results including the lubrication correction are shown by the open symbols . . . . . . . 48 35 Settling of a sphere (a 4.8) onto a horizontal wall. The gap between the particle surface and wall, h, relative to the hydrodynamic radius, is plot ted as a function of the nondimensional time (open circles). . . 49 36 Illustration of the algorithm to determine the list of clusters. . . 50 37 The maximum cluster size as a function of the cluster cutoff gap hs/a at varying volume fractions, . . . . . . 51 38 Number of clusters as a function of the cluster cutoff gap hs/a at varying volume fractions, . . . . . . . .... 51 41 Settling of a horizontal pair of particles at different Reynolds numbers; Re ~ 0.1 (circles), 0.05 (squares), 0.025 (triangles). . . ..... 58 42 Settling of a pair particles at various orientations: 0 = 0 (circles), 450 (squares), 680 (triangles), 900 (diamonds). The Reynolds number Re = 0.1 in each case . . . . . . . . 59 43 Particle velocity fluctuations as a function of time with a with W = 48a. 60 44 Particle velocity fluctuations as a function of time for 'Box" and "Bounded" geometries at various widths: W/a = 16 (circles), 32 (squares), and 48 (triangles) respectively. . . . . . . .... 61 45 Mean settling velocity, (U\), and fluctuations in settling velocity, AU 2, as a function of time for different system heights, H. . . . 63 46 Particle velocity fluctuations at steady state. Profiles are shown for the "Box" system (solid symbols) at different widths: W/a 1= 6 (circles), 32 (squares), and 48 (triangles) . . . . . . 64 47 Steady state particle velocity fluctuations as a function of system width for simulation versus experimental fit (solid line). . . . 65 48 Effect of inertia on the mean settling velocity and fluctuations in settling velocity . . . . . . . . . 66 49 Effect of grid resolution on the mean settling velocity and fluctuations in settling velocity . . . . . . . . 67 410 Effect of compressibility on the spatial variation of the mean settling ve locity. The data is for Reynolds number of Re = 0.06. . . . 68 4 11 Particle volume fraction as a function of height, z. The dots indicate the instantaneous average, and the steadystate (1000 < t/ts < 1400) density profile in the viewing window is shown in the .,.li 'ent plots. . . 69 4 12 Structure factor describing horizontal density fluctuations, S(k1), for var ious system sizes: W 16a, 32a, and 48a. ....... . . 71 4 13 Structure factor at different angles; vertical [0,0,1] direction (circles), 450 [1,0,1] direction (squares), and horizontal [1,0,0] direction (triangles). 73 414 Time evolution of structure factor for horizontal density fluctuations. 74 51 Profiles of the particle volume fractions as a function of the system height. 78 52 Comparison of particle volume fraction profiles for monodisperse and poly disperse suspensions at steadystate, 1000 < t/ts < 1400. . . 79 53 Particle velocity fluctuation as a function of time for a 10'. polydisperse suspension. The data is taken for three system widths: W/a 1= 6 (cir cles), 32 (squares), 48 (triangles). . . . . . ... 80 ix 54 Comparison of particle velocity fluctuations for monodisperse and 10'. polydisperse suspensions . . . . . . 81 55 Comparison of structure factors for different degrees of polydispersity: a) monodisperse (circles) and 1(,'. polydisperse (squares) suspensions (W=48a, N=72000); b) 2'., (circles), 5'U (squares), and 10't. (triangles) polydis perse suspensions (W 32a, N 32000). . . . . . 82 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy SEDIMENTATION OF HARD SPHERE SUSPENSIONS AT LOW REYNOLDS NUMBER USING A LATTICEBOLTZMANN METHOD By NhanQuyen Nguyen August 2004 Chair: Anthony J. Ladd Major Department: Chemical Engineering Particle motion in a settling suspension is directly influenced by hydrodynamic interactions, such that the particle velocities will exhibit a random component, char acterized by a hydrodynamically induced diffusion. Statistically, this means that the particles tend to be uniformly distributed within the bulk. Theoretical predictions for the sedimentation velocity in homogeneous suspensions are in good agreement with ex periments. On the other hand, a calculation of the particle velocity fluctuations, based on the same assumptions, is in conflict with experimental results. Our main objective was to understand this discrepancy. Numerical simulations have shown that dynamics within the bulk are not indepen dent of the size of the container, even when the boundaries are located far away from the observation window. We hypothesized that density fluctuations in the bulk drain away to the suspensionsediment and suspensionsupernatent interfaces. This generates a continual decay in velocity fluctuations, which is replenished by hydrodynamic diffu sion due to short range multiparticle interactions. A balance between convective and diffusive transport of density fluctuations leads to a correlation length beyond which the velocity fluctuations are screened. Numerical results were found to support this hypothesis. An examination of the microstructure as the suspension settled showed that it evolves toward a statistically nonrandom state, so that particles are distributed very uniformly at large length scales. We identified this result based on the observed damping of the structure factor, S(k) k2 as k goes to zero. Our numerical simulations show that the longrange correlations found in monodis perse suspensions are destroyed by small amounts of polydispersity, which are found in typical laboratory experiments. The variation in particle size also leads to a segregation of different sizes as the suspension settles. The net effect is a pronounced stratification of particle concentration at the supernatent front that persists into the bulk. Because of the stratification, the velocity fluctuations within the bulk do not have to convect toward the density gradients at the front; rather, the gradient is strong enough inside the bulk to create a sink for the velocity fluctuations. CHAPTER 1 INTRODUCTION We studied gravity driven nonBrownian particles within a viscous fluid. We deal with lowReynolds number flow, such that the dynamics of particles in the dispersion is strongly affected by hydrodynamic interactions generated by the relative motion of particles and fluid. Owing to the complexity of suspension dynamics, in 1985, a paradox for steady state sedimentation was pointed out by Caflisch and Luke (CL) (1). Consider a steadily sedimenting suspension of hard spheres. A concentration fluctuation gives rise to a velocity field decaying as 1/r with distance r from the origin of the fluctuations. The linearity of Stokes flow implies that the velocity field from many spatially distributed concentration fluctuations is simply a sum Ecu, of the individual contributions. If these concentration fluctuations take place in a random, spatially uncorrelated manner throughout the suspension, the resulting variance Eu2 in the velocity at any point in the suspension would be the sum of the squares of the individual contributions. This sum Eiuf has N ~ L3 terms with N solute particles in a container of linear dimension L in all directions. The 1/r contribution mentioned above leads to a velocity fluctuation of < ui >~ L2, so that ESu2 ~ L. The CL prediction is based on two simple physical arguments: The fluid flow is slow so that the Reynolds number Re The concentration fluctuations are statistically independent from one point to another in space. However, most experiments do not find the size dependence predicted by CL; the reason for the disagreement is unclear. One hypothesis is that a sufficiently strong anticorrelation of density fluctuations develops at large length scales, which suppress the CL divergence. We present findings from latticeBoltzmann simulations to verify this hypothesis. Simulations of the sedimentation problem require largescale numerical calcula tions, for which the latticeBoltzmann method is a useful tool. The simulation is based 1 on the development of microscopic interactions that lead to macroscopic fluid flow. This is done by applying the Boltzmann equation developed from gas kinetic theory. The simulation uses a spatial grid for the fluid phase, making it an inherently parallizeable algorithm where simulations of largescale systems spread over multiple processors can be performed. With boundary conditions at the surface of the particles, the lattice Boltzmann method becomes applicable to multiphase flows. We are able follow the trajectories of each particle in the suspension over time, giving a detailed picture of the particle dynamics. 1.1 Properties of a Sedimenting Suspension Sedimentation deals with systems that exist in a nonequilibrium steady state. The particles are driven through the solvent by a density mismatch, and have a downward velocity resulting from a balance between gravity and viscous drag force. This leads to a phase separation over time, which develops a region of pure fluid on top, and a 1. r of sediment at the bottom. The shrinking domain of bulk suspension is separated from the pure fluid by a horizontal interface, moving downward, and from the sediment by another horizontal interface, moving upward. Ultimately, these interfaces meet, the motion ceases, and separation is complete. Further observation of a settling suspension shows that particles exhibit random motion even when the particles are large and Brownian motion is negligible. The flow developed by each particle influencing the other makes the dynamics chaotic, and highly sensitive to initial conditions. The resulting chaos implies that the time evolution of coarsegrained quantities must be described using diffusion coefficients and noise sources, even though the microscopic dynamics in the absence of thermal Brown ian motion are entirely deterministic (2, 3). This induced large scale diffusive behavior, in the absence of thermal noise, is called hydrodynamic diffusion or hydrodynamic dispersion. Important macroscopic questions arise in the study of this ":ipl!, prototype system. For instance: what is the time of separation? What is the distribution of velocities and concentration? How is this process affected by the container boundaries and dimensions? 1.2 Fluid Dynamic Equations The equations of motion of a fluid are in essence a reformulation of Newton's second law (which states that the force exerted on a body of mass m equals m times the accelerations a) (4). For a fluid volume V, the mass density and acceleration can vary from place to place. The volume V is therefore divided into infinitesimally small volumes dV over which p and a are constant. Hence for an arbitrary volume AV F ma= padV. (1.1) JAV The force F exerted on a fluid volume AV consists of two distinct contributions: A body force exerted on the total volume. A surface force exerted on the boundary A of volume AV. Hence F = Fbody + Fsurf. (1.2) Common body forces are gravity and external electric or magnetic fields. Defining feat as the body force per unit volume, we can write Fbody j fxtdV. (1.3) JAv The surface forces can be found by considering the stresses exerted on the boundary A of the volume AV: F,,,f I dA (1.4) where II is the stress tensor that specifies the local normal and tangential forces per unit area exerted on a fluid element; and dA is a vector pointing perpendicular to the surface, its magnitude determined by the surface area of the infinitesimal element dA. Using the divergence theorem, Eq. 1.4 can be written as F8f= V HdV. (1.5) J.2, the equation of motion for a fluid can be written as From Eqs. 1.1, 1.3, and 1.5 into 1.2, the equation of motion for a fluid can be written as pa = V H + f(6t. (1.6) In specifying the acceleration term, an Eulerian description of the velocity u(r, t) is presented. Hence, Du au a + u Vu. (1.7) Dt at The two terms on the right hand side can be interpreted as the local rate of change in velocity (Ou/Ot) due to changes at position r and a convective rate of change in velocity (u Vu) due to the transport of the element to a different position (4). The value Vu is a secondorder tensor describing the gradient of the velocity field. It can be split into symmetric and antisymmetric parts as Vu 1 (Vu + Vu) + (Vu VuT) (1.8) where VuT is the transpose of Vu. Equation 1.8 can then be rewritten as Vu S + A (1.9) where S is the rateofstrain tensor that describes fluid deformation, and A is the vorticity tensor that corresponds to a solidbody rotation (4). The stress tensor H must still be defined. In the absence of fluid flow, the stresses in a fluid are caused by the static fluid pressure, which is isotropic (i.e., H PI (1.10) where P is the hydrostatic pressure and I is the unit tensor). However, when fluid flow occurs, the stress tensor becomes anisotropic, and can be expressed as the sum of an isotropic pressure P (usually different from the static pressure) and a nonisotropic contribution. For incompressible Newtonian liquids, this nonisotropic contribution is given by 2pS, where p is the viscosity of the fluid. This proportionality between the stress and the rate of strain is a phenomenological relationship, which many fluids have been found to obey (4). Hence H PI + 2pS. (1.11) Taking the divergence (and using VuT 0 for an incompressible fluid) gives V. VP + pV2u. (1.12) Substituting Eq. 1.7 and 1.12 into Eq. 1.6 gives the equation of motion for the velocity field p o + u Vu ) VP + V2+ ft, (1.13) which together with the continuity equation for an incompressible fluid V u 0, (1.14) are referred to as the NavierStokes equations. The NavierStokes equations can be written in nondimensional form by introducing a characteristic length scale I and velocity uo. Defining the dimensionless variables i uot/1, r~ r/1l, ii u/uo, P IP/puo and fet = fetl2/PUo (1.15) then substituting them into Eq. 1.13 gives S+ i Vu) VP + V2U e + t, (1.16) where the tildes represent nondimensional parameters and operators. The constant developed on the LHS is defined as the Reynolds number, Re = lUop (1.17) /t and can be interpreted as a measure of the ratio of inertial to viscous forces. 1.3 Hydrodynamics at Low Reynolds Number The standard problem in hydrodynamics is to find the velocity of a particle on which forces are acting. The simple example of a sphere moving under the action of a constant external force through a quiescent fluid can be considered as a model problem. As the isolated object settles in the fluid, it creates a disturbance flow in its surrounding neighborhood which becomes zero at large distances. After a short time, the sphere will reach a steadystate velocity, U. To determine U, we must find the fluid flow field u around the sphere, and calculate the stress exerted on the sphere surface by the fluid. When the Re NavierStokes equation reduces to the Stokes or creepingflow equation. In dimensional form, Eq. 1.13 can be rewritten as VP + pV2u ft (1.18) with V u 0. (1.19) A solution to the Stokes equation, Eq. 1.18, can be developed using Fourier transforms (5, 6). First, the Fouriertransform pair is defined in real space r and Fourier space k as u(k) J eikru(r)dr (1.20) u(r) J e ikrf(k)dk. (1.21) Taking the Fourier transform of Eq. 1.18 yields ikP + pk2 f F(k), (1.22) ku 0. (1.23) Taking the dot product of k on the transformed equation of motion, Eq. 1.22, elimi nates the u term to give P(k) i k (1.24) Substituting P into Eq. 1.22 gives u(k) = 2 F F k) (1.25) To invert the transformed solutions (Eqs. 1.24 and 1.25) we use the following relation f27 Fk A F V (1.26) (27 )3 j k2 47r 1 F eik.rdk = F 1 (1.27) (27)3 k2 4 7r 1 Fkk F.(1 rr \ and (2 ikdk F (1.28) (27)3 k4 0y 8 yrr3 From Eqs. 1.24 and 1.26, the pressure field is determined to be 1 P(r) 3r fet, (1.29) 47r3 with fext representing the force from a point source. The velocity field is determined by taking the inverse Fourier transform of Eq. 1.25 with Eqs. 1.27 and 1.28 to give u(r) (I + 3 fet. (1.30) 87/1 r r 3 Hence, the Oseen tensor is given as 0(r) rr. (1.31) O(r) (I + 3t 87 \r r3 We now look at the case when the full velocity field is not needed, but only the force on a particle is desired. Consider a rigid sphere of radius a and surface S immersed in a fluid with velocity field u,(r) in the absence of the sphere. The sphere is centered at ro, translates with velocity Uo and rotates about its center with angular velocity f. Owing to its motion, the sphere affects the fluid through forces distributed over its surface, H(r) per unit area. From Eq. 1.30 the velocity field at r is now uo(r) plus that generated by the superposition of forces H(r')dr' (5), or u(r) u,(r) f 0(r r') H(r')dr'. (1.32) JS At the surface of the sphere, this must match the sphere velocity given by U = Uo + Qx (r ro). Thus, Uo + Q x (r ro) u,(r) 0(r r') H( lr')dr'. (1.33) JS Integration of Eq. 1.33 over the sphere surface results in IH(r')dr' = F [uo(r) Uo] dr. (1.34) Js 2a Js Expanding u, (r) about the sphere center for creeping flow, and integrating over S yields Faxen's first law F 67/a (uu + ua 2oU)0 o (1.35) The subscript zero indicates evaluation at ro, the particle center. This is an important result, showing that in a nonuniform flow, the particle does not move with the mean fluid velocity (uo) in the absence of an external force, but that there is a correction proportional to a27V2 U. Faxen's second law can be derived by taking the vector product of Eq. 1.33 with (r ro) before carrying out the second integration, and then separating the equation into antisymmetric and symmetric parts, which must be satisfied individually (5). In this way the torque on a sphere (7) is obtained as T 87pa3 r(V x u)o 2O (1.36) and the force dipole (stresslet) as S a3 + 2V) (VuC + Vu). (1.37) From Eq. 1.35, Stokes's law for the force exerted on a sphere settling in a infinite medium is obtained. By setting Uo = 0 or u. = 0 we get F 6= 6wau. or = 67raUo, (1.38) and the torque is T 87/Sa3Q. (1.39) 1.3.1 LongRange Interactions and the Divergence Problem To obtain the flow disturbance caused by the presence of a sphere, Eq. 1.32 must be evaluated. By setting II = F/47ra2 3pUo/2a and expanding the Oseen tensor about the center of the sphere (5) leaves u(r) 2a Uo 0. O(R) + (r' ro) VO(R) + (r' ro)(r' ro) VVO(R) +... dr' = 6paUo ( + a2V2 0(R) (1.40) with R = r ro, and Uo representing the Stokes settling velocity for an isolated sphere. Equation 1.40 is derived from the identities / dr' = 47a2 (1.41) (r' ro)(r' ro)dr' = a4 I, (1.42) Js 3 and the remaining terms integrate to zero, since i(r' ro)'dr' = 0 for n odd, (1.43) V 0 0, and V2V20 0. Substituting the Oseen tensor (Eq. 1.31) into Eq. 1.40 gives the velocity field (5) 3a a2 3a a2) r Uor 44) u ) = 1+ Uo 1+ 1+ (144) 4R 3R2 ) 4R R2 R2 which satisfies the boundary conditions u Uo, R a (1.45) u 0, R + oc. (1.46) Therefore, the disturbance is long range, decaying as 1/R. Now we add one more sphere at a position r2 from ri so that iri r2 12 > a. To first order from Eq. 1.44, Sphere 1 feels the flow caused by Sphere 2 with the following strength Ul12 3 (1 + R12R12) Uo + ... (1.47) 4R12 where the vector R12 is the unit vector R12/R12. The velocity of Sphere 1 is now increased with respect to the single sphere sedimentation velocity Ui Uo + (1 + R12R12) U+ (1.48) 4R12 If the same reasoning is applied to the sedimentation of a cloud of N identical spheres, then to first order in the inverse of the interparticle distance, we get N U Uo + (1 + RjRi) U + ... (1.49) Because the disturbance in the fluid due to particle interactions is long range, the summation in Eq. 1.49 will diverge as the volume of the system becomes large. The divergence problem was solved by Batchelor (8), who pointed out that there is an upcurrent of displaced fluid that had not been considered. This backflow cancels the longrange hydrodynamic interactions, and causes the settling velocity of a suspension to be less than that of an isolated particle. 1.3.2 MultiParticle Interactions With the linearity of the Stokes equation presented in Sect. 1.3, the development of the resistance and mobility relations for a single particle moving in a fluid can be obtained, where F j (TT. U, (1.50) T =(.RR Q. (1.51) The forces F and torques T exerted on the particle by the fluid are related to the particle velocities U and angular velocities Q via the configuration dependent friction coefficients ( (6, 9), which are developed from the solution of the Stokes equations. For multiparticle interactions, the above equations become N sF < . Upu + C . ] (1.52) j1 N Tj [C U, + (1.53) j 1 where the summation extends over all N particles in the suspension, and C is the resistance matrix. The superscripts T and R refer to the translational and rotational components of C, with TR representing the coupling between the two components. Equations 1.52 and 1.53 can be inverted to find the corresponding mobility matrices tp, relating the particle velocities to the hydrodynamic forces and torques, N U [> T. Fj + pT T] (1.54) j1 N o. L[ Fy + pR Tj] (1.55) j 1 Thus, the resistance and mobility matrices completely characterize the linear relations between the force and torque with the translational and rotational velocities for particle motion at low Reynolds number flow. 1.3.3 Particles Near Contact For the case when particles come in close proximity with one another, lubrication theory is applied to provide the leading order terms in the resistance and mobility matrix. Lubrication theory is calculated based on the motion of two surfaces near contact and the forces are pairwise additive in this limit, F1 All A12 B11 B12 Ui F2 A21 A22 B21 B22 U2 p (1.56) TI Bil B12 C1 C12 1I T2 / 21 B22 C21 C22 2 The square resistance matrix contains second rank tensors A, B, and C. The elements of the resistance matrix obey a number of symmetry conditions. The reciprocal theorem states that the resistance matrix is symmetric (6, 9), so that y = ji (1.57) Bi Bj, (1.58) C3 C (1.59) With e = R/R being the unit vector along the line of centers, we can write (6) S = + Y(J ee1), (1.60) Y B =Y kek, (1.61) CJ = X + Y i , .). (1.62) 12 The form of these scalar functions are developed, which is a function of the dimension less separation a/R. The dimensionless functions for diagonal elements of the resistance matrix will approach unity for large R because the scales were chosen by considering the single sphere result (6). From Eqs. 1.57 1.62, we obtain 10 independent non dimensional scalar functions to describe the resistance matrix. Analytic approximations to these scalar functions can then be derived from lubrication theory for particles in near contact. CHAPTER 2 DESCRIPTION OF THE LATTICEBOLTZMANN MODEL The latticeBoltzmann model is based on the Boltzmann equation developed from gas kinetic theory. Kinetic theory is the branch of statistical physics dealing with the dynamics of nonequilibrium processes and their relaxation to thermodynamic equilibrium. The Boltzmann equation, established by Ludwig Boltzmann in 1872, is its cornerstone. 2.1 Molecular Theory for Fluid Flow According to the molecular theory of matter, a macroscopic volume of gas ( I 1 cm3) is a system of a very large number ( Iv, 1020) of molecules moving in a rather irregular way. In principle, we may assume that the molecules are particles obeying the laws of classical mechanics. It is also assumed that the laws of interaction between the molecules are perfectly known so that, in principle, the evolution of the system is computable, provided suitable initial data are given. However, solving the initial value problem for a number of particles of a realistic order of magnitude ( '*, N 1020) is an impossible task. As a result, only averages can be computed and are related to macroscopic quantities as pressure, temperature, stresses, heat flow, etc.; this is the basic idea of statistical mechanics. From statistical mechanics, we talk about probabilities instead of certainties; that is, a given particle will not have a definite position and velocity, but only probabilities of having different positions and velocities. Under suitable assumptions, the information required to compute averages for these systems can be reduced to the solution of an equation, the socalled Boltzmann equation. With a solution to the Boltzmann equation, we can then deduce macroscopic fluid flow from the microscopic model. 2.1.1 Boltzmann Equation The kinematic equation for the onebody distribution function (10) reads as follows: Dtf [= t + P r + F. p] f(r, p,t) = C12, (2.1) where the lefthand side represents the streaming motion of the molecules along the trajectories associated with the external force field F and C12 represents the effect of intermolecular (twobody) collisions taking molecules in/out of the streaming trajectory. Here r is the position of the particle, p = mc its linear momentum, and F is the force experienced by the particle as a result of intermolecular interactions and possibly external fields. The collision operator involves the twobody distribution function f12 that ex presses the probability of finding a molecule, t 1, around ri with velocity Cl and a second molecule, 2, around r2 with velocity c2, both at time t. The dynamic equation for C12 can be developed, but because the collisions are coupled this equation calls into p1 i, the threebody function C123, which in turn depends on C1234 and so on down an endless line known as the BBGKY hierarchy, after Bogoliubov, Born, Green, Kirkwood, and Yvon (11). To simplify Eq. 2.1, Boltzmann assumed that the system is a dilute gas of point like structureless molecules interacting via a shortrange twobody potential. Under such conditions, intermolecular interactions can be described solely in terms of localized binary collisions, where molecules spends most of their time on free trajectories perfectly unaware of each other. With this picture, the collision term splits into gain (G) and loss (L) components: C12 G L = (fl'2' f12)gg, Q)ddp2 (2.2) corresponding to direct/inverse collisions taking molecules out/in the volume ele ment (10). The parameter a is the collision cross section, which expresses the number of molecules with relative speed gQ around the solid angle Qf. Then comes Boltzmann's closure assumption: f12 f2. (2.3) This closure assumes no correlations between molecules entering a collision (molecular chaos or Stosszahlansatz) and is the essential component that describes the Boltzmann equation. The Boltzmann equation then takes the form: Lat + P B Or + Fext p] f(r, p, t) = (f, f2 fif2)g(g, Q)dQdp2. (2.4) The left hand side is Newton's equation for single particle dynamics, while the right hand side describes intermolecular collisions under the Stosszahlansatz approximation. 2.1.2 Connection between Microscopic and Macroscopic Dynamics The mass density p(r, t) is the integral of the density in the oneparticle phase space f(r, c, t) with respect to all possible velocities p(r, t)= Jmfdc, (2.5) where m is the molecular mass. Here, the probability distribution function f(r, c, t) describes the number of molecules at time t, with position r lying within the velocity space c. Because of the probabilistic meaning of f, the density p is the expected mass per unit volume at (r, t) or the product of the molecular mass m by the probability density of finding a molecule at (r, t), which turns out to be the number density n(r, t) n(r, t) = p(r,t)/m. (2.6) From Eq. 2.5, the momentum pu can be written as follows pu = cmfdc (2.7) or, using the following components: ', = cimfdc. (2.8) The velocity u can be perceived as the macroscopic observation developed from molecular motion; it is zero for the steady state of a gas enclosed in a box at rest. Each molecule has its own velocity c which can be decomposed into the sum of u and another velocity v c u, (2.9) which describes the random deviation of the molecule from the ordered motion with velocity u. The velocity v is usually called the peculiar velocity or the random velocity; it coincides with c when the gas is macroscopically at rest (12). Thus, we have from Eqs. 2.5, 2.8, and 2.9, S, ,, fdc = cimfdc i mfdc ,", /', = 0. (2.10) Another quantity needed for developing macroscopic hydrodynamics is the mo mentum flux. Since momentum is a vector, we therefore consider the flow of the jth component of momentum in the ith direction; given by I ci(cjmf)dc = cicnjmfdc. (2.11) Eq. 2.11 shows that the momentum flux is described by a symmetric tensor of second order. It is known that in a macroscopic description only a part of the microscopically evaluated momentum flux will be identified as convective, because the integral in Eq. 2.11 will be in general different from zero even if the gas is macroscopically at rest. In order to find out how the momentum flux will appear in a macroscopic description, we have to use the splitting of c into mean velocity u and peculiar velocity v. This is written as where Eq. 2.5 and 2.10 have been applied. The momentum flux therefore breaks into two parts, one of which is recognized as the macroscopic momentum flux, while the second part is a hidden momentum flux due to the random motion of the molecules. For example, in a fixed region of gas, we observe that the change in momentum is only attributed to the matter that enters and leaves the region, which means the second part has no macroscopic explanation unless it is attributed to the action of a force exerted / \ \\ I ^ / 0^  S/ Figure 21: The 18 velocity vectors (ci) in the latticeBoltzmann model. on the boundary of the region of interest by the contiguous region of the gas. In other words, the integral of f, .,,, fdc appears as the contribution to the stress tensor (12) pij ,= .mfdc. (2.13) Here, the trace pai = ljpjj is the isotropic part of the stress tensor. More specifically, the fluid pressure can be written as P pii/3 for the case at equilibrium. 2.2 LatticeBoltzmann Method The computational utility of latticeBoltzmann models depends on the fact that only a small set of velocities are necessary to simulate the N i, rStokes equations (13). Thus the more general kinetic equations in Sec. 2.1 are simplified by considering a discrete set of velocities only. The lattice used in this simulation is a three dimensional cubic lattice for which there are 18 velocities, Fig. 21. This model uses the [100] and [110] directions, with twice the density of particles moving in the [100] directions as in [110] directions. This gives the symmetry condition necessary to ensure that the hydrodynamic transport coefficients are isotropic. In this work the 18velocity model is augmented with stationary particles, which enables it to account for small deviations from the incompressible limit, although in simulations of stationary flows it is found that the numerical differences are small (14). In the latticeBoltzmann method (LBM), probabilities substitute for discrete particles in the particle velocity function ni(r, t). The LBM then describes the time evolution of the discretized velocity distribution ni(r + ciAt, t + At) = ni(r, t) + A, [n(r, t)], (2.14) where Ai is the change in ni due to molecular collisions at the lattice nodes, and At is the time step. This oneparticle distribution function describes the mass density of particles with velocity ci, at a lattice node r, and at a time t; r, t, and ci are discrete, whereas ni is continuous. The hydrodynamic fields, mass density p, momentum density j = pu (where u is the fluid flow velocity), and momentum flux H, are moments of this velocity distribution: P = ii, j =pu = niCi, H = niCiCi. (2.15) i i i The postcollision distribution, nu(r, t) + Ai [n(r, t)], is propagated for one time step, in the direction ci. The collision operator Ai(n) depends on all the ni's at the node represented collectively by n(r, t), with the constraints of mass and momentum conservation. Exact expressions for the Boltzmann collision operator have been derived for several latticegas models (15, 16) by making use of the ii,,. I. ii (iI ... assumption (that the effect of the collision operator depends only on the instantaneous state of a node). However, such collision operators are complex and illsuited to numerical simulation. A computationally more useful form for the collision operator Ai(n) can be constructed by linearizing about the local equilibrium(17) neq Ai(n) A(neq) + inj, (2.16) where Lij are the matrix elements of the linearized collision operator, and Aj(neq) = 0. L is calculated by considering the general principles of conservation and symmetry and the construction of the eigenvalues and eigenvectors of L. The population density associated with each velocity has a weight ac' that de scribes the fraction of particles with velocity ci in a system at rest; these weights depend only on the speed ci and are normalized so that .,"" 1. The velocities ci are chosen such that all particles move from node to node simultaneously. For any cubic lattice, ac Cicic C2c2I, (2.17) where c = Ax/At, Ax is the grid p' ii: and C2 is a numerical coefficient that depends on the choice of weights. However, in order for the viscous stresses to be independent of direction, the velocities must also satisfy the isotropy condition; ac CiaCis3CyCi6 = C 4C4 {a36 + 6ay66 + 6a63}. (2.18) The isotropy condition requires that a multispeed model be applied, which in this case is the 18velocity model for a simple cubic lattice. As a first step, the correct form of the equilibrium distribution must be established. To determine the equilibrium distribution, the velocity distribution function is split into a local equilibrium part and a nonequilibrium part i nq + neq. (2.19) The equilibrium distribution is a collisional invariant (i.e., A,(nq) = 0). The form of the equilibrium distribution is constrained by the moment conditions required to reproduce the inviscid (Euler) equations with nondissipative hydrodynamics on large length and time scales. In particular, the second moment of the equilibrium distribution should be equal to the inviscid momentum flux PI + puu: p = nq (2.20) j =ci c pu (2.21) 11eq __E c = pc I + puu (2.22) The equilibrium distribution can be used in Eqs. 2.20 and 2.21 (c.f. Eq. 2.15) because mass and momentum are conserved during the collision process, thus nnne neqci _c 0. (2.23) i i The pressure in Eq. 2.22, P= pc takes the form of an ideal gas equation of state with adiabatic sound speed cs. This is also adequate for the liquid phase if the density fluctuations are small (i.e. the Mach number M = u/c8 < 1), so that VP = c Vp. For small Mach numbers, the equilibrium distribution can be expanded as a power series in u ci/c and it then follows from Eq. 2.22 that the weights must be chosen so that C2 (c/c)2; i.e. from Eq. 2.17 Saccc c2I. (2.24) A suitable form for the equilibrium distribution of the 19velocity model that satisfies Eqs. 2.20 2.22, as well as the isotropy condition (18) (Eq. 2.18), is S j ac pPuu:+ (2c I) (2.25) where c = c c = Ax/At, P = pc2, and the coefficients of the three speeds are a 1 1 .1 (2.26) 3 18' a 36(2.26) In this case the coefficient in Eq. 2.18 is C4 = (c/c)4. From Eqs. 2.25 and 2.26, the inviscid hydrodynamic equations are correctly repro duced. The viscous stresses come from moments of the nonequilibrium distribution, as in the (C'!i pn ,Enskog solution of the Boltzmann equation. The fundamental limitation of this class of latticeBoltzmann model is that the Mach number be small, less than 0.3; our suspension simulations ahv keep M < 0.1. Having constructed the local equilibrium distribution, we come to the relaxation of the nonequilibrium component of the local distribution. Beyond the conservation of mass and momentum inherent in neq, the eigenvalue equations of the collision operator for neq develop the isotropic relaxation of the stress tensor. The linearized collision operator must satisfy the following eigenvalue equations; Lj =0, Zc4 = 0, c=c4 Acc, c Lj AcA, (2.27) i i i i where cFc indicates the traceless part of cici. The first two equations follow from conservation of mass and momentum (c.f. Eq. 2.23), and the last two equations describe the isotropic relaxation of the stress tensor; the eigenvalues A and A, are related to the shear and bulk viscosities and lie in the range 2 < A < 0. Equation 2.27 accounts for only 10 of the eigenvectors of L. The remaining 9 modes are higherorder eigenvectors of L that are not relevant to the N ',1i rStokes equations, but which do affect the boundary conditions at the solidfluid interfaces. In general the eigenvalues of these kinetic modes are set to 1, which both simplifies the simulation and ensures a rapid relaxation of the nonhydrodynamic modes (19). The collision operator can be further simplified by taking a single eigenvalue for both the viscous and kinetic modes (20, 18). This exponential relaxation time (ERT) approximation, A = nq /T, has become the most popular form for the collision operator because of its simplicity and computational efficiency. However, the absence of a clear time scale separation between the kinetic and hydrodynamic modes can sometimes cause significant errors at solidfluid boundaries, and thus the more flexible collision operator described by Eq. 2.27 is preferred. A 3parameter collision operator is used in the suspension simulations, allowing for separate relaxation of the 5 shear modes, 1 bulk mode, and 9 kinetic modes. The postcollision distribution n* = ni + Ai is written as a series of moments (Eq. 2.15), as seen with the equilibrium distribution (Eq. 2.25), *a( + j.Ci_ + (puu + flneq,*) 2:( Cc (2.28) The zeroth (p) and first (j = pu) moments are the same as in the equilibrium distri bution (Eq. 2.25), but the nonequilibrium second moment HIeq is modified by the collision process, according to Eq. 2.27: neq,* (1 + A)H + ( + A,)(Ineq : 1)1, (2.29) 3 where nueq II Ieq. The kinetic modes can also contribute to the postcollision distribution, but the eigenvalues of these modes are set at 1, so that they have no effect on rn. Equation 2.29 with A = A, = 1 is equivalent to the ERT model with T = 1; for A < 1, the kinetic modes relax more rapidly than the viscous modes, which is the proper limit for hydrodynamics. The macrodynamical behavior arises from the latticeBoltzmann equation with the multitimescale analysis applied (15, 14). The NavierStokes equations, tp + V. (pu) 0(2.30) t(pu) + V (puu) + Vpc pV2u + (v+ V(V. ), are recovered in the low velocity limit, with viscosities p pc2At( + and / pc At ( + 1). (2.31) The factors of 1/2 and 1/3, which appear in the definitions of the shear and bulk viscosities serve to correct for numerical diffusion, so that viscous momentum diffuses at the expected speed. In the presence of an externally imposed force density f, for example a pressure gradient or a gravitational field, the time evolution of the latticeBoltzmann model includes an additional contribution f/(r,t), ni(r + ciAt, t + At) i(r, t) + Ai [n(r, t)] + fi(r, t). (2.32) This forcing term can be expanded in a power series in the velocity, f f'JQcr i + (uf + fu) :(cce cI)j At( fi ac c2 2c4 )t) (2.33) The expression relating the force density to the change in velocity distribution was determined by matching the macroscopic dynamics from Eq. 2.32 to the NavierStokes equations (14). More accurate solutions to the velocity field are obtained if it includes a portion of the momentum (14) added to each node, j' pu' = nbce + fAt. (2.34) 2.2.1 Gas Kinetics to Fluid Dynamics: The ChapmanEnskog Expansion In order to solve the Boltzmann equation introduced in Sec. 2.1 for realistic nonequilibrium situations, we must rely upon approximation methods, in particular, perturbation procedures. In order to do this, we look for a parameter C which we consider to be small. Therefore, in obtaining macroscopic hydrodynamic equations, we introduce the nondimensional Knudsen number denoted by Kn Ku = l/ 1, (2.35) where I is the mean free path between molecular collisions, and d is a typical length scale. Thus, Kn > 0 corresponds to a fairly dense gas and Kn o0 to a free molecular flow (i.e., a flow where molecules have negligible interactions with each other). The C'!i p1in ,Enskog expansion, then, considers multiple time scales, of orders e and c2, to recover the N ,',i, rStokes equation. In this sense, the N ,',i, rStokes equation describes two kinds of processes, convection and diffusion, which act on two different time scales (21). If we consider only the first scale, we obtain the compressible Euler equation; if we move on to the second one we obtain the NavierStokes equation while losing the compressibility effect. If we want both compressibility and diffusion, we have to keep both scales at the same time and think of n in terms of time and space derivatives n + 2a (2.36) at at1 Ot2 On Ce (2.37) Or OrI This enables us to introduce two different time variables t1 = et, t2 = C2t and a new space variable ri = er, such that the fluid dynamical variables are functions of ri, t1, t2,. The compatibility conditions at the first order give that the time derivatives of the fluid dynamic variables with respect to t1 is determined by the Euler equation but the derivatives with respect to t2 are determined only at the next level and are given by the terms of the NavierStokes equation describing the effect of viscosity and heat conductivity (21). The two contributions are to be added as indicated in Eq. 2.36 to obtain the full time derivative and thus write the NavierStokes equation. While the C'!i Ipn :Enskog expansion was originally applied to the continuous Boltzmann equation in Eq. 2.1, here we derive the NavierStokes equation using the simplified lattice Boltzmann equation represented as follows: n1(r + c1At, t + At) n1(r, t) = Q(r, t), (2.38) with Q(r, t) [ (r, t) n (r, t)]. (2.39) T The collision operator here is from the single relaxation time BhatnagarGrossKrook (BGK) model where naq represents the equilibrium distribution and 7 is a typical timescale associated with collisional relaxation to the local equilibrium. To evaluate the LBGK model, first a Taylor series expansion for ni(r + cAt, t + At) is performed, ni(r + ciAt, t + At) = n(r, t) + At + c At + at ar, (At)2 a2 2n, 02n, t)2 a + 2caa, + ciGcli9 (2.40) Next we introduce the two time scales and one spatial scale as noted in Eq. 2.36 and 2.37. The distribution function is also expanded accordingly: n nio) + Cni + C2 n2) + 0(c3). (2.41) Substitute Eq. 2.40 into the LHS of Eq. 2.38 with Eq. 2.39 gives ani an (At)2 a2n, a2ni _2ni At + CiaAt + + 2ci, + CiaCi3 at ar, 2 at2 ata9rar3 1 +(ni nq) 0. (2.42) 7 Substitute Eqs. 2.36, 2.37, and 2.41 into Eq. 2.42 gives CAt at + Ca + C2At [i + a + Cia a + (2.43) at, O8ia at, Ot2 9rla 2 2 (A2 (0) 2 (0) 2n() e + 2ci + ciciClaa + (2.44) 2 9t 8tj OtjO8Ia 8 9rI9rj S(no) n ) + Cn + 2 2) 0. (2.45) T T T Arranging terms in ascending orders of c gives 0(e) o) 0(e1): (1) 0(2) (2) ATt E at1 TAt [ant [ at1 (At)2 2 + cia Orl,, (+1) Cia  9rla (0) Ot2 a2n(O) atlat1 2 (o) + 2ci ia Oti1r1 2 (0) ] + CiIC 'Orl 8r08ri/3 It is noted that the equilibrium distribution, along with Eq. 2.46, should satisfy the following constraints: S m p = i n i i=0 i0 m PUa = Ciani i=0 m can i=0 m E (0) in0 i=o) i= (2.49) (2.50) m E (0) Cian i, i=0 where m is the total number of directions in the lattice model, and the molecular mass is set to unity. With these constraints applied to Eq. 2.41, the constraints for the 1st and 2nd order components in c become 0n "UO mT 0 i 0an i=0 Tnl i (2) i U (2.51) (2.52) m i (2) Cian i i=0 Using the above relation in Eq. 2.47 gives the following reduction + Cia O a m 0 o) i=0 Op at, 9pua + 0rla a m (0) i 0 0. Furthermore, the second moment is obtained by Z Ci3 at, + Cia Ol] 0 0+ > i=0 9t, ,(0) + a (9r l o (2.46) (2.47) (2.48) [an40 (2.53) ,(0) 0 (2.54) Recall from the derivation in Eq. 2.12 using the peculiar velocity, the viscous stress in terms of the distribution function n can be determined to be i=0o P( = PuU3 + PSa3. (2.55) Thus, to first order in c, from Eq. 2.53 the continuity equation is obtained and from Eqs. 2.54 and 2.55 the Euler equation is determined au'3 au,3 OP pt +p ua  (2.56) at rla rla Next, the 0(c2) terms are analyzed. From Eq. 2.48, the relations in Eqs. 2.49, 2.51, 2.53, and 2.54 are applied to get a relation in terms of t2, and rn (2) (1) (0) (1)' 0 + a + cia ao + (2.57) i =0 i =0 m A 2 (0) n2 (0) 2 (0) TnAt 2n(O) 2n(O) 12n' I ___0lT _+ 2Cia + Cia Ci ail 2 talar Otarri rar3 Then i=0 a m 9t2 z 0) ra cm Substitution of Eqs. 2.53 and 2.54 into the 4th, ates and overall value of zero. Thus, the zeroth (2.58) (2.59) (2.60) 5th, and 6th terms of Eq. 2.57 gener moment of 0(c2) becomes ap 0. Ot2 For the momentum conservation, Eq. 2.52 is applied to give m m 1) Tn (0), (1) S(2) 0 __, + + + i=0 i=0 m CAt n2 (0) a2n 0) 2n (0) 0 2 ia t1it1 + 2ci t + ciaCi3 i  i0 2 atl8rl OtiOrio rlaOrl (2.61) (2.62) Here, applying Eqs. 2.50 and 2.52 gives at2  aput at2 Applying Eq. 2.47 into the 3rd component of Eq. 2.62 gives a T Or i0 . a m arla i=0 IE 1 at1 + cn  O7 r ,, The 4th component in Eq. 2.62 can be evaluated further by substitution to give At m 2  iO a2n(7) ci at at8 81 At 02pUa 2 at1at1 Ata ( 2 at1 OPri Br~i ) At a a + Sat, ar (PUaU + pc~ ) At 2 0 Op 2 r ri3 ati where 0(u2) terms are neglected due to low Mach number flows. Here c, represents the speed of sound which is characterized by the pressure at equilibrium (i.e P pRT pc2). Substitution from Eq. 2.53 gives At 2 a a c2 Orla r(pu) 2 90 r^ a ~ A2 (2.69) The 5th component of Eq. 2.62 is considered in conjunction with the 1st term in Eq. 2.65 to get (1 raa 1 ()u Otria Ori3 (1 T)At ap(O) t Or l 3 a/3 (t T) AtV7,7. (PU) (2.70) (2.63) (2.64) (2.65) (2.66) (2.67) (2.68) and the 6th component of Eq. 2.62 with the second term in Eq. 2.65 gives T Ata a icpc ) (using defn. from Eq. 2.25) \~ 2} dr^GOrjir3 _0 i=0 ( 1 a a m) c c^ c c(cispuK) (2.71) t iAlc a(pu3) + 2 (pu3) 2 21 099r1 0aria 0rl/3 ) Atc, [V2(pu) + 2VV (pu)] (2.72) The new term in Eq. 2.71 is taken from the istropic condition for this lattice model. Finally combining Eqs. 2.64, 2.69, 2.70, and 2.72 into Eq. 2.62 gives t(pu) ( ) Atc [V2(pu) + VV. (pu)] (2.73) Thus the dynamic shear and bulk viscosities for the lattice BGK model are PS P= =B 7 t) Atc. (2.74) The end result is to sum over all orders of c. With Eq. 2.53 added to Eq. 2.61, the continuity equation is determined: a + a = 0. (2.75) at Ora Eqs. 2.54, with Eq. 2.55, added to Eq. 2.73 gives a + u r3 OrP + t Atc [V2(pu) + VV. (pu)] (2.76) In the incompressible limit, Ou3 0 (2.77) 89r3 8pu a 9pua SP au pu3 + pV2(pu). (2.78) at ar3 Ora Thus, the continuum N ,'. i, Stokes equation has been derived. In the end, the basic result of the Ch '!pi iiEnskog procedure is that the macroscopic description of a fluid is recovered by a suitable expansion of the Boltzmann equation. 0 0 0 o * O 0 O 0 a 0 o 0 0 o o o o o Da b Figure 2 2: Location of boundary nodes for a curved surface (a) and two surfaces in near contact (b). 2.2.2 SolidFluid Boundary Conditions Boundary conditions in the latticeBoltzmann model are straightforward to implement, even for nonplanar surfaces (19). Solid particles are defined by a surface, which cuts some of the links between lattice nodes. Fluid particles moving along these links interact with the solid surface at boundary nodes placed halfway along the links. Thus we obtain a discrete representation of the particle surface, which becomes more and more precise as the particle gets larger (Fig. 22). The motion of a particle is determined by the forces and torques acting on it. These forces consist of external, nonhydrodynamic forces, and hydrodynamic forces. Nonhydrodynamic forces, such as attraction and/or repulsion between surfaces due to an interparticle potential, can be calculated separately and incorporated into the equations of motion. In contrast, the hydrodynamic forces are calculated from the resultant velocity distributions, ni. Stationary solid objects are introduced into the simulation by the "bounce ., .I: collision rule at the specified boundary nodes, in which incoming distributions are reflected back towards the nodes they came from. Surface forces are calculated from the momentum transfer at each boundary node and summed to give the force and torque on each object (19). While more accurate methods that require knowledge of the shape of the particle surface have been investigated (22, 23, 24, 25, 26, 27), the bounceback rule can be applied to surfaces of arbitrary shape, without additional complications. However, for secondorder convergence the bounceback rule requires a calibration of the hydrodynamic radius, which is not alvb, convenient. The boundary nodes are located midway between interior (solid) and exterior (fluid) nodes (28, 29). The normal collision rules are carried out at all fluid nodes, and augmented by bounceback rules at the midpoints of links connecting the lattice nodes. This is depicted by the arrows in Fig. 22a linking the normal distribution between fluid (circles) and boundary nodes (squares). In Poiseuille flow the "link bounce '.,. 1 rule gives velocity fields that deviate from the exact solution by a constant slip velocity, Us 8 ULBE UExact, proportional to L2 (30) (L being the channel width); u8/ue Ax2/L2, (2.79) Ub 0 t+ 01  Figure 23: Moving boundary node update. The solid circles are fluid nodes, squares are boundary nodes with the associated velocity vector, and open circles are interior nodes. where c = L2 Vp /8v is the exact velocity at the center of the channel and the coefficient 3 depends on the eigenvalues of the collision operator. In the past the lattice nodes were treated on either side of the boundary surface in an identical fashion, so that fluid fills the whole volume of space, both inside and outside the solid particles. Because of the relatively small volume inside each particle, the interior fluid quickly relaxes to rigidbody motion, characterized by the particle velocity and angular velocity, and closely approximates a rigid body of the same mass (14). However, at short times the inertial lag of the fluid is noticeable, and the contribution of the interior fluid to the particle force and torque reduces the stability of the particle velocity update. For these reasons, with the si:. I .i' in Ref. (31), the interior fluid is removed from the particle representation. A somewhat different implementation of the same idea is described in Ref. (32). The moving boundary condition (19) without interior fluid (31) is implemented as shown in Fig. 23. We take the set of fluid nodes r just outside the particle surface, and for each node all the velocities Cb such that r + cbAt lies inside the particle surface. Each of the corresponding population densities is then updated according to a simple rule which takes into account the motion of the particle surface (19); 2acb pub. Cb nb(rt+At) n(r,t) (2.80) Cs where n4 (r, t) is the postcollision distribution at (r, t) in the direction cb, and cb' Cb. For a stationary boundary the population is simply reflected back in the direction it came from with ub = 0 (15, 33). The local velocity of the particle surface, Ub U U+ Qx(rb R), (2.81) is determined by the particle velocity U, angular velocity fl, and center of mass R; rb r + CbhAt is the location of the boundary node. The mean density po in Eq. 2.80 is used instead of the local density p(r, t) since it simplifies the update procedure. The differences between the two methods are small, of the same order (pu3) as the error terms in the latticeBoltzmann model. Test calculations show that even large variations in fluid density (up to 10'.) have a negligible effect on the force (less than 1 part in 104). As a result of the boundary node updates, momentum is exchanged locally between the fluid and the solid particle, but the combined momentum of solid and fluid is conserved. The forces exerted at the boundary nodes can be calculated from the momentum transferred in Eq. 2.80, f(rb,t + 'At) A3 2n(r, t) 2a oUb Cb, (2.82) The particle forces and torques are then obtained by summing f(rb) and rb x f(rb) over all the boundary nodes associated with a particular particle. It can be shown analyt ically that the force on a planar wall in a linear shear flow is exact (19), and several numerical examples of latticeBoltzmann simulations of hydrodynamic interactions are given in Ref. (34). To understand the physics of the moving boundary condition, one can imagine an ensemble of particles, moving at constant speed Cb, impinging on a massive wall oriented perpendicular to the particle motion. The wall itself is moving with velocity Ub < Cb. The velocity of the particles after collision with the wall is cb + 2ub and the force exerted on the wall is proportional to Cb Ub. Since the velocities in the latticeBoltzmann model are discrete, the desired boundary condition cannot be implemented directly, but we can instead modify the density of returning particles so that the momentum transferred to the wall is the same as in the continuous velocity case. It can be seen that this implementation of the noslip boundary condition leads to a small mass transfer across a moving solidfluid interface. This is physically correct and arises from the discrete motion of the solid surface. Thus during a time step At the fluid is flowing continuously, while the solid particle is fixed in space. If the fluid cannot flow across the surface there will be large artificial pressure gradients, arising from the compression and expansion of fluid near the surface. For a uniformly moving particle, it is straightforward to show that the mass transfer across the surface in a time step At (Eq. 2.80) is exactly recovered when the particle moves to its new position. For example, each fluid node .,dli ,. ent to a planar wall has 5 links intersecting the wall. If the wall is advancing into the fluid with a velocity U, then the mass flux across the interface (from Eq. 2.80) is poU. Apart from small compressibility effects, this is exactly the rate at which fluid mass is absorbed by the moving wall. For sliding motion, Eq. 2.80 correctly predicts no net mass transfer across the interface. 2.2.3 Hydrodynamic Radius of a Particle Although the linkbounceback rule at the solidfluid interface is second order accurate for planar walls oriented along lattice symmetry directions, it is only first order accurate for channels oriented at arbitrary angles (35, 36). Thus for large channels, the hydrodynamic boundary is displaced by an amount A from the physical boundary, where A is independent of channel width but depends on the orientation of the channel with respect to the underlying lattice. Convex bodies sample a variety of boundary orientations, so that it is not possible to make an analytical determination of the displacement of the hydrodynamic boundary from the solid particle surface. Nevertheless, the displacement of the boundary can be determined numerically from simulations of flow around isolated particles. By considering the size of the particles to be the hydrodynamic radius, ahy = a + A, rather than the physical radius a, approximate secondorder convergence can be obtained, even for dense suspensions (34). The hydrodynamic radius can be determined from the drag on a fixed sphere in a pressuredriven flow (34). Galilean invariance can be demonstrated by showing that the same force is obtained for a moving particle in a stationary fluid (34). Since the 100 120 140 160 180 200 1.420 a = 8.2.7 a = 8.5 1.41500 100 120 140 160 180 200 t/t, Figure 24: Drag force F as a function of time, normalized by the drag force on an isolated sphere, Fo = 67/aU. Time is measured in units of the Stokes time, ts = a/U. particle samples different boundary node maps as it moves on the grid, it is important to sample different particle positions when determining the hydrodynamic radius, especially when the particle radius is small (< 5Ax). Rather than averaging over many fixed configurations, we chose to have the particle move slowly across the grid, at constant velocity, sampling different boundary node maps as it goes. The changing boundary node maps lead to fluctuations in the drag force, as shown in Fig. 24. The force has been averaged over a Stokes time, ts, so that the relative fluctuations in force are comparable to the relative fluctuations in velocity of a neutrally buovint particle in a constant force field. The force fluctuations, 6F = /(< F2 > < F >2)/ < F >, are of the order of one percent for particles of radius 2 3Ax, and are considerably smaller for larger particles (Table 2 1). More sophisticated boundary conditions have been developed using finitevolume methods (37, 38) and interpolation (23, 25, 39). Both methods reduce the force fluctuations by at least an order of magnitude from those observed here, but even with the simple bounceback scheme, the fluctuations in force can be reduced by an appropriate choice of particle radius. We have noticed 35 Table 21: Variance in the computed drag force 6F = V< F2 > < F >2/ < F > for a particle of radius a moving along a random orientation with respect to the grid. a/Ax 2.7 2.5 8.2 8.5 6p 5.738e03 1.208e02 4.332e04 5.674e04 that fluctuations in particle force are strongly correlated with fluctuations in particle volume. Thus we choose the radius of the boundary node map so as to minimize fluctuations in particle volume for random locations on the grid. It can be seen from Table 21 that a two fold reduction in the force fluctuations is possible by this procedure, although for sufficiently large particles the difference is minimal. A set of optimal particle radii is given in Table 22. The bounceback rule leads to a velocity field in the region of the boundary that is firstorder accurate in the grid spacing Ax. The hydrodynamic boundary (the surface where the fluid velocity field matches the velocity of the particle) is displaced from the particle surface by a constant, A (Fig. 25), that depends on the viscosity of the fluid (34). For the range of kinematic viscosities used in this work, 1/6 < v* < 1/1200, A varies from 0 to 0.5Ax (Table 22); the dimensionless kinematic viscosity v* = vAt/Ax2. For small particles (a < 5Ax), A also depends weakly on the particle radius (Table 22). Although the difference between the hydrodynamic boundary and the physical boundary is small, it is important in obtaining accurate results with computationally useful particle sizes. Calibration of the hydrodynamic radius leads to approximately secondorder convergence from an inherently firstorder boundary condition; it will not be necessary when more accurate boundary conditions are implemented. The hydrodynamic radii (ahy) in Table 22 were determined from the drag force on a single sphere in a periodic cubic cell (34). The Reynolds number in these simulations Table 22: Hydrodynamic radius ahy (in units of Ax) for various fluid viscosities; a is the input particle radius. a/Ax 1.24 2.7 4.8 6.2 8.2 v* = 1/6 1.10 2.66 4.80 6.23 8.23 v* = 1/100 1.42 2.92 5.04 6.45 8.46 v* = 1/1200 1.64 3.18 5.31 6.72 8.75 Figure 25: Actual (solid lines) and hydrodynamic (dashed lines) surfaces for a particle settling onto a wall. (< 0.2) was sufficiently small to have a negligible effect on the drag force. The time averaged force was used to determine the effective hydrodynamic radius for three different kinematic viscosities: v* = 1/6, v* = 1/100, and v* = 1/1200. In each case the cell dimensions were 10 times the particle radius and the corrections for periodic boundary conditions (about !n' .) were made as described in Ref. (34). The difference between the hydrodynamic radius and the input radius, A = ahy a, is independent of particle size for radii a > 5Ax, and its magnitude increases with decreasing kinematic viscosity (14). The kinematic viscosity v* = 1/6 gives a hydrodynamic radius that is the same as the input radius (for sufficiently large particles), and is the natural choice for simulations at low Reynolds number. At higher Reynolds numbers a lower viscosity is necessary to maintain incompressibility (14, 34), and for accurate results it is then essential to use the calibrated hydrodynamic radius. In order to implement the hydrodynamic radius in a multiparticle suspension, all distance calculations are based on the hydrodynamic radius (as shown in Fig. 25); the input radius a is only used to determine the location of the boundary nodes. It should be noted that not all combinations of particle radius and viscosity can be used. Table 22 indicates that particle radii less than 3Ax cannot be used with a kinematic viscosity v* = 1/6, since the hydrodynamic radius is then less than the input radius. 2.2.4 Particle Motion An explicit update of the particle velocity U(t + At) U(t) + F(t) (2.83) [(t + At) Q(t) + T(t) (2.84) I has been found to be unstable (34) unless the particle radius is large or the particle mass density is much higher than the surrounding fluid. In previous work (34) the instability was reduced, but not eliminated, by averaging the forces and torques over two successive time steps. Subsequently, an implicit update of the particle velocity was proposed (40) as a means of ensuring stability. Here we present a generalized version of that idea, which can be adapted to situations where two particles are in near contact. The particle force and torque can be separated into a component that depends on the incoming velocity distribution and a component that depends, via ub, on the particle velocity and angular velocity (Eqs. 2.81 and 2.82); F Fo (FU F (2.85) T To CTU. U CT T (2.86) The velocity independent forces and torques are determined at the halftime step Fo(t + At) At Z 2n(r,t)cb (2.87) b To(t + tAt) At 2n(r, t)(rb x Cb), (2.88) b where the sum is over all the boundary nodes, b, describing the particle surface, with cb representing the velocity associated with the boundary node b and pointing towards the particle center. The matrices FU 2p x3 acbcbCb (2.89) c t b F 2po Ax3 c 2At aCbcb(rb X Cb) (2.90) b TU 2po Ax3 acb(rb X Cb)Cb (2.91) b cAt 3 cb(rb X Cb)(rb X Cb) (2.92) b are highfrequency friction coefficients, and describe the instantaneous force on a particle in response to a sudden change in velocity. The magnitude of the friction coefficients can be readily estimated, thereby establishing bounds on the stability of an explicit update. Apart from irregularities in the discretized surface, (FU and TQo are diagonal matrices, while (FQ TU = 0. For a node .'Id ',ent to a planar wall Z, aci2 c2 where the sum is over the 5 directions that cross the wall. The number of such nodes is approximately 477a2/Ax2, so that (FU 207 poAxa2 (2.93) 9 At Similarly, TQ 87 poAxa4 (2.94) 9 At These estimates of the translational and rotational friction coefficients are within 21 i. and 50'. of numerically computed values, respectively. The stability criterion for an explicit update (FUAt/m < 2 then reduces to a simple condition involving the particle radius and mass density: 5 < 2. (2.95) 3 psa The corresponding condition for the torque leads to the same stability criterion S 5 pfAx < 2, (2.96) I 3 psa whereas with interior fluid the numerical factors were 6 and 10 respectively (34), showing that interior fluid destabilizes an explicit update. The friction coefficients in Eqs. 2.85 and 2.86 are essentially constant, fluctuating slowly in time as the particle moves on the underlying grid; thus the particle velocities can be updated assuming these friction matrices are constant. The equations of motion can then be written in matrix form as 1 U(t + At) U(t) + aC(F aCF S(t + At) (t) aC" 1+a cT Fo FU() CFQ(t) (2.97) To (TuU(t) (CT (t) where a is a parameter that controls the degree of implicitness. An explicit update (34) corresponds to a = 0, an implicit update (40) corresponds to a = 1, and a second order semiimplicit update corresponds to a = The explicit, implicit, and semiimplicit updates evaluate the velocitydependent force at t, t + At, and t + At respectively. In practice we find only small differences between semiimplicit and fully implicit methods and we use the fully implicit method (a = 1) in this work. The boundary node map is updated infrequently (every 10100 time steps) and the 6 x 6 matrix inversion need only be done when the map is updated. We note that in the limit of (FUAt/m >> 1 only the fully implicit (a = 1) version of Eq. 2.97 reduces to the steady state force and torque balance, Fo FU(t + At) To CT (t + At) 0. The semiimplicit method (a = ) produces an oscillating solution and the explicit method (a = 0) a diverging solution. An implicit update of the particle velocities requires two passes through the boundary nodes. On the first pass the population densities are used to calculate Fo and To. The implicit equations (Eq. 2.97) are then solved for U(t + At) and [2(t + At) for the given implicit parameter a. These velocities are then used to calculate the new population densities in a second sweep through the boundary nodes. The computational overhead incurred by the boundary node updates is typically less than 10i i. even at high concentrations. CHAPTER 3 INCORPORATING LUBRICATION INTO THE LATTICEBOLTZMANN METHOD 3.1 Introduction LatticeBoltzmann simulations (19, 34) are being increasingly used to calculate hydrodynamic interactions (31, 41, 42, 43, 44, 45, 46), by direct numerical simulation of the fluid motion generated by the moving interfaces. However, as two particles approach each other the calculation of the hydrodynamic force breaks down in the gap between the particles, typically when the minimum distance between the two surfaces is of the order of the grid spacing. For example, it is impractical to use sufficient mesh resolution to resolve the flow in the small gaps that can result from an imposed shear flow. At high shear rates the rheology of a suspension of spheres is qualitatively affected by lubrication forces between particles separated by gaps less than O.Ola, where a is the particle radius (47, 48). A simulation at this resolution (t 107 grid points per particle) is unfeasible for more than a few particles. The number of grid points can be reduced by using body fitted coordinates (49) or adaptive meshes (23, 25), but the small zone size in the gap reduces the time step that can be used to integrate the flow field (50). The problem is exacerbated by the uniform grid used in latticeBoltzmann simulations, but it should be noted that similar techniques, using particles embedded in regular grids, are becoming increasingly popular in computational fluid dynamics (51, 52). Some aspects of this work may therefore be applicable to such methods as well. Simulations of hydrodynamically interacting particles typically use multiple expansions of the Stokes equations (53, 54). Again the calculations break down when the particles are close to contact, unless an impractical number of multiple moments are used (55). A solution to this problem is to calculate the lubrication forces between pairs of particles for a range of small interparticle gaps, and then patch these solutions on to the friction coefficients calculated from the multiple expansion. The method exploits the fact that lubrication forces can be added pairbypair, and has been shown to work well in practice (53). A simplified version of this approach has already been adopted to include normal lubrication forces in latticeBoltzmann simulations (56). Here we extend our previous work to include all components of the lubrication force and torque, and test the numerical scheme for the interactions between a spherical particle and a planar wall. We find that the hydrodynamic interactions can be well represented over the entire range of separations by patching only the most singular components of the lubrication force onto the force calculated from the lattice Boltzmann model. This is simpler than the Stokesian dynamics approach, where the patch is calculated at every separation. In this chapter, we propose a comprehensive solution to the technical difficulties that arise when particles are close to contact, in particular the loss of mass conserva tion. The bulk of our numerical results are a series of detailed tests of the hydrody namic interactions between two surfaces in near contact. We demonstrate that after including corrections for the lubrication forces we obtain accurate results over a wide range of fluid viscosities. Finally, we describe an efficient implicit algorithm for updat ing the particle velocities even in the presence of stiff lubrication forces. An explicit solution of these differential equations requires either that the particles are much denser than the surrounding fluid (34), or that very small time steps are used to update the particle velocities. On the other hand, if the particle velocities are updated implicitly, the resulting system of linear equations severely limits the number of particles that can be simulated. We describe what we call a "clusterimplicit" method, whereby groups of closely interacting particles are grouped in individual clusters and interactions be tween particles in the same cluster are updated implicitly, while all other lubrication forces are updated explicitly. In fluidized suspensions clusters typically contain less than 10 particles, even at high concentrations, so that the implicit updates are a small overhead. Our simulations efficiently cope with clusters of several hundred particles by using conjugategradient methods, and only slow down if the number of particles in the cluster exceeds 103. 3.1.1 Surfaces Near Contact When two particle surfaces come within 1 grid spacing, fluid nodes are excluded from regions between the solid surfaces (Fig 22b), leading to a loss of mass conserva tion. This happens because boundary updates at each link produces a mass transfer (2acbpoUbCb /C)Ax3 across the solidfluid interface, which is necessary to accommodate the discrete motion of the particle surface (see Section 2.2.2). The total mass transfer in or out of an isolated particle, AM = 2Axo U c abb C+ acbrb X Cb 0, (3.1) 8s b b regardless of the particle's size or shape. Although the sums Lb acbCb and Lb acbrb X Cb are zero for any closed surface, when two particles are close to contact some of the boundary nodes are missing and the surfaces are no longer closed. In this case AM / 0 and mass conservation is no longer ensured. Two particles that remain in close proximity never reach a steady state, no matter how slowly they move, since fluid is constantly being added or removed, depending on the particle positions and velocities. If the two particles move as a rigid body mass conservation is restored, but in general this is not the case. The accumulation or loss of mass occurs slowly, and in many dynamical simula tions it fluctuates with changing particle configuration but shows no longterm drift. However, we enforce mass conservation, particlebyparticle, by redistributing the excess mass among the boundary nodes (c.f. Eq. 2.80) n, (r, t + At) n (r, t) 2 abpU A (3.2) where A = AZ3 L acbpo. b The force and torque arising from this redistribution of mass is small, but not exactly zero; Ax3 0 AM b AF A [ A accb (3.3) At A b b Ax3po AM AT At A acbrb X Cb (3.4) b They can be compactly included by a redefinition of the friction coefficients, Eqs. 2.89 2.92, replacing Cb and rb x Cb by their deviation from the mean, a> bCb acbrb X Cb c Zb Yab and rb X Cb =acb (3.5) b b so that Cb > Cb c and rb X Cb rb X Cb C rb X Cb. Then the force and torque are correctly calculated from Eqs. 2.85 and 2.86, even when mass is redistributed. 3.1.2 Lubrication Forces When two particles are in near contact, the fluid flow in the gap cannot be resolved. For particle sizes that are typically used in multiparticle simulations (a < 5Ax), the lubrication breakdown in the calculation of the hydrodynamic interaction occurs at gaps less than O.la. However, in some flows, notably the shearing of a dense suspension, qualitatively important physics occurs at smaller separations, typically down to O.Ola. Here we describe a method to implement lubrication corrections into a latticeBoltzmann simulation. For particles close to contact, the lubrication force, torque, and stresslet can be calculated from a sum of pairwiseadditive contributions (53), and if we consider only singular terms, they can be calculated from the particle velocities alone (57). In the grandresistancematrix (58) formulation FI An Bil B22 T Bi C11 C12 U12 T2 B22 C12 C22 1Q (3.6) S GH Hi H12 f2 S2 G22 H21 H22 where F2 F1, U12 U1 U2, and the friction matrices are as given in Ref. (58). We have made full use of the symmetries inherent in the friction matrices, but without assuming that the particle radii are the same. Most importantly, the external flow field does not enter into the lubrication calculation, which removes the need for estimates of the local flow field. We have noted in previous latticeBoltzmann simulations (14, 56) that the calcu lated forces follow the Stokes flow results down to a fixed separation, around 0.5Ax, and remain roughly constant thereafter. The solid symbols in Fig. 31, for example, show this behavior for the normal force between a spherical particle and a plane wall. This ,. 1 a lubrication correction of the form of a difference between the lubrication force at h and the force at some cut off distance hN; i.e. Flub 67 l 1 U12 12 h = 0 h> hN, where U12 E UI U2, h R12 01 02 is the gap between the two surfaces, and the unit vector R12 R12/IR121 The friction coefficients in Eq. 3.6 are all products of a scalar function of the gap between the particles, either 1/h or In(1/h), multiplied by a polynomial of the unit vector connecting the particle centers (58). For two spheres of arbitrary size, there are a total of 15 independent scalar coefficients, which fall into 3 classes. Again using the notation of Ref. (58) these are XA XG XG (normal force); YA YB Y, Y2G2 (tangential force); and Ycl, Yc, yC, YH1, yH, H, yH (rotation). We implement our lubrication correction by calculating a modified form of each scalar coefficient as in Eq. 3.7; for example XA ) XA hXA N Xll(h) Xl(h) xIAI(hN) h < hN (3.8) XAl (h) 0 h>hN which vanishes at the cutoff distance h = hN. We allow for different cut off distances for each of the 3 types of lubrication interaction. 3.1.3 Particle Wall Lubrication The hydrodynamic interactions between two moving surfaces have been calculated for the simplest geometry, namely a spherical particle .,Ii :ent to a planar wall. We used 3 different particle sizes, with input radii a/Ax 2.7, 4.5, and 8.2, chosen to minimize volume fluctuations (see Section 2.2.2) with the exception of the results for a/Ax = 4.5, which were generated before the optimum radius (4.8Ax) was determined. The hydrodynamic radii, ashy(a, v*), which are used to determine the 1000 100 10 1 0.0 1000 100 10 1 0.0 1000 100 10 1 a =2.7 )01 0.01 0.1 1 a =4.5 a0.1 0.01 0.1A .. . 0.001 0.01 0.1 1 h/ah Figure 31: Normal force on a particle of input radius a settling onto a horizontal planar surface at zero Reynolds number, with Fo=67y/.,, U. Solid symbols indicate: v* = 1/6 (circles), v* = 1/100 (triangles), and v* = 1/1200 (squares). Results including the lubrication correction are shown by the open symbols. positions of the lubricating surfaces, were taken from Table 22. The location of the planar wall was shifted by A(v*), corresponding to the a > o limit in Table 22 (Fig. 25). In this way we ensure that the lubricating surfaces are in the same position as the hydrodynamic boundaries in the latticeBoltzmann simulations. The unit cell is periodic in four directions with a top and bottom wall, and cell length of 10 times the particle radius, which is sufficiently large that the effects of periodic images were negligible. The simulation determines the steady state force and torque on the particle for a given velocity and angular velocity, which was then used to calculate the friction coefficient as a function of the gap h from the wall. Simulation results for the frictional forces and torques are shown in Figs. 31 34 for three different fluid viscosities v* = 1/6, 1/100, and 1/1200. The normal force shows the trend discussed in section 3.1.2 for each particle size and fluid viscosity (Fig. 31). The simulated force (solid symbols) follows the exact result (solid line) down to an interparticle gap, hN < Ax, that is independent of particle C C C 5 4 3 2 1 0. 5 4 3 2 1 0. 5 4 3 2 1 a 2.7 001 0.01 0.1 1  a=4.5 001 0.01 0.1 1 a = 8.2 A* 0.001 0.01 0.1 1 h/a Figure 32: Tangential force on a particle settling next to a vertical planar surface at zero Reynolds number. Simulation data for drag force divided by Fo0 Cr,/. U, is indi cated by solid symbols. Results including the lubrication correction are shown by the open symbols. size. For larger particles the latticeBoltzmann method captures progressively more of the lubrication force, but even for a = 8.2Ax there are noticeable deviations for h/ahy < 0.01. The simulation reproduces more of the lubrication force at smaller viscosities because the shift in the hydrodynamic radius d, , the contact of the par ticle surfaces. The data obtained for a particle radius of 8.2Ax was used to determine the lubrication cutoff hN(v*) for each viscosity, and the numerical values are recorded in Table 31. These lubrication corrections bring the simulated normal force into agreement with theory for all the particle sizes, interparticle gaps, and fluid viscosities studied (open symbols in Fig. 31). The corresponding result for the force and torque Table 31: Lubrication ranges for various kinematic viscosities, determined for a sphere of radius a = 8.2Ax. The ranges are determined separately for the normal, hN, tangen tial, hT, and rotational, hR, motions. hN/Ax hr/Ax hR/Ax v* = 1/6 0.67 0.50 0.43 v* = 1/100 0.24 0.50 0.15 v* = 1/1200 0.10 0.50 0.00 0.6 0.5 a = 2.7 0.4  0.3 0.2 0.1 A 0 .. A . 0.001 0.01 0.1 1 0.6 0.5 a = 4.5  0.4 0.3 O 0.2  0.1 AU AA 0  0.001 0.01 0.1 1 0.6 ............. 0.5 a sphere sliding along the wall is shown in Figs. 32 & 33. Aga we see 8.2that the 0.4 0.3  0.2 0.1 SA 0.001 0.01 0.1 1 h/ah Figure 3 3: Torque on a particle settling next to a vertical planar surface at zero Reynolds number. Simulation data for torque divided by To 87,/.? U, is indicated by solid symbols. Results including the lubrication correction are shown by the open symbols. on a sphere sliding along the wall is shown in Figs. 3 2 & 3 3. Again we see that the lubrication correction gives consistently accurate forces and torques, though not quite as accurate as the normal force. The lubrication ranges for tangential motion were found to be independent of the fluid viscosity, as shown in Table 3 1. We also noticed that the reciprocal relations are obeyed; the force on a rotating sphere is similar to the data in Fig. 3 3. The calculated torque on a rotating sphere (Fig. 3 4) is in agreement with theory for the higher viscosities v* 1,/6 and v* 1/100, but not at the lowest viscosity V* 1/1200. Here the latticeBoltzmann method over predicts the torque on a rotating sphere by 20; '". We think that the error is caused by the large difference between the hydrodynamic and input radii (A 0.55Ax) and it implies that viscosities v* less than 0.01 are not suitable for suspension simulations, at least with bounceback boundary conditions. In practice this is not a serious limitation: a viscosity i* 0.01 4 ... a 2.7 1 0.001 0.01 0.1 1 4 a 45 0.001 0.01 0.1 1 4 . a 8.2 0.001 0.01 0.1 1 h/a Figure 34: Torque on a particle rotating next to a vertical planar surface at zero Reynolds number. Simulation data for torque divided by T 87',/., Q, is indicated by solid symbols. Results including the lubrication correction are shown by the open symbols. allows simulations with a Reynolds number up to 10 per grid point (with a Mach num ber ~ 0.1), which is at or beyond the limit of resolution of the flow. In other words, there is little practical value in viscosities less than 0.01. Finally, we examined the dynamic motion of a particle (a = 4.8Ax) settling onto a solid wall (Fig. 35). The lubrication force was calculated using the ranges given in Table 31. The particle was given a finite mass and placed with an initial gap of h = 0.2ahy between the particle and wall. The simulations were performed at a Reynolds number Re ~ 0.02 by applying a constant force to the particle. The results show the expected exponential decay of the gap between the particle and wall over time, in quantitative agreement with lubrication theory (59) (shown by the solid line). 3.1.4 Cluster Implicit Method The lubrication forces complicate the update of the particle velocity because they involve interactions between many particles, especially at higher concentrations. For simplicity we update the particle velocities in two steps; first the latticeBoltzmann 10 0 10 102 103 v*= 1/6 10 4 .. ... .. 0.1 1 10 100 1000 10 0 10 102 103 V*= 1/100 104 .. .. .. .. .. ... 0.1 1 10 100 1000 vt/a2 Figure 35: Settling of a sphere (a 4.8) onto a horizontal wall. The gap between the particle surface and wall, h, relative to the hydrodynamic radius, is plotted as a func tion of the nondimensional time (open circles). forces and torques (Eq. 2.97), followed by the lubrication forces. The overall procedure is still first order accurate, but the lubrication forces can cause instabilities whenever the particles are in near contact. The instability is driven by the normal forces, and the stability criteria (At 67Ta2At 9 rTAt 2p < 2 (3.9) m i7psa3h 2 psah is violated when h is less than ~ O.lAx. It is impractical to solve all the equations implicitly, so we implemented an algorithm which uses an implicit update only where necessary. In our simulations we used the conservative criteria (At/m < 0.1. Schematically, we solve the coupled differential equations c Ax+b (3.10) by splitting the dissipative matrix A into regular and singular components, A = AR + As. In our context As only contains components of the normal friction coefficient a % 9 9 32 b C Figure 36: Illustration of the algorithm to determine the list of clusters. when the gap between particles is less than the stability cut off, hs, determined from Eq. 3.9. Thus AR contains all the nonzero components of the lubrication correction but with the interparticle separation in the normal force regularized by h8 so that the larger of hij and h, is used to calculate the force. The remaining normal force is included in As, with the form of Eq. 3.7, but with hN replaced by h,. Using a mixed explicitimplicit differencing, (t + At) (t) AR. x(t) As x(t + At) + b, (3.11) At we obtain the first order update (1 + AsAt) x(t + At) x(t) ARAt x(t) + bAt. (3.12) The important point is that, by a suitable relabeling of the particle indices, As can be cast into a block diagonal form, with the potential for an enormous reduction in the computation time for the matrix inversion. The relabeling is achieved by a 1000 , 100 = 0.50 )= 0.25 S= 0.10 10 1 I I I I 0 0.02 0.04 0.06 0.08 0.1 0.12 hs/a Figure 37: The maximum cluster size as a function of the cluster cutoff gap hs/a at varying volume fractions, Q. cluster analysis. First, all pairs of particles that are closer than the stability cut off are identified, and a list is made of all such pairs. An illustration is shown in Fig. 36a, where pairs of particles with separations less than h8 are indicated by the solid lines. The cluster labels are initialized to the particle index; each particle is then relabeled by giving it the smallest label of all the particles it is connected to. After 1 pass the labels are as shown in Fig. 36b and after 2 passes 3 distinct clusters have been identified, each with a unique label (Fig. 36c). The algorithm stops when no further relabeling takes place. Although more efficient schemes are possible, this simple scheme is more than adequate for our purposes. Once the clusters have been identified, the implicit equations can be solved for each cluster. We use conjugate gradients to exploit the sparseness of As, which is extensive even within each diagonal block. 1000 i 100 S10 = 0.50 = 0.25 E B 4= 0.10 1 0 0.02 0.04 0.06 0.08 0.1 0.12 hs/a Figure 38: Number of clusters as a function of the cluster cutoff gap hs/a at varying volume fractions, Q. The computational cost of the clusterimplicit algorithm depends primarily on the maximum cluster size, which is shown in Fig. 37 as a function of the cluster cutoff gap hs. A random distribution of 1000 particles was sampled in a periodic box at volume fractions 0.10, 0.25, and 0.50. At low and moderate volume fractions the cluster size is only weakly dependent on hs/a, ranging from 27, and clusters of this size impose a negligible computational overhead. However, at high volume fractions there is a percolation threshold at hs/a ~ 0.02 beyond which a single cluster more or less spans the whole volume. In this case the cluster will grow to encompass almost all the particles in the system. Thus at high densities computational efficiency requires that h8/a < 0.02. When combined with the stability criteria, which implies hs/a w 1/a2, we find a minimum radius of a = 10Ax to keep iAt/m < 0.5. A simulation of several hundred such particles is possible on a personal computer or desktop workstation. In Fig. 38 we show the corresponding number of clusters. In general there is a steep rise in the number of clusters with increasing hs/a, leveling off to around 100 clusters. The sharp drop in the number of clusters at the highest volume fraction is associated with the percolation transition, as seen in Fig. 37. 3.2 Conclusions In this work we have described and tested several extensions to the lattice Boltzmann method for particle suspensions, which enable reasonably accurate force calculations to be made even for particles in near contact. In particular we have shown how to deal with problems of mass conservation when two particles are in near contact, and how to account for the lubrication forces between closelyspaced particles. Numerical tests show that the forces and torques between a particle and a plane wall can be computed to within a few percent of the exact result for Stokes flow. We note that the torque on a rotating sphere .,i d :ent to a plane wall is seriously in error (31i,) when the fluid viscosity is very small (1/1200). This ,i. I that the calibration procedure may break down when the hydrodynamic boundary is displaced by more than Ax/2 from the physical one. Inclusion of the lubrication forces leads to large forces and stiff differential equa tions for the particle velocities. We have developed a mixed explicitimplicit velocity 53 update, which minimizes the number of linear equations that must be solved, while maintaining absolute stability. CHAPTER 4 SEDIMENTATION OF HARDSPHERE SUSPENSIONS AT LOW REYNOLDS NUMBER 4.1 Introduction Particles larger than a few microns tend to settle out of suspension, because grav itational forces then dominate over the diffusive flux arising from gradients in particle concentration. The detailed dynamics of the most idealized flow, the sedimentation of hard spheres in the absence of inertia and Brownian motion is still controversial. When a suspension settles, each particle experiences a different shielding of the fluid di ,. due to the fluctuating arrangements of its neighbors. These hydrodynamic interactions are long range, decaying .,,iiii:l ,tically as 1/R, and drive large fluctuations in particle velocity, which, for particles more than about 10pm in diameter, completely dominate the thermal Brownian motion. A relatively simple calculation shows that, if the par ticle positions are independently and uniformly distributed, the velocity fluctuations are proportional to the linear dimension of the container (1). However, two different sets of experiments found that the velocity fluctuations converge to a fixed value for sufficiently large systems (60, 61). The qualitative discrepancy between theory and experiment has generated considerable attention, focusing on possible mechanisms for screening the hydrodynamic interactions (62) and thereby saturating the velocity fluctuations when the system size is larger than the hydrodynamic screening length. The key theoretical idea (62) is that hydrodynamic interactions can be screened by changes in suspension microstructure, analogous to the screening of electrostatic inter actions in charged systems. Hydrodynamic screening occurs when a test particle and its neighbors are, collectively, neutrally buo,,t with respect to the bulk suspension; in other words, when density fluctuations at length scales larger than the screening length are suppressed. This requires that a rather delicate longrange correlation develop in the distribution of particle pairs, which is resistant to the randomizing effects of hydro dynamic dispersion. Several different bulk mechanisms for the microstructural changes have been proposed, including threebody hydrodynamic interactions (62), a convective instability (63), and a coupled convectiondiffusion model (2, 64). However, these the ories are all inconsistent with the results of numerical simulations (42, 56, 65), which have shown that, in homogeneous suspensions (with periodic boundary conditions), particles are distributed randomly at separations beyond a few particle diameters and the velocity fluctuations remain proportional to container size. More recently, the idea that the container boundaries pl1 ,i a critical role in determining the amplitude of the velocity fluctuations has been explored in de tail (46, 66, 67). The primary goal has been to discover if container walls can quali tatively change the hydrodynamic interactions in the bulk suspension. Brenner (66) analyzed the effects of vertical container walls on the particle velocity fluctuations, but found that it does not eliminate the dependence on system size. However, numerical simulations (46) subsequently showed that rigid boundaries at the top and bottom of the vessel cause a strong timedependent damping of the velocity fluctuations even in bulk regions far from the walls. These boundaries introduce interfaces between sedi ment, suspension, and supernatent fluid, which are absent in homogenous systems with periodic boundary conditions. The timedependent damping of the velocity fluctuations was explained by convective draining of largescale fluctuations in particle density to these interfaces (46), a si. I i. o made earlier by Hinch (68). A different picture of the effect of container walls was proposed independently, and more or less simultaneously by Tee (67), where it was ir.. . 1. that hydrodynamic dispersion at the suspension supernatent interface leads to a stratification in the particle concentration. A stratified suspension introduces another length scale, beyond which the hydrodynamic interac tions are screened (69). However, a model based on convection of density fluctuations leads to qualitatively different conclusions from a stratified suspension, and we will compare these ideas to the results of our simulations in Secs. 4.3.1 and 5.1.1. Recent experimental results have cast doubt on the conclusion that the velocity fluctuations are necessarily independent of container size (67, 70). Brenner (70) found that the velocity fluctuations in very dilute suspensions (volume fraction Q < 1 .) did not saturate, even for container dimensions as large as 800a, where a is the particle radius. Tee et al. (67) found that for large cells and dilute suspensions they could not even obtain a steady state. Instead the velocity fluctuations d. .1iv for the duration of the experiment, no matter what the height of the container; these observations could be explained by an increasing stratification of the suspension with time. On the other hand some experiments found that the velocity fluctuations were steady and independent of system size, even in very large vessels (71). To make further progress, it will be necessary to reconcile the apparently conflicting experimental results, and to develop the correct physical picture underlying the screening of the hydrodynamic interactions in settling suspensions. The focus of the present work is on the microstructure of a settling suspension, as characterized by the distribution of pairs of particles in the bulk suspension. At low Reynolds number the instantaneous particle velocities are completely determined by the particle positions, and it can be shown that the dominant contribution to the velocity fluctuations can be calculated from the structure factor N S(k)= N 1 exp(ik rij), (4.1) ij 1 which is the Fourier transform of the pair correlation function (64, 65). The effects of hydrodynamic screening show up in the longwavelength (small k) behavior of the S(k); an ..imptotic k2 dependence of the structure factor indicates screening and an eventual saturation of the velocity fluctuations with increasing container size. It has not proven possible to measure S(k) directly by light I I ii:. due to the large size of the particles, but direct imaging has been used to calculate related spatial fluctuations in particle concentration (72). In our simulations we have been able to determine the structure factor directly. Here we present a detailed account of our results for both the velocity fluctuations and the microstructure, including effects of polydispersity, which may be important in interpreting experimental results. 4.1.1 Sedimentation Simulations Simulations of sedimentation are computationally demanding because of the large number of particles necessary to capture the length scales involved in the development of hydrodynamic screening. The simulation parameters are therefore a carefully chosen compromise between several competing factors, limited by the requirement that a calculation complete in a reasonable amount of time, of the order of one month. Experi ments (61) ,r.. I that the characteristic length scale is the mean interparticle spacing I = aQ1/3 and that the screening length is of the order of 101. The computational time required for our simulations is, to a first approximation, proportional to the number of fluid grid points, and therefore the total volume is an important limitation. We have limited our calculations to a single volume fraction 0.13 (1 I 2a), which was chosen as a reasonable compromise between keeping I small and avoiding the additional complications of a highly concentrated suspension. We used a cell with a square cross section, and studied a range of widths from W = 81 to W = 241. In most instances the cells were bounded at the top and bottom by rigid impermeable walls and we have found from experience that a cell height H = 1000a is necessary to allow time for the suspension to reach a steady state. We have varied the cell height in some instances, as reported in Sec. 4.2.2. At the chosen volume fraction the suspensions contains between 8000 and 72000 solid particles. A key component of this work is to assess the effects of the macroscopic boundary conditions on the suspension microstructure and dynamics. A noslip boundary breaks the macroscopic translational invariance in the direction normal to the boundary surface, and because of the longrange character of hydrodynamic interactions, this symmetry breaking can have a significant effect far from the boundary. We have therefore compared results for three different sets of boundary conditions: systems bounded in all directions by rigid walls, which we denote as "Box", systems bounded by walls at the top and bottom, while the vertical faces are periodic ("Bounded"), and systems with periodic boundary conditions on all faces ("PBC"). For systems with periodic boundaries on all faces, the height to width ratio was set to 4:1, as previous work showed no change occurred with taller cells (56). In the simulations with noslip boundaries at the top and bottom ("Box" and "Bounded"), data was taken in a window, typically of height 200a, centered around 400a from the bottom of the vessel, as is typical in experimental measurements. We have checked that the results are insensitive to the position of the viewing window as long as it is located in a region of uniform particle concentration. Fully periodic systems ("PBC") are spatially homogeneous, and data is averaged over the entire volume in this case. When considering the system dimensions, it should be noted that an extra 2a has been added to all dimensions bounded by noslip walls. This is to allow for the excluded volume with the walls and keep the particle concentration as close as possible to 1;:'. We will identify the system dimensions without any excluded volume correction. A particle radius of two grid spacings, a = 2Ax, was used in all these simulations, as compared with a = 1.5Ax or a = 1.25Ax in previous work (46, 56). Test calculations with small number of particles (see Sec. 4.2.3) i... i that there are no in i,i". changes in the results if larger particles are used. The largest systems studied in this work contains 72000 particles and approximately 20 million grid points. The simulations were run for 2000 Stokes times, where the Stokes time ts = a/Uo is the time it takes an isolated particle to settle one radius. The gravitational force was set so that the Stokes time corresponded to 250 time steps and the complete settling simulation therefore required half a million steps. The calculation takes about 1 month on a cluster of 8 processors. The typical particle Reynolds number, based on the mean settling velocity, was Re = 2p (U) a/p = 0.06, whereas in laboratory experiments it is usually one to two orders of magnitude smaller. Although the residual inertia is a potential source of complication, tests described in Sec. 4.2.3 and laboratory experiments (73) ~i, 1 effects of inertia are negligible at Reynolds numbers Re < 1. 4.1.2 Settling of Pairs of Particles Here we examine the settling of pairs of particles to assess the effects of grid artifacts and the small residual fluid inertia. Two spheres settling sidebyside repel each other when the spacing is of the order of the particle diameter (74, 75). The rotation of the particles gives rise to a lift force at finite Reynolds number, which tends to drive the particles apart. At very small Reynolds numbers, Re < 0.1, the repulsion between the particles become negligible and the separation distance does not change appreciably over time (76). We examine the settling of two particles at a range of particle Reynolds numbers 0.0025 < Re < 0.1. Fig. 41 shows that the particles maintain their orientation, and slowly drift apart with a velocity proportional to Re. * I I I ' 4  2 0 200 400 600 800 100 , 80 60 40 20 200 400 600 800 t/t, S Figure 41: Settling of a horizontal pair of particles at different Reynolds numbers; Re ~ 0.1 (circles), 0.05 (squares), 0.025 (triangles). Fig. 42 shows the separation distance Ar/a between particle centers as they settle with different initial orientations to the horizontal. Particles oriented in the horizontal (0 = 0) and vertical (0 = 90) directions maintain their initial orientations over time, although the particles drift away from each other (horizontal) or towards each other (vertical), due to the residual fluid inertia. At low Reynolds numbers, a pair of exactly spherical particles oriented at some intermediate angle should settle as a unit, maintaining the separation and orientation of the pair. Our results show that particles at offaxis orientations have slightly different velocities so the pair rotates with time, even at very low Re. This is a numerical artifact caused by the discrete lattice. It can only be reduced by using larger particles with more grid points to represent the particle surface, or by a more complicated boundary condition that smooths out the i ,.. d representation of the surface shown in Fig. 22. Nevertheless the grid artifacts cause much smaller variations in particle velocity than the polydispersity inherent in laboratory experiments. 4.2 Velocity Fluctuations In this section we examine the transient and steadystate fluctuations in particle velocity under different macroscopic boundary conditions. This issue has already been S. AAAAA 2 **** 0 200 400 600 800 100 I I 60A A A 40 U AA A A A A A 20 200 400 600 800 t/t S Figure 42: Settling of a pair particles at various orientations: 0 = 0 (circles), 450 (squares), 680 (triangles), 900 (diamonds). The Reynolds number Re = 0.1 in each case. examined from several different perspectives; here we present a more complete and uptodate account of our own investigations. 4.2.1 TimeDependent Velocity Fluctuations The initial particle configuration is sampled from the mostprobable (or maximum entropy) distribution of hard spheres, which was generated by an equilibrium Monte Carlo simulation. At moderate concentrations this distribution deviates from a random distribution because of the excluded volume. Thus the initial distribution does not obey Poisson statistics, but the pair probability is nevertheless random for separations of more than 510 particle radii. The particle velocities are initialized by a run of 200ts where the particle positions are frozen. Thus the initial configuration for the dynamical simulation (t = 0) is an idealized wellstirred distribution, without largescale spatial inhomogeneities. Simulations have shown that there is a qualitative difference in the time evolution of the velocity fluctuations depending on the macroscopic boundary conditions (46). The velocity fluctuations shown in Fig. 43 illustrate the different time dependence of "PBC" and "Box" geometries. With periodic boundary conditions the velocity fluctuations initially increase with time, reaching a plateau after approximately 200ts. I I I I I I 0.5 Periodic 0.5 Box ) 0.4  "A 0.3 0 0.2  0.1 E 0 0 500 1000 1500 2000 0.08  ^ 0.06 A t 0.04 E V 0.02 MENMONEE 0 0 500 1000 1500 2000 t/t Figure 43: Particle velocity fluctuations as a function of time with a with W = 48a. This has previously been explained in terms of an increasing number of particle contacts as the settling proceeds (56, 65). This is an example of a change in suspension microstructure, i.e. how pairs of particles are distributed, even though the particle concentration is exactly the same. However, when the container is bounded by no slip walls, the fluctuations decay with time, over a period of approximately 1OOOts. The initial fluctuations are comparable to the values in the homogeneous suspension ("PBC"), but in the "Box" geometry they decay over time to a much smaller value. A similar decay in the velocity fluctuations has been observed experimentally (70, 71), over comparable time scales. Simulations also demonstrate that the boundary conditions on the vertical faces of the container do not piw a qualitatively important role in determining the magnitude of the velocity fluctuations. The data illustrated in Fig. 44 show that velocity fluctua tions decay in a similar way for both "Box" and "Bounded" geometries, and therefore the amplitude of the velocity fluctuations is largely controlled by the boundary con ditions at the top and bottom of the container (46). As an additional test, we ran a few simulations with noslip vertical walls and periodic boundary conditions at the top and bottom; these results were similar to the fully periodic case. Vertical noslip boundaries may 1'1 i, an important role in determining the fluctuations in containers 0.4 I I A Box A" A 0.2 *OA A 0 500 1000 1500 2000 0.4 A 1 Ac Bounded EA A 1A 0.2 1 1 ** A*****k A 0 500 1000 1500 2000 tt ,S Figure 44: Particle velocity fluctuations as a function of time for 'Box" and "Bounded" geometries at various widths: W/a = 16 (circles), 32 (squares), and 48 (triangles) respectively. with large aspect ratios in the horizontal plane (66), but our results show that there is no qualitative effect in containers with a square crosssection. The decay of the velocity fluctuations in the bulk of the suspension cannot be explained by direct screening of the longrange interactions by the noslip boundaries. Although a noslip boundary does screen the hydrodynamic interactions of particles near the boundary, the screening does not extend to particles that are farther from the wall than they are to each other (66). Thus in the measurement window, which is typically ~ 5W from the bottom, the hydrodynamic interactions are not significantly affected by the container walls. Moreover, the timedependent decay of the fluctuations ~,. 1i that there is a rearrangement of the typical particle configurations during the settling process. The most likely mechanism is that the container bottom and the suspensionsupernatent interface act as sinks of fluctuation energy, as i.. 1. 1 earlier by Hinch (68). A horizontal density fluctuation is idealized as two regions (or blobs) side by side, one slightly heavier than the average and the other slightly lighter. These density fluctuations convect to one of these two interfaces and are absorbed by the density gradient at the interface. Scaling arguments ii. I that a horizontal density fluctuation of length I convects with a velocity ,. = Uo(01/a)1/2 with respect to the mean (68), so that fluctuations drain away on a time scale t,(1) l/t, a (1/aO)1/2S. (4.2) In the absence of noslip boundaries at the top and bottom of the container, largescale density fluctuations recirculate through the suspension, and the fluctuations in density and velocity are therefore time independent. To account for the velocity fluctuations reaching steady state, the model must include some mechanism for replenishing the largescale density fluctuations; otherwise the velocity fluctuations will continually decrease. We will assume that smallscale (of order a) density fluctuations are generated by conversion of gravitational potential energy, and then spread out by hydrodynamic diffusion resulting from shortrange multiparticle interactions. The characteristic length of the density fluctuation will grow to order I in a time of order td(l) I 12/D, (4.3) where D e Uoa is the hydrodynamic dispersion coefficient. By balancing the convection time tc with the diffusion time td we obtain a critical length scale, lc = af 1/3, beyond which fluctuations drain away more rapidly than they can be replenished by diffusion. Thus the system can reach a steady state with a correlation length that is independent of system width and proportional to the mean interparticle spacing a 1/3, as observed experimentally (61, 71). This is the physical idea underlying the convectiondiffusion model proposed by Levine et al.(64), which will be discussed in detail in Sec. 4.3.2. An alternative explanation, based on the development of a stratified density profile in the suspension (77), will be discussed in Sec. 4.3.1. Within a viewing window far from the sediment and supernatent interfaces the velocity fluctuations reach a quasisteady state after a time period of the order of 1OOOts. The duration of the steady state, prior to the supernatent interface reaching the viewing window, depends on the system height, as shown in Fig. 45. However, the mean settling velocity and the velocity fluctuations within the viewing window are otherwise unaffected by system height. For the height used in most of our simulations Sa 0.4 S 0.2 Ha = 1000 V m H/a = 1500 A H/a = 2000 0 1 0 500 1000 1500 I S 0.08 A A ^A  0.04 * 0 I 0 500 1000 1500 tt S Figure 45: Mean settling velocity, (U\\), and fluctuations in settling velocity, NAUX2), as a function of time for different system heights, H. (H = 1OOOa) we typically have an upper time limit to the simulations of about 1500ts, before the supernatent interface starts to interfere with the measurements in the viewing window. We find that the time taken to reach steady state depends on the container width (Fig. 44), ranging from 600 700ts for the smallest system (W 16a) to about 1000ts for the largest system (W = 48a). 4.2.2 SteadyState Settling A settling particle produces a disturbance to the fluid flow that decays as r1 at low Reynolds number, where r is the distance from the disturbance. This disturbance produces a fluctuation in the velocity of another particle that decays as R2, where R is the interparticle separation. Assuming the particles are distributed uniformly, then integration over volume leads to the wellknown result: (U2) x L, (4.4) where L is the container dimension (1). Precise results have been worked out in the lowconcentration limit for several different geometries (65, 77, 78). The data in Fig. 4 4 shows that the initial velocity fluctuations follow the Caflisch and Luke (1) scaling Cl 0 A Cl Cl 0 A Cl H U U S 6 I I I I  AA A AAAA A AA 0.06 A A A _A AA m o .**A A A  0 10 20 30 40 5( 0.04 0 A 0.03 AA AAA 0.02 SA A 0 AA 00 A A A 0.01 A li g' __ ei A 0 I I I L& A _AAAAA_ 10 20 30 40 x/a Figure 46: Particle velocity fluctuations at steady state. Profiles are shown for the "Box" system (solid symbols) at different widths: W/a = 16 (circles), 32 (squares), and 48 (triangles). (Eq. 4.4), while the steady state velocity fluctuations grow more slowly with system size. Table 41: Steady state particle velocity fluctuations in a monodisperse suspension as a function of system width. W/a Box Center Bounded Time (ts) Space (a) < AUIf > 7/U 16 0.033 0.036 0.027 Box 6001000 150550 32 0.048 0.052 0.053 8001200 200500 48 0.059 0.063 0.076 11001500 250450 < AUJ > U02 16 0.012 0.016 0.007 Bounded 7001100 200650 32 0.018 0.022 0.013 8001200 250600 48 0.022 0.028 0.018 10001400 250500 The steady state fluctuations across the viewing window are shown in Fig. 46. The time window for averaging was determined separately for each system, depending on the time to reach a quasisteady state (Fig. 44); the duration of the time window was 400ts in all cases. The results were also averaged over a range of vertical position, chosen so that the density and velocity fluctuations were constant apart from statistical errors (Fig. 54). The spatial and temporal windows used in Fig. 46 are listed in Table 41. For the "Bounded" geometry the velocity fluctuations were averaged 0 over the whole horizontal plane, while for the "Box" geometry the fluctuations were measured in a thin slice of width 2a located in the center of the container. Figure 46 shows that even though the velocity fluctuations grow less than linearly with width, they do not converge to a widthindependent value for the range of system sizes studied. 2 HH 0 Mono : Bounded Mono: Box 0 2 4 6 8 10 12 14 1 0.5 0 2 4 6 8 10 12 14 W ol/3/a Figure 47: Steady state particle velocity fluctuations as a function of system width for simulation versus experimental fit (solid line). The steady state velocity fluctuations are summarized in Table 41. Results for the "Box" geometry are averaged across the horizontal plane and along the centerline of the system (the central position in Fig. 46); results for the "Bounded" geometry are shown averaged over the horizontal plane. The results for the "Box" geometry show more tendency to saturation than the "Bounded" geometry, but it is not clear that this is a qualitative as opposed to a quantitative difference. It will take larger system sizes than are computationally feasible at present to resolve this issue. However, Fig. 47 shows that the simulation data agrees rather well with experimental measurements for comparable container sizes. The solid lines are fits (61) to a v iii. I j of experiments (61, 79, 80). The experimental data ,ir:. i that somewhat larger system sizes are needed for the velocity fluctuations to saturate and become independent of system size. We note that inclusion of polydispersity improves the agreement with experiment, but we 67 will see in Sec. 5.1.3 that a different screening mechanism is operating in polydisperse suspensions. It should be noted that in using the analytic fit to experimental data in Fig. 47, we have taken the smallest container dimension as the controlling length, whereas Segre (61) used the larger lateral dimension in fitting their experiments. They claimed a better fit from this choice, but it cannot be justified physically or theoretically. 4.2.3 Numerical Errors 0.5 I I Re =0.015 * Re =0.03 o 0.4 A Re =0.06 SRe =0.12 A .. Itlmam.. V 0.3 * 0.2 I I 0 500 1000 1500 o 0.08 . 0.04 a 0 I 0 500 1000 1500 t/t S Figure 48: Effect of inertia on the mean settling velocity and fluctuations in settling velocity. The simulations described in this paper are not expected to be a completely quantitative description of Stokesflow hydrodynamics. In particularly, we cannot disregard the possibility that inertial effects p1 i, a role on scales larger than the particle radius. At the volume fraction used in this work, we can expect that particle velocity correlations persist for distances of the order of 40a (61), and the Reynolds number based on this distance is 1.2. Nevertheless the data in Fig. 48 shows that, in a system of width W 16a (N = 8000), inertia p1l ,, no role for particle Reynolds numbers Re < 0.1. We have verified that this conclusion holds in a larger systems as well (W = 32a, N = 32000); no significant difference in velocity fluctuations was observed between Re = 0.06 and 0.03. Again computational limitations prevent a study * 0 4 0 V 0.2 a= 1.25 a 2.0 A a = 2.7 0 I 0 500 1000 1500 A 0.1 A < 0.05 0 I , 0 500 1000 1500 tt S Figure 49: Effect of grid resolution on the mean settling velocity and fluctuations in settling velocity. of the effects of inertia in larger systems, but experimental results (73) i.. 1 that it is small whenever Re < 1. An important concern with latticeBoltzmann simulations are artifacts introduced by the motion of particles across the grid. It was shown in Sec. 4.1.2 that these grid artifacts cause a weak but measurable dispersion in the trajectories of pairs of particle, and the same random noise may serve to weaken the microstructural changes that lead to suppression of velocity fluctuations. However, simulations with larger particles, which have smoother trajectories (81), are not noticeable different. The effects of grid resolution are shown in Fig. 49; the larger particles correspond to a more refined mesh. There are quantitative differences in the velocity fluctuations calculated with a = 1.25Ax (46), and the present results (a = 2Ax) over the time range 250 < t/ts < 1000. A further increase in particle size (a = 2.7Ax) does not lead to statistically significant changes. Finally, we have noticed that the mean settling velocity is not alvi independent of vertical position in the bulk region at Re = 0.06, although at lower Re the mean settling velocity was essentially constant. It turns out that this is an artifact of the finite compressibility of the latticeBoltzmann model, rather than an inertial effect. 0.4 Ma 0.004 Ma= 0.001 0.35 0.25 0."00 400 600 z/a Figure 410: Effect of compressibility on the spatial variation of the mean settling velocity. The data is for Reynolds number of Re = 0.06. Under the conditions of the simulations, with a kinematic viscosity v = O.167Ax2/At, a Reynolds number Re a 0.06 corresponds to a Mach number Ma = 0.004. By reducing the viscosity, we can reduce the Mach number, while keeping the Reynolds number constant. Figure 410 shows that at constant Reynolds number, Re ~ 0.06, the mean velocity is independent of position at sufficiently small Mach numbers. However, we could not detect significant differences in the fluctuations at either smaller Re or smaller Ma, although the nonuniform settling velocity does lead to a small reduction in particle concentration (~ 0.01) during the course of the simulation. 4.3 Microstructure The dynamics of suspensions at low Reynolds numbers are controlled by the distribution of particle positions, which is sufficient to determine the particle velocities at any instant of time. In a settling suspension, subtle shifts in the pair correlation function can have a dramatic effect on the macroscopic behavior. Specifically, it has been established that the amplitude of the velocity fluctuations are largely determined by the structure factor (64, 65), and in particular its lowk limiting behavior: xUr2 f S(k) F kokpk2] < U Ua 4 U6aJ k4 dk, (4.5) where S(k) was defined in Eq. 4.1 and X is a numerical factor of order 1. It has not been possible to measure the structure factor of a settling suspension of non Brownian particles experimentally, since the particles are too large for light scattering 70 measurements. Fluctuations in particle concentration have been measured within a cylindrical or rectangular window (72), but this only gives an angle average of the pair distribution. In numerical simulations it is possible to calculate S(k) as a function of both wavelength and direction. This is important since some theories predict that the structure factor becomes highly anisotropic at longwavelengths (64), with the horizontal fluctuations vanishing as k2 while the vertical fluctuations remain finite at all wavelengths. Here we report on the structure factor in a steadily settling suspension, and investigate the possibility of stratification of the particle concentration (77). 4.3.1 Particle Concentration S0.125 0.6 I I 0. 5 Bounded 0.4 t 0.120 S0.33 0.2 0.115 0.1 0 I 0.110 0 200 400 600 200 300 400 500 S0.125 S.6I I I I ' 0 .6 Box 0.5 0.4 0 .120 0.3 0.2 0.115 0.1 0 0 .110 0 200 400 600 200 300 400 500 z/a z/a Figure 411: Particle volume fraction as a function of height, z. The dots indicate the instantaneous average, and the steadystate (1000 < t/ts < 1400) density profile in the viewing window is shown in the .,.1i ,ent plots. Figure 411 shows that the density profile is uniform in the bulk, with a sharp interface between the suspension and supernatent fluid. A sharp interface has also been observed experimentally (82), but recently it has been ir..I 1 that hydrodynamic dispersion at the suspensionsupernatent interface could lead to a weak stratification of the particles and a nonuniform density in the bulk (3, 67, 77). The time evolution of the concentration profile, Q(z, t), can be described by a convectiondiffusion equation, which, in a Lagrangian frame moving with the mean settling velocity U(Q), can be written as + U'a = D92; (4.6) U' = QdU/dQ is the velocity of a density perturbation with respect to the mean settling speed and D is the hydrodynamic diffusion coefficient. The solution to Eq. 4.6 is qualitatively similar to the profiles shown in Fig. 411; in particular there is a bulk region where the concentration is constant. However, if the diffusion coefficient is very large, the suspensionsupernatent interface spreads into the bulk leading to a weakly stratified suspension, with a small density gradient even in the bulk (77). Stratification has been si:. 1 as a mechanism for suppressing velocity fluctu ations (69), by making it possible for fluctuations in particle concentration to reach a neutrally buovi,nt position in the suspension without draining to the interfaces. How ever, the density profiles shown in Fig. 411 are not stratified; instead the concentration profile is uniform in the bulk and the suspensionsupernatent interface is sharp. The plots of average concentration show that any mean variation in density is well below the statistical noise, even after averaging over several hundred Stokes times. Our data  ir:. I that any residual density gradient must have a characteristic length of at least 104a, or 10 times the container height. Hydrodynamic dispersion does cause a spreading of the interface, but this is com pensated by hindered settling, which convects the less dense regions at a higher velocity than the high density regions and thereby sets up a convective flux in opposition to the diffusive flux. Balancing the convective and diffusive concentration fluxes with respect to a frame moving with the mean settling velocity, we estimate an interface thickness, D\l/U' w 5a, at this volume fraction. Here the vertical dispersion coefficient, D = 4 K(UH) a, was taken from our simulations and is comparable to experimental measurements (80). Our simulations do not rule out the possible significance of in terfacial diffusion in very dilute suspensions (0 < 1 ), such as are commonly used in ParticleImageVelocimetry measurements (67, 70, 71), but we do exclude it as a general mechanism for hydrodynamic screening. 72 1.5 0. 5 [  1.5 0 1 2 S0.5 0 0 1 2 k a/l Figure 412: Structure factor describing horizontal density fluctuations, S(k), for various system sizes: W 16a, 32a, and 48a. 4.3.2 Structure Factor The structure factors were calculated directly from particle positions within the viewing window in simulations with the "Bounded" geometry. Periodic boundary conditions in the horizontal plane allow for more accurate Fourier analysis and this is reflected in the varying quality of data for fluctuations in vertical and horizontal directions shown in Fig. 413. Figure 412 shows that, in comparison with the initial equilibrium distribution, horizontal density fluctuations are strongly suppressed by no slip boundaries at the top and bottom of the container. On the other hand it is already known that suspensions with periodic boundary conditions do not undergo significant changes in microstructure during settling (56, 83). In particular the structure factor remains finite at all wavelengths. The significance of this result is that it demonstrates that the macroscopic boundaries have a profound effect on the distribution of particles in the bulk suspension. It may also be the mechanism by which hydrodynamic interac tions are screened during the settling process. A physical interpretation of the data in Fig. 412 is that, if the viewing window was divided into vertical slices of sufficient size, there would be essentially the same number of particles in each slice, rather than the expected variation of order VrN8, where N, is the number of particles in the slice. The structure factor in a settling suspension develops a strong anisotropy, as shown in Fig. 413. Although the horizontal density fluctuations are proportional to kI_, the vertical fluctuations tend to a nonzero constant at small kH. Damping of horizontal density fluctuations (at least as fast as k{) is the minimum requirement for hydrodynamics I, as was shown by Levine et al. (64). They proposed that there are two qualitatively distinct nonequilibrium phases for settling suspensions, an unscreened phase characterized by a random microstructure and a screened phase where the horizontal density fluctuations are damped out at long wavelengths. They derived an expression for the nonequilibrium structurefactor, S(k) k + .7) D i + Dlk2 + 7k/k2' (4.7) which is consistent in functional form with the structure factor obtained in our numer ical simulations (Fig 413). The renormalized parameters, the fluctuations in particle flux Ni and the diffusion coefficients Di, together with the damping coefficient, 7, were calculated from coupled field equations describing the evolution of the particle concen tration and fluid velocity. According to the theory, the phase boundary is determined by the anisotropy in the renormalized noise, NI/NH, and diffusivity, DI/D\\. The structure factor data can be used to extract ratios of the parameters that appear in Eq. 4.7; namely NL/7 = 0.4a2, and NII/D = 0.17. We obtained NL/7 and NII/DII from the lowk behavior of the horizontal and vertical density fluctuations, but since our data is rather noisy, it is impossible to extract a meaningful value of D/7, which appears as a quartic correction to the ..imptotic k2 dependence of the structure factor. Nevertheless, for the sake of completeness we will use our best estimate, Di/7 = 0.5a2, to determine the ratio NI/NH 0.7. When combined with tracerdiffusion measurements of D/DII = 0.16, this sI r' 1 we are near the transition between screened and unscreened phases (64). Unfortunately our data is not sufficiently precise to enable a definitive conclusion to be drawn. Significantly larger systems sizes ka/K Figure 413: Structure factor at different angles; vertical [0,0,1] direction (circles), 450 [1,0,1] direction (squares), and horizontal [1,0,0] direction (triangles). will be necessary for a quantitative comparison with the predictions of the theory, with greatly increased computational requirements. Density fluctuations have been measured at other lowindex directions; for example Fig. 413 also shows data for the 450 [1,0,1] direction. Unfortunately the system is not periodic for any kvectors that lie outside the horizontal plane, so this data is inherently noisier. Further complicating the analysis is the effect of renormalization, which produces a shoulder that is best seen for the horizontal [1,0,0] direction at ka 7r/4. For ka > 7/4 the structure factor is similar to the equilibrium distribution (Fig. 4 12) but the convective effects described qualitatively in Sec. 4.2.1 produce significant damping when ka < 7/4. The data for the 450 may also be showing a shoulder at a smaller k, ka 7 point, because hydrodynamic screening requires that the density fluctuations are damped in all directions except the vertical. In practice this means the shoulder is expected to move to smaller and smaller k as the direction rotates from horizontal to vertical. Although the structure factors measured in the simulations are consistent with the predictions of Levine (64), their theory postulates that all the important dynamics occurs in the bulk, independent of the macroscopic boundary conditions. Although this is a logical assumption, our numerical simulations have shown that it is incorrect. Simulations with periodic boundary conditions (42, 56) do not exhibit any damping of the horizontal density fluctuations, as would be expected if the model were correct in all essentials. Instead, our simulations t:: 1 that the container bottom and the suspensionsupernatent interface act as sinks of fluctuation energy, as . 1. 1 earlier by Hinch (68). Random density fluctuations convect to one of these two interfaces and are absorbed by the density gradient at the interface. The data shown in Fig. 414 supports this conclusion, albeit not conclusively. Here we show the structure factor in the viewing window during its evolution from an equilibrium state to the steady state. Despite the limited time averaging (a total of 4 x 200ts = 800ts for each plot), the implication is that the long wavelength fluctuations decay fastest. If so, this is evidence for the immediate convection of largescale density fluctuations (68), rather than the establishment of a density gradient by hydrodynamic diffusion (77). The screening by Mucha (77) is generated by transient vertical variations in particle concentration, rather than by a steadystate change in pair correlations. The key difference is that the theory by Levine (64) adds strongly anisotropic concentration fluctuations, which are an empirical representation of the supposed effects of the ii ,' body hydrodynamic interactions. To obtain a hydrodynamically screened phase, the theory requires that smallscale concentration fluctuations are largest in the horizontal plane, N% > N (64). By contrast, in a random suspension the density fluctuations are isotropic on all length scales. Our simulations show that there is a pronounced change 1.5 Equilibrium t/t = 0200 t/t = 200400 A t/t 400600 106 50 4 Figure 414: Time evolution of structure factor for horizontal density fluctuations. in the anisotropy of the density fluctuations as the settling proceeds, so the character of the noise may change as the suspension evolves to steady state. 4.3.3 Mean Settling Speed Our results i.. I that the mean settling velocity decreases during the settling process by about 25'. as shown in Fig. 48. This is a generic result whenever there is a noslip wall bounding the top and bottom of the container, and is independent of system height (Fig. 45), Reynolds number (Fig. 48) and grid resolution (Fig. 49). To our knowledge a systematic variation of settling velocity with time has not been reported previously, but it is consistent with a timedependent reduction in number density fluctuations, which has been observed experimentally by Lei et al.(72). In that work it was shown that the numberdensity fluctuations in a fixed volume decreases as the suspension settles. We have observed a similar decrease in numberdensity fluctuations in our simulations, but changes in microstructure show up more clearly in the structure factor (Sec. 4.3.2). If the change in settling velocity is due to changes in the microstructure at long wavelengths, then it can be estimated from the relation between the settling velocity and structure factor (65), which is quite precise at long wavelengths: < AU > 6 [S(k) Seq(k)] k4 dk, (4.8) The steadystate structure factor was taken from Eq. 4.7 with the parameters deter mined from the simulation data (Sec. 4.3.2), and the isotropic equilibrium structure factor was fitted with a quadratic polynomial; the integration was taken up to a maxi mum ka = 7/2, beyond which it was assumed that the structure factors were similar. This gave a predicted decrease in the mean settling speed < AU >= 0.20Uo, while the simulation result was around 0.13U0. 4.4 Conclusions The focus of this work has been to address the role of macroscopic boundary conditions on the microstructure of a settling suspension. In monodisperse suspensions we have found that there is a rearrangement of the pair distribution in the bulk region of the suspension, which suppresses the longwavelength density fluctuations, especially in the horizontal plane where a clear k2 dependence was observed. Simulations with different macroscopic boundary conditions on the container walls have demonstrated that the boundary conditions at the top and bottom of the container plh a crucial role in determining the distribution of particle pairs. The measured structure factor is consistent with the key qualitative features predicted theoretically by Levine et al. (64) using a renormalized convectiondiffusion model. However, we note that this theory cannot yet explain why the macroscopic boundary conditions should plhJ a crucial role. Our simulations predict that the mean settling velocity decreases during the settling process. While this has been shown to be consistent with the observed changes in suspension microstructure there is no experimental confirmation at the present time. We have found no evidence for stratification in monodisperse suspensions at moderate concentrations. Our results show that the interface is sharp and that the concentration in the bulk is very uniform. The , ..Ii. i by Mucha (77) that our observations (46) could be explained in terms of stratification is incorrect. We con clude that there is a mechanism for microstructural rearrangement in the bulk of a monodisperse suspension during settling, which leads to a substantial reduction in the amplitude of the velocity fluctuations from the predictions for randomly distributed particles. Based on our simulations and the experimental measurements of Bernard Michel (70), we believe that the question of saturation of the velocity fluctuations in monodisperse suspensions is still open. Both our simulations and these experi ments show the steadystate velocity fluctuations growing less rapidly than linearly in container size, W, but not yet fully saturating. CHAPTER 5 SEDIMENTATION OF A POLYDISPERSE SUSPENSION 5.1 Introduction Particles used in laboratory measurements have a typical polydispersity in the range 5 10'. (67, 71, 80) although some experiments use significantly more monodis perse particles (61, 70). In a fluidized bed it is well known that particles segregate according to size, with a denser suspension of larger particles at the bottom and less dense suspension of lighter particles at the top. The larger particles have an inherently higher settling velocity and therefore match the fluidization velocity at a higher con centration than the smaller particles. Thus a fluidized bed of polydisperse particles is inevitably stratified in both concentration and volume fraction. We have studied the settling of polydisperse suspensions to find out if there is significant segregation during the settling process, and what effects that may have on the velocity fluctuations and microstructure. 5.1.1 Stratification Figure 51 shows the instantaneous density profiles in a 10'(. polydisperse suspen sion at steady state (t = 1200ts). In contrast to a monodisperse suspension (Fig. 411), the suspension supernatent interface is considerably broader, as shown in Fig. 5la, and there is a weak underlying stratification of the particle volume fraction. The volume fraction profile is broken down in Fig. 5 lb into contributions from three different size ranges. It can be seen that there is considerable segregation by this time and the suspensionsupernatent interface is dominated by the smallest particles. Small particles tend to drift into the upper region of the front, while the large particles settle faster and are dominant in the region nearest the dense pack; medium size particles are dis tributed throughout the suspension. We emphasize that the interface spreading seen in Fig. 5la is a result of differential convection of different particle sizes rather than hydrodynamic diffusion (77), which we have seen does not produce a broad interface at 0.6 a) 0.3 b) 0.5 Large 0.4 0  Medium 0.0.2 Small 0.3  0.2 0. 0.1 , 200 400 600 0 200 400 600 z/a z/a Figure 51: Profiles of the particle volume fractions as a function of the system height. this volume fraction. Even at concentrations less than 1 convective spreading due to polydispersity is often more important than hydrodynamic diffusion (84). Our simulations, Fig. 52, show that there is a stratification of the particle mass density (or volume fraction) in polydisperse suspensions, which persists deep into the bulk. In comparison with the monodisperse suspension, here there is a clear decrease of the particle mass density (or volume fraction) with height that is visible over the statistical noise. Segregation and stratification can be explained by a onedimensional convectiondiffusion equation analogous to Eq. 4.6: Ot i + 0 (U ) = D2. (5.1) where the particle sizes have been divided into a number of different size ranges, with volume fraction Oi and settling velocity Ui Ui,o f(b). The settling velocity was constructed from the usual approximation that takes the settling velocity of the isolated particle and a hindrance function, f(y), based on the total volume fraction y= Yi. In the absence of diffusion, D = 0, we find an interface spreading and segregation that is qualitatively similar to the simulation data, confirming that polydispersity is the dominant mechanism for stratification at this volume fraction. More complicated explanations of stratification based on interface diffusion (77) or convection of density fluctuations (85) have not considered the effects of polydispersity, which may partly or even fully account for the experimental observations (see Sec. 5.1.4). 0.125 0.125 Mono Poly 0.120 0.120 4 0.115 0.115 0.110 0.110 0.105 _______ 0.105 ______,,, 200 300 400 500 200 300 400 500 z/a z/a Figure 52: Comparison of particle volume fraction profiles for monodisperse and poly disperse suspensions at steadystate, 1000 < t/ts < 1400. 5.1.2 Velocity fluctuations Our simulations indicate that the concentration profile at the interface between a polydisperse suspension and the supernatent fluid is continually evolving during settling, so the designation of a steadystate may be more arbitrary than in the monodisperse case. Nevertheless there are time windows of about 400ts duration when the velocity fluctuations are essentially stationary, as shown in Fig. 53. The time dependence and steadystate values of the velocity fluctuations are comparable to the monodisperse case (Fig. 47), although the anisotropy between vertical and horizontal velocity fluctuations is about 4 for the polydisperse suspension, comparable to experiment, while it is somewhat less for monodisperse suspensions. The velocity fluctuations in the different geometries are summarized in Table 51 in the same format as in Table 41. Table 51: Steadystate particle velocity fluctuations in a 10'. polydisperse suspension as a function of system width. W/a Box Center Bounded Time (ts) Space (a) < AU~f > /Uo 16 0.053 0.061 0.053 Box 6001000 150400 32 0.083 0.094 0.095 8001200 200400 48 0.104 0.115 0.144 11001500 250400 < AU_2 > U02 16 0.011 0.015 0.009 Bounded 8001200 200500 32 0.021 0.028 0.018 8001200 250500 48 0.029 0.037 0.029 10001400 200400 81 0.4  0.3 0 A 0.2  *A S0.12 0 I I I 0 500 1000 1500 2000 0.08 I 1 0.06 A 0.04  r A i A A V 0.02 *II *I I : : I I i A 0 00000000 A 0 500 1000 1500 2000 t/t Figure 5 3: Particle velocity fluctuation as a function of time for a 1C0. polydis perse suspension. The data is taken for three system widths: W/a 16 (circles), 32 (squares), 48 (triangles). The profiles of the velocity fluctuations as a function of height are qualitatively different in monodisperse and polydisperse suspensions. Figure 54 shows that the velocity fluctuations in monodisperse suspensions are independent of vertical position in the bulk suspension, while in polydisperse suspensions they decay with increasing height. It has been pointed out by Mucha et al. (77) that a stratified suspension is expected to lead to a reduction of the velocity fluctuations near the suspension supernatent interface, because the density gradient is larger there (Fig. 5la) and so is more effective in damping the velocity fluctuations. The decreasing fluctuations shown in Fig. 54 is further evidence that the polydisperse suspension is stratified, while the monodisperse suspension is not. The data in Figs. 52 and 54 ,., 1 that the velocity fluctuations are controlled by stratification for z > 400a. The gradient in volume fraction at this point (z = 400a), 3a = d(logQ)/dz 2.5 x 104, is consistent with a scaling argument (77) that si. I velocity fluctuations are controlled by stratification when 3 > critical, which for our system is roughly 2 x 104a1. 0.1 1 1 0.1 1 1 Mono Poly 0.08 0.08 0.06 0.06 A  0.04 0.04 0.02 .A:, : 0.02 0 200 400 600 800 0 200 400 600 800 z/a z/a Figure 54: Comparison of particle velocity fluctuations for monodisperse and 10', polydisperse suspensions. 5.1.3 Structure Factor Our simulations show that small amounts of polydispersity destroy the delicate mi crostructural rearrangements observed in strictly monodisperse suspensions. Horizontal density fluctuations in 1('. polydisperse suspensions are not damped at long wave lengths as they are for monodisperse suspensions, as shown in Fig. 55a. This ir.. I that different mechanisms for damping the velocity fluctuations exist. In monodisperse suspensions, the distribution of particle pairs adjusts itself during settling so that the density fluctuations are damped as k2 at long wavelengths. This leads to a time dependent damping of the velocity fluctuations and the possibility of a sizeindependent saturation. In polydisperse suspensions the microstructure is apparently randomized by the varying settling speeds and this mechanism no longer holds. However in this case the particle velocity fluctuations may be damped by stratification due to segregation of the different particle sizes. The structure factors for different degrees of polydispersity are compared in Fig. 5 5b. The systems are smaller in this case and so there are fewer kvectors. The structure factor for 2'". polydispersity is similar to the monodisperse case, apparently vanishing as k_ at long wavelengths. For higher degrees of polydispersity, 5'. and 10(. the structure factor tends to a nonzero value at low k1. It appears that only very small degrees of polydispersity can be tolerated if laboratory experiments are to mimic the properties of a monodisperse suspension. a) b) a 0.4 0.4 0 A A* A ** U 0.2 0.2 U0 0.2 0.4 0 0.2 0.4 k a/c ka/Ic Figure 55: Comparison of structure factors for different degrees of polydispersity: a) monodisperse (circles) and 10(. polydisperse (squares) suspensions (W 48a, N 72000); b) (circles), 5'. (squares), and 1C. (triangles) polydisperse suspensions (W 32a, N 32000). 5.1.4 Stratification in Laboratory Experiments Our investigations ir.. that there may be more than one mechanism at work in determining the velocity fluctuations in a settling suspension. In moderately dense suspensions, interface diffusion is prevented by the stronger convective effects of hindered settling, leading to a sharp interface and a uniform particle concentration in the bulk. On the other hand, moder ate to high concentration experiments (79, 80) have used suspensions with significant polydispersity, in excess of 5'. In this case the particle density becomes stratified over the duration of the experiment, Fig. 52, even deep into the bulk where measurements are made. The development of longrange correlations (Fig. 4 12) are hindered by even small amounts of polydispersity, so that stratification may be the dominant mechanism for controlling the velocity fluctuations in these experiments. We have solved the convectiondiffusion equation (Eq. 5.1) at 1t;: concentration and 1A. polydispersity; we find a stratification f3a 3 x 104 at a point H/3 of the way up the container and at a time when the suspension front has fallen approximately H/2. The calculated stratification is comparable to what was found in our polydisperse simulations (H lOOOa) and is sufficiently strong to control the velocity fluctuations (Fig. 5 4). Many experiments (61, 67, 70, 71, 84) are done at very low particle concentra tions, to permit the use of ParticleImageVelocimetry measurements. In this case hindered settling is negligible and the interface spreads by both diffusion and con vective segregation of different particle sizes. It is not computationally feasible to simulate suspensions at these low concentrations by the latticeBoltzmann method (Sec. 4.1.1), but it is interesting to consider whether stratification can be important in these systems. Mucha (77) has considered this question in detail for monodisperse suspensions; here we compare the contributions of segregation and dispersion, via solu tions of Eq. 5.1. Again we consider the stratification /3 at a location H/3 after a time H/2Uo as a reference point. Hindered settling leads to the development of a traveling concentration profile of constant shape (82), but in dilute suspensions the shape of the profile is constantly evolving. In the absence of hindered settling, the natural length scale is the height of the container and in these units there are only two parameters; the degree of polydispersity, a, and the dimensionless diffusion coefficient, D* = D/UoH; the dimensionless stratification is now 3H. In the absence of polydispersity (a = 0), large diffusion coefficients D* > 0.01 are needed to produce significant stratification. In our opinion this is a weakness of the analysis of Mucha (77); in order for interface diffusion to produce significant stratification, numerical values of the diffusion coeffi cient must exceed all experimental measurements. However, even modest amounts of polydispersity can produce significant stratification. For a = 0.1, stratification is in fact dominated by segregation; a reduced diffusion coefficient D* 0.01 is insufficient to contribute significantly to the overall stratification. Even at 5'. polydispersity, segre gation produces significant stratification and can be assisted by a small hydrodynamic dispersion, D* = 0.001. At 2'. polydispersity or less, segregation is insignificant, and the convectiondiffusion model .:: I0 that only interface diffusion can produce stratification. Of the experiments listed at the beginning of the previous paragraph, only one set (70), systematically examines the dependence on container size using nearly monodisperse particles, with polydispersity as low as 1 In these experiments a clear saturation with container size was not observed, particularly at the lowest volume fraction. It also took significantly longer for the suspension to reach a steadystate than in other experiments. The polydispersity was about 7'. in some experiments (67, 71, 84), which is sufficient for significant stratification by segregation. Segre et al. (61) used nearly monodisperse particles (o = 0.01), but primarily varied the volume fraction rather than the container size. We conclude that the issue of the saturation of the velocity fluctuations in a strictly monodisperse suspension is still an open question. 5.2 Conclusions Most experiments use particles with significant polydispersity in size. In this case our results show that the pair correlations are noticeably more random than in the monodisperse case and for a > 0.02, it seems likely that the driving force for microstructural rearrangements is no longer sufficient to combat the additional randomization from variations in particle velocity. However increasing polydispersity leads to stratification of the particle concentration even at moderate concentrations. Since in this case the stratification is driven by convective segregation of different particle sizes it can spread over large distances, even in the presence of hindered settling. We ~i. 1 that segregationinduced stratification may have a profound effect on the interpretation of experimental measurements, even in dense suspensions. Several recent experiments have been carried out at very low particle concentra tions, while our simulations are limited by computational efficiency to much larger volume fractions, in this case 0.13. Our results do not therefore discount the pos sibility that interface diffusion causes significant stratification at low volume fractions as claimed by Mucha (77). Nevertheless the convectiondiffusion model shows that segregation is often more important than interface diffusion in promoting stratification. Only in dilute suspensions of nearly monodisperse particles would we expect interface diffusion to be important, although stratification by segregation can occur under a broader range of circumstances. CHAPTER 6 SUMMARY We have applied the latticeBoltzmann method to examine particle dynamics in a gravity driven suspension. Initially, the particles exist in a randomly distributed state that exhibits large particle velocity fluctuations and no long range pair correlations within the microstructure. With periodic boundaries applied, the fluctuations are found to increase over short times and persist thereafter. However, noslip boundaries intro duced at the top and bottom of the vessel lead to a decay in the velocity fluctuation over time within the bulk, which is comparable with experimental findings. Here, the steadystate microstructure exhibits an ordered configuration at long wavelengths, which signifies that screening is present when noslip boundaries are applied. The numerical results are shown to compare well with experimental results, but due to computational limitations, a complete analysis of the sedimentation problem has not yet been attained. In experiments, the velocity fluctuation are found to become independent of system size beyond a certain point (61, 80). Our results do not unambiguously show this effect. The velocity fluctuations are found to grow less than linearly with system width, but they do not converge to a width independent value for the system sizes evaluated. However, quantitative results for the systems studied agree well with a fit to experimental data (61). Because experiments were carried out with larger containers before the screening effect became pronounced, simulations of larger system sizes should be evaluated for a more conclusive comparison. Luke (69) proposed a mechanism due to stratification by hydrodynamic diffusion at the suspensionsupernatent front leading to a damping of the velocity fluctuations within the bulk. We studied the concentration profile for a monodisperse suspension at 1;;''. volume fraction, but found a sharp front develop at the interface. Here, hindered settling, because of particle crowding interfering with the motion of the individual particles, has a stronger net effect than a diffusive spreading. The concentration was found then to be uniformly distributed within the bulk. This is contrary to 86 simulations (77) and experiments (61, 71) which were performed at low volume fractions, where stratification from diffusive spreading occurs more readily. As we see it, the stratification effect does suppress velocity fluctuations, but only at low volume fractions where hindered settling is not present. Without stratification, the mechanism for density fluctuation decay is then attributed to convective transport. We hypothesize that density fluctuations in the bulk can drain away at the suspensionsediment and suspensionsupernatent interfaces. This generates a continual decay in velocity fluctuations, which is replenished by hydrodynamic diffusion due to short range multiparticle interactions. A balance between convective and diffusive transport of density fluctuations leads to a correlation length scale beyond which the velocity fluctuations are screened. Further examination of the suspension microstructure showed that an evolution towards a statistically nonrandom state developed, where particles are distributed uniformly at large length scales. We identified this result based on the damping of the long range pair correlations in the structure factor, S(k) k as k1 goes to zero. Moreover, the pair correlations are found to be anisotropic between the vertical and horizontal planes, such that the vertical fluctuations were finite at all wavelengths. This correlates with a convectiondiffusion model (64), verifying that the transport of fluctuations due to convection is the dominant mechanism present leading to the decay in particle velocity fluctuations. Monodisperse suspensions represent an idealized situation where there is no varia tions in particle size. Typical experiments at moderate to high volume fractions (79, 80) exhibit size polydispersity of 5'. or greater. We find that for polydispersity in excess of 5'. a different mechanism for the decay in velocity fluctuations is observed. The variation in particle size leads to a segregation of different sizes as the suspension set tled. Smaller particles tend to settle slower and stay closer to the supernatent front, while large particles settle faster accumulating more into the bulk. The net effect was a pronounced stratification of particle concentration at the suspensionsupernatent front that persists into the bulk. Because of the stratification, the velocity fluctuations within the bulk do not have to convect towards the density gradients at the front, rather the 88 gradient is strong enough inside the bulk to create a sink for the velocity fluctuations to decay. The structure factor in the polydisperse suspension does not exhibit the same damping at long wavelengths seen in monodisperse systems, indicating a more ran domly distributed microstructure. However, velocity fluctuations are still damped, but this is due to stratification induced by segregation of different particle sizes. We i::. 1 then that polydispersity may have a larger effect in experimental measurements than previously supposed. 