Computationally efficient coarse-graining XDEM/CFD modeling of fixed-bed combustion of biomass

(1)

ContentslistsavailableatScienceDirect

Combustion and Flame

journalhomepage:www.elsevier.com/locate/combustflame

Computationally eﬃcient coarse-graining XDEM/CFD modeling of ﬁxe d-b e d combustion of biomass

Jingyuan Zhang

^a

, Tian Li

^a^,^b^,^∗

, Henrik Ström

^c^,^a

, Terese Løvås

^a

aDepartment of Energy and Process Engineering, Faculty of Engineering, NTNU — Norwegian University of Science and Technology, Trondheim, Norway

bRISE Fire Research, Tiller 7092, Norway

cDivision of Fluid Dynamics, Department of Mechanics and Maritime Sciences, Chalmers University of Technology, SE-412 96 Gothenburg, Sweden

a rt i c l e i nf o

Article history:

Received 16 November 2020 Revised 10 November 2021 Accepted 11 November 2021

Keywords:

CFD Biomass

Fixed-bed combustion Coupling

a b s t r a c t

Inthe multi-scalemodeling ofadense particlesystem, the particlephase and thegas phase canbe modeledonvastlydifferentscales.Thecouplingbetweenthetwomodelshasacriticalinfluenceonthe predictionsobtainedfromthecombinedframeworkbutcanbeaccomplished inavarietyofwaysun- derdifferentassumptions.Inthiswork,atransient3Dmodelusinganewcouplingapproachforfixed- bedcombustionofbiomassispresented.ThedevelopedmodelisformulatedasanEulerian-Lagrangian framework.Aparticlegrid,generatedbasedonthefluidgrid,isappliedasatransfergrid,andadiffu- sionoperationisimplementedtosmooththeinteractionsbetweenthegasphaseandtheparticles.The interactionsbetweengasand solidphasesaswell asthe radiativeheattransferbetweenparticlesare considered.Theparticlemotionisresolvedbythesoft-spheremodel,whereastheconversioniscalcu- latedbasedonathermallythickparticlemodel.Allsub-modelsareoptimizedtoenhancecomputational efficiency.The3Dmodelisvalidatedbycomparingthesimulationswithlaboratory-scaleexperimentsfor afixed-bedoperatedincounter-currentcombustionmode.Thekeysimulationparametersareconfigured bysensitivity analysis.The simulationresults areingoodagreement withthe experimentalmeasure- ments,and thecombustion regimeswithdifferentairinletconditionsarewell captured.Thecoupling effects arediscussedindetail.The particlegridsizeinfluencesthe predictionofthe transientresults, andtheinterplaybetweentheheattransfermechanismsinsidethefixed-bedandthecouplingscheme isthoroughlyanalyzed.Bothinter-particleradiationandgas-to-particleconvectionplayessentialrolesin theheattransferinsidethefuelbed,whiletheinter-particleheatconductioncanbeneglected.

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

1. Introduction

Asacarbon-basedrenewableenergyresource,biomassisasub- stituteforfossilfuels suchascoal.Powergeneratedfrombiomass is mainly from thermal conversion processes. Most plants use grate-fired systems for their simplicity and versatility. The grate furnacehasalowconstructioncostandcanbeoperatedwithvar- ious feed materials with a wide range of moisture contents [1]. However, the emissions are not only due to the feed but highly influenced bytheoperatingparameters [2].Inthisrespect,abet- ter designofthegratefurnacecanimprovecombustionefficiency andreducetheemissions.Thedetailedinformationfrominsidethe furnacethat isrequiredforoptimizationcanbe hardtomeasure.

Complementarytoexperimentalmethods,computationalﬂuiddy-

∗Corresponding author.

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

namics(CFD)simulationisanalternativetooltooptimizethecom- bustionprocessofagrate-ﬁredfurnaceowingtoitscosteffective- nessandfastturnovertime.

To capture the thermochemical conversion of biomass in the fuel bed and the subsequent gas phase conversion in the free- boardregion,atransientmulti-phase simulationisneeded. There are two main approaches in conventional CFD regarding multi- phasemodels.TheﬁrstoneistheEulerian-Eulerianapproach.One exampleusingthisprincipleistheso-calledporousmedia model (PMM).TheNavierStokesequations(includingtheconservationof mass, species and energy) for the gas phase are solved using a phase fraction approach, and related source terms that describe the momentum, mass and energy transfer between the gas and solid phases. The particle sub-modelsare registered to the com- putationalgridstocalculatethesourceterms.Onecellcancontain severalidenticalparticlesorevenacertainfractionofoneparticle [3].The inter-particle momentum andmass transfers are usually ignored.Theinter-particleheattransfercanberesolvedbythepar-

https://doi.org/10.1016/j.combustﬂame.2021.111876

(2)

ticle sub-modelsonthe surfaceheatbalance, orby resolvingthe energycontinuityequation forthesolidphase,whichregardsthe particlesasacontinuum[4].PMMcanbecomputationallyeﬃcient todealwithalargenumberofparticles,butitrequiresadditional sub-modelstogetsolidpropertiesintermsofacontinuousphase fromindividualparticles.Thismodelapproachalsocannotdescribe themovementofparticlesfromﬁrstprinciples,astheinformation aboutindividualparticleshasbeenaveragedout.Insomestudies, extramodelsareemployedtopredictbeddeformationduetopar- ticleshrinkage[5].

The other approach is the Eulerian-Lagrangian approach. The gasphaseisstillresolvedbytheNavierStokesequations,whilethe particles are describedby the extended discrete elementmethod (XDEM). In theclassical discrete elementmethod (DEM) method, thepositionandmotionintimeandspaceofeachparticleorpar- celaretracked.XDEM extendstheDEMbyaddingpropertiesand various processesattachedtotheparticle.Forexample,thesingle particlethermalconversionmodelwiththeconsiderationofintra- particle heat transfer canbe coupled to theDEM usingthe soft- sphere modeltocalculatetheparticlecollisions.Plentyofstudies of such XDEM/CFD simulations (or DEM/CFD in some literature) have been conductedin thepast decade. PetersandBruch et al.

[6,7] implementeda thermallythickparticlemodelinto theDEM framework.Theinter-particleheattransferthroughconductionand radiation ismodeled, andthe simulations ofbiomasscombustion were conductedunderﬁxed-bedconditions.Wieseetal.[8]modeled particle collisions withcylindrically shaped particles in the DEM/CFD framework for pellets. In their simulations, a transfer gridwasemployedtothedatatransferbetweengasandparticles.

Mahmoudi etal. [9] modeled the beech woodchips combustion in a ﬁxed-bed in the DEM/CFD method. A good agreement between measured andsimulatedbed massloss andbedtempera- turewasobtained.Furthermore,asimulationofaplant-scalemu- nicipal solid waste (MSW) incineration on grates has been pre- sentedbyWissingetal[10].

The XDEM/CFD approach can however capture multi-physics phenomena and provide more detailed information. It is usually quitecomputationally costlyastheparticlesearchinthecollision model andthe thermally thick particleconversion model require massivecalculationtimewhentheparticlenumberbecomeslarge.

Mehrabianetal.[11]implemented anumericallyeﬃcientparticle conversion model and a particle position rearrangement strategy to reducecomputational costina gratefurnacesimulation.How- ever,themotionoftheparticlesderivedby theactual mechanical motionalongwiththegrateandshrinkage-inducedcollapseisim- portanttopredictthefuelbedlocalignitionandﬂamestability.

Apartfromtheconcernaboutthecomputationaleﬃciency,an- otherissueabouttheXDEM/CFDapproachisthecouplingbetween the two models. When the particle size is larger than the ﬂow scale that needs to be resolved for the gas phase, the particles will be larger than one computational cell andthe data transfer between the CFDmodel andXDEM is nolonger straightforward.

In conventional DEM/CFD modeling, the particle data is directly transferred betweenthecellinwhichtheparticle’scentroidislo- cated, because,forthe gasphase the particleis still viewedasa Lagrangian point particle (LPP). Therefore, the inter-phase transfers are modeled by 0-dimensionalmodels, and forparticle conversion sub-modelstheassumptionisthatthecellshouldbeable to provide the gas ﬁelds that are inﬁnitely far fromthe particle surface.Therefore,forthecouplingbetweenthegasandlargepar- ticles, special schemes forthe datatransfer between theparticle sub-modelsandthesub-particle scalecomputationalgridsarere- quired.

In fact, the coupling of ǣlarge particlesǥ with a continuous phase is a common issue in various research ﬁelds. Sun et al.

[12] made a review over the so-called ǣcoarse-grainingǥ meth-

odsandsummarizedfourmainmethods:particlecentroidmethod (PCM),thedividedparticlevolumemethod(DPVM),thestatistical kernel method(SKM) and the two-grid method (TGM).All these methodsareaimingtomaptheparticlepropertiesgeneratedfrom aLagrangianpointtoanEulerianﬁeldtosmooththecouplingbe- tween thesolid phase andthe gasphase. Thissmoothingis typ- ically intended to increase the robustness of the solver and de- creasethedependenceoftheobtainedsolutionon thenumerical gridsandsettingsused.However,mostofthesestudiesfocusedex- clusivelyonhydrodynamics,whereas one mayexpectthe various couplingissues tobe all themorepronounced innon-isothermal reactiveﬂows.

Itis alsovery importantthat thecouplingmethod canrecon- struct the gasphase information required by the assumptions of the particle sub-model, when it is used to describe the particle conversionprocessunderreactingconditions.Itisessentialtoen- surethe applicabilityof theseextended sub-modelsinthe simulation. Inour previous study[13], the effects on the simulations ofthesingleparticlethermalconversionprocessbyusingdifferent coarse-graining methods were evaluated. The studyshowed that thecouplingschemeisanintrinsicpropertyofthecombinedmod- elingframeworkandthatitisextremelydifficult(ifnotpractically impossible) to choose coupling parameters that donot influence theresultsinanyway.Thechallengesincouplingwillbemorese- vereinfixed-beds,duetotheclosepackingattainedinsuchbeds.

However, this issue isnot studied carefully inmost of the work onthesimulationoftheconversionofbiomassinﬁxed-bedsusing theXDEM/CFDapproach.

Theobjectiveofthisworkistodevelop acomputationallyeﬃ- cientmodelbased ontheXDEM/CFDapproachforthesimulation of the combustion of biomass particles. The model will be val- idatedagainst experimental measurements. Specialattention will bepaidtothecouplingbetweentheparticlesandthegasphasein the combinedmodeling framework. The thermalconversion pro- cessinalaboratory-scaleﬁxed-bedwillbeanalyzed byusingthe simulationresults.

2. Mathematicalmodeling

The model used in this work is in an Eulerian-Lagrangian frameworkdevelopedusingOpenFOAM®.Thegasphaseissolved usingtheNavierStokes equationsandtheheat andspeciestrans- port equations. The presence ofthe particle is considered by in- troducingthephasevolumefractionandthesourcetermsintothe equationswithspecificcouplingmethods,whichwillbedescribed indetail inSection 2.2.The governingequations andthe calcula- tionofthermophysicalpropertiesofthegasmixturearepresented inourpreviouswork[13].Inthiswork,theinletsuperficialgasve- locityofthefixed-bedisbelow0.42m/s,soonlythelaminarflow regime is considered. A turbulencemodel can howeverbe easily employed inthe developedmodeling framework.Forthe particle sub-models,thecorrelationsthatfitinthecorrespondingturbulent flowregimesarepreferred.Usually,such correlationsonlyrequire generalflow informationratherthandetailedturbulenceinforma- tionasinput.Theoptimizationofthecomputationalefficiencyde- scribedinSection2.4islimitedtofixedbeds.Iftheparticleshave fiercemotions,forexamplein afluidized bed, thecurrentmodel mayfailinpredictingthemovementoftheparticles.

2.1. XDEMmodel

2.1.1. Particlemotionmodel

The particle momentum isgoverned by Newton’s second law.

Thepotentialashlossduetobreakageduringcollapsesorbedmo- tion is not considered in the current model, as well as the ash transformations, which maychange the structure of the ashand

(3)

create larger agglomerates. The conventional drag force models, such as the Ergun [14] and Wen&Yu [15] drag correlations, will verylikelyover-predictthedragforceactingoncharparticleswith an integral thick andporous ashlayer, resulting inan underesti- matedminimumﬂuidizationvelocity.Togetridoftheinterference from the drag model,and considering that in the ﬁxed-bed,the motionoftheparticlesisagravity-drivenprocess,onlythegravity forceandthecontactforcesareconsidered:

m_idU_i^p dt =

n

j=1

F^c_{i j}+F^g_i, (1)

Ii

d

ω

i^p

dt = n

j=1

Mi j, (2)

wherem_iandI_iarethemassandmomentofinertia ofparticlei, respectively, U_i^p is the velocity of theparticlei,

ω

_i^pis the angular velocityofparticlei,F^c_{i j} andM_{i j} arethecontactforce andcontact torqueactingonparticleibyits jthcontact(eitherwithaparticle ora wall),respectively,nisthe numberoftotalcontactsforpar- ticle i, F^g_i is the gravity force actingon paticle i. The soft-sphere collision modelisimplementedtocalculatethecontactforce and torque.FurtherdetailscanbefoundinFernandesetal.[16].Inthis study,themotionofcylindricalparticlesarecalculatedasspherical particleswithavolume-equivalentdiameter.

2.1.2. Particleconversionmodel

Alayer-basedthermallythickparticlemodelisemployed[3,17]. The intra-particle heat transfer in the radial directionis considered. Itisassumedthat theconversionsub-processes:drying,devolatilization andcharburnout,occurininfinitelythinfronts,and these fronts divide the particle into fourlayers: wet woodlayer, dry woodlayer, char layerand ashlayer. The conversion ofeach sub-processiscalculatedbyseparatesub-models.Theparticlesub- modelsareprovidedintheAppendix,andmoredetailedmodelde- scriptionsarepresentedinourpreviouswork[13].Still,compared withthemodelin[13],thereare some adjustmentsmadeinthis work.Firstly,athresholdof1%oftheparticle’soriginalmassisde- fined onevery layer todetermine whetherto calculateheat con- ductioninsidethelayer.Ifthelayermassislessthanthethreshold, thentheheatconductionequationforthislayerwillnotbesolved andthelayer temperaturewill beassignedasthetemperatureat the outer front. Bydoing this, solving fortheheat conduction in a very thinlayer isavoided, anda muchlarger time step incal- culating theparticle heattransfercould beapplied withoutcaus- ingnumericaldiffusionandtemperatureoscillations.Italsogreatly enhancestheoverall computational efficiencyof theparticlesub- models. It should be noted that a proper thermal-mass thresh- oldmayfurtheroptimizethecalculation,whichrequiresfuturein- vestigations.Secondly,sincethecharburnoutsub-modeldoesnot consider the competition for oxygen with the combustible gases released fromthe devolatilization process that resultsin a radial flow (Stefan flow)in thiswork,therefore charburnoutis treated as a sequential process that may be initiatedonly after the devolatilization hasbeen accomplished. Anotheradjustment is that inthepreviouswork,theporosityoftheashlayerisafixedvalue, whileinthiswork,theporosityiscalculatedbyassumingtheash intrinsicdensityisaconstant.Thisisduetotheinsightsfromthe parameterstudyabouttheshrinkagefactor

η

^for^the^char^burnout

process, which is the volume ratio between the converted char layerandthegeneratedashlayer[13],thatisconductedinalater section. Thecalculatedporosityismorerealistic whenaccounting forvaryingshrinkagefactors.Alltheotheradjustmentsofthesub- modelswillbepresentedinthefollowingsections.

2.1.3. Surfaceheatandmasstransfermodel

Inthiswork,becausetheadiabatic wallboundaryconditionis used,theheatconductionandradiationbetweenthewallandpar- ticlearecurrentlyneglected.Theheatbalanceequationatparticle surfaceisthusformulatedas:

qin,i=qconv,i+ n

j=1

q_cond_,j_i+q_rad,i, (3)

where q_in_,_i is the heat ﬂow rate into particle i through the surface,andthispartof theheat causestheparticle temperatureto increase,andisconsumedby endothermicsub-processes, suchas drying.q_con_v_,_iistheheat exchangeratebetweenparticlei andits localsurroundinggasbyconvection,whichcanbecalculatedas:

q_con_v_,i=hcA_p,i

(

^T^g⁻^Tp,i

)

^, ⁽⁴⁾

whereA_p_,_iistheparticlesurfacearea,TgandT_p_,_iarethelocalgas temperatureandparticlesurfacetemperature, respectively,andhc

istheheattransfercoeﬃcient. Forcylindricalparticlesina ﬁxed- bed,theNusseltnumber(Nu)correlationproposedbySinghaletal.

[18]is used in this work.

Nup=hcdp

kg =1.77+0.29

ε

⁻⁰^.⁸¹^Re⁰p^.⁷³Pr⁰_f^.⁵⁰, (5) wheredpistheparticlediameter,kgistheheatconductivityofthe gasphase,

ε

^is^the^local^bed^porosity.

q_cond_,_j_iistheheatexchange ratebetweenparticles jandidue toconduction.Inafixed-bed,conductionismainlyduetoparticle- particle static contact, and the equation proposed by Batchelor etal.[19] andmodifiedby Zhou etal.[20] is adopted.The sim- plifiedconductionq_cond_,_j_icanbecalculatedas:

q_cond_,_j_i= 4rc,i j

Tp,j−Tp,i

1/k_p,i+1/k_p,j

, (6)

wherek_p_,_i andk_p_,_j aretheheat conductivityforparticlesi and j, respectively,r_c_,_{i j} isthe particle-particle contactradius,which can be obtainedfromthe DEM simulationbased onthe Hertzelastic contacttheory.HeretheHeron’sformulaisemployedtocalculate r_c_,_{i j}:

r_{c,i j}=2A

d_{i j}, (7)

A=

s

(

^s−rj

)(

^s−ri

)(

^s−di j

)

, (8)

s=

(

^rj+ri+di j

)

2 , (9)

wherer_jandr_iaretheradiusofparticles jandi,respectively,and d_{i j}is the distance between the centroids of particles jandi.

ThelastterminEq.(3),q_rad_,_iistheheatexchangeratebetween particleianditssurroundingenvironmentbyradiation.Acommon practiceistocreatea radiationcontrolvolume,based onparticle i’s position, andsearch forthe particle’s immediateneighbors to calculateanequivalentlocalradiationtemperatureT_eq_,_i[8,11].Then q_rad_,_icanbecalculatedas:

qrad,i=

σ

^f

(

^Teq,i4−Ti4

)

, (10) where

^,

σ

^and ^f ^are ^the emissivity, the StefanBoltzmann con- stant and the view factor, respectively. To avoid a large number ofcomputationallyexpensiveparticlesearches,here, theimmedi- ateneighborsareapproximatedbythecontactparticlesforparti- clei,which arealreadycalculated intheDEM step.Assuming all thecontactscanbetreatedequaltoparticlei,theequivalentlocal radiationtemperatureforparticleicanbecalculatedas:

T_eq_,_i⁴=

α

^Tg⁴+

(

¹−

α )

¹_nⁿ

j=1

T_j⁴, (11)

(4)

Fig. 1. Two-grid method.

where

α

îs ^the ^gas^phase ^fraction,ând ^the^quartic ^root ^meanîs

used.

The masstransfer ofreactantstothe particlesurfaceiscalcu- latedthrough the masstransfer coeﬃcienthm,whichis obtained bytheSherwoodnumber(Sh).TheanalogoftheNucorrelation is employed:

Shp=hmdp

De,g =1.77+0.29

ε

⁻⁰^.⁸¹^Re⁰p^.⁷³Sc⁰_f^.⁵⁰, (12) whereDe,gistheeffectivediffusivityofthegasphase.

2.2. Couplingbetweengasandparticle

The couplingbetweenthe XDEM andthe CFDmodelincludes twoparts.Thefirstistosamplethegasphasepropertiesrequired by the particle sub-models. Because the particles are treated as point particles in thefluid model, thedetails of theflow around theparticles arenot fullyresolved. Asaresult, point-particlecor- relations arerequiredtocalculatethemass,momentumandheat exchangebetweentheparticleandthesurroundingfluid.Thethe- oreticalrequirementsonthesamplingaregovernedbythederiva- tion of the sub-model itself. For example, when calculating the heat transfer using the Nu number correlation, one wouldwant thesamplingtoprovidethefar-fieldgasphasetemperature.Inthis step,theTGM[21]conceptisadopted.TGMisanefficientcoarse- grainingmethodtodealwithalargeparticlenumber,asthecom- putational cost is independent from the particle number [13]. A particlegridisgeneratedfromthefluidgrid.Theparticlegridcell size islarger than thatof thefluid grid,andalsolarger than the particle diameter. There are no governing equations resolved on the particle grid. It is only used to transfer andaverage the gas information to theparticles. Since the particlegrid andthe fluid grid are bothfixed, oneefficient waytomap thegas datato the particle grid isto createa grid cellID matching table beforethe simulation,asshowninFig.1.Whenthefluidgridcellispartially overlappedwithdifferentparticlegridcells,itwillbe assignedto the closest particle grid cell according to its centroid. In such a situation,theparticlegridcellsrepresentedbytheIDmatchingta- bleareapproximationsoftheparticlegridcellsgeneratedinitially.

Thismethodisevenmoreconvenientwhenunstructuredmeshis usedfortheﬂuidgrid.Aftereveryﬂuidtimestep,theparticlegrid willcollectthemeanvaluesaccordingtothetable.

Thevolumetricmeanisusedinthedatamapping,andthegas phasepropertiescanbecalculatedas:

φ

pg,i=

∀^j∈pg,i

Vf g,j

φ

f g,j

∀^j∈pg,i

V_{f g}_,_j , (13)

where

φ

pg,i and

φ

f g,j represent anygiven gas phase property in ithparticlecelland jthﬂuidcell,respectively,V_{f g}_,_j isthecellvolume of jth ﬂuid cell. The gas phase properties includes temperature,speciesconcentration,densityandthermal-physical properties,suchasheatcapacityandheatconductivity.

The second partofthe couplingis totransfer particledata to theCFDmodel.Thedatathatneedstobetransferredincludesthe mass and heat source terms (the momentum exchange between the gas and particle is neglected in this work), and the particle properties.Theseare,forexample,theparticlevolume neededto obtain the phase fraction, the equivalent particle absorption coefficient andthe equivalentparticle scatteringfactor that are re- quiredbythegasphaseradiationmodel,whichwillbepresented inSection2.3.1.Becausetheparticlecanbe muchlargerthanthe fluid grid cell, the particle data need to be redistributed on the fluidgrid.Thisredistribution shouldconsidertheparticle’s physical presenceinspaceandtherobustnessofthe solver,aswellas introduceaslittleadditionalspatialdiffusionofthedataaspossi- ble.Consequently,theredistribution oftheparticle dataneeds to beconfinedtoarelativelynarrowregionaroundtheparticlecen- ter.Inthis work,the diffusion-basedmethod(DBM) proposed by Capecelatroetal.[22] andfurtherdevelopedbySunetal.[12]is implemented.Theparticledatawillbetransferredtothefluidgrid cellwhereitscentroidislocated,then,beforethegoverningequa- tions ofthe gasphase are solved, a diffusionoperation isimple- mentedbysolving:

∂

^S

∂τ

⁼

∇

²S, (14)

where

τ

îs^the^diffusion^time^variable,ând^S^representsâny^parti-

clescalardata.Theno-flux boundaryconditionshouldbe applied tothecomputationaldomaintoguaranteetheconservationofthe passivescalar Sduringsuch a diffusionoperation. When thedif- fusion has been accomplished, the particle data will result in a distributionthat isequivalenttothe distributionwithaGaussian kernel[12,23]. Forexample,the scalar Scan be the value of the particle’s volume. After the diffusion operation, the redistributed S field divided by the volume field of the fluid grid cell repre- sentsthe smoothedvolume fraction field ofthe solid phase. The length scaleb forthis diffusionoperation, whichis equivalentto thebandwidthinGaussiandistribution,isdeterminedbythetotal diffusion

τ

total byb=

4

τ

total.Fordifferentparticledata,adiffer- entlengthscalecanbeassignedbyusingdifferentvaluesof

τ

total. The DBMsmears particledata ateach particle’sposition without calculationofmassivegridsearching,whilethetotalparticlenum- berhaslittleimpactonthecomputationalexpenses[13].Theover- allcouplingmethodisillustratedinFig.2.

There are two main beneﬁts of combining the TGM and the DBM. The ﬁrst is that by choosing a proper

τ

total, the DBM can avoid over-smearing the particle data instead of using the TGM also for the particle data redistribution. The second is that al- though acorresponding DBMsamplingmethodforthegas phase propertiesistheoreticallypossible(e.g.todiffusethegasdatathat aretemporarilystoredonaseparatemeshateachtime),manyfac- torsmayinfluencethesamplingresults,suchasmeshresolutions, meshtypes,anddiffusionoperationparameters. Itsnumericalac- curacyandefficiencyhave notbeenwell studied. Previousworks [23,24]usingtheDBMonlyappliedthesamesamplingmethodsas thePCM.The TGMcansampletheaveraged gasphaseproperties witha muchhighercomputational efficiency[13],which areless disturbedbytheGaussiandistributionofthesourceterms.

Itshould benotedthat thecurrentapproachliesbetweenthe traditionalEulerian-Lagrangianmethods(e.g.particlesaresmallin relationto thefluid cell sizes)andthe particle-resolved methods (e.g.boundary-fitted grid orimmersed-boundarymethods andso on).Thisapproachavoidsthecomputationallyheavyfullyresolved calculationbutisnotabletofulfillthetheoreticalrequirementsof theEulerian-Lagrangianapproach.Acertain degreeofapproxima- tion hasto be introduced, such asthe particle grid size andthe lengthscalebintheTGMandtheDBM,respectively.Therelation betweenthecouplingparametersandtheireffectsonthesimula-

(5)

Fig. 2. Coupling method. The solid black lines are the particle grid, whereas the grey lines indicate the ﬂuid grid.

tionneedstobeanalyzed.Byidentifyingtheintermediatecoupling length as an independent parameter of the coupling model, the couplingeffectsandthegridindependencefortheCFDsolvercan be studied separately in different respects. As mentioned above, the relevantconsiderations whenderiving theTGMandtheDBM formulations areindependent andpossiblyeven partlycontradic- tory. Thus, there is no reason to expect a priori that the same lengthscaleshouldbeusedinthesetwoaspectsofthecoupling.

2.3. Gasphasesub-models

2.3.1. Gasphaseradiationmodel

Thediscreteordinates(DO)model[25]isused tosolvethera- diative transferequation (RTE)forthe gasphase.Consideringthe particleeffects,assumingtherefractiveindexofthegasphaseis1, theRTEcanbewrittenas:

∇

·

(

^I^s

)

⁺

(

^a^g⁺^a^p⁺

σ

p

)

^I

(

^r^,^s

)

⁼^a^g

σ

^Tg4

π

⁺^E^p

+

σ

^p

4

π

4π 0

I

r,s

s·s

d

^,⁽¹⁵⁾

where Iis theradiation intensity,agis theabsorption coefficient ofthegasphase,sandrarethedirectionandthepositionvectors, respectively,îsthe scatteringphase function,îs^{a solid}ângle, ap,Epand

σ

^p are the equivalent absorption coefficient,equivalent emissionandequivalentparticlescatteringfactorduetothepres- enceofparticles.Thoseequivalentcoefficientswillhaveinfluences onthetransportoftheradiationintensity,whichcanberegarded as the resultofthe blocking effectof particles on the gasphase radiative heat transfer. The finite-volume method (FVM) is used in OpenFOAM®to discretizeEq. (15), andthe particles equivalent propertiescanbecalculatedby:

ap= N

i=1

pi

Api

V , (16)

A_pi=

π

^d²pi

4 , (17)

Ep= N

i=1

piA_pi

σ

^Tpi⁴

π

^V ^, ⁽¹⁸⁾

σ

^p⁼^N

i=1

1−f_pi

1−

pi

A_pi

V , (19)

where

pi,A_pi,T_piandf_piaretheemissivity,projectedarea,surface temperatureandscatteringfactorofparticlei,V isthecellvolume, andNisthenumberoftheparticlesinthecell.

SinceH₂OandCO₂ arethemainproductsinthegasphase,the absorptioncoeﬃcienthasbeenassumedtobe[26]:

ag=0.1

(

^wH2O+w_CO₂

)

^, ⁽²⁰⁾

wherew_H

2O andw_CO

2 are the mass fractionof H₂O andCO₂,respectively.

2.3.2. Gasphasereactions

Inourpreviouswork[13],globalreactionswereused.However, inthisstudy,thecombustionislaminarandanaccurateprediction of the ignition front is important.For this reason, a detailed re- actionmechanismisusedforthehomogeneousreactions,namely the 32species and255 reactions set developed by Liet al.[27]. Thisdetailedmechanismdoeshowevernotincludetar decompo- sition,andthereforereactionNo.4inTable4from[13]isaddedto the mechanism. Toincrease the simulationeﬃciency, thetabula- tionofdynamicadaptivechemistry(TDAC)method[28]isused.

2.4. Solutionstrategy

The particle collision sub-model requires particle searching, whichisquitecomputationallycostly.However,inafixed-bedsim- ulation,theparticlesarealmostconstantlyinstagnantcontacts.It isthusnotnecessarythattheparticlecollisionmodelissolvedat everytime step.Hence,onlyataspecifiedfrequencyistheparti- clemotionsub-modelthattakesadvantageofthispseudo-steady- statebehaviorapplied.Ineverymfluidtimesteps,theparticlemo- tion sub-model is solved at the first n time steps (m and n are arbitraryintegers,wheren<m).Forexample,inthisstudy,inev- ery100fluidtime steps,theparticlemotioncanbe resolvedonly atthefirst 20time steps,then thecomputational expenseofthe particlemotionsub-model is reducedtoone-fifth. The fluidtime steps forwhich the motionis not calculated,the contactsof the particle andthe particle-particle contact radius r_{c,i j} are kept the sameasthe lastpredictions fromtheparticle motionsub-model, whileother particle properties,such as,the particletemperature, willstillbeupdatedbysolvingtheparticleconversionsub-models.

The valuesused inthisstudy arebased on oursimulation expe- rience. It isworth mentioningthat such a calculation strategy is specificallydesignedforthe solid phasewithslow motion.Ifthe particlesmovementisintense,forexample,inafluidizedbed,this strategyshouldnotbeadopted.TheXDEMcouldalsouseasmaller timestepthanthefluidtime step.Theoverallcalculationscheme isshownasfollows.

1. ResolvetheCFDmodel,updategasphasepropertiesatthecur- rentﬂuidtime.

2. Map the gas data to the particle grid and calculate the gas phasepropertiesasseenbyeachparticle.

3. Resolve particle conversion sub-models in one particle time step.

4. Resolveparticlemotionsub-model inoneparticletime step,if inthisﬂuidtimestepitisrequired.

5. Advancetothenextparticletimestep,andrepeatsteps 2and 3,untilthe accumulatedparticletime steps equal thecurrent ﬂuidtimestep,thentransferparticledatatotheﬂuidgrid.

6. Diffuseparticledatafromt=0sto

τ

totalsbyDBMontheﬂuid grid.

7. AdvancetheCFDsolvertothenextﬂuidtime.

(6)

Fig. 3. Computational domain and boundary conditions.

Table 1

Particle physical properties.

Value Unit

Moisture 9.8 w %

Ash content 3.3 w %(dry basis)

Ash intrinsic density 2000 kg ·m ⁻³ Equivalent diameter, d p 0.0084 m

Cylinder aspect ratio 1.67 -

Density of the particle 1030 kg ·m ⁻³ Shrinkage factor for char burnout 0.05-0.35 - Ash layer porosity calculated -

3. Validationandsensitivitystudy

Inordertovalidate thenumericalmodel,importantproperties predictedby simulation arecompared withthe measurements of Porteiro etal.[29].Theexperimentissetup asafixed-bedcombustion inatube reactor.Theignitionfront propagatesinthedi- rectionoppositetotheairflow,whichisknownascounter-current mode.Thefront ofconversionandignitionpropagatedownwards, which can be quantified by the measured ignited massflux. The simulationdomainisasectionoftheexperimentalsetup.Inaccor- dancewiththeexperiment,thegeometryandboundaryconditions areillustratedinFig.3.

Threedifferentsetsofﬂuidgridsaregeneratedforagridinde- pendencetest.Thethreegridshave79,920(ﬁne),23,040(medium) and6,144(coarse)cells,respectively.Thecorresponding meancell side lengthsare2.6mm,3.9mmand6.1mm,forthethreegrids.

Additionally, a coarser grid (PCMGrid) is generated for the PCM simulation.BecausePCMnormallyrequiresthatthecellshouldbe muchlargerthantheparticlesize,thePCMGridhas672cellswith ameancellsidelengthof12.7mm.

The thermophysicalpropertiesofthegasphase,aswell asthe calculation ofreaction heat, are evaluated by the standard NASA polynomials [30]. The cases using poplarpellets as fuel [29] are chosen forvalidation. Some relevantparticlepropertiesare listed inTable1.Otherpropertiesarethesameaspresentedinourpre- viouswork(Table5)[13].

In the simulation, 2000 poplar pellets were injected into the domainfromthetop,andthefuelbedwasthusformedbymeans ofa randomstackingofparticlesthrough DEMcalculation.Inthe experiments, the bed wasignited fromthe top. Correspondingly, theparticlesinthetoplayerofthefuelbedareinitializedaschar particles(particlesthathaveaccomplisheddryinganddevolatiliza- tion)withatemperatureof1250K,whichissufficienttostartthe charoxidation.Foralltheparticlesub-models,atimestepof10⁻⁴ sisused,whilefortheCFDsolverthetimestepis10⁻³s.Thepar- ticlemotionsub-modelissolvedinthefirst20timefluidstepsfor every100timesteps.

For the DBM, the diffusiontime

τ

total wassub-divided into6 shortertimesteps.Differentlengthscales(b)forthediffusionop- eration are adopted for different particle data. For the exchange source terms,b=dpis applied,because theinteractions between

Fig. 4. Instantaneous ignited mass ﬂux calculated using PCM and the proposed method with different grid resolution (data for Exp. is the mean ignited mass ﬂux).

The airﬂow rate is 0.276 kg ·m ⁻²·s ⁻¹, and the particle grid size is 1 . 5 d p.

the particleand thegas phase are expectedto occur in the particle’s local region, without being overlapped too much withits neighbors.While fortheparticle properties,likethe particlevol- ume,b=2d_pisapplied,becausetheDBMisexpectedtosmearthe particlepropertiesmoresmoothlyinspacetoincreasetherobust- nessofthesolver. TGMrequiresthat theconﬁgurationofthesize oftheparticlegridisdeﬁnedinthecouplingmethod.Asensitivity studywasconducted, andthe averageside length ofthe particle grid cell of1dp, 1.5dp and 2dp are compared. The defaultvalues fortheparticlegrid sizeandfortheshrinkagefactorduringchar combustionare1.5dpand0.3,respectively.

3.1. Gridsensitivityanalysis

The combined TGM-DBM coupling method presented in Section2.2onlyworkswhenthesizeoftheparticlegridislarger thanthefluidgrid.Thethreefluidgridswithdifferentresolutions aretestedagainstthePCMonthePCMgrid.ThegridsizeforPCM is veryclose to theparticle grid size.The transient ignitedmass flux isshowninFig. 4.The time-averagedignitedmassflux pre- dictedbyusingfine,mediumandcoarsegirdsare0.056,0.063and 0.064kg·m⁻²·s⁻¹,respectively,whichareingoodagreementwith experimentalobservation.Fortheinstantaneousignitedmassflux, thecoarsergridpredictedalargeroscillation.

As showninFig. 4,thePCM isfound to fail inpredictingthe ﬁxed-bed combustion. In the experiments, there will be an ignition front propagating downwards. Particles above the ignition front will release heat intensively due to the char burnout and the combustionof volatilesreleased in thedevolatilization stage.

InthePCM simulation,thecomputational cellsaretoolarge. The source termsin heat andspecies equationswill be evenly redis- tributedintoacorrespondingly largespace,resultinginanunder- estimationofthegasphasetemperature.Then,theparticlewillbe cooleddownbyconvection,andthebedwillbeunabletosustain thepropagationoftheconversionfront.Itshouldbe notedthatif theTGMisusedforboththesamplingandredistribution,similar resultswillbe obtained. Suchdeviationsindicate thenecessityof apropercouplingmethodinthesimulation.

Thepredictionsusingthe newcouplingmethodwithdifferent ﬂuid resolutions show a clear periodictrend in theignited mass ﬂux asshownin Fig.4. The resultshave a similar period,which isrelatedtotheparticlegridsizeaswillbediscussedinthelater sections.Since the TGM coupling is implemented withthe same couplingparameters for thedifferent simulations,the differences are mainly due to that the DBM coupling has different perfor-

(7)

Fig. 5. Performance of DBM on different grids. The particle grid is marked in gray, and the ﬂuid grid is marked in blue. (For interpretation of the references to colour in this ﬁgure legend, the reader is referred to the web version of this article.)

Fig. 6. The probability distribution of initial heights of particle centroids.

mances onthe grids withdifferentresolutions. Fig.5 showsthat thegridresolutionhaseffectsontheDBMcouplingtakingtheen- ergy sourceterm(S_h) asanexample.Aparticlewhosecentroidis locatedattheedgeofacellboundaryoftheparticlegrid,willhave quite a different impacton the neighboring particlecells than if its centroidis placedcloserto thecenterofthegrid cell.Witha finer fluid grid, theparticle data are smearedmore smoothly after the DBM has been applied, while with a coarser fluid grid, the particle’s influence isstill mainlyon the particle grid,which holds theparticle’s centroid.Ingeneral, a finerfluid gridwillre- duce thegradientofgasphasepropertiesontheparticlegrid.So, when the fluid grid becomes finer, the curves in Fig. 4 become more smooth. It is worth noting that a mesh-independent solu- tion doesnot meanthat thefluctuationsareeliminated. Thepar- ticle mass lossrate isvery uneven, especially regarding different conversionsub-processes.So,suchfluctuationscouldarisebecause of the discrete nature of the XDEM representation of the particles in the bed, andthus reflect a correct physicalphenomenon.

Forexample,whenthewidthofthefixed-bedisquitelimited,the packing of the particles with uniformed diameterwill show the layering structure,asshowninFig.6.Then,thequasi-steadycon- version willshow a periodiccharacter astheignitionfront prop- agatinglayerbylayer oftheparticles. Inthesimulation,the finer grid also requires a larger computational cost, especially for the DBMcoupling.Consideringalltheseeffects,themediumfluidgrid ischosenfortheremainderofthisstudy.

Fig. 7. Instantaneous ignited mass flux predicted by different sizes of particle grid (data for Exp. is the mean ignited mass flux). The airflow rate is 0.276 kg ·m ⁻²·s ⁻¹, and the fluid grid is medium.

3.2. Couplingparametersstudy

FortheDBMcoupling,theparticleisalreadymuchlargerthan the ﬂuid grid cell and the major concern is therefore that the smoothingofthesourcetermsshouldnotcauseexcessivenumer- ical diffusion. Large positive source terms will exacerbate solver robustness [31],thus there are usually strong arguments for ex- tensivesmoothingfromastabilityperspective.Asdiscussedabove, duetothattheparticlesareindirectcontactwitheachother,the particle’sdirectinﬂuencewillbe limitedtotheregionclosetoits surface.Thereisnotmuchfreevoidspacetomotivatechangingto adifferentlengthscaleforthesourcetermsotherthanb=dp.As longastheparticlepropertiesarechangingsmoothlyinrelationto thegasphase,itwillalsoenhancetherobustnessofthesolver.

The primary use of the particle grid is used to calculate the means of thegas data asseen by theparticles. In the fixed-bed simulation,iftheparticlegridcell islargerthan2timesthe par- ticlediameter,theinfluencefromtheneighboringparticleswillbe largerthanwhatistypicallypossibletomotivatefromphysicalar- gumentsforthetargetedparticle.Thegasphasebeyondtheparti- cle’sfirstneighborsisunlikelytohavedirectinteractionwiththe particle.

Thepredictedtransientignitedmassﬂuxeswithdifferentparti- clegridsizesareshowninFig.7.Intheﬁrst300s,theresultusing 1dp particle grid sizehas 3periods, while theresult using1.5dp

particlegridhas2 periods.Theperiods areexactlyinverselypro- portionaltotheparticlegridsize,showingthatthechoiceofpar- ticlegridspacingwillsignificantlyinfluencethestatesoftransient simulation.However,thetime-averagedignitedmassfluxfromthe two predictions are quite close (0.63 and 0.60 kg·m⁻²·s⁻¹, for 1.5dpand1dpcase,respectively),whileforthepredictionwith2dp

particlegrid,itseemsthattheconvectivecoolingisoverestimated inthesamewayasinthePCMpredictions.Inthisstudy,theparticle grid size of 1.5d_p was selected for further investigation, as its predictionagrees with theexperimental measurement. At the same time, thechosen grid size ismore in correspondencewith themodelassumptionthat theparticlegridsize shouldbe larger thantheparticlesize.

3.3. Sensitivitystudyonparticleshrinkagefactor

Unlike the particle shrinkage during the drying and devolatilizationprocess,thechangingvolumerelatedtocharconver- sionisassociatedwithalargeruncertainty.Itisnotonlyrelatedto theashcontentbutisalsoaffectedbythecombustiontemperature

(8)

Fig. 8. Mean ignited mass ﬂux predicted with different shrinkage factors for the char conversion process. The particle grid size is 1 . 5 d pand the ﬂuid grid is medium.

Fig. 9. The influence of shrinkage factors ( η= 0.05, 0.1, and 0.15) on instantaneous ignited mas flux at an airflow rate of 0.452 kg ·m ⁻²·s ⁻¹.

andeventheencapsulationofthematerial[32].Theshrinkagefac- tor forthe charconversion process (

η

⁾ ^will^eventually ^determine

the volume of the ashlayer andthe particle’s outermost surface area, whichis criticalinthe calculationofthe heattransfer with thegasphaseandneighboringparticles.Aseriesofshrinkagefac- tors forthecharconversionprocessare testedunderfouroperat- ing conditions,whichareoxygenlimited,reactionlimited,cooling byconvectionandblown-off conversionregimes.

The predictionsofthemeanignitedmassﬂuxestogether with the experimental measurements are shown in Fig.8. The overall trend isthat thebedconversionrateincreaseswiththedecrease of

η

^.^Because

η

ônlyâffects^the^volumeôf^theâsh^layer,ît^will^have

noinﬂuenceonthecharlayersurfacearea.Bydeﬁnition,asmaller

η

^means^a^smaller^eventual^particle^volume.^During^char^burnout,

a smalleroutermost surfacearea willreduce the coolingthrough convective heat transfer with the gas phase. Then, the particle’s surface temperature will be overestimated. In Fig. 8, for

η

^equal

to 0.05, 0.1 and0.15, thepredictions at the highestairﬂow rates are absent,becauseinthesethreecasesthesimulationsarequite unstable,whichmeansthattheperiodicalchangingoftheignited massﬂuxisnotestablished,asshowninFig.9.

For

η

=0.15,thebedtemperaturesat50s,150sand250sare shown inFig. 10.The bedis not totallyextinguished atthe high airﬂowrate.Alongthewallregionwherethegasvelocityisrather low,thebedconversionisstillself-sustained.Asmall

η

^also^helps

thecharparticle,whichisatthelastconversionstagewithasmall volume,ﬁllintothegapsinthenext particlelayer.Moredetailed

Fig. 10. Particle temperature at 50, 150, and 250 s after the ignition of the bed ( η

= 0.15).

experimental observationisneededto validate thisphenomenon.

Johanssonetal.[33] alsopointedout that theash propertiesare critical in determining the burning velocity and extinction phe- nomenaathighairﬂowrates,byaffectingtheradiativeheattrans- ferinthebedintheirsimulationofﬁxed-bedcombustionthrough a PMMapproach. The sensitivityanalysison

η

^indicates ^that ^the

modelofthemorphologyoftheashlayercouldplayanimportant roleinsuch XDEM/CFDmodelingontheﬁxed-bedcombustion.A shrinkage factor between0.25 and0.35 results in a good agreementwiththeexperimentalmeasurements.Inthisstudy,

η

=0.3 isadoptedinthelatersimulations,andtheﬁnalvolumechanging isinareasonablerangeconsideringtheashcontent[34].

4. Resultsanddiscussion 4.1. Ignitedmassﬂux

Fromtheabovesections,thefollowingsimulationswereconﬁg- uredwiththemedium ﬂuidmesh,1.5dpparticlegrid,and

η

=0.3.

The predicted mean ignited mass fluxes at all measured airflow ratesareshowninFig.11.Theexcessairratiosarecalculatedwith regard to thestoichiometric air(5.59 kgof dryairper kg offuel burnt) which is determined by the chemical equivalent formula ofthe poplar pellets(CH₁_.₅₅O₀_.₇₅, majorelements only) [29]. The predictions by the developedXDEM/CFD modelare in very good agreementwiththeexperimentalmeasurements.Theslopesinthe three combustion zones, which include oxygen limited, reaction limited andcooling by convection, are well captured, as well as theextinctionatahighairflowrate.

In the reaction limited zone, simulations predict a small in- creaseinthemeanignitedmassflux,whichissignificantlylower thanintheoxygen-limitedregime,butsomewhathigherthanthat observedintheexperiments.Welooked furtherintothesedevia- tionsandobservedthattheyoriginatefromtheTGMcouplingand theparticlesurfacemasstransfersub-model.Firstly, themeanO₂ concentration calculatedontheparticlegrid,whichisneededfor the particle conversion sub-model, is not equal to that infinitely faraway fromtheparticlesurface, asisassumedinthetheoreti- cal underlyingderivation ofthe sub-model[13].This factimplies thatthefar-fieldconcentrationemployedinthesub-modelwillbe lowerthan thatattheinlet,andwillincrease withthe inletflow rate. Secondly, compared to the mass diffusion through the ash layer,themassdiffusionfromthegasphasetotheparticlesurface isthedominatingprocess.Themasstransfercoefficienthm calcu- latedfromEq.(12)isincreasing asthegasvelocityincreases.The charconversionrateisdirectlydeterminedbytheabovetwovari- ables [4]. Still,it isconcluded that thedeviations inthe reaction limitedzoneareacceptable.

4.2. Thecouplingeffects

The instantaneous ignited mass ﬂux by process contribution predictedinthecasewheretheairﬂowrateis0.327kg·m⁻²·s⁻¹ (referredtoasthereactionlimitedcaselater)isshowninFig.12.

(9)

Fig. 11. Mean ignited mass ﬂux at different airﬂow rates (a) and excess air ratios (b).

Fig. 12. Instantaneous ignited mass ﬂux including the contribution from each sub- process at an airﬂow rate of 0.327 kg ·m ⁻²·s ⁻¹(reaction limited case).

The majorityofthemasslossisthrough thedevolatilizationpro- cess. Ineach givenperiod,theconversion ratesofthedryingand char burnout processes change oppositely to the rate of the devolatilization process. Such periodic oscillationis strongly related to theparticle gridconﬁguration. Fromthe bottom ofthebedto the height of 0.06 m, there are 5 layers of particle grid cells in theverticaldirection.Thethermalconversionprocessinthesegrid layers isillustratedinFig.13intermsofparticletemperatureand conversionratio,respectively.

ThecoloredregionsinFig.13marktheboundariesoftheparti- cle’sstatusdistributionatcertainpointsintime.Thelayersarere-

ferredtoas1stto5thlayerparticlecellsforconvenience.At50s, theconversionrateofdevolatilizationisincreasing,whilethedry- ing and charburnout rates are decreasing. Fig.13 showsthat at 50s,inthe4thparticlecellsattheheightaround0.04-0.05m,the temperaturedifferencebetweentheparticles canindeedbe quite large. Except forthe charparticles, mostparticles inthese particlecells justcompleted drying orare still inthe drying process.

Atthisstage, theparticles containmainlythe drywood.There is acleardryingfrontattheheightofaround0.04mthattemporar- ilystoppedbetweenthe3rdand4thlayeroftheparticlecells.At 100s,theparticlesinthe3rdlayeroftheparticlecellsstartedthe dryingprocess.Thedrywoodinthe4thlayeroftheparticlecells isstill beingconsumed.When thedevolatilizationratedecreases, moreO₂ couldjointhecharburnoutinsteadofthehomogeneous gas phase reactions with the released pyrolysis gases. At 200 s, Fig.13 showsaclearstagnation ofthepropagationof thedrying frontagain.

Thesizeoftheparticlegridcellsisdesignedtobe largerthan thelength scaleoftheconversionfront, andasa resultthefront will eitherbe smoothed out ormove discretely ata length scale dictatedbythemeshsize.Thegasphasetemperatureiscriticalin thepredictionofthedryingprocess,andthenumericaleffectsin- troducedby the couplingschemeconsequently produce feedback effectsinthemodelingpredictionofthebedconversionprocesses.

Asa summary,thecoupling effectoftheTGM showsthatan ac- curate,mesh-independentpredictionofthelocationofthedrying frontishardtoachieve.

4.3. Heatcontributioninﬁxed-bedcombustionofbiomass

The particles whose initial centroid positions are located be- tween0.039mand0.042m(inthe4thlayeroftheparticlecells, 340 particles) are tracked for their entire conversion process in an effort to quantifyparticle-level variations, something that the XDEM/CFDapproachintrinsicallyincorporatesbutwhichisﬁltered outinEulerianmethods.Foreveryparticle,theparticleconversion time count startsonly atthe onset ofthermal conversion(when conversionequals 0.005).Forthereaction limitedcase,the parti- cleconversionratioandsurfacetemperatureversustheconversion timeareshowninFigs.14and15.Theshadowedregionshowsthe boundariesofthedistribution,andtheevolutionofthreerepresen- tativeparticlesisalsoplottedassolidlinesintheﬁgures.

Particles A, B and C are identiﬁed to show different conversion histories. They can be described as a fast converting particle, medium converting particleandslow convertingparticle, respectively. The maindifference is observed inthe dryingprocess because for the devolatilization and char burnout processes, the curves inFig. 14 show similar slopes. Fig. 15 reveals that in the dryingprocessdifferentparticlesexperienceratherdifferentheat- ingrates.Afterthedevolatilization, whenthecharburnoutstarts, there is a sudden temperature jump. The temperatures remain similarduringthecharburnoutphaseforthedifferentparticles.

TheheatcontributiontermsinEq.(3),togetherwiththecom- bustionheat,whichisreleasedfromthecharburnoutprocess,for particlesA,BandCareshowninFigs.16,17and18,respectively.

Thefastconvertingparticleisdirectlycontactedtothe hightem- perature charparticles. The particleis heatedup mainly through particleto particle radiation,whiledueto the freshairﬂowfrom the bottom, the convection by the gas acts ascooling the particle. Forthemedium converting particle,itis not directlyincon- tact withany charburning particles, butit belongs to a particle cellthatincludeshigh-temperatureregions.Convectionandradia- tionhavesimilarcontributions totheparticletemperaturerise.It isshownin Fig.17,particleB experienced areignition ataround 160 s. Thisis because, atthis moment,particle B moved from a particle grid cell in the 4th particle cell layer to another one in

(10)

Fig. 13. Particle temperature (a) and conversion (b) along with the bed height at selected points in time for the reaction limited case. Two vertical dashed lines approximately separate the conversion into three stages: drying, devolatilization, and char burnout. The left vertical dashed line is based on the moisture content, and the right vertical dashed line is based on the average char yield.

Fig. 14. Particle conversion versus conversion time for the reaction limited case.

Fig. 15. Particle surface temperature versus conversion time for the reaction limited case.

the 3rd layer whichhas a higher O₂ concentration. For the slow converting particle, it is very likely that at the beginning ofthe period, it is in contactwithparticles that are going through the devolatilizationprocess. Thetemperatureoftheparticle undergo- ing devolatilizationismainlyintherangeof600-800K.Theheat- ing processisratherslow.However, whenitsneighborscomplete devolatilization, andstart char burnout, the slowconverting particle will start to behave like the fast converting particle.In the charburnoutstage,bothconvectionandradiationcontributetothe cooling ofthe particle.As showninFigs. 16,17 and18,theheat transferduetotheparticletoparticleconductioncanbeneglected forall threeparticles.Other studiesalsoindicate similarﬁndings,

Fig. 16. Heat ﬂow rate of different mechanisms for particle A.

Fig. 17. Heat ﬂow rate of different mechanisms for particle B.

asforexamplebyShinetal.[34].However,itisstillpossiblethat inthepresentstudytheparticletoparticleheatconductionisun- derestimatedbycurrentsub-models.

Thecasewiththeairﬂowrateof0.126kg·m⁻²·s⁻¹ isreferred toasrepresentativeoftheoxygenlimitedcase.Acomparisonbe- tween the oxygen limited caseand the reaction limited case on theparticleconversionprocessisshownonFig.19.Alargepartof theparticleconversionprocessesareoverlapping,however,during thecharburnoutstage,theoxygenlimitedcasehasalowerignited massﬂux,andtheparticlesconversionratesaremuchlower.

The probability distributions of the heat contribution for the tracked particle are shown inFigs. 20–23. Figs. 20 and 22 show

(11)

Fig. 18. Heat ﬂow rate of different mechanisms for particle C.

Fig. 19. Particle conversion versus conversion time for the reaction limited case and the oxygen limited case.

Fig. 20. Contribution from different heat exchange mechanisms during drying and devolatilization for the reaction limited case.

theratiobetweenconvectionandradiationheatﬂowratesduring the drying anddevolatilizationprocesses forthe reaction limited caseandtheoxygenlimitedcase,respectively.Inbothcases,radiation isthemain contributorto heatsourceterms.Nearlyhalf of theparticleslostheatthroughconvection(whereqconv/q_radisneg- ative while q_rad is positive). Compared with the reaction limited case, theconvective heatingcontributesmoretotheparticletem- perature increase intheoxygen limitedcases.Forthe charburnout stage, the heattransfer fromthe particleto theenvironment

Fig. 21. Contribution from different heat exchange mechanisms during combustion for the reaction limited case.

Fig. 22. Contribution from different heat exchange mechanisms during drying and devolatilization for the oxygen limited case.

Fig. 23. Contribution from different heat exchange mechanisms during combustion for the oxygen limited case.

is mainly archived by convection and radiation. For the reaction limitedcase,theconvectivecoolingismuchstrongerthanthera- diativecooling,whichisquiteasexpected.

5. Conclusion

An XDEM/CFD model for the simulation of biomass combus- tionin aﬁxed-bed hasbeendeveloped,witha speciﬁc emphasis