Table Of Content

Steady state particle distribution of a dilute sedimenting suspension Bogdan Cichocki and Krzysztof Sadlej∗ Warsaw University, Institute of Theoretical Physics Hoz˙a 69, 00-681 Warsaw, Poland (Dated: 2nd February2008) Sedimentation of a non-Brownian suspension of hard particles is studied. It is shown that in 5 the low concentration limit a two-particle distribution function ensuring finite particle correlation 0 length can be found and explicitly calculated. The sedimentation coefficient is computed. Results 0 are compared with experiment. 2 n I. INTRODUCTION be achieved through a screening mechanism similar in a spirit to that introduced by Koch&Shaqfeh [6]. J 8 In recent years the familiar problem of sedimentation 2 of a non-Browniansuspension of particles has gained in- terest as new insight into the phenomenon is gathered. II. BBGKY HIERARCHY ] t Despite vivid research in the field, some questions still f o remainintriguingpuzzles. Evenwhen consideringobser- The system under consideration is a suspension of N s vationtimescalesanddimensionregimes,whichsanction spheres of equal radius a immersed in an incompressible . t the neglect of Brownian motion, the dynamics show to fluid of shear viscosity η. In the regime studied the par- a be many-body and highly chaotic leading to astonishing ticle Reynolds number is assumed to be small enough to m consequences. treat inertial effects as virtually absent. Moreover the - d One of the main, still open problems, is the determi- Pecletnumber describingthe relativeinfluence of hydro- n nation of the particle probability distribution function dynamic effects to Brownian motion, is presumed to be o for a steadily sedimenting suspension. It was first chal- large. On those assumptions and on the laboratory time c lenged by Batchelor [1]. In order to derive an equation scales imposed the fluid motion is governed by the sta- [ forthepair-distributionfunctionheconsideredapolidis- tionary Stokes equation. The dynamics of the particles 2 perse suspension. Taking then the limit of identical par- are therefore given by the equations of motion v ticles, a particle distribution function for mono-disperse 5 suspensions could be derived. But these results where dR N 9 i =U = µ (X)F (1) 4 shown to be ambiguous (depending on the way the limit dt i ij j 1 is taken) and therefore doubtful ([1], [2], [3]). The prob- Xj=1 0 lem of particle distribution is often skimmed overby the where µ is the translational part of the mobility ma- 5 assumption, that the particles are randomly distributed trix, principally a matrix function of the configuration 0 in the fluid. Indeed this will be the case, if Brownian / motionsplayasignificantroleinthedynamicsofthesys- X = (R1,...,RN) of all the particles (if not indicated t a tem, creating a randomising mechanism. But whenever otherwise, dimensionless - normalised to 2a - distances m will be used). It describes the hydrodynamic interac- this type of stochastic motion is absent, there seems to tionsbetweenparticlesandconnectstheforcesactingon - be no justification for such assumptions. Further, it has d particles to the velocities they acquire in a given config- been [4] pointed out, that a equilibrium (equal probabil- n uration. The forces F will be assumed to be constant o ity)distributionofparticlesresultsincharacteristicsofa j and equal F for all particles. c non-Browniansuspension, which have not been detected : inexperiments[5]. Moreoveritcanactuallybeexplicitly The particle distribution function satisfies the Liou- v ville equation, which for the system considered, has the shownthatsuchadistributioncannotbeasolutionofthe i X Liouville equation incorporating full multi-body interac- following structure r tions, therefore falling completely out of consideration. a ∂ρ(X;t) Inthisworkweproposeawaytotacklethisissue. The ∂t + ∇i· µij(X)Fρ(X;t) =0 (2) derivation of the statistical properties of the system will i,j X (cid:2) (cid:3) be based on the conjecture that correlations in steady state must be of finite length. This assumption, orig- Themobilitymatrixcanbefoundfollowingaformalpro- inating in the basic ideas of statistical mechanics, has cedure presented for example in ref.[7]. Summarising, profound consequences. It, as will be shown, leads to in order to calculate the hydrodynamic interaction, the an effective equation for the pair distribution function. boundary conditions on the particles surfaces are sub- Truncationof long rangehydrodynamic interactions will stituted by induced force densities introduced for each particle. These force densities, induced by fluid flow, produce subsequent fluid velocity fields which are prop- agated and reflect off other particles and thereby cause ∗Electronicaddress: [email protected] interactions. In consequence, the mobility matrix can 2 be expressedas a scattering sequence equivalent to a su- where in both equations, terms depending on higher- perposition of all these reflections. This sequence con- order distribution functions have been omitted. tains one-particle scattering operators, which construct Thereduceddistributionfunctionsfactorizeforgroups induced force densities on particles given velocity fields ofparticles thatarefar awayfromeachother. Therefore around them, and Green operators which propagate the these functions can be expanded in terms of correlation influence of a force density concentrated on a given par- functions according to ticle, resulting in fluid velocity fields aroundother particles. Inan unbounded, infinite system, the Greenopera- n(12) = h(1)h(2)+h(12) tor is the Oseen tensor. n(123) = h(1)h(2)h(3)+h(12)h(3) The full N-body mobility matrix can therefore be ex- +h(13)h(2)+h(23)h(1)+h(123) (8) panded in terms of the number of particles which enter eachtermofthescatteringsequence. Thissocalledclus- where h(1...s) is a s-particle correlationfunction which ter expansion has the structure vanisheswheneveranysubsetofparticlesisdraggedaway from the rest. If the correlations in a system are to be N µ (X)=µ(1)(1)+ µ(2)(1i)+... of finite length, which will be assumed, the correlation 11 11 11 functions must decayfasterthen the inverseofthe inter- i=2 X particle distance cubed. N µ (X)=µ(2)(12)+ µ(3)(12i)+... (3) Usingtheexpansion(8),theBBGKYhierarchyforthe 12 12 12 reduced distribution functions can be transcribed into a i=3 X hierarchy of equations for the time evolution of the suc- whereµ(s)(1...,s)denotesallthetermsofthescattering cessive correlation functions. An analysis of these equa- ij tionsshowsthatifsystemconfigurationissuchthatpar- sequence of µ which involve only, but all the particles ticlesformtwouncorrelatedclusters,someexpressionson {1,...,s}. the r.h.s,formulatedinterms of the scattering sequence, Next we introduce reduced distribution functions contain single Green operators (so called connection or N! articulationlines [7],[9]) joining these two groups of par- n(1,2,...,s)= d(s+1)···dNρ(X) (4) ticles. Such terms result in long range contributions be- (N −s)! Z cause a solitary Green operator includes slowly decaying These functions represent the average densities of pairs, parts proportional to 1/rγ where r is the relative dis- triplets, etc. of particles in given configurations. tance between the groups and γ =1,2,3. The evolution AnintegrationoftheLiouvilleequationinallthevari- ofacorrelationfunctionisthereforegivenbyanequation ables except the first s, together with the cluster expan- which contains long range terms and is therefore non- sion of the hydrodynamic interactions, leads to an infi- local in space. This structure of the hierarchy equations nite hierarchy of equations governing the time evolution is inconsistentwith the finite correlationlength assumed of the reduced distribution functions. This is the ana- for the correlation functions. As will be shown in the logueoftheBBGKY(Bogolubov-Born-Green-Kirkwood- next section, considering the low concentration limit a Yvon)hierarchyintroducedinthekinetictheoryofgases particle distribution can be found, which cancels out all [8]. It should be considered in the thermodynamic limit long range terms in the hierarchy, thereby saving finite when the number of particles and the volume of the sys- correlationlength. tem go to infinity, while the density is kept constant. Theequationsforthe twoandthree particledistribution functions have the structure III. LOW CONCENTRATION LIMIT PAIR DISTRIBUTION FUNCTION ∂n(12) = − ∇ ·µ (12)Fn(12) ∂t i ij Inthelowconcentrationlimitthedominatingpartofa i,j=1,2 X s-particlereduceddistributionfunctionisproportionalto − d3 ∇ ·[µF] (12;3)n(123)+...(5) ns. Consequently, the first terms of the hierarchy equa- i i Z i=1,2 tions have a well established order in density if the con- X centration of particles is low. where for example Considerthe equationgivingthe time evolutionof the two-particle correlationfunction. It is derived through a 3 substitutionof(8)into(5). Inthelowestorderofdensity [µF] (12;3)=[µ(2)(13)+µ(2)(13)+ µ(3)(123)]F (6) 1 11 13 1j it vanishes due to the relation [1] j=1 X ∇ ·µ (12)=0. (9) and i ij i,j=1,2 X 3 ∂n(123) =− ∇ ·µ (123)Fn(123)+... (7) which is a consequence of the isotropic nature of the ∂t i ij hydrodynamic interactions. The two-particle equation i,j=1 X 3 yields therefore in the dilute limit no condition for the Figure1: Thetwo-particlereduceddistributionfunctiong(R) two-particle correlation function. This peculiarity was satisfying thefinite correlation length criterion also encounteredby Bathelor[1]. His solutionwas based g(R) on the analysis of a polidisperse suspension where the two-particle equation yields a condition, but was shown 17.5 to be ambiguous ([1] ,[2],[3]). We turn our attention to 15 thenexthierarchyequation,i.e. thethirdequationwhich 12.5 emerges when (8) is addressed to (7), truncated to the terms proportional to the density cubed. 10 Aspointedout,ther.h.softhisequationcontainslong 7.5 rangeterms. If allthree particles are far awayfromeach other,bothsides ofthe equationcontainonlyshotrange 5 terms. But consider the case when one particle (e.g. the 2.5 particle with index 1) is far away from a close pair of R particles (consequently particles with indices 2 and 3). 1.02 1.04 1.06 1.08 1.1 All the terms, which in such a configurationlead to long range contributions can be identified and singled out resulting in the expression whereαtd andβtd arefunctionsappearingintheexplicit ij ij n∇ · (µ (2|1) + µ (23|1)+µ (2|3|1))n(23) + formofthe generalisedtwo-particlemobilitymatrixcon- 2 21 21 21 necting the translational velocity of the particles to the n∇3·(cid:2)(µ31(3|1) + µ31(32|1)+µ31(3|2|1))n(23)(cid:3) (10) stresslet [11]. Since limR→∞g(R)=1 where e(cid:2)ach vertical line stands for a single Gree(cid:3)n op- 1 ∞ 3(B−A) erator connecting uncorrelated particles. For example g(R)= exp dR (12) 1−A R(1−A) µ (23|1) stands for the sum of all terms of the scatter- (cid:20)ZR (cid:21) 21 ingsequenceofµ(3)(123)thatcontainonlyasingleGreen Notethatthissolutionisindependentofthedirectiondis- 21 tinguishedbytheexternalforcefield,theeffectsofwhich operator connecting particle 1 and the group (2,3). are in the low concentration limit completely dispersed In what follows, we will show that there exists a two- by the isotropic hydrodynamic interactions. Further, to particle probability distribution function which cancels verifytheconsistencyofthereasoningimposed,itcanbe out these long range terms and therefore leads to finite shownthatthetwo-particledistributionfunctionderived correlation length (12) guarantees convergence of all terms up to order n3 In equation (10) the long range contributions appear in equation (5) governing the time evolution of the two- in the Greenoperatorsbinding particle 1 with the group particle reduced distribution function. consisting of particles 2 and 3. Their explicit form is retrievedthrough an analysis of the multipole matrix el- ements of the Green operators, which can be found e. g. IV. RESULTS in [9]. When taking into account the symmetry relations A. Stationary particle distribution function upon exchange of particles, and a change of variables where r becomes the vector joining the distant particle and the center of the group of close particles, while R The functions A(R) and B(R) can be expressed as the vector denoting the relative position of the two par- series in inverse powers of inter-particle distance [11]. ticleswithinthegroup,wereachtheconclusionthatmost In numerical approximations of the integral (12) these of the long range contributions cancel out automatically sequences are truncated at 1000 and 800 terms respec- - only the multipoles proportional to the inverse of the tively. Lubrication corrections and far limit asymptotic inter-particle distance survive. with g(R)−1≈0.1953/R6 are taken into account. The A further expansion in terms of R = |R| (R ≪ r) computed function g(R) is plotted atfig. 1. Further, the shows that only terms proportionalto 1/r2 remain. The structurefactoratk=0iscalculated: S(0)=1−1.64φ. equationwhichemergeshasanidenticalstructure to the one derived by Batchelor&Green[10] and describing the relativemotionofa pairofparticlesinashearflow. The resultingparticledistributionfunctionn(R)=n2g(R)is B. Sedimentation coefficient isotropic. It is given by the solution of the differential equation In the low concentration limit the average sedimentation velocity U of a particle in the suspension measured 1 dg 3(A−B) dA (1−A) − − =0 (11) relative to the Stokes velocity U0 may be expanded in a gdR R dR series in the powers of volume fraction φ=4πa3n/3. (cid:26) (cid:27) A and B are functions of the scalar distance R given by U therelationsA=2(αt1d1+αt1d2)/RandB =2(β1td1+β1td2)/R U0 =1+φK +O(φ2) (13) 4 µ is the mobility of a single sphere. Dimensionless dis- 0 Figure 2: Dimensionless sedimentation velocity as a function tance R normalisedto 2a is adopted. The integral can of volume fraction - experimental results after Hanratty et 12 beinterpretedasacorrectionemergingduetothechange al.[14]. Pointsrepresentexperimentaldataandthesolidcurve of the distribution away from equilibrium. A numerical isafit. Thesolidlinerepresentsthesedimentationcoefficient K = −3.87, whereas the dotted line is the coefficient for a calculation yields the result K = −3.87. It is compared suspension in equilibrium with K =−6.55. with experimental data of Hanratty et al. on fig. 2. A fit to the experimental data of Ham et al. [15] leads to the sedimentation coefficient equal to −3.9, in excellent agreement with the calculated value. V. CONCLUDING REMARKS We have shown that a mono-disperse suspension of hard, non-Brownian particles can develop a micro- structure which insures finite correlation length. A scheme based on the formal schemes of non-equilibrium statistical mechanics and a detailed analysis of hydrodynamic interaction was used in the low-concentration regime to derive an equation for the two-particle distribution function. The solution shows to be isotropic due to the domination of the isotropic hydrodynamic interactions. The mechanism which insures the cutoff of long range correlations can be interpreted as hydrodynamic screening in the sense that the micro-structure of the suspension arranges itself in a way which truncates long range parts of the interactions. A close pair of particles feels the influence of a third distant particle through an The linear coefficient K can be expressed in terms of a effectiveshearingflow. Screeningleadstoaconfiguration microscopicexpression [13] involvingtwo-body mobility of the close pair, which hinders the effect of the distant matrices and the pair distribution function. particle. Theidea,thatasuspensionmightdevelopami- crostucture which changes the characteristics of it, was 2 2 introducedbyKoch&Shaqfeh[6]upontheanalysisofthe K =K + (g(12)−1)Tr µ(2)(12) dR 0 πµ 1i 12 problem of diverging velocity fluctuation. 0 ZR12≥1 "i=1 # X (14) Further, the sedimentation coefficient was calculated. whereK =−6.546([13],[12])isthesedimentationcoef- The value −3.87 was found. It agrees very well with the 0 ficient for a equilibrium distribution of hard spheres and experimental results cited. [1] G. K. Batchelor, J. Fluid Mech. 119,379, (1982) Mechanics, Wiley,New York,(1975) [2] G. K. Batchelor, C. S. Wen J. Fluid Mech. 124,495, [9] B. Cichocki, M. L.Ekiel-Jez˙ewska, P.Szymczak,E. Wa- (1982) jnryb, J. Chem. Phys. 117, 1231 (2002). [3] F. Feuillebois, Some Theoretical Results for the Motion [10] G. K. Batchelor, J. T. Green, J. Fluid Mech. 56, of Solid Spherical Particles in a Viscous Fluid, in Mul- 401 (1972) tiphase Science and Technology, edited by J. P. Hansen, [11] B.Cichocki,B.U.Felderhof,R.SchmitzPhys.Chem.Hy- D. Levesque, and J. Zinn-Justin (Elsevier, New York, drodyn. 10, 383 (1988). 1991), p763. [12] G. K.Batchelor, J. Fluid Mech. 52, 245 (1972) [4] R.E.Caflisch,J.H.C.Luke,Phys.Fluids28,259(1985). [13] B. Cichocki, B. U. Felderhof, J. Chem. Phys. 89, [5] P.N.Segré,E.Herbolzheimer,P.M.ChaikinPhys. Rev. 1049 (1988). Lett 79, 2574 (1997). [14] T.J. Hanratty, A. Bandukwala, AIChE J. 3, 293- [6] D. L. Koch, E. S. G. Shaqfeh, J. Fluid Mech. 224, 6 (1957). (275) (1991). [15] J. M. Ham and G. M. Homsy, Int. J. Multiphase Flow [7] B. Cichocki, P. Szymczak J. Chem. Phys. 121, 14, 533 (1988). 3329 (2004). [8] R. Balescu, Equilibrium and Nonequilibrium Statistical