The effect of Stefan flow on Nusselt number and drag coefficient of spherical particles in non-isothermal gas flow

ContentslistsavailableatScienceDirect

International Journal of Multiphase Flow

journalhomepage:www.elsevier.com/locate/ijmulflow

The effect of Stefan flow on Nusselt number and drag coefficient of spherical particles in non-isothermal gas flow

Thamali R. Jayawickrama

^a,∗

, Nils Erland L. Haugen

^b,c

, Matthaus U. Babler

^d

, M.A. Chishty

^a

, Kentaro Umeki

^a

aEnergy Engineering, Div. Energy Science, Luleå University of Technology, Luleå 971 87, Sweden

bDepartment of Energy and Process Engineering, Norwegian University of Science and Technology, Kolbjørn Hejes vei 1 B, 7491 Trondheim, Norway

cDepartment of Thermal Energy, SINTEF Energy Research, Kolbjørn Hejes vei 1 A, 7491 Trondheim, Norway

dDepartment of Chemical Engineering, KTH Royal Institute of Technology, SE-10044 Stockholm, Sweden

a rt i c l e i nf o

Article history:

Received 17 August 2020 Revised 21 December 2020 Accepted 25 March 2021 Available online 29 March 2021 Keywords:

Drag coefficient Nusselt number Stefan flow Boundary layer multiphase reactive flow

a b s t r a c t

AStefanflowcanbegeneratedduringaphasechangeorreactionsofaparticleimmersedinafluid.This study investigates theeffectofStefanflowonthe exchangeofmomentum(dragcoefficient(CD))and heattransfer(Nusseltnumber(Nu))betweentheparticleandbulk-fluid.Fullyresolvedsimulationswere carriedoutforaflownearasphericalparticleimmersedinauniformbulkflow.Theimmersedboundary methodisusedforimplementingfluid-solidinteractionsandtheparticleisconsideredasastaticbound- arywithfixedboundaryconditions.Inanon-isothermalflow,thechangesinthermophysicalproperties attheboundarylayerplayedaroleinthevariationofCDandNubyaStefanflowfurther.Thepreviously developedmodelforthedragcoefficientofasphericalparticleinauniformisothermalflowwasmodi- fiedforauniformnon-isothermalflow.Themodelisdevelopedbasedonphysicalinterpretation.Anew modelisdevelopedfortheNusseltnumberforasphericalparticlewithauniformStefanflowcombining availablemodelsinliterature.ThemodelsarevalidatedforStefanReynoldsnumber−8Re_{s f,p}25and particleReynoldsnumberof2Re_f30ingasflow(i.e.Pr≈0.7).

© 2021TheAuthor(s).PublishedbyElsevierLtd.

ThisisanopenaccessarticleundertheCCBYlicense(http://creativecommons.org/licenses/by/4.0/)

1. Introduction

Particle-laden flows have many complexities due to e.g. flow separation, particle wakes, multi-particle effects, Stefan flow ef- fects and reactions. Such flows are associated with physical ef- fects that have a wide range of length and time scales. For ex- ample, the largest length scale in pulverized boilers (reactor) is O(^101m) ^{while the smallest physical scale (particle boundary} layer) isO(^10−5m)^{andthesmallestchemicalscaleis}O(^10−10m)^. Therefore, it is currently impossible to resolve all scales in any numerical setup of practical relevance. This gap can be bridged by developingmodels describing the effects occurring atsmaller scales (smallest physical/chemical scales). The smallest physical scales(O(^{10−5m))canbestudiedthroughdetailednumericalsim-} ulations. In contrast to experimental data, numerical simulations createavirtualenvironmentthat ismuch moreversatiletoeluci- date the relevanttransport phenomenaandthat can be used for

∗Corresponding author.

E-mail address: [email protected] (T.R. Jayawickrama).

developingmodels.Inthecurrentstudy,weinvestigatetheStefan floweffectsinparticle-ladenflowsusingnumericalsimulations.

AStefan flow iscreated when thereis a netflow of gas/fluid towardsorawayfromasolid surfacethat isreactingorundergo- ing aphase change(MurphyandShaddix,2003). Some examples are:evaporation,condensationandcombustionofdropletsaswell aspulverizedfuelcombustionandgasification.TheStefanflowcan affecttheexchangeofmass,momentumandheatbetweenthesur- faceandthebulk fluidinparticle-ladenflows.Models forNusselt number(Nu),Sherwoodnumber(Sh)andthedragcoefficient(C_D) areusedtocalculateheat,massandmomentumtransferbetween the particle and the fluid, respectively. However, this study will onlyconsidertheNusseltnumberandthedragcoefficient.

In the past, the Stefan flow effectwas considered for droplet evaporation and combustion (Renksizbulut and Yuen, 1983b;

1983a; Abramzon andSirignano, 1989; Harpole, 1981). Lately,an interest for the effect of the Stefan flow has emerged for coal combustion applications due to high reactive gas concentration inOxy-fuelcombustion(O₂/CO₂)comparedtoair-fuelcombustion (N₂/O₂).Theimportance ofStefanflow inOxy-fuelcombustionof coal isemphasized byYu etal.(2013). Accordingto them,aSte-

https://doi.org/10.1016/j.ijmultiphaseflow.2021.103650

0301-9322/© 2021 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ )

Nomenclature

RomanSymbols

Symbol Description(Units) A crosssectionarea(m²)

cp specificheatcapacity(Jkg⁻¹ K⁻¹) D diameteroftheparticle(m) F force(N)

h heattransfercoefficient(Wm⁻²K⁻¹)

−

→I identitymatrix(1)

L latentheatofevaporation(Jkg⁻¹)

−

→n unitnormalvector(1) p pressure(Pa)

R radius(m) S surfacearea(m²) T temperature(K)

t weightingfactor(between0to1)(-) U velocity(ms⁻¹)

−

→u velocityvector(ms⁻¹) V volume(m³)

GreekSymbols

δ

^{boundarylayerthickness(m)}

μ

^{viscosity(Pas)}

ρ

^{density(kgm−3)}

τ

^timescale(s)

λ

^thermalconductivity(Wm⁻¹K⁻¹) Subscripts

b boilingpoint(-) B boundarylayer(-)

∞ parameterscalculatedatthefar-fieldcondition(-) f parameters calculated atthe film condition (when

t=0.5)(-) l liquid(-)

s f withStefanflowconditions(-)

p parameterscalculatedattheparticlesurface(-) Dimensionlessnumbers

B_T Spaldingheattransfernumber(B_T=^cp(TbL^−T∞)) C_D Dragcoefficient(C_D= _0.5ρ^F_U2A)

Nu Nusseltnumber(Nu= ^hd_λ) Pe Pecletnumber(Pe=Re×Pr) Pr Prandtlnumber(Pr= ^cp_λ^μ) Re Reynoldsnumber(Re=^ρ_μ^UD)

fan flowhasa stronginfluenceonthemasstransfer rateinZone II conversion(kineticallyanddiffusioncontrolled)whiletheeffect isinsignificantinZoneIII(diffusioncontrolled)duringburnoutpe- riod.StillitisnotclearfromtheirresultswhenStefanflowcanbe neglected.

The main objectiveofthe currentpaperisto studythe effect of Stefan flow on Nusselt number and drag coefficient for non- isothermal conditions(i.e.whenthereisatemperaturedifference betweenparticleandgasfield).Eventhoughthemodelisgeneric andmeantto beapplicable fora varietyofconditions,it wasde- veloped andvalidated witha primary interest on entrained-flow biomassgasification.Assummarizedinthenextsection,weaimto fillagapinknowledgeandmodels,especiallyunderthepresence oflarge temperaturedifferences(i.e.>100K). Hereafter’tempera- ture difference (^{T)’ means the} temperaturedifference between the solid particle (sphere) surface and the far-field of the fluid.

Simulationsresolvingtheboundarylayerarecarriedoutforalam- inarflowsurroundingastaticsphericalparticle.Multi-component

effects were avoided forthe simplicityof work. The applicability of our modelfor the drag coefficient, developed in ourprevious work under isothermal conditions (Jayawickrama etal., 2019), is assessed andextended to non-isothermal conditions.In addition, a newmodel describing the effect ofStefan flow on the Nusselt numberisdeveloped.

2. Previousstudies

2.1. Nusseltnumberathightemperaturedifference

The Nusselt number (Nu=hD/

λ

^{) is usually expressed as a}

function of Reynolds number (Re=

ρ

^UD/

μ

^{) and Prandtl number}

(Pr=cp

μ

/

λ

^{). A Nusselt number formula that is applicable for}

both high and low temperature difference conditions is hard to find in the literature. Two popular models are the models of Whitaker(1972) and themodel of Ranz-MarshallRanz andMar- shall(1952).Theformerreadsas:

Nu=2+

(

⁰.4Re¹²+0.06Re²³

)

^Pr0.4

μ

∞

μ

p

¹4

, (1)

wherethermophysicalproperties(i.e.

λ

,

ρ

,

μ

,andcp) forthecal- culationofNusseltnumber,Reynoldsnumber,andPrandtlnumber are based on far-field conditions,

μ

∞ is the viscosityat far-field conditionand

μ

^{pistheviscosityatparticlesurfacecondition.The}

Ranz-MarshallmodelRanzandMarshall(1952)isgivenas:

Nu=2+0.6Re¹²Pr¹³, (2)

wherethermophysicalpropertiesatfilmconditionareusedtocal- culateNusseltnumber,ReynoldsnumberandPrandtlnumbers,in- stead of those at far-field conditions. Film condition is defined as the average between the far-field condition and the surface condition, i.e.T_f=(^T∞+Tp)/2 where T_∞ andTp are the far-field and surface temperatures, respectively. At low temperature dif- ferences and Reynolds numbers (≈0<Re<100), the Whitaker model(Eq.1)typicallygivespredictionsthatare closertotheac- tual values(NikrityukandMeyer, 2014), whilethe Ranz-Marshall model (Eq. (2)) can be applied for high temperature differences (1<Re<130)(Ellendtetal.,2018).

Thereare numerousworkson developingmodelsfortheNus- selt number associated with droplet evaporation. Evaporation at hightemperature differences requires consideration of the varia- tionofthermophysicalproperties,suchasthermalconductivity(

λ

⁾

and specific heat capacity (cp). This effect can be accounted for throughacorrectionfactorfortheNusseltnumber(Harpole,1981), orby introducinga referencetemperature(Naraslmhan andGau- vin,1967;Downingm,1966; YuenandChen,1978).Thereference temperatureisthencalculatedasfollows:

Tt=tT_∞+

(

^1−t

)

^{Tp, (3)}

wheret isweightfactor.

2.2. EffectofStefanflowonNusseltnumber

DifferentmodelsfortheNusseltnumberdevelopedforevapora- tionofsingledropletsaresummarizedbyZhifuetal.(2013).They have categorizedthe available models into theoretical, numerical andexperimental models.Accordingto their comparisons,all the modelsaredeviatingfromexperimentalresultswhentheevapora- tionrates arehigh. Therefore,they havedevelopeda modelwith acorrection factorthat isapplicableforhighevaporationratesas well.Inthismodel,theNusseltnumberisgivenas:

Nu_Zh=f_TNu, (4)

where

f_T=

(

^1+BTp

)

^{−23, (5)}

and

Nu=2+0.552Re¹²Pr¹³. (6)

HeretheSpaldingheattransfernumber(BTp)isdefinedas:

BTp=cp,p

(

^T∞−T_b

)

L , (7)

where Lis latent heat of evaporationandT_b isthe boiling point temperature. The Reynolds number is calculated based on prop- erties at theparticle surface, while the Prandtl numberis calcu- latedbasedonfar-fieldcondition.TheNusseltnumberiscalculated basedon propertiesattheparticlesurfacecondition.Itisnoticed that themodel ofZhifu etal. (2013) hasno explicit dependence ontheStefanflow.TheeffectofStefanflowisaccountedforindi- rectly throughtheevaporationrate,characterizedby theSpalding heattransfernumber.

Niazmand and Renksizbulut (2003) used the model devel- oped for droplet evaporation by Renksizbulut & Yuen (Nu_RY) Renksizbulut andYuen(1983a) forthegeneralized caseofa Ste- fanflow:

NuRY=2+0.57Re^1/2Pr^1/3

(

^1+BTf

)

^{0.7 (8)}

wheretheReynoldsnumberiscalculatedbasedonparticlesurface conditions, Prandtl numberis calculated based on film condition andtheSpaldingheattransfernumberisdefinedas:

B_Tf= PrRes f

Nu , (9)

where Re_{s f}=

ρ

^Us fD

μ

^{, (10)}

is the Reynoldsnumber basedon the Stefan flowvelocity (here- after, called Stefan Reynolds number). The variation of thermo- physical properties are neglected forNiazmand and Renksizbulut Niazmand andRenksizbulut(2003)andtheselectedrangeofSte- fan flows wasbasedon droplet evaporation(0.01 ^UU^{s f}_∞ 0.04).

Murphy&ShaddixMurphyandShaddix(2003)haveformulateda Nusselt number(Nu_M) correlation, forStefan flow in a quiescent environment.Assumingconstantproperties,theirexpressionreads as:

Nu_M=Nu

(

^PrRes f

)

/Nu

e⁽PrRes f)/Nu−1, (11) where Nu=2 istheNusselt numberina quiescentflowwithout Stefanflow.Recently,KestelKestel(2016)developedanewempir- ical model applicable for the convective flow environment based on his simulation data that gave better accuracy than the other available models.Inthismodel,whichisapplicable forRe<200, Re_{s f}<20and0.744<Pr<1.5,theNusseltnumber(Nu_K)isgiven as:

NuK=Nuexp

−0.54PrRe¹_{s f}^.126 Nu^1.052

, (12)

where

Nu=2+0.39Re^0.56Pr^0.45. (13) In Eq.12 and13,allpropertiesarecalculated basedontherefer- encetemperatureasdefinedinEq. (3)whentheweightfactoris t=0.9. Thismodel hasa largenumber offittingparameters and it doesagreebetterwithsimulationresults.However, itdoesnot necessarilyrepresentthephysicalphenomena.

Insummary,mostofthecurrentlyavailablemodelsforNusselt numberforparticleswithStefanflowinaconvectiveenvironment are empirical. Oneof thevery few theoretical models (ofEq.11)

(MurphyandShaddix,2003)developedfortheNusseltnumberof particles with Stefan flow is for a quiescent environment and is basedona constantpropertyassumption.Therefore,thereare no modelsforStefanflowinaconvectiveenvironmentbasedonphys- icalinterpretationwhileconsideringvariationofproperties.

2.3. Dragcoefficientsathightemperaturedifferences.

The drag coefficient is defined asC_D=F/(⁰.5

ρ

^U2A), where F is the drag force, A is the cross-sectional area of the particle,

ρ

is the density of the fluid and U is the velocity difference be- tweentheparticleandthefluid.Therearemanycorrelationsavail- able to calculate fluid dragon a solid spherical object. However, mostofthesemodelshavebeendevelopedforisothermalorclose to isothermal conditions. This makes these models fail at high temperaturedifferences, since variationsofpropertieshave tobe considered inorderto accuratelycalculate thedrag. The Schiller- NaumannmodelSchillerandNaumann(1935)forthedragcoeffi- cient,givenas:

C_D=24

Re

(

^{1+0.15Re0.687}

)

^{, (14)}

is a widely used drag model. Recently, Ellendt et al.

Ellendt et al. (2018) have suggested a correction factor (

φ

⁾

for the Schiller-Naumann correlation considering non-isothermal effects:

CD=24

Re

(

^{1+0.15Re0.687}

) φ

^;

φ

^=0.273

(

¹−0.883Re

)

ρ

∞

ρ

^{p −1}

+1, (15)

when1<Re<130.Here,theReynoldsnumberisevaluatedatthe surfacetemperatureofthe sphere,

ρ

∞ isthe densityofthe fluid in the far-fieldand

ρ

^{p is the densityof the fluid at the particle}

surface.Thefluid densityenteringtheexpressionforthedragco- efficient(CD=F/(⁰.5

ρ

^U2A)^{)isatfar-field}conditions.

2.4. EffectsofStefanflowondragcoefficients.

Similar to the Nusseltnumber, the models developed for the combustion and evaporation of sprays are available for the drag coefficient under the influence ofa Stefanflow (Yuen andChen, 1976; Eisenkalam et al., 1967; Renksizbulut and Yuen, 1983b).

One common approach is the so-called one-third rule proposed by Yuen and Chen Yuen and Chen (1976). The one-third rule uses ordinary drag models, for example the one of Schiller- NaumannSchillerandNaumann(1935)(seeEq. (14)),foranevap- oratingdroplet,butwiththeReynoldsnumbercalculatedas:

Re=

ρ

∞UD

μ

^{t , (16)}

where

μ

^{t isthedynamicviscosityobtainedatthereferencetem-}

perature, as given by Eq. (3), with a weight factor of t=1/3. This model is applicable in the range of 1<Re<2000 and 0<

B_T<3.Thesameresultwasconfirmedby RenksizbulutandYuen Renksizbulut and Yuen (1983b) for an evaporating droplet from theirsimulations.However,theapproachdescribedabovedoesnot includea dependencyon theStefan velocity andistherefore not expected to be suitable unless the Stefan flow velocity is small comparedtothevelocityofthemeanflow.

Studiesof theeffect ofStefanflow on thedragcoefficient for generalized cases have always assumed isothermal conditions as pertheauthorsknowledge.MostrecentworksaredonebyJayaw- ickramaetal.Jayawickramaetal.(2019),KestelKestel(2016)and Miller&BellanMillerandBellan(1999).Thelattertwohavedevel- opedempiricalmodelsforthedragcoefficientofasphericalobject

Fig. 1. Computational domain for the simulations. i, i = 1 to 5 representing the coarsest mesh to finest mesh. D _−x,iis the distance from the centre of the sphere to negative x-direction and D ₊x,iis the distance from the centre of the sphere to positive x-direction in level i (See the Table 2 ).

witha Stefanflow. Kestel’s modelisapplicable forawider range of Stefan flows (0<Re_sf20 and Re200). Both models have severalfittingparameters. Jayawickramaetal.(2019) developeda modelbasedonaphysicalinterpretationofthedragthatrequired only one fittingparameter. This modelwas validatedagainst nu- merical simulationsinthe rangeof−1Re_sf3(anegativeRe_sf means inward Stefanflow)andRe<14.All threemodels areap- plicable for isothermal conditionsonly. Therefore, it is important tostudytheeffectofaStefanflowonthedragcoefficientinclud- ingthermaleffectsaswell.

3. Methodology

In thecurrentwork,numerical simulationsare carriedout for a flowaround astatic, sphericalparticle withconstantsize using OpenFOAM. The simulationdomain andboundariesareshownin Fig. 1. The incoming gas flow to the simulation domain is uni- formanditstemperatureiskeptat1400K.AuniformStefanflow is givenasa boundaryconditionattheparticle surface.Different casesare simulatedbyvaryingspheresurfacetemperature,diam- eterandincomingflowvelocity,resultinginavarietyofReynolds numbers. Variation ofproperties with temperature is considered (SeeAppendixAformoredetails.).TheReynoldsnumberiswithin the limit of steady, axi-symmetric flow (Re<210) (Johnson and Patel, 2017) andthe Mach numberofthe flow is well below0.1.

Therefore,theflowisessentiallyin-compressible.Theintra-particle heattransferisnotconsideredandtheparticletemperatureiskept uniformbothinspaceandtime.Radiativeheattransferisalsone- glected.Thefluidisgovernedbythesteady,incompressible,lami- narflowequations,wheremassconservationyieldsthecontinuity equationas:

∇

^·

( ρ

^−→u

)

^{=0, (17)}

whilemomentumconservationgives:

( ρ

^−→u ·

∇ )

^−→u =−

∇

^p+

∇

·

μ

^[

∇

^−→u +

∇

^−→uT−2

3

( ∇

· −→u

)

^−→I]. (18) Finally,fromenergyconservationweget:

∇

^·

( ρ

^cp−→uT

)

⁼⁻

∇

^·

λ∇

^{T. (19)}

Eqs.17,18and19werediscretizedusingsecond-orderschemes withthefinitevolumemethod.

3.1. Boundaryconditions.

The temperatureofthe inletboundary is keptat 1400 K.The exitof the domain is considered asan outflow boundary, where thegradientsofthevelocityandtemperaturearesettozero.The boundariesattheside ofthedomainare treatedasslipwalls. In theslip wallboundary condition, thevelocity componentnormal tothe wall iszero.Inaddition, thegradients oftemperatureand theothervelocitycomponentsinthenormaldirectiontothewall are alsosetto bezero.Along theaxisof symmetry,a symmetric boundaryconditionisapplied.IntheSymmetricboundarycondi- tion, thevelocity component normalto the symmetry plane and the gradients ofall the other properties normalto the plane are settozero.Onlyaquarterofthedomainissimulatedastheflow isaxisymmetric.

A Cartesian mesh is used for the simulation. The immersed boundary method (IBM) was applied for the implementation of the solid boundary. In this work, the discrete forcing approach (Mittal andIaccarino, 2005), which directlyapplies the presence of a solid body through boundary conditions (Jasak et al., 2014), isused.Thevalue ofanyparameter ofacellthat crossestheim- mersedboundaryiscalculatedbyinterpolatingvaluesbetweenthe immersed boundary and neighboring cells (Fadlun et al., 2000).

The Stefan velocity is considered as a uniform velocity normal

Table 1

Conditions maintained for far-field velocity, particle diameter and particle temperature. Far-field temperature was kept at T ∞= 1400 K.

Condition Inlet velocity ( m / s ) Diameter ( mm ) T p(K ) Re f

1 0.5 1.0 400 4.88

1200 2.66 1600 2.10

2 3.0 0.5 400 14.64

1200 7.98 1600 6.31

3 3.0 1.0 400 29.29

1200 15.98 1600 13.74

to theimmersedboundary (Dirichletboundary condition).Foran outwardly directed Stefanflow, thetemperatureoftheoutflowis equaltothesurfacetemperatureoftheparticle.

Thepressuregradientissettozeroatthesolidboundary(Neu- mann boundary conditions). Treatment of Neumann and Dirich- let boundaryconditionsintheimmersedboundarymethodisex- plainedinJayawickramaetal.(2019).

3.2. Simulationconditionsandprocedure

For all the simulations in this work, the fluid (including the fluidoftheStefanflow)wasassumedtobepurenitrogen.Thein- letvelocity, diameteroftheparticleandtemperaturerangeofthe fluidandthespherewereselectedbasedonpulverizedcombustion andgasificationapplicationsatatmosphericpressure.Thevelocity at theinlet varied between0.5-3 ms⁻¹ andthe diameterofthe particle is between 0.5-1.0 mm. The rangeof Stefan flow veloci- tieswasselectedbasedonresultsfromKreitzbergetal.(2016)and Umeki et al. (2012) for devolatilization and char conversion of biomass.Thechoiceofbulkfluidtemperature(1400K)isbasedon therangeoftypicalbulkfluidtemperaturesobservedinpilotscale experiments of entrained-flow gasification (Sepman et al., 2017).

Fuel particles in entrained flow gasifiers are usually colder than the surroundinggasbecause ofpredominantlyendothermicreac- tions andthe lack ofan oxygen rich atmosphere, except forthe near burner zone. The particletemperaturecan, however,exceed thegastemperaturebyca.200Kinpulverizedcombustion,where oxygen isavailable forcharcombustion reactions(Lietal., 2018).

Therefore, we selected three different fuel particle temperatures (Tp=400,1200,and1600K),each representingdrying,chargasi- fication, andcharoxidation stages,respectively. Theparticle tem- peratures andfar-fieldconditionsstudiedin thisworkare shown inTable1.

We used the OpenFOAM environment foam-extend-4.0 (Welleretal., 1998) forthesimulations. Theimmersedboundary solverforincompressible,steady-stateconditionswasmodifiedto accountfornon-isothermal,variabledensityandvariableproperty conditions. The solver uses quadratic interpolation (Jasak et al., 2014) forthe reconstruction of the solid phase boundary condi- tionsintotheclosestfluidcells.

The preliminary domain size and mesh resolution was se- lected basedon previousstudies (Jayawickrama etal.,2019;Con- stantetal.,2017;Richter andNikrityuk,2012)forisothermalflow around a sphere. The inlet conditions and Stefan flow velocities are similar to the ones used for the isothermal simulations in our previous work (Jayawickrama etal., 2019). Therefore,the do- mainsizeisunchangedforthecurrentnon-isothermalsimulations (64D×32D×32D). There are, however, two main differences in thenon-isothermalcasescomparedtotheisothermalcases.

The first difference is that a reduction (increase) of particle temperature increases (decreases) the Reynolds number (Re), re- sulting in a thinner (thicker) boundary layer for non-isothermal

Table 2

Distance from the centre of the particle in diameters ( D ) in the computational domain (See Fig. 1 ).

(a) Mesh I

i D _−x,i D _+x,i D _y,i, D _z,i i/D

1 16 48 16 0.32

2 3 6 3 0.16

3 2 5 2 0.08

4 1.5 3 1.5 0.04

5 1.2 2 1.2 0.02

(b) Mesh II

i D _−x,i D _+x,i D _y,i, D _z,i Delta i/D

1 16 48 16 0.32

2 6.5 12 6.5 0.16

3 5.5 10 5.5 0.08

4 4.5 6 4.5 0.04

5 3.5 4 3.5 0.02

(c) Mesh III

i D _−x,i D _+x,i D _y,i, D _z,i i/D

1 16 48 16 0.16

2 3 6 3 0.08

3 2 5 2 0.04

4 1.5 3 1.5 0.02

5 1.2 2 1.2 0.01

conditions.Meshrefinementteststhereforehadtobecarriedout.

ThetestswerecarriedoutwiththehighestReynoldsnumbercon- ditions(condition3ofTable1withparticletemperature400K)and withthesmallestpossibleboundarylayerthickness(inwardStefan flowcondition).Twomeshrefinementlevelsweretested,asshown inTable2(MeshIandMeshIII).

The other difference between the isothermal and non- isothermal cases is due to the difference between the thermal (

δ

th)andtheviscousboundarylayerthickness(

δ

v^is).AsthePrandtl number(Pr)islessthan1,thethermalboundarylayerthicknessis largerthantheviscousboundarylayerthickness(

δ

th>

δ

vis).There- fore,thesizeofthemeshrefinementregionshavetobeexamined.

This was carried out for the lowest Reynolds number condition (condition1ofTable1withparticletemperature1600K) withthe largestpossible boundarylayerthickness (highestoutward Stefan flow).Tworefinementregionsizesweretested,whichisshownin Table 2(Mesh I andMesh II).Difference betweenMesh I,II and IIIwereverysmallinCD andNuandthevelocityandtemperature fields around the boundary layer were also identicalwhen com- paringallthemeshes.Therefore,MeshIIIwasusedforthesimu- lations.

Table3showstheselectionofmeshrefinementlevelsandsize ofrefinement regions usedforthe simulations inthispaper. The finalmeshforalltheconditionswasMeshIIIwiththehighestre- finement0.01D.

3.3. EstimationofdragcoefficientandNusseltnumber

Thedragcoefficientisadimensionlessquantity usedtorepre- sentforcesactingonthesurfaceofabodyimmersedinafluid.For asphericalbodywithradiusR,itcanbecalculatedas:

C_D,f=

−

→F_P,x+−→ F_visc,x

1

2

ρ

fU_∞²

( π

^R2

)

^{, (20)}

where

ρ

f is thefluid densityoffilm condition.The pressureand viscousforcesaregivenas

−

→F p=

S

pp−→nds, (21)

and

−

→F_visc=−

S

μ

^p

( ∇

^−→u +

∇

^−→uT

)

^−→nds, (22)

Table 3

Mesh refinement results and refinement domain size results as explained in section 3.2 . The drag ( C D) and Nusselt number ( Nu ) calculated at far-field conditions and Stefan Reynolds number ( Re s f) calculated at particle surface condition.

Re sf Mesh C D Error (% of mesh III or II) Nu ∞ Error (% of mesh III or II)

-7.98 mesh I 3.01 10.12 5.40 0.15

mesh III 3.36 - 5.39 -

2.36 mesh I 10.52 0.25 2.16 2.44

mesh II 10.55 2.11 -

mesh III 10.94 2.16 -

Fig. 2. Drag coefficient ( C D) at film condition for the case where there is no Stefan flow. Lines: Correlations of Ellendt et al. Ellendt et al. (2018) at different particle temperatures (400 K,1200 K,1600 K), symbols: results from our numerical simula- tions. Green: isothermal. Cyan: T p= 400 K. Red: T p= 1200 K. Blue: T p= 1600 K.

respectively.Here,theintegrationisoverthesurfaceSofthepar- ticle. In the above, pp is the extrapolated pressure at the parti- clesurface.Onlythecomponents−→

Fpand−→

F_visc inthedirectionof themeanflowareaccountedforwhencalculatingthedragcoeffi- cient, sincethe othercomponents arecanceled duetosymmetry.

The Nusselt numberis calculatedbased on theoverall difference inenthalpyfluxattheboundariesofthesimulationdomain.Here, thefar-fieldbasedNusseltnumberiscalculatedasfollows:

Nu_∞= (

ρ

^−→ucpTS)in+−→u_{s f}(

ρ

^cpTS)sph−((

ρ

^−→ucpT)^−→ndS)out

S_sph(^Tp−T_∞) ^×

2R

λ

∞, (23) wheresubscriptsin,outandsphreferstotheconditionsatthein- letboundary,theoutletboundaryandtheparticlesurface,respec- tively,andSisthesurfaceareaoftherelevantboundary.

3.4. Validation

In order tovalidate the code, simulations were carriedout to examine if the code reproduces known resultsboth for the drag coefficientandtheNusseltnumber.

For the validation of the code with respect to the drag coef- ficient, non-isothermal simulations without Stefan flow were car- riedout.Thedragcoefficientsobtainedfromthesimulationsbased on Eq. 20 were compared with the model suggested by Ellendt etal.Ellendtetal.(2018)(seeEq.15).AsshowninFig.2,thenu- merical results show good agreement withthe model of Ellendt et al.Please note that, when determining the modelpredictions, theReynoldsnumberisbasedonfilmconditions.

We are interested in the Nusselt number at strongly non- isothermal conditions, i.e., where the temperaturedifference be- tween theparticle surface andthe far-fieldis high(>100 K). In

Fig. 3. Normalized drag C D,s f/C D,0at film condition and Normalized Stefan flow ve- locity U s f/U _∞. C _D,0is the drag coefficient without Stefan flow. Simulation conditions:

U ∞= 3.0 m s ⁻¹, T ∞= 1400 K, and D = 1.0 mm.

order to validate the code with respect to the Nusselt number, simulations were carried out with a strong temperature differ- ence, but without Stefan flow. The results were compared with theRanz-Marshall model(Eq.2),which isapplicable forstrongly non-isothermalconditions(see section2.1).Table4showsagood agreementbetweenthenumericalresultsandthemodeldata.

4. ResultsandDiscussion

4.1. TheeffectofStefanflowonthedragcoefficientunder non-isothermalconditions

Bycomparingthesimulationresultsobtainedatisothermaland non-isothermal conditions, it is possible to isolate the physical effects of the Stefan flow (e.g. due to the change in boundary layerthickness)fromthermaleffects(e.g.variationofthermophys- icalpropertiesduetothechangeintemperature).Figure3shows thenormalizeddragcoefficient(C_{D,s f}/C_D,0)againstnormalizedSte- fanflowvelocity(U_{s f}/U_∞)forbothisothermalandnon-isothermal conditions(condition3ofTable1).

Thedragcoefficientisnormalizedbythecorrespondingdragas obtainedwithoutaStefanflow(C_D,0).Here,C_D,0andC_{D,s f} arecalcu- latedbasedonfilmcondition(SeeEq.20).Ascanbeseenfromthe figure, the temperature difference has a significant effect on the slopeofthecurve,especiallyforhightemperaturedifferences.The dragreductionbytheStefanflowismoresignificantwhenthepar- ticletemperatureislowerthanthesurroundinggas(T_p<T_∞)and viceversa. Thesame behavior canbe observed (notshownhere) forconditions 1and2 (see Table1) as well. Itmeans that apart fromthephysicaleffectsoftheStefanflow,thethermaleffecthas tobe considered todescribe thechangeofC_D fornon-isothermal conditions.

Table 4

Comparison of Nusselt numbers ( Nu ) and the drag coefficient ( C D) without Stefan flow from simulations and the Ranz-Marshall model ( Eq. 2 )), respectively the model of El- lendt et al. ( Eq. 15 ). Far-field temperature ( T ∞) is 1400 K for all the cases. Conditions 1-3 are listed in Table 1 , while conditions 4-5 are presented in the following: condition 4: D = 1 . 0 mm and U ∞= 5 . 94 m s ⁻¹, Condition 5: D = 1 . 0 mm and U ∞ = 11 . 88 m s ⁻¹

Condition T p Re f Nu Error C _D,f Error

K - Sim Model % Sim Model %

1 400 4.88 3.16 3.32 4.8 5.90 6.48 8.95

1200 2.66 2.84 2.89 1.7 11.39 11.53 1.22 1600 2.10 2.74 2.78 1.4 13.67 14.38 4.98

2 400 14.65 4.04 4.28 5.6 2.76 2.67 1.34

1200 7.99 3.57 3.55 0.6 5.01 4.76 2.15 1600 6.31 3.41 3.35 1.8 5.90 5.95 0.30

3 400 29.29 5.16 5.23 1.3 1.74 1.68 0.92

1200 15.98 4.34 4.20 3.3 3.09 2.92 1.55 1600 12.63 4.06 3.91 3.8 3.60 3.64 0.24

4 1600 25.0 4.95 4.69 5.5

5 1600 50.0 6.22 5.80 7.2

ThethermaleffectsoftheStefanflowcanbestudiedbyinves- tigatingFig.4,whichshowsthevariationofthevelocityandtem- peraturefieldsintheboundarylayer.WithouttheStefanflow(blue lines), thevelocity gradient ofthe non-isothermal case (Tp<T_∞) is slightly larger than that of the isothermal case. Nevertheless, we canseefromFig.2that thedragcoefficient forTp<T_∞ (non- isothermal case)islowerthanfortheisothermalcase.Thisisbe- causethecontributionfromthechangeinthermophysicalparam- etersismoresignificant thanthechangeinboundarylayerthick- ness (Eq.20). Tobemorespecific: onewouldexpectthedragco- efficienttoincreasewhentheboundarylayergetsthinner(higher velocity gradients), butthis effect is more than compensated by thedecreaseinviscosityduetothelowertemperature.Inessence, thelocalReynoldsnumberisincreasedwhentheparticletemper- aturebecomeslowerthanthefar-fieldtemperature,anditisclear fromFig.2thatthedragcoefficientdecreasewithincreasingRe_f.

In contrast,the same non-isothermal caseshowsa more pro- nounced expansion of the velocity boundary layer with an out- ward Stefan flow (red lines) than does the isothermal case. This pronounced changeinthevelocityisduetotheexpansion ofthe gasfromtheStefanflow asitisheated.Since itisthevelocityof the Stefanflowthat iskept constant betweendifferentcases,the totalmassfluxduetotheStefanflowismuchhigherforthenon- isothermal case(sincethe fluid densityismorethan three times higher at400Kthan at1400 K).Thismeans that astheinitially coldgasemittedfromtheparticleat400Kisheatedup,itacceler- atesandpushestheboundarylayeroutwards.Infact,thenormal- ized temperature plotin Fig. 4bshows the decreasein gas tem- peratureneartheparticlesurfacewithoutwardStefanflow.Asfor the inward Stefan flow, both velocity and thermal boundarylay- ersshowedexactlyoppositetrendsfromtheoutwardStefanflow, i.e.steeper velocity gradient andthinnerthermal boundarylayer.

These observations imply the importance to consider thechange inthermophysicalparameterswhenmodellingthedragcoefficient under non-isothermalconditions.Therefore, themodeldeveloped inourpreviouspaper(Jayawickramaetal.,2019),whichwasbased on isothermal simulations, needs tobe extended to considerthe effectofthevariationofthermo-physicalproperties.

Our previous study under isothermal conditions (Jayawickrama et al., 2019) showed that the drag coefficient changesduetoaStefanflow.Thischangeisprimarilycausedbya modification oftheviscous forcesduetothechangeinboundary layer thickness. Followingthe idea in Jayawickrama etal.(2019), the current studyuses a simple model forthe effect ofa Stefan flow on the drag coefficient. It is related to the change in the volume of the boundary layer due to the Stefan flow, and is

Fig. 4. (a) Normalized velocity in the mean flow direction ( U x/U _∞); (b) Normal- ized temperature ( T /T ∞). Both figures are drawn as functions of the normalized distance from the centre of the sphere ( y/R ) along the y -axis ( θ= 90 ^◦). Simulation conditions: U ∞= 3.0 m s ⁻¹, T ∞ = 1400 K, and D = 1.0 mm. Solid lines: isothermal and Dashed lines: non-isothermal results ( T p= 400 K).

proposedas:

C_{D,s f}=CD,0×CD,r, (24)

whenCD,0 isthedragcoefficientundernon-isothermalconditions without aStefan flow(see e.g. Eq. (15)), andC_D,r is acorrection termthat accountsfortheeffects ofa Stefanflow, inadditionto anythermaleffectsofthisStefanflow. Thiscorrectiontermtakes intoaccounttwoeffects:oneisduetothetemperaturedifference between the particle surface andthe far-field while the other is duetothevariationofthetemperaturefieldduetotheStefanflow.

Botheffectscanbeaccountedforbyusingamodifiedtemperature (T˜) basedon thevolumetric contributionof theStefan flow (V_{s f}) and its temperature (T_{s f}=Tp), and the volume of the boundary layerwithoutStefanflow(VB)anditstemperature(T_f= ^T∞+2^Tp);

T˜=VBT_f+V_{s f}Tp

VB+V_{s f} , (25)

where

V_{s f}=4

π

^R2Us f

τ

⁽²⁶⁾

is the added volume dueto the Stefan flow,with the flow time- scalegivenas:

τ

= 2

(

^R+

δ )

U_∞ . (27)

Furthermore,thevolumeoftheboundarylayerisgivenas:

V_B=4

3

π (

^R+

δ )

³−4

3

π

^R3, (28)

when

δ

= 2AR

Re_f, (29) istheclassicalboundarylayerthickness,where

Re_f=

ρ

fU_∞d

μ

f

(30)

andAisamodelconstant.BysubstitutingV_{s f} andV_BinEq.25with the correspondingexpressions foundinEq.26and28weobtain:

T˜=Tf+^U_U^{s f}_∞f

(

^Ref

)

^Ts f

1+^U_U^{s f}_∞f

(

^Ref

)

^{, (31)}

where

f

(

^Ref

)

=3

(

¹+ 2A

Re_f

)

¹

(

^3A

Ref

+6

(

^A Ref

)

²⁺⁴

(

^A Ref

)

³

)

^{. (32)}

Now,T˜willbeusedtocalculatethedragcoefficientwithoutStefan flow(C˜_D,0)suchthat thenon-isothermal modelforC_{D,s f} becomes:

C_{D,s f}=C˜D,0×CD,r, (33)

whereC˜_D iscalculatedfromthemodifiedSchiller-Naumannequa- tion(Eq.15)fornon-isothermalconditions:

C˜D,0=24

Re˜

(

¹+0.15Re˜^0.687

) φ

^;

φ

^=0.273

(

¹−0.883Re^˜

)( ρ

∞

ρ

^{p −1}

⁾

^+1,

(34) whereRe˜ istheReynoldsnumbercalculatedwithpropertiesatT˜. C_D,r iscalculated based onthe modeldeveloped fromisothermal simulations(Jayawickramaetal.,2019)where:

CD,r= V_B

V_{s f}+VB = 1

1+^U_U^{s f}_∞f

(

^Ref

)

^{. (35)}

Intheabove,thetildeoverC_D isusedtohighlightthatitisbased on properties calculated atT˜. The constant A is calculated using non-linearleast-squaresregressiontominimizetheerrorbetween themodelandthesimulationresults(nlinfitinMATLAB).Thefinal valueofAis2.93.

Fig.5,whichshowsthedragcoefficientasafunctionoftheSte- fanflowReynoldsnumber,comparestheabove modelwithsimu- lationresults.Themodelisanextensionofthepreviousisothermal model presented in Jayawickrama et al. (2019). This new model capturestheeffectsofnon-isothermal,uniformbulkflowanduni- formStefanflow.Modeldataandsimulationresultsarematching well and it has only one fitting parameter (A). The model has a good qualitative performance forboth negativeand positive Ste- fanflowconditions,anditisbasedonaphysicalinterpretationof thermaleffectsduetopropertyvariationsandtheStefanflow,and physicaleffectsduetopressure,viscosityandStefanflow.

The(relative)root-mean-squareerror(Eq.36)withallthedata in Fig. 5 was 9.6%. The error was relatively high for Tp=400 K (Fig. 5a), with the maximum value reaching 28%. When only considering the data from the temperature difference of 200 K (Fig. 5b-c), the maximum relative error of the model was 6% andtheroot-mean-squareerrorwas4.6%.Root-mean-squareerror (RMSE_C

d)iscalculatedasfollows:

RMSECd=100×

[

(

^Cd,model−Cd,simulations Cd,simulations

)

^2]

n , (36)

whereC_d,modelisthevaluepredictedbythemodelEq.33-(35)and C_d,simulations is thevalue calculated from the simulations andn is thenumberofsimulationsconsidered.

The models aretested andvalidatedfor theparticle Reynolds number range of 2Re_f30, Stefan Reynolds number range of−8Re_{s f,p}25,andtemperaturerangeof400KT_p1600K withuniform Stefanflow. The developed modelshould be appli- cable for the valid temperature ranges of the modified Schiller- Naumann model (Eq. 34). However, one should be careful when extrapolatingtheapplicabilitybeyondtherangeofvalidationcon- ditions.Forexample,themodelmightnot bevalidathigherpar- ticle Reynolds number due to flow separation or the change in therelativemagnitudebetweenthepressureforceandtheviscous force.

4.2. NusseltnumberwithStefanflow

Murphy&Shaddix MurphyandShaddix (2003)hasdeveloped a theoretical model that accountsfor the effect ofa Stefan flow when calculatingthe Nusseltnumber ofa sphere immersedin a quiescentfluid(SeeEq.11).Intheir model,theNusseltnumberis calculatedasNu_M=Nu₀fcorr,whereNu₀=2istheNusseltnum- berofasphericalparticlewithnoStefanflowinaquiescentfluid and fcorr isa correction termthat accountsfor the effectof the Stefan flow. One way to apply this model directly for the cases withconvectiveflows istoreplacetheNusseltnumber,Nu₀,with theonewithaconvectiveflow,asgivenbye.g.theRanz-Marshall model. However, the prediction withthis approach doesnot de- scribethesimulationresults.Thesameobservationwasdiscussed by Kestel Kestel (2016), who proceeded to develop an empirical modelwithseveralfittingparametersEqs.12-(13).

As discussed in the previous section, the temperature in the boundarylayerchangesduetotheStefanflow,especiallywhenthe temperaturedifferencesaresignificant. Thischangeshouldbe re- flectedinthecharacteristictemperaturewhencalculatingtheNus- selt number. In this work, we apply a multiplication law to de- scribe the effectof a Stefan flow (by Eq.11) and the effectof a convective flow (Eq.2), but considering the change in character- istictemperature.Thisapproach inpracticecalculatestheNusselt

Fig. 5. Comparison of the drag coefficient from the model, i.e. Eqs. 31 - 35 , (lines) and the simulations based on T (symbols). ˜ Particle temperature ( T p) is (a) 400 K (b) 1200 K or (c) 1600 K. Condition 1: U _∞= 0.5 m s ⁻¹and D = 1.0 mm. Condition 2: U _∞= 3.0 m s ⁻¹and D = 0.5 mm. Condition 3: U _∞= 3.0 m s ⁻¹and D = 1.0 mm.

numberbasedonthethermophysicalpropertiesusingthevolume averaged temperature derived earlier(Eq.31).The model forthe Nusseltnumberneedstobeapplicableforconvectiveflowsaround aspherewithhightemperaturedifferences. Here,wehaveuseda Ranz-Marshall typemodelby parameterfittingtheoriginal Ranz- Marshall model withsimulation datawithout Stefan flow, to ob- tain:

Nu˜ =2+0.570Re˜^0.537Pr˜^1/3, (37) where Re and Prwere calculatedbased on the volume averaged temperature, T˜,asgiveninEq.31.Now wecan replacethe Nus- seltnumberwithoutStefanflow(Nu)inMurphy&Shaddixmodel (Eq.11) with the model presented in Eq.37,such that the final model for Nu, accounting for non-isothermal effects and Stefan flowreadsas:

Nus f,f=Nu˜ q

e^q−1, (38)

whereq= Pr_fRe_{s f,p}

Nu˜ andNu_{s f,f} calculated basedonfilmcondition forthethermalconductivity(

λ

f).TheStefanflow Reynoldsnum- ber (Re_{s f}) is calculatedbased on particlesurfaceconditionwhile thePrandtlnumber(Pr)iscalculatedbasedonfilmcondition.

It is clearthat the volumeaveraged temperaturemust lie be- tween theparticle temperature(Tp) andthe far-fieldtemperature

(T_∞).From thedefinitionofthe volumeaveraged temperature,as givenbyEq.25(respectivelyEq.31),itcanbe shownthatthisis notthe casewhenV_{s f}/V_B<−0.5.(Thiscorresponds toasituation wherethereisaverystronginward Stefanflow.)Thismeans that theexpressiongivenbyEq.25cannot beusedtodefinethevol- umeaveragedtemperatureforsuchacondition.Therefore,thevol- umeaveraged temperatureisassumedtobeequaltothefar-field temperaturewhenV_{s f}/V_B<−0.5.Thismeansthat,

T˜=

Tf+^UU^{s f}∞f(Ref)Ts f

1+^UU^{s f}∞f(Ref)

(

^Eq.31

)

^forVs f/VB−0.5,

(

³⁹

)

T_∞ forV_{s f}/VB<−0.5,

(

⁴⁰

)

where f(^Ref)^{iscalculatedfromEq.32andA}=0.4.

Tovalidatethemodel,theNusseltnumberwascalculatedfrom thesimulationwiththeconditions1,2,and3(seeTable1)includ- ingonenegativeStefanflowcasewithV_{s f}/V_B<−0.5.Fig.6depicts thecomparisonbetweensimulationresults(symbols)andthepre- dictionsobtainedwiththemodelpresentedinEq. (38)(lines).

The(relative)root-mean-squareerror(Eq.36afterreplacingthe termC_d withNu) with all the data in Fig.6 was 12.6%.The er- rorwasrelativelyhighforTp=400K(Fig.6a),withthemaximum valuereaching73%.Whenonlyconsideringthedatafromthetem-

Fig. 6. The Nusselt number comparison between the model ( Eq. 38 -lines) and simulation (symbols) data with Stefan flow. Particle temperature ( T p) is (a) 400 K (b) 1200 K or (c) 1600 K. Condition 1: U ∞= 0.5 m s ⁻¹and D = 1.0 mm. Condition 2: U ∞= 3.0 m s ⁻¹and D = 0.5 mm. Condition 3: U ∞= 3.0 m s ⁻¹and D = 1.0 mm.

peraturedifferenceof200K(Fig.6b-c),themaximumrelativeer- rorofthemodelwas9%andtheroot-mean-squareerrorwas3.8%. ThemodelisdevelopedforcalculatingtheNusseltnumberfor a sphericalparticlewithuniformStefanflow, immersedinauni- form convective flow. It was validated for the Reynolds number (Re_f) 2Re_f30,StefanReynoldsnumber(Re_{s f,p}) −8Re_{s f,p} 25 and temperature range 400 K-1600 K for nitrogen gas atmo- sphere.TheparametersfortheNusseltnumberwithoutStefanflow (Eq.37)wereestimatedbyfittingthesimulationdatapresentedin thisstudy.Differentsetsofparametersmightbeapplicablefordif- ferentReynoldsnumberandtemperatureranges.

ThemodelforthedragcoefficientandtheNusseltnumberwere both developed by assuming that the change in temperature in- side the boundary-layer occurs due to variations in Stefan flow velocity, Stefan flow temperatureandfar-fieldtemperaturealone.

Thiswouldnotbethecasewhenthereareotherphenomenathat affectthe boundary-layertemperature,such ase.g.,homogeneous reactions. The model is based on the assumption that the pres-

sureforce andtheviscous force areofthe sameorderofmagni- tudeandthatonlytheviscous forceisaffectedbytheStefanflow (seeJayawickramaetal.(2019)).Thismightnotbetrueforhigher Reynoldsnumbers.

5. Conclusions

Theeffectofa Stefanflowon thedragcoefficientandNusselt numberwasstudiedforauniformflowaroundasphericalparticle.

Theeffectwasinvestigatedatnon-isothermalconditionsusingre- solvednumericalsimulations.Particlediameter,slipvelocity,parti- cletemperature,andStefanflowvelocityfrom/totheparticlehave beenvaried duringthesimulations. Therangeof StefanReynolds numberof−8Re_{s f,p}25,Reynoldsnumberof2Re_f30and particletemperatures(Tp)of400K,1200 Kand1600Kwerecon- sideredinthesimulations.Thefar-fieldtemperature(T_∞)waskept constantat1400K.