Investigation of granular batch sedimentation via DEM–CFD coupling

ORIGINALPAPER

InvestigationofgranularbatchsedimentationviaDEM–CFDcoupling

T.Zhao·G.T.Houlsby·S.Utili

Received:3February2014/Publishedonline:12November2014©Springer-VerlagBerlinHeidelberg2014

AbstractThispaperpresentsthreedimensionalnumeri-calinvestigationsofbatchsedimentationofsphericalparti-clesinwater,byanalysesperformedbythediscreteelementmethod(DEM)coupledwithcomputationalfluiddynam-ics(CFD).Byemployingthismodel,thefeaturesofbothmechanicalandhydraulicbehaviourofthefluid-solidmix-turesystemarecaptured.Firstly,theDEM–CFDmodelisvalidatedbythesimulationofthesedimentationofasinglesphericalparticle,forwhichananalyticalsolutionisavail-able.Thenumericalmodelcanreplicateaccuratelytheset-tlingbehaviourofparticlesaslongasthemeshsizeratio(Dmesh/d)andmodelsizeratio(W/Dmesh)arebothlargerthanagiventhreshold.Duringgranularbatchsedimentation,segregationofparticlesisobservedatdifferentlocationsinthemodel.Coarsegrainscontinuouslyaccumulateatthebot-tom,leavingthefinergrainsdepositedintheupperpartofthegranularassembly.Duringthisprocess,theexcessporewaterpressureinitiallyincreasesrapidlytoapeakvalue,andthendissipatesgraduallytozero.Meanwhile,thecompressibilityofthesedimentsdecreasesslowlyasasoillayerbuildsupatthebottom.Consolidationofthedepositedlayeriscausedbytheself-weightofgrains,whilethecompressibilityofthesampledecreasesprogressively.

ElectronicsupplementarymaterialTheonlineversionofthisarticle(doi:10.1007/s10035-014-0534-0)containssupplementarymaterial,whichisavailabletoauthorizedusers.

T.Zhao(B)·G.T.Houlsby

DepartmentofEngineeringScience,UniversityofOxford,OxfordOX13PJ,UK

e-mail:tzhaooxford@gmail.com

S.Utili

SchoolofEngineering,UniversityofWarwick,CoventryCV47AL,UK

KeywordsGranularbatchsedimentation·DEM–CFDcoupling·Segregation·Porewaterpressure·Effectivestress·Compressibility

1Introduction

Thesedimentationandconsolidationprocessesofgranu-larmaterialsarecommoninbothterrestrialandsubmergedenvironments.Inthenaturalenvironment,granularmate-rialsoftensettlecontinuouslytowardstheseabed,lakeorriverfloortoformaloosesedimentlayer.Astheskeletonofthesedimentlayerisextremelycompressible,itunder-goesrelativelylargestrainunderadditionalloads.Sedimen-tationprocessesareveryimportantforsolid-liquidsepa-ration,ascanbefoundinthefieldsofchemical,mining,wastewater,food,pharmaceuticalandotherindustries[1–3].Therelatedresearchincludestheoretical,experimentalandnumericalinvestigationsonavarietyofmaterials[3–10].

Duringgranularsedimentation,theaveragesettlingveloc-ityofasuspensionisthemostnotableparameterquantifyingthedynamicbehaviourofthefluid-solidmixturesystem[11–15].AsfirstproposedbyKynch[4],theaveragesettlingvelocityofasuspensiondependsonlyonthelocalconcen-trationofsolidmaterials(inadditiontothecharacteristicsoftheparticlesandfluid).Thisstatementhasbeenvalidatedbyexperimentalmeasurementsofthegrainsettlingratesingran-ularsuspensionsystems[6,8,16].Numericalsimulationsofsedimentationmainlyfocusonthesettlinganddepositionalbehaviourofparticlesineithermonodisperse[7,12,17]orpolydisperse[10,18]systems.Thisresearchrevealsthreedistinctzonesinagrainsettlingsystem,enumeratedfromtoptobottomas:1.thehinderedsettlingzone,wherethesurfacesettlingvelocityisapproximatelyconstant;2.the

123

922transitionzonewherethesettlingratedecreasesgraduallytowardszero;and3.thecompressionzonewhereasoillayerisformedatthebottomandconsolidationoccursduetotheself-weightofsediments.However,asreportedbyRichard-sonandZaki[7],zone1isabsentforveryfineflocculatedpulps,andthesettlingratedecreasesprogressively.

Formanyyears,analyticalandnumericalinvestiga-tionsofsedimentationemployingempiricalcorrelationsofthemixturepropertiesbasedonlaboratoryexperimentshavebeenreported[3,7,19–22].However,afullysystem-aticstudyofsedimentationisstilllackingasthefluid-solidmixturepresentsahighlyheterogeneousstructurefeaturingaspatiallynon-uniformdistributionofparti-cles[6,23].Thenumericalinvestigationreportedhereisbasedontheconceptthatthemotionofparticlesiscom-pletelygovernedbytheNewtonianequationsofmotion,andinterparticlecollisionsaremodelledbythesoftparti-cleapproach[3].ThefluidflowiscalculatedbytheNavier–Stokesequations[3,11,24].Consequently,thediscreteele-mentmethod(DEM)[25]andcomputationalfluiddynamics(CFD)[26]techniquescanbeusedtostudythemechanicalandhydraulicbehaviourofparticlesandfluidflow,respec-tively.Bycouplingthesetwomethods,acompleteanalysisofthesedimentationofafluid-solidmixturesystemcanbeachieved[27].

TheworkpresentedhereispartofaresearchefforttoextendtheDEMmodellingoflandslides[28–30]tothesub-mergedcase.Thepaperisorganisedasfollows.WeexplainthetheoryandmethodologyoftheDEMandCFDinSect.2.Thenumericalinvestigationofgranularsedimentationispre-sentedinSect.3.WefirstvalidatetheDEM–CFDcouplingmodelbysimulatingthesedimentationofasinglesphericalparticle(Sect.3.1).Then,weanalysethebatchsedimen-tationofparticles,examiningthesegregationofparticles,generationanddissipationofexcessporewaterpressures,andevolutionofcontactforcechains(Sect.3.2).Section4summarizestheresultsandmainconclusionsachievedbythework.

2Theoryandmethodology

Theequationsgoverningafluid-solidmixturesystemarederivedfromthetheoryofmultiphaseflow[27,31].Figure1illustratesthecomponentsofthemixture,whichconsistsoffluid,andparticles.Thefluiddensity(ρf),velocity(U)andparticlepackingporosity(n)arefunctionsofspatialposi-tionandtime.TheDEMandCFDopensourcecodesESyS-Particle[32,33]andOpenFOAM[34]wereemployedforthesimulationspresentedhere.ThecouplingalgorithmfromChenetal.[27]originallywritteninYADE[35]wasimple-mentedinESyS-Particlebytheauthors.

123

T.Zhaoetal.

Fig.1Aninfinitesimallysmallelementofthefluid-solidsystem

2.1Equationsgoverningparticlemotion

AccordingtoNewton’ssecondlawofmotion,theequationgoverningthetranslationalmotionofasphericalparticleisexpressedas:

md2󰀍󰀂i−dt2−→xi=mi−→g+→fnc+−→f󰀇tc+−f−f−luid→(1)

c

wheremiisthemassofparticlei;→−xiisthetroid;−→gisthegravitationalacceleration;−position→ofitscen-fncand−→ftcarethe

normalandtangentialparticle-particlecontactforcesexertedbytheneighbouringparticlesonparticlei,whicharecalcu-latedusingthelinear-springandrollingresistancemodelasdetailedbyIwashita[36],Jiangetal.[37]andBelheineetal.[38];thesummationofthecontactforcesisoverallparticlesincontactwithparticlei;−f−f−luid→

the

istheforceexertedbyfluidontheparticle,whichwillbedefinedinSect.2.2.TherotationalmotionofasphericalparticleisgovernedbyEq.(2),as:

Id󰀍idt−→ωi=

−→rc×−→ftc+−M→r(2)cwhere−→Iiisthemomentofinertiaabouttheωiistheangularvelocity;→−graincentroid;

rcisthevectorfromtheparticle

masscentretothecontactpoint;−M→

ristherollingresistantmoment,whichinhibitsparticlerotationoverotherparti-cles.TherollingresistantmomentusedintheDEMmodelaccountsapproximatelyfortheangularshapeandinterlock-ingeffectsbetweenparticles.2.2Fluid-particleinteraction

Theinteractionforcebetweenfluidandparticles(−f−f−luid→

)consistsoftwoparts:hydrostaticandhydrodynamicforces[39].Thehydrostaticforceaccountsforthefluidpressure

Investigationofgranularbatchsedimentationgradientaroundanindividualparticle(i.e.buoyancy)[27,40,41],expressedas:−→fbi=−vpi∇p

(3)

where−→

fiisthehydrostaticbuoyantforceactingonparticlei,vb

piisthevolumeofparticlei;pisthefluidpressure.Thehydrodynamicforcesactingonaparticlearethedrag,liftandvirtualmassforces.Thedragforceiscausedbytheviscousshearingeffectoffluidontheparticle;theliftforceiscausedbythehighfluidvelocitygradient-inducedpressuredifferenceonthesurfaceoftheparticleandthevirtualmassforceiscausedbyrelativeaccelerationbetweenparticleandfluid[42–44].ThelattertwoforcesarenormallyverysmallwhencomparedtothedragforceinsimulatingfluidflowatrelativelylowReynoldsnumbers[44].Therefore,theliftandvirtualmassforcesareneglectedinthecurrentDEM–CFDcouplingmodel.Inthisprocess,thedragforceoccurswhenthereisanon-zerorelativevelocitybetweenfluidandparti-cles.Itactsattheparticlecentreinadirectionoppositetotheparticlemotionrelativetothefluid[45].Experimentalcorre-lations[20,21,46]andnumericalsimulations[47–49]forthedragforcearereportedintheliterature.Inthisresearch,thedragforce(Fdi)actingonanindividualparticleiscalculatedusingtheempiricalcorrelationproposedbyDiFelice[22],as:

Fdi=1

πd22Cdρf

4

|U−V|(U−V)n−χ+1(4)

whereCdisthedragforcecoefficient;dandVareparticlediameterandvelocity,respectively.

Theporositycorrectionfunctionn−(χ+1)inEq.(4)rep-resentstheinfluenceofthepackingconcentrationofgrainsonthedragforce.Theexpressionforthetermχis[22]:󰀋󰀆χ=3.7−0.65exp−

1.5−log󰀁2󰀅

10Rep

2

(5)

whereRep=ρfdn|U−V|/μistheReynoldsnumberdefinedattheparticlesizelevel,withμbeingthefluidvis-cosity.Inthecurrentanalyses,χrangesfrom3.4to3.7.Thereareseveraldefinitionsofdragforcecoefficientreportedintheliterature[49–51].AcomparisonbetweenthesecorrelationsandexperimentaldataisshowninTable1andFig.2.AccordingtoFig.2,itcanbeobservedthatthecor-relationofBrownandLawler[51]matchestheexperimentaldatawellinthewholerangeofReynoldsnumbersconsid-eredinthisresearch.Therefore,thedragforcecoefficientisimplementedintheDEM–CFDcouplingmodelas:Cd=24󰀂0󰀇0.407Re1+0.150Re.681

+1+8710(6)

Re923

Table1Somecorrelationsofdragforcecoefficientforsphericalpar-ticlesStokes[42]Re24

Newton(Inertial)[3]0.44

󰀉24󰀆

SchillerandNaumann[50]1.0+0.15Re0.687󰀁Re≤1000DallaValle[46]󰀎

0Re.44Re0.63+√4.8

󰀃>1000

2ReBrownandLawler[47]

Re24

󰀆1+0.150Re0.681󰀁+0.4071+8710

Re1000010.810000.6)dC( tn1000.4eiciffeoC10103104105 Experimental data ga Brown and Lawler (2003)rD1DallaValle (1948)Schiller and Nauman (1935)Stokes0.1Newton10-410-2100102104106Reynolds Number (Re)Fig.2Dragcoefficientofvariouscorrelationsandexperimentaldataforsphericalparticles

Thetotalforceexertedbyfluidonasingleparticleisthereforeexpressedas:

−f−f−luid→

=−vpi∇p+1πd22Cdρf

4

|U−V|(U−V)n−χ+1(7)2.3Governingequationsoffluidflow

Thecontinuityandmomentumequationsgoverningthemotionoffluidflowinafluid-solidmixturesystemarederivedfromthetheoryofmultiphaseflow[31],as:∂󰀆ρfn

󰀁

∂󰀆

∂ρt+󰀆󰀁fnU󰀁∇·ρfnU=0(8)∂t+∇·󰀆ρfnUU󰀁=−n∇p+n∇·τ+nρf→−g+fd

(9)wherefdis󰀌theaveragedragforceperunitfluidvolume,definedaswithiniN=1Fdi/Vmesh,withNbeingthenumberofparticlesthefluidmeshcell,Vmeshisthevolumeofthefluidmeshcell;τisthefluidviscousstresstensor.

123

924T.Zhaoetal.

malizedsettlingtimeisdefinedas:

󰀆󰀁gd2ρs−ρftt

[T]==

tc18μH

(11)

Theverticalposition(h)isnormalizedbytheinitialparal-lelepipedheight(H),as:[H]=h/H

(12)

Fig.3Configurationofthegranularsedimentationmodel

3NumericalinvestigationofgrainsedimentationThisresearchexaminesthesettlingbehaviourofparticlesin

fluidwithinaparallelepiped,employingtheDEMcoupledwithCFD.Thefundamentalparametersgoverningtheset-tlingprocessarethefluiddensity(ρf)andviscosity(μ),thewidthoftheparallelepiped(W),thediameterofparticle(d)andtheporosityofthegranularpacking(n).BasedonthedimensionalanalysisperformedbyRichardsonandZaki[7],afunctionrelatingalltheseparametersis:

󰀈󰀊

WρdUrUr

,n,(10)=f

U0µdwhereUristherelativesettlingvelocitybetweenparticleandfluid,U0=gd2(ρs−ρ)/18μistheterminalsettlingvelocityofasinglesphericalparticleinfluid,calculatedbyStokes’lawofsedimentation[46].Intheanalysis,thenor-Table2Inputparametersofthesedimentationsimulation

InEq.(10),thefirsttwodimensionlessgroupsontherighthandsidecorrespondtotheReynoldsnumberandtheporos-itycorrectiontermsofthegoverningequationsintheDEM–CFDcouplingmodel(i.e.Eqs.(4)–(6)).Thesizeratio(W/d)representstheinfluenceofmodelsizeonthesettlingveloc-ityofparticles.Asperiodicboundariesareusedinthelateraldirectionsofthefluidmodel(seeFig.3),theinfluenceofwallfriction(non-slipeffect)onthesettlingbehaviourofparticlescanbeneglectedinthesimulations.

ThenumericalmodelconfigurationisshowninFig.3.Theparallelepipedhascross-sectionaldimensionsof0.025m×0.025mandaheightof1.0m.Theparticlesofvarioussizesareinitiallyrandomlygeneratedwithintheparallelepipedandthensettledownwardsundergravity.Theinputparame-tersofthesimulationsarelistedinTable2.3.1Sedimentationofasingleparticle

AsimulationofthesedimentationofasingleparticlehasbeenusedtovalidatetheDEM–CFDcouplingcode.Asphericalparticlewithradius1mmsettlesfromaposition9cmbelowtheuppersurfaceofthefluidmodel.Theporosityofthefluidmeshcellinwhichtheparticleisplacediscalculatedas0.99.Themotionofthesphericalparticleisgovernedby:󰀆󰀁4143∂Ur

πrρs=πr3ρs−ρfg−πr2ρfCdUr23∂t32

(13)

wherer=d/2istheparticleradius.

InEq.(13),thedragforcecoefficient(Cd)iscalculatedfromEq.(6).SinceUalsoappearsintheexpressionofthedragcoefficient,itisnotstraightforwardtoobtainasolu-

DEMparametersParticlediameter,d(mm)Granulardensity,ρs(kg/m3)Normalstiffness,Kn(N/m)Shearstiffness,Ks(N/m)Particlefrictionangle,φμ(◦)

Value[1.8,3.8]2,6503.0×1072.7×107301.00.1

CFDparametersFluiddensity,ρf(kg/m3)Viscosity,μ(Pa·s)SimulationparametersGravity,g(m/s2)

DEMtimestepsize,󰀍ts(s)CFDtimestepsize,󰀍tf(s)Couplingfrequency∗,α

Value1,0000.001Value−9.8110−710−5100

∗

ThecouplingfrequencyisthenumberofiterationstepsusedintheDEMinonecouplinginterval

Coefficientofrollingstiffness,βCoefficientofplasticmoment,η

123

Investigationofgranularbatchsedimentation

0.30)0.25s/m( y0.20 Analytical resultst Numerical resultsicolev0.15 gniltte0.10S0.050.000.00.10.20.30.4Time (s)Fig.4Thesettlingvelocityofaparticlewithdiameterbeing1mm

tionforthesettlingvelocityfromEq.(13)explicitly.Thus,aforwardfinitedifferencenumericaltechniqueisusedtocalculatetherelativesettlingvelocityatdifferenttimes.InFig.4,thecalculatedanalyticalresultsforsettlingvelocityarecomparedwiththenumericalones.Inthistest,thefluiddomainismeshedinthex-,y-andz-directionswith5×5×200equalsizedparallelepipedcells.Thepar-ticlesettlesfromaninitialstaticstateandthenacceleratesuntiltheterminalvelocityisreached.ThenumericalresultsmatchanalyticalonescalculatedbyEq.(13)well.Theter-minalvelocityoftheparticleis0.28m/s.

Asthecouplingmethodologyonlydescribestheaver-ageparameters(e.g.dragforce,flowvelocity,pressure)ofthefluid-solidmixture,thefluidflowaroundtheparticlesisnotexplicitlyrepresented.Duringthecalculation,thelocalporosityisassumedtobeevenlydistributedwithinonefluidmeshelement[52].Inordertogetaccurateresults,severalDEMparticlesarerequiredtofitinsideoneCFDmeshele-ment,whichmeansthatthesizeratiobetweenfluidmeshdimension(Dmesh)andparticlediameter(d)shouldbelargerthansomecriticalvalues.Aseriesofnumericalsimulationswithdifferentmeshsizeratioswereconductedtoexplorethissizeeffect:thesizeoffluidmeshwasvariedfrom0.0025to0.025m,whiletheparticlediameterwasheldconstantfordifferentsimulations).Theterminalvelocityofgrainsisnor-malizedbytheanalyticalsettlingvelocityofasingleparticle(U0),asshowninFig.5.

Thevalueofmeshsizeratio(Dmesh/d)reflectstheaccu-racyoftheaveragingprocessusedintheDEM–CFDcou-plingmodel.Ifitistoosmall,fluctuationofsettlingvelocityoccursacrossthedomainduetotheheterogeneouspackingofthefluid-solidmixture.AccordingtoFig.5,theparticlesettlingvelocitymatchesthetheoreticalvaluewhenthesizeratio(Dmesh/d)islargerthan5.Thisconclusionagreeswellwiththecriticalmesh-grainsizeratiosuggestedbyItasca[52]forthesettlingbehaviourofasingleparticle.

925

1.0)U00.8/U (ytic0.6olev g0.4niltteS0.20.001234567Size ratio (Dmesh/d)Fig.5Theeffectofsizeratio(Dmesh/d)onthenormalizedsettlingvelocity

100Sample Height (m)800.15 (base)n0.50 (middle)ah0.85 (top)t re60nif egat40necreP2001.52.02.53.03.54.0Particle diameter (mm)Fig.6Theparticlesizedistributioncurveofthesampleatthreedif-ferentlocations

Ingeneral,toconsidertheresolutionofCFDcalcula-tionandpossibleboundarywallfrictioneffects,oneneedstostudytheinfluenceofthemodelsizeratio(i.e.theratiobetweenthewidthofthefluidmodel(W)andthemeshsize(Dmesh))onthegranularsettlingbehaviour.Itasca[52]sug-gestthatthisratioshouldbenolessthan5andthattheimple-mentationofperiodicboundariesintheCFDcaneffectivelyreducetheboundarywallfrictioneffects.TheItascarecom-mendationsareimplementedinthecurrentmodel,andthenumericalresultsobtainedarethereforethoughttobeinde-pendentoftheboundarywallfrictioneffects.

3.2NumericalsimulationofbatchgranularsedimentationFornumericalsimulationsofgranularbatchsedimentation,6000polydispersedparticlesarerandomlygeneratedwithinaparallelepiped(seeFig.3).Theparticlesizedistribution(PSD)ischeckedalongtheparallelepiped.AsshowninFig.6,thePSDcurvesatthreelocationsalongtheparal-lelepiped(base,middleandtop)overlapeachother,indi-

123

926

100T.Zhaoetal.

100Percentage smaller thanPercentage smaller than8060402001.5 Base Middle TopInitial PSD curve[T] = 0.0806040200 Base Middle TopInitial PSD curve[T] = 1.82.02.53.03.54.01.52.02.53.03.54.0Particle diameter (mm)100100Particle diameter (mm)Percentage smaller than8060402001.5Percentage smaller than Base Middle TopInitial PSD curve[T] = 3.6806040200 Base Middle TopInitial PSD curve[T] = 5.42.02.53.03.54.01.52.02.53.03.54.0Particle diameter (mm)Particle diameter (mm)Fig.7Particlesizedistributionsinsuspensionsandsedimentsatdifferenttimes

catingthattheinitialDEMsampleisuniform.Theaverageporosityis0..Inthefollowinganalyses,thefluidmodelismeshedinx-,y-andz-directionswith5×5×200equalsizedfixed-gridparallelepipedcells,asabove.3.3Observationsfromnumericalanalysis3.3.1Segregationofparticles

AccordingtoStokes’lawofsedimentation[46],duringthesettlingprocess,coarsegrainssettlefasterthanthefinerones,leadingtosegregationofgrains.Toillustratethis,wepresenttheparticlesizedistributionatdifferentlocationsinthesampleinFig.7.Thesystemhasbeendividedintothreelayers(base,middleandtop),witheachlayercon-tainingthesameamountofparticles.Atthebeginningofsedimentation,segregationisnotsignificant,withthethreePSDcurvesalmostthesameasthoseoftheinitialassem-bly.At[T]=1.8,thePSDcurveofthebaseregiondevi-atesslightlybelowtheinitialPSDduetotheaccumulationofcoarsegrainsatthebottom.Asthesimulationcontin-ues,thePSDcurvesoftheupperandbottomregionsdevi-ategraduallyfromtheinitialPSDcurve.ThePSDcurveforthemiddleregionremainsalmostthesameastheini-tialPSDcurveuntil[T]=3.6,andthendeviatesgradu-allytowardsthePSDcurveofgrainsinthebaseregion.

Theseobservationsindicatethatthecoarsegrainsmovefasterthanthefinerones,leadingtoahigherconcentra-tionofcoarsegrainsatthebaseandfinergrainsatthetop.

Initially,coarsegrainsinthetoplayersettleintothemiddlelayer,whilethecoarsegrainsinthemiddlelayersettleintothebaselayer.Asaresult,onlytheportionofcoarsegrainsatthetopandbaselayershavechanged,whilethePSDinthemiddlelayerremainsunchanged.After[T]=5.4,arelativelythickdepositisformedatthebase,andthePSDcurvesofthethreelayersremainunchanged.Eventhoughsomefinegrainsarestillsuspendedinthefluid,theysettleveryslowlyontopofthedeposit.Attheendofthesimulation,onlythePSDcurvefortheupperregionliesabovetheinitialPSDcurve,whilethePSDcurvesofgrainsinthemiddleandbaselayersoverlapeachotherandlieslightlybelowtheinitialPSDcurve.Thisindicatesthatinsmallscalesedimentationsimulations,thesegregationofgrainsmanifestsitselfprin-cipallyasahigherconcentrationoffinegrainsatthetopofthesedimentlayer.3.3.2Densityprofile

Animportantfeatureofbatchsedimentationisthegradualchangeofbulkdensityofthesuspensionduetotheseg-regationofparticles.Thebulkdensitycanbecalculated

123

Investigationofgranularbatchsedimentation927

as:

ρb=nρf+(1−n)ρs

(14)

whereρbisthebulkdensityρsisthedensityofsolidmate-rials.

Duringthesimulation,theporosityineachfluidcellisrecorded,whichisusedtocalculatethebulkdensity,Eq.(14).Figure8illustratestheevolutionofthebulkdensityprofilesofthefluid-solidmixture.Itcanbeobservedthatthebulkdensityofthesuspensionincreasesgraduallyalongtheparallelepipedheighttowardsthebottomduetotheseg-regationofgrains.Inthebottomregion,grainsaccumulatetoformadensesedimentlayer,whichconsolidateprogres-

Fig.8Thedensityprofileofthefluid-solidmixtureatdifferentnor-malizedtimes

sivelyundertheself-weightoftheoverlyinggrains.After[T]=10.8,thebulkdensityatthebottomreachesanapprox-imatelyconstantvalueof1,965kg/m3.Avideoofthesimu-lationshowingtheparticlessettlingin3Disprovidedinthesupplementarymaterial.

ItispossibletoextractthepositionsofgrainsduringtheDEMsimulation.Figure9illustratesthegranularcolumnatfoursettlingtimes,inwhichagridissuperimposedonthesampletoshowtherelativelocationofthegrains.Tovisu-alizethegrainmotion,particlesatdifferentheightsoftheparallelepipedarecolouredredandgreen.Bytrackingthepositionsoftheuppermostgrains,thepositionofthefluid–suspensioninterfaceisobtained.However,asmanygrainsdecelerateandconsolidateatthebottom,thereisnowell-definedinterfacebetweensuspensionandsediment.Thus,aclosecomparisonbetweenseveralsuccessivesnapshotsofthemodeliscarriedouttofindastablepositionofthesuspension-sedimentinterfaceataspecifictime.Theevo-lutionofthesuspension-sedimentinterfaceisplottedasadashedcurveinFig.10.

BasedonFig.8,aseriesofconstantdensitypointsalongtheparallelepipedatdifferentsettlingtimeswasmappedontoFig.10.Inthisfigureitcanbeobservedthatthegeneralpat-ternsofconstantdensitycurvesissimilartothetheoreticalresultsofsedimentationbyKynch[4].Morespecifically,theuppersurfaceinitiallyacceleratesdownwardstoreachacon-stantvelocitywithinaveryshorttimeandthensettlesataconstantvelocitytowardsthebottom.Astheinitialsuspen-sionisuniformlygenerated,theheightagainsttimecurveisastraightlineinthehinderedsettlingzone.Whenthesedi-mentsapproachthebottom,theydecelerateduetograininter-actionsuntilthesettlingvelocitiesbecomezero.Inthesus-pension,thebulkdensityrangesfrom1000to1500kg/m3.

Fig.9Measurementofsedimentsheight

123

928Fig.10Heightagainsttimeplotforconstantdensitycurves

Theobservationsobtainedinthisresearchmatchtheanalyt-icalresultsbyTiller[5]well.

Asthegrainsaccumulateatthebottom,thebulkdensitychangesgraduallyfromtheintermediatetodensepackingstate,asrepresentedbythedensitychangebetweenthesus-pension,zonesofintermediatedensityandthetopofrela-tivelydensegrainlayers.Inaddition,thestablesuspension-sedimentinterfacecurvepassesthroughtheintermedi-atedensityregionandisveryclosetothe1,700kg/m3densitycurve,whichsuggeststhatastablesoilstructurecanbeformedwithdensitybeingequaltoorlargerthan1,700kg/m3.

3.3.3ExcessporewaterpressureandeffectivestressDuringthesimulation,theexcessporewaterpressure(u)alongtheparallelepipedisrecordedintheCFDmodel,basedonwhich,theevolutionoftheexcessporewaterpressureatthebottomoftheparallelepipedisobtained(seeFig.11).InFig.11,themeasuredexcessporewaterpressuresarenor-malizedbythehydrostaticpressureofwateratthebottomoftheparallelepiped(p0=ρfgH).Astheparticlesaccelerateimmediatelyaftertheapplicationofgravity,thenormalizedexcessporewaterpressureincreasesquicklytoreachthepeakvalueof0.17at[T]=0.54.Afterthattime,aloosesoilstructurebuildsupatthebottom,andtheexcessporewaterpressuredissipatesgradually.Inthisprocess,alineardissipationperiodisobservedbetween[T]=0.54and[T]=7.2.After[T]=7.2,astablesedimentlayerisformedatthebottomregion,andtheexcessporewaterpressuredissipatesslowlyduetoconsolidationofthesediment.

Consideringthefactthatsomeparticlesaccumulateatthebottomoftheparallelepipedbefore[T]=0.54inFig.11,themeasuredvalueofmaximumnormalisedexcessporewaterpressureshouldbesmallerthantheanalyticalvalue.Aratio-nalestimationofthepeaknormalisedexcessporewaterpres-

123

T.Zhaoetal.

Fig.11Evolutionofporepressureatthebottomofparallelepiped(p0isthehydrostaticpressure)

sureinFig.11canbetheinterceptofthelinearpressurelineontheverticalaxisatpointA,withanormalizedexcessporewaterpressureof0.178.Ontheotherhand,theana-lyticalmaximumnormalisedexcessporewaterpressurecanbecalculatedasthedifferencebetweentotalstressandthehydrostaticwaterpressureofthesystematthebeginningofthesimulation(seeEq.(15)).

󰀏

unρ_max/p0=

f+(1−n)ρs󰀄gH−ρfgH

e󰀊ρfgH󰀈=(1−n)ρs

ρf−1(15)

Basedontheinputparameters,theanalyticalmaximumnor-malisedexcessporewaterpressureofthefluid-solidmixturesystemiscalculatedas0.18whichcanmatchtheestimatedpeakvaluefromFig.11(≈0.178)verywell.

TheisochronesofexcessporewaterpressureatdifferenttimesareillustratedinFig.12a.Asonlywaterexistsintheregionabovethewater-suspensioninterface,theexcessporewaterpressureremainsnilinthatregion.From[T]=0.54,theexcessporewaterpressurebuildsupandvarieslinearlyalongtheparallelepiped.Astimepasses,alineardistributionofexcessporewaterpressureremainsintheupperpartoftheparallelepiped,whileatthebottom,thereisaregionofcon-stantexcessporewaterpressure.Thiscanbeexplainedbythefactthat,oncethegrainsstopmoving,thereisnorelativemotionbetweenparticlesandfluid,andthereforenoexcessporewaterpressuregradientexistsinthesediment.However,inthesuspension,theparticlescansettledownwardscontin-uously,andthus,thelineardistributionpatternremainsunal-teredthere.Aqualitativecomparisoncanbemadebetweenthenumericalresultsofourstudy(Fig.12a)andtheexper-imentalresultsfromBeenandSills[8](Fig.12b).Eventhoughthematerialsusedinthenumericalandexperimen-

InvestigationofgranularbatchsedimentationFig.12Isochronesofexcessporewaterpressurealongtheparal-lelepipedatdifferenttimes.aTheDEM–CFDsimulationresults.bResultsfromBeen[8](thetimehasbeennormalizedbythecharacter-isticsettlingtimeofagrainwithdiameterbeing5μm)

talmodelsarequitedifferent,thegeneralfeaturesoftheisochronesofexcessporewaterpressurepresentedasdimen-sionlessparametersarequalitativelythesame.

Duringthesimulation,thetotalstressactingonthebot-tomoftheparallelepipediscalculatedfromtheporewaterpressureandcontactforces,as:σ=ps+u+

FS

(16)

wherepsisthehydrostaticpressure;uistheexcessporewaterpressure;Fisthecontactforceexertedbygrainsonthebottomoftheparallelepiped;Sisthecrosssectionareaoftheparallelepiped.

Figure13showsthatthetotalstressmeasuredatthebot-tomoftheparallelepipedisalmostconstantthroughoutthesimulation.Theconstanttotalstressindicatesthattheperi-

929

Fig.13Totalstressatthebaseofthebottom

Fig.14Stressdistributionalongtheparallelepipedat[T]=3.6

odicboundaryconditionemployedinthesimulationhaveeffectivelyreducedtheinfluenceofwallfrictionontheover-allsettlingbehaviourofparticles.

AccordingtoBeen[16],thetotalstressofthesuspensionalongtheparallelepipedcanbeevaluatedbyintegratingthedensityprofiles,whiletheeffectivestressiscalculatedbysubtractingthemeasuredporewaterpressurefromthetotalstress.Figure14illustratesthedistributionofporewaterpres-sureandtotalstressalongtheparallelepipedat[T]=3.6.Atthistime,theparticlesinitiallylocatedintheupperregionoftheparallelepiped(abovepointA)havealreadysettled,suchthatonlywaterexiststhere.Thus,theporewaterpressureisequaltothetotalstress.Inthesuspension(betweenthepointsAandB),grainscaneithersettleataconstantveloc-ity(nearPointA)orcollidewitheachother,formingaloosesoilstructure(nearPointB).Theweightofparticlesispartlyorwhollysupportedbythefluidviscousdragandhydrosta-

123

930T.Zhaoetal.

ticforces.Thus,theprofileofporewaterpressuredeviatesgraduallyfromthetotalstresscurve,indicatingthateffectivestressoccurswithinthesample.Inthesedimentzone(belowpointB),thedensesoilstructurecansustaintheoverlyingloadsandarelativelylargeeffectivestressoccursthere.Therelationshipbetweenvoidratioandnormalizedeffec-tivestress(σ󰀈)withinthesoilsampleatdifferentsettlingtimesisshowninFig.15.Inthisfigure,threedistinctzonescanbeidentified:settling,transitionandconsolidationzones.Thesettlingzoneappearswhenthenormalizedeffectivestressissmallerthan0.01.Inthiszone,thevoidratiovarieswidelyfrom4.0to14,indicatingthatthesoilisextremelycompressible.Inthetransitionzone,aloosesoilstruc-turebeginstobuildupandthenormalisedeffectivestressincreasesfrom0.01to0.05,whilethevoidratiodecreasesfrom5.0to0.5gradually.However,nouniquerelationshipbetweenvoidratioandeffectivestressisobserved,whichis

closetotheobservationreportedinBeenandSills[8].Intheconsolidationzone,thecompressibilityofsoilislowandthevoidratiovariesverylittle.3.3.4Forcechains

Duringthesedimentation,thesedimentsaccumulategradu-allyatthebottomoftheparallelepipedtoformastructuredsoillayer,whichcanbevisualizedbyplottingthecontactforcechainsofthewholegranularassembly.Straightlinesareusedtoconnectthecentresofeachpairofparticlesincon-tact.Thewidthoftheselinesisproportionaltothemagnitudeofcontactforces,whiletheorientationalignswiththecon-tactforcevectors.Byplottingallthecontactforcesasstraightlines,agraphofforcechainsofthesamplecanbeobtained.AsshowninFig.16,thecontactforceisnormalizedbythecharacteristichydrostaticforceactingonthecrosssectionofaparticle(e.g.thediameterisD)atthebottomoftheparallelepiped:

F󰀁[F]=󰀆ρs−ρfgHiD2

(17)

Atthebeginningofthesimulation,thesedimentlayerisverythin,sothatthecontactforcesbetweengrainsaresmall.Asthesoilstructurebuildsup,themagnitudeofcontactforceismainlycontrolledbytheself-weightofparticles.Duetoconsolidation,thecontactforcesatthebottomregionofthemodelincreasegradually.Thestrongforcechainsofthefinalstablesoilstructurepreferablyorientvertically,indicatingthattheself-weightofsoilhasasignificantinfluenceonthedistributionofcontactforces.

4Conclusions

Fig.15Thevoidratioagainstnormalizedeffectivestressatvarioussimulationtime

Fig.16Forcechainofthesedimentsatdifferenttimes

ThispapersetsouttoinvestigategranularsedimentationviaDEM–CFDcoupledsimulations.Abenchmarksimulation

123

Investigationofgranularbatchsedimentationofthesedimentationofasinglesphericalparticleinfluidwasusedtovalidatethenumericalmodel.Inthistest,thesettlingvelocitymatchestheavailableanalyticalsolutionswell.CriteriaforassessingtheaccuracyoftheDEM–CFDcouplingcalculationsareprovided,i.e.thesizeratiobetweenfluidmeshsizeandparticlediametershouldbelargerthan5.Itisalsoshownthattheuseofperiodicboundariesinthelateraldirectionsofthefluidmodelcaneffectivelyreducethefluidboundarywallfrictioneffects.

Duringbatchsedimentation,progressivesegregationofparticlescanbevisualizedbycurvesoftheparticlesizedis-tributionatdifferentlocations.Segregationissignificantneartheupperregionofthemodel,whileitisnotsoevidentinthemiddleandbottomregions.Inthisprocess,thecoarsegrainsaccumulateatthebottom,leavingthefineronestosettleontothesurfaceofthedeposit.Asaresult,thebulkdensityofthefluid-solidmixturedecreasesgraduallywiththeheight.Attheendofthesimulation,thebulkdensityofsedimentsisaconstantvalue(i.e.1,965kg/m3).Thecurvesdescribingtheverticaldownwardtrajectoryofthefluid–suspensioninter-faceandtheincreaseofsuspension–sedimentinterfaceareinagreementwiththetheoreticalresultsproposedbyKynch[4].Asparticlescontinuouslysettledownwards,formingalooselayeratthebottomoftheparallelepipeddomain,sed-imentsconsolidateslowlyundertheweightoftheoverlyinggrains.Duringthisprocess,excessporewaterpressurebuildsupthendissipatesslowly.Thenormalizedmaximumexcessporewaterpressureofthesuspensionis0.178,whichisveryclosetotheanalyticalvalue.TheisochronesoftheexcessporewaterpressureexhibitaqualitativeagreementwiththeexperimentalresultsbyBeenandSills[8].

Inthispaper,aDEM–CFDcouplingformulationispre-sentedfortheinvestigationofgrainsedimentationinflu-ids.Manyotherapplicationsoftheformulationarepossible,e.g.submarinelandslides,mudflowsandriverscouring.ThecomputationalefficiencyofthenumericalmodeldependshighlyonthenumbersofparticlespresentandthesizeofthemeshcellsintheCFDmodel.

AcknowledgmentsThisworkofthefirstauthorissupportedbyMarieCurieActions-InternationalResearchStaffExchangeScheme(IRSES).“geohazardsandgeomechanics”,GrantNo.294976.

References

1.Bürger,R.,Wendland,W.L.:Sedimentationandsuspensionflows:historicalperspectiveandsomerecentdevelopments.J.Eng.Math.41(2),101–116(2001)

2.DiFelice,R.:Thesedimentationvelocityofdilutesuspensionsofnearlymonosizedspheres.Int.J.Multiph.Flow25(4),559–574(1999)

3.Komiwes,V.,Mege,P.,Meimon,Y.,Herrmann,H.:Simulationofgranularflowinafluidappliedtosedimentation.Granul.Matter8(1),41–54(2006)

931

4.Kynch,G.J.:Atheoryofsedimentation.Trans.FaradaySoc.48,166–176(1952)

5.Tiller,F.M.:Revisionofkynchsedimentationtheory.AIChEJ.27(5),823–829(1981)

6.Font,R.:Compressionzoneeffectinbatchsedimentation.AIChEJ.34(2),229–238(1988)

7.Richardson,J.F.,Zaki,W.N.:Sedimentationandfluidisation:partI.Chem.Eng.Res.Des.75(Supplement(0)),S82–S100(1954)8.Been,K.,Sills,G.C.:Self-weightconsolidationofsoftsoils:anexperimentalandtheoreticalstudy.Géotechnique31(4),519–535(1981)

9.Tsuji,Y.,Kawaguchi,T.,Tanaka,T.:Discreteparticlesimulationoftwo-dimensionalfluidizedbed.PowderTechnol.77(1),79–87(1993)

10.Shih,Y.T.,Gidaspow,D.,Wasan,D.T.:Hydrodynamicsofsedi-mentationofmultisizedparticles.PowderTechnol.50(3),201–215(1987)

11.Wachmann,B.,Kalthoff,W.,Schwarzer,S.,Herrmann,H.J.:Col-lectivedragandsedimentation:comparisonofsimulationandexperimentintwoandthreedimensions.Granul.Matter1(2),75–82(1998)

12.Batchelor,G.K.,Green,J.T.:Thehydrodynamicinteractionoftwo

smallfreely-movingspheresinalinearflowfield.J.FluidMech.56,375–400(1972)

13.Hinch,E.J.:Anaveraged-equationapproachtoparticleinteractions

inafluidsuspension.J.FluidMech.83,695–720(1977)

14.Batchelor,G.K.:Sedimentationinadilutepolydispersesystemof

interactingspheres.Part1.Generaltheory.J.FluidMech.119,379–408(1982)

15.Batchelor,G.K.,Wen,C.S.:Sedimentationinadilutepolydisperse

systemofinteractingspheres.Part2.Numericalresults.J.FluidMech.124,495–528(1982)

16.Been,K.:StressStrainBehaviourofaCohesiveSoilDeposited

UnderWater.UniversityofOxford,Oxford(1980)

17.Mills,P.,Snabre,P.:Settlingofasuspensionofhardspheres.Euro-phys.Lett.25(9),651(1994)

18.Davis,R.H.,Gecol,H.:Hinderedsettlingfunctionwithnoempiri-calparametersforpolydispersesuspensions.AIChEJ.40(3),570–575(1994)

19.Othmer,D.F.:Fluidization.ReinholdPublishingCorporation,New

York(1956)

20.Ergun,S.:Fluidflowthroughpackedcolumns.Chem.Eng.Prog.

48(2),–94(1952)

21.Wen,C.Y.,Yu,Y.H.:Mechanicsoffluidization.Chem.Eng.Prog.

Symp.Ser.62(62),100–111(1966)

22.DiFelice,R.:Thevoidagefunctionforfluid-particleinteraction

systems.Int.J.Multiph.Flow20(1),153–159(1994)

23.Fitch,B.:Kynchtheoryandcompressionzones.AIChEJ.29(6),

940–947(1983)

24.Kalthoff,W.,Schwarzer,S.,Ristow,G.,Herrmann,H.J.:Onthe

applicationofanovelalgorithmtohydrodynamicdiffusionandvelocityfluctuationsinsedimentingsystems.Int.J.Mod.Phys.C07(04),543–561(1996)

25.Cundall,P.A.,Strack,O.D.L.:Adiscretenumericalmodelforgran-ularassemblies.Géotechnique29(1),47–65(1979)

26.Anderson,J.D.:ComputationalFluidDynamics:TheBasicswith

Applications.McGraw-Hill,NY(1995)

27.Chen,F.,Drumm,E.C.,Guiochon,G.:Coupleddiscreteelement

andfinitevolumesolutionoftwoclassicalsoilmechanicsprob-lems.Comput.Geotech.38(5),638–7(2011)

28.Utili,S.,Crosta,G.B.:Modelingtheevolutionofnaturalcliffs

subjecttoweathering:1.Limitanalysisapproach.J.Geophys.Res.EarthSurf.116(F1),F01016(2011)

29.Boon,C.W.,Houlsby,G.T.,Utili,S.:Newinsightsinthe1963

Vajontslideusing2Dand3Ddistinctelementmethodanalyses.Geotechnique(inpress)(2014)

123

932

30.Utili,S.,Crosta,G.B.:Modelingtheevolutionofnaturalcliffs

subjecttoweathering:2.Discreteelementapproach.J.Geophys.Res.EarthSurf.116(F1),F01017(2011)

31.Brennen,C.E.:FundamentalsofMultiphaseFlow.CambridgeUni-versityPress,Cambridge(2005)

32.Weatherley,D.,Boros,V.,Hancock,W.:ESyS-ParticleTutorialand

User’sGuideVersion2.1.EarthSystemsScienceComputationalCentre,TheUniversityofQueensland(2011)

33.Abe,S.,Place,D.,Mora,P.:Aparallelimplementationofthelattice

solidmodelforthesimulationofrockmechanicsandearthquakedynamics.PureAppl.Geophys.161(11–12),2265–2277(2004)34.OpenCFDOpenFOAM—TheopensourceCFDtoolbox,http://

www.openfoam.com/

35.Šmilauer,V.,Catalano,E.,Chareyre,B.,Dorofeenko,S.,Duriez,

J.,Gladky,A.,Kozicki,J.,Modenese,C.,Scholtès,L.,Sibille,L.,Stránský,J.,Thoeni,K.:YadeDocumentation,(TheYADEProject;http://yade-dem.org/doc/)(2010)

36.Iwashita,K.:Rollingresistanceatcontactsinsimulationofshear

banddevelopmentbyDEM.J.Eng.Mech.124(285),285–292(1998)

37.Jiang,M.J.,Yu,H.S.,Harris,D.:Anoveldiscretemodelforgran-ularmaterialincorporatingrollingresistance.Comput.Geotech.32(5),340–357(2005)

38.Belheine,N.,Plassiard,J.P.,Donzé,F.V.,Darve,F.,Seridi,A.:

Numericalsimulationofdrainedtriaxialtestusing3Ddiscreteele-mentmodeling.Comput.Geotech.36(1–2),320–331(2009)

39.Shafipour,R.,Soroush,A.:Fluidcoupled-DEMmodellingof

undrainedbehaviorofgranularmedia.Comput.Geotech.35(5),673–685(2008)

40.Kafui,D.K.,Johnson,S.,Thornton,C.,Seville,J.P.K.:Paralleliza-tionofaLagrangian–EulerianDEM/CFDcodeforapplicationtofluidizedbeds.PowderTechnol.207(1–3),270–278(2011)

123

T.Zhaoetal.

41.Zeghal,M.,ElShamy,U.:Acontinuum-discretehydromechani-calanalysisofgranulardepositliquefaction.Int.J.Numer.Anal.MethodsGeomech.28(14),1361–1383(2004)

42.Drew,D.A.,Lahey,R.T.:Somesupplementalanalysisconcerning

thevirtualmassandliftforceonasphereinarotatingandstrainingflow.Int.J.Multiph.Flow16(6),1127–1130(1990)

43.DEMSolutions:EDEM-CFDCouplingforFLUENT:UserGuide

(2010)

44.Kafui,K.D.,Thornton,C.,Adams,M.J.:Discreteparticle-continuumfluidmodellingofgas-solidfluidisedbeds.Chem.Eng.Sci.57(13),2395–2410(2002)

45.Guo,Y.:AcoupledDEM/CFDanalysisofdiefillingprocess—PhD

Thesis.TheUniversityofBirmingham(2010)

46.Stokes,G.G.:MathematicalandPhysicalPapers.CambridgeUni-versityPress,Cambridge(1901)

47.Choi,H.G.,Joseph,D.D.:Fluidizationbyliftof300circularpar-ticlesinplanePoiseuilleflowbydirectnumericalsimulation.J.FluidMech.438,101–128(2001)

48.Zhang,J.,Fan,L.S.,Zhu,C.,Pfeffer,R.,Qi,D.:Dynamicbehavior

ofcollisionofelasticspheresinviscousfluids.PowderTechnol.106(1–2),98–109(1999)

49.Beetstra,R.,vanderHoef,M.A.,Kuipers,J.A.M.:Dragforce

ofintermediateReynoldsnumberflowpastmono-andbidispersearraysofspheres.AIChEJ.53(2),4–501(2007)

50.DallaValle,J.M.:Micromeritics:TheTechnologyofFineParticles,

2ndedn.PitmanPub.Corp(1948)

51.Brown,P.,Lawler,D.:Spheredragandsettlingvelocityrevisited.

J.Environ.Eng.129(3),222–231(2003)

52.Itasca:CCFD-PFCDocumentaton,1.0vols.ItascaConsulting

GroupInc,Minnesota(2008)