Inter-compartment interaction in multi-impeller mixing
Part I. Experiments and multiple reference frame CFD
Haringa, Cees; Vandewijer, Ruben; Mudde, Robert F.
DOI
10.1016/j.cherd.2018.06.005
Publication date
2018
Document Version
Final published version
Published in
Chemical Engineering Research and Design
Citation (APA)
Haringa, C., Vandewijer, R., & Mudde, R. F. (2018). Inter-compartment interaction in multi-impeller mixing:
Part I. Experiments and multiple reference frame CFD. Chemical Engineering Research and Design, 136,
870-885. https://doi.org/10.1016/j.cherd.2018.06.005
Important note
To cite this publication, please use the final published version (if applicable).
Please check the document version above.
Copyright
Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy
Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.
This work is downloaded from Delft University of Technology.
‘You share, we take care!’ – Taverne project
https://www.openaccess.nl/en/you-share-we-take-care
Otherwise as indicated in the copyright section: the publisher
is the copyright holder of this work and the author uses the
Dutch legislation to make this work public.
ContentslistsavailableatScienceDirect
Chemical
Engineering
Research
and
Design
j o u r n a l ho me p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c h e r d
Inter-compartment
interaction
in
multi-impeller
mixing:
Part
I.
Experiments
and
multiple
reference
frame
CFD
Cees
Haringa
∗,
Ruben
Vandewijer,
Robert
F.
Mudde
1TransportPhenomena,DepartmentofChemicalEngineering,DelftUniversityofTechnology,vanderMaasweg9, 2629HZDelft,theNetherlands
a
r
t
i
c
l
e
i
n
f
o
Articlehistory:
Received19November2017
Receivedinrevisedform24April
2018
Accepted1June2018
Availableonline15June2018
MSC: 00-01 99-00Keywords: Mixing Schmidtnumber CFD Multipleimpellers Macroinstability
a
b
s
t
r
a
c
t
CFDsimulationsofmixinginsingle-phasemulti-RushtonstirredtanksbasedontheRANS
methodologyfrequentlyshowanover-predictionofthemixingtime.Thishintsatan
under-predictionofthemassexchangebetweenthecompartmentsformedaroundtheindividual
impellers.SomestudiesrecommendtuningtheturbulentSchmidtnumbertoaddressthis
issue,butthisappearstobeanad-hoccorrectionratherthanphysicaladjustment,thereby
compromisingthepredictivevalueofthemethod.Inthiswork,westudytheflowprofile
inbetweentwoRushtonimpellersinstirredtank.Thedatahintsatthepresenceof
macro-instabilities,andapeakinturbulentkineticenergyintheregionofconvergentflow,which
bothmaypromoteinter-compartmentmassexchange.CFDstudiesusingthesteady-state
multiplereferenceframemodel(unsteadysimulationsaretreatedinpartII)inherentlyfail
toincludethemacro-instability,andunderestimatetheturbulentkineticenergy,thereby
stronglyover-estimatingmixingtime.Furthermore,theresultsarehighlymesh-sensitive,
withincreasingmeshdensityleadingtoapoorerpredictionofthemixingtime.Despite
properresultsfor1-impellerstudies,wedonotdeemMRF-RANSmodelssuitableformixing
studiesinmulti-impellergeometries.
©2018InstitutionofChemicalEngineers.PublishedbyElsevierB.V.Allrightsreserved.
1.
Introduction
Computationalfluiddynamics(CFD)simulationsoffera
rela-tivelycheapandfastapproachtowardsevaluatingthemixing
performanceofa rangeofimpellerconfigurations,without
theneedofalengthyexperimentalcampaign.Thisrequires
that CFD simulations sufficiently capture the true mixing
behavior,whichgaverisetoasignificantbodyofvalidation
literature.Areviewofliteraturefocusingonsinglephaseflows
withRushtonturbines(Section2)showsthatReynolds
Aver-agedNavierStokes(RANS)simulationsyielddecentresults
Abbreviations:MI,macro-instability;CoM,coefficientofmixing;LDA,laserDoppleranemometry;PIV,particleimagevelocimetry;RANS,
ReynoldsaveragedNavierStokes;MRF,multiplereferenceframes;SM,slidingmesh;(S/R)KE,standard/realizablek−;RSM,Reynolds
stressmodel;(D/L)ES,detached/largeEddysimulation.
∗ Correspondingauthor.
E-mailaddress:r.f.mudde@tudelft.nl(C.Haringa).
1 Currentaddress:DSMBiotechnologyCenter,AlexanderFleminglaan1,2613AXDelft,theNetherlands.
for single-impeller geometries (Coroneo et al., 2011), but
appeartoconsistently overestimatemixingtimes in
multi-impellergeometries(Montanteetal.,2005;Moˇst ˇeketal.,2005;
Kukukováetal.,2005;Jahodaetal.,2007);thisover-estimation
mayholdforotherimpellertypesatlargespacing(Montante
etal.,2005).LargeEddysimulations(LES)appeartoadequately
capturemixingbehavior,buttherequiredcomputationtime
prohibits routine application. We assess this
overestima-tion ofthemixingtimeinRushton-stirredtankswithlarge
impellerspacingbyRANSsimulations.Ouroriginal
hypothe-siswasthatofRANSsimulationsunder-estimatingtinthe
https://doi.org/10.1016/j.cherd.2018.06.005
Nomenclature Roman
an dampingcoefficient,Eq.(3),–
Ct tracerconcentration,kg/m3
C off-bottomclearanceimpeller,m
cn strengthcoefficient,Eq.(3),–
C inter-impellerclearance,m
D impellerdiameter,m
D diffusioncoefficient,m2/s
Dt diffusioncoefficient,turbulent,m2/s
Ep/Et fractionofperiodicenergyinspectrumofu,–
f frequency,s−1 f1.1 basefrequency1,s−1 f2.1 basefrequency2,s−1 f1.2 harmonicfrequency1,s−1 f2.2 harmonicfrequency2,s−1 FQ pumpingnumber,– h axialcoordinate,m H tankheight,m
k slotnumber(autocorrelation),–
kt turbulentkineticenergy,m2/s2
kMI macro-instabilitykineticenergy,m2/s2
kt,* kt+kMI,m2/s2
kt,fit ktcomputedviaauto-correlationfit,m2/s2
kt,SI ktcomputedviaspectralintegration,m2/s2
kMI,fit kMIcomputedviaauto-correlationfit,m2/s2 kMI,SI kMIcomputedviaspectralintegration,m2/s2
N impellerrevolutions,1/s
Nc no.gridcells,–
Qax axialflowrate,L/s
r radialposition,m
R tankradius,m
Sct turbulentSchmidtnumber,–
T tankdiameter,m
t timestepsize,s
u fluctuatingvelocity,m/s
U meanvelocity,m/s
Utip impellertipspeed,m/s
V tankvolume,m3
Vi gridcellvolume,m3
w Tukey–Hanningwindow,−
Greek
laserwavelength,nm
viscosity,dynamic,Pas
t viscosity,dynamic,turbulent,Pas
turbulentdissipationrate,m2s3
ˆ auto-correlationcoefficient,–
l density,kg/m3
lag lagtime,s
noisecomponent,–
95 mixingtime,s
95 mixingnumber,N·95−
horizontal plane segregating the compartments formed
around the individual impellers, hence under-estimating
mass exchange between them. During our investigation,
we additionally observed a macro-instability in the
inter-compartmentregion,whichsignificantlycontributedtomass
exchange. Inseveral prior publications,inter-compartment
mass exchange was boosted by tuning of the turbulent
Schmidtnumber,Sct=t/(Dt)(Montanteetal.,2005;Gunyol
etal.,2009;Haringaetal.,2016).Weposethatthistuningis
notbasedonphysicalreasoning.Instead,it isapatchwork
solutioncoveringexchangebyvirtueofthemacro-instability,
aswellasunder-predictedturbulentexchange,asTominaga
andStathopoulos(2007)similarlydiscussedformispredicted
scalarspreadinginjets.Suchanad-hoccorrection
compro-misesthepredictivecapabilitiesofRANS-basedstirredtank
simulations. Whereas other publications on multi-impeller
mixingfocusedprimarilyonoverallmixingbehavior,westudy
theflowinregionbetweentheimpellersindetail,usinga
com-bination ofLaser DopplerAnemometry (LDA) experiments,
RANS andLESsimulations.Wecombinetheseresultswith
both steady and unsteady CFD simulations, to answer as
to why RANS simulations poorly capturemixing in
multi-impeller tanks. In this first part of the work, wefocus on
theexperimentalassessmentandsteadyRANSsimulations,
unsteady simulations including LES are considered in part
two.
2.
Literature
overview:
RANS
simulations
of
Rushton
turbines
Inthisliteratureoverview,wediscusstheuseofRANS
simula-tionstostudymixinginvesselsstirredwithRusthonturbines.
LESandDESmodelsarethesubjectofpartIIofthiswork.
2.1. Theimpellerdischargestream
Mosthydrodynamicresearchinstirredtanksfocusedonthe
impeller discharge stream and trailing vortices, being the
regions wheremuchofthe hydrodynamic actionhappens.
Experimentstypicallymeasuredtheturbulentkineticenergy
(kt)anddischargestreamvelocitybyLDA(WuandPatterson,
1989;Ng et al.,1998;Ng andYianneskis, 2000;Murthy and Joshi, 2008; Schäfer et al., 1997; Vennekeret al., 2010; Lee andYianneskis,1998;Derksenetal.,1999)orparticleimage
velocimetry (PIV)(Khopkaret al., 2004;Ranadeetal., 2001;
Escudié and Liné, 2003; Escudié et al., 2004; Huchet et al., 2009; Baldi and Yianneskis, 2004; Baldi et al., 2004; Ducci andYianneskis,2005;Linéetal.,2013).Earlystudiestypically
computed the turbulent dissipation rate on dimensional
grounds(WuandPatterson,1989),orviatheenergybalance
(EscudiéandLiné, 2003;Escudiéetal., 2004). Direct
mea-surements have been conducted with high-resolution LDA
andPIV(Huchetetal.,2009;BaldiandYianneskis,2004;Baldi et al., 2004; Ducciand Yianneskis, 2005). These show that
estimationsoffromdimensionalanalysisarereasonablein
thebulk,butunderestimateclosetotheblade,where
tur-bulenceisstronglyanisotropic(Hartmannetal., 2004).The
bulkofCFDstudiesreportconsensusregardingtheimpeller
dischargevelocity(GunyolandMudde,2009;Ngetal.,1998;
ZakrzewskaandJaworski,2004;Ammaretal.,2011;Laneetal., 2000;Brucatoetal.,1998;JenneandReuss,1999),while
turbu-lenceparametersyieldlessuniversalagreement.Brucatoetal.
(1998)underestimatedktusingthestandardk−(SKE)model,
usinginner–outer(IO)andslidingmesh(SM)impeller
model-ing.GunyolandMudde(2009)reportedpooragreementusing
RNG–KEbutfoundgoodresultsforktwithSKEandthe
realiz-ablek−(RKE)models,withmultiple-referenceframe(MRF)
impeller modeling. TheReynolds stress model (RSM)
com-paredlessfavorably,butoutperformedRNG–KE.Singhetal.
(2011)notedktinthedischargestreamisreasonablywell
poorer,withthetrailingvorticesbeingtoosmall,and
inward-curvedcomparedtoexperimentaldata(Derksenetal.,1999;
Ranadeet al., 2001;Escudiéet al., 2004). Underestimations
inofup to50% werereportedbyNgetal.(1998),Ngand
Yianneskis(2000),whomodeledthefull3Dimpeller
geome-try.With1–3cellsacrossthebladethickness,poorlyresolved
flowaroundtheblademaycausethisunderestimation.
Stud-iesusingsheet-bodyimpellersandbafflesperformedbetter
(GunyolandMudde,2009;Coroneoetal.,2011;Singhetal.,
2011), although high mesh densities (106+ cells for a one
impeller,360◦ domain) are required formesh independent
results(Coroneoetal.,2011).Forthese,resultsforboththe
SKEandRKEmodelagreewellwiththedimensional
assess-mentofbyWuandPatterson(GunyolandMudde,2009).This
isunsurprising,asthek−methodisbuiltlargelyonthesame
assumptions (isotropy and a single turbulent length-scale
Lres=Ak3/2t /).Singhet al.(2011)(SM–SKE)showedthatthe
agreementinbreaksdownnearthebladewhencomparedto
directmeasurements.Furtherawayfromtheimpeller,where
thek−modelassumptionsreasonablyhold(Hartmannetal.,
2004),decentagreementinwasobservedbySinghetal.
2.2. MixinginRushton-stirredtanks
2.2.1. Definitionofmixingtime
We report the 95% mixing time in dimensionless form;
95=N·95 with 95 the mixing time ins and N the
agita-tionrateins−1.Forpoint-measurementsweconsider95as
thetimewherethenormalizedtracerconcentrationCt/ ¯Ct is
between0.95and1.05( ¯Ctthevesselaverage)(Kukukováetal.,
2005;Jahoda etal., 2007; Bujalski etal., 2002a,b,c; Jaworski etal., 2000a,b;Magelli etal., 2013). Probeshavethe
advan-tageofrecordingdynamicsthatare directly comparableto
CFDsimulations,buttheresultsmaybelocationspecific,and
invasiveprobesmaydisturbthelocalflow.Thecoefficientof
mixing(CoM)(Kukukováet al.,2008;Kukukova etal., 2009;
Hartmannetal.,2006)(alternativelycalledCoefficientof
Vari-ation(CoV))isusedtoquantify95 basedonCFDfielddata.
AsCFDmeshesaretypicallynon-uniform,volume-weighing
shouldbeappliedinthedeterminationoftheCoM(Hartmann
etal.,2006): CoM=
⎛
⎜
⎝
i Ct,i− ¯Ct ¯ Ct2 Vi iVi
⎞
⎟
⎠
(1)Kukukováet al.(2008) set CoM=0.05 asthe 95%mixing
limit, while Hartmann et al. (2006) report CoM=0.0283 as
95% of the volume being 95% mixed, based on numerical
experiments. We apply the limit of Hartmann, but do not
considereitherindicatortobesuperiorprovidedconsistency
isapplied.Coroneoetal.(2011)alsousedtheCoMto
quan-tifymixingbyplanelaserinducedfluorescence(P-LIF).Some
studieswithalternativemixingindicatorsarefoundin
litera-ture,suchas95%ofthevolumebeingdecolored(Moo-Young
etal.,1972).LeeandYianneskis(1997)suspendedthermally
sensitive liquid-crystal particles and applied a heat pulse,
measuring 95when95%oftheparticleshadthesamehue.
2.2.2. Singleimpellersystems
The mixing time correlation of Yeoh et al. (2005) yields
95=30.7forasingleRushtonturbine,withD=T/3,H=T.The
Ruszkowski–Grenvillecorrelation(Yeohetal.,2005),validfor
variousimpellertypes(Nienow,1997),predicts 95=27.5.
Typ-icalexperimentalandCFDresultsgive 95≈30–40(Table1),
showingthecorrelationsarereasonable.Thereisquitesome
variationintheresultshowever.Inexperiments,thismaybe
duetoprobelocation(RaghavRaoandJoshi,1988;Kukuková
et al., 2005; Jahodaet al., 2007), or dueto different
meth-odsusedto quantifythemixing time(Leeand Yianneskis,
1997;Moo-Youngetal.,1972),thatarenotmutually
compa-rable.SimilarobservationscanbemadefortheCFDresults.
LESmethods –discussed inmoredetail inpart II–
gener-allyperformwell,butdifferentmixingquantificationmethods
yielddifferent results.Especially 3D-CoM(Hartmannet al.,
2006))yieldsahigher95,asitconsiderstheentirevolume,
includingpoorlymixedregions.Thiscallscautionwhen
com-paringmixingtimesbetweendifferentstudies,experimental
ornumerical.FocusingonRANSmethods,slidingmesh
sim-ulationswerereportedtoover-predict 95(Jahodaetal.,2007;
Zadghaffarietal.,2010).Incontrast,theMRF–SKEapproach
ofJahodaetal.andKukukovaetal.underestimated 95
com-pared with their measurements. The agreement between
experimentalprobemeasurementsandnon-invasiveLIF
mea-surements(Distelhoffet al.,1997;Javed etal., 2006)makes
itunlikelythesedifferencesareduetotheprobeinfluencing
flow.Morelikely,themeshresolutionwasinsufficient.Results
byCoroneoetal.(differentgeometry:D=T/3,C=T/2,H=T)do
showgoodagreementbetweenPLIFdataandMRF–SKE
simu-lations(2D-CoM)qualitativelyandquantitatively.Theirmesh
of2000kcellswasmuchfinerthanusedinotherMRF
stud-ies.Coroneo et al.foundthe discharge velocityprofileand
torque-basedpowernumberPo,oftenusedformesh
depen-dencystudies,areratherinsensitivetomeshdetails(Coroneo
etal.,2011).ThedissipationbasedpowernumberPoismore
sensitive,andinbetteragreementwiththemeshsensitivity
of 95.
2.2.3. Multi-impellersystems
Theflowinmulti-impellersystemsdependsstronglyonthe
inter-impeller clearance (Rutherford et al., 1996). We focus
on Rushton-stirred systems with a spacing C=T, which
exhibitparallelflow(MishraandJoshi,1994;Rutherfordetal.,
1996),formingseparatecompartmentsaroundtheindividual
impellers.Werefertotheregionwherethesecompartments
meetastheinter-compartmentregion.withtheexact
hori-zontalseparationbetweenthecompartmentsreferredtoas
the inter-compartmentplane.In thisplane, themean flow
isdominantlyradial,andlittlemassexchangebetweenthe
compartmentsisexpectedbythemeanconvectiveflow.This
formsarate-limitingstepinmixing;comparingthe
experi-mental 95between1-,2-and4impellersystems(Table1–3,
respectively)shows 95morethandoublesupondoublingthe
numberofimpellers.Allsummarizedmulti-impellerstudies
report 95 basedon invasiveprobes,withtracerinjectedat
thevesseltop.Incaseofmultipleprobes,thebottomprobe
reading isreported, givingthehighest 95.Early studiesby
Jaworskietal.(2000a,b)andBujalskietal.(2002a,c),withvery
limitedmeshresolution,reporteduptoafactor2
overestima-tionof 95,bothwithSM–SKEandMRF–SKE.Kukukováetal.
(2005)reportedan8.6%over-estimationforMRF–SKE,while
Jahodaetal.(2007)observed20%and26%over-estimations
forMRF–SKEandSM–SKE,respectively.TheLESsimulations
ofJahodaetal.(2007) showedhighlysimilarprobe
dynam-icscomparedtotheirexperiments;thehigherdegreeofnoise
in experimental data did in an 11.4% under-estimationof
4-Table1–Anoverviewofmixingtimesforasingleimpellergeometry(D=T/3,C=T/3,H=T):experimental(t95).Measure
indicatesthemethodusedtoquantify95:x-P:xprobes,T:thermal,C:decolorization,3D-CoM:coefficientofmixing,
volume.MM:min–maxmethod.
Author 95 Measure Re Method Sct Mesh Ref.
Kukukovaetal. 33 1-P 46,667 Conductivity − − (Kukukováetal.,2005)
Jahodaetal. 34 1-P 46,667 Conductivity − − (Jahodaetal.,2007)
Leeetal. 21.6 T 40,000 Thermal − − (LeeandYianneskis,1997)
Distelhoffetal 36 8-P 24,000 LIF − − (Distelhoffetal.,1997)
Javedetal 32 32-P 24,000 LIF − − (Javedetal.,2006)
Moo-Youngetal. 36 C 24,000 Decolorization − − (Moo-Youngetal.,1972)
RaghavRaoetal. 48.2(avg.) 1-P >100,000 Conductivity − − (RaghavRaoandJoshi,1988)
Rewatkaretal. 39.5(avg.) 1-P >200,000 Conductivity − − (RewatkarandJoshi,1991)
Nereetal. 35.3 MM 46,722 B.C.k− Notreported 280k (Nereetal.,2003)
Kukukovaetal. 23.5 1-P 46,722 MRFk− 0.7 330k (Kukukováetal.,2005)
Jahodaetal. 13.5 1-P 46,722 MRFk− 0.7 615k (Jahodaetal.,2007)
Jahodaetal. 50 1-P 46,722 SMk− 0.7 615k (Jahodaetal.,2007)
Javedetal 27 32-P 24,000 SMk− 0.7 112k (Javedetal.,2006)
Zadghaffarietal 40–80 1-P 41,667 SMLES 0.7 971k (Zadghaffarietal.,2010)
Jahodaetal. 26 1-P 46,722 SM–LES 0.7 615k (Jahodaetal.,2007)
Yeohetal. 33.2 48-P 4000 SM–LES 0.7 330k (Yeohetal.,2005)
Hartmannetal. 39.8–54.7 3D-CoM 24,000 SM–LES 0.7 13,824k (Hartmannetal.,2006)
Table2–Anoverviewofmixingresultswith2impellergeometries:CFDandexperimentaldata.Unlessotherwise stated,D=T/3andC=T.Allmixingtimemeasurementswereprobe-based.ForturbulentSchmidtnumberSct,n.r.=not reported.
Author Geometry 95 Re Method Sct Mesh Ref.
Kukukovaetal. C=T/3 104.5 46,722 Conductivity − − (Kukukováetal.,2005)
Kukukovaetal. C=T/3 113.5 46,722 MRFk− 0.7 286k (Kukukováetal.,2005)
Jahodaetal. C=T/3 92 46,722 Conductivity 0.7 − (Jahodaetal.,2007)
Jahodaetal. C=T/3 110 46,722 MRFk− 0.7 1230k (Jahodaetal.,2007)
Jahodaetal. C=T/3 116 46,722 SMk− 0.7 1230k (Jahodaetal.,2007)
Jahodaetal. C=T/3 81.5 46,722 SMLES 0.7 1230k (Jahodaetal.,2007)
Jaworskietal. C=T/4,D=T/2 46.6 162,000 Conductivity − − (Jaworskietal.,2000a,b)
Jaworskietal. C=T/4,D=T/2 48.9 216,000 Conductivity − − (Jaworskietal.,2000a,b)
Jaworskietal. C=T/4,D=T/2 137.5 162,000 SMk− n.r. 70k/180◦ (Jaworskietal.,2000a,b)
Jaworskietal. C=T/4,D=T/2 125.2 216,000 SMk− n.r. 70k/180◦ (Jaworskietal.,2000a,b)
Jaworskietal. C=T/4,D=T/2 140 324,000 SMk− n.r. 70k/180◦ (Jaworskietal.,2000a,b)
Bujalskietal. C=T/4,D=T/2 176 216,000 MRFk− n.r. 115k (Bujalskietal.,2002a,c)
Table3–Anoverviewofmixingresultsfor3and4impellergeometries.Inallcases,H=n·Twherenisthenumberof impellers,C=TandD=T/3.
Author Geometry 95 Re Method Sct Mesh Ref.
Jahodaetal. 3Imp.,C=T/2 197 21,000 Conductivity − − (JahodaandMachon,1994)
Jahodaetal. 3Imp.,C=T/2 186 25,210 Conductivity − − (JahodaandMachon,1994)
Jahodaetal. 3Imp.,C=T/2 212 33,610 Conductivity − − (JahodaandMachon,1994)
Jahodaetal. 3Imp.,C=T/2 210 42,015 Conductivity − − (JahodaandMachon,1994)
Jahodaetal. 3Imp.,C=T/2 204 50,415 Conductivity − − (JahodaandMachon,1994)
Mosteketal. 3Imp.,C=T/2 239 46,722 Conductivity − − (Moˇst ˇeketal.,2005)
Mosteketal. 3Imp.,C=T/2 286 46,722 MRFk− 0.7 1537k (Moˇst ˇeketal.,2005)
Mosteketal. 3Imp.,C=T/3 236 46,722 Conductivity − − (Moˇst ˇeketal.,2005)
Mosteketal. 3Imp.,C=T/3 273 46,722 MRFk− 0.7 1537k (Moˇst ˇeketal.,2005)
Jahodaetal.,1994 4Imp.,C=T/2 400 21,000 Conductivity − − (JahodaandMachon,1994)
Jahodaetal.,1994 4Imp.,C=T/2 387 25,210 Conductivity − − (JahodaandMachon,1994)
Jahodaetal.,1994 4Imp.,C=T/2 392 33,610 Conductivity − − (JahodaandMachon,1994)
Jahodaetal.,1994 4Imp.,C=T/2 405 420,15 Conductivity − − (JahodaandMachon,1994)
Jahodaetal.,1994 4Imp.,C=T/2 402 50,415 Conductivity − − (JahodaandMachon,1994)
Montanteetal. 4Imp.,C=T/2 461 37,500 Conductivity − − (Montanteetal.,2005)
Mosteketal. 4Imp.,C=T/2 437 46,722 Conductivity − − (Moˇst ˇeketal.,2005)
Montanteetal. 4Imp.,C=T/2 476 37,500 SMk− 0.1 350k (Montanteetal.,2005)
Mosteketal. 4Imp.,C=T/2 466 46,722 MRFk− 0.7 1984k (Moˇst ˇeketal.,2005)
Mosteketal. 4Imp.,C=T/3 439 46,722 Conductivity − − (Moˇst ˇeketal.,2005)
impellersystems.Montanteetal.requiredastronglylower
turbulentSchmidtnumber,Sct=0.1(Montanteetal.,2005),to
yieldagreementwithexperimentsina4impellertank;
simi-larobservationsweremadeinotherstudies(Delafosseetal.,
2014;Gunyoletal.,2009).Moˇst ˇeketal.(2005)didnottuneSct,
andreportedanapproximately20%overestimationin 95
com-paredtotheirmeasurementswith3impellers(40%compared
toJahodaandMachon(1994)).Theyfounddecentagreement
in4impellersystemscomparedtotheirownmeasurements
(theydidover-estimate 95comparedtoJahodaandMachon
(1994)).AnimportantobservationintheworkofMosteketal.
for3and4impellers,isthatincreasingmeshdensityledtoan
increasein 95,withmeshindependencenotyetreachedatthe
finestmeshes.Thenumberofgridcellsinallmulti-impeller
studiesreportedabovewaslowerthaninthesingle-impeller
studyofCoroneoetal.(2011),implyingthatfurtherrefinement
mayleadtoaconsiderablystrongerover-estimationof 95in
multi-impellersystems;thedecentresultsofMosteketal.for
4impellersmaybeduetoinsufficientresolution,ratherthan
inherentlygoodmodelperformance.Hydro-dynamically,the
inter-compartmentregionisunder-studied;studiestypically
focusontheimpellerdischargestreamratherthanthe
qui-escentbulkregion(Delafosseetal.,2014).Micaleetal.(1999)
reportedasmallpeakinktintheinter-compartmentregion
withLDAdata, possiblyduetothe generationofturbulent
kinetic energy by the colliding circulation loops. This
tur-bulencemayenhancemixingacrosstheinter-compartment
plane,andhencedeservesmoreattention:itmaybean
under-estimation ofturbulence inthe compartment segreregion
leadstopoormixingpredictionswithRANSmodels,masked
bytheeffectsofinsufficientmeshresolution.
A second phenomenon that may impact
inter-compartmentmixing isthe presence ofMacro-instabilities
(MIs).Thesehavebeenextensivelystudiedinsingle-impeller
geometries,bothexperimentallyandnumerically(Roussinova
et al., 2003; Hartmann et al., 2004; Nikiforaki et al., 2004).
Nikiforakietal.(2003)suggestjetinstabilitiesorinstabilities
by precessing vortices are the dominant cause. Paglianti
et al. (2006, 2008) report boththese effects are present in
Rushton-stirredtanks.As forbothphenomenafrequencyf
scaleslinearwithagitationrateN,MIfrequenciesarereported
hereinf/N.Inthedischargestreamofa2impellersystem,
frequenciesoff/N≈0.02andf/N≈0.055arefoundfor
precess-ingvorticesandjetinstabilities,respectively(Pagliantietal.,
2008). Guillard et al. similarly identified instabilities with
f/N≈0.05−0.08intheregionabovethetopimpeller(Guillard
etal.,2000a,b).Ifthesemacro-oscillationsleadtooscillations
intheinter-compartmentplane, theymaypromotemixing
betweenthecompartments.Byconstruction,frozen-flowfield
MRFsimulationsarenotabletocapturetheseoscillations.
2.3. Hypothesis
Noconsistent under- or over-estimationin 95 is observed
for RANS simulations with 1 impeller, and Coroneo et al.
(2011)showedquantitativeandqualitativeagreementin
mix-ingbehaviorcanbeobtainedwithMRF–SKE,givensufficient
meshresolution.Incontrast,multi-RushtonRANS-CFD
stud-ies showa consistentover-predictionof 95, which incase
ofMRFsimulationworsenswithincreasingmeshresolution.
Thishintsthattheoverestimationof 95insuchsystems
orig-inatesfrom poorassessmentinter-compartmentexchange,
whileintra-compartment can be properly captured. In the
inter-compartment plane, the flow is dominantly parallel
Table4–Meshesusedinthiswork.2IFrepresentsa 360◦domainand2IPrepresentsthe60◦domain.Thelast letter(s)representthemeshquality(C=crude,
M=medium,F=fine,SF=super-fine).Ncisthenumber ofgridcells.
Name Nc Domain Methods
2IP-C 94k 60◦ RKE,RSM 2IP-M 506k 60◦ RKE,RSM 2IP-F 812k 60◦ RKE,RSM, SM-RKE 2IF-C 648k 360◦ RKE,RSM, SM-RKE,LES 2IF-M 1997k 360◦ RKE,RSM, SM-RKE, SM-RSM,LES 2IF-F 5884k 360◦ RKE,RSM, SM-RKE
2IF-SF 10,584k 360◦ RKE,RSM,LES
meaning that turbulence and/or MIs control mass trans-portbetweenthecompartments.InlinewithJaworskietal. (2000b), we hence pose the hypothesis that RANS models
inherentlyover-estimate 95 by(a)anunder-predictionoft
duetoincorrectassessmentofturbulenceintheparallelflow
regionand/or(b)failingtocapturetheeffectofMIson
inter-compartmentmixing.
3.
Materials
and
methods
3.1. CFDsetup
We focus on the realizable k− (RKE) and Reynolds stress
model(RSM) withlinear-pressurestrainformulation.These
modelsarewellestablishedinstirredtankmodeling;a
sum-marycanbefoundinGunyolandMudde(2009),forexample.
A2-impellerstirredtankwithheightH=2T,impeller
diame-terD=T/3,off-bottomclearanceC=T/3andimpellerclearance
C=Twasmodeled.Here,tankdiameterT=0.29m,asused
byJahodaetal.(2007).Thetankcontains4baffleswithwidth
T/10.Allinternalsweremodeledassheetbodies(Gunyoland Mudde,2009;Coroneoetal.,2011).Botha360◦ domainand
60◦domainweremodeledtoassesstheinfluenceofa
symme-tryassumption(GunyolandMudde,2009;Gunyoletal.,2009;
Haringaetal.,2016,2017).Usinga60◦sectionintroduces6
baf-fles;theeffectofthisadditionalbafflingwasfoundtobesmall
(GunyolandMudde,2009).AsinJahodaetal.(2007),thebottom
tracerconcentrationprobeisplacedatT/4fromthebottom,
betweenthebaffles,atT/20fromthewall.Thetopprobeis
placedat1.25T,atthesameradialandangularposition.The
agitationrateissettoN=5s−1.InordertousetheMRFmethod,
the domain isseparated inastationaryzone, androtating
zonesaroundtheindividualimpellers.Theserotatingzones
donotoverlapwiththeinter-compartmentplane.Structured
hexahedralmesheswereused;thegridsizesNcarereported
inTable4,withthe2IFseriesbeingthefulldomain,andthe2IP
seriesrepresentinga60◦slice.Meshdetailsarefurther
spec-ifiedinSupplementarymaterialB.Spatialdiscretizationwas
setto2ndorderupwind(GunyolandMudde,2009;Coroneo
etal.,2011)andstandardwallfunctionswereemployed.The
vessel topwas setto no-shear tomimicafreesurface, all
other walls were no-slip. Convergence was declared when
theresiduals<10−5andthemeanvelocityremainedwithin
0.1%over1000iterations.Thevelocityandturbulencefields
pas-Fig.1–MeasurementgridfortheLDA-experimentsinthe mid-compartmentplane.Axialpositionsgivenare
non-dimensionalizedwiththetankdiameterT,radialwith tankradiusR=D/2.y=0isthecentralpositionbetweenthe impeller(themostlikelylocationoftheinter-compartment plane),atheight5T/6fromthebottom,andr=0the impellershaft.
sivescalarinasphericalvolume(witharadiusof0.0125m),
atheighth=0.551mfromthebottom,atr=0.0725m,inthe
baffleplane.Atimestepsizeoft=0.005swasusedwith
sec-ondorder implicittime discretization.Thetracer and bulk
fluidhad equalproperties, =1000kg/m3 and =0.001Pas,
suchthatthetracerwillnotdisturbtheflowfield.TheFLUENT
defaultsimplegradientdiffusionhypothesis(SGDH)wasused
todetermineturbulentscalardiffusion,modelingthediffusive
fluxasJi=−(lD+t/Sct)∇CtwithD=10−9m2/sthe
molecu-lardiffusioncoefficient,ttheturbulentviscosity,Ctthescalar
concentration.Sct=0.7istheturbulentSchmidtnumber.An
unpublishedsimulationwiththegeneralizedgradient
diffu-sionhypothesis(GGDH)yieldednosignificantchangein 95,
bothforRKEandRSMsimulations.
TheaxialflowrateQaxthroughthemid-impellerplane(y=0
inFig.1)wascomputedbyiso-clippingthissurface,retaining
onlylocationswithuax<0(downwardflow),andreportingthe
totalmass-flowratethroughtheremainingfacesinFLUENT.
Duetomassconservation,anequalvaluewasfoundwhen
iso-clippingtoretainonlyupwardflow.
3.2. LDAsetup
Weuseda2-probeTSIpowersightlasersystem(150mW)with
1=561nm(axial)and2=532nm(radial),operatedin
back-wardscatteringmode.Aglasstank(T=0.26m,H=2T,C=T/3,
C=T,4aluminumbafflesofT/10)wasplacedinawater-filled
rectangularglassboxforrefractiveindexmatching.Two
stan-dardRushtonturbinesweremountedonacenteredshaftof
diameterds=0.02m.Agraphicallayoutofthesetupisgiven
inSupplementarymaterialA.Theexperimentsreported in
thisworkwereconductedaftertheCFDsimulations.Dueto
the equipment available inour lab,the experimental tank
diameterisslightlysmallerthanthe diameterusedinCFD
simulations. All geometric ratios are equal, however. The
experimentswereconductedatN=5.78s−1duetolimitations
ofthemotor,givingRe=4.34×104,comparedto4.67×104(at
N=5.0s−1)byJahodaetal.(2007)andinourCFDwork.Wedo
notexpectthesedifferencestoinfluencetheresults,asboth
simulationsandexperimentsareconductedinthefully
turbu-lentregime,andresultsarepresentedindimensionlessform.
Hollowglassseeding(dp=8–12m,StokesnumberSt≈0.03)
wasused.
Verificationmeasurementswereconductedintheimpeller
outflowat6radialpositions(r=[50,60,70,77,90,102]mm),at
theheightoftheimpellerdiscofbothimpellers,inthe
baf-fle plane.Theaveragedataratewas 340–500Hz(axial) and
830–850Hz(radial).Dataintheinter-compartmentregionwas
recordedat7radialand11axialpositions,whichare
speci-fiedinFig.1,withy=0locatedatheight5T/6,measuredfrom
thebottom.Eachmeasurement50,000datapointswere
col-lected,measurementsweredonein5-fold.Theaveragedata
rateswere289Hz(axial)and343Hz(radial).Theseare
insuf-ficienttoresolvethefullturbulencespectrum,butsufficeto
measuretheenergy-carryingmotionswhichareexpectedto
berelevantforinter-compartmentmassexchange.
3.2.1. Dataprocessing
LDA data processing was conductedin MATLAB8.6.0. The
setupcontainedaburstcountertoremovefalseregistrations,
thevelocitybiaswascorrectedusinggate-timeweighing.The
meanvelocityisretrievedbytime-averagingthevelocity
sig-nal;subtractingthemeanyieldsthefluctuatingvelocity.The
contribution of noise and periodic flow components were
assessedviatheslotted auto-correlationmethodwithlocal
variancenormalization(TummersandPasschier,1772c;van
Maanenetal.,1999),withtheauto-correlationˆcomputedas:
ˆ(klag)=
u(ti)u(tj)
u2(ti) u2(tj)
(2)
withkanintegerslotnumber.Ideally,atzerolag-timeˆ=1;
(white)noisecausesalowervalueinpractice,thedifference
isusedtoassessthenoisecontribution.Tosegregate
peri-odicflowcontributions,vanMaanen (1999)proposed fitting
theauto-correlationfunctionwithadampedcosines,Eq.(3):
ˆfit()=b+c0e−˛0+ in=1cne−˛ncos(2·nf ) (3)
Summingtheamplitudes cnthengivesthecontribution
ofperiodicmotionstothefluctuatingkineticenergy.
Alterna-tively,thekineticenergycontributionofperiodicflowcanbe
estimatedfromthespectraldensityfunctionEq.(4):
S(f)= t
12ˆ(klag)w(klag)cos(kflag)(4)
WithwaTukey–Hanningwindow(TummersandPasschier,
1996, 1772c) with width tm(f)=tmax f0/f, where lower
decreases low-frequencyresolution, but high can induce
modulation We found =5provided a good trade-off.The
relative contribution of periodic components follows from
integratingEq.(4)overtherelevantfrequencies,whichwas
doneusingtrapezoidalintegration.Withthesemethods,the
RMSvelocityurms=
¯
uu canbecorrectedfortheperiodic
components,byEq.(5)whenperiodicfittingisusedandEq.
(6)whenusingspectralintegration:
ut=u·
1−− (cn) (5)
ut=u·
1−−(Ep/Et) (6)
Here,Ep/Etisthefractionofenergycontainedinthe
peri-odicspectralcomponentsofu. Wereportkt astheMI-free
fluctu-Fig.2–LDAresults:currentdata(topandbottomimpelleroutflow)comparedwithearlierstudiesbyWuandPatterson (1989),MurthyandJoshi(2008).Thepositionisscaledwithtankandimpellertipradius,UradandktwiththetipspeedUtip.
atingenergy basedon u, and kMI as the kinetic energyof
macro-instabilities,kt*−kt.Sincethetangentialvelocityisnot
measured,isotropyisassumedtoyieldkt=3/4( ¯u2t,ax+ ¯u2t,rad).
Experimentally, we use the subscripts fit and SI to
distin-guishbetweenkt and kMI computedwithcorrelationfitting
andspectralintegration,respectively.Intheimpelleroutflow
ˆ wasevaluated forlag=250mswith500 slots,which
suf-ficedtocapturenoiseandbladepassages.Thewell-defined
impeller-inducedperiodicitygavepreferencetoEq.(5)for
peri-odicestimation,taking the first 3harmonicsinto account.
In the inter-compartment region, ˆ was evaluated first for
lag=250mswith500slotstoestimatenoise(noimpeller
influ-encewasobserved),andnextatlag=40swith1600slotsto
evaluateMI-inducedperiodicity.Asahigh-resolution,longlag
timecomputationwasunfeasible,weestimatedthetotal
spec-traldensityfunctionbymatchingthespectrumofbothlag
timerangesintheoverlappingfrequencyrangetoconstruct
asinglespectrum.Intheinter-compartmentregion,ˆwas
fit-tedwith4cosines,at2frequencies(f1.1andf2.1)andtheirfirst
harmonics(f1.2andf2.2).Spectraldensityintegrationrequiresa
cut-offbetweenturbulentandMIcomponents.Weconsidered
therangef/N<0.1(10impellerrevolutions)asMI-components,
basedonvisualassessmentofthespectra.Ofcourse,the
cut-offpointcanbedebated,butasweusethesamedefinition
betweenexperimentandCFDchangingthecut-offshouldnot
affectthecomparison.
4.
Results
and
discussion
4.1. Experimental
4.1.1. LDAverification
RanadeandJoshi(1990)comparedarangeimpelleroutflow
measurements,summarizedinthegraybandofFig.2. Our
datacoincideswiththeupperbound,weobserveasomewhat
slowerdeclineinUradwithradialposition.Thedifferencesmay
arisebecauseofvariationsinsignalprocessing,LDAsetup,
andsensitivityofUrad totheaxial position.Therapiddrop
atr=102mmisattributedtothevicinityofthebaffle.Agood
agreementinturbulentkineticenergyktisobservedcompared
toWuandPatterson(1989),MurthyandJoshi(2008)reported
ahigherktnearthetip,likelysincetheymeasuredallthree
velocitycomponentsandavoidedtheisotropyassumption.
4.1.2. inter-compartmentdynamics
Wenow focus on the dynamicsin the inter-compartment
region. LDA results forthis region are presented inFig. 3,
incombinationwithCFDforMRF–RKE.TheCFDresultswill
bediscussed later,inSection4.4.1. Itcanbeobservedthat
Uradexhibitsmirror-symmetryintheplaney=0(Fig.3A),with
somedeviationsnearthebaffle,andadecreasingmagnitude
ofUradnearingthestirrershaft.Thetypicalflowpatternsof
discturbinesmeantheflowmovestowardsy=0inthewall
region(convergentpattern),thisswitchestoflowawayfrom
y=0inthecenter(divergentpattern)betweenr/R=0.694and
r/R=0.566(Fig.3B).
Thetotalfluctuatingkineticenergykt*reachesamaximum
aty=0atevery radialposition,exceptclosest tothebaffle
(Fig.3C).Thispeakinkt* couldoriginatebothfromMIsand
from turbulent kineticenergygeneratedbythe collisionof
theflowloopsofthetop-andbottomimpeller.Todiscriminate
betweenthesepossibleorigins,thepresenceoflow-frequency
periodiccomponentsintheLDAsignalwasassessedby
gen-eratingtheauto-correlationfunctionˆofthesignalandfitting
withEq.(3).ThefitsarevisualizedinFig.4,withthefit
param-etersgiveninTable5.
As noted inSection 3.2.1, two frequencies(f1.1 and f2.1)
and their first harmonics (f1.2 and f2.2) are fitted. This
givesf1.1/N=0.020±0.001andf2.1/N=0.061±0.003axially,and
f1.1/N=0.015±0.007 and f2.1/N=0.052±0.015 radially. These
valuesagreewell withthe precessingvortex and jet
insta-bility frequenciesas reportedbyPagliantietal. (2008).The
totalMIkineticenergy,kMI,fit,reachesamaximumaty=0atall
radialpositions.Anoverviewofallfittingresultsisprovided
inTable5.
The radial oscillations are comparatively weak:
Cn,rad<0.075 at all locations except at r/R=0.694. The
resultatr/R=0.694isafittingerror,nostrongperiodic
com-ponentsareobservableinFig.4,andthedampingcoefficients
an (Eq.(3))are attheconstraintvaluean=0.33,indicatinga
poor fit.Axial oscillationshavea highermagnitude; Cn,ax
variesfrom0.1attheshaftto0.22atr/R=0.694.Closetothe
shaftf1.1 istwicethestrength ofitsharmonicf1.2/N=0.040.
For r/R>0.358, f1.2/N is by far the stronger contribution.
Overall, f2.1/N=0.061hasthe highestmagnitude,exceeding
f1.2/Nbyafactor1.5−2,exceptatr/R=0.2. Thenoise
com-ponent is around 1–6% at all points, except r/R=0.765.
This is consistent at all axial locations; the correlation
function shows a much more rapiddecay at this position
(not shown). Wehave notbeen able to pinpoint the exact
reasonforthishighernoiselevel;it couldbethatthe
pres-ence ofthe baffle leadsto relatively dominantsmall-scale
turbulence, meaning the data-rate was insufficient.
Fur-thermore, imperfections in refractive-index matching are
Fig.3–Axialprofiles(baffleplane)ofA:Urad,B:Uax,C:kt,comparingLDAresults(symbols)withCFDdata(lines)for
MRF–RKE,360◦domain.ByconstructionkMI=0forMRFsimulations,hencekt=kt*.Infig.C,theblackrectanglesrepresent
thetotalkineticenergykt*,thebluediamondstheturbulentkineticenergykt,andthereddiamondstheMIenergykMI.
Lines:CFDresultsatdifferentmeshdensities(dotted:2IF−C,dash-dot:2IF−M,dashed:2IF−F,solid:2IF−SF).Row(D) showsthetangentialaverageoftacquiredfromtheCFDsimulations.(Forinterpretationofthereferencestocolorinthis
legend,thereaderisreferredtothewebversionofthearticle.)
baffle may lead to higher background noise by reflection
effects.
The frequencies fitted by Eq. (3) contribute 8–20% to
thefluctuating kineticenergyaty=0.MIs havealess well
definedfrequencythanstirrerbladepassages,whichsuggests
auto-correlation fitting under-estimates their contribution
byexcessivedampening and/ormissingout additional
fre-quencies. This makes it likely kMI,fit is underestimated by
Eq. (5). Using the spectral density function and
determin-ing kMI,SI usingEq. (6), the kineticenergycontained inthe
rangef/N<0.1liesbetween28and49% aty=0,asreported
inTable6.However,this methoddoes notdiscriminateMI
frequenciesfromotherlarge-scalecontributions,andis
con-sidered over-estimative. The non-discriminatory nature of
the spectral density method does make comparison with
LES spectra more straightforward (part II), as LES
simula-tions showed less well-defined frequencies. Hence, the MI
energies reported in the CFD comparison come from the
spectral density method, withthe notion that the true MI
energy islikely tolie in between the estimates ofEqs. (5)
and (6). Fig. 5shows the spectral density function of ˆuax.
Inagreementwiththeauto-correlationfitsdiscussed
previ-ously,aclearbi-modalpeakrepresentingf1.2/Nandf2.1/Nis
visible,withincreasingprominenceattheouterradial
posi-tions.
Since we are interested in inter-compartment mixing,
uax,MIisofparticularinterest;inthisdirectionthejet
insta-bilityfrequencycontributesmost.WhetherMIswillstrongly
affectmixingdependsonwhethertheobservedoscillations
representacross-flowbetweenthecompartments,orthe
seg-regationbetweenthecompartmentsmovingupanddownas
Fig.4–Fittedauto-correlationfunctionsofuax(top)andurad (bottom)at3radiallocations,intheplaney=0.Blue:LDAdata. Red:fittedfunctionwith4dampedcosines,at2basefrequenciesandtheirfirstharmonic.(Forinterpretationofthe referencestocolorinthislegend,thereaderisreferredtothewebversionofthearticle.)
Table5–RMS-velocity,noisefactor,MI-correctionfactor,oscillationfrequenciesandkineticenergyintheplaney=0at differentradialpositions,usingtheperiodic-fittingapproach(Eq.(5)).Thecoefficientcrepresentthecontributiontothe Reynoldsstress ¯uu,withci.jthejthharmonicoffrequencycomponenti.
r/R 0.200 0.358 0.435 0.512 0.566 0.694 0.765 urad(m/s) 0.149 0.135 0.127 0.115 0.107 0.091 0.083 rad 0.015 0.017 0.022 0.032 0.037 0.064 0.317 c1.1,rad 0 0.006 0.016 0.017 0.006 0.032 0.010 c1.2,rad 0.037 0.048 0.040 0.007 0.003 0.002 0.001 c2.1,rad 0 0.015 0.075 0.009 0.012 0.045 0.006 c2.2,rad 0.032 0.005 0.028 0.016 0.011 0.039 0.012 f1.1,rad/N 0.020 0.001 0.015 0.015 0.022 0.017 0.017 f2.1,rad/N 0.049 0.058 0.056 0.053 0.050 0.043 0.062 uax(m/s) 0.082 0.098 0.112 0.134 0.154 0.194 0.151 ax 0.063 0.054 0.045 0.041 0.033 0.027 0.192 c1.1,ax 0.049 0.022 0.029 0.031 0 0 0 c1.2,ax 0.022 0.005 0.052 0.095 0.064 0.070 0.055 c2.1,ax 0.027 0.044 0.043 0.102 0.113 0.135 0.073 c2.2,ax 0.003 0.007 0.009 0.014 0 0.011 0 f1.1,ax/N 0.020 0.017 0.021 0.020 0.021 0.020 0.020 f2.1,ax/N 0.053 0.062 0.065 0.062 0.062 0.061 0.061 103×k t,fit/U2tip 7.73 7.63 7.49 7.54 8.85 10.80 6.17 103×k MI,fit/U2tip 0.62 0.41 0.78 1.36 1.23 2.40 0.70
Table6–MI-correctionfactorandkineticenergyintheplaney=0atdifferentradialpositions,usingthespectral integrationapproach(Eq.(6)).TheRMSvelocityandnoisefactorareasinTable5.
r/R 0.200 0.2358 0.3435 0.512 0.566 0.694 0.765 (Ep/Et)rad(Eq.(6)) 0.409 0.380 0.337 0.267 0.222 0.140 0.029 (Ep/Et)ax(Eq.(6)) 0.164 0.140 0.200 0.288 0.391 0.546 0.318 103·kt,2/U2tip 5.38 5.57 5.91 6.32 6.60 6.77 4.67 103·k MI,2/Utip2 3.04 2.47 2.35 2.58 3.52 6.45 2.22
papersuggestcompartmentcross-flowpredominantlyoccurs
closetotheshaft.However,wefirstfocusontheperformance
ofMRFsimulationswithsteady-statehydrodynamics(frozen
flowfield),whichareinherentlyunabletocaptureMIsbytheir
steadystatenature.
4.2. CFD:multiplereferenceframes
4.2.1. CFDvalidation:theimpelleroutflow
Fig.6showsUrad,ktandinthetopimpelleroutflow(360◦,
flow-Fig.5–Spectraldensityfunctionsofˆuaxatvariousradialpositionsandy=0.Dashedline:S(f/N)∝(f/N)−5/3.Dash-dotline:
S(f/N)∝(f/N)−1.Dottedline:cut-offfrequencybetweenMIandturbulence.
Fig.6–TangentiallyaveragedprofilesofUrad,ktandinthetopimpelleroutflowcomparedwithphase-averagedLDAdata,
forMRFsimulations.Toprow:realizablek−,360◦domain.Bottomrow:RSM360◦domain.LinesrepresentCFDdata;dotted line:2IF−C,dash-dotline:2IF−M,dashedline:2IF−F,solidline:2IF−SF.Symbolsrepresentexperimentaldata.
Abbreviations:W.P.WuandPatterson(1989),M.J.MurthyandJoshi(2008),D.Y.DucciandYianneskis(2005).Theblue crossesintheUradplotrepresenttheupper-andlowerboundofthestudiesreviewedbyRanadeandJoshi(1990).(For
interpretationofthereferencestocolorinthislegend,thereaderisreferredtothewebversionofthearticle.)
fieldinMRF,kMI=0byconstruction,andkt=kt*.Theprofiles
weretangentiallyaveragedinorderaccountforallimpeller
angles. In our view, this is a fairer comparison with LDA
datathantheinstantaneousoutflowprofiledirectlyfromthe
impellertip,althoughit reducesagreementnearthe baffle
wheretheflowisnotrotation-dominated.Noexperimental
profilesforhavebeenreportedforthecurrentstudy,asthe
equipment didnothavesufficientresolutionforsuch
mea-surements.
BothturbulencemodelspredictUradandktwellwithinthe
LDAdatarange;ktdoesdeviateclosetotheimpeller,wherethe
isotropicturbulenceassumptionbreaksdown.Interestingly,
thenon-isotropicRSMpredictsaktdecreaseatthetip,which
isinbetteragreementwithourdataandWuandPatterson
Table7–Comparisonofdimensionlessmixingtimes95
fortheMRFmethod:modelsandmeshdependency. Methods:P=bottomprobe,C=coefficientofmixing.The workofHartmannisfollowedtosettheCoM-boundary.
Mesh RKE,bot/top RKE,CoM RSM,bot/top RSM,CoM
2IF-C 133.5/110.5 148.5 138.1/107.1 152.0 2IF-M 146.1/120.0 162.0 157.8/127.9 173.0 2IF-F 145.2/120.3 160.5 145.2/113.1 159.5 2IF-SF 180.9/151.4 200.0 n.m. 172.5 2IP-C 96.2/n.m. 106.5 113.0/n.m 129.0 2IP-M 117.1/95.0 129.0 166.4/n.m 161.5 2IP-F 135.9/113.1 152.0 176.5/n.m 161.0 n.m.=notmeasured.
andJoshi,2008),whileonlythelattermeasuredall
(fluctuat-ing)velocitycomponents.RSMyieldsasuperiorassessment
for; atthehighestmeshdensitythepeakdissipationrate
recorded byBaldiand Yianneskis(2004),Baldi et al.(2004),
DucciandYianneskis(2005)isreasonablycaptured.TheRKE
modelagreeswiththedimensional assessmentofWuand
Patterson(WuandPatterson,1989),likelyarisingfromthe
sim-ilarunderlyingassumptions.Exceptforatthebladetipwith
2IF−C,nomeshdependencyisobservedintheRKEoutflow
profiles.TheRSMmodelismoremeshsensitive,but
indepen-denceisreachedwith2IF−F.Furtherparametersrelatedto
meshdependencearediscussedinSupplementarymaterial
B.Basedonimpelleroutflow data, bothturbulencemodels
performsatisfactory. Due tosimilarity thebottomimpeller
profilesandprofilesforthe60◦meshareomitted.
4.3. Mixingtimes
Wequantify 95bytwomethods;aprobeatthesamelocation
asusedbyJahodaetal.(2007)fordirectcomparison,andthe CoMasusedbyHartmannetal.(2006)toquantifymixingin
theentirevolume.Lackingexperimentaldataforthelatter,it
isreportedformutualcomparisonbetweensimulationsonly.
Theresultsare reportedin Table7. Experimentally,Jahoda
etal.reportabottomprobemixingtime 95,bot=92.Fortheir
top-probenomixingtimeisexplicitlyreported, 95,top≈75is
estimatedfromtheirproberesponseprofiles.
Forthe360◦-meshes,theprobe-based 95forbothRKEand
RSMexceedstheexperimentalforbothprobessubstantially.
Theover-predictionincreaseswithincreasingmeshdensity,
with a factor 2 over-estimation at the finest mesh, while
stillbeingpossiblymesh-dependent.Weconsideritunlikely
suchadegreeofover-estimationisaresultofthe(invasive)
probe influencing the flow experimentally, which is
sup-portedbysimilarresultsbetweeninvasiveandnon-invasive
point measurementsin singleimpeller systems (Distelhoff
etal., 1997;Javedetal.,2006).Probereadingsmay be
loca-tiondependent,implyingtheconcludedover-estimationof 95
maybelocation-specific. Toexcludethis,additionalprobes
havebeenmonitoredinthe bottomcompartment, yielding
very similar 95, both forMRF (Supplementary material C)
andSM–LESsimulations(partIIofthisstudy).Thisislogical
consequenceofinter-compartmentmixingbeingtherate
lim-itingstep,leadingtoarelativelyhomogeneousconcentration
withinthebottomcompartment.Thesimilar 95between
sev-eralbottom-compartmentprobes,combinedwiththesimilar
degree of over-estimation of 95 by the top- and
bottom-probecomparedtoexperimentaldata,givesconfidencethe
Table8–Axialflow-rateQax(downward)throughthe mid-impellerplaney=0inL/s.Forthe60◦mesh,the valueismultipliedby6.
Mesh RKE RSM RKE RSM
360◦ 360◦ 60◦ 60◦
X-C 1.36 1.58 3.04 2.45
X-M 1.14 1.28 2.72 1.34
X-F 1.48 1.31 1.10 2.25
X-SF 0.74 1.46 – –
concludedover-estimationin 95 isnotanartifactofprobe
location.Hence,thelogicalconclusionisthattheoffsetin 95
resultsfrominherentlimitationsinthemodelingapproach. WedonotesomeoutliersintheCFDresults;bothforRKE andRSM,95deviatesfromthetrendfor2IF−F.Thisappears
to bea mesheffect, which is furtherdiscussed in Section
4.4.1.The60◦meshyieldsasystematicallylower 95withRKE.
Thiswasnotaconsequenceofthe6-foldsymmetry
introduc-ingadditionalinjectionpoints,butratherofincreasedaxial
exchangebetweenthecompartments(seeSection4.4.1).With
theRSMmodel,the60◦ meshyieldssimilarCoMfiguresas
inthe 360◦ mesh,butwhereas theCoM becomesconstant,
theprobevaluekeepsincreasing.Thisdoesindicatedifferent
mixingbehaviorwithinthecompartmentsfortheRSMmodel.
4.4. Inter-compartmentflow
4.4.1. Multiplereferenceframes
360◦ domain.Theprofiles forUrad and Uax compare
reason-ablywellwithexperimentaldatausingRKE(Fig.3AandB),
althoughbothareover-predictedclosetothebaffle,probably
duetotheuseofa2Dbafflegeometry.Thechangefrom
con-vergingtodivergingaxialflowiswellcaptured.Somemesh
differencescanbedetected,butasidefrom2IF−Fnomesh
performsparticularlybetterorworse.
For 2IF−F, the peak in Urad and inflection point of Uax
lieslightlybelowy=0.Thisindicatestheinter-compartment
planehasshiftedcomparedtoy=0.InTable8wereportthe
integraldownwardflowrateQaxaty=0.Thehigherthisvalue,
themoreconvectiveexchangebetweenthestirrer
compart-mentsbythemeanflow(bymassconservation,theupward
flowratewasalwaysequal).Qaxtypicallydecreaseswithmesh
density,but2IF−F−RKEisaclearoutlier,whichagreeswith
the shiftinUrad andUaxobservedinFig. 3.Taken together,
thisexplainswhy 95iscomparativelylowin2IF−F,despite
kt andtangentially-averagedturbulent viscosity¯t (Fig. 3D)
beingsimilartotheothercases:thehigherQaxpromotes
mix-ing betweenthe stirrercompartments.Neartheshaft,kt is
reasonablycaptured,butthisagreementalreadybreaksdown
forr/R=0.358.Thepeakinktaty=0,mostprominentinthe
convergingflowregion,isunderestimatedbyaboutafactor
two usingRKE; ontopofthat theenergycontainedinMIs
isnotcapturedbyvirtueofthesteadystatenatureofMRF
simulations.Thetangentialaveragesoft,Fig.3D,showlittle
differencebetweenthemeshes;thevaluefor2IF−SFnarrowly
exceedsthe others.Thisindicates¯t isnotresponsiblefor
theincreasing 95 uponincreasingmeshdensity.Itis,
how-ever,likelythat¯t isstructurallyunder-predictedduetothe
poorassessmentofkt,leadingtoastructuralover-prediction
of 95asQaxconvergestoalow,mesh-independentvalue.It
isthistrendofdecreasingQaxwithanincreasinggrid
reso-lutionthat causesthe sensitivity 95 tomeshdensity.Note
Fig.7–Axialprofiles(baffleplane)ofA:Urad,B:Uax,C:kt,comparingLDAresults(symbols)withCFDdata(lines)for
MRF–RSM,360◦domain.ByconstructionkMI=0forMRFsimulations,hencekt=kt*.Infig.C,theblackrectanglesrepresent
thetotalkineticenergykt*,thebluediamondstheturbulentkineticenergykt,andthereddiamondstheMIenergykMI.
Lines:CFDresultsatdifferentmeshdensities(dotted:2IF−C,dash-dot:2IF−M,dashed:2IF−F,solid:2IF−SF).Row(D) showsthetangentialaverageoftacquiredfromtheCFDsimulations.(Forinterpretationofthereferencestocolorinthis
legend,thereaderisreferredtothewebversionofthearticle.)
increasedexchangebetweenthecompartments,asthe
inter-compartmentplane and y=0 may notexactly overlap, but
theredoesappeartobeacorrelationbetweenthetwo.TheUrad
andUaxprofilesforRSM(Fig.7)areinlesseragreementwith
experimentaldata.For2IF−Mthevelocityprofilesfavorably
matchtheLDAmeasurements,butforboth2IF−Fand2IF−SF
theUradpeakshowsstrangeasymmetricbehavior.Asforthe
RKEcase,theskewedinter-compartmentplaneincreasesQax
inthey=0plane(Table8),explainingwhy 95isnearlysimilar
for2IF−M,2IF−Fand2IF−SFwithMRF–RSM.
60◦ domain,MRFIncontrasttothefulldomain,Urad and
Uaxarepoorlypredictedinthe60◦domain(Fig.8,rowsA,B).
ThemagnitudeofUradisstronglyunder-estimatedeverywhere
exceptclosetothebaffle,whilereversalfromconvergingto
divergingaxialflowisnotatallcapturedwithinthemeasured
region.ThepeakUradispredictedataloweraxialpositionthan
measuredexperimentally.TheMRF–RKEmodeldoesyielda
comparatively good estimation ofkt for the 2IP−SF mesh,
resultinginahighertangentiallyaveragedt(Fig.8,rowE).
Furthermore Qax is significantly higher in the 60◦ domain
(Table8).Combined,thehightandaxialflowrateexplainthe
low 95observedforRKEinthe60◦domain.
ComparedtoRKE,RSMshows muchmoremesh
depen-dency,withunexpectedUradpeaks.ktisreasonablyassessed
forr/R<0.5.Overall,tandQaxaresimilarbetweenthe360◦
and60◦domain,explainingthesimilarityin 95betweenthem.
Aswasthecasewiththefull-domain,thereseemstobeashift
inplanepositioncomparedtotheexperimentaldataforthe
Fig.8–Axialprofiles(baffleplane)ofA:Urad,B:Uax,C:kt,comparingLDAresults(symbols)withCFDdata(lines)for
MRF–RKE,60◦domain.ByconstructionkMI=0forMRFsimulations,hencekt=kt*.Infig.C,theblackrectanglesrepresentthe
totalkineticenergykt*,thebluediamondstheturbulentkineticenergykt,andthereddiamondstheMIenergykMI.Lines:
CFDresultsatdifferentmeshdensities(dotted:2IF−C,dash-dot:2IF−M,dashed:2IF−F,solid:2IF−SF).Row(D)showsthe tangentialaverageoftacquiredfromtheCFDsimulations.(Forinterpretationofthereferencestocolorinthislegend,the
readerisreferredtothewebversionofthearticle.)
5.
Concluding
remarks
Multi-Rushtonstirredtanksarepronetocompartment
forma-tionaroundtheimpellers, withinter-compartmenttransfer
becoming the rate limiting factor for mixing. The mixing
time 95 appears to depend strongly on 2 factors:
turbu-lent exchange between the compartments, and the effect
of macro-instabilities (MIs) on inter-compartment
dynam-ics. Using LDA, we find evidence ofMIs atf/N=0.020 and
f/N=0.062,whichareinagreementwithearlierreported
val-uesfortheprecessingvortexandjetinstability(Pagliantietal.,
2006,2008).Dependingonthemethodofestimation,theMIs
containatbetween15%(fromauto-correlation)and30%(from
spectralintegration)ofthetotalfluctuatingkineticenergykt*,
althoughclosetothe bafflethismay beashighas20–50%,
mostsignificantlyintheaxialdirection.Thismakesitlikely
thatMIshaveaconsiderableeffectonmasstransferbetween
thecompartments.
Multiple reference frame(MRF) CFD modeling has been
employedinpriorliteraturetomulti-impellerstirredtanks,
exploitingitssteady-statenaturetoreducecomputationtime.
Aliteratureanalysisrevealedthat(1) 95 typicallyincreased
withmeshdensity,(2)thesimulationsreportedinliterature
wereunlikelymesh-independent,and(3)mesh-independent
results are expected to over-predict 95 strongly. Our own
MRFsimulations,focusedontheinter-compartmentregion,
confirmstheincreasingover-predictionsof 95with
increas-ingmeshdensity.Asignificantunder-estimationofkt*inthe
inter-compartmentregionisobserved.Inpartthisisdueto
the complete absenceofmacro-instabilities insteady-state
energy kt is lower than the LDA measurements; this is
attributed to Reynolds-averaged MRF simulations
predict-ingalargelyshear-freeinter-compartmentplane,leadingto
anunder-predictioninthegenerationofkt bythecolliding
compartmentflow-loops. Theturbulent viscosity t in the
inter-compartmentregion,whichviatheturbulentSchmidt
number Sct controls turbulent mixing, is nearly
indepen-dent of mesh density. While the inherent over-prediction
of 95 may be attributedto the under-predicted kt*, via an
under-estimatedt),theincreasing 95withmeshdensityhas
anotherorigin.Itwasfoundthat, despitetheradialflowin
theinter-compartmentplane,theaxialflowrateQaxthrough
thehorizontalplaney=0isnotexactly0,inducingsome
con-vectiveaxialmixing.ThisQaxreducedwithincreasingmesh
resolution,corresponding totheincreasing 95.There were
someoutlierswhereacomparativelylow 95andhighQaxwas
observed; further scrutiny revealed the inter-compartment
planewasslightlytiltedinter-compartmentplanewas
pre-dicted here. Thisimplies the presenceofMIs allowssome
ambiguityintheinter-compartmentplanepositionpredicted
insteady-stateMRF,whichcanmakethesimulationshighly
sensitiveto themesh layout andinitial conditions.In this
study we exclusively made use of structured hexahedral
meshes.Withtetrahedralmeshes,ahigherdegreeof
numeri-caldiffusionmayalsolower 95.InthelimitofNc→∞,strong
over-predictionsof 95 are still expected, regardless ofthe
meshtype.Usinga60◦sectionofthedomainyieldedgoodflow
agreementintheimpellerregion,butverypooragreementin
theinter-compartmentregion.Typically,Qaxwashigher
yield-inglower 95,butnotasaresultofbetterflowprediction.More
likely,forcingtheinherentlyunstablesegregationplaneinto
6-fold periodicity introduces artificial axial exchange. This
means the assumptionofperiodicityshould beavoided in
multi-impellermixingstudies,notbecauseitartificially
multi-pliesfeedpoints,butbecauseitpoorlycapturesflowbetween
thecompartments.Toconclude, theMRFmodelisnotable
toassesstherelevantdynamicsformixinginmulti-Rushton
tanks,andquitelikelysoinothercaseswhere
compartmen-talizedflowisobserved.TuningofSctissometimesproposed
inliteraturetoimproveagreementin 95.Thisad-hoc
solu-tionislikelytonegativelyimpactintra-compartmentmixing,
which does appear to be assessed properly, judging from
single-impellermixingsimulations(Coroneoetal.,2011).
Fur-thermore,theexact Sct correctiondependsonthe strength
ofMIs, the degree ofunderestimation of t, and the
pre-dictedpositionoftheinter-compartmentplane,whichmay
all vary on a case-by-case basis, thereby losing predictive
value.Inapplicationswheretransientmethodsmixing
simu-lationsarecomputationallyunfeasible,andcrudeagreement
inmixingbehaviorsuffices,Scttuningmaybetolerable,but
thelimitationsoftheMRFmodelmustbekeptinmind.For
detailedmixingstudies,suchasthequantitativeassessment
ofimpellerconfigurationsinmulti-impellermixingtanks,the
observedlimitationsmaketheMRFmodelwithafrozen
flow-fieldinherentlyunsuitable.Inthesecondpartofthiswork,we
willdiscusshowslidingmeshRANSandLESperforminthis
respect.
Acknowledgements
Wewish to thank our colleagues at ECUST, Shanghai and
DSMSinochem pharmaceuticals forourongoing
collabora-tion.WearegratefultoProf.H.vandenAkker,Dr.S.Kenjeres,
Dr. L. Portela and Prof.H.J. Noorman for their advice and
insights.CHacknowledgesallstudentsfollowingMSccourse
CH3421, who contributed tothe simulation resultsas part
oftheircoursereports.EvertWagner,JosThieme,Christiaan
Schinkel,StefantenHagen,YoupvanGoozenandRuudvan
Tolaregratefullyacknowledgedfortheirtechnicalassistance.
Thisworkhasbeenconductedwithinamulti-partyresearch
project,betweenDSM-SinochemPharmaceuticals,TUDelft,
East ChinaUniversity ofScienceandTechnology and
Guo-jia,subsidizedbyNWOandMoST(NWO-MoSTJointprogram
2013DFG32630).Allsponsorsaregratefullyacknowledged.
Appendix
A.
Supplementary
data
Supplementary data associated with this
arti-cle can be found, in the online version, at
https://doi.org/10.1016/j.cherd.2018.06.005.
References
Ammar,M.,Chtourou,W.,Driss,Z.,Abid,M.,2011.Numerical
investigationofturbulentflowgeneratedinbaffledstirred
vesselsequippedwiththreedifferentturbinesinoneand
two-stagesystem.Energy36(8),5081–5093,
http://dx.doi.org/10.1016/j.energy.2011.06.002.
Baldi,S.,Yianneskis,M.,2004.Onthequantificationofenergy
dissipationintheimpellerstreamofastirredvesselfrom
fluctuatingvelocitygradientmeasurements.Chem.Eng.Sci.
59(13),2659–2671,http://dx.doi.org/10.1016/j.ces.2004.03.021.
Baldi,S.,Ducci,A.,Yianneskis,M.,2004.Determinationof
dissipationrateinstirredvesselsthroughdirectmeasurement
offluctuatingvelocitygradients.Chem.Eng.Technol.27(3),
275–281,http://dx.doi.org/10.1002/ceat.200401979.
Brucato,A.,Ciofalo,M.,Grisafi,F.,Micale,G.,1998.Numerical
predictionofflowfieldsinbaffledstirredvessels:a
comparisonofalternativemodellingapproaches.Chem.Eng.
Sci.53(21),3653–3684,
http://dx.doi.org/10.1016/S0009-2509(98)00149-3.
Bujalski,W.,Jaworski,Z.,Nienow,A.,2002a.CFDstudyof
homogenizationwithdualRushtonturbines-comparison
withexperimentalresults.Chem.Eng.Res.Des.80(1),97–104,
http://dx.doi.org/10.1205/026387602753393402.
Bujalski,J.,Jaworski,Z.,Bujalski,W.,Nienow,A.,2002b.The
influenceoftheadditionpositionofatraceronCFD
simulatedmixingtimesinavesselagitatedbyaRushton
turbine.Chem.Eng.Res.Des.80(8),824–831,
http://dx.doi.org/10.1205/026387602321143354.
Bujalski,W.,Jaworski,Z.,Nienow,A.,2002c.CFDstudyof
homogenizationwithdualRushtonturbines–comparison
withexperimentalresults:PartII.Themultiplereference
frame.Chem.Eng.Res.Des.80(1),97–104,
http://dx.doi.org/10.1205/026387602753393402.
Coroneo,M.,Montante,G.,Paglianti,A.,Magelli,F.,2011.CFD
predictionoffluidflowandmixinginstirredtanks:numerical
issuesabouttheRANSsimulations.Comput.Chem.Eng.35
(10),1959–1968,
http://dx.doi.org/10.1016/j.compchemeng.2010.12.007.
Delafosse,A.,Collignon,M.-L.,Calvo,S.,Delvigne,F.,Crine,M.,
Thonart,P.,Toye,D.,2014.CFD-basedcompartmentmodelfor
descriptionofmixinginbioreactors.Chem.Eng.Sci.106, 76–85.
Derksen,J.,Doelman,M.,VandenAkker,H.,1999.
Three-dimensionalLDAmeasurementsintheimpellerregion
ofaturbulentlystirredtank.Exp.Fluids27(6),522–532,
http://dx.doi.org/10.1007/s003480050376.
Distelhoff,M.F.W.,Marquis,A.J.,Nouri,J.M.,Whitelaw,J.H.,1997.
Scalarmixingmeasurementsinbatchoperatedstirredtanks.
Can.J.Chem.Eng.75(4),641–652,
!["Benedykt Dybowski", Kazimierz Stołyhwo, Wrocław 1957 : [recenzja]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACH5BAEAAAAALAAAAAABAAEAAAICRAEAOw==)