Eulerian modeling of inertial and diffusional aerosol deposition in bent pipes

University of Groningen

Eulerian modeling of inertial and diffusional aerosol deposition in bent pipes

Frederix, E. M. A.; Kuczaj, A. K.; Nordlund, M.; Veldman, A. E. P.; Geurts, B. J.

Published in:

Computers & fluids

DOI:

10.1016/j.compfluid.2017.09.018

IMPORTANT NOTE: You are advised to consult the publisher's version (publisher's PDF) if you wish to cite from

it. Please check the document version below.

Document Version

Publisher's PDF, also known as Version of record

Publication date:

2017

Link to publication in University of Groningen/UMCG research database

Citation for published version (APA):

Frederix, E. M. A., Kuczaj, A. K., Nordlund, M., Veldman, A. E. P., & Geurts, B. J. (2017). Eulerian

modeling of inertial and diffusional aerosol deposition in bent pipes. Computers & fluids, 159, 217-231.

https://doi.org/10.1016/j.compfluid.2017.09.018

Copyright

Other than for strictly personal use, it is not permitted to download or to forward/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 (like Creative Commons).

Take-down policy

If you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediately and investigate your claim.

Downloaded from the University of Groningen/UMCG research database (Pure): http://www.rug.nl/research/portal. For technical reasons the number of authors shown on this cover page is limited to 10 maximum.

ContentslistsavailableatScienceDirect

Computers

and

Fluids

journalhomepage:www.elsevier.com/locate/compfluid

Eulerian

modeling

of

inertial

and

diffusional

aerosol

deposition

in

bent

pipes

E.M.A.

Frederix

a,∗

_,

_A.K.

_Kuczaj

d,a

_,

_M.

_Nordlund

d

_,

_A.E.P.

_Veldman

c,a

_,

_B.J.

_Geurts

a,b

a Multiscale Modeling and Simulation, Faculty EEMCS, University of Twente, P.O. Box 217, Enschede 7500 AE, The Netherlands

b Anisotropic Turbulence, Fluid Dynamics Laboratory, Faculty of Applied Physics, Eindhoven University of Technology, P.O. Box 513, Eindhoven 5600 MB, The Netherlands

c Institute of Mathematics and Computing Science, University of Groningen, P.O. Box 407, Groningen 9700 AK, The Netherlands d Philip Morris International R&D, Philip Morris Products S.A., Quai Jeanrenaud 5, Neuchatel CH-20 0 0, Switzerland

a

r

t

i

c

l

e

i

n

f

o

Article history: Received 19 August 2016 Revised 31 May 2017 Accepted 23 September 2017 Available online 25 September 2017 Keywords: Aerosol Deposition Drift Diffusion Transport Bent pipe

Droplet size distribution

a

b

s

t

r

a

c

t

ThispaperpresentsasectionalEulerianaerosolmodelforsize-dependentdropletdepositionatwallsof thedomain,drivenbybothdiffusionandinertia.Themodelisbasedontheinternallymixedassumption andemploystheformulationforcompressibleaerosols.Itisvalidatedinabentpipegeometry against modelsandexperimentalandnumericaldatafromliterature.Goodagreementisfoundinboththe diffu-sionandinertialdepositionregimes.Toimprovetheoverpredictionofinertialdepositionbyaboundary treatmentthatadoptszero-gradientdropletwallvelocity,weuseacorrectedwallvelocity,basedonan analyticalsolutionofthedropletmotionnearthewall.Inthebentpipesettingthecorrectedwallvelocity isfoundtoreducetheoverpredictionofdepositionandislesssensitivetogridrefinement.Wealsoshow thatrefiningthecomputational meshnearthe pipewallimprovesthepredicted deposition efficiency, significantly.Finally,wepresentaparameterstudyvaryingtheReynoldsnumberandthebendcurvature. Thedeposition efficiencycurveis recordedfordropletdiametersrangingfromthenanometer scaleto beyondthemicrometerscale,whichisauniquecontributionofthispaper. Thecompletesizerangeis simulatedinonlyonesimulation,duetothesectional approach.Inthediffusion-dominated regimean increaseinReynoldsnumberleadstoagradualenhancementofdeposition.Intheinertialregime,where dropletdriftdominatesdeposition,amuchstrongerdependenceontheReynoldsnumberisfound.The dependenceofthedepositiononthebendcurvatureislesspronounced.Theresultsshowninthispaper establishtheroleofEulerian simulationinpredictingdepositionofaerosoldropletsand areusefulfor understandingsize-dependentaerosoldepositioninothermorecomplexconfinedgeometries.

© 2017TheAuthors.PublishedbyElsevierLtd. ThisisanopenaccessarticleundertheCCBY-NC-NDlicense. (http://creativecommons.org/licenses/by-nc-nd/4.0/)

1. Introduction

Two of the main mechanisms of aerosol deposition are iner-tial impaction anddiffusion [1]. Bothprocesses are strongly size dependent; smallaerosoldroplets depositdue tohighdiffusivity, large droplets due to large momentum and intermediately-sized droplets deposit more scarcely. In a setting where the Reynolds number issufficiently highso thatgravitational settling becomes negligible,thisleadstothewell-knowndepositionefficiencycurve which has a characteristic ‘V’shape. This was observed in many kindsofgeometrythatinvolveaerosoldeposition,e.g.,in

respira-tory flow [2–4] or flow around an object [5,6]. The exact shape

∗ _{Corresponding author.}

E-mail address: e.m.a.frederix@utwente.nl (E.M.A. Frederix).

ofthisdepositioncurvecharacterizesthefiltrationefficiencyofan object orgeometry, and isvery usefulfor understanding aerosol deposition.

A common way to study aerosol deposition is to consider aerosolflow through a bent pipe.The bent pipe geometryoffers a simple setting inwhich the mechanisms behind aerosol depo-sitioncanbe systematicallystudied. Infact,thebent pipecanbe usedasahighlyidealizedmouth–throatmodeltoemulateaerosol dropletdepositioninthehumanairways,see[7].Bystudyingthe bentpipeaqualitativeimpressioncanbeformedofboththeflow andaerosoldepositionpatterns.

Thebentpipeproblemhasbeenstudiedby manyauthorsand thereforeoffersareliable andwell-understood pointofreference. Earlier theoretical works studied particle trajectories in the bent pipegivenanapproximateflow field[8–10].Morerecently,many https://doi.org/10.1016/j.compfluid.2017.09.018

authorshave publishedCFD simulations of particledeposition in bends using Lagrangian (e.g., [11–14]) or Eulerian methods (e.g.,

[15–20]). Others have studied aerosol deposition in pipe bends

experimentally, e.g., [7,12,21,22]. The seminal work of Pui et al.

[21] iskey in theliterature concerning bentpipe aerosol deposi-tionefficiencyandaclearsourceforvalidation.Forexample,Pilou etal.[15] found good agreementforReynolds numbers Re ₌100 and Re =1000andVasquezetal. [20] found good agreementfor

Re =1000whileanoverpredictionofthedepositionefficiencyfor

Re ₌100wasobserved.

Most bent pipe studies focused on the ‘inertial deposition regime’,lookingataerosoldropletsorparticleswithaStokes num-bertypicallylarger than0.01.Forthesedropletsitis theirinertia thatleadstoa collisionwithageometrywall.However, asnoted before, sufficiently small droplets may also deposit by diffusion. Infact,inmany applicationstheaerosoldroplet sizeissuch that droplet diffusion and inertia are two important effects, e.g., see

[2,23,24]. Inthis paper,we consider aerosoldeposition in abent

pipefordropletsizesrangingfromthenanometerscaletobeyond themicrometerscale.

Recently, wedevelopedanEulerian,sectional,internally-mixed aerosol model [25,26] capable of predicting the evolution of thedropletsize distributionundergoingnucleation, condensation, evaporationandcoagulation.Weformulatedacompressiblemodel inwhichthemixturedensityisconstitutedbyanumberof chem-ical species, either present as vapor or in the form of liquid droplets.Buildingonthatfoundation,inthispaperweextendthis modeltoinclude dropletdrift,diffusion andwall deposition.The mainobjectiveof thispaperisto presentthemodel andto vali-dateourEulerianapproachagainstdatafromliterature,inboththe diffusionandtheinertialregime.Moreover,westudyhow predic-tionsdepend onthe chosengrid andboundarytreatment forthe dropletvelocity.

The sectional Eulerian model retains a compressible formula-tioninwhichthemixturedensityiscomposedofbothvaporsand liquids, mitigating a passive scalar approach asis done in many otherworks,e.g.,[15,23].Thiscouplestheaerosolprocesses,such asdropletdriftanddiffusionbutalsonucleationandcondensation

(see[25,26])tothetransport equationsformass,momentumand

energy. This may be relevant in cases where mixture compress-ibilityisimportant,orwheretemperaturechangesarelarge. How-ever,alsoinsystemsnotexhibitingthesefeaturesthecompressible fluid framework is beneficial for obtaining general and accurate modelsasreliableconstitutiverelationscanbeformulated explic-itly.Incombinationwithapressure-basedapproach[27] this com-bines consistency inthe physical modelwith computational effi-ciency.Wedevelop a schemewhich,byconstruction,implements twoconstraintequations ensuring (1)that speciesmassfractions alwaysadd up tounityand(2)that thefirst momentofthesize distributionisalsoreflectedintheliquidmassconcentration solu-tion.

For large droplets we compare the predictions of our model against aforementioned experimental and numerical bent pipe studies. For small droplets we compare against the analytical straightpipediffusionaldepositionmodelofIngham[28].Inboth regimeswe findgoodagreement,providedsufficient spatial reso-lutionofthesolutionisadopted.

In this paper we present a detailed numerical study of our modelfor droplet diffusion, drift andsubsequent deposition. We studythetwocasespresentedbyPuietal.[21],forReynolds num-ber Re =100 and Re =1000, on five different meshes in which wecompareresultsobtainedwithorwithoutgridrefinementnear thewall.Weusetwoboundarytreatmentsforthedropletvelocity atthe wall,i.e., a‘zero-gradient’ boundarycondition, keepingthe dropletvelocity fromcell center tothe wall constant, anda cor-rectedboundaryconditionasproposedin[23],employingthe

an-alyticalsolutionofthedropletequationofmotionnearthewallin alinearizedflowfield.Thecorrectedboundaryconditionisshown to decrease the overestimation of the deposition efficiency, and generallyislessresolutionsensitive.Weshow thatthewall grid-refined meshes improve the predictions of the deposition curve significantly. We also present a parameter study for the depen-dence of the deposition efficiency on the Reynolds number and the bend curvature.An enhancement ofboth diffusion and iner-tialdeposition isshownforincreasing Reynoldsnumber whereas thedependenceonthebendcurvatureissmall.

In the model validation presented in this paper, we consider aerosol deposition in a bent pipe fordroplet sizesranging from the nanometer scale to beyond the micrometer scale. This enor-mous size-range is the unique feature ofour model: within one formulation the corresponding deposition efficiency curve span-ning the complete size domain is predicted. Moreover, the sec-tionalformulationspanningmanydecades indropletsizesallows forastraight-forwardextension toincludeaerosolprocesses such asnucleation,condensation,evaporationandcoagulationor break-up,aswasdonebefore(see[25,26,29]).Thecombinationofthese capabilitiesformsauniqueandquitecompleteaerosolmodel.

The layout of this paper is as follows. In Section 2 we will, starting from the equation ofmotion for the droplet size distri-bution,constructasetofequationsdescribinga compressible ‘in-ternally mixed’ [30] multi-speciesaerosolin an Eulerianway, in-cludinganewdrift fluxterm.Next,inSection 3,wewilladopta finite volume method and discretize the transport equations ac-cordingly. Again, special attention is paid to retaining the con-sistency among the equations, also in their discrete forms. Two boundary treatments forthe dropletvelocity at thewall are dis-cussed. In Section 4.1 we present the bend pipe geometry and

inSection 4.2 the fluid velocitysolution isvalidatedagainst data

from[15].Next,inSection4.3 thegridsensitivityofthesolutionis shownusingboththecorrected anduncorrected wall treatments. We validate the inertial and diffusion regime of the deposition curveinSections4.4 and4.5,respectively.Finally,inSection4.6 we presentaparameterstudy.

2. An internally mixed Eulerian aerosol model with droplet drift and diffusion

In this section, we will discuss the ‘internally mixed’ multi-species Eulerian model that we adopt for the description of the evolutionofanaerosolmixture.Subsequently,we extenditto in-corporatedriftfluxanddiffusionterms,basedonasize-dependent driftvelocityandBrowniandiffusivity.Theconsiderationstakenin arrivingatthedriftfluxmodelwillbediscussedhere.

2.1. Mass and droplet concentration transport equations

LetusconsideravolumeinwhichwehaveN species,present as vapor and liquid,where the liquid phase is contained in dis-perseddroplets.Withrespecttothetotalmassinthisvolume,the

j thspecieshasamassfraction Y jpresentasvaporand Z j present

asliquid.Bydefinition,wehave 

j

(

Yj+Zj

)

=1. (1)

Vapors areassumedtobe idealgasesandliquids areassumedto beincompressible.UsingAmagat’slaw[31],themixturedensity

ρ

canbe relatedtothespecies-specificvaporcompressibilityratios, species-specific liquiddensities, pressure andtemperature, giving anequationofstateintheformof

with

ψ

(p, T )the‘effective’compressibilityratio.FollowingFrederix etal.[26],weintroduce

ψ

−1

₍

_p,T

₎

₌ j Yj

ψ

 j

(

T

)

+p j Zj

ρ

 j

(

T

)

, (3)

with

ψ

_j

(

T

)

the vapor compressibility ratio and

ρ

_j

(

T

)

the liquid density,bothfor j -speciespureconstituentsofspecies j only.These quantitiesarebothtemperature-dependent.Since(2) togetherwith

(3) isbasedonAmagat’slawithastheimplicationthatvolumeis anextensivequantity,meaningthatthespecificvolumesare inde-pendentofmixturecomposition.

Fortheevolutionofthemixturedensity,assumingthatdroplets and vaporsare moving with the same velocity u ,the continuity equationholds:

∂

t

ρ

+

∇

·

(

ρ

u

)

=0 (4)

This can be expanded for all species in both phases, giving the transportequationsforthe j thspeciesvaporandliquidmass con-centrations:

∂

t

(

ρ

Yj

)

+

∇

·

(

ρ

Yju

)

=0 (5a)

∂

t

(

ρ

Zj

)

+

∇

·

(

ρ

Zju

)

=0. (5b)

Eq.(5b) provides informationabouttheevolutionoftheliquid mass concentration for each species, under the assumption that the overall liquid motionis equal to themixture motion, requir-ingtheliquiddropletstobesufficientlysmall.Wemayreadily re-lax thisassumption byincorporating thedropletsize distribution instead. In fact, the liquid mass is presentin the form of many dispersed droplets, suggesting that the liquid phase may also be describedbythedropletsizedistribution n (s , x , t )(s , x , t ),where s

isthemass ofadroplet.Thisisusefulforthe descriptionof pro-cessesthatdependonthesizeofadroplet,suchasdropletdriftor diffusion.Fortheevolutionofthedropletsizedistribution,we in-troducetheGeneralDynamicEquation(GDE),see[26,30]. Extend-ingthemodelasgivenby(5)toincludedropletmotionthatdiffers fromthemixturemotion,wereplacethemixturevelocity u inthe GDEwiththesize-dependent dropletvelocity v (s ), whichis writ-tenas

v

(

s

)

=u+u

(

s

)

, (6)

where u (s ) isthe liquiddrift velocity ofan s -sized droplet with respecttothemotionofthecarriergas, u .ThecorrespondingGDE for n (s , x , t )(s , x , t )canbeexpressedas

∂

tn

(

s,x,t

)(

s,x,t

)

+

∇

·

(

u

(

s,x,t

))

+

∇

·

(

u

(

s

)

n

(

s,x,t

)(

s,x,t

))

=

∇

·

(

D

(

s

)

∇

n

(

s,x,t

)(

s,x,t

))

, (7) where D

(

s

)

isthe droplet diffusivity, to be definedmomentarily.

Eq.(7) containstwoconvectivefluxes:one withrespect to u and a second one withrespect to u (s ). The latter flux term, aswell astheright-hand side diffusionterm,expresses thatdroplets can moveindependentlyofthemixture.Thisextensionrequirestoalso modify(5)toadequately reflecttheseadditionaldynamics. More-over,thecontinuityequationreflectingmassconservationmustbe augmentedwithanextraconservativedivergenceterm incorporat-ing thelocalappearance orremovalofdropletsby drift.Weturn tothistasknext.

Letus firstconsider the droplet sizedistribution for an inter-nallymixedaerosol,aswasdiscussedin[26].Thefirstmomentof thissizedistributionisrequiredtobeequaltothetotalmass con-centrationofdroplets,i.e.,

ρ

Z=

_∞

0

sn

(

s,x,t

)(

s,x,t

)

ds, (8)

where Z =jZ j.Thisequationimpliesthat

ρ

Z and n (s , x , t )(s , x , t )

aremutuallyconsistentanditallows ustorelatethedropletdrift tothe rateof changeof thetotal dropletmass concentration. By multiplying(7) by s ,andthentakingtheintegralfrom0to_∞in s , wefind

∂

t

(

ρ

Z

)

+

∇

·

(

ρ

uZ

)

+

∇

·

(

Z−1fZ

)

=0, (9)

wherewe haveintroduced theproduct Z −1Z (which isunity) into thedriftdivergenceterms,forlateruse,andwhere

f=fdrift_{− f}di f f ₍₁₀₎ and fdrift₌  _∞ 0 su

(

s

)

n

(

s,x,t

)(

s,x,t

)

ds (11a) fdiff₌  _∞ 0 sD

(

s

)

∇

n

(

s,x,t

)(

s,x,t

)

ds (11b) Theflux f canbeconsideredasthetotalfluxofliquid concen-trationmovingaway fromthemixturemotionduetodrift or dif-fusion.Recallingthat Z =jZ jandintroducingthisin(9),wecan

expand(9) in j -spaceas

∂

t

(

ρ

Zj

)

+

∇

·

(

ρ

uZj

)

+

∇

·

(

Z−1fZj

)

=0, (12)

Eq.(12) isconsistent with(9).Inthe caseofno driftandno dif-fusion,thethirdtermoftheleft-handsideiszeroandwerecover theoriginaltransportequation for Z jinwhichliquidisconvected

by u only,Eq.(5b).Satisfyingseparately(12) foreach j automati-cally satisfies(9) aswell.The product Z −1Z j asfoundin(12) can

beconsideredasthemassfractionofliquidspecies j withrespect tothetotal liquid massconcentration.Therefore,theflux Z −1f Z jin

(12) canbe interpretedasthetotalfluxof j-species liquid concen-trationdriftingawayfromthemixture.

The transport equation of j -species liquid massconcentration,

Eq. (12), is fully consistent with that of the size distribution,

Eq.(7),imposedbytheconsistencyrelation(8).Wesetoutto for-mulateanewaugmentedcontinuityequation,butforthiswemust still consider the transport of vapor concentration,

ρ

Y j. From a

physicalpointofviewweassumethatwhendropletsdrift,the vol-umethey‘vacate’isreplenishedbyacounterflowofvapormixture ofequalvolume.Inthisway,givenourcompressiblemixture equa-tionofstate(2) and(3) whichisbasedontheideathatvolumeis an extensive quantity, pressure remains uniform. We now intro-duceacompensatingvapordriftterm,toaccountforthis counter-flow. Without the explicit compensation of volume, the pressure wouldlocallychange duetodroplet drift,which leadsto a pres-suregradient-induced flow in u . This would affect the complete mixture.Let

ρ

v _{denote the localmean vapor densityand}

_ρ

 _the

localmeanliquiddensity.Thenthetotalvapormassconcentration driftflux,compensatingforthedropletmassdriftflux,isgivenby

h=

γ

f, (13)

with

γ

=

ρ

v_/

_ρ

_{,andthe j -species-specificcontributionistakenas}

Y −1h Y _j,i.e.,weighedwiththerelativecontributionofspecies j to thetotallocalvapormass,where Y =jY j.Subtractingthis

com-pensatingflux(itisinoppositedirectionrelativetotheliquidflux) fromtheleft-handsideofEq.(5a),wefind

∂

t

(

ρ

Yj

)

+

∇

·

(

ρ

uYj

)

−

∇

·

(

Y−1hYj

)

=0 (14)

Sincewehavechosentocompensatethesupplyorremovalof vol-ume asa resultof droplet drift by an equal vapor mixture vol-ume,inevitablythe localdensitywilldecreasewhen dropletsare movingaway,becausegenerally

γ

_{ 1.}Thisperturbationofdensity shouldbeaccountedforinthecontinuityequation,andinthe cor-respondingpressure equation. Byadding(12)–(14), andsumming over j wefind:

When

γ

=1, thecompensatingvaporvolume hasthe samemass asthe driftingdroplets, leadingto zeronetmasschangeby drift. Indeed,inthat casethe third termin Eq.(15) becomeszero, re-ducingittotheusualcontinuityequation.

Tosummarize,wenowestablishedasetoftransportequations formassconcentration, Eq.(15),liquidandvapor species-specific massconcentrations, Eqs. (12) and (14), and forthe droplet size distribution, Eq.(7), for the compressible Eulerian description of amulti-speciesinternallymixedaerosol.Consistencyamongthese equationsisenforcedbythemassfractionunityconstraint,Eq.(1), andthefirstmomentmassconcentrationconstraint,Eq.(8).

2.2. Drift velocity model

Intheprevioussectionweintroducedthedropletvelocity v (s ),

Eq.(6), differing fromthe mixture velocity u by a drift velocity. Wemustformulateamodelforthedescriptionof v (s ),asa func-tionofthedropletsize.Generally,withinthedriftfluxframework thereare two popular choices toestablish such a model. First,a so-calledlocal equilibrium approximation framework (see [6,32]) maybeadoptedinwhich v (s )maybefoundusinganalgebraic re-lationwhichstems fromthedropletequation ofmotioninwhich theaccelerationof thedroplet is assumedtobe equal to the lo-calcarrierfluid acceleration.Second, wemayapply thecomplete equation ofmotion ofa droplet retaining theacceleration ofthe droplet.Inthispaperweadoptthis‘fullmodel’.

Following Manninen et al.[32] the motion ofa single s -sized dropletcan,inaLagrangianway,bedescribedbyNewton’ssecond lawofmotion,i.e.,

zdtv=fD+fG, (16)

whereweonlyincludethedragforce f Dandthegravitationalforce

f Gasthedominantcontributionsfordropletswith

γ

 1[33].The

dragforcecanbewrittenas[34]

fD=−

1

2Ad

ρ

CD

|

u

|

u, (17)

withcross-sectional droplet area A d and dragcoefficient C D.The

gravityforceactingonthedropletisgivenby

fG=

(

ρ

−

ρ

v

)

Vdg, (18)

withdropletvolume V dandgravitationalaccelerationvector g . C D

isafunctionofthedropletReynoldsnumber,givenby Red=

d

ρ

v

|

u

|

μ

, (19)

withdroplet diameter d and vapor mixture viscosity

μ

.For suf-ficientlylow Re d (see [34] formore detail), C D takes theform of

Stokesdrag, CD,St=_Re24

d

. (20)

Eq.(16) constitutesapartial differentialequation (PDE)inthe Euleriancontext.Wemayexpressthematerialtimederivative as:

dtv=

∂

tv+

(

v·

∇)

v. (21)

UndertheassumptionthatStokesdragapplies,thePDEfor v ,given amass s andcorrespondingdiameter d ,becomes

∂

tv+

(

v·

∇

)

v=−

1

τ

(

v− u

)

+

(

1−

γ

)

g, (22)

withthedropletrelaxationtime

τ

=

ρ

d2

18

μ

. (23)

Eq.(22) can be solved numericallyfor v . Note that v = v

(

s,x ,t

)

,

i.e.,aEulerianvelocityfieldexpressingthetime-,space-and size-dependentdropletvelocity.

2.3. The Stokes–Einstein model for Brownian diffusion

Thediffusiveflux f diff_{dependsonthesize-dependentdiffusivity}

D

(

s

)

.ThisdiffusivityalsoappearsinthediffusiontermintheGDE. WeadopttheStokes–Einstein equation,providingamodelforthe Browniandiffusivityofasphericalbody.Itimplies[1]

D

(

s

)

= kBTCc

3

πμ

d, (24)

with k B theBoltzmannconstant,temperature T ,Cunningham

cor-rection factor C c, mixture viscosity

μ

and droplet diameter d =

d

(

s

)

relatedtothedropletmass s as s=1

6

ρ



π

d

3_{. (25)}

WecomputetheCunninghamcorrectionfactor,accountingfor sur-faceslipofsmalldroplets,as[1]

Cc=1+

λ

d



2.34+1.05exp



−0.39d

λ



, (26)

with

λ

the mean free path of the vapor mixture. For d 

λ

this functionbecomesunity.For d 

λ

, C cbecomesproportionalto d −1,

making the diffusivity (24) increase quadratically as d becomes smaller.

3. A sectional method retaining mass fractions and liquid mass

In the previous section, we motivated a set of equations for thecompressibleEuleriandescriptionofaninternallymixed multi-species aerosol. At the analytical levelwe showed that two con-sistency relations,one for themass fractions(1) andone forthe firstmomentofthedropletsizedistribution(8),weresatisfied.In thissection,wewillintroducethesectionalformulationto approx-imatea solutionforthesetofequations.Furthermore,atthe dis-cretelevelwe developourmethodsuch thatthe twoconsistency relations,(1) and(8),aresatisfiedbyconstruction.

3.1. The sectional formulation

Following Frederixetal.[25],we approximate thecontinuous size distributionfunction n (s , x , t )(s , x , t ) by asectional formula-tion,inwhichdropletsaredividedin_{P so-called}‘sections’labeled bytheir size s i.Mathematically, n (s , x , t )(s , x , t )isrepresentedby

n

(

s,x,t

)(

s,x,t

)

=

ρ



i

Mi

δ

(

s− si

)

, (27)

where N i=

ρ

M i isthe total number ofdroplets per unit volume

havingsize s _i,i.e.,

ρ

Mi=

yi+1

yi

n

(

s,x,t

)(

s,x,t

)

ds, (28) where y i and y i+1 arethelower andupperboundaryofsection i .

Theformulation(27) impliesthat,with

δ

(

s − si

)

havingasunitone

overmassand

ρ

M inumberper unitvolume,thesizedistribution

hasunitnumberperunitmassperunitvolume.

Bytakingtheintegralovertheinterval[y i, y i+1]of(7),and

us-ing(28),wefind

∂

t

(

ρ

Mi

)

+

∇

·

(

ρ

uMi

)

+

∇

·

(

ρ

uiMi

)

=0, (29)

wherethe driftvelocity u _i= u 

(

s i

)

isonly evaluatedat s i dueto

theDiracdeltafunctionrepresentationof n (s , x , t )(s , x , t )in(27). In terms of the sectional discretization the consistency relation

(8) becomes Z=

i

3.2. Finite volume discretization

OpenFOAM®_{offersacell-centeredcollocatedfinitevolume(FV)}

framework,inwhichwesolveoursystemofequations.Adetailed analysis of the integration ofthe Pressure-Implicit withSplitting ofOperatorsPISO methodinOpenFOAMwaspresentedin[27] for a two-moment description ofan aerosol, andextended to a sec-tional modelin[26].In thisframework, next,we include aerosol droplet drift. Inthe finite-volumemethod we consider a compu-tational volume V with F faces,labeledby faceindex1_{≤ f≤ F}and

A ftheoutward normalfacevectorwith

|

A f

|

=A f, i.e.,thesurface

areaofface f .Integrating(29) overthevolume V ,wefind, approx-imating

ρ

M iasconstantin V ,

∂

t

(

ρ

Mi

)

+D

φ

fMi, f

+D

φ

i, fMi, f

=0, (31) withfluxes

φ

f =

(

ρ

u

)

f· Af,

φ

idrift, f =

(

ρ

ui

)

f· Af and

φ

diff i, f =Di, f

∇

f

(

ρ

Mi

)

· Af (32) with

φ

 i, f=

φ

idrift, f −

φ

idiff, f (33)

and M _{i, f} =

(

M i

)

f, where (· )f denotes interpolation to the face f ,

whichweaddressmomentarily.Also,in(32) thereappearsaterm

∇

f(

ρ

M i). This is the numerical representation of the gradient of

ρ

M _iatface f .Notethatthisisgenerallydifferentthan (

∇

(

ρ

M _i))_f, i.e., thegradient of

ρ

M i computed atcell centers and only then

interpolatedtoface f .Finally,in(31) weadoptthenotation D

af

= 1 V  f af, (34)

summing a f overthefaces f enclosing V .Whatremainsisthe

dis-cretizationofthetimederivative.Weillustratetheremainingpart of the developmentof our method by adopting the notationally compact implicit Eulerscheme. We stress, however,that our ap-proachremainsapplicabletoanyothertimediscretizationscheme. UsingtheimplicitEulerscheme,(31) maybetimeintegratedfrom discrete time t m _{at time level m to t} m+1_{= t} m₊

_{t at time level}

m +1,as

(

ρ

Mi

)

m+1−

(

ρ

Mi

)

m

t =−D

φ

fMi, f

m+1 − D

φ

 i, fMi, f

m+1 . (35) ThisequationcanbesolvedusingthecompressiblePISOalgorithm, as shown in [26,27]. The choice of interpolation scheme for the computationof M i,fisessentialforthepreservationofpositivityof

M i.Forexample,thelinearinterpolationschemeisknownto

pro-duce oscillationsnearsharpgradients;itisnotmonotonicity pre-serving. Theupwindscheme,on theother hand,takes M i,fequal

to M i coming from the upwind direction, where the upwind

di-rectionisdeterminedbythesignoftheflux.Theupwindscheme isTVD(TotalVariationDiminishing).TVDschemeswereshownto havethemonotonicity property,i.e., thenumberoflocalextrema inthesolutiondoesnotincreaseandlocalminimaare nondecreas-ingandlocalmaximaarenonincreasing[35].Thisalsomeansthat when startingwith a positive solution, a TVDscheme preserves positivity.

Generally, a TVD schemedetermines its interpolation weights based ona limiter function. Thislimiter, inturn, isa function of thetransportedsolutionvariableitself,aswellasthefaceflux.In the case of(35), both divergenceterms may be easily combined into one flux,

φ

f+

φ

i, f, on which the limiter can be based. For

now,weleavethechoiceofinterpolationschemefor M i,f

undeter-mined (this choice willbe addressed inSection 3.4), butassume thatanappropriateinterpolationisusedthatguaranteespositivity of M m_i+1.

The set of P Eq. (35) forms the discrete counterpart of the droplet size distributiontransport Eq.(7).While for (7) we then enforcedtheanalyticalconstraint(8) tofindtheconsistent trans-portequationfor Z ,wewillnowfollowasimilarrouteforthe nu-mericalmodelstartingfrom(35),andapplyingthediscrete coun-terpartof(8),i.e.,condition(30),tothetransportequationfor M i.

Multiplying(35) by s _iandsummingover i ,wefind

(

ρ

Z

)

m+1₋

₍

_ρ

_Z

₎

m

t =−D

φ

fZf

m+1 − D

Z−1_f

φ

_fZf

m+1 + Jm+1_, (36) where

φ

 f=  i si

φ

_{i, f} Mi, f (37)

Byfollowingthe samesteps asin Section2, butnowatthe dis-crete level,we mayguarantee that the first moment consistency relation as expressed by (8) also holds discretely. Following this strategy,weexpand(36) in j -space,whichgives

(

ρ

Zj

)

m+1−

(

ρ

Zj

)

m

t =−D

φ

fZj, f

m+1 − D

Z−1_f

φ

_fZj, f

m+1 . (38) Toensure that thisequation isconsistent with (36) in thesense thatsumming(38) over j yields(36),werequirethat

Zf=



j

Zj, f. (39)

Thisrelationhasan importantnumericalconsequence,i.e., it im-plies that Z f must be computed from the individual Z j,f

inter-polants,andnot byfirstcomputing Z atcell centersandthen in-terpolatingthistothefaces. Whileathighspatialresolutions the differencesbetween



jZ j,fand(



jZ j)fmaybesmall,weadhereto

thealternativedefinitionof Z atthefaces,following(39).

Forconvenienceofimplementation,we combineboth convec-tivetermsin(38) intoasingleone,containingone fluxonwhich theinterpolationschemecanbebased.Thisgives

(

ρ

Zj

)

m+1−

(

ρ

Zj

)

m

t =−D



φ

f+Z−1f

φ

f



Zj, f

m+1 . (40)

The term betweensquare brackets on the right-hand side forms thefluxwithwhich Z jistransportedatface f .Atthelevelof

im-plementation,a difficulty withthisformisthat the fluxcontains

Z −1_{f ,}whichisundefinedfor Z _f→0.Fornumericalstability,the fol-lowingform, wherewemultiply theflux by Z _fanddivide Z _j,fby

Z f,improvesthis:

(

ρ

Zj

)

m+1−

(

ρ

Zj

)

m

t =−D



φ

fZ˜f+

φ

f



Z_{j, f} Zf+





m+1 , (41)

wherethe termbetweenthe square brackets isthe total flux on whichtheinterpolationschemeusedfor Z _j,fand Z _fisbased.This formhasthreeconsequences:

1.Thenumber



isaverysmallnumber.Incaseof Z _f=0as com-putedaccordingto(39),theintroductionof



preventsdivision by zero. For Z f→0,which, dueto positivity, alsoimpliesthat

Z _j,f→0 for all j , the term Z _{j, f/}

(

Z _f+



)

also goes to zero; no liquidispresentandnoliquidisconvected.

2.Theconvectedscalarbecomes Z −1_f Z j, f,whichisnon-linearin Z j.

We computeitexplicitly, basedonthelatestiterativesolution inthePISOalgorithm,see[27].

3.Insidethefluxthereappearsa Z fwhichreceivedanadditional

tilde in its notation. The reason forthis is that Z ˜f inside the

flux cannot be computedfrom a limiter interpolation scheme which is based on the flux, as this Z f is part of the flux

andtherefore denoteit asZ ˜f, indicating this. Dueto thefirst

momentconsistencyrelation(30) wecancomputeZ ˜_f as ˜

Zf=



i

siMi, f, (42)

whereeach M i,fiscomputedbyaTVDinterpolationscheme.

Eq. (41) is the discrete equivalentof(12). Bythe same token, wecandiscretizethe Y _j-Eq.(14) inananalogousway:

(

ρ

Yj

)

m+1−

(

ρ

Yj

)

m

t =−D



φ

fY˜f−

γ φ

f



Y_{j, f} Yf+





m+1 , (43)

where,asbefore,wehaveintroduced



forrobustness.Thetermin thesquare bracketsintheright-handsideisidentifiedastheflux withwhich Y j,f/Y f isconvected. Also,we introducedY ˜f forwhich

wederiveanexpressionnext.

Adding(41)and(43) toeachotherandsummingtheresultover

j givesthediscreteformofthecontinuityEq.(15),i.e.,

ρ

m+1

₍

_Y+Z

₎

m+1₋

_ρ

m

₍

_Y+Z

₎

m

t =−D

φ

f



_˜ Yf+Z˜f



+

φ

 f[1−

γ

]

m+1 , (44) providedthatwhensumming(43) over j wehave



jYj, f

Yf =

1, (45)

whichisanalogousto(39).Aftercomparisonofthisdiscreteform ofthecontinuity equation withits exact counterpart,(15), it be-comes clearthat we must requireY ˜f+Z ˜f=1. Toguarantee this,

wecomputeY ˜f as ˜ Yf=1− ˜Zf=1−  i ziMi, f. (46)

Thefinalformofthediscretedensityequationbecomes

ρ

m+1₋

_ρ

m

t =−D

φ

f+

φ

f[1−

γ

]

m+1 . (47)

3.3. Deposition boundary treatment

Whenweuseano-slipboundaryconditionforthemixture ve-locity u then u iszeroatthe wall.The onlytwo mechanisms al-lowingfordropletdepositiononsuchwallsarethen,seeEq.(7),a non-zero drift velocity u (s ), or a non-zero gradient of n (s , x , t ) perpendicular to the wall, i.e.,

∇

n (s , x , t )· n with n the general wallnormal.At the discretelevel, thistranslatesintosaying that

φ

drift

i, f =0toenableinertialdepositionor

φ

i, fdiff =0toenable

diffu-sionaldeposition,with f afaceatthewall,seeEq.(32).Wediscuss thesecontributionsnext.

3.3.1. Inertial deposition

First, we considerinertialdeposition. Weconsidera computa-tionalvolume V locatednexttoawallwithcellcenterposition x c

andtheface centeratthewall locatedat x f,seeFig.1 (left).The

cell-outwardunitnormalofthewallface f isdefinedas nf=

Af

Af

, (48)

andisshowninFig.1 (left).Atthecellcenterwehaveafluid ve-locity u and i thsection droplet velocity v _i=v

(

s _i

)

.Forvery small dropletsthese two vectorsbecome identical,andthe drift veloc-ity,i.e., u _i= v i− u tends towards zero.Giventhe no-slip

bound-aryconditionfor u atthewall,thedropletvelocityat x c,forvery

smalldroplets,will alsotendto zeroasthecomputationalgrid is refined. We now assume that droplets maintain their velocity v i

closetothewall.Thisleadstothefollowing uncorrected deposition boundary condition : BC

(

vi

)

atthewall →

∇

vi· n=0, ifvi· n>0 vi=0, otherwise (49a) BC

(

n

(

s,x,t

))

atthewall →

∇

n

(

s,x,t

)

· n=0, (49b) which,forthedriftvelocity,isamixedDirichlet–Neumann bound-aryconditionbasedonthecondition v _{i· n>}0.Ifthisconditionis satisfiedthenthedropletsclosetothewallhaveavelocity point-ing into the wall, andwe set the droplet velocity at the wall to be equal to this velocity (i.e., the face-normal gradient of v _i is zero). If thiscondition isnot satisfied,we set thedroplet veloc-ityatthewall tozero.The boundaryconditionfor n (s , x , t )is of typeNeumann,enforcingazerogradientsolutionperpendicularto thewall. Atthe discretelevel weimplementthe uncorrected de-positionboundaryconditionasfollows:

vi, f =max

(

vi,c· nf,0

)

nf (50a)

M_{i, f} =Mi,c (50b)

ρ

f =

ρ

c (50c)

wheresubscript(· )cdenotesthecellcenteredvaluecorresponding

the control volume of which f isa face. Eq. (50a) sets v i,f equal

totheface-normalcomponentof v _i,cifthiscomponentispositive, orotherwisetozero.Notethat(50a) alwayssetsavectoratface f

inlinewith n f,discardingthecomponentsof v corthogonal tothe

facenormal.Thiscanbe doneasthedrift dropletflux,expressed by Eq.(11), containsthe inner product with A f,making only the

wall-normalcomponentof v relevant.

The uncorrected deposition boundary condition can lead to a significantoverpredictionofthedropletvelocityatthewall[23].In practice,asdroplets travelfrom x c to x f,they decelerate.Longest

and Oldham[23] proposed a corrected deposition boundary treat- ment ,employingtheanalyticalsolutionofthedroplettrajectoryin betweencellcenterandfacecenter,assumingalineardecreaseof thewall-normalfluid velocityto zeroatthewall.Followingtheir work, the equation of motion for the wall-normal velocity of a dropletinsection i ,undergoingonlyStokesdrag,isgivenby[34]

dt

v

i, f =

τ

−1

(

uf−

v

i, f

)

, (51)

where

v

i, f =v i· nf, u f=u · nf and

τ

=

ρ

d2Cc

18

μ

(52)

the droplet relaxation time. The mixture velocity, at the interval [x c, x f],canbeassumedtolinearlydecreasetozero,i.e.,

u=−

_δ

x

f

uc, (53)

with x =

(

x − xf

)

· nf the wall-normal coordinate,

δ

f=

|

x f− xc

|

thecell-to-facedistanceand u c= u c· nf thewall-normal fluid

ve-locity at cell center. Inserting (53) into (51) yields a second or-der‘spring-damper’ordinary differentialequationforwall-normal droplet position x , to whichthe well-known solution isgiven by Longest and Oldham[23]. We can now determine the time t f at

which x _=−d/2_,i.e.,thefirsttimeatwhichthedropletintercepts the wall. At thistime the ‘interception dropletvelocity’ v i,f(t f) is

computed,andusedasthevaluefor v iatthewall:

Fig. 1. Left: a schematic overview of a computational volume with a wall face and the velocity boundary conditions. The cell center is located at x c and the wall face at x f . At the cell center the fluid has velocity u and the droplets of size s i velocity v i . The wall has outward unit normal n f . The wall-inward component of the droplet velocity,

vi · n f , is shown in gray. Right: M i as a function of x near the wall. For drift the concentration M i at the wall is approximated as cell-centered value M i,c . The gradient ∇f M i at the wall is approximated by the gray line working with M i = 0 at the wall.

Table 1

Chosen numerical schemes for the discretization of the indicated terms.

Term Scheme OpenFOAM

∂ t Implicit Euler [39] Euler

uf , T f Linear upwind scheme linearUpwind linear Yj,f , Z j,f , M i,f van Leer limiter [40] vanLeer

vi,f Upwind upwind

∇ (cell-center gradient) Central linear

∇2 (Laplacian) _{Central linear uncorrected} ∇f (face-normal gradients) Central uncorrected

3.3.2. Diffusional deposition

We modela perfectly-absorbingboundary toapproximate dif-fusional deposition. When droplets hit a wall they are absorbed instantly,andremovedfromthedomain.Diffusionaldeposition,as concludedbefore,isdrivenby

∇

n (s , x , t )_{· n}beingnon-zeroatthe wall.However,inthetwopreviouslydefineddepositionboundary conditionsweset

∇

n (s , x , t )· ntobezero:thenumericaltreatment ofthedriftdeposition preventsdiffusionaldeposition. Alsoatthe discretelevelweencounterthisproblem.Fig.1 (left) sketchesthe solutionof M inearthewall.Atthewall, M iattainsthevalueof M i

at thecell center, indeedshowingno gradient atthe wall. How-ever,duetotheperfectlyabsorbingboundarycondition,we com-putethe gradientof M iattheface in the diffusion flux asif M i is

zeroatthewall toinducediffusiontransport towardsthewall. If weretain(50c),thenthistranslatesintosaying

∇

f

(

ρ

Mi

)

=

ρ∇

fMi=−

ρ

c

Mi,c

xf− xc

at xf, (55)

whichisalsoschematicallyshowninFig.1 (left)bythegrayline. Even though the discrete implementations are conflicting, we

use (50a) or(54) forthe computation ofthe drift deposition flux

and(55) forthecomputationofthe diffusional deposition flux .This approach was also adopted in the drift flux models of Xi and Longest[36] andLongestandOldham[23],andprovedtobe suc-cessful.

3.4. Schemes and methods

We implementour modelin OpenFOAM, usingthe compress-iblePISO(PressureImplicitwithSplittingofOperator)algorithmas documentedin[27].Togetherwiththetransportequationsfor n (s ,

x , t )(discretizedintermsof M i), Y j, Z jand

ρ

aspresentedbefore,

we solve the Navier–Stokes equations forthe mixture velocity u

andtheenergyequationfortemperature T ,see[27].InthePISO al-gorithmthecontinuityEq.(15) isrewrittenintothepressure

equa-tion,usinganequationofstate,see[37,38].InOpenFOAMwemust selectsuitable spatialandtemporalschemesforthediscretization of the equations. Table 1 shows an overview of our choices. All discretizationschemeslisted inTable1 are thestandardschemes implementedinOpenFOAM.

4. Validation of inertial and diffusion aerosol deposition in bent pipes

In this section, we use our model and method to simulate aerosoldropletdepositioninbentpipegeometries.Forthe descrip-tionofthebentpipegeometry,Section 4.1,wecloselyfollowthe works by the authors in [10,15,20,21], all of which studied iner-tialdeposition intwogeometricallydifferent90° bentpipes,each operatedat adifferent Reynoldsnumber. In thissection, we will discussthebentpipegeometryandstudythenumericalflowand deposition solutions interms of temporal and spatial sensitivity. Moreover, in both the diffusional and inertial deposition regime we compareourresults againstdata fromliterature.Finally, bent pipedepositionforalargerangeofReynoldsandcurvatureratios isshown.

4.1. The bent pipe setup

Fig. 2 showsschematically the bent pipe geometry. In agree-mentwith[15],theinletofthebentsectionisextendedbya dis-tance D ,andthe outletbya distance2D ,with D thediameterof thepipe,asdepicted.We retainthischoice ingeometryto allow direct comparisonwith the results ofPilou et al.[15]. Following Puietal.[21] weset

R= r

Table 2

Simulation parameters for water droplets carried by air at room temperature and atmospheric conditions. The fluid velocity through the pipe is based on the Reynolds number.

Parameter Value Unit μ 1 . 81 × 10 −5 _m 2 /s ρ 1.1898 kg/m 3 ρ ₁₀ 3 _kg/m 3 T0 293.15 K D(Re = 100) 0 . 93 × 10 −3 _m D(Re = 10 0 0) 0 . 93 × 10 −3 _m r R D /2 _m U Reμ/( ρD ) m/s p0 10 5 Pa g 0 m/s 2

asthecurvatureratio,i.e.,theratioofbendradius r andpipe ra-dius R definedas R =D/ 2.TheReynoldsnumberisdefinedas Re=

ρ

UD

μ

, (57)

withbulk velocity U .TheStokes number,expressingtheratio be-tween droplet inertial time scale and the fluid convective time scale,isgivenby

St=

ρ

d2UCc

18

μ

R . (58)

where

ρ

istheliquiddropletdensityand d thedropletdiameter. Notethat theStokesnumberisbasedon R whereas theReynolds numberisbasedon D , inaccordanceto [21].ThePeclet number, introducedhereforlateruse,expressestheratioofconvectiveand diffusivedroplettransport.Itisdefinedas

Pe=UD

D , (59)

withDthesize-dependentdiffusivity,asgivenby(24).Inthebent pipegeometryanimportantquantityistheDeannumber,defined as

De=√Re

R. (60)

Theflow field inbends ofcircularcrosssection only dependson this number, expressing the ratio of the centrifugal and inertial forcestotheviscousforces.

Table2 listsallvaluesforthesimulationparameters,basedon

waterdropletsimmersedinairatroomtemperature T 0and

atmo-sphericpressure p 0.InallsimulationswesettheReynoldsnumber

byspecifying U ,whilekeeping

ρ

and

μ

constant andtakingpipe diameter D from [15]. Simulations are essentially governed only by Re and R  [41].We set gravityto zero,such that gravitational settlingof droplets isneglected. In [15,20] itwas concludedthat theeffectsofgravitywerefoundtoberelativelysmallforthetwo studiedcasesintheirwork.However,in[21] itwasconcludedthat forthe Re =100casetheorientationofthebentpipewithrespect togravityhadan effectonthedeposition resultsandthe experi-mentalreproducibility.Therefore,inthe Re =100casewewillalso includegravityinoursimulations fortheinertialdeposition vali-dation.Thiswillbeclearlyindicated.

Forthemesh,weapplyan‘O-typemultiblockstructure’aswas donein[15].Thepipecrosssectioncontainsan internalsquareof size D /5, to avoidvery small computational volumes atthe pipe centerline,aswould appearin apurely cylindrical mesh.The in-ternalsquare isextended tothe pipe wallby an enclosing cylin-dricalmesh,see Fig.2 (right). Towardsthe pipe wallwe apply a boundarylayergridrefinement,such that inradialdirectioncells atthewallareafactor5smallerthancellsattheinternalsquare.

Fig. 2. Left: schematic overview of the bent pipe geometry. The inlet extension has length D , the outlet extension length 2 D , with D the pipe diameter. The bent has radius r . Right: cross section of the O-type mesh, shown for Mesh C, with the in- ternal block of width D /5, number of internal square cells N s and number of radial cells N r . The color indicates the magnitude of the parabolic inlet velocity profile.

Table 3

Definition of the used meshes with N s the number of cells for a side of the internal square, N r the number of cells in radial direction, G the linear grading factor in radial direction (indicating the ratio between the smallest and largest cell, with the smallest cell at the wall) and N t the total number of computational cells for R  = 7 . The ‘no-bl’ suffix implies that no grid refinement to the wall is applied, i.e., G = 1 . Mesh N s N r G N t xc / D xw / D A-no-bl 5 9 1 41.820 0.040 0 0 0.035858 B-no-bl 7 12 1 104.720 0.02857 0.02988 C-no-bl 9 16 1 234.549 0.02222 0.02241 D-no-bl 11 20 1 459.459 0.01818 0.01792 E-no-bl 14 25 1 922.488 0.01429 0.01434 A 5 20 5 86.700 0.03568 0.007136 B 7 26 5 211.344 0.02752 0.005503 C 9 34 5 465.885 0.02108 0.004217 D 11 44 5 944.163 0.01632 0.003263 E 14 55 5 1.893.528 0.01307 0.002613

We introduce five differentgrids, A, B, C, Dand E,each being a refined version of its parent with Abeing the coarsest.Also, we definefivenon-gradedgrids,identicaltothefivegradedgridsbut without a boundarylayer grid refinement. Inthe pipe cross sec-tion,themeshisdefinedby N sand N r, beingthenumberofcells

spanningonesideoftheinternalsquareandthenumberofcellsin radial directionfromtheinternal squareto thepipewall, respec-tively.ThesetwonumbersarealsoindicatedinFig.2 (right).Along thepipeaxiscellshaveauniformspacing,chosen suchthat their widthinaxialdirectionisequaltothewidthincross-sectional di-rectionofacellwithintheinternalsquare.Thenumberofcellsin axialdirectionisa functionof R .Table3 givesthevaluesforNs,

NrandNt foreachgrid.Thetotalnumberofcells,Nt,isindicated

inTable3 for R =7.Themeshesarechosen suchthat N t roughly

doubles from one mesh to the next. This means that the typical celldimension

x decreasesroughlybyafactor21/3_{.Inthatsense,}

MeshEcontainscellsonlyafactor24/3_{smallerintypicalsizewith}

respecttoMeshA.

For the boundary conditions we apply the ones shown in

Table 4 forinlet, outlet andwall, inaddition to the uncorrected

orcorrecteddepositionboundarytreatmentforinertialdeposition andtheperfectlyabsorbingboundaryconditioninthecomputation ofthediffusionterms.Thefluxesfortheliquidconcentrations fol-low directlyfrom(11),andarecomputedexplicitly.At time t =0 we startfroma quiescentstatein whichnoaerosolispresentin thesystem.Weuseadistributionofsectionalrepresentativesizes

s iwhichisevenlyspacedinlogd -space,andrangesfrom

Table 4

List of boundary conditions for each variable and boundary patch. ZG stands for ‘zero gradient’.

Variable Inlet Outlet Wall Z (droplet) 10 −5 _{ZG ZG} Y (air) 1 − Z ZG ZG Mi M 0 ZG Mixed  p ZG 0 ZG T 293.15 ZG 293.15 u U  _{ZG 0} v u (Un)corrected ZG  As described in _{Section 3.3 .}

 Parabolic velocity profile with its maximum at 2 U such that the mean bulk flow is U .

size.Wesetaparabolicinletvelocityprofilewithitsmaximumat 2U ,forboth u and v i.Thisprofileisrampedstartingfroma

quies-centstate u ₌0overthetimeinterval t _∈[0_,

τ

˜]_,with ˜

τ

= L

U (61)

where L isthecenterlinelengthofthepipe.Thetime

τ

˜isthebulk flow-throughtime ofthesystem, aquantitywhichwe canuseto definethenon-dimensionalflow-throughtime t  as

t= t ˜

τ

. (62)

For t >

τ

˜ theparabolic inletfluid anddroplet velocity profile re-mainsconstantwithitsmaximumat2U .

4.2. Flow solution

In the work of Puiet al.[21] inertial deposition of aerosolis experimentallystudiedusingtwocases: Re =100, R _{=7and Re =}

1000, R =5.7.Thesetwocaseswerealsonumericallyinvestigated by Puiet al.[15], Tsai and Pui [10] and Vasquez etal. [20]. We compareourflowsolutionsagainst thewell-establishedresultsof Pilouetal.[15].

Fig. 2(left) shows pipe cross section lines C 1–C1 and C 2–C2.

Alongthesetwo lineswe compute thescaled axialvelocity mag-nitude,i.e., |u |/U .Fig. 3 showsthisquantity along C 1–C1 and C 2–

C 2 for the bend inlet (I), halfway the pipe bend (H) and at the

bend outlet(O) forthe Re ₌100 and Re ₌1000 cases, computed onthefiverefinedmeshes,A–E.Also,thedataofPilouetal.[15] is shown. For most meshes the results are very similar. Generally, we find goodagreement with[15], forall meshes. We see a no-table change in the solution towards the literature data forboth casesasthemeshisrefined.Thisisinagreementwiththe obser-vations of [20], where the Re =100 solution becomes effectively grid-independent beyond 106 _{cells while the Re =1000 solution}

becomesnearlygrid-independentfor3.3× 106_cells.

Inthispaper,we are inparticularfocusing ontheaccuracyof the solutionintermsof thediffusionandinertialdeposition. We willpresenttheassessmentofthegriddependenceofthose quan-tities inthe next section. We concludethat our flow predictions agreewellwithdatafromliterature.

4.3. Convergence of the deposition efficiency

Arelevantquantityinthestudyofaerosoldepositionisthe de- position efficiency

η

,definedasthe ratioof‘that whatdepositsin the system’ over ‘that what enters the system’. At thenumerical level,wedefinethedepositionefficiencyfordropletsection i as

η

i=  f∈W



i, f  f∈I



i, f , (63)

with W and I the setoffaces belongingtothewallandinlet, re-spectively,and



i, f=

φ

f+

φ

idrift, f

M_{i, f}+

φ

diff

i, f (64)

the total flux of s i-sized droplets crossing face f . Note that the

fluxescontain theface f surfaceareathrough theircorresponding definitions(11), such that summingthem over thesets W andI

represents the numerical equivalent of a surface integral. In the sectional formulation we can now, given a sectional distribution spanningsomespacein s ,computehow

η

dependsonthesize s .

Fig.4 (left)showsthelogarithmofthedepositionefficiencyas a functionof thelogarithm ofdroplet diameter d .In thislog-log space we recover a V-shaped deposition efficiency curve, where theleftsideoftheVisgovernedbydiffusionaldepositionandthe rightsidebyinertialdeposition.Forsmall d thedropletdiffusivity increases,increasingdiffusionaldepositioninturn.Forlarge d the dropletinertialtimescalebecomeslargeranddropletsdriftaway fromthecarriergastrajectory,increasinginertialdeposition.

Ourtime-resolvingalgorithm allowsto obtain

η

asatransient quantity. It will take some time before

η

has, given the initial condition of quiescence, converged to its steady state solution.

Fig. 4 (left) shows

η

as a function of d for six non-dimensional flow-through times t , computed using Mesh E. As t  increases, the

η

-curvesareseentoapproachasteadystate,forbothmeshes.

Fig.4 (right)showshowtheminimumofthe

η

-curve settlesasa functionofnon-dimensionalflow-throughtime t .Forthe Re =100 case, having a smaller Reynolds and Dean number, steady state isonlyattainedaround t =30, whereas the Re =1000 case con-vergesmorerapidly.Also,anincreaseingriddensityincreasesthe transienttime.ThismeansthatifthesolutionforMeshEis well-developed,soarethesolutionsonthecoarsermeshes,motivating thechoiceto onlyshow resultsforMeshEinFig.4.Forthe ‘no-bl’ meshes(see Table3) the transienttime amounts to lessthan 5 flow-through times(not shown in Fig.4). In the remainder of thispaperall presentedresultsaretakenat t =30,i.e.,the bulk flowwasallowedtopassthecompletelengthofthepipe30times, beforerecordingthesteady-statedepositionefficiency

η

.Thiswas foundadequate,asatthistimeallcurvesinFig.4 (right)have be-comereasonably‘flat’,indicatingsteadystate.

Fig.5 shows, for Re ₌100 and Re ₌1000_, thedeposition effi-ciency computed on all five meshes, with or without boundary layer, using the corrected and uncorrected wall velocity bound-arytreatment.The‘no-bl’ meshesproduceadeposition efficiency curve very dissimilar to that of the refined meshes witha large grid dependence, clearly indicating the relevance of a boundary layer mesh. Considering the solutions for Mesh A to E, we see that the deposition efficiency is only sensitive to the mesh den-sityinthebottomoftheV-shape,where

η

hasdropped toabout 10−4.Inthediffusion-drivenleftarmandinertia-drivenrightarm of the deposition efficiency curve the results are close to grid-independent. The uncorrected boundary condition gives, as ex-pected,asignificantoverpredictionofthedepositionefficiency,in particularfor the‘no-bl’ meshesandfor thecoarser wall-refined meshes.Forthecorrectedboundaryconditionthesolutionismuch lessgridsensitive,although forthe‘no-bl’meshesthedeposition efficiencyfor Re ₌1000aroundthebottomoftheVvariesnotably. Boththeuncorrectedandcorrectedboundaryconditionsappearto produce almost the same efficiency curve on the finest mesh E, bothfor Re ₌100and Re ₌1000.

Tostudythegriddependenceofthesolutionquantitativelywe introducethedepositioncurveintegral

e= 

log

η

dlogd, (65)

in which



is the domain in d shown inFig. 5. This relation is approximatednumericallybyasimpleRiemannsumoverall

sec-Fig. 3. Bent pipe steady state axial velocity profiles for Re = 100 (top row) and Re = 10 0 0 (bottom row), along pipe cross sections C 1 –C 1 (left column) and C 2 –C 2 (right

column), at the bend inlet (I), halfway the bend (H) and at the bend outlet (O), each horizontally offset by 2 for better visibility of the profiles. From gray to black (—) mesh A–E are shown. The cross sections C 1 –C 1 and C 2 –C 2 are defined in Fig. 2 (left). The data of Pilou et al. [15] ( •) is also shown.

Fig. 4. Left: the Re = 100 deposition efficiency ηfor a range of droplet diameters d , shown at six non-dimensional flow-through times (as indicated) for Mesh E. Right: the convergence of the scaled minimum of the logarithm of ηin time, for Re = 100 ( •) and Re = 10 0 0 ( ), computed on mesh E. The value for η∞ is taken at t  = 30 . Both

figures contain results using the uncorrected boundary treatment.

tions.Table5 gives e forthefivemeshes.Alsoshowninthistable isconvergencemeasure

e ,definedas:

e2= 



(

log

η

− log

ηparent

)

2_{dlogd (66)}

where

η

parent belongs to the coarser parent mesh.

e indicates

how,witheachrefinementstep,thesolutionchangesinlog-space. Generally,forfiner meshesit isshown that

e becomes smaller inboththe Re ₌100and Re ₌1000case,forbothboundary condi-tions,indicatingthat thesolutionconverges.Thisinformationcan alsobe visually distilled fromFig. 5 wherethe lines forMesh D andEareclosertoeachotherthanthoseofMeshAandB.

Fig.6 showsscaled depositioncurve integral e/e ref asfunction

of

x w asgivenby Table3.Thereferencevalue e ref istakenas e

fortheuncorrectedsimulationwithboundaryrefinementonMesh E.FromMeshAtoE

x wdecreasesroughlyonedecade.The

quan-tity e appears to convergence to unity as

x w becomes smaller,

sincethecurvesappeartobecomeflatterforsmall

x w.

Figs.5 and6 andTable5 indicatethatinmostpracticalcasesa computationperformedontherefinedMeshAusingthecorrected boundary conditiongivesa numericallyadequate estimate of the deposition efficiency in the diffusion, intermediate (i.e., the bot-tomofthedepositioncurve)andinertialregime,whichisalready within10%ofthedepositioncurvecomputedontherefinedMesh E.Havingestablishedanimpressionofthenumericalrobustnessof thesolutionintermsofthevelocityfieldanddepositionefficiency bothasafunctionofgriddensity,thequestionremainshow

phys-Fig. 5. The grid dependence of the η-curve for Re = 100 (left) and Re = 10 0 0 (right), for all meshes as defined by Table 3 , using the uncorrected boundary condition (top) and corrected boundary condition (bottom).

Table 5

Deposition curve integral e , Eq. (65) , and convergence measure e , Eq. (66) , for Re = 100 and Re = 10 0 0 , uncorrected and corrected deposition velocity and for all five meshes.

Uncorrected Corrected Re = 100 Re = 10 0 0 Re = 100 Re = 10 0 0 Mesh e e e e e e e e A −8.0321 −8.6685 −8.2704 −9.7996 B −8.1047 0.1003 −8.6045 0.1311 −8.2500 0.1027 −9.5519 0.4312 C −8.2084 0.1039 −8.7509 0.1261 −8.2710 0.1101 −9.5210 0.1603 D −8.3003 0.0862 −8.9282 0.1393 −8.3223 0.0632 −9.5393 0.1117 E −8.3676 0.0613 −9.0832 0.1235 −8.3749 0.0515 −9.5682 0.0822 A-no-bl −6.3214 −3.7847 −7.7941 −7.0371 B-no-bl −6.6076 0.2059 −5.0040 0.6980 −7.8413 0.1403 −7.8298 0.4130 C-no-bl −6.9363 0.2355 −6.0595 0.6248 −7.9229 0.1646 −8.4605 0.3314 D-no-bl −7.1788 0.1788 −6.6810 0.3861 −7.9759 0.1648 −8.8041 0.1944 E-no-bl −7.3997 0.1689 −7.1492 0.3029 −8.0362 0.1073 −9.0306 0.1420

Fig. 6. Deposition curve integral e , Eq. (65) , as function of the typical scaled wall grid cell size xw / D , for Re = 100 (left) and Re = 10 0 0 (right), with ( ◦, •) and without ( , ) boundary layer and with the corrected ( ◦, ) and uncorrected ( •, ) boundary treatment for the droplet wall velocity.