<%BANNER%>

Sedimentation of hard sphere suspensions at low Reynolds number ssing a lattice-Boltzmann method

University of Florida Institutional Repository

PAGE 4

Iwouldliketothankmyadvisor,TonyLadd,forhisguidanceandunderstandingingivingmetheopportunitytoworkonthesedimentationproblem.Theworkwaslongandtedious,buthisprofessionalismandhardworkwereaninspirationthatkeptmegoing.IwouldalsoliketothankmyfellowcolleaguesintheChemicalEngineeringdepartment.Withsuchgoodpeoplearound,itmadethelonghoursofdoingresearchbearable.IthankmylabmembersJonghoonLee,foralwaysbeinghelpful;ByoungjinChun,forhiscongeniality;BerkUsta,forhisgenerosity;andMuraliRangarayan,forbeingagoodguy.ThanksalsogotoPiotrSzymczakandRolfVerberg,twopostdocswhoseinniteknowledgeandwisdomIwasfortunatetolearnfrom.Finally,IthankmyMom,Dad,brothers,andsisterforbeingapartofmylife. iv

PAGE 5

Page ACKNOWLEDGMENTS ................................ iv LISTOFTABLES .................................... vii LISTOFFIGURES ................................... viii ABSTRACT ....................................... xi CHAPTER 1INTRODUCTION ................................. 1 1.1PropertiesofaSedimentingSuspension .................. 2 1.2FluidDynamicEquations ......................... 3 1.3HydrodynamicsatLowReynoldsNumber ................ 5 1.3.1Long-RangeInteractionsandtheDivergenceProblem ...... 8 1.3.2Multi-ParticleInteractions ..................... 10 1.3.3ParticlesNearContact ....................... 11 2DESCRIPTIONOFTHELATTICE-BOLTZMANNMODEL ......... 13 2.1MolecularTheoryforFluidFlow ..................... 13 2.1.1BoltzmannEquation ........................ 14 2.1.2ConnectionbetweenMicroscopicandMacroscopicDynamics .. 15 2.2Lattice-BoltzmannMethod ......................... 17 2.2.1GasKineticstoFluidDynamics:TheChapman-EnskogEx-pansion ................................ 22 2.2.2Solid-FluidBoundaryConditions .................. 30 2.2.3HydrodynamicRadiusofaParticle ................ 33 2.2.4ParticleMotion ........................... 37 3INCORPORATINGLUBRICATIONINTOTHELATTICE-BOLTZMANNMETHOD ....................... 40 3.1Introduction ................................. 40 3.1.1SurfacesNearContact ........................ 42 3.1.2LubricationForces .......................... 43 3.1.3ParticleWallLubrication ...................... 44 3.1.4ClusterImplicitMethod ...................... 48 3.2Conclusions ................................. 52 4SEDIMENTATIONOFHARD-SPHERESUSPENSIONSATLOWREYNOLDSNUMBER ......................... 53 4.1Introduction ................................. 53 4.1.1SedimentationSimulations ..................... 56 v

PAGE 6

.................... 57 4.2VelocityFluctuations ............................ 59 4.2.1Time-DependentVelocityFluctuations .............. 59 4.2.2Steady-StateSettling ........................ 63 4.2.3NumericalErrors .......................... 66 4.3Microstructure ................................ 68 4.3.1ParticleConcentration ....................... 69 4.3.2StructureFactor ........................... 71 4.3.3MeanSettlingSpeed ........................ 75 4.4Conclusions ................................. 76 5SEDIMENTATIONOFAPOLYDISPERSESUSPENSION .......... 77 5.1Introduction ................................. 77 5.1.1Stratication ............................. 77 5.1.2Velocityuctuations ......................... 79 5.1.3StructureFactor ........................... 81 5.1.4StraticationinLaboratoryExperiments ............. 82 5.2Conclusions ................................. 84 6SUMMARY ..................................... 85 REFERENCES ...................................... 88 BIOGRAPHICALSKETCH ............................... 94 vi

PAGE 7

Page 2{1VarianceinthecomputeddragforceF=p .................................... 35 2{2Hydrodynamicradiusahy(inunitsofx)forvariousuidviscosities;aistheinputparticleradius. ............................ 35 3{1Lubricationrangesforvariouskinematicviscosities,determinedforasphereofradiusa=8:2x.Therangesaredeterminedseparatelyforthenor-mal,hN,tangential,hT,androtational,hR,motions. ............. 46 4{1Steadystateparticlevelocityuctuationsinamonodispersesuspensionasafunctionofsystemwidth. .......................... 64 5{1Steady-stateparticlevelocityuctuationsina10%polydispersesuspen-sionasafunctionofsystemwidth. ...................... 79 vii

PAGE 8

Page 2{1The18velocityvectors(ci)inthelattice-Boltzmannmodel. ........ 17 2{2Locationofboundarynodesforacurvedsurface(a)andtwosurfacesinnearcontact(b). ................................. 29 2{3Movingboundarynodeupdate.Thesolidcirclesareuidnodes,squaresareboundarynodeswiththeassociatedvelocityvector,andopencirclesareinteriornodes. ................................ 31 2{4DragforceFasafunctionoftime,normalizedbythedragforceonanisolatedsphere,F0=6aU.TimeismeasuredinunitsoftheStokestime,tS=a=U. ................................. 34 2{5Actual(solidlines)andhydrodynamic(dashedlines)surfacesforaparti-clesettlingontoawall. ............................ 36 3{1Normalforceonaparticleofinputradiusasettlingontoahorizontalpla-narsurfaceatzeroReynoldsnumber,withF0=6ahyU.Solidsymbolsindicate:=1=6(circles),=1=100(triangles),and=1=1200(squares).Resultsincludingthelubricationcorrectionareshownbytheopensymbols. .................................. 45 3{2TangentialforceonaparticlesettlingnexttoaverticalplanarsurfaceatzeroReynoldsnumber.SimulationdatafordragforcedividedbyF0=6ahyU,isindicatedbysolidsymbols.Resultsincludingthelubricationcorrectionareshownbytheopensymbols. ........................ 46 3{3TorqueonaparticlesettlingnexttoaverticalplanarsurfaceatzeroReynoldsnumber.SimulationdatafortorquedividedbyT0=8a2hyU,isindicatedbysolidsymbols.Resultsincludingthelubricationcorrectionareshownbytheopensymbols. .............................. 47 3{4TorqueonaparticlerotatingnexttoaverticalplanarsurfaceatzeroReynoldsnumber.SimulationdatafortorquedividedbyT0=8a3hy,isindicatedbysolidsymbols.Resultsincludingthelubricationcorrectionareshownbytheopensymbols. .............................. 48 3{5Settlingofasphere(a=4:8)ontoahorizontalwall.Thegapbetweentheparticlesurfaceandwall,h,relativetothehydrodynamicradius,isplot-tedasafunctionofthenon-dimensionaltime(opencircles). ........ 49 3{6Illustrationofthealgorithmtodeterminethelistofclusters. ........ 50 3{7Themaximumclustersizeasafunctionoftheclustercutogaphs=aatvaryingvolumefractions,. .......................... 51 viii

PAGE 9

............................... 51 4{1SettlingofahorizontalpairofparticlesatdierentReynoldsnumbers;Re0:1(circles),0.05(squares),0.025(triangles). ............. 58 4{2Settlingofapairparticlesatvariousorientations:=0(circles),45(squares),68(triangles),90(diamonds).TheReynoldsnumberRe=0:1ineachcase. ................................. 59 4{3ParticlevelocityuctuationsasafunctionoftimewithawithW=48a. 60 4{4Particlevelocityuctuationsasafunctionoftimefor`Box"and\Bounded"geometriesatvariouswidths:W=a=16(circles),32(squares),and48(triangles)respectively. ............................. 61 4{5Meansettlingvelocity,Uk,anductuationsinsettlingvelocity,DU2kE,asafunctionoftimefordierentsystemheights,H. ............. 63 4{6Particlevelocityuctuationsatsteadystate.Prolesareshownforthe\Box"system(solidsymbols)atdierentwidths:W=a=16(circles),32(squares),and48(triangles). ......................... 64 4{7Steadystateparticlevelocityuctuationsasafunctionofsystemwidthforsimulationversusexperimentalt(solidline). .............. 65 4{8Eectofinertiaonthemeansettlingvelocityanductuationsinsettlingvelocity. ..................................... 66 4{9Eectofgridresolutiononthemeansettlingvelocityanductuationsinsettlingvelocity. ................................. 67 4{10Eectofcompressibilityonthespatialvariationofthemeansettlingve-locity.ThedataisforReynoldsnumberofRe=0:06. ............ 68 4{11Particlevolumefractionasafunctionofheight,z.Thedotsindicatetheinstantaneousaverage,andthesteady-state(1000
PAGE 10

........................... 81 5{5Comparisonofstructurefactorsfordierentdegreesofpolydispersity:a)monodisperse(circles)and10%polydisperse(squares)suspensions(W=48a,N=72000);b)2%(circles),5%(squares),and10%(triangles)polydis-persesuspensions(W=32a,N=32000). .................... 82 x

PAGE 11

Particlemotioninasettlingsuspensionisdirectlyinuencedbyhydrodynamicinteractions,suchthattheparticlevelocitieswillexhibitarandomcomponent,char-acterizedbyahydrodynamicallyinduceddiusion.Statistically,thismeansthattheparticlestendtobeuniformlydistributedwithinthebulk.Theoreticalpredictionsforthesedimentationvelocityinhomogeneoussuspensionsareingoodagreementwithex-periments.Ontheotherhand,acalculationoftheparticlevelocityuctuations,basedonthesameassumptions,isinconictwithexperimentalresults.Ourmainobjectivewastounderstandthisdiscrepancy. Numericalsimulationshaveshownthatdynamicswithinthebulkarenotindepen-dentofthesizeofthecontainer,evenwhentheboundariesarelocatedfarawayfromtheobservationwindow.Wehypothesizedthatdensityuctuationsinthebulkdrainawaytothesuspension-sedimentandsuspension-supernatentinterfaces.Thisgeneratesacontinualdecayinvelocityuctuations,whichisreplenishedbyhydrodynamicdiu-sionduetoshortrangemultiparticleinteractions.Abalancebetweenconvectiveanddiusivetransportofdensityuctuationsleadstoacorrelationlengthbeyondwhichthevelocityuctuationsarescreened.Numericalresultswerefoundtosupportthishypothesis. Anexaminationofthemicrostructureasthesuspensionsettledshowedthatitevolvestowardastatisticallynonrandomstate,sothatparticlesaredistributed xi

PAGE 12

Ournumericalsimulationsshowthatthelong-rangecorrelationsfoundinmonodis-persesuspensionsaredestroyedbysmallamountsofpolydispersity,whicharefoundintypicallaboratoryexperiments.Thevariationinparticlesizealsoleadstoasegregationofdierentsizesasthesuspensionsettles.Theneteectisapronouncedstraticationofparticleconcentrationatthesupernatentfrontthatpersistsintothebulk.Becauseofthestratication,thevelocityuctuationswithinthebulkdonothavetoconvecttowardthedensitygradientsatthefront;rather,thegradientisstrongenoughinsidethebulktocreateasinkforthevelocityuctuations. xii

PAGE 13

Westudiedgravitydrivennon-Brownianparticleswithinaviscousuid.Wedealwithlow-Reynoldsnumberow,suchthatthedynamicsofparticlesinthedispersionisstronglyaectedbyhydrodynamicinteractionsgeneratedbytherelativemotionofparticlesanduid.Owingtothecomplexityofsuspensiondynamics,in1985,aparadoxforsteadystatesedimentationwaspointedoutbyCaischandLuke(CL)( 1 ).Considerasteadilysedimentingsuspensionofhardspheres.Aconcentrationuctuationgivesrisetoavelocityelddecayingas1=rwithdistancerfromtheoriginoftheuctuations.ThelinearityofStokesowimpliesthatthevelocityeldfrommanyspatiallydistributedconcentrationuctuationsissimplyasumiuioftheindividualcontributions.Iftheseconcentrationuctuationstakeplaceinarandom,spatiallyuncorrelatedmannerthroughoutthesuspension,theresultingvarianceu2inthevelocityatanypointinthesuspensionwouldbethesumofthesquaresoftheindividualcontributions.Thissumiu2ihasNL3termswithNsoluteparticlesinacontaineroflineardimensionLinalldirections.The1=rcontributionmentionedaboveleadstoavelocityuctuationofL2,sothatiu2L. TheCLpredictionisbasedontwosimplephysicalarguments: However,mostexperimentsdonotndthesizedependencepredictedbyCL;thereasonforthedisagreementisunclear.Onehypothesisisthatasucientlystronganticorrelationofdensityuctuationsdevelopsatlargelengthscales,whichsuppresstheCLdivergence.Wepresentndingsfromlattice-Boltzmannsimulationstoverifythishypothesis. Simulationsofthesedimentationproblemrequirelarge-scalenumericalcalcula-tions,forwhichthelattice-Boltzmannmethodisausefultool.Thesimulationisbased 1

PAGE 14

onthedevelopmentofmicroscopicinteractionsthatleadtomacroscopicuidow.ThisisdonebyapplyingtheBoltzmannequationdevelopedfromgaskinetictheory.Thesimulationusesaspatialgridfortheuidphase,makingitaninherentlyparallizeablealgorithmwheresimulationsoflarge-scalesystemsspreadovermultipleprocessorscanbeperformed.Withboundaryconditionsatthesurfaceoftheparticles,thelattice-Boltzmannmethodbecomesapplicabletomulti-phaseows.Weareablefollowthetrajectoriesofeachparticleinthesuspensionovertime,givingadetailedpictureoftheparticledynamics. FurtherobservationofasettlingsuspensionshowsthatparticlesexhibitrandommotionevenwhentheparticlesarelargeandBrownianmotionisnegligible.Theowdevelopedbyeachparticleinuencingtheothermakesthedynamicschaotic,andhighlysensitivetoinitialconditions.Theresultingchaosimpliesthatthetime-evolutionofcoarse-grainedquantitiesmustbedescribedusingdiusioncoecientsandnoisesources,eventhoughthemicroscopicdynamicsintheabsenceofthermalBrown-ianmotionareentirelydeterministic( 2 3 ).Thisinducedlarge-scalediusivebehavior,intheabsenceofthermalnoise,iscalledhydrodynamicdiusionorhydrodynamicdispersion. Importantmacroscopicquestionsariseinthestudyofthis"simple"prototypesystem.Forinstance:whatisthetimeofseparation?Whatisthedistributionofvelocitiesandconcentration?Howisthisprocessaectedbythecontainerboundariesanddimensions?

PAGE 15

4 ).ForauidvolumeV,themassdensityandaccelerationcanvaryfromplacetoplace.ThevolumeVisthereforedividedintoinnitesimallysmallvolumesdVoverwhichandaareconstant.HenceforanarbitraryvolumeV TheforceFexertedonauidvolumeVconsistsoftwodistinctcontributions: Hence Commonbodyforcesaregravityandexternalelectricormagneticelds.Deningfextasthebodyforceperunitvolume,wecanwrite ThesurfaceforcescanbefoundbyconsideringthestressesexertedontheboundaryAofthevolumeV: whereisthestresstensorthatspeciesthelocalnormalandtangentialforcesperunitareaexertedonauidelement;anddAisavectorpointingperpendiculartothesurface,itsmagnitudedeterminedbythesurfaceareaoftheinnitesimalelementdA.Usingthedivergencetheorem,Eq. 1.4 canbewrittenas FromEqs. 1.1 1.3 ,and 1.5 into 1.2 ,theequationofmotionforauidcanbewrittenas

PAGE 16

Inspecifyingtheaccelerationterm,anEuleriandescriptionofthevelocityu(r;t)ispresented.Hence, Thetwotermsontheright-handsidecanbeinterpretedasthelocalrateofchangeinvelocity(@u=@t)duetochangesatpositionrandaconvectiverateofchangeinvelocity(uru)duetothetransportoftheelementtoadierentposition( 4 ).Thevalueruisasecond-ordertensordescribingthegradientofthevelocityeld.Itcanbesplitintosymmetricandantisymmetricpartsas 2ru+ruT+1 2ruruT whereruTisthetransposeofru.Equation 1.8 canthenberewrittenas (1.9) whereSistherate-of-straintensorthatdescribesuiddeformation,andisthevorticitytensorthatcorrespondstoasolid-bodyrotation( 4 ). Thestresstensormuststillbedened.Intheabsenceofuidow,thestressesinauidarecausedbythestaticuidpressure,whichisisotropic(i:e:, wherePisthehydrostaticpressureandIistheunittensor).However,whenuidowoccurs,thestresstensorbecomesanisotropic,andcanbeexpressedasthesumofanisotropicpressureP(usuallydierentfromthestaticpressure)andanonisotropiccontribution.ForincompressibleNewtonianliquids,thisnonisotropiccontributionisgivenby2S,whereistheviscosityoftheuid.Thisproportionalitybetweenthestressandtherateofstrainisaphenomenologicalrelationship,whichmanyuidshavebeenfoundtoobey( 4 ).Hence

PAGE 17

Takingthedivergence(andusingruT=0foranincompressibleuid)gives SubstitutingEq. 1.7 and 1.12 intoEq. 1.6 givestheequationofmotionforthevelocityeld whichtogetherwiththecontinuityequationforanincompressibleuid arereferredtoastheNavier-Stokesequations. TheNavier-Stokesequationscanbewritteninnondimensionalformbyintroducingacharacteristiclengthscalelandvelocityu0.Deningthedimensionlessvariables ~t=u0t=l;~r=r=l;~u=u=u0;~P=lP=u0and~fext=fextl2=u0 thensubstitutingthemintoEq. 1.13 gives @~u wherethetildesrepresentnondimensionalparametersandoperators.TheconstantdevelopedontheLHSisdenedastheReynoldsnumber, Re=lu0 ; andcanbeinterpretedasameasureoftheratioofinertialtoviscousforces.

PAGE 18

owelduaroundthesphere,andcalculatethestressexertedonthespheresurfacebytheuid.WhentheRe1,thetermsontheLHSofEq. 1.16 areneglected,andtheNavier-StokesequationreducestotheStokesorcreeping-owequation.Indimensionalform,Eq. 1.13 canbere-writtenas withru=0: AsolutiontotheStokesequation,Eq. 1.18 ,canbedevelopedusingFouriertransforms( 5 6 ).First,theFourier-transformpairisdenedinrealspacerandFourierspacekas ^u(k)=Zeikru(r)dr (2)3Zeikr^u(k)dk: TakingtheFouriertransformofEq. 1.18 yields Takingthedotproductofkonthetransformedequationofmotion,Eq. 1.22 ,elimi-natesthe^utermtogive ^P(k)=iFk Substituting^PintoEq. 1.22 gives ^u(k)=1 Toinvertthetransformedsolutions(Eqs. 1.24 and 1.25 )weusethefollowingrelation 4r;

PAGE 19

1 (2)3ZF 4r and1 (2)3ZFkk 8rrr FromEqs. 1.24 and 1.26 ,thepressureeldisdeterminedtobe P(r)=1 4r3rfext; withfextrepresentingtheforcefromapointsource.ThevelocityeldisdeterminedbytakingtheinverseFouriertransformofEq. 1.25 withEqs. 1.27 and 1.28 togive 8I Hence,theOseentensorisgivenas 8I Wenowlookatthecasewhenthefullvelocityeldisnotneeded,butonlytheforceonaparticleisdesired.ConsiderarigidsphereofradiusaandsurfaceSimmersedinauidwithvelocityeldu1(r)intheabsenceofthesphere.Thesphereiscenteredatr0,translateswithvelocityU0androtatesaboutitscenterwithangularvelocity.Owingtoitsmotion,thesphereaectstheuidthroughforcesdistributedoveritssurface,(r)perunitarea.FromEq. 1.30 thevelocityeldatrisnowu1(r)plusthatgeneratedbythesuperpositionofforces(r0)dr0( 5 ),or Atthesurfaceofthesphere,thismustmatchthespherevelocitygivenbyU=U0+(rr0).Thus, IntegrationofEq. 1.33 overthespheresurfaceresultsin

PAGE 20

Expandingu1(r)aboutthespherecenterforcreepingow,andintegratingoverSyieldsFaxen'srstlaw 6a2r2u1)0U0: Thesubscriptzeroindicatesevaluationatr0,theparticlecenter.Thisisanimportantresult,showingthatinanonuniformow,theparticledoesnotmovewiththemeanuidvelocity(u1)intheabsenceofanexternalforce,butthatthereisacorrectionproportionaltoa2r2u1. Faxen'ssecondlawcanbederivedbytakingthevectorproductofEq. 1.33 with(rr0)beforecarryingoutthesecondintegration,andthenseparatingtheequationintoantisymmetricandsymmetricparts,whichmustbesatisedindividually( 5 ).Inthiswaythetorqueonasphere( 7 )isobtainedas 2(ru1)0 andtheforcedipole(stresslet)as 3a31+1 10a2r2(ru1+ruT1): FromEq. 1.35 ,Stokes'slawfortheforceexertedonaspheresettlinginainnitemediumisobtained.BysettingU0=0oru1=0weget andthetorqueis 1.32 mustbeevaluated.Bysetting=F=4a2=3U0=2aandexpandingtheOseentensor

PAGE 21

aboutthecenterofthesphere( 5 )leaves 2(r0r0)(r0r0):rrO(R)+:::dr0=6aU01+1 6a2r2O(R) (1.40) withR=rr0,andU0representingtheStokessettlingvelocityforanisolatedsphere.Equation 1.40 isderivedfromtheidentities andtheremainingtermsintegratetozero,since 1.31 )intoEq. 1.40 givesthevelocityeld( 5 ) whichsatisestheboundaryconditions Therefore,thedisturbanceislongrange,decayingas1=R. Nowweaddonemoresphereatapositionr2fromr1sothatjr1r2j=R12a.TorstorderfromEq. 1.44 ,Sphere1feelstheowcausedbySphere2withthefollowingstrength wherethevector^R12istheunitvectorR12=R12.ThevelocityofSphere1isnowincreasedwithrespecttothesinglespheresedimentationvelocity

PAGE 22

IfthesamereasoningisappliedtothesedimentationofacloudofNidenticalspheres,thentorstorderintheinverseoftheinterparticledistance,weget Becausethedisturbanceintheuidduetoparticleinteractionsislongrange,thesummationinEq. 1.49 willdivergeasthevolumeofthesystembecomeslarge. ThedivergenceproblemwassolvedbyBatchelor( 8 ),whopointedoutthatthereisanupcurrentofdisplaceduidthathadnotbeenconsidered.Thisbackowcancelsthelong-rangehydrodynamicinteractions,andcausesthesettlingvelocityofasuspensiontobelessthanthatofanisolatedparticle. 1.3 ,thedevelopmentoftheresistanceandmobilityrelationsforasingleparticlemovinginauidcanbeobtained,where TheforcesFandtorquesTexertedontheparticlebytheuidarerelatedtotheparticlevelocitiesUandangularvelocitiesviathecongurationdependentfriction-coecients( 6 9 ),whicharedevelopedfromthesolutionoftheStokesequations. Formulti-particleinteractions,theaboveequationsbecome wherethesummationextendsoverallNparticlesinthesuspension,andistheresistancematrix.ThesuperscriptsTandRrefertothetranslationalandrotationalcomponentsof,withTRrepresentingthecouplingbetweenthetwocomponents.

PAGE 23

Equations 1.52 and 1.53 canbeinvertedtondthecorrespondingmobilitymatrices,relatingtheparticlevelocitiestothehydrodynamicforcesandtorques, Thus,theresistanceandmobilitymatricescompletelycharacterizethelinearrelationsbetweentheforceandtorquewiththetranslationalandrotationalvelocitiesforparticlemotionatlowReynoldsnumberow. ThesquareresistancematrixcontainssecondranktensorsA;B;andC. Theelementsoftheresistancematrixobeyanumberofsymmetryconditions.Thereciprocaltheoremstatesthattheresistancematrixissymmetric( 6 9 ),sothat ~Bij=Bji; Withe=R=Rbeingtheunitvectoralongthelineofcenters,wecanwrite( 6 )

PAGE 24

Theformofthesescalarfunctionsaredeveloped,whichisafunctionofthedimension-lessseparationa=R.ThedimensionlessfunctionsfordiagonalelementsoftheresistancematrixwillapproachunityforlargeRbecausethescaleswerechosenbyconsideringthesinglesphereresult( 6 ).FromEqs. 1.57 1.62 ,weobtain10independentnon-dimensionalscalarfunctionstodescribetheresistancematrix.Analyticapproximationstothesescalarfunctionscanthenbederivedfromlubricationtheoryforparticlesinnearcontact.

PAGE 25

However,solvingtheinitialvalueproblemforanumberofparticlesofarealisticorderofmagnitude(say,N'1020)isanimpossibletask.Asaresult,onlyaveragescanbecomputedandarerelatedtomacroscopicquantitiesaspressure,temperature,stresses,heatow,etc.;thisisthebasicideaofstatisticalmechanics. Fromstatisticalmechanics,wetalkaboutprobabilitiesinsteadofcertainties;thatis,agivenparticlewillnothaveadenitepositionandvelocity,butonlyprobabilitiesofhavingdierentpositionsandvelocities.Undersuitableassumptions,theinformationrequiredtocomputeaveragesforthesesystemscanbereducedtothesolutionofanequation,theso-calledBoltzmannequation.WithasolutiontotheBoltzmannequation,wecanthendeducemacroscopicuidowfromthemicroscopicmodel. 13

PAGE 26

10 )readsasfollows: wheretheleft-handsiderepresentsthestreamingmotionofthemoleculesalongthetrajectoriesassociatedwiththeexternalforceeldFandC12representstheeectofintermolecular(two-body)collisionstakingmoleculesin/outofthestreamingtrajectory.Hereristhepositionoftheparticle,pmcitslinearmomentum,andFistheforceexperiencedbytheparticleasaresultofintermolecularinteractionsandpossiblyexternalelds. Thecollisionoperatorinvolvesthetwo-bodydistributionfunctionf12thatex-pressestheprobabilityofndingamolecule,say1,aroundr1withvelocityc1andasecondmolecule,say2,aroundr2withvelocityc2,bothattimet.ThedynamicequationforC12canbedeveloped,butbecausethecollisionsarecoupledthisequationcallsintoplaythethree-bodyfunctionC123,whichinturndependsonC1234andsoondownanendlesslineknownastheBBGKYhierarchy,afterBogoliubov,Born,Green,Kirkwood,andYvon( 11 ). TosimplifyEq. 2.1 ,Boltzmannassumedthatthesystemisadilutegasofpoint-likestructurelessmoleculesinteractingviaashort-rangetwo-bodypotential.Undersuchconditions,intermolecularinteractionscanbedescribedsolelyintermsoflocalizedbinarycollisions,wheremoleculesspendsmostoftheirtimeonfreetrajectoriesperfectlyunawareofeachother.Withthispicture,thecollisiontermsplitsintogain(G)andloss(L)components: correspondingtodirect/inversecollisionstakingmoleculesout/inthevolumeele-ment( 10 ).Theparameteristhecollisioncrosssection,whichexpressesthenumberofmoleculeswithrelativespeed~g=g~aroundthesolidangle~.

PAGE 27

ThencomesBoltzmann'sclosureassumption: Thisclosureassumesnocorrelationsbetweenmoleculesenteringacollision(molecularchaosorStosszahlansatz)andistheessentialcomponentthatdescribestheBoltzmannequation.TheBoltzmannequationthentakestheform: ThelefthandsideisNewton'sequationforsingleparticledynamics,whiletherighthandsidedescribesintermolecularcollisionsundertheStosszahlansatzapproximation. wheremisthemolecularmass.Here,theprobabilitydistributionfunctionf(r;c;t)describesthenumberofmoleculesattimet,withpositionrlyingwithinthevelocityspacec.Becauseoftheprobabilisticmeaningoff,thedensityistheexpectedmassperunitvolumeat(r;t)ortheproductofthemolecularmassmbytheprobabilitydensityofndingamoleculeat(r;t),whichturnsouttobethenumberdensityn(r;t) FromEq. 2.5 ,themomentumucanbewrittenasfollows or,usingthefollowingcomponents: Thevelocityucanbeperceivedasthemacroscopicobservationdevelopedfrommolecularmotion;itiszeroforthesteadystateofagasenclosedinaboxatrest.

PAGE 28

Eachmoleculehasitsownvelocitycwhichcanbedecomposedintothesumofuandanothervelocity whichdescribestherandomdeviationofthemoleculefromtheorderedmotionwithvelocityu.Thevelocityvisusuallycalledthepeculiarvelocityortherandomvelocity;itcoincideswithcwhenthegasismacroscopicallyatrest( 12 ).Thus,wehavefromEqs. 2.5 2.8 ,and 2.9 Anotherquantityneededfordevelopingmacroscopichydrodynamicsisthemo-mentumux.Sincemomentumisavector,wethereforeconsidertheowofthejthcomponentofmomentumintheithdirection;givenby Eq. 2.11 showsthatthemomentumuxisdescribedbyasymmetrictensorofsecondorder.Itisknownthatinamacroscopicdescriptiononlyapartofthemicroscopicallyevaluatedmomentumuxwillbeidentiedasconvective,becausetheintegralinEq. 2.11 willbeingeneraldierentfromzeroevenifthegasismacroscopicallyatrest.Inordertondouthowthemomentumuxwillappearinamacroscopicdescription,wehavetousethesplittingofcintomeanvelocityuandpeculiarvelocityv.Thisiswrittenas whereEq. 2.5 and 2.10 havebeenapplied.Themomentumuxthereforebreaksintotwoparts,oneofwhichisrecognizedasthemacroscopicmomentumux,whilethesecondpartisahiddenmomentumuxduetotherandommotionofthemolecules.Forexample,inaxedregionofgas,weobservethatthechangeinmomentumisonlyattributedtothematterthatentersandleavestheregion,whichmeansthesecondparthasnomacroscopicexplanationunlessitisattributedtotheactionofaforceexerted

PAGE 29

Figure2{1: The18velocityvectors(ci)inthelattice-Boltzmannmodel. ontheboundaryoftheregionofinterestbythecontiguousregionofthegas.Inotherwords,theintegralofRvivjmfdcappearsasthecontributiontothestresstensor( 12 ) Here,thetracepii3i=1piiistheisotropicpartofthestresstensor.Morespecically,theuidpressurecanbewrittenasP=pii=3forthecaseatequilibrium. 13 ).ThusthemoregeneralkineticequationsinSec. 2.1 aresimpliedbyconsideringadiscretesetofvelocitiesonly.Thelatticeusedinthissimulationisathreedimensionalcubiclatticeforwhichthereare18velocities,Fig. 2{1 .Thismodelusesthe[100]and[110]directions,withtwicethedensityofparticlesmovinginthe[100]directionsasin[110]directions.Thisgivesthesymmetryconditionnecessarytoensurethatthehydrodynamictransportcoecientsareisotropic.Inthisworkthe18-velocitymodelisaugmentedwithstationaryparticles,whichenablesittoaccountforsmalldeviationsfromtheincompressiblelimit,althoughinsimulationsofstationaryowsitisfoundthatthenumericaldierencesaresmall( 14 ).

PAGE 30

Inthelattice-Boltzmannmethod(LBM),probabilitiessubstitutefordiscreteparticlesintheparticlevelocityfunctionni(r;t).TheLBMthendescribesthetimeevolutionofthediscretizedvelocitydistribution whereiisthechangeinniduetomolecularcollisionsatthelatticenodes,andtisthetimestep.Thisone-particledistributionfunctiondescribesthemassdensityofparticleswithvelocityci,atalatticenoder,andatatimet;r,t,andciarediscrete,whereasniiscontinuous.Thehydrodynamicelds,massdensity,momentumdensityj=u(whereuistheuidowvelocity),andmomentumux,aremomentsofthisvelocitydistribution: Thepost-collisiondistribution,ni(r;t)+i[n(r;t)],ispropagatedforonetime-step,inthedirectionci.Thecollisionoperatori(n)dependsonalltheni'satthenoderepresentedcollectivelybyn(r;t),withtheconstraintsofmassandmomentumconservation.ExactexpressionsfortheBoltzmanncollisionoperatorhavebeenderivedforseverallattice-gasmodels( 15 16 )bymakinguseofthe\molecularchaos"assumption(thattheeectofthecollisionoperatordependsonlyontheinstantaneousstateofanode).However,suchcollisionoperatorsarecomplexandill-suitedtonumericalsimulation.Acomputationallymoreusefulformforthecollisionoperatori(n)canbeconstructedbylinearizingaboutthelocalequilibrium( 17 )neq whereLijarethematrixelementsofthelinearizedcollisionoperator,andi(neq)=0.LiscalculatedbyconsideringthegeneralprinciplesofconservationandsymmetryandtheconstructionoftheeigenvaluesandeigenvectorsofL. Thepopulationdensityassociatedwitheachvelocityhasaweightacithatde-scribesthefractionofparticleswithvelocityciinasystematrest;theseweightsdependonlyonthespeedciandarenormalizedsothatiaci=1.Thevelocitiesciare

PAGE 31

chosensuchthatallparticlesmovefromnodetonodesimultaneously.Foranycubiclattice, wherec=x=t,xisthegridspacing,andC2isanumericalcoecientthatdependsonthechoiceofweights.However,inorderfortheviscousstressestobeindependentofdirection,thevelocitiesmustalsosatisfytheisotropycondition; Theisotropyconditionrequiresthatamulti-speedmodelbeapplied,whichinthiscaseisthe18-velocitymodelforasimplecubiclattice. Asarststep,thecorrectformoftheequilibriumdistributionmustbeestablished.Todeterminetheequilibriumdistribution,thevelocitydistributionfunctionissplitintoalocalequilibriumpartandanon-equilibriumpart Theequilibriumdistributionisacollisionalinvariant(i:e:,i(neq)=0).Theformoftheequilibriumdistributionisconstrainedbythemomentconditionsrequiredtoreproducetheinviscid(Euler)equationswithnon-dissipativehydrodynamicsonlargelengthandtimescales.Inparticular,thesecondmomentoftheequilibriumdistributionshouldbeequaltotheinviscidmomentumuxPI+uu: TheequilibriumdistributioncanbeusedinEqs. 2.20 and 2.21 (c.f.Eq. 2.15 )becausemassandmomentumareconservedduringthecollisionprocess,thus

PAGE 32

ThepressureinEq. 2.22 ,P=c2s,takestheformofanidealgasequationofstatewithadiabaticsoundspeedcs.Thisisalsoadequatefortheliquidphaseifthedensityuctuationsaresmall(i.e.theMachnumberM=u=cs1),sothatrP=c2sr. ForsmallMachnumbers,theequilibriumdistributioncanbeexpandedasapowerseriesinuci=c2s,anditthenfollowsfromEq. 2.22 thattheweightsmustbechosensothatC2=(cs=c)2;i.e.fromEq. 2.17 Asuitableformfortheequilibriumdistributionofthe19-velocitymodelthatsatisesEqs. 2.20 2.22 ,aswellastheisotropycondition( 18 )(Eq. 2.18 ),is 2c4s;(2.25) wherecs=p 3;a1=1 18;ap 36:(2.26) InthiscasethecoecientinEq. 2.18 isC4=(cs=c)4. FromEqs. 2.25 and 2.26 ,theinviscidhydrodynamicequationsarecorrectlyrepro-duced.Theviscousstressescomefrommomentsofthenon-equilibriumdistribution,asintheChapman-EnskogsolutionoftheBoltzmannequation.Thefundamentallimitationofthisclassoflattice-BoltzmannmodelisthattheMachnumberbesmall,lessthan0.3;oursuspensionsimulationsalwayskeepM<0:1. Havingconstructedthelocalequilibriumdistribution,wecometotherelaxationofthenon-equilibriumcomponentofthelocaldistribution.Beyondtheconservationofmassandmomentuminherentinneq,theeigenvalueequationsofthecollisionoperatorfornneqdeveloptheisotropicrelaxationofthestresstensor. Thelinearizedcollisionoperatormustsatisfythefollowingeigenvalueequations; where 2.23 ),andthelasttwoequations

PAGE 33

describetheisotropicrelaxationofthestresstensor;theeigenvaluesandvarerelatedtotheshearandbulkviscositiesandlieintherange2<<0.Equation 2.27 accountsforonly10oftheeigenvectorsofL.Theremaining9modesarehigher-ordereigenvectorsofLthatarenotrelevanttotheNavier-Stokesequations,butwhichdoaecttheboundaryconditionsatthesolid-uidinterfaces.Ingeneraltheeigenvaluesofthesekineticmodesaresetto-1,whichbothsimpliesthesimulationandensuresarapidrelaxationofthenon-hydrodynamicmodes( 19 ). Thecollisionoperatorcanbefurthersimpliedbytakingasingleeigenvalueforboththeviscousandkineticmodes( 20 18 ).Thisexponentialrelaxationtime(ERT)approximation,i=nneqi=,hasbecomethemostpopularformforthecollisionoperatorbecauseofitssimplicityandcomputationaleciency.However,theabsenceofacleartimescaleseparationbetweenthekineticandhydrodynamicmodescansometimescausesignicanterrorsatsolid-uidboundaries,andthusthemoreexiblecollisionoperatordescribedbyEq. 2.27 ispreferred. A3-parametercollisionoperatorisusedinthesuspensionsimulations,allowingforseparaterelaxationofthe5shearmodes,1bulkmode,and9kineticmodes.Thepost-collisiondistributionn?i=ni+iiswrittenasaseriesofmoments(Eq. 2.15 ),asseenwiththeequilibriumdistribution(Eq. 2.25 ), 2c4s:(2.28) Thezeroth()andrst(j=u)momentsarethesameasintheequilibriumdistri-bution(Eq. 2.25 ),butthenon-equilibriumsecondmomentneqismodiedbythecollisionprocess,accordingtoEq. 2.27 : 3(1+v)(neq:I)I;(2.29) whereneq=eq.Thekineticmodescanalsocontributetothepost-collisiondistribution,buttheeigenvaluesofthesemodesaresetat1,sothattheyhavenoeectonn?i.Equation 2.29 with=v=1isequivalenttotheERTmodelwith=1;for<1,thekineticmodesrelaxmorerapidlythantheviscousmodes,whichistheproperlimitforhydrodynamics.

PAGE 34

Themacrodynamicalbehaviorarisesfromthelattice-Boltzmannequationwiththemulti-time-scaleanalysisapplied( 15 14 ).TheNavier-Stokesequations, arerecoveredinthelowvelocitylimit,withviscosities 2andv=c2st2 3v+1 3:(2.31) Thefactorsof1=2and1=3,whichappearinthedenitionsoftheshearandbulkviscositiesservetocorrectfornumericaldiusion,sothatviscousmomentumdiusesattheexpectedspeed. Inthepresenceofanexternallyimposedforcedensityf,forexampleapressuregradientoragravitationaleld,thetimeevolutionofthelattice-Boltzmannmodelincludesanadditionalcontributionfi(r;t), Thisforcingtermcanbeexpandedinapowerseriesinthevelocity, 2c4st:(2.33) TheexpressionrelatingtheforcedensitytothechangeinvelocitydistributionwasdeterminedbymatchingthemacroscopicdynamicsfromEq. 2.32 totheNavier-Stokesequations( 14 ).Moreaccuratesolutionstothevelocityeldareobtainedifitincludesaportionofthemomentum( 14 )addedtoeachnode, 2ft:(2.34) 2.1 forrealisticnonequilibriumsituations,wemustrelyuponapproximationmethods,inparticular,perturbationprocedures.Inordertodothis,welookforaparameterwhichweconsidertobesmall.Therefore,inobtainingmacroscopichydrodynamicequations,we

PAGE 35

introducethenondimensionalKnudsennumberdenotedbyKn Kn==l=d;(2.35) wherelisthemeanfreepathbetweenmolecularcollisions,anddisatypicallengthscale.Thus,Kn!0correspondstoafairlydensegasandKn!1toafreemolecularow(i:e:,aowwheremoleculeshavenegligibleinteractionswitheachother). TheChapman-Enskogexpansion,then,considersmultipletimescales,ofordersand2,torecovertheNavier-Stokesequation.Inthissense,theNavier-Stokesequationdescribestwokindsofprocesses,convectionanddiusion,whichactontwodierenttimescales( 21 ).Ifweconsideronlytherstscale,weobtainthecompressibleEulerequation;ifwemoveontothesecondoneweobtaintheNavier-Stokesequationwhilelosingthecompressibilityeect.Ifwewantbothcompressibilityanddiusion,wehavetokeepbothscalesatthesametimeandthinkofnintermsoftimeandspacederivatives @t=@n @t1+2@n @t2 @r=@n @r1: Thisenablesustointroducetwodierenttimevariablest1=t;t2=2tandanewspacevariabler1=r,suchthattheuiddynamicalvariablesarefunctionsofr1;t1;t2;. Thecompatibilityconditionsattherstordergivethatthetimederivativesoftheuiddynamicvariableswithrespecttot1isdeterminedbytheEulerequationbutthederivativeswithrespecttot2aredeterminedonlyatthenextlevelandaregivenbythetermsoftheNavier-Stokesequationdescribingtheeectofviscosityandheatconductivity( 21 ).ThetwocontributionsaretobeaddedasindicatedinEq. 2.36 toobtainthefulltimederivativeandthuswritetheNavier-Stokesequation. WhiletheChapman-EnskogexpansionwasoriginallyappliedtothecontinuousBoltzmannequationinEq. 2.1 ,herewederivetheNavier-Stokesequationusingthe

PAGE 36

simpliedlatticeBoltzmannequationrepresentedasfollows: with i(r;t)=1 ThecollisionoperatorhereisfromthesinglerelaxationtimeBhatnagar-Gross-Krook(BGK)modelwhereneqirepresentstheequilibriumdistributionandisatypicaltime-scaleassociatedwithcollisionalrelaxationtothelocalequilibrium. ToevaluatetheLBGKmodel,rstaTaylorseriesexpansionforni(r+cit;t+t)isperformed, NextweintroducethetwotimescalesandonespatialscaleasnotedinEq. 2.36 and 2.37 .Thedistributionfunctionisalsoexpandedaccordingly: SubstituteEq. 2.40 intotheLHSofEq. 2.38 withEq. 2.39 gives t@ni SubstituteEqs. 2.36 2.37 ,and 2.41 intoEq. 2.42 gives (2.43) (2.44) 1

PAGE 37

Arrangingtermsinascendingordersofgives Itisnotedthattheequilibriumdistribution,alongwithEq. 2.46 ,shouldsatisfythefollowingconstraints: wheremisthetotalnumberofdirectionsinthelatticemodel,andthemolecularmassissettounity.WiththeseconstraintsappliedtoEq. 2.41 ,theconstraintsforthe1stand2ndordercomponentsinbecome 0=mXi=0n(1)i=mXi=0n(2)i 0=mXi=0cin(1)i=mXi=0cin(2)i: UsingtheaboverelationinEq. 2.47 givesthefollowingreduction @t1mXi=0n(0)i+@ @r1mXi=0cin(0)i=0!@ @t1+@u Furthermore,thesecondmomentisobtainedby @t1mXi=0cin(0)i+@ @r1mXi=0cicin(0)i=0!@u

PAGE 38

RecallfromthederivationinEq. 2.12 usingthepeculiarvelocity,theviscousstressintermsofthedistributionfunctionncanbedeterminedtobe Thus,torstorderin,fromEq. 2.53 thecontinuityequationisobtainedandfromEqs. 2.54 and 2.55 theEulerequationisdetermined Next,theO(2)termsareanalyzed.FromEq. 2.48 ,therelationsinEqs. 2.49 2.51 2.53 ,and 2.54 areappliedtogetarelationintermsoft2,and (2.57) @t1mXi=0n(1)i=0; @t2mXi=0n(0)i=@ @t2; @r1mXi=0cin(1)i=0: SubstitutionofEqs. 2.53 and 2.54 intothe4th,5th,and6thtermsofEq. 2.57 gener-atesandoverallvalueofzero.Thus,thezerothmomentofO(2)becomes @t2=0: Forthemomentumconservation,Eq. 2.52 isappliedtogive (2.62)

PAGE 39

Here,applyingEqs. 2.50 and 2.52 gives @t1mXi=0cin(1)i=0; @t2mXi=0cin(0)i=@u ApplyingEq. 2.47 intothe3rdcomponentofEq. 2.62 gives @r1mXi=0cicin(1)i=t@ @r1mXi=0cici"@n(0)i The4thcomponentinEq. 2.62 canbeevaluatedfurtherbysubstitutiontogive t @t1@P(0) =t @t1@ @r1uu+c2s @r1@ @t1; whereO(u2)termsareneglectedduetolowMachnumberows.Herecsrespresentsthespeedofsoundwhichischaracterizedbythepressureatequilibrium(i:eP=RT=c2s).SubstitutionfromEq. 2.53 gives =t @r1@ @r1(u)=t The5thcomponentofEq. 2.62 isconsideredinconjunctionwiththe1stterminEq. 2.65 toget (1)t@ @t1@ @r1mXi=0cicin(0)1=(1)t@ @t1@ @r1P(0)(1)tc2s@ @r1@ @r1(u)=(1)tc2srr(u) (2.70)

PAGE 40

andthe6thcomponentofEq. 2.62 withthesecondterminEq. 2.65 gives 2t@ @r1@ @r1mXi=0cicicin(0)i(usingdefn:fromEq: )1 2t c2s@ @r1@ @r1mXi=0acicicici(ciu) (2.71) 2tc2s@ @r1@ @r1(u)+2@ @r1@ @r1(u)=1 2tc2sr2(u)+2rr(u): ThenewterminEq. 2.71 istakenfromtheistropicconditionforthislatticemodel.FinallycombiningEqs. 2.64 2.69 2.70 ,and 2.72 intoEq. 2.62 gives @t2(u)=1 2tc2sr2(u)+rr(u): ThusthedynamicshearandbulkviscositiesforthelatticeBGKmodelare 2tc2s: Theendresultistosumoverallordersof.WithEq. 2.53 addedtoEq. 2.61 ,thecontinuityequationisdetermined: @t+@u Eqs. 2.54 ,withEq. 2.55 ,addedtoEq. 2.73 gives 2tc2sr2(u)+rr(u): Intheincompressiblelimit, (2.77) Thus,thecontinuumNavier-Stokesequationhasbeenderived.Intheend,thebasicresultoftheChapman-EnskogprocedureisthatthemacroscopicdescriptionofauidisrecoveredbyasuitableexpansionoftheBoltzmannequation.

PAGE 41

Figure2{2: Locationofboundarynodesforacurvedsurface(a)andtwosurfacesinnearcontact(b).

PAGE 42

19 ).Solidparticlesaredenedbyasurface,whichcutssomeofthelinksbetweenlatticenodes.Fluidparticlesmovingalongtheselinksinteractwiththesolidsurfaceatboundarynodesplacedhalfwayalongthelinks.Thusweobtainadiscreterepresentationoftheparticlesurface,whichbecomesmoreandmorepreciseastheparticlegetslarger(Fig. 2{2 ). Themotionofaparticleisdeterminedbytheforcesandtorquesactingonit.Theseforcesconsistofexternal,non-hydrodynamicforces,andhydrodynamicforces.Non-hydrodynamicforces,suchasattractionand/orrepulsionbetweensurfacesduetoaninter-particlepotential,canbecalculatedseparatelyandincorporatedintotheequationsofmotion.Incontrast,thehydrodynamicforcesarecalculatedfromtheresultantvelocitydistributions,ni. Stationarysolidobjectsareintroducedintothesimulationbythe\bounce-back"collisionruleatthespeciedboundarynodes,inwhichincomingdistributionsarereectedbacktowardsthenodestheycamefrom.Surfaceforcesarecalculatedfromthemomentumtransferateachboundarynodeandsummedtogivetheforceandtorqueoneachobject( 19 ).Whilemoreaccuratemethodsthatrequireknowledgeoftheshapeoftheparticlesurfacehavebeeninvestigated( 22 23 24 25 26 27 ),thebounce-backrulecanbeappliedtosurfacesofarbitraryshape,withoutadditionalcomplications.However,forsecond-orderconvergencethebounce-backrulerequiresacalibrationofthehydrodynamicradius,whichisnotalwaysconvenient. Theboundarynodesarelocatedmidwaybetweeninterior(solid)andexterior(uid)nodes( 28 29 ).Thenormalcollisionrulesarecarriedoutatalluidnodes,andaugmentedbybounce-backrulesatthemidpointsoflinksconnectingthelatticenodes.ThisisdepictedbythearrowsinFig. 2{2 alinkingthenormaldistributionbetweenuid(circles)andboundarynodes(squares).InPoiseuilleowthe\linkbounce-back"rulegivesvelocityeldsthatdeviatefromtheexactsolutionbyaconstantslipvelocity,us=uLBEuExact,proportionaltoL2( 30 )(Lbeingthechannelwidth);

PAGE 43

Figure2{3: Movingboundarynodeupdate.Thesolidcirclesareuidnodes,squaresareboundarynodeswiththeassociatedvelocityvector,andopencirclesareinteriornodes. whereuc=L2jrpj=8istheexactvelocityatthecenterofthechannelandthecoecientdependsontheeigenvaluesofthecollisionoperator. Inthepastthelatticenodesweretreatedoneithersideoftheboundarysurfaceinanidenticalfashion,sothatuidllsthewholevolumeofspace,bothinsideandoutsidethesolidparticles.Becauseoftherelativelysmallvolumeinsideeachparticle,theinterioruidquicklyrelaxestorigid-bodymotion,characterizedbytheparticlevelocityandangularvelocity,andcloselyapproximatesarigidbodyofthesamemass( 14 ).However,atshorttimestheinertiallagoftheuidisnoticeable,andthecontributionoftheinterioruidtotheparticleforceandtorquereducesthestabilityoftheparticlevelocityupdate.Forthesereasons,withthesuggestionsinRef.( 31 ),theinterioruidisremovedfromtheparticlerepresentation.AsomewhatdierentimplementationofthesameideaisdescribedinRef.( 32 ). Themovingboundarycondition( 19 )withoutinterioruid( 31 )isimplementedasshowninFig. 2{3 .Wetakethesetofuidnodesrjustoutsidetheparticlesurface,andforeachnodeallthevelocitiescbsuchthatr+cbtliesinsidetheparticlesurface.Eachofthecorrespondingpopulationdensitiesisthenupdatedaccordingtoasimplerulewhichtakesintoaccountthemotionoftheparticlesurface( 19 );

PAGE 44

wherenb(r;t)isthepost-collisiondistributionat(r;t)inthedirectioncb,andcb0=cb.Forastationaryboundarythepopulationissimplyreectedbackinthedirectionitcamefromwithub=0( 15 33 ).Thelocalvelocityoftheparticlesurface, isdeterminedbytheparticlevelocityU,angularvelocity,andcenterofmassR;rb=r+1 2cbtisthelocationoftheboundarynode.Themeandensity0inEq. 2.80 isusedinsteadofthelocaldensity(r;t)sinceitsimpliestheupdateprocedure.Thedierencesbetweenthetwomethodsaresmall,ofthesameorder(u3)astheerrortermsinthelattice-Boltzmannmodel.Testcalculationsshowthatevenlargevariationsinuiddensity(upto10%)haveanegligibleeectontheforce(lessthan1partin104). Asaresultoftheboundarynodeupdates,momentumisexchangedlocallybetweentheuidandthesolidparticle,butthecombinedmomentumofsolidanduidisconserved.TheforcesexertedattheboundarynodescanbecalculatedfromthemomentumtransferredinEq. 2.80 2t)=x3 Theparticleforcesandtorquesarethenobtainedbysummingf(rb)andrbf(rb)overalltheboundarynodesassociatedwithaparticularparticle.Itcanbeshownanalyt-icallythattheforceonaplanarwallinalinearshearowisexact( 19 ),andseveralnumericalexamplesoflattice-BoltzmannsimulationsofhydrodynamicinteractionsaregiveninRef.( 34 ). Tounderstandthephysicsofthemovingboundarycondition,onecanimagineanensembleofparticles,movingatconstantspeedcb,impingingonamassivewallorientedperpendiculartotheparticlemotion.Thewallitselfismovingwithvelocityubcb.Thevelocityoftheparticlesaftercollisionwiththewalliscb+2ubandtheforceexertedonthewallisproportionaltocbub.Sincethevelocitiesinthelattice-Boltzmannmodelarediscrete,thedesiredboundaryconditioncannotbeimplementeddirectly,butwecaninsteadmodifythedensityofreturningparticlesso

PAGE 45

thatthemomentumtransferredtothewallisthesameasinthecontinuousvelocitycase.Itcanbeseenthatthisimplementationoftheno-slipboundaryconditionleadstoasmallmasstransferacrossamovingsolid-uidinterface.Thisisphysicallycorrectandarisesfromthediscretemotionofthesolidsurface.Thusduringatimestepttheuidisowingcontinuously,whilethesolidparticleisxedinspace.Iftheuidcannotowacrossthesurfacetherewillbelargearticialpressuregradients,arisingfromthecompressionandexpansionofuidnearthesurface.Forauniformlymovingparticle,itisstraightforwardtoshowthatthemasstransferacrossthesurfaceinatimestept(Eq. 2.80 )isexactlyrecoveredwhentheparticlemovestoitsnewposition.Forexample,eachuidnodeadjacenttoaplanarwallhas5linksintersectingthewall.IfthewallisadvancingintotheuidwithavelocityU,thenthemassuxacrosstheinterface(fromEq. 2.80 )is0U.Apartfromsmallcompressibilityeects,thisisexactlytherateatwhichuidmassisabsorbedbythemovingwall.Forslidingmotion,Eq. 2.80 correctlypredictsnonetmasstransferacrosstheinterface. 35 36 ).Thusforlargechannels,thehydrodynamicboundaryisdisplacedbyanamountfromthephysicalboundary,whereisindependentofchannelwidthbutdependsontheorientationofthechannelwithrespecttotheunderlyinglattice.Convexbodiessampleavarietyofboundaryorientations,sothatitisnotpossibletomakeananalyticaldeterminationofthedisplacementofthehydrodynamicboundaryfromthesolidparticlesurface.Nevertheless,thedisplacementoftheboundarycanbedeterminednumericallyfromsimulationsofowaroundisolatedparticles.Byconsideringthesizeoftheparticlestobethehydrodynamicradius,ahy=a+,ratherthanthephysicalradiusa,approximatesecond-orderconvergencecanbeobtained,evenfordensesuspensions( 34 ). Thehydrodynamicradiuscanbedeterminedfromthedragonaxedsphereinapressure-drivenow( 34 ).Galileaninvariancecanbedemonstratedbyshowingthatthesameforceisobtainedforamovingparticleinastationaryuid( 34 ).Sincethe

PAGE 46

Figure2{4: DragforceFasafunctionoftime,normalizedbythedragforceonanisolatedsphere,F0=6aU.TimeismeasuredinunitsoftheStokestime,tS=a=U. particlesamplesdierentboundarynodemapsasitmovesonthegrid,itisimportanttosampledierentparticlepositionswhendeterminingthehydrodynamicradius,especiallywhentheparticleradiusissmall(<5x).Ratherthanaveragingovermanyxedcongurations,wechosetohavetheparticlemoveslowlyacrossthegrid,atconstantvelocity,samplingdierentboundarynodemapsasitgoes.Thechangingboundarynodemapsleadtouctuationsinthedragforce,asshowninFig. 2{4 .TheforcehasbeenaveragedoveraStokestime,tS,sothattherelativeuctuationsinforcearecomparabletotherelativeuctuationsinvelocityofaneutrallybuoyantparticleinaconstantforceeld.Theforceuctuations,F=p 2{1 ).Moresophisticatedboundaryconditionshavebeendevelopedusingnite-volumemethods( 37 38 )andinterpolation( 23 25 39 ).Bothmethodsreducetheforceuctuationsbyatleastanorderofmagnitudefromthoseobservedhere,butevenwiththesimplebounce-backscheme,theuctuationsinforcecanbereducedbyanappropriatechoiceofparticleradius.Wehavenoticed

PAGE 47

Table2{1: VarianceinthecomputeddragforceF=p thatuctuationsinparticleforcearestronglycorrelatedwithuctuationsinparticlevolume.Thuswechoosetheradiusoftheboundarynodemapsoastominimizeuctuationsinparticlevolumeforrandomlocationsonthegrid.ItcanbeseenfromTable 2{1 thatatwofoldreductionintheforceuctuationsispossiblebythisprocedure,althoughforsucientlylargeparticlesthedierenceisminimal.AsetofoptimalparticleradiiisgiveninTable 2{2 Thebounce-backruleleadstoavelocityeldintheregionoftheboundarythatisrst-orderaccurateinthegridspacingx.Thehydrodynamicboundary(thesurfacewheretheuidvelocityeldmatchesthevelocityoftheparticle)isdisplacedfromtheparticlesurfacebyaconstant,(Fig. 2{5 ),thatdependsontheviscosityoftheuid( 34 ).Fortherangeofkinematicviscositiesusedinthiswork,1=61=1200,variesfrom0to0:5x(Table 2{2 );thedimensionlesskinematicviscosity=t=x2.Forsmallparticles(a<5x),alsodependsweaklyontheparticleradius(Table 2{2 ).Althoughthedierencebetweenthehydrodynamicboundaryandthephysicalboundaryissmall,itisimportantinobtainingaccurateresultswithcomputationallyusefulparticlesizes.Calibrationofthehydrodynamicradiusleadstoapproximatelysecond-orderconvergencefromaninherentlyrst-orderboundarycondition;itwillnotbenecessarywhenmoreaccurateboundaryconditionsareimplemented. Thehydrodynamicradii(ahy)inTable 2{2 weredeterminedfromthedragforceonasinglesphereinaperiodiccubiccell( 34 ).TheReynoldsnumberinthesesimulations Table2{2: Hydrodynamicradiusahy(inunitsofx)forvariousuidviscosities;aistheinputparticleradius.

PAGE 48

Figure2{5: Actual(solidlines)andhydrodynamic(dashedlines)surfacesforaparticlesettlingontoawall. (<0:2)wassucientlysmalltohaveanegligibleeectonthedragforce.Thetimeaveragedforcewasusedtodeterminetheeectivehydrodynamicradiusforthreedierentkinematicviscosities:=1=6,=1=100,and=1=1200.Ineachcasethecelldimensionswere10timestheparticleradiusandthecorrectionsforperiodicboundaryconditions(about40%)weremadeasdescribedinRef.( 34 ). Thedierencebetweenthehydrodynamicradiusandtheinputradius,=ahya,isindependentofparticlesizeforradiia>5x,anditsmagnitudeincreaseswithdecreasingkinematicviscosity( 14 ).Thekinematicviscosity=1=6givesahydrodynamicradiusthatisthesameastheinputradius(forsucientlylargeparticles),andisthenaturalchoiceforsimulationsatlowReynoldsnumber.AthigherReynoldsnumbersalowerviscosityisnecessarytomaintainincompressibility( 14 34 ),andforaccurateresultsitisthenessentialtousethecalibratedhydrodynamicradius. Inordertoimplementthehydrodynamicradiusinamulti-particlesuspension,alldistancecalculationsarebasedonthehydrodynamicradius(asshowninFig. 2{5 );theinputradiusaisonlyusedtodeterminethelocationoftheboundarynodes.Itshouldbenotedthatnotallcombinationsofparticleradiusandviscositycanbeused.Table 2{2 indicatesthatparticleradiilessthan3xcannotbeusedwithakinematicviscosity=1=6,sincethehydrodynamicradiusisthenlessthantheinputradius.

PAGE 49

mF(t) (2.83) IT(t) (2.84) hasbeenfoundtobeunstable( 34 )unlesstheparticleradiusislargeortheparticlemassdensityismuchhigherthanthesurroundinguid.Inpreviouswork( 34 )theinstabilitywasreduced,butnoteliminated,byaveragingtheforcesandtorquesovertwosuccessivetimesteps.Subsequently,animplicitupdateoftheparticlevelocitywasproposed( 40 )asameansofensuringstability.Herewepresentageneralizedversionofthatidea,whichcanbeadaptedtosituationswheretwoparticlesareinnearcontact. Theparticleforceandtorquecanbeseparatedintoacomponentthatdependsontheincomingvelocitydistributionandacomponentthatdepends,viaub,ontheparticlevelocityandangularvelocity(Eqs. 2.81 and 2.82 ); Thevelocityindependentforcesandtorquesaredeterminedatthehalf-timestep 2t)=x3 2t)=x3 wherethesumisoveralltheboundarynodes,b,describingtheparticlesurface,withcbrepresentingthevelocityassociatedwiththeboundarynodebandpointingtowardsthe

PAGE 50

particlecenter.Thematrices (2.90) (2.92) arehigh-frequencyfrictioncoecients,anddescribetheinstantaneousforceonaparticleinresponsetoasuddenchangeinvelocity. Themagnitudeofthefrictioncoecientscanbereadilyestimated,therebyestablishingboundsonthestabilityofanexplicitupdate.Apartfromirregularitiesinthediscretizedsurface,FUandTarediagonalmatrices,whileF=TU=0.ForanodeadjacenttoaplanarwallPiacic2i=5 18c2wherethesumisoverthe5directionsthatcrossthewall.Thenumberofsuchnodesisapproximately4a2=x2,sothat Similarly, Theseestimatesofthetranslationalandrotationalfrictioncoecientsarewithin20%and50%ofnumericallycomputedvalues,respectively.ThestabilitycriterionforanexplicitupdateFUt=m<2thenreducestoasimpleconditioninvolvingtheparticleradiusandmassdensity: 5 3fx sa<2:(2.95) Thecorrespondingconditionforthetorqueleadstothesamestabilitycriterion I5 3fx sa<2;(2.96) whereaswithinterioruidthenumericalfactorswere6and10respectively( 34 ),showingthatinterioruiddestabilizesanexplicitupdate.

PAGE 51

ThefrictioncoecientsinEqs. 2.85 and 2.86 areessentiallyconstant,uctuatingslowlyintimeastheparticlemovesontheunderlyinggrid;thustheparticlevelocitiescanbeupdatedassumingthesefrictionmatricesareconstant.Theequationsofmotioncanthenbewritteninmatrixformas whereisaparameterthatcontrolsthedegreeofimplicitness.Anexplicitupdate( 34 )correspondsto=0,animplicitupdate( 40 )correspondsto=1,andasecondordersemi-implicitupdatecorrespondsto=1 2.Theexplicit,implicit,andsemi-implicitupdatesevaluatethevelocity-dependentforceatt,t+t,andt+1 2trespectively.Inpracticewendonlysmalldierencesbetweensemi-implicitandfullyimplicitmethodsandweusethefullyimplicitmethod(=1)inthiswork.Theboundarynodemapisupdatedinfrequently(every10-100timesteps)andthe66matrixinversionneedonlybedonewhenthemapisupdated.WenotethatinthelimitofFUt=m>>1onlythefullyimplicit(=1)versionofEq. 2.97 reducestothesteadystateforceandtorquebalance,F0FUU(t+t)=T0T(t+t)=0.Thesemi-implicitmethod(=1 2)producesanoscillatingsolutionandtheexplicitmethod(=0)adivergingsolution. Animplicitupdateoftheparticlevelocitiesrequirestwopassesthroughtheboundarynodes.OntherstpassthepopulationdensitiesareusedtocalculateF0andT0.Theimplicitequations(Eq. 2.97 )arethensolvedforU(t+t)and(t+t)forthegivenimplicitparameter.Thesevelocitiesarethenusedtocalculatethenewpopulationdensitiesinasecondsweepthroughtheboundarynodes.Thecomputationaloverheadincurredbytheboundarynodeupdatesistypicallylessthan100%,evenathighconcentrations.

PAGE 52

19 34 )arebeingincreasinglyusedtocalculatehydrodynamicinteractions( 31 41 42 43 44 45 46 ),bydirectnumericalsimulationoftheuidmotiongeneratedbythemovinginterfaces.However,astwoparticlesapproacheachotherthecalculationofthehydrodynamicforcebreaksdowninthegapbetweentheparticles,typicallywhentheminimumdistancebetweenthetwosurfacesisoftheorderofthegridspacing.Forexample,itisimpracticaltousesucientmeshresolutiontoresolvetheowinthesmallgapsthatcanresultfromanimposedshearow.Athighshearratestherheologyofasuspensionofspheresisqualitativelyaectedbylubricationforcesbetweenparticlesseparatedbygapslessthan0:01a,whereaistheparticleradius( 47 48 ).Asimulationatthisresolution(107gridpointsperparticle)isunfeasibleformorethanafewparticles.Thenumberofgridpointscanbereducedbyusingbodyttedcoordinates( 49 )oradaptivemeshes( 23 25 ),butthesmallzonesizeinthegapreducesthetimestepthatcanbeusedtointegratetheoweld( 50 ).Theproblemisexacerbatedbytheuniformgridusedinlattice-Boltzmannsimulations,butitshouldbenotedthatsimilartechniques,usingparticlesembeddedinregulargrids,arebecomingincreasinglypopularincomputationaluiddynamics( 51 52 ).Someaspectsofthisworkmaythereforebeapplicabletosuchmethodsaswell. SimulationsofhydrodynamicallyinteractingparticlestypicallyusemultipoleexpansionsoftheStokesequations( 53 54 ).Againthecalculationsbreakdownwhentheparticlesareclosetocontact,unlessanimpracticalnumberofmultipolemomentsareused( 55 ).Asolutiontothisproblemistocalculatethelubricationforcesbetweenpairsofparticlesforarangeofsmallinterparticlegaps,andthenpatchthesesolutionsontothefrictioncoecientscalculatedfromthemultipoleexpansion.Themethodexploitsthefactthatlubricationforcescanbeaddedpair-by-pair,and 40

PAGE 53

hasbeenshowntoworkwellinpractice( 53 ).Asimpliedversionofthisapproachhasalreadybeenadoptedtoincludenormallubricationforcesinlattice-Boltzmannsimulations( 56 ).Hereweextendourpreviousworktoincludeallcomponentsofthelubricationforceandtorque,andtestthenumericalschemefortheinteractionsbetweenasphericalparticleandaplanarwall.Wendthatthehydrodynamicinteractionscanbewellrepresentedovertheentirerangeofseparationsbypatchingonlythemostsingularcomponentsofthelubricationforceontotheforcecalculatedfromthelattice-Boltzmannmodel.ThisissimplerthantheStokesiandynamicsapproach,wherethepatchiscalculatedateveryseparation. Inthischapter,weproposeacomprehensivesolutiontothetechnicaldicultiesthatarisewhenparticlesareclosetocontact,inparticularthelossofmassconserva-tion.Thebulkofournumericalresultsareaseriesofdetailedtestsofthehydrody-namicinteractionsbetweentwosurfacesinnearcontact.Wedemonstratethatafterincludingcorrectionsforthelubricationforcesweobtainaccurateresultsoverawiderangeofuidviscosities.Finally,wedescribeanecientimplicitalgorithmforupdat-ingtheparticlevelocitieseveninthepresenceofstilubricationforces.Anexplicitsolutionofthesedierentialequationsrequireseitherthattheparticlesaremuchdenserthanthesurroundinguid( 34 ),orthatverysmalltimestepsareusedtoupdatetheparticlevelocities.Ontheotherhand,iftheparticlevelocitiesareupdatedimplicitly,theresultingsystemoflinearequationsseverelylimitsthenumberofparticlesthatcanbesimulated.Wedescribewhatwecalla\cluster-implicit"method,wherebygroupsofcloselyinteractingparticlesaregroupedinindividualclustersandinteractionsbe-tweenparticlesinthesameclusterareupdatedimplicitly,whileallotherlubricationforcesareupdatedexplicitly.Inuidizedsuspensionsclusterstypicallycontainlessthan10particles,evenathighconcentrations,sothattheimplicitupdatesareasmalloverhead.Oursimulationsecientlycopewithclustersofseveralhundredparticlesbyusingconjugate-gradientmethods,andonlyslowdownifthenumberofparticlesintheclusterexceeds103.

PAGE 54

2{2 b),leadingtoalossofmassconserva-tion.Thishappensbecauseboundaryupdatesateachlinkproducesamasstransfer(2acb0ubcb=c2s)x3acrossthesolid-uidinterface,whichisnecessarytoaccommodatethediscretemotionoftheparticlesurface(seeSection 2.2.2 ).Thetotalmasstransferinoroutofanisolatedparticle, M=2x30 regardlessoftheparticle'ssizeorshape. AlthoughthesumsPbacbcbandPbacbrbcbarezeroforanyclosedsurface,whentwoparticlesareclosetocontactsomeoftheboundarynodesaremissingandthesurfacesarenolongerclosed.InthiscaseM6=0andmassconservationisnolongerensured.Twoparticlesthatremainincloseproximityneverreachasteadystate,nomatterhowslowlytheymove,sinceuidisconstantlybeingaddedorremoved,dependingontheparticlepositionsandvelocities.Ifthetwoparticlesmoveasarigidbodymassconservationisrestored,butingeneralthisisnotthecase. Theaccumulationorlossofmassoccursslowly,andinmanydynamicalsimula-tionsituctuateswithchangingparticlecongurationbutshowsnolong-termdrift.However,weenforcemassconservation,particle-by-particle,byredistributingtheexcessmassamongtheboundarynodes(c:f:Eq. 2.80 ) A(3.2) whereA=x3Pbacb0. Theforceandtorquearisingfromthisredistributionofmassissmall,butnotexactlyzero; F=x30 AXbacbcb# T=x30 AXbacbrbcb#:

PAGE 55

Theycanbecompactlyincludedbyaredenitionofthefrictioncoecients,Eqs. 2.89 2.92 ,replacingcbandrbcbbytheirdeviationfromthemean, sothatcb!cb 2.85 and 2.86 ,evenwhenmassisredistributed. Forparticlesclosetocontact,thelubricationforce,torque,andstressletcanbecalculatedfromasumofpairwise-additivecontributions( 53 ),andifweconsideronlysingularterms,theycanbecalculatedfromtheparticlevelocitiesalone( 57 ).Inthegrand-resistance-matrix( 58 )formulation whereF2=F1,U12=U1U2,andthefrictionmatricesareasgiveninRef.( 58 ).Wehavemadefulluseofthesymmetriesinherentinthefrictionmatrices,butwithoutassumingthattheparticleradiiarethesame.Mostimportantly,theexternalowelddoesnotenterintothelubricationcalculation,whichremovestheneedforestimatesofthelocaloweld.

PAGE 56

Wehavenotedinpreviouslattice-Boltzmannsimulations( 14 56 )thatthecalcu-latedforcesfollowtheStokesowresultsdowntoaxedseparation,around0:5x,andremainroughlyconstantthereafter.ThesolidsymbolsinFig. 3{1 ,forexample,showthisbehaviorforthenormalforcebetweenasphericalparticleandaplanewall.ThissuggestsalubricationcorrectionoftheformofadierencebetweenthelubricationforceathandtheforceatsomecutodistancehN;i.e. whereU12=U1U2,h=jR12ja1a2isthegapbetweenthetwosurfaces,andtheunitvector^R12=R12=jR12j. ThefrictioncoecientsinEq. 3.6 areallproductsofascalarfunctionofthegapbetweentheparticles,either1=horln(1=h),multipliedbyapolynomialoftheunitvectorconnectingtheparticlecenters( 58 ).Fortwospheresofarbitrarysize,thereareatotalof15independentscalarcoecients,whichfallinto3classes.AgainusingthenotationofRef.( 58 )theseareXA11,XG11,XG22(normalforce);YA11,YB11,YB22,YG11,YG22(tangentialforce);andYC11,YC12,YC22,YH11,YH12,YH21,YH22(rotation).WeimplementourlubricationcorrectionbycalculatingamodiedformofeachscalarcoecientasinEq. 3.7 ;forexample whichvanishesatthecut-odistanceh=hN.Weallowfordierentcutodistancesforeachofthe3typesoflubricationinteraction. 2.2.2 )withtheexceptionoftheresultsfora=x=4:5,whichweregeneratedbeforetheoptimumradius(4:8x)wasdetermined.Thehydrodynamicradii,ahy(a;),whichareusedtodeterminethe

PAGE 57

Figure3{1: NormalforceonaparticleofinputradiusasettlingontoahorizontalplanarsurfaceatzeroReynoldsnumber,withF0=6ahyU.Solidsymbolsindicate:=1=6(circles),=1=100(triangles),and=1=1200(squares).Resultsincludingthelubricationcorrectionareshownbytheopensymbols. positionsofthelubricatingsurfaces,weretakenfromTable 2{2 .Thelocationoftheplanarwallwasshiftedby(),correspondingtothea!1limitinTable 2{2 (Fig. 2{5 ).Inthiswayweensurethatthelubricatingsurfacesareinthesamepositionasthehydrodynamicboundariesinthelattice-Boltzmannsimulations.Theunitcellisperiodicinfourdirectionswithatopandbottomwall,andcelllengthof10timestheparticleradius,whichissucientlylargethattheeectsofperiodicimageswerenegligible.Thesimulationdeterminesthesteadystateforceandtorqueontheparticleforagivenvelocityandangularvelocity,whichwasthenusedtocalculatethefrictioncoecientasafunctionofthegaphfromthewall.SimulationresultsforthefrictionalforcesandtorquesareshowninFigs. 3{1 3{4 forthreedierentuidviscosities=1=6,1=100,and1=1200. Thenormalforceshowsthetrenddiscussedinsection 3.1.2 foreachparticlesizeanduidviscosity(Fig. 3{1 ).Thesimulatedforce(solidsymbols)followstheexactresult(solidline)downtoaninterparticlegap,hN
PAGE 58

Figure3{2: TangentialforceonaparticlesettlingnexttoaverticalplanarsurfaceatzeroReynoldsnumber.SimulationdatafordragforcedividedbyF0=6ahyU,isindi-catedbysolidsymbols.Resultsincludingthelubricationcorrectionareshownbytheopensymbols. size.Forlargerparticlesthelattice-Boltzmannmethodcapturesprogressivelymoreofthelubricationforce,butevenfora=8:2xtherearenoticeabledeviationsforh=ahy<0:01.Thesimulationreproducesmoreofthelubricationforceatsmallerviscositiesbecausetheshiftinthehydrodynamicradiusdelaysthecontactofthepar-ticlesurfaces.Thedataobtainedforaparticleradiusof8:2xwasusedtodeterminethelubricationcutohN()foreachviscosity,andthenumericalvaluesarerecordedinTable 3{1 .Theselubricationcorrectionsbringthesimulatednormalforceintoagreementwiththeoryforalltheparticlesizes,interparticlegaps,anduidviscositiesstudied(opensymbolsinFig. 3{1 ).Thecorrespondingresultfortheforceandtorque Table3{1: Lubricationrangesforvariouskinematicviscosities,determinedforasphereofradiusa=8:2x.Therangesaredeterminedseparatelyforthenormal,hN,tangen-tial,hT,androtational,hR,motions. =1=60.670.500.43=1=1000.240.500.15=1=12000.100.500.00

PAGE 59

Figure3{3: TorqueonaparticlesettlingnexttoaverticalplanarsurfaceatzeroReynoldsnumber.SimulationdatafortorquedividedbyT0=8a2hyU,isindicatedbysolidsymbols.Resultsincludingthelubricationcorrectionareshownbytheopensymbols. onasphereslidingalongthewallisshowninFigs. 3{2 & 3{3 .Againweseethatthelubricationcorrectiongivesconsistentlyaccurateforcesandtorques,thoughnotquiteasaccurateasthenormalforce.Thelubricationrangesfortangentialmotionwerefoundtobeindependentoftheuidviscosity,asshowninTable 3{1 .Wealsonoticedthatthereciprocalrelationsareobeyed;theforceonarotatingsphereissimilartothedatainFig. 3{3 .Thecalculatedtorqueonarotatingsphere(Fig. 3{4 )isinagreementwiththeoryforthehigherviscosities=1=6and=1=100,butnotatthelowestviscosity=1=1200.Herethelattice-Boltzmannmethodoverpredictsthetorqueonarotatingsphereby20-30%.Wethinkthattheerroriscausedbythelargedierencebetweenthehydrodynamicandinputradii(=0:55x)anditimpliesthatviscositieslessthan0.01arenotsuitableforsuspensionsimulations,atleastwithbounce-backboundaryconditions.Inpracticethisisnotaseriouslimitation:aviscosity=0:01

PAGE 60

Figure3{4: TorqueonaparticlerotatingnexttoaverticalplanarsurfaceatzeroReynoldsnumber.SimulationdatafortorquedividedbyT0=8a3hy,isindicatedbysolidsymbols.Resultsincludingthelubricationcorrectionareshownbytheopensymbols. allowssimulationswithaReynoldsnumberupto10pergridpoint(withaMachnum-ber0:1),whichisatorbeyondthelimitofresolutionoftheow.Inotherwords,thereislittlepracticalvalueinviscositieslessthan0.01. Finally,weexaminedthedynamicmotionofaparticle(a=4:8x)settlingontoasolidwall(Fig. 3{5 ).ThelubricationforcewascalculatedusingtherangesgiveninTable 3{1 .Theparticlewasgivenanitemassandplacedwithaninitialgapofh=0:2ahybetweentheparticleandwall.ThesimulationswereperformedataReynoldsnumberRe0:02byapplyingaconstantforcetotheparticle.Theresultsshowtheexpectedexponentialdecayofthegapbetweentheparticleandwallovertime,inquantitativeagreementwithlubricationtheory( 59 )(shownbythesolidline).

PAGE 61

Figure3{5: Settlingofasphere(a=4:8)ontoahorizontalwall.Thegapbetweentheparticlesurfaceandwall,h,relativetothehydrodynamicradius,isplottedasafunc-tionofthenon-dimensionaltime(opencircles). forcesandtorques(Eq. 2.97 ),followedbythelubricationforces.Theoverallprocedureisstillrstorderaccurate,butthelubricationforcescancauseinstabilitieswhenevertheparticlesareinnearcontact.Theinstabilityisdrivenbythenormalforces,andthestabilitycriteria m=6a2t 3sa3h=9 2t sah<2 (3.9) isviolatedwhenhislessthan0:1x. Itisimpracticaltosolvealltheequationsimplicitly,soweimplementedanalgorithmwhichusesanimplicitupdateonlywherenecessary.Inoursimulationsweusedtheconservativecriteriat=m<0:1.Schematically,wesolvethecoupleddierentialequations _x=Ax+b bysplittingthedissipativematrixAintoregularandsingularcomponents,A=AR+AS.InourcontextASonlycontainscomponentsofthenormalfrictioncoecient

PAGE 62

Figure3{6: Illustrationofthealgorithmtodeterminethelistofclusters. whenthegapbetweenparticlesislessthanthestabilitycuto,hs,determinedfromEq. 3.9 .ThusARcontainsallthenon-zerocomponentsofthelubricationcorrectionbutwiththeinterparticleseparationinthenormalforceregularizedbyhssothatthelargerofhijandhsisusedtocalculatetheforce.TheremainingnormalforceisincludedinAS,withtheformofEq. 3.7 ,butwithhNreplacedbyhs.Usingamixedexplicit-implicitdierencing, t=ARx(t)ASx(t+t)+b;(3.11) weobtaintherstorderupdate (1+ASt)x(t+t)=x(t)ARtx(t)+bt:(3.12) Theimportantpointisthat,byasuitablerelabelingoftheparticleindices,AScanbecastintoablockdiagonalform,withthepotentialforanenormousreductioninthecomputationtimeforthematrixinversion.Therelabelingisachievedbya

PAGE 63

Figure3{7: Themaximumclustersizeasafunctionoftheclustercutogaphs=aatvaryingvolumefractions,. clusteranalysis.First,allpairsofparticlesthatarecloserthanthestabilitycutoareidentied,andalistismadeofallsuchpairs.AnillustrationisshowninFig. 3{6 a,wherepairsofparticleswithseparationslessthanhsareindicatedbythesolidlines.Theclusterlabelsareinitializedtotheparticleindex;eachparticleisthenrelabeledbygivingitthesmallestlabelofalltheparticlesitisconnectedto.After1passthelabelsareasshowninFig. 3{6 bandafter2passes3distinctclustershavebeenidentied,eachwithauniquelabel(Fig. 3{6 c).Thealgorithmstopswhennofurtherrelabelingtakesplace.Althoughmoreecientschemesarepossible,thissimpleschemeismorethanadequateforourpurposes.Oncetheclustershavebeenidentied,theimplicitequationscanbesolvedforeachcluster.WeuseconjugategradientstoexploitthesparsenessofAS,whichisextensiveevenwithineachdiagonalblock. Figure3{8: Numberofclustersasafunctionoftheclustercutogaphs=aatvaryingvolumefractions,.

PAGE 64

Thecomputationalcostofthecluster-implicitalgorithmdependsprimarilyonthemaximumclustersize,whichisshowninFig. 3{7 asafunctionoftheclustercutogaphs.Arandomdistributionof1000particleswassampledinaperiodicboxatvolumefractions0.10,0.25,and0.50.Atlowandmoderatevolumefractionstheclustersizeisonlyweaklydependentonhs=a,rangingfrom2-7,andclustersofthissizeimposeanegligiblecomputationaloverhead.However,athighvolumefractionsthereisapercolationthresholdaths=a0:02beyondwhichasingleclustermoreorlessspansthewholevolume.Inthiscasetheclusterwillgrowtoencompassalmostalltheparticlesinthesystem.Thusathighdensitiescomputationaleciencyrequiresthaths=a<0:02.Whencombinedwiththestabilitycriteria,whichimplieshs=a1=a2,wendaminimumradiusofa=10xtokeept=m<0:5.Asimulationofseveralhundredsuchparticlesispossibleonapersonalcomputerordesktopworkstation. InFig. 3{8 weshowthecorrespondingnumberofclusters.Ingeneralthereisasteepriseinthenumberofclusterswithincreasinghs=a,levelingotoaround100clusters.Thesharpdropinthenumberofclustersatthehighestvolumefractionisassociatedwiththepercolationtransition,asseeninFig. 3{7 Inclusionofthelubricationforcesleadstolargeforcesandstidierentialequa-tionsfortheparticlevelocities.Wehavedevelopedamixedexplicit-implicitvelocity

PAGE 65

update,whichminimizesthenumberoflinearequationsthatmustbesolved,whilemaintainingabsolutestability.

PAGE 66

1 ).However,twodierentsetsofexperimentsfoundthatthevelocityuctuationsconvergetoaxedvalueforsucientlylargesystems( 60 61 ).Thequalitativediscrepancybetweentheoryandexperimenthasgeneratedconsiderableattention,focusingonpossiblemechanismsforscreeningthehydrodynamicinteractions( 62 )andtherebysaturatingthevelocityuctuationswhenthesystemsizeislargerthanthehydrodynamicscreeninglength. Thekeytheoreticalidea( 62 )isthathydrodynamicinteractionscanbescreenedbychangesinsuspensionmicrostructure,analogoustothescreeningofelectrostaticinter-actionsinchargedsystems.Hydrodynamicscreeningoccurswhenatestparticleanditsneighborsare,collectively,neutrallybuoyantwithrespecttothebulksuspension;inotherwords,whendensityuctuationsatlengthscaleslargerthanthescreeninglengtharesuppressed.Thisrequiresthataratherdelicatelong-rangecorrelationdevelopinthedistributionofparticlepairs,whichisresistanttotherandomizingeectsofhydro-dynamicdispersion.Severaldierentbulkmechanismsforthemicrostructuralchanges 54

PAGE 67

havebeenproposed,includingthree-bodyhydrodynamicinteractions( 62 ),aconvectiveinstability( 63 ),andacoupledconvection-diusionmodel( 2 64 ).However,thesethe-oriesareallinconsistentwiththeresultsofnumericalsimulations( 42 56 65 ),whichhaveshownthat,inhomogeneoussuspensions(withperiodicboundaryconditions),particlesaredistributedrandomlyatseparationsbeyondafewparticlediametersandthevelocityuctuationsremainproportionaltocontainersize. Morerecently,theideathatthecontainerboundariesplayacriticalroleindeterminingtheamplitudeofthevelocityuctuationshasbeenexploredinde-tail( 46 66 67 ).Theprimarygoalhasbeentodiscoverifcontainerwallscanquali-tativelychangethehydrodynamicinteractionsinthebulksuspension.Brenner( 66 )analyzedtheeectsofverticalcontainerwallsontheparticlevelocityuctuations,butfoundthatitdoesnoteliminatethedependenceonsystemsize.However,numericalsimulations( 46 )subsequentlyshowedthatrigidboundariesatthetopandbottomofthevesselcauseastrongtime-dependentdampingofthevelocityuctuationseveninbulkregionsfarfromthewalls.Theseboundariesintroduceinterfacesbetweensedi-ment,suspension,andsupernatentuid,whichareabsentinhomogenoussystemswithperiodicboundaryconditions.Thetime-dependentdampingofthevelocityuctuationswasexplainedbyconvectivedrainingoflarge-scaleuctuationsinparticledensitytotheseinterfaces( 46 ),asuggestionmadeearlierbyHinch( 68 ).Adierentpictureoftheeectofcontainerwallswasproposedindependently,andmoreorlesssimultaneouslybyTee( 67 ),whereitwassuggestedthathydrodynamicdispersionatthesuspension-supernatentinterfaceleadstoastraticationintheparticleconcentration.Astratiedsuspensionintroducesanotherlengthscale,beyondwhichthehydrodynamicinterac-tionsarescreened( 69 ).However,amodelbasedonconvectionofdensityuctuationsleadstoqualitativelydierentconclusionsfromastratiedsuspension,andwewillcomparetheseideastotheresultsofoursimulationsinSecs. 4.3.1 and 5.1.1 Recentexperimentalresultshavecastdoubtontheconclusionthatthevelocityuctuationsarenecessarilyindependentofcontainersize( 67 70 ).Brenner( 70 )foundthatthevelocityuctuationsinverydilutesuspensions(volumefraction<1%)didnotsaturate,evenforcontainerdimensionsaslargeas800a,whereaistheparticle

PAGE 68

radius.Teeetal:( 67 )foundthatforlargecellsanddilutesuspensionstheycouldnotevenobtainasteadystate.Insteadthevelocityuctuationsdecayedforthedurationoftheexperiment,nomatterwhattheheightofthecontainer;theseobservationscouldbeexplainedbyanincreasingstraticationofthesuspensionwithtime.Ontheotherhandsomeexperimentsfoundthatthevelocityuctuationsweresteadyandindependentofsystemsize,eveninverylargevessels( 71 ).Tomakefurtherprogress,itwillbenecessarytoreconciletheapparentlyconictingexperimentalresults,andtodevelopthecorrectphysicalpictureunderlyingthescreeningofthehydrodynamicinteractionsinsettlingsuspensions. Thefocusofthepresentworkisonthemicrostructureofasettlingsuspension,ascharacterizedbythedistributionofpairsofparticlesinthebulksuspension.AtlowReynoldsnumbertheinstantaneousparticlevelocitiesarecompletelydeterminedbytheparticlepositions,anditcanbeshownthatthedominantcontributiontothevelocityuctuationscanbecalculatedfromthestructurefactor whichistheFouriertransformofthepaircorrelationfunction( 64 65 ).Theeectsofhydrodynamicscreeningshowupinthelong-wavelength(smallk)behavioroftheS(k);anasymptotick2dependenceofthestructurefactorindicatesscreeningandaneventualsaturationofthevelocityuctuationswithincreasingcontainersize.IthasnotprovenpossibletomeasureS(k)directlybylightscattering,duetothelargesizeoftheparticles,butdirectimaginghasbeenusedtocalculaterelatedspatialuctuationsinparticleconcentration( 72 ).Inoursimulationswehavebeenabletodeterminethestructurefactordirectly.Herewepresentadetailedaccountofourresultsforboththevelocityuctuationsandthemicrostructure,includingeectsofpolydispersity,whichmaybeimportantininterpretingexperimentalresults.

PAGE 69

compromisebetweenseveralcompetingfactors,limitedbytherequirementthatacalculationcompleteinareasonableamountoftime,oftheorderofonemonth.Experi-ments( 61 )suggeststhatthecharacteristiclengthscaleisthemeaninterparticlespacingl=a1=3andthatthescreeninglengthisoftheorderof10l.Thecomputationaltimerequiredforoursimulationsis,toarstapproximation,proportionaltothenumberofuidgridpoints,andthereforethetotalvolumeisanimportantlimitation.Wehavelimitedourcalculationstoasinglevolumefraction=0:13(l2a),whichwaschosenasareasonablecompromisebetweenkeepinglsmallandavoidingtheadditionalcomplicationsofahighlyconcentratedsuspension.Weusedacellwithasquarecrosssection,andstudiedarangeofwidthsfromW=8ltoW=24l.InmostinstancesthecellswereboundedatthetopandbottombyrigidimpermeablewallsandwehavefoundfromexperiencethatacellheightH=1000aisnecessarytoallowtimeforthesuspensiontoreachasteadystate.Wehavevariedthecellheightinsomeinstances,asreportedinSec. 4.2.2 .Atthechosenvolumefractionthesuspensionscontainsbetween8000and72000solidparticles. Akeycomponentofthisworkistoassesstheeectsofthemacroscopicboundaryconditionsonthesuspensionmicrostructureanddynamics.Ano-slipboundarybreaksthemacroscopictranslationalinvarianceinthedirectionnormaltotheboundarysurface,andbecauseofthelong-rangecharacterofhydrodynamicinteractions,thissymmetrybreakingcanhaveasignicanteectfarfromtheboundary.Wehavethereforecomparedresultsforthreedierentsetsofboundaryconditions:systemsboundedinalldirectionsbyrigidwalls,whichwedenoteas\Box",systemsboundedbywallsatthetopandbottom,whiletheverticalfacesareperiodic(\Bounded"),andsystemswithperiodicboundaryconditionsonallfaces(\PBC").Forsystemswithperiodicboundariesonallfaces,theheighttowidthratiowassetto4:1,aspreviousworkshowednochangeoccurredwithtallercells( 56 ).Inthesimulationswithno-slipboundariesatthetopandbottom(\Box"and\Bounded"),datawastakeninawindow,typicallyofheight200a,centeredaround400afromthebottomofthevessel,asistypicalinexperimentalmeasurements.Wehavecheckedthattheresultsareinsensitivetothepositionoftheviewingwindowaslongasitislocated

PAGE 70

inaregionofuniformparticleconcentration.Fullyperiodicsystems(\PBC")arespatiallyhomogeneous,anddataisaveragedovertheentirevolumeinthiscase.Whenconsideringthesystemdimensions,itshouldbenotedthatanextra2ahasbeenaddedtoalldimensionsboundedbyno-slipwalls.Thisistoallowfortheexcludedvolumewiththewallsandkeeptheparticleconcentrationascloseaspossibleto13%.Wewillidentifythesystemdimensionswithoutanyexcludedvolumecorrection. Aparticleradiusoftwogridspacings,a=2x,wasusedinallthesesimulations,ascomparedwitha=1:5xora=1:25xinpreviouswork( 46 56 ).Testcalculationswithsmallnumberofparticles(seeSec. 4.2.3 )suggestthattherearenomajorchangesintheresultsiflargerparticlesareused.Thelargestsystemsstudiedinthisworkcontains72000particlesandapproximately20milliongridpoints.Thesimulationswererunfor2000Stokestimes,wheretheStokestimetS=a=U0isthetimeittakesanisolatedparticletosettleoneradius.ThegravitationalforcewassetsothattheStokestimecorrespondedto250timestepsandthecompletesettlingsimulationthereforerequiredhalfamillionsteps.Thecalculationtakesabout1monthonaclusterof8processors.ThetypicalparticleReynoldsnumber,basedonthemeansettlingvelocity,wasRe=2hUia==0:06,whereasinlaboratoryexperimentsitisusuallyonetotwoordersofmagnitudesmaller.Althoughtheresidualinertiaisapotentialsourceofcomplication,testsdescribedinSec. 4.2.3 andlaboratoryexperiments( 73 )suggesteectsofinertiaarenegligibleatReynoldsnumbersRe<1. 74 75 ).TherotationoftheparticlesgivesrisetoaliftforceatniteReynoldsnumber,whichtendstodrivetheparticlesapart.AtverysmallReynoldsnumbers,Re<0:1,therepulsionbetweentheparticlesbecomenegligibleandtheseparationdistancedoesnotchangeappreciablyovertime( 76 ).WeexaminethesettlingoftwoparticlesatarangeofparticleReynoldsnumbers0:0025
PAGE 71

Figure4{1: SettlingofahorizontalpairofparticlesatdierentReynoldsnumbers;Re0:1(circles),0.05(squares),0.025(triangles). Fig. 4{2 showstheseparationdistancer=abetweenparticlecentersastheysettlewithdierentinitialorientationstothehorizontal.Particlesorientedinthehorizontal(=0)andvertical(=90)directionsmaintaintheirinitialorientationsovertime,althoughtheparticlesdriftawayfromeachother(horizontal)ortowardseachother(vertical),duetotheresidualuidinertia.AtlowReynoldsnumbers,apairofexactlysphericalparticlesorientedatsomeintermediateangleshouldsettleasaunit,maintainingtheseparationandorientationofthepair.Ourresultsshowthatparticlesato-axisorientationshaveslightlydierentvelocitiessothepairrotateswithtime,evenatverylowRe.Thisisanumericalartifactcausedbythediscretelattice.Itcanonlybereducedbyusinglargerparticleswithmoregridpointstorepresenttheparticlesurface,orbyamorecomplicatedboundaryconditionthatsmoothsouttheraggedrepresentationofthesurfaceshowninFig. 2{2 .Neverthelessthegridartifactscausemuchsmallervariationsinparticlevelocitythanthepolydispersityinherentinlaboratoryexperiments.

PAGE 72

Figure4{2: Settlingofapairparticlesatvariousorientations:=0(circles),45(squares),68(triangles),90(diamonds).TheReynoldsnumberRe=0:1ineachcase. examinedfromseveraldierentperspectives;herewepresentamorecompleteandup-to-dateaccountofourowninvestigations. Simulationshaveshownthatthereisaqualitativedierenceinthetimeevolutionofthevelocityuctuationsdependingonthemacroscopicboundaryconditions( 46 ).ThevelocityuctuationsshowninFig. 4{3 illustratethedierenttimedependenceof\PBC"and\Box"geometries.Withperiodicboundaryconditionsthevelocityuctuationsinitiallyincreasewithtime,reachingaplateauafterapproximately200tS.

PAGE 73

Figure4{3: ParticlevelocityuctuationsasafunctionoftimewithawithW=48a. Thishaspreviouslybeenexplainedintermsofanincreasingnumberofparticlecontactsasthesettlingproceeds( 56 65 ).Thisisanexampleofachangeinsuspensionmicrostructure,i.e.howpairsofparticlesaredistributed,eventhoughtheparticleconcentrationisexactlythesame.However,whenthecontainerisboundedbyno-slipwalls,theuctuationsdecaywithtime,overaperiodofapproximately1000tS.Theinitialuctuationsarecomparabletothevaluesinthehomogeneoussuspension(\PBC"),butinthe\Box"geometrytheydecayovertimetoamuchsmallervalue.Asimilardecayinthevelocityuctuationshasbeenobservedexperimentally( 70 71 ),overcomparabletimescales. Simulationsalsodemonstratethattheboundaryconditionsontheverticalfacesofthecontainerdonotplayaqualitativelyimportantroleindeterminingthemagnitudeofthevelocityuctuations.ThedataillustratedinFig. 4{4 showthatvelocityuctua-tionsdecayinasimilarwayforboth\Box"and\Bounded"geometries,andthereforetheamplitudeofthevelocityuctuationsislargelycontrolledbytheboundarycon-ditionsatthetopandbottomofthecontainer( 46 ).Asanadditionaltest,weranafewsimulationswithno-slipverticalwallsandperiodicboundaryconditionsatthetopandbottom;theseresultsweresimilartothefullyperiodiccase.Verticalno-slipboundariesmayplayanimportantroleindeterminingtheuctuationsincontainers

PAGE 74

Figure4{4: Particlevelocityuctuationsasafunctionoftimefor`Box"and\Bounded"geometriesatvariouswidths:W=a=16(circles),32(squares),and48(triangles)respectively. withlargeaspectratiosinthehorizontalplane( 66 ),butourresultsshowthatthereisnoqualitativeeectincontainerswithasquarecross-section. Thedecayofthevelocityuctuationsinthebulkofthesuspensioncannotbeexplainedbydirectscreeningofthelong-rangeinteractionsbytheno-slipboundaries.Althoughano-slipboundarydoesscreenthehydrodynamicinteractionsofparticlesneartheboundary,thescreeningdoesnotextendtoparticlesthatarefartherfromthewallthantheyaretoeachother( 66 ).Thusinthemeasurementwindow,whichistypically5Wfromthebottom,thehydrodynamicinteractionsarenotsignicantlyaectedbythecontainerwalls.Moreover,thetime-dependentdecayoftheuctuationssuggeststhatthereisarearrangementofthetypicalparticlecongurationsduringthesettlingprocess.Themostlikelymechanismisthatthecontainerbottomandthesuspension-supernatentinterfaceactassinksofuctuationenergy,assuggestedearlierbyHinch( 68 ).Ahorizontaldensityuctuationisidealizedastworegions(orblobs)sidebyside,oneslightlyheavierthantheaverageandtheotherslightlylighter.Thesedensityuctuationsconvecttooneofthesetwointerfacesandareabsorbedbythedensitygradientattheinterface.Scalingargumentssuggestthatahorizontaldensity

PAGE 75

uctuationoflengthlconvectswithavelocityvlU0(l=a)1=2withrespecttothemean( 68 ),sothatuctuationsdrainawayonatimescale Intheabsenceofno-slipboundariesatthetopandbottomofthecontainer,large-scaledensityuctuationsrecirculatethroughthesuspension,andtheuctuationsindensityandvelocityarethereforetimeindependent. Toaccountforthevelocityuctuationsreachingsteadystate,themodelmustincludesomemechanismforreplenishingthelarge-scaledensityuctuations;otherwisethevelocityuctuationswillcontinuallydecrease.Wewillassumethatsmall-scale(ofordera)densityuctuationsaregeneratedbyconversionofgravitationalpotentialenergy,andthenspreadoutbyhydrodynamicdiusionresultingfromshort-rangemulti-particleinteractions.Thecharacteristiclengthofthedensityuctuationwillgrowtoorderlinatimeoforder whereDU0aisthehydrodynamicdispersioncoecient.Bybalancingtheconvectiontimetcwiththediusiontimetdweobtainacriticallengthscale,lca1=3,beyondwhichuctuationsdrainawaymorerapidlythantheycanbereplenishedbydiusion.Thusthesystemcanreachasteadystatewithacorrelationlengththatisindependentofsystemwidthandproportionaltothemeaninterparticlespacinga1=3,asobservedexperimentally( 61 71 ).Thisisthephysicalideaunderlyingtheconvection-diusionmodelproposedbyLevineetal:( 64 ),whichwillbediscussedindetailinSec. 4.3.2 .Analternativeexplanation,basedonthedevelopmentofastratieddensityproleinthesuspension( 77 ),willbediscussedinSec. 4.3.1 Withinaviewingwindowfarfromthesedimentandsupernatentinterfacesthevelocityuctuationsreachaquasi-steadystateafteratimeperiodoftheorderof1000tS.Thedurationofthesteadystate,priortothesupernatentinterfacereachingtheviewingwindow,dependsonthesystemheight,asshowninFig. 4{5 .However,themeansettlingvelocityandthevelocityuctuationswithintheviewingwindowareotherwiseunaectedbysystemheight.Fortheheightusedinmostofoursimulations

PAGE 76

Figure4{5: Meansettlingvelocity,Uk,anductuationsinsettlingvelocity,DU2kE,asafunctionoftimefordierentsystemheights,H. (H=1000a)wetypicallyhaveanuppertimelimittothesimulationsofabout1500tS,beforethesupernatentinterfacestartstointerferewiththemeasurementsintheviewingwindow.Wendthatthetimetakentoreachsteadystatedependsonthecontainerwidth(Fig. 4{4 ),rangingfrom600700tSforthesmallestsystem(W=16a)toabout1000tSforthelargestsystem(W=48a). whereListhecontainerdimension( 1 ).Preciseresultshavebeenworkedoutinthelow-concentrationlimitforseveraldierentgeometries( 65 77 78 ).ThedatainFig. 4{4 showsthattheinitialvelocityuctuationsfollowtheCaischandLuke( 1 )scaling

PAGE 77

Figure4{6: Particlevelocityuctuationsatsteadystate.Prolesareshownforthe\Box"system(solidsymbols)atdierentwidths:W=a=16(circles),32(squares),and48(triangles). (Eq. 4.4 ),whilethesteadystatevelocityuctuationsgrowmoreslowlywithsystemsize. Table4{1: Steadystateparticlevelocityuctuationsinamonodispersesuspensionasafunctionofsystemwidth. ThesteadystateuctuationsacrosstheviewingwindowareshowninFig. 4{6 .Thetimewindowforaveragingwasdeterminedseparatelyforeachsystem,dependingonthetimetoreachaquasi-steadystate(Fig. 4{4 );thedurationofthetimewindowwas400tSinallcases.Theresultswerealsoaveragedoverarangeofverticalposition,chosensothatthedensityandvelocityuctuationswereconstantapartfromstatisticalerrors(Fig. 5{4 ).ThespatialandtemporalwindowsusedinFig. 4{6 arelistedinTable 4{1 .Forthe\Bounded"geometrythevelocityuctuationswereaveraged

PAGE 78

overthewholehorizontalplane,whileforthe\Box"geometrytheuctuationsweremeasuredinathinsliceofwidth2alocatedinthecenterofthecontainer.Figure 4{6 showsthateventhoughthevelocityuctuationsgrowlessthanlinearlywithwidth,theydonotconvergetoawidth-independentvaluefortherangeofsystemsizesstudied. Figure4{7: Steadystateparticlevelocityuctuationsasafunctionofsystemwidthforsimulationversusexperimentalt(solidline). ThesteadystatevelocityuctuationsaresummarizedinTable 4{1 .Resultsforthe\Box"geometryareaveragedacrossthehorizontalplaneandalongthecenterlineofthesystem(thecentralpositioninFig. 4{6 );resultsforthe\Bounded"geometryareshownaveragedoverthehorizontalplane.Theresultsforthe\Box"geometryshowmoretendencytosaturationthanthe\Bounded"geometry,butitisnotclearthatthisisaqualitativeasopposedtoaquantitativedierence.Itwilltakelargersystemsizesthanarecomputationallyfeasibleatpresenttoresolvethisissue.However,Fig. 4{7 showsthatthesimulationdataagreesratherwellwithexperimentalmeasurementsforcomparablecontainersizes.Thesolidlinesarets( 61 )toavarietyofexperiments( 61 79 80 ).Theexperimentaldatasuggeststhatsomewhatlargersystemsizesareneededforthevelocityuctuationstosaturateandbecomeindependentofsystemsize.Wenotethatinclusionofpolydispersityimprovestheagreementwithexperiment,butwe

PAGE 79

willseeinSec. 5.1.3 thatadierentscreeningmechanismisoperatinginpolydispersesuspensions.ItshouldbenotedthatinusingtheanalyticttoexperimentaldatainFig. 4{7 ,wehavetakenthesmallestcontainerdimensionasthecontrollinglength,whereasSegre( 61 )usedthelargerlateraldimensioninttingtheirexperiments.Theyclaimedabettertfromthischoice,butitcannotbejustiedphysicallyortheoretically. Eectofinertiaonthemeansettlingvelocityanductuationsinsettlingvelocity. ThesimulationsdescribedinthispaperarenotexpectedtobeacompletelyquantitativedescriptionofStokes-owhydrodynamics.Inparticularly,wecannotdisregardthepossibilitythatinertialeectsplayaroleonscaleslargerthantheparticleradius.Atthevolumefractionusedinthiswork,wecanexpectthatparticlevelocitycorrelationspersistfordistancesoftheorderof40a( 61 ),andtheReynoldsnumberbasedonthisdistanceis1.2.NeverthelessthedatainFig. 4{8 showsthat,inasystemofwidthW=16a(N=8000),inertiaplaysnoroleforparticleReynoldsnumbersRe<0:1.Wehaveveriedthatthisconclusionholdsinalargersystemsaswell(W=32a,N=32000);nosignicantdierenceinvelocityuctuationswasobservedbetweenRe=0:06and0:03.Againcomputationallimitationspreventastudy

PAGE 80

Figure4{9: Eectofgridresolutiononthemeansettlingvelocityanductuationsinsettlingvelocity. oftheeectsofinertiainlargersystems,butexperimentalresults( 73 )suggestthatitissmallwheneverRe<1. Animportantconcernwithlattice-Boltzmannsimulationsareartifactsintroducedbythemotionofparticlesacrossthegrid.ItwasshowninSec. 4.1.2 thatthesegridartifactscauseaweakbutmeasurabledispersioninthetrajectoriesofpairsofparticle,andthesamerandomnoisemayservetoweakenthemicrostructuralchangesthatleadtosuppressionofvelocityuctuations.However,simulationswithlargerparticles,whichhavesmoothertrajectories( 81 ),arenotnoticeabledierent.TheeectsofgridresolutionareshowninFig. 4{9 ;thelargerparticlescorrespondtoamorerenedmesh.Therearequantitativedierencesinthevelocityuctuationscalculatedwitha=1:25x( 46 ),andthepresentresults(a=2x)overthetimerange250
PAGE 81

Figure4{10: Eectofcompressibilityonthespatialvariationofthemeansettlingvelocity.ThedataisforReynoldsnumberofRe=0:06. Undertheconditionsofthesimulations,withakinematicviscosity=0:167x2=t,aReynoldsnumberRe0:06correspondstoaMachnumberMa=0:004.Byreducingtheviscosity,wecanreducetheMachnumber,whilekeepingtheReynoldsnumberconstant.Figure 4{10 showsthatatconstantReynoldsnumber,Re0:06,themeanvelocityisindependentofpositionatsucientlysmallMachnumbers.However,wecouldnotdetectsignicantdierencesintheuctuationsateithersmallerReorsmallerMa,althoughthenon-uniformsettlingvelocitydoesleadtoasmallreductioninparticleconcentration(0:01)duringthecourseofthesimulation. 64 65 ),andinparticularitslow-klimitingbehavior: aZS(k) whereS(k)wasdenedinEq. 4.1 andisanumericalfactoroforder1.Ithasnotbeenpossibletomeasurethestructurefactorofasettlingsuspensionofnon-Brownianparticlesexperimentally,sincetheparticlesaretoolargeforlightscattering

PAGE 82

measurements.Fluctuationsinparticleconcentrationhavebeenmeasuredwithinacylindricalorrectangularwindow( 72 ),butthisonlygivesanangleaverageofthepairdistribution.InnumericalsimulationsitispossibletocalculateS(k)asafunctionofbothwavelengthanddirection.Thisisimportantsincesometheoriespredictthatthestructurefactorbecomeshighlyanisotropicatlong-wavelengths( 64 ),withthehorizontaluctuationsvanishingask2whiletheverticaluctuationsremainniteatallwavelengths.Herewereportonthestructurefactorinasteadilysettlingsuspension,andinvestigatethepossibilityofstraticationoftheparticleconcentration( 77 ). Particlevolumefractionasafunctionofheight,z.Thedotsindicatetheinstantaneousaverage,andthesteady-state(1000
PAGE 83

writtenas dt+U0@z=D@2z;(4.6)U0=dU=disthevelocityofadensityperturbationwithrespecttothemeansettlingspeedandDisthehydrodynamicdiusioncoecient.ThesolutiontoEq. 4.6 isqualitativelysimilartotheprolesshowninFig. 4{11 ;inparticularthereisabulkregionwheretheconcentrationisconstant.However,ifthediusioncoecientisverylarge,thesuspension-supernatentinterfacespreadsintothebulkleadingtoaweaklystratiedsuspension,withasmalldensitygradienteveninthebulk( 77 ). Straticationhasbeensuggestedasamechanismforsuppressingvelocityuctu-ations( 69 ),bymakingitpossiblefoructuationsinparticleconcentrationtoreachaneutrallybuoyantpositioninthesuspensionwithoutdrainingtotheinterfaces.How-ever,thedensityprolesshowninFig. 4{11 arenotstratied;insteadtheconcentrationproleisuniforminthebulkandthesuspension-supernatentinterfaceissharp.Theplotsofaverageconcentrationshowthatanymeanvariationindensityiswellbelowthestatisticalnoise,evenafteraveragingoverseveralhundredStokestimes.Ourdatasuggeststhatanyresidualdensitygradientmusthaveacharacteristiclengthofatleast104a,or10timesthecontainerheight. Hydrodynamicdispersiondoescauseaspreadingoftheinterface,butthisiscom-pensatedbyhinderedsettling,whichconvectsthelessdenseregionsatahighervelocitythanthehighdensityregionsandtherebysetsupaconvectiveuxinoppositiontothediusiveux.Balancingtheconvectiveanddiusiveconcentrationuxeswithrespecttoaframemovingwiththemeansettlingvelocity,weestimateaninterfacethickness,Dk=U05a,atthisvolumefraction.Heretheverticaldispersioncoecient,Dk=4Uka,wastakenfromoursimulationsandiscomparabletoexperimentalmeasurements( 80 ).Oursimulationsdonotruleoutthepossiblesignicanceofin-terfacialdiusioninverydilutesuspensions(<1%),suchasarecommonlyusedinParticle-Image-Velocimetrymeasurements( 67 70 71 ),butwedoexcludeitasageneralmechanismforhydrodynamicscreening.

PAGE 84

Figure4{12: Structurefactordescribinghorizontaldensityuctuations,S(k?),forvarioussystemsizes:W=16a,32a,and48a. 4{13 .Figure 4{12 showsthat,incomparisonwiththeinitialequilibriumdistribution,horizontaldensityuctuationsarestronglysuppressedbyno-slipboundariesatthetopandbottomofthecontainer.Ontheotherhanditisalreadyknownthatsuspensionswithperiodicboundaryconditionsdonotundergosignicantchangesinmicrostructureduringsettling( 56 83 ).Inparticularthestructurefactorremainsniteatallwavelengths.Thesignicanceofthisresultisthatitdemonstratesthatthemacroscopicboundarieshaveaprofoundeectonthedistributionofparticlesinthebulksuspension.Itmayalsobethemechanismbywhichhydrodynamicinterac-tionsarescreenedduringthesettlingprocess.AphysicalinterpretationofthedatainFig. 4{12 isthat,iftheviewingwindowwasdividedintoverticalslicesofsucientsize,

PAGE 85

therewouldbeessentiallythesamenumberofparticlesineachslice,ratherthantheexpectedvariationoforderp Thestructurefactorinasettlingsuspensiondevelopsastronganisotropy,asshowninFig. 4{13 .Althoughthehorizontaldensityuctuationsareproportionaltok2?,theverticaluctuationstendtoanon-zeroconstantatsmallkk.Dampingofhorizontaldensityuctuations(atleastasfastask2?)istheminimumrequirementforhydrodynamicscreening,aswasshownbyLevineetal:( 64 ).Theyproposedthattherearetwoqualitativelydistinctnon-equilibriumphasesforsettlingsuspensions,anunscreenedphasecharacterizedbyarandommicrostructureandascreenedphasewherethehorizontaldensityuctuationsaredampedoutatlongwavelengths.Theyderivedanexpressionforthenon-equilibriumstructure-factor, whichisconsistentinfunctionalformwiththestructurefactorobtainedinournumer-icalsimulations(Fig 4{13 ).Therenormalizedparameters,theuctuationsinparticleuxNiandthediusioncoecientsDi,togetherwiththedampingcoecient,,werecalculatedfromcoupledeldequationsdescribingtheevolutionoftheparticleconcen-trationanduidvelocity.Accordingtothetheory,thephaseboundaryisdeterminedbytheanisotropyintherenormalizednoise,N?=Nk,anddiusivity,D?=Dk. ThestructurefactordatacanbeusedtoextractratiosoftheparametersthatappearinEq. 4.7 ;namelyN?==0:4a2,andNk=Dk=0:17.WeobtainedN?=andNk=Dkfromthelow-kbehaviorofthehorizontalandverticaldensityuctuations,butsinceourdataisrathernoisy,itisimpossibletoextractameaningfulvalueofD?=,whichappearsasaquarticcorrectiontotheasymptotick2dependenceofthestructurefactor.Nevertheless,forthesakeofcompletenesswewilluseourbestestimate,D?==0:5a2,todeterminetheratioN?=Nk0:7.Whencombinedwithtracer-diusionmeasurementsofD?=Dk=0:16,thissuggestswearenearthetransitionbetweenscreenedandunscreenedphases( 64 ).Unfortunatelyourdataisnotsucientlyprecisetoenableadenitiveconclusiontobedrawn.Signicantlylargersystemssizes

PAGE 86

Figure4{13: Structurefactoratdierentangles;vertical[0,0,1]direction(circles),45o[1,0,1]direction(squares),andhorizontal[1,0,0]direction(triangles). willbenecessaryforaquantitativecomparisonwiththepredictionsofthetheory,withgreatlyincreasedcomputationalrequirements. Densityuctuationshavebeenmeasuredatotherlow-indexdirections;forexampleFig. 4{13 alsoshowsdataforthe45o[1,0,1]direction.Unfortunatelythesystemisnotperiodicforanyk-vectorsthatlieoutsidethehorizontalplane,sothisdataisinherentlynoisier.Furthercomplicatingtheanalysisistheeectofrenormalization,whichproducesashoulderthatisbestseenforthehorizontal[1,0,0]directionatka=4.Forka>=4thestructurefactorissimilartotheequilibriumdistribution(Fig. 4{12 )buttheconvectiveeectsdescribedqualitativelyinSec. 4.2.1 producesignicantdampingwhenka<=4.Thedataforthe45omayalsobeshowingashoulderatasmallerk,ka=6,butawidercellisnecessarytoconrmthis.Thisisanimportantpoint,becausehydrodynamicscreeningrequiresthatthedensityuctuationsaredampedinalldirectionsexceptthevertical.Inpracticethismeanstheshoulderisexpectedtomovetosmallerandsmallerkasthedirectionrotatesfromhorizontaltovertical. AlthoughthestructurefactorsmeasuredinthesimulationsareconsistentwiththepredictionsofLevine( 64 ),theirtheorypostulatesthatalltheimportantdynamicsoccursinthebulk,independentofthemacroscopicboundaryconditions.Althoughthisisalogicalassumption,ournumericalsimulationshaveshownthatitisincorrect.Simulationswithperiodicboundaryconditions( 42 56 )donotexhibitanydamping

PAGE 87

ofthehorizontaldensityuctuations,aswouldbeexpectedifthemodelwerecorrectinallessentials.Instead,oursimulationssuggestthatthecontainerbottomandthesuspension-supernatentinterfaceactassinksofuctuationenergy,assuggestedearlierbyHinch( 68 ).Randomdensityuctuationsconvecttooneofthesetwointerfacesandareabsorbedbythedensitygradientattheinterface.ThedatashowninFig. 4{14 supportsthisconclusion,albeitnotconclusively.Hereweshowthestructurefactorintheviewingwindowduringitsevolutionfromanequilibriumstatetothesteadystate.Despitethelimitedtimeaveraging(atotalof4200tS=800tSforeachplot),theimplicationisthatthelongwavelengthuctuationsdecayfastest.Ifso,thisisevidencefortheimmediateconvectionoflarge-scaledensityuctuations( 68 ),ratherthantheestablishmentofadensitygradientbyhydrodynamicdiusion( 77 ). ThescreeningbyMucha( 77 )isgeneratedbytransientverticalvariationsinparticleconcentration,ratherthanbyasteady-statechangeinpaircorrelations.ThekeydierenceisthatthetheorybyLevine( 64 )addsstronglyanisotropicconcentrationuctuations,whichareanempiricalrepresentationofthesupposedeectsofthemany-bodyhydrodynamicinteractions.Toobtainahydrodynamicallyscreenedphase,thetheoryrequiresthatsmall-scaleconcentrationuctuationsarelargestinthehorizontalplane,N0?N0k( 64 ).Bycontrast,inarandomsuspensionthedensityuctuationsareisotropiconalllengthscales.Oursimulationsshowthatthereisapronouncedchange Figure4{14: Timeevolutionofstructurefactorforhorizontaldensityuctuations.

PAGE 88

intheanisotropyofthedensityuctuationsasthesettlingproceeds,sothecharacterofthenoisemaychangeasthesuspensionevolvestosteadystate. 4{8 .Thisisagenericresultwheneverthereisano-slipwallboundingthetopandbottomofthecontainer,andisindependentofsystemheight(Fig. 4{5 ),Reynoldsnumber(Fig. 4{8 )andgridresolution(Fig. 4{9 ).Toourknowledgeasystematicvariationofsettlingvelocitywithtimehasnotbeenreportedpreviously,butitisconsistentwithatime-dependentreductioninnumber-densityuctuations,whichhasbeenobservedexperimentallybyLeietal:( 72 ).Inthatworkitwasshownthatthenumber-densityuctuationsinaxedvolumedecreasesasthesuspensionsettles.Wehaveobservedasimilardecreaseinnumber-densityuctuationsinoursimulations,butchangesinmicrostructureshowupmoreclearlyinthestructurefactor(Sec. 4.3.2 ). Ifthechangeinsettlingvelocityisduetochangesinthemicrostructureatlongwavelengths,thenitcanbeestimatedfromtherelationbetweenthesettlingvelocityandstructurefactor( 65 ),whichisquitepreciseatlongwavelengths: Thesteady-statestructurefactorwastakenfromEq. 4.7 withtheparametersdeter-minedfromthesimulationdata(Sec. 4.3.2 ),andtheisotropicequilibriumstructurefactorwasttedwithaquadraticpolynomial;theintegrationwastakenuptoamaxi-mumka==2,beyondwhichitwasassumedthatthestructurefactorsweresimilar.Thisgaveapredicteddecreaseinthemeansettlingspeed=0:20U0,whilethesimulationresultwasaround0:13U0.

PAGE 89

inthehorizontalplanewhereacleark2dependencewasobserved.Simulationswithdierentmacroscopicboundaryconditionsonthecontainerwallshavedemonstratedthattheboundaryconditionsatthetopandbottomofthecontainerplayacrucialroleindeterminingthedistributionofparticlepairs.ThemeasuredstructurefactorisconsistentwiththekeyqualitativefeaturespredictedtheoreticallybyLevineetal:( 64 )usingarenormalizedconvection-diusionmodel.However,wenotethatthistheorycannotyetexplainwhythemacroscopicboundaryconditionsshouldplayacrucialrole. Oursimulationspredictthatthemeansettlingvelocitydecreasesduringthesettlingprocess.Whilethishasbeenshowntobeconsistentwiththeobservedchangesinsuspensionmicrostructurethereisnoexperimentalconrmationatthepresenttime. Wehavefoundnoevidenceforstraticationinmonodispersesuspensionsatmoderateconcentrations.Ourresultsshowthattheinterfaceissharpandthattheconcentrationinthebulkisveryuniform.ThesuggestionbyMucha( 77 )thatourobservations( 46 )couldbeexplainedintermsofstraticationisincorrect.Wecon-cludethatthereisamechanismformicrostructuralrearrangementinthebulkofamonodispersesuspensionduringsettling,whichleadstoasubstantialreductionintheamplitudeofthevelocityuctuationsfromthepredictionsforrandomlydistributedparticles.BasedonoursimulationsandtheexperimentalmeasurementsofBernard-Michel( 70 ),webelievethatthequestionofsaturationofthevelocityuctuationsinmonodispersesuspensionsisstillopen.Bothoursimulationsandtheseexperi-mentsshowthesteady-statevelocityuctuationsgrowinglessrapidlythanlinearlyincontainersize,W,butnotyetfullysaturating.

PAGE 90

67 71 80 )althoughsomeexperimentsusesignicantlymoremonodis-perseparticles( 61 70 ).Inauidizedbeditiswellknownthatparticlessegregateaccordingtosize,withadensersuspensionoflargerparticlesatthebottomandlessdensesuspensionoflighterparticlesatthetop.Thelargerparticleshaveaninherentlyhighersettlingvelocityandthereforematchtheuidizationvelocityatahighercon-centrationthanthesmallerparticles.Thusauidizedbedofpolydisperseparticlesisinevitablystratiedinbothconcentrationandvolumefraction.Wehavestudiedthesettlingofpolydispersesuspensionstondoutifthereissignicantsegregationduringthesettlingprocess,andwhateectsthatmayhaveonthevelocityuctuationsandmicrostructure. 5{1 showstheinstantaneousdensityprolesina10%polydispersesuspen-sionatsteadystate(t=1200tS).Incontrasttoamonodispersesuspension(Fig. 4{11 ),thesuspensionsupernatentinterfaceisconsiderablybroader,asshowninFig. 5{1 a,andthereisaweakunderlyingstraticationoftheparticlevolumefraction.ThevolumefractionproleisbrokendowninFig. 5{1 bintocontributionsfromthreedierentsizeranges.Itcanbeseenthatthereisconsiderablesegregationbythistimeandthesuspension-supernatentinterfaceisdominatedbythesmallestparticles.Smallparticlestendtodriftintotheupperregionofthefront,whilethelargeparticlessettlefasterandaredominantintheregionnearestthedensepack;mediumsizeparticlesaredis-tributedthroughoutthesuspension.WeemphasizethattheinterfacespreadingseeninFig. 5{1 aisaresultofdierentialconvectionofdierentparticlesizesratherthanhydrodynamicdiusion( 77 ),whichwehaveseendoesnotproduceabroadinterfaceat 78

PAGE 91

Figure5{1: Prolesoftheparticlevolumefractionsasafunctionofthesystemheight. thisvolumefraction.Evenatconcentrationslessthan1%,convectivespreadingduetopolydispersityisoftenmoreimportantthanhydrodynamicdiusion( 84 ). Oursimulations,Fig. 5{2 ,showthatthereisastraticationoftheparticlemassdensity(orvolumefraction)inpolydispersesuspensions,whichpersistsdeepintothebulk.Incomparisonwiththemonodispersesuspension,herethereisacleardecreaseoftheparticlemassdensity(orvolumefraction)withheightthatisvisibleoverthestatisticalnoise.Segregationandstraticationcanbeexplainedbyaone-dimensionalconvection-diusionequationanalogoustoEq. 4.6 : wheretheparticlesizeshavebeendividedintoanumberofdierentsizeranges,withvolumefractioniandsettlingvelocityUi=Ui;0f().Thesettlingvelocitywasconstructedfromtheusualapproximationthattakesthesettlingvelocityoftheisolatedparticleandahindrancefunction,f(),basedonthetotalvolumefraction=Pii.Intheabsenceofdiusion,D=0,wendaninterfacespreadingandsegregationthatisqualitativelysimilartothesimulationdata,conrmingthatpolydispersityisthedominantmechanismforstraticationatthisvolumefraction.Morecomplicatedexplanationsofstraticationbasedoninterfacediusion( 77 )orconvectionofdensityuctuations( 85 )havenotconsideredtheeectsofpolydispersity,whichmaypartlyorevenfullyaccountfortheexperimentalobservations(seeSec. 5.1.4 ).

PAGE 92

Figure5{2: Comparisonofparticlevolumefractionprolesformonodisperseandpoly-dispersesuspensionsatsteady-state,1000
PAGE 93

Figure5{3: Particlevelocityuctuationasafunctionoftimefora10%polydis-persesuspension.Thedataistakenforthreesystemwidths:W=a=16(circles),32(squares),48(triangles). Theprolesofthevelocityuctuationsasafunctionofheightarequalitativelydierentinmonodisperseandpolydispersesuspensions.Figure 5{4 showsthatthevelocityuctuationsinmonodispersesuspensionsareindependentofverticalpositioninthebulksuspension,whileinpolydispersesuspensionstheydecaywithincreasingheight.IthasbeenpointedoutbyMuchaetal:( 77 )thatastratiedsuspensionisexpectedtoleadtoareductionofthevelocityuctuationsnearthesuspension-supernatentinterface,becausethedensitygradientislargerthere(Fig. 5{1 a)andsoismoreeectiveindampingthevelocityuctuations.ThedecreasinguctuationsshowninFig. 5{4 isfurtherevidencethatthepolydispersesuspensionisstratied,whilethemonodispersesuspensionisnot.ThedatainFigs. 5{2 and 5{4 suggestthatthevelocityuctuationsarecontrolledbystraticationforz>400a.Thegradientinvolumefractionatthispoint(z=400a),a=d(log)=dz2:5104,isconsistentwithascalingargument( 77 )thatsuggestsvelocityuctuationsarecontrolledbystraticationwhen>critical,whichforoursystemisroughly2104a1.

PAGE 94

Figure5{4: Comparisonofparticlevelocityuctuationsformonodisperseand10%polydispersesuspensions. 5{5 a.Thissuggeststhatdierentmechanismsfordampingthevelocityuctuationsexist.Inmonodispersesuspensions,thedistributionofparticlepairsadjustsitselfduringsettlingsothatthedensityuctuationsaredampedask2atlongwavelengths.Thisleadstoatime-dependentdampingofthevelocityuctuationsandthepossibilityofasize-independentsaturation.Inpolydispersesuspensionsthemicrostructureisapparentlyrandomizedbythevaryingsettlingspeedsandthismechanismnolongerholds.Howeverinthiscasetheparticlevelocityuctuationsmaybedampedbystraticationduetosegregationofthedierentparticlesizes. ThestructurefactorsfordierentdegreesofpolydispersityarecomparedinFig. 5{5 b.Thesystemsaresmallerinthiscaseandsotherearefewerk-vectors.Thestructurefactorfor2%polydispersityissimilartothemonodispersecase,apparentlyvanishingask2?atlongwavelengths.Forhigherdegreesofpolydispersity,5%and10%,thestructurefactortendstoanon-zerovalueatlowk?.Itappearsthatonlyverysmalldegreesofpolydispersitycanbetoleratediflaboratoryexperimentsaretomimicthepropertiesofamonodispersesuspension.

PAGE 95

Figure5{5: Comparisonofstructurefactorsfordierentdegreesofpolydispersity:a)monodisperse(circles)and10%polydisperse(squares)suspensions(W=48a,N=72000);b)2%(circles),5%(squares),and10%(triangles)polydispersesuspensions(W=32a,N=32000). 79 80 )haveusedsuspensionswithsignicantpolydispersity,inexcessof5%.Inthiscasetheparticledensitybecomesstratiedoverthedurationoftheexperiment,Fig. 5{2 ,evendeepintothebulkwheremeasurementsaremade.Thedevelopmentoflong-rangecorrelations(Fig. 4{12 )arehinderedbyevensmallamountsofpolydispersity,sothatstraticationmaybethedominantmechanismforcontrollingthevelocityuctuationsintheseexperiments.Wehavesolvedtheconvection-diusionequation(Eq. 5.1 )at13%concentrationand10%polydispersity;wendastraticationa=3104atapointH=3ofthewayupthecontainerandatatimewhenthesuspensionfronthasfallenapproximatelyH=2.Thecalculatedstraticationiscomparabletowhatwasfoundinourpolydispersesimulations(H=1000a)andissucientlystrongtocontrolthevelocityuctuations(Fig. 5{4 ). Manyexperiments( 61 67 70 71 84 )aredoneatverylowparticleconcentra-tions,topermittheuseofParticle-Image-Velocimetrymeasurements.Inthiscase

PAGE 96

hinderedsettlingisnegligibleandtheinterfacespreadsbybothdiusionandcon-vectivesegregationofdierentparticlesizes.Itisnotcomputationallyfeasibletosimulatesuspensionsattheselowconcentrationsbythelattice-Boltzmannmethod(Sec. 4.1.1 ),butitisinterestingtoconsiderwhetherstraticationcanbeimportantinthesesystems.Mucha( 77 )hasconsideredthisquestionindetailformonodispersesuspensions;herewecomparethecontributionsofsegregationanddispersion,viasolu-tionsofEq. 5.1 .AgainweconsiderthestraticationatalocationH=3afteratimeH=2U0asareferencepoint.Hinderedsettlingleadstothedevelopmentofatravelingconcentrationproleofconstantshape( 82 ),butindilutesuspensionstheshapeoftheproleisconstantlyevolving.Intheabsenceofhinderedsettling,thenaturallengthscaleistheheightofthecontainerandintheseunitsthereareonlytwoparameters;thedegreeofpolydispersity,,andthedimensionlessdiusioncoecient,D=D=U0H;thedimensionlessstraticationisnowH.Intheabsenceofpolydispersity(=0),largediusioncoecientsD>0:01areneededtoproducesignicantstratication.InouropinionthisisaweaknessoftheanalysisofMucha( 77 );inorderforinterfacediusiontoproducesignicantstratication,numericalvaluesofthediusioncoe-cientmustexceedallexperimentalmeasurements.However,evenmodestamountsofpolydispersitycanproducesignicantstratication.For=0:1,straticationisinfactdominatedbysegregation;areduceddiusioncoecientD=0:01isinsucienttocontributesignicantlytotheoverallstratication.Evenat5%polydispersity,segre-gationproducessignicantstraticationandcanbeassistedbyasmallhydrodynamicdispersion,D=0:001.At2%polydispersityorless,segregationisinsignicant,andtheconvection-diusionmodelsuggeststhatonlyinterfacediusioncanproducestratication. Oftheexperimentslistedatthebeginningofthepreviousparagraph,onlyoneset( 70 ),systematicallyexaminesthedependenceoncontainersizeusingnearlymonodisperseparticles,withpolydispersityaslowas1%.Intheseexperimentsaclearsaturationwithcontainersizewasnotobserved,particularlyatthelowestvolumefraction.Italsotooksignicantlylongerforthesuspensiontoreachasteady-statethaninotherexperiments.Thepolydispersitywasabout7%insomeexperiments( 67 71

PAGE 97

84 ),whichissucientforsignicantstraticationbysegregation.Segreetal:( 61 )usednearlymonodisperseparticles(=0:01),butprimarilyvariedthevolumefractionratherthanthecontainersize.Weconcludethattheissueofthesaturationofthevelocityuctuationsinastrictlymonodispersesuspensionisstillanopenquestion. Severalrecentexperimentshavebeencarriedoutatverylowparticleconcentra-tions,whileoursimulationsarelimitedbycomputationaleciencytomuchlargervolumefractions,inthiscase=0:13.Ourresultsdonotthereforediscountthepos-sibilitythatinterfacediusioncausessignicantstraticationatlowvolumefractionsasclaimedbyMucha( 77 ).Neverthelesstheconvection-diusionmodelshowsthatsegregationisoftenmoreimportantthaninterfacediusioninpromotingstratication.Onlyindilutesuspensionsofnearlymonodisperseparticleswouldweexpectinterfacediusiontobeimportant,althoughstraticationbysegregationcanoccurunderabroaderrangeofcircumstances.

PAGE 98

Wehaveappliedthelattice-Boltzmannmethodtoexamineparticledynamicsinagravitydrivensuspension.Initially,theparticlesexistinarandomlydistributedstatethatexhibitslargeparticlevelocityuctuationsandnolongrangepaircorrelationswithinthemicrostructure.Withperiodicboundariesapplied,theuctuationsarefoundtoincreaseovershorttimesandpersistthereafter.However,no-slipboundariesintro-ducedatthetopandbottomofthevesselleadtoadecayinthevelocityuctuationovertimewithinthebulk,whichiscomparablewithexperimentalndings.Here,thesteady-statemicrostructureexhibitsanorderedcongurationatlongwavelengths,whichsigniesthatscreeningispresentwhenno-slipboundariesareapplied. Thenumericalresultsareshowntocomparewellwithexperimentalresults,butduetocomputationallimitations,acompleteanalysisofthesedimentationproblemhasnotyetbeenattained.Inexperiments,thevelocityuctuationarefoundtobecomeindependentofsystemsizebeyondacertainpoint( 61 80 ).Ourresultsdonotunambiguouslyshowthiseect.Thevelocityuctuationsarefoundtogrowlessthanlinearlywithsystemwidth,buttheydonotconvergetoawidthindependentvalueforthesystemsizesevaluated.However,quantitativeresultsforthesystemsstudiedagreewellwithattoexperimentaldata( 61 ).Becauseexperimentswerecarriedoutwithlargercontainersbeforethescreeningeectbecamepronounced,simulationsoflargersystemsizesshouldbeevaluatedforamoreconclusivecomparison. Luke( 69 )proposedamechanismduetostraticationbyhydrodynamicdiusionatthesuspension-supernatentfrontleadingtoadampingofthevelocityuctuationswithinthebulk.Westudiedtheconcentrationproleforamonodispersesuspensionat13%volumefraction,butfoundasharpfrontdevelopattheinterface.Here,hinderedsettling,becauseofparticlecrowdinginterferingwiththemotionoftheindividualparticles,hasastrongerneteectthanadiusivespreading.Theconcentrationwasfoundthentobeuniformlydistributedwithinthebulk.Thisiscontraryto 86

PAGE 99

simulations( 77 )andexperiments( 61 71 )whichwereperformedatlowvolumefractions,wherestraticationfromdiusivespreadingoccursmorereadily.Asweseeit,thestraticationeectdoessuppressvelocityuctuations,butonlyatlowvolumefractionswherehinderedsettlingisnotpresent. Withoutstratication,themechanismfordensityuctuationdecayisthenattributedtoconvectivetransport.Wehypothesizethatdensityuctuationsinthebulkcandrainawayatthesuspension-sedimentandsuspension-supernatentinterfaces.Thisgeneratesacontinualdecayinvelocityuctuations,whichisreplenishedbyhydrodynamicdiusionduetoshortrangemultiparticleinteractions.Abalancebetweenconvectiveanddiusivetransportofdensityuctuationsleadstoacorrelationlengthscalebeyondwhichthevelocityuctuationsarescreened. Furtherexaminationofthesuspensionmicrostructureshowedthatanevolutiontowardsastatisticallynon-randomstatedeveloped,whereparticlesaredistributeduniformlyatlargelengthscales.Weidentiedthisresultbasedonthedampingofthelongrangepaircorrelationsinthestructurefactor,S(k?)!k2?ask?goestozero.Moreover,thepaircorrelationsarefoundtobeanisotropicbetweentheverticalandhorizontalplanes,suchthattheverticaluctuationswereniteatallwavelengths.Thiscorrelateswithaconvection-diusionmodel( 64 ),verifyingthatthetransportofuctuationsduetoconvectionisthedominantmechanismpresentleadingtothedecayinparticlevelocityuctuations. Monodispersesuspensionsrepresentanidealizedsituationwherethereisnovaria-tionsinparticlesize.Typicalexperimentsatmoderatetohighvolumefractions( 79 80 )exhibitsizepolydispersityof5%orgreater.Wendthatforpolydispersityinexcessof5%,adierentmechanismforthedecayinvelocityuctuationsisobserved.Thevariationinparticlesizeleadstoasegregationofdierentsizesasthesuspensionset-tled.Smallerparticlestendtosettleslowerandstayclosertothesupernatentfront,whilelargeparticlessettlefasteraccumulatingmoreintothebulk.Theneteectwasapronouncedstraticationofparticleconcentrationatthesuspension-supernatentfrontthatpersistsintothebulk.Becauseofthestratication,thevelocityuctuationswithinthebulkdonothavetoconvecttowardsthedensitygradientsatthefront,ratherthe

PAGE 100

gradientisstrongenoughinsidethebulktocreateasinkforthevelocityuctuationstodecay. Thestructurefactorinthepolydispersesuspensiondoesnotexhibitthesamedampingatlongwavelengthsseeninmonodispersesystems,indicatingamoreran-domlydistributedmicrostructure.However,velocityuctuationsarestilldamped,butthisisduetostraticationinducedbysegregationofdierentparticlesizes.Wesuggestthenthatpolydispersitymayhavealargereectinexperimentalmeasurementsthanpreviouslysupposed.

PAGE 101

[1] R.E.CaischandJ.H.C.Luke,\Varianceinthesedimentationspeedofasuspension,"Phys.Fluids,vol.28,pp.759,1985. [2] SriramRamaswamy,\Issuesinthestatisticalmechanicsofsteadysedimentation,"AdvancesinPhysics,vol.50,pp.297{341,2001. [3] P.J.MuchaandM.P.Brenner,\DiusivitiesandFrontPropagationinSedimen-tation,"Phys.ofFluids,vol.15,pp.1305,2003. [4] TheoG.M.VanDeVen,ColloidalHydrodynamics,AcademicPress,NewYork,1989. [5] W.B.Russel,D.A.Saville,andW.R.Schowalter,ColloidalDispersions,CambridgeUniversityPress,1989. [6] S.KimandR.T.Miin,\Theresistanceandmobilityfunctionsoftwoequalspheresinlow-Reynolds-numberow,"Phys.Fluids,vol.28,pp.2033,1985. [7] G.K.BatchelorandJ.T.Green,\Thehydrodynamicinteractionoftwosmallfreely-movingspheresinalinearoweld,"J.FluidMech.,vol.56,pp.401,1972. [8] G.K.Batchelor,\Sedimentationinadilutedispersionofspheres,"J.FluidMech.,vol.52,pp.245,1972. [9] J.HappelandH.Brenner,Low-ReynoldsNumberHydrodynamics,MartinusNijho,Dordrecht,1986. [10] SauroSucci,TheLatticeBoltzmannEquationforFluidDynamicsandBeyond,ClarendonPress,Oxford,2001. [11] KersonHuang,StatisticalMechanics,Wiley,NewYork,1987. [12] CarloCercignani,TheBoltzmannEquationandItsApplications,Springer-Verlag,NewYork,1975. [13] U.Frisch,B.Hasslacher,andY.Pomeau,\LatticegasautomatafortheNavier-Stokesequation,"Phys.Rev.Lett.,vol.56,pp.1505,1986. [14] A.J.C.LaddandR.Verberg,\Lattice-Boltzmannsimulationsofparticle-uidsuspensions,"J.Stat.Phys.,vol.104,pp.1191{1251,2001. [15] U.Frisch,D.d'Humieres,B.Hasslacher,P.Lallemand,Y.Pomeau,andJ-P.Rivet,\Latticegashydrodynamicsintwoandthreedimensions,"ComplexSystems,vol.1,pp.649,1987. [16] G.R.McNamaraandG.Zanetti,\UseoftheBoltzmannequationtosimulatelattice-gasautomata,"Phys.Rev.Lett.,vol.61,pp.2332,1988. 89

PAGE 102

[17] F.Higuera,S.Succi,andR.Benzi,\Latticegasdynamicswithenhancedcolli-sions,"Europhys.Lett.,vol.9,pp.345,1989. [18] Y.H.Qian,D.d'Humieres,andP.Lallemand,\LatticeBGKmodelsfortheNavier-Stokesequation,"Europhys.Lett.,vol.17,pp.479{484,1992. [19] A.J.C.Ladd,\NumericalsimulationsofparticulatesuspensionsviaadiscretizedBoltzmannequationPartI.Theoreticalfoundation,"J.FluidMech.,vol.271,pp.285,1994. [20] H.Chen,S.Chen,andW.H.Matthaeus,\Recoveryofthenavier-stokesequationsusingalattice-gasBoltzmannmethod,"Phys.Rev.A,vol.45,pp.R5339{5342,1992. [21] CarloCercignaniandGilbertoMedeirosKremer,TheRelativisticBoltzmannEquation:TheoryandApplications,BirkhauserVerlag,Boston,2002. [22] S.Y.Chen,D.Martinez,andR.W.Mei,\OnboundaryconditionsinlatticeBoltzmannmethods,"Phys.Fluids,vol.8,pp.2527{2536,1996. [23] O.FilippovaandD.Hanel,\Grid-renementforlattice-BGKmodels,"J.Comput.Phys.,vol.147,pp.219,1998. [24] R.S.Maier,R.S.Bernard,andD.W.Grunau,\BoundaryconditionsforthelatticeBoltzmannmethod,"Phys.Fluids,vol.8,pp.1788{1801,1996. [25] R.W.Mei,L.S.Luo,andW.Shyy,\Anaccuratecurvedboundarytreatmentinthelatticeboltzmannmethod,"J.Comput.Phys.,vol.155,pp.307{330,1999. [26] D.R.Noble,S.Y.Chen,J.G.Georgiadis,andR.O.Buckius,\Aconsistenthydrodynamicboundary-conditionforthelatticeBoltzmannmethod,"Phys.Fluids,vol.7,pp.203{209,1995. [27] P.A.Skordos,\InitialandboundaryconditionsforthelatticeBoltzmannmethod,"Phys.Rev.E,vol.48,pp.4823{4842,1993. [28] A.J.C.LaddandD.Frenkel,\Dynamicsofcolloidaldispersionsvialattice-gasmodelsofanincompressibleuid,"inCellularAutomataandModelingofComplexPhysicalSystems,P.Manneville,N.Boccara,G.Y.Vichniac,andR.Bidaux,Eds.,Berlin-Heidelberg,1989,number46inSpringerProceedingsinPhysics,pp.242{245,Springer-Verlag. [29] J.A.SomersandP.C.Rem,,"inShellConferenceonParallelComputing,G.A.vanderZee,Ed.1988,LectureNotesonComputerScience. [30] X.He,Q.Zou,L-S.Luo,andM.Dembo,\AnalyticsolutionsofsimpleowsandanalysisofnonslipboundaryconditionsforthelatticeBoltzmannBGKmodel,"J.Stat.Phys.,vol.87,pp.115{136,1997. [31] C.K.Aidun,Y.N.Lu,andE.Ding,\DirectanalysisofparticulatesuspensionswithinertiausingthediscreteBoltzmannequation,"J.FluidMech.,vol.373,pp.287{311,1998.

PAGE 103

[32] M.W.Heemels,M.H.J.Hagen,andC.P.Lowe,\Simulatingsolidcolloidalparticlesusingthelattice-Boltzmannequation,"J.Comput.Phys.,vol.164,pp.48{61,2000. [33] R.Cornubert,D.d'Humieres,andC.D.Levermore,\AKnudsenlayertheoryforlatticegases,"PhysicaD,vol.47,pp.241,1991. [34] A.J.C.Ladd,\NumericalsimulationsofparticulatesuspensionsviaadiscretizedBoltzmannequationPartII.Numericalresults,"J.FluidMech.,vol.271,pp.311,1994. [35] I.GinzbourgandP.M.Adler,\Boundaryconditionanalysisforthethree-dimensionallattice-Boltzmannmodel,"J.Phys.IIFrance,vol.4,pp.191,1994. [36] I.GinzbourgandD.d'Humieres,\Localsecond-orderboundarymethodsforlattice-Boltzmannmodels,"J.Stat.Phys.,vol.84,pp.927,1996. [37] H.D.Chen,C.Teixeira,andK.Molvig,\Realizationofuidboundaryconditionsviadiscreteboltzmanndynamics,"Int.J.Mod.Phys.C,vol.9,pp.1281{1292,1998. [38] H.Chen,\Volumetricformulationofthelattice-Boltzmannmethodforuiddynamics:Basicconcept,"Phys.Rev.E,vol.58,pp.3955{3963,1998. [39] M'hamedBouzidi,MouaouiaFirdaouss,andPierreLallemand,\MomentumtransferofaBoltzmann-latticeuidwithboundaries,"Phys.Fluids,vol.13,pp.3452,2001. [40] C.P.Lowe,D.Frenkel,andA.J.Masters,\Long-timetailsinangularmomentumcorrelations,"J.Chem.Phys.,vol.103,pp.1582{1587,1995. [41] A.J.C.Ladd,H.Gang,J.X.Zhu,andD.A.Weitz,\Time-dependentcollectivediusionofcolloidalparticles,"Phys.Rev.Lett.,vol.74,pp.318,1995. [42] A.J.C.Ladd,\Hydrodynamicscreeninginsedimentingsuspensionsofnon-Brownianspheres,"Phys.Rev.Lett.,vol.76,pp.1392,1996. [43] P.N.Segre,O.P.Behrend,andP.N.Pusey,\Short-timeBrownianmotionincolloidalsuspensions-Experimentandsimulation,"Phys.Rev.E,vol.52,pp.5070{5083,1995. [44] D.W.Qi,\LatticeBoltzmannsimulationsofparticlesinnonzeroReynoldsnumberows,"J.FluidMech.,vol.385,pp.41{62,1999. [45] P.Raiskinmaki,A.Shakib-Manesh,A.Koponen,A.Jasberg,M.Kataja,andJ.Timonen,\Simulationsofnon-sphericalparticlessuspendedinashearow,"Comput.Phys.Commun.,vol.129,pp.185,2000. [46] A.J.C.Ladd,\Eectsofcontainerwallsonthevelocityuctuationsofsediment-ingspheres,"Phys.Rev.Lett.,vol.88,pp.048301,2002. [47] J.F.Brady,\Rheologyofconcentratedcolloidaldispersions,"J.Chem.Phys.,vol.99,pp.567{581,1993.

PAGE 104

[48] J.R.MelroseandR.C.Ball,\Thepathologicalbehaviorofshearedhard-sphereswithhydrodynamicinteractions,"Europhys.Lett.,vol.32,pp.535{540,1995. [49] RenweiMeiandWeiShyy,\OntheFiniteDierence-BasedLatticeBoltzmannMethodinCurvilinearCoordinates,"J.Comp.Phys,vol.143,pp.426,1998. [50] C.K.Ghadder,\Onthepermeabilityofunidirectionalbrousmedia:Aparallelcomputationalapproach,"Phys.Fluids,vol.7,pp.2563,1995. [51] R.Glowinski,T.W.Pan,T.I.Hesla,D.D.Joseph,andJ.Periaux,\AFictitiousDomainApproachtotheDirectNumericalSimulationofIncompressibleViscousFlowpastMovingRigidBodies:ApplicationtoParticulateFlow,"J.Comput.Phys.,vol.169,pp.363,2001. [52] S.O.UnverdiandG.Tryggvason,\Afront-trackingmethodforviscous,incom-pressible,multi-uidows,"J.Comput.Phys.,vol.100,pp.25{37,1992. [53] J.F.BradyandG.Bossis,\Stokesiandynamics,"Ann.Rev.Fluid.Mech.,vol.20,pp.111,1988. [54] A.J.C.Ladd,\Hydrodynamictransportcoecientsofrandomdispersionsofhardspheres,"J.Chem.Phys.,vol.93,pp.3484,1990. [55] B.CichockiandB.U.Felderhof,\Short-timediusioncoecientsandhighfrequencyviscosityofdilutesuspensionsofsphericalBrownianparticles,"J.Chem.Phys.,vol.89,pp.1049,1988. [56] A.J.C.Ladd,\Sedimentationofhomogeneoussuspensionsofnon-Brownianspheres,"Phys.Fluids,vol.9,pp.491{499,1997. [57] I.L.ClaeysandJ.F.Brady,\Lubricationsingularitiesofthegrandresistancetensorfortwoarbitraryparticles,"PhysicoChem.Hydrodyn.,vol.11,pp.261,1989. [58] SangtaeKimandSeppoJ.Karrila,Microhydrodynamics:PrinciplesandSelectedApplications,Butterworth-Heinemann,Boston,MA,1991. [59] B.CichockiandR.B.Jones,\Imagerepresentationofasphericalparticlenearahardwall,"PhysicaA,vol.258,pp.273{302,1998. [60] H.NicolaiandE.Guazzelli,\Eectofthevesselsizeonthehydrodynamicdiusionofsedimentingspheres,"Phys.Fluids,vol.7,pp.3,1995. [61] P.N.Segre,E.Herbolzheimer,andP.M.Chaikin,\Long-rangecorrelationsinsedimentation,"Phys.Rev.Lett.,vol.79,pp.2574,1997. [62] D.L.KochandE.S.G.Shaqfeh,\Screeninginsedimentingsuspensions,"J.FluidMech.,vol.224,pp.275,1991. [63] P.TongandB.J.Ackerson,\AnalogiesbetweencolloidalsedimentationandturbulentconvectionathighPrandtlnumbers,"Phys.Rev.E,vol.58,pp.R6931{R6934,1998. [64] A.Levine,S.Ramaswamy,E.Frey,andR.Bruinsma,\Screenedandunscreenedphasesinsedimentingsuspensions,"Phys.Rev.Lett.,vol.81,pp.5944,1998.

PAGE 105

[65] A.J.C.Ladd,\Dynamicalsimulationsofsedimentingspheres,"Phys.FluidsA,vol.5,pp.299,1993. [66] M.P.Brenner,\Screeningmechanismsinsedimentation,"Phys.Fluids,vol.11,pp.754{772,1999. [67] S.Y.Tee,P.J.Mucha,L.Cipelletti,S.Manley,M.P.Brenner,P.N.Segre,andD.A.Weitz,\NonuniversalVelocityFluctuationsofSedimentingParticles,"Phys.Rev.Lett.,vol.89,pp.054501{1,2002. [68] E.J.Hinch,\Sedimentationofsmallparticles,"inDisorderandMixing,E.Guyon,Y.Pomeau,andJ.P.Nadal,Eds.,Dordrecht,1988,pp.153{161,KluwerAca-demic. [69] J.H.C.Luke,\Decayofvelocityuctuationsinastablystratiedsuspension,"Phys.Fluids.,vol.12,pp.1619{1621,2000. [70] G.Bernard-Michel,A.Monavon,D.Lhuillier,D.Abdo,andH.Simon,\ParticleVelocityFluctuationsandCorrelationLengthsinDiluteSedimentingSuspensions,"Phys.Fluids,vol.14,pp.2339{2349,2002. [71] E.Guazzelli,\Evolutionofparticle-velocitycorrelationsinsedimentation,"Phys.Fluids,vol.13,pp.1537{1540,2001. [72] X.Lei,B.J.Ackerson,andP.Tong,\SettlingStatisticsofHardSphereParticles,"Phys.Rev.Lett.,vol.86,pp.3300{3303,2001. [73] M.L.Cowan,J.H.Page,andD.A.Weitz,\Velocityuctuationsinuidizedsuspensionsprobedbyultrasoniccorrelationspectroscopy,"Phys.Rev.Lett.,p.453456,2000. [74] InchulKim,SaidElghobashi,andWilliamA.Sirignano,\Three-dimensionalowovertwospheresplacedsidebyside,"J.FluidMech.,vol.246,pp.465,1993. [75] K.O.L.F.Jayaweera,B.J.Mason,andG.W.Slack,\Behaviorofclustersofspheresfallinginaviscousuid,"J.FluidMech.,vol.20,pp.121,1964. [76] J.WuandR.Manasseh,\DynamicsofDual-ParticlesSettlingUnderGravity,"Int.J.MultiphaseFlow),vol.24,pp.1343{1358,1998. [77] PeterJ.Mucha,Shang-YouTee,DavidA.Weitz,BorisI.Shraiman,andMichaelP.Brenner,\AModelforVelocityFluctuationsinSedimentation,"J.FluidMech.,vol.501,pp.71{104,2004. [78] D.L.Koch,\Hydrodynamicdiusioninasuspensionofsedimentingpointparticleswithperiodicboundaryconditions,"Phys.FluidsA,vol.6,pp.2894,1994. [79] J.M.HamandG.M.Homsy,\Hinderedsettlingandhydrodynamicdispersioninquiescentsedimentingsuspensions,"Int.J.MultiphaseFlow,vol.14,pp.533,1988. [80] H.Nicolai,B.Herzhaft,E.J.Hinch,L.Oger,andE.Guazzelli,\Particleveloc-ityuctuationsandhydrodynamicself-diusionofsedimentingnon-Brownianspheres,"Phys.Fluids,vol.7,pp.12,1995.

PAGE 106

[81] N.-Q.NguyenandA.J.C.Ladd,\Lubricationcorrectionsforinlattice-Boltzmannsimulationsofparticlesuspensions,"Phys.Rev.E,vol.66,pp.046708,2002. [82] R.H.DavisandM.A.Hassen,\Spreadingoftheinterfaceatthetopofaslightlypolydispersesedimentingsuspension,"J.FluidMech.,vol.196,pp.107,1988. [83] N.-Q.NguyenandA.J.C.Ladd,\Microstructureinasettlingsuspensionofhardspheres,"Phys.Rev.E,vol.69,pp.050401(R),2004. [84] L.Bergougnoux,S.Ghicini,E.Guazzelli,andJ.Hinch,\SpreadingFrontsandFluctuationsinSedimentation,"Phys.Fluids,vol.15,pp.1875{1887,2002. [85] P.N.Segre,\OriginofStabilityinSedimentation,"Phys.Rev.Lett.,vol.89,pp.254503,2002.

PAGE 107

Nhan-QuyenNguyenwasbornintheBinhThuanprovinceofVietnamonAugust10,1977.HecametoAmericawhenhewas3yearsold,livinginNewOrleans,LA.untilheturned22.HeattendedTulaneUniversityinNewOrleans,wherehereceivedhisbachelor'sdegreeinchemicalengineeringandmathinMay1999.HewentontograduateschoolattheUniversityofFlorida,chemicalengineeringdepartment,wherehereceivedhisdoctoratedegreeinAugust2004. 95


Permanent Link: http://ufdc.ufl.edu/UFE0006486/00001

Material Information

Title: Sedimentation of hard sphere suspensions at low Reynolds number ssing a lattice-Boltzmann method
Physical Description: Book
Language: English
Creator: Nguyen, Nhan Quyen ( Dissertant )
Ladd, Anthony J. ( Thesis advisor )
Publisher: University of Florida
Place of Publication: Gainesville, Fla.
Publication Date: 2004
Copyright Date: 2004

Subjects

Subjects / Keywords: Chemical Engineering thesis, Ph.D   ( local )
Dissertations, Academic -- UF -- Chemical Engineering   ( local )

Notes

Abstract: Particle motion in a settling suspension is directly influenced by hydrodynamic interactions, such that the particle velocities will exhibit a random component, characterized 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 experiments. 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 independent 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 suspension-sediment and suspension-supernatent 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 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) --> k² as k goes to zero. Our numerical simulations show that the long-range correlations found in monodisperse 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.
Subject: Boltzmann, hydrodynamics, LBE, lubrication, sedimentation, suspensions
General Note: Title from title page of source document.
General Note: Document formatted into pages; contains 107 pages.
General Note: Includes vita.
Thesis: Thesis (Ph.D.)--University of Florida, 2004.
Bibliography: Includes bibliographical references.
General Note: Text (Electronic thesis) in PDF format.

Record Information

Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
System ID: UFE0006486:00001

Permanent Link: http://ufdc.ufl.edu/UFE0006486/00001

Material Information

Title: Sedimentation of hard sphere suspensions at low Reynolds number ssing a lattice-Boltzmann method
Physical Description: Book
Language: English
Creator: Nguyen, Nhan Quyen ( Dissertant )
Ladd, Anthony J. ( Thesis advisor )
Publisher: University of Florida
Place of Publication: Gainesville, Fla.
Publication Date: 2004
Copyright Date: 2004

Subjects

Subjects / Keywords: Chemical Engineering thesis, Ph.D   ( local )
Dissertations, Academic -- UF -- Chemical Engineering   ( local )

Notes

Abstract: Particle motion in a settling suspension is directly influenced by hydrodynamic interactions, such that the particle velocities will exhibit a random component, characterized 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 experiments. 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 independent 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 suspension-sediment and suspension-supernatent 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 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) --> k² as k goes to zero. Our numerical simulations show that the long-range correlations found in monodisperse 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.
Subject: Boltzmann, hydrodynamics, LBE, lubrication, sedimentation, suspensions
General Note: Title from title page of source document.
General Note: Document formatted into pages; contains 107 pages.
General Note: Includes vita.
Thesis: Thesis (Ph.D.)--University of Florida, 2004.
Bibliography: Includes bibliographical references.
General Note: Text (Electronic thesis) in PDF format.

Record Information

Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
System ID: UFE0006486:00001


This item has the following downloads:


Full Text










SEDIMENTATION OF HARD SPHERE SUSPENSIONS
AT LOW REYNOLDS NUMBER USING A LATTICE-BOLTZMANN METHOD



















By

NHAN-QUYEN 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

Nhan-Quyen 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 Long-Range Interactions and the Divergence Problem . 8
1.3.2 Multi-Particle Interactions . . . . . 10
1.3.3 Particles Near Contact . . . . . 11

2 DESCRIPTION OF THE LATTICE-BOLTZMANN 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 Lattice-Boltzmann Method ........ . . . 17
2.2.1 Gas Kinetics to Fluid Dynamics: The C!i inp ,i-Enskog Ex-
pansion . . . . . ... 22
2.2.2 Solid-Fluid Boundary Conditions ..... . . 30
2.2.3 Hydrodynamic Radius of a Particle .... . . 33
2.2.4 Particle M otion . . . . . . 37

3 INCORPORATING LUBRICATION INTO THE
LATTICE-BOLTZMANN 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 HARD-SPHERE 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 Time-Dependent Velocity Fluctuations . . ... 59
4.2.2 Steady-State 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

2-1 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

2-2 Hydrodynamic radius ahy (in units of Ax) for various fluid viscosities; a is
the input particle radius . . . . . . . 35

3-1 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

4-1 Steady state particle velocity fluctuations in a monodisperse suspension as
a function of system width . . . . . . 64

5-1 Steady-state particle velocity fluctuations in a 10W'. polydisperse suspen-
sion as a function of system width. . . . . . 79















LIST OF FIGURES
Page

2-1 The 18 velocity vectors (ci) in the lattice-Boltzmann model. . .... 17

2-2 Location of boundary nodes for a curved surface (a) and two surfaces in
near contact (b). . . . . . . . .... 29

2-3 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

2-4 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

2-5 Actual (solid lines) and hydrodynamic (dashed lines) surfaces for a parti-
cle settling onto a wall . . . . . . . 36

3-1 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

3-2 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

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 . . . . . . . 47

3-4 Torque on a particle rotating next to a vertical planar surface at zero Reynolds
number. Simulation data for torque divided by To-87/'~., Q, is indicated
by solid symbols. Results including the lubrication correction are shown
by the open symbols . . . . . . . 48

3-5 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 non-dimensional time (open circles). . . 49

3-6 Illustration of the algorithm to determine the list of clusters. . . 50

3-7 The maximum cluster size as a function of the cluster cutoff gap hs/a at
varying volume fractions, . . . . . . 51









3-8 Number of clusters as a function of the cluster cutoff gap hs/a at varying
volume fractions, . . . . . . . .... 51

4-1 Settling of a horizontal pair of particles at different Reynolds numbers;
Re ~ 0.1 (circles), 0.05 (squares), 0.025 (triangles). . . ..... 58

4-2 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

4-3 Particle velocity fluctuations as a function of time with a with W = 48a. 60

4-4 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

4-5 Mean settling velocity, (U\), and fluctuations in settling velocity, AU 2,
as a function of time for different system heights, H. . . . 63

4-6 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

4-7 Steady state particle velocity fluctuations as a function of system width
for simulation versus experimental fit (solid line). . . . 65

4-8 Effect of inertia on the mean settling velocity and fluctuations in settling
velocity . . . . . . . . . 66

4-9 Effect of grid resolution on the mean settling velocity and fluctuations in
settling velocity . . . . . . . . 67

4-10 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 steady-state (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

4-14 Time evolution of structure factor for horizontal density fluctuations. 74

5-1 Profiles of the particle volume fractions as a function of the system height. 78

5-2 Comparison of particle volume fraction profiles for monodisperse and poly-
disperse suspensions at steady-state, 1000 < t/ts < 1400. . . 79

5-3 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









5-4 Comparison of particle velocity fluctuations for monodisperse and 10'-.
polydisperse suspensions . . . . . . 81

5-5 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 LATTICE-BOLTZMANN METHOD

By

Nhan-Quyen 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 suspension-sediment and suspension-supernatent 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 long-range 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 non-Brownian particles within a viscous fluid. We

deal with low-Reynolds 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 >~ L-2, 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 lattice-Boltzmann simulations to verify

this hypothesis.

Simulations of the sedimentation problem require large-scale numerical calcula-

tions, for which the lattice-Boltzmann 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 large-scale systems spread over multiple processors can

be performed. With boundary conditions at the surface of the particles, the lattice-

Boltzmann method becomes applicable to multi-phase 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 non-equilibrium 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 coarse-grained 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 second-order 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 rate-of-strain tensor that describes fluid deformation, and A is the

vorticity tensor that corresponds to a solid-body 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 Navier-Stokes equations.

The Navier-Stokes 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 steady-state 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
Navier-Stokes equation reduces to the Stokes or creeping-flow equation. In dimensional

form, Eq. 1.13 can be re-written 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 Fourier-transform 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)

k-u 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 Long-Range 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

long-range hydrodynamic interactions, and causes the settling velocity of a suspension

to be less than that of an isolated particle.

1.3.2 Multi-Particle 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 multi-particle interactions, the above equations become
N
sF < -. Upu + C -. ] (1.52)
j-1
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)
j-1
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 LATTICE-BOLTZMANN MODEL

The lattice-Boltzmann 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 non-equilibrium 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 so-called 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 one-body distribution function (10) reads as

follows:

Dtf [= t + P- r + F. p] f(r, p,t) = C12, (2.1)

where the left-hand 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 (two-body) 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 two-body 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 three-body 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 short-range two-body 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 one-particle 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 2-1: The 18 velocity vectors (ci) in the lattice-Boltzmann 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 Lattice-Boltzmann Method

The computational utility of lattice-Boltzmann models depends on the fact that

only a small set of velocities are necessary to simulate the N i, r-Stokes 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. 2-1. 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 18-velocity 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 lattice-Boltzmann 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 one-particle 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 post-collision 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 lattice-gas 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 ill-suited 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 multi-speed model be applied, which in this case

is the 18-velocity 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 non-equilibrium part

i -nq + neq. (2.19)

The equilibrium distribution is a collisional invariant (i.e., A,(n-q) = 0). The form

of the equilibrium distribution is constrained by the moment conditions required to

reproduce the inviscid (Euler) equations with non-dissipative 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 19-velocity 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 non-equilibrium distribution,

as in the (C'!i -pn ,-Enskog solution of the Boltzmann equation. The fundamental

limitation of this class of lattice-Boltzmann 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 non-equilibrium 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 A-cA, (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 higher-order

eigenvectors of L that are not relevant to the N ',1i r-Stokes equations, but which do

affect the boundary conditions at the solid-fluid 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 non-hydrodynamic 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 solid-fluid boundaries, and thus the more flexible

collision operator described by Eq. 2.27 is preferred.

A 3-parameter 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

post-collision 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 non-equilibrium 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 post-collision

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 lattice-Boltzmann equation with the
multi-time-scale analysis applied (15, 14). The Navier-Stokes 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 lattice-Boltzmann 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 Navier-Stokes
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 Chapman-Enskog 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, r-Stokes equation. In this sense, the N ,',i, r-Stokes 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 Navier-Stokes 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 Navier-Stokes 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 Navier-Stokes equation.

While the C'!i Ipn :-Enskog expansion was originally applied to the continuous

Boltzmann equation in Eq. 2.1, here we derive the Navier-Stokes 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 Bhatnagar-Gross-Krook

(BGK) model where naq represents the equilibrium distribution and 7 is a typical

time-scale 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 r-la 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 -
i-O


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 -ra-a 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
a-u--- 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 ii-Enskog 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 Solid-Fluid Boundary Conditions

Boundary conditions in the lattice-Boltzmann model are straightforward to

implement, even for non-planar 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. 2-2).

The motion of a particle is determined by the forces and torques acting on it.

These forces consist of external, non-hydrodynamic forces, and hydrodynamic forces.

Non-hydrodynamic forces, such as attraction and/or repulsion between surfaces due

to an inter-particle 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 bounce-back

rule can be applied to surfaces of arbitrary shape, without additional complications.

However, for second-order convergence the bounce-back 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 bounce-back rules at the midpoints of links connecting the lattice nodes.

This is depicted by the arrows in Fig. 2-2a 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 L-2 (30) (L being the channel width);


u8/ue Ax2/L2,


(2.79)









Ub






0




t+ 01 -

Figure 2-3: 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 rigid-body 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. 2-3. 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 post-collision 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 lattice-Boltzmann 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 lattice-Boltzmann 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 lattice-Boltzmann 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 no-slip boundary condition leads

to a small mass transfer across a moving solid-fluid 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 link-bounce-back rule at the solid-fluid 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 second-order convergence can be obtained, even for dense suspensions (34).

The hydrodynamic radius can be determined from the drag on a fixed sphere in

a pressure-driven 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 2-4: 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. 2-4. 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 buov-int 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 finite-volume 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 bounce-back scheme, the fluctuations

in force can be reduced by an appropriate choice of particle radius. We have noticed







35

Table 2-1: 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.738e-03 1.208e-02 4.332e-04 5.674e-04


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 2-1 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 2-2.

The bounce-back rule leads to a velocity field in the region of the boundary

that is first-order 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. 2-5), 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 2-2); the dimensionless kinematic

viscosity v* = vAt/Ax2. For small particles (a < 5Ax), A also depends weakly on

the particle radius (Table 2-2). 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 second-order convergence from an inherently first-order

boundary condition; it will not be necessary when more accurate boundary conditions

are implemented.

The hydrodynamic radii (ahy) in Table 2-2 were determined from the drag force on

a single sphere in a periodic cubic cell (34). The Reynolds number in these simulations

Table 2-2: 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 2-5: 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 multi-particle suspension, all

distance calculations are based on the hydrodynamic radius (as shown in Fig. 2-5);

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 2-2 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 half-time 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 high-frequency 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

semi-implicit update corresponds to a = The explicit, implicit, and semi-implicit

updates evaluate the velocity-dependent force at t, t + At, and t + -At respectively. In

practice we find only small differences between semi-implicit and fully implicit methods

and we use the fully implicit method (a = 1) in this work. The boundary node map

is updated infrequently (every 10-100 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 semi-implicit 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
LATTICE-BOLTZMANN METHOD

3.1 Introduction

Lattice-Boltzmann 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 lattice-Boltzmann 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 pair-by-pair, 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 lattice-Boltzmann

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 "cluster-implicit" 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 conjugate-gradient 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 2-2b), 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 solid-fluid 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 long-term drift.

However, we enforce mass conservation, particle-by-particle, 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

lattice-Boltzmann simulation.

For particles close to contact, the lubrication force, torque, and stresslet can be

calculated from a sum of pairwise-additive contributions (53), and if we consider only

singular terms, they can be calculated from the particle velocities alone (57). In the

grand-resistance-matrix (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 lattice-Boltzmann 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. 3-1, 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 (al+a2)2 hhN (3.7)

= 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 cut-off 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 3-1: 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 2-2. The location of the

planar wall was shifted by A(v*), corresponding to the a -> o limit in Table 2-2

(Fig. 2-5). In this way we ensure that the lubricating surfaces are in the same position

as the hydrodynamic boundaries in the lattice-Boltzmann 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. 3-1 3-4 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. 3-1). 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 3-2: 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 lattice-Boltzmann 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 3-1. 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. 3-1). The corresponding result for the force and torque


Table 3-1: 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. 3-2 & 3-3. 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 lattice-Boltzmann 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 bounce-back

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 3-4: 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. 3-5). The lubrication force was calculated using the ranges

given in Table 3-1. 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 lattice-Boltzmann










10 0

10-

10-2

10-3 v*= 1/6

10 -4 .. ... ..
0.1 1 10 100 1000

10 0
10


10-2

10-3 V*= 1/100

10-4 .. .. .. .. .. ...
0.1 1 10 100 1000
vt/a2

Figure 3-5: 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 non-dimensional 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 3-6: 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 non-zero 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
explicit-implicit 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 3-7: 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. 3-6a,

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. 3-6b and after 2 passes 3 distinct clusters have been identified,

each with a unique label (Fig. 3-6c). 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 3-8: Number of clusters as a function of the cluster cutoff gap hs/a at varying
volume fractions, Q.









The computational cost of the cluster-implicit algorithm depends primarily on

the maximum cluster size, which is shown in Fig. 3-7 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 2-7, 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. 3-8 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. 3-7.

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 closely-spaced 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 explicit-implicit velocity







53

update, which minimizes the number of linear equations that must be solved, while

maintaining absolute stability.















CHAPTER 4
SEDIMENTATION OF HARD-SPHERE 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 long-range 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 three-body hydrodynamic interactions (62), a convective

instability (63), and a coupled convection-diffusion 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 time-dependent 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 time-dependent damping of the velocity fluctuations

was explained by convective draining of large-scale 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 long-wavelength (small k) behavior of the

S(k); an ..i-mptotic 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 = aQ-1/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 no-slip boundary breaks

the macroscopic translational invariance in the direction normal to the boundary

surface, and because of the long-range 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 no-slip 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 no-slip 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 side-by-side 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. 4-1 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 4-1: Settling of a horizontal pair of particles at different Reynolds numbers;
Re ~ 0.1 (circles), 0.05 (squares), 0.025 (triangles).


Fig. 4-2 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 off-axis 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. 2-2. 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 steady-state 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


60-A A A
40 U AA A A A A A
20

200 400 600 800
t/t
S

Figure 4-2: 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

up-to-date account of our own investigations.

4.2.1 Time-Dependent Velocity Fluctuations

The initial particle configuration is sampled from the most-probable (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 5-10 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 well-stirred distribution, without large-scale 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. 4-3 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 4-3: 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. 4-4 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 no-slip vertical walls and periodic boundary conditions at the

top and bottom; these results were similar to the fully periodic case. Vertical no-slip

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 4-4: 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 cross-section.

The decay of the velocity fluctuations in the bulk of the suspension cannot be

explained by direct screening of the long-range interactions by the no-slip boundaries.

Although a no-slip 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 time-dependent 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

suspension-supernatent 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 no-slip boundaries at the top and bottom of the container, large-scale

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 large-scale density fluctuations; otherwise

the velocity fluctuations will continually decrease. We will assume that small-scale (of

order a) density fluctuations are generated by conversion of gravitational potential

energy, and then spread out by hydrodynamic diffusion resulting from short-range

multi-particle 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 convection-diffusion

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 quasi-steady 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. 4-5. 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 4-5: 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. 4-4), ranging from 600 700ts for the smallest system (W 16a)

to about 1000ts for the largest system (W = 48a).

4.2.2 Steady-State Settling

A settling particle produces a disturbance to the fluid flow that decays as r-1 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 R-2, where R

is the interparticle separation. Assuming the particles are distributed uniformly, then

integration over volume leads to the well-known result:


(U2) x L, (4.4)


where L is the container dimension (1). Precise results have been worked out in the

low-concentration 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 4-6: 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 4-1: 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 600-1000 150-550
32 0.048 0.052 0.053 800-1200 200-500
48 0.059 0.063 0.076 1100-1500 250-450
< AUJ > U02 16 0.012 0.016 0.007 Bounded 700-1100 200-650
32 0.018 0.022 0.013 800-1200 250-600
48 0.022 0.028 0.018 1000-1400 250-500



The steady state fluctuations across the viewing window are shown in Fig. 4-6.

The time window for averaging was determined separately for each system, depending

on the time to reach a quasi-steady state (Fig. 4-4); 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. 5-4). The spatial and temporal windows used in Fig. 4-6 are listed in

Table 4-1. 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 4-6

shows that even though the velocity fluctuations grow less than linearly with width,

they do not converge to a width-independent 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 4-7: 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 4-1. Results for the

"Box" geometry are averaged across the horizontal plane and along the centerline of

the system (the central position in Fig. 4-6); 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. 4-7

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. 4-7, 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 4-8: 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 Stokes-flow 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. 4-8 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 4-9: 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 lattice-Boltzmann 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. 4-9; 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 lattice-Boltzmann 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 4-10: 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 4-10 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 non-uniform 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 low-k 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 long-wavelengths (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 4-11: Particle volume fraction as a function of height, z. The dots indicate the
instantaneous average, and the steady-state (1000 < t/ts < 1400) density profile in the
viewing window is shown in the .,.1i ,ent plots.


Figure 4-11 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 suspension-supernatent interface could lead to a weak stratification of

the particles and a non-uniform density in the bulk (3, 67, 77). The time evolution of

the concentration profile, Q(z, t), can be described by a convection-diffusion 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. 4-11; in particular there is a bulk

region where the concentration is constant. However, if the diffusion coefficient is very

large, the suspension-supernatent 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. 4-11 are not stratified; instead the concentration

profile is uniform in the bulk and the suspension-supernatent 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 Particle-Image-Velocimetry 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 4-12: 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. 4-13. Figure 4-12 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. 4-12 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. 4-13. Although the horizontal density fluctuations are proportional

to kI_, the vertical fluctuations tend to a non-zero 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 non-equilibrium 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 non-equilibrium structure-factor,


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 4-13). 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 low-k 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 ..i-mptotic 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

tracer-diffusion 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 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).

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 low-index directions; for example

Fig. 4-13 also shows data for the 450 [1,0,1] direction. Unfortunately the system is not

periodic for any k-vectors 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

suspension-supernatent 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. 4-14

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 large-scale 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 steady-state 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 small-scale 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 = 0-200
t/t = 200-400
A t/t 400-600


106-

50 4


Figure 4-14: 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. 4-8. This is a generic result whenever there

is a no-slip wall bounding the top and bottom of the container, and is independent of

system height (Fig. 4-5), Reynolds number (Fig. 4-8) and grid resolution (Fig. 4-9).

To our knowledge a systematic variation of settling velocity with time has not been

reported previously, but it is consistent with a time-dependent reduction in number-

density fluctuations, which has been observed experimentally by Lei et al.(72). In that

work it was shown that the number-density fluctuations in a fixed volume decreases

as the suspension settles. We have observed a similar decrease in number-density

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 steady-state 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 long-wavelength 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 convection-diffusion 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 steady-state 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 5-1 shows the instantaneous density profiles in a 10'(. polydisperse suspen-

sion at steady state (t = 1200ts). In contrast to a monodisperse suspension (Fig. 4-11),

the suspension supernatent interface is considerably broader, as shown in Fig. 5-la, 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

suspension-supernatent 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. 5-la 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 5-1: 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. 5-2, 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 one-dimensional

convection-diffusion 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 5-2: Comparison of particle volume fraction profiles for monodisperse and poly-
disperse suspensions at steady-state, 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 steady-state 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. 5-3. The

time dependence and steady-state values of the velocity fluctuations are comparable

to the monodisperse case (Fig. 4-7), 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 5-1 in the same format

as in Table 4-1.

Table 5-1: Steady-state 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 600-1000 150-400
32 0.083 0.094 0.095 800-1200 200-400
48 0.104 0.115 0.144 1100-1500 250-400
< AU_2 > U02 16 0.011 0.015 0.009 Bounded 800-1200 200-500
32 0.021 0.028 0.018 800-1200 250-500
48 0.029 0.037 0.029 1000-1400 200-400







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 5-4 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. 5-la) and so is

more effective in damping the velocity fluctuations. The decreasing fluctuations shown

in Fig. 5-4 is further evidence that the polydisperse suspension is stratified, while the

monodisperse suspension is not. The data in Figs. 5-2 and 5-4 -,--.-, -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 10-4, 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 10-4a-1.









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 5-4: 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. 5-5a. 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 size-independent

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 k-vectors. 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 non-zero 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 5-5: 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. 5-2, even

deep into the bulk where measurements are made. The development of long-range

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 convection-diffusion equation (Eq. 5.1) at

1t;: concentration and 1A. polydispersity; we find a stratification f3a 3 x 10-4 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 Particle-Image-Velocimetry 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 lattice-Boltzmann 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 convection-diffusion 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 steady-state 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 segregation-induced 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 convection-diffusion 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 lattice-Boltzmann 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, no-slip 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

steady-state microstructure exhibits an ordered configuration at long wavelengths,

which signifies that screening is present when no-slip 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 suspension-supernatent 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 suspension-sediment and suspension-supernatent 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 non-random 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 convection-diffusion 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 suspension-supernatent 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.