CFD modeling of biomass combustion and gasification in fluidized bed reactors using a distribution kernel method

ContentslistsavailableatScienceDirect

Combustion and Flame

journalhomepage:www.elsevier.com/locate/combustflame

CFD modeling of biomass combustion and gasification in fluidized bed reactors using a distribution kernel method

Miao Yang

^a

, Jingyuan Zhang

^b

, Shenghui Zhong

^a,c

, Tian Li

^b,d

, Terese Løvås

^b

, Hesammedin Fatehi

^a

, Xue-Song Bai

^a,∗

aDepartment of Energy Sciences, Lund University, Lund 22100, Sweden

bDepartment of Energy and Process Engineering, Faculty of Engineering, NTNU - Norwegian University of Science and Technology, Trondheim, Norway

cState Key Laboratory of Engines, Tianjin University, 135 Yaguan Rd, Tianjin 300350, China

dRISE Fire Research, Tiller, Norway

a rt i c l e i nf o

Article history:

Received 5 May 2021 Revised 7 September 2021 Accepted 8 September 2021

Keywords:

Biomass combustion and gasification CFD simulation

Fluidized bed furnace MP-PIC

Distribution kernel method

a b s t r a c t

A three-dimensional reactive multi-phase particle-in-cell (MP-PIC) model is employed to investigate biomasscombustion and gasificationinfluidizedbed furnaces.The MP-PICmodel considered here is basedonacoarsegrainmethod(CGM)whichclustersfuelandsandparticlesintoparcels.CGMiscompu- tationallyefficient,however,itcancausenumericalinstabilityiftheclusteredparcelsarepassingthrough smallcomputationalcells,resultinginover-loadingofsolidparticlesinthecells.Toovercomethisprob- lem,inthisstudy,adistributionkernelmethod(DKM)isproposedandimplementedinanopen-source CFDcode,OpenFOAM.InDKM, aredistributionprocedureisemployedtospreadthesolidvolumeand sourcetermsoftheparticlesintheparceltothedomaininwhichtheparticlesareclustered.Thenu- mericalstiffnessproblemcausedbytheCGMclusteringcanberemediedbythismethod.Validationof themodelwasperformedusingdatafromdifferentlab-scalereactors.Themodelwasshowntobeable tocapturethetransientheattransferprocessinalab-scalebubblingfluidizedbedreactorundervarying fluidizationvelocitiesandloadsofsand.Then,themodelwasusedtostudythecombustion/gasification processinabubbling fluidizedbedreactorunder varyingambient temperatures,equivalentair ratios, andsteam-to-biomassratios.TheperformanceofDKMwasshowntoimprovetheaccuracyandthero- bustnessofthemodel.

© 2021TheAuthor(s).PublishedbyElsevierInc.onbehalfofTheCombustionInstitute.

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

1. Introduction

Biomassisconsideredapromisingenergysourcebecauseofits worldwide availability, ease of access, and renewable generation within a short period[1,2]. Biomasscan be usedin differentap- plications, such asheat orpower generation, chemical synthesis, andproductionofnanomaterials.Biomasscanbetransformedinto liquid, gaseous, and solid fuels through different chemical, phys- ical, and biological conversion processes [3]. As an alternative to fossil fuelforpower generation andheating, biomass energycan contribute significantly towards the objectives of the UN’s Paris Agreement in reducing greenhouse gas emissions. Fluidizedbeds (FB) havebeenadopted incoal/biomassgasificationandcombus- tion duetoits highefficiencyforgas-solidcontactandadvantage

∗Corresponding author.

E-mail address: Xue-Song.Bai@energy.lth.se (X.-S. Bai).

ofcontinuousoperation[4].Multi-scaleandmulti-physicochemical processes,suchascomplexhydrodynamicsofdensegas-solidflow, particlecollision, heatandmass transfer,radiation,homogeneous andheterogeneouschemical reactions,andturbulent combustion, occur simultaneously in a FB furnace [5]. Many experimental methods have been used to reveal the working mechanisms of FB[6,7].Theexperimentalmethodshavethedisadvantageofhigh cost and long research cycles. Thus, the computational fluid dy- namic(CFD)approachisconsideredasanefficientmethodforin- vestigationofthecomplexgas-solid two-phaseflowandcombus- tionprocess[8].

Twomain approachesare adoptedtomodelgas-solidflowus- ing CFD, the Euler–Euler approach and the Euler–Lagrange ap- proach.IntheEuler–Eulerapproach,bothgasandsolidphasesare considered as the continuous phase while in the Euler–Lagrange approach,thegas phaseis consideredasa continuousphaseand solid phase asa discrete phase [9]. Inmodeling fluid-solidinter- action,thetwo-fluid model(TFM),developedbasedontheEuler–

https://doi.org/10.1016/j.combustflame.2021.111744

0010-2180/© 2021 The Author(s). Published by Elsevier Inc. on behalf of The Combustion Institute. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ )

Euler framework [10], has been extensively employed due to its low computationalcost.Themodelishoweverunabletodescribe the properties ofthe solid-phaseon the individual particle level, which makes itdifficult to considerthe distributionof theparti- cletypesandsizessince separatetransport equationsneed tobe solved foreachsizeandtype.Havingmanyvolumeaveragingand interpolationmethodsforbothgasandsolidphasecanleadtoac- cumulatednumericalerrorandcanresultinamesh-dependentso- lution.

In the framework of Euler–Lagrange approach, two groups of models are developed, discrete element method (DEM) [11] and multi-phaseparticle-in-cell(MP-PIC)method[12].Comparedwith theEuler–Eulerframework,theEuler–Lagrangeapproachcantrack eachparticleindividuallyandthepropertiesoftheparticles(diam- eter,density,velocity,temperature,andchemicalcomposition,etc.) can be takeninto account with a high accuracy [13]. This, how- ever,leadstoasharpincreaseinthecomputationalload,especially whenthecollisionbetweentheparticlesisconsidered.IntheDEM model, the collision forces foreach individual particle are calcu- latedbasedontheinteractionbetweenalltheindividualparticles inthesystem.Themodelhasahighaccuracy,butitisbeyondthe computationalcapacitytomodelindustrialscalefluidizedbedfur- naceswherequadrillionparticlesormoreareinvolved.TheMP-PIC model eliminatesthedifficulties ofcalculating interparticle inter- actions by mapping particlepropertiesto an Eulergrid andthen mapping back computedstress tensors to particle positions [14]. This method was first developed by Andrewsand O’Rourke [15]. MP-PIC is computationallymuch moreefficient than theclassical DEM methodsincemoreefficientcollisionmodelsareusedinthe simulationsandamuchlargertimestepisallowed[9].Recently,a generalmodelforthenumericalcalculationofcollisionalexchange ofmass,momentum,andenergybetweenparticleswaspresented [16],whichusesaBhatnager,Gross,andKrook(BGK)model[17]in a transport equation forthe particle distribution function, inor- der toapproximate the ratesat which the collisionsbring about thelocalequilibriumofparticleproperties.IntheBGKmodel,the effect ofcollisionsis representedby a simple relaxationtermon the right-hand side of the Boltzmann equation. A series of arti- cleshavebeenpublishedreportingimprovedBKG-modelcollision dampingtime forparticlevelocity fluctuations[18]andtheaddi- tional effect that drives the particle velocity distribution toward isotropy[19].

In an industrial-scale fluidized bed furnace, it is particularly important to reduce the computational loaddue to thepresence of a large number of particles in the system. The coarse grain method (CGM), in which thevirtual particles (known asparcels) are used to represent real particles, is widely employed in the Euler–Lagrangeframework[20].Thismethodwasdiscussedinde- tail in our previous work [21] that employed the CGM to simu- latelab-scalefluidizedbeds.OneofthecharacteristicsoftheCGM methodisthehighlocalcarrierloadincertaincomputationalcells since the fuel and sand particles are clusteredinto parcels. This maycauseproblemsinsimulatingindustrialfluidizedbedfurnace withverycomplexgeometrieswhichleadstothepresenceofsmall computational cells,andawiderangeofparticlesizedistribution (PSD). Numerical instability andnonphysical over-load can occur when large virtual particles pass through small size cells. Large particles contribute to large source terms in the gas-phase gov- erning equationswhich causeconvergence problemornumerical artifactsdueto largespatialgradients occurringinthesmallsize cells[22].Thisissue isparticularlyimportantin themostwidely used and low-cost particle centroid method(PCM), in whichthe entirebodyoftheparticleisassumedtobeinthelocalcellwhere the particlecentroidis.Athighparticleloadconditions,thesolid volume fractioninthesmalllocalcell canbelargerthanthatthe cellcanphysicallypermit.

Different methods have been proposed to deal withthe situ- ation when the ratioof themesh size to the particlesize is too small,includingthedivided particlevolume method(DPVM), the cube averaging method (CAM), the two-grid method (TGM), and the diffusion-basedmethod (DBM) [22,23]. In ourprevious work [24],threeofthesemethods,i.e.,CAM,TGM,andDBMwereinves- tigatedandtheirimpactonthesourcetermdistributionwaseval- uatedinasingle-particlecombustioncase.IntheDBM,thesource termsofaLagrangianparticlearedistributedintoanEulerianfield accordingtoastatisticalkernelfunction.Thesourcetermsarecal- culatedusingthePCMmodelbeforethedistribution,sincethegas propertiesrequired by theparticle sub-modelsare sampledfrom the particles’ local cell. The DBM is more robust compared with the other investigatedmethods; however, thecomputational effi- ciencydecreasesrapidlywiththeincreaseofthemeshresolution.

IntheTGM,avirtualcoarsegridiscreatedbasedonthefinegrid tosolve theparticlepropertiesandthesourcetermsare mapped tothefinegrid[25].TheTGM,whichtreatstheparticlesinalocal grid, maycausesignificant errors for thesimulations of particle- laden flows. In the CAM, a virtual cubic region is created as an interaction media between the particle andthe gas phase. Com- paredwiththeTGMandtheDBM,theCAMhasanobviousadvan- tage in computational efficiency for dense multiphase flow with unstructuredmeshes. However, two independentmeshesneed to beconstructedintheCAM,whichincreasesthecomplexityofthe implementationandparallelcomputation.

The purpose ofthis work isto develop a comprehensiveMP- PIC model taking into account chemical reactions to investigate theperformanceofdensegas-solidfluidized bedreactors.Thefo- cus is on developinga robust method that can handle the local over-loadprobleminsmall-sizegridcells.InSection2wepresent thebasic mathematicalframework ofthe Euler–Lagrangegovern- ing equationsandthe mathematical formulationforthe DKM. In Section 3 we introduce the numerical methodology to solve the Euler–Lagrangegoverning equations,aswell asthecell searching strategyandparallelcomputationmethodfortheDKM.Thevalida- tionofthemodelandthe performanceoftheDKMare discussed inSection4andtheconclusionsarepresentedinSection5.

2. Mathematicalformulation

In the MP-PIC approach, the governing equations of the con- tinuous phase and the discrete phase are describedin the Euler and Lagrange frameworks, respectively. The interaction between thediscreteandthecontinuousphaseismodeledusingthemass, momentumandenergy sourceterms.The massconservationand energyconservationequationsareadoptedfromRefs.[26–28]and the momentum equations are adopted from Refs. [28,29]. The physicalinterpretationofthe termsinthe governingequationsis not discussedindetail,butthe readerisreferred totheprevious publicationsforfurtherdetails.

2.1. Gasphasegoverningequations

Innumericalsimulationsoffluidizedbedfurnaces,Reynoldsav- eragedNavier–Stokes(RANS)turbulenceclosureisoftenused.The gas phase governing equations consist of the Reynolds-averaged continuity, momentum, energy and species transport equations.

Thecontinuityequationis

∂ α

^g

ρ

g

∂

^{t +}

∇

^·

α

^g

ρ

gug

=Sm, (1)

wheretheoverbardenotesthatthequantityisReynoldsaveraged, andthetildedenotesFavreaveraged.

α

^{gisthe gas}volume fraction,

ρ

^{gis the gas density,ug}is the velocity vector of the gas,andSm

represents thegas formationratedueto thermochemicalconver- sionofthefuelparticles.Themomentumequationsare

∂ α

g

ρ

gug

∂

^{t +}

∇

^·

α

^g

ρ

gugug

=

α

^g

ρ

gg−

α

^g

∇

^pg+

∇

^·

( α

^g

τ

^g

)

^+Su,

(2) where g is the gravitationalacceleration, pg is the gas pressure,

τ

^{gisthe sum}of viscous stressand Reynolds stress, andSu is the sourcetermduetomomentumexchangebetweenthegasandthe solidphase.Theenergyequationis

∂

(

^αgρg

(

^h+K

) )

∂^t +

∇

^·

α

^g

ρ

gug

h+K

=

α

^g

ρ

gug·g+^∂α_∂^g_t^pg +

∇

·

α

^g

ρ

g

ef f

∇

^h

+

α

^gQ˙r+

α

^gQ˙com+Sq, (3) where h is thespecific enthalpy ofthe gas, andK is the kinetic energyofthegasflow.e f f isthesumofmolecularandturbulent heatdiffusioncoefficients,e f f=g+

μ

^t/(

ρ

gPrt),wherePrt isthe turbulent Prandtl numberand

μ

^{t isthe turbulent eddyviscosity.}

Q˙_risthemeansourcetermduetoradiativeheattransfer,Q˙_comis the mean source termdue to volatilechemical reactions, and Sq

isthemeansourcetermduetothermochemicalconversionofthe solidfuel.Thespeciestransportequationis

∂ α

g

ρ

gYg,k

∂

^{t +}

∇

·

α

g

ρ

gugY_g,k

=

∇

·

α

g

ρ

gD_{e f f}

∇

^Yg,k

+

α

g

ω

˙g,k

+SYk, (4) where Y_g,k is the mass fraction of species k in the gas mixture, and

ω

˙g,k isthe meanchemical reaction rateof speciesk. D_{e f f} is the effective mass diffusion coefficient for species k taking both the viscous andturbulentcontributions intoaccount, D_{e f f}=Dg+

μ

^t/(

ρ

gSct)^,whereSct istheturbulent Schmidtnumber.S_Y

k isthe mean formationrateof speciesk dueto thermochemicalconver- sionofthesolidfuelparticles.

A Partially Stirred Reactor (PaSR) model is used to account for turbulence chemistry interaction when computingthe source termsduetogasphasechemicalreactions(

ω

˙g,k,Q˙_com)[30].Inthe PaSRmodel,themeanreactionratesaremodeledas

ω

˙g,k=

κ ω

^˙g,k

(

^Y,T,p

)

^{, (5)}

where

κ

^{isthevolumefractionofthereactivemixture,}

κ

=_τc+^τc_τm.

τ

^{c and}

τ

^{m denote the local} chemical reaction time andthe local mixingtime,respectively.Thechemicalreactiontime,

τ

^c,isdeter-

mined fromthemeanreactionratesofthefuel(

ω

˙f)andtheoxi- dizerorthegasificationagents(

ω

˙^o),

1

τ

^{c =max}

{

⁻

ω

^˙f

Y_f ,−

ω

˙^o

Yo

}

^{, (6)}

wheresubscripts f andodenotethefuelandoxidizerorthegasi- ficationagents,respectively.Themixingtime

τ

^mismodeledas

τ

^m=Cmix

ν

ε

^{, (7)}

whereC_mix,isamodelconstant(C_mix=1.0inthisstudy).

ν

^and

ε

denotethekinematicviscosityandthedissipationrateofturbulent kineticenergy,respectively.

Thestresstensor

τ

^{gintheEq.(2)isthesumoftheviscousand}

Reynoldsstressesandcanbewrittenas

τ

g=

τ

l+

τ

t. (8)

ThestresstensorforaNewtonianfluidisexpressedas

τ

l=

μ

g

(( ∇

^ug

)

+

( ∇

^ug

)

^T−2

3

( ∇

·ug

)

^I

)

, (9)

andtheReynoldsstressesaremodeledaccordingto

τ

^t=

μ

^t

(( ∇

^ug

)

+

( ∇

^ug

)

^T−2

3

( ∇

·ug

)

^I

)

−2

3

ρ

gkI, (10) where

μ

gisthedynamicviscosity,andIisthesecond orderunit tensor.Standardk−

ε

^{modelisusedtodeterminetheeddyviscos-}

ity,

μ

^t=

ρ

gC_μk²/

ε

^{,wherek istheturbulentkineticenergy.kand}

ε

^{aremodeledusingthefollowingtransportequations:}

∂

(

^αgρgk

)

∂^t +

∇

^·

α

^g

ρ

gugk

=

∇

^·

α

^g

μ

^g+μ_σk^t

∇

^k

+

α

gPk−

α

g

ρ

g

ε

+Sk, (11)

∂

(

^αgρgε

)

∂^t +

∇

^·

α

^g

ρ

gug

ε

=

∇

^·

α

^g

μ

^g+μ_σε^t

∇ε

+

α

^gεk

C_ε1P_k−C_ε2

ρ

g

ε

+S_ε, (12) whereP_k=

τ

^t:

∇

^{ugistheproductionrateofturbulentkineticen-}

ergy.S_k andS_ε are thesourcetermsduetogas−solidinteraction.

Standardvaluesofmodelconstantsareused,C_μ=0.09,C_ε1=1.44, C_ε2 =1.92,C_σk=1.0andC_{σ ε} =1.3[5,31].

2.2. Solidphasegoverningequations

Two types of particles, biomass and sand, exist in the solid phase.During thethermochemical conversionprocess of biomass particles, the sand particles are assumed to be chemically inert.

The interactionbetween theparticles andthe surroundinggas is through mass and momentum exchange and heat transfer. The mass, momentum, and energy conservation equations for solid phaseintheMP-PICmodelaredescribedinthefollowing.

2.2.1. Massconservationequationandpyrolysismodels

Biomass particle conversion (direct combustion/pyrolysis and gasification)can beseenasatwo-stageprocess:pyrolysis(or de- volatilization) andheterogeneous conversion of char [32,33]. The pyrolysis product consistsof heavy hydrocarbon species (suchas tar), light hydrocarbon species (such as methane), water, carbon monoxide,andcarbondioxide,etc.Pyrolysismodelsincludingde- tailed tar specieshavebeen reported(e.g.,[34]) andused inour previousstudiesofsingleparticlecombustion[35].Theaimofthis work is to develop a robust MP-PIC model that can be used to studya wide rangeof particle loads.Hence, a one-step pyrolysis model and a simplified homogeneous volatile gas and heteroge- neouscharreactionmechanismareemployedinthepresentwork, followingtheliterature[5,8,21,31,36–40].However,moreadvanced and complex chemical kinetic models (e.g., including tar chem- istry) [34,41,42]can beemployed in thedevelopedmodel frame- work.Theone-steppyrolysisreactionmodelis

Biomass→x1CO+x2CO2+x3H2O+x4H2+x5CH4

+x6Ash_(s)+x7C_(s), (13) wherex_j arethe stoichiometric constants. Ash_(s) andC_(s) denote respectivelyashandcharthatareinsolidphase.

Heterogeneous reactions of char with the surrounding gas species (such as O2, CO2, H2O) is a complex process, involving char-O₂,char-CO₂ and char-H₂Oreactions. High molecular weight hydrocarbons(tar) are treatedasunstable products andreactions withsulfurandnitrogenarenottakenintoaccount[5,36].Table1 lists the homogeneous and heterogeneous reactions considered, wherethe chemicalkinetic rateconstantsare takenfromthelit- erature.

Themassconservationequation fortheithbiomass particleis writtenas

dmi

dt =m˙i=m˙_vapor,i+m˙dev^ol,i+m˙char,i, (14)

Table 1

Homogeneous and heterogeneous reactions considered in biomass combustion and gasifica- tion. Note: C ₍s)is solid phase char. C krepresents molar concentration of gas species k .

Reference Homogeneous reactions Kinetic rate [Kmol/m ³/s]

R1 [31,36] CH 4+ H 2O → CO + 3H 2 R 1= 0 . 312 exp (−15 , 098 /T g)C CH4C H2O

R2 [31,36] CO + H 2O → CO 2+ H 2 R 2= 2 . 5 ×10 ⁸exp (−16 , 597 /T g)C COC H2O

R3 [31,36] CO 2+ H 2→ CO + H 2O R 3= 9 . 43 ×10 ⁹exp (−20 , 563 /T g)C CO2C H2

R4 [31,36] CH 4+ 2O 2→ CO 2+ 2H 2O R 4= 2 . 119 ×10 ¹¹exp (−24 , 379 /T g)C ⁰_CH^.2₄C ¹_O^.₂³ R5 [31,36] CO + 0.5O 2 → CO 2 R 5 = 1 . 0 ×10 ¹⁰exp (−15 , 154 /T g)C COC ⁰O^.2⁵C ⁰H^.2⁵O

R6 [31,36] H 2+ 0.5O 2 →H 2O R 6 =2 . 2 ×10 ⁹exp (−13 , 109 /T g)C H2C O2

Reference Heterogeneous reactions Kinetic rate [s/m]

R7 [40,43] C ₍s)+ 0.5O 2 → CO R 7 = 0 . 046 ×10 ⁷exp (−13 , 523 / (R uT i)) R8 [40,43] C (s)+ H 2O→CO + H 2 R 8 =1 . 71 ×10 ⁷exp (−211 , 0 0 0 / (R uT i)) R9 [40,43] C ₍s)+ CO 2→ 2CO R 9= 9 . 1 ×10 ⁶exp (−166 , 00 / (R uT i))

wherem˙_vapor,i,m˙_devol,iandm˙_char,idenotetheevaporationrate,the devolatilizationrateandthecharconversionrate,respectively.The moistureevaporationrateismodeledas[31],

˙

m_vapor,i= ShD_{di f f,va} di

P_sat,Ti

RuTi,s−X_v pg

RuTi,s

AsiM_v, (15)

whereShistheSherwoodnumber,whichismodeledusingRanz–

Marshall correlation,Sh=(²+0.6Re¹_i^/2Sc^1/3)^[44].Reiis Reynolds number based on the particle size andrelative velocity between the ith particleandthe surroundinggas.Sc is theSchmidt num- ber of the surrounding vapor. D_{di f f,va}, P_sat,T

i, As_i, T_i,s,X_v, and M_v represent respectively the vapor diffusion coefficient, the satura- tion pressure, the surface area of particle, the surface tempera- tureoftheparticle,themolarfractionofvaporinthesurrounding gas, andthe molar weight ofvapor. Ru is the universal gas con- stant.d_iisanequivalentsphereparticlediametercomputedbased on the particle mass m_i and a constant particle density

ρ

i, i.e., d_i=(^6mi/

πρ

i)^1/3.

Therateofdevolatilizationiscomputedbasedonthepyrolysis reactionmodel(Eq.(13)),

˙

mdev^ol,i=−Adexp

− E_d RuTi

m_volat,i, (16)

where A_d=5.0×10⁶ [s⁻¹] and E_d=1.2×10⁸ [J/Kmol] are rate constants[5],m_volat,i isthemass ofthevolatileremaining inthe particle.These rateconstantshavebeen validatedunderdifferent biomassdevolatilizationconditions[5,8,21].

The rateofcharconversion iscomputedbasedon thehetero- geneousreactions(R₇,R₈andR₉)listedinTable1,

˙ m_char,i=

3

j=1

˙

m_{char,i j}, (17)

where m˙_{char,i j} representthe char consumption rates by reactions withO₂(j=1,reactionR7),H₂O(j=2,reactionR8),andCO₂(j= 3,reactionR9),

˙

m_{char,i j}=−Asip_j Rdi f f,jRkin,j

R_{di f f,j}+R_kin,j, (18)

R_{di f f,j}=C_j[0.5

(

^Tg+T_i

)

^]0.75

d_i , (19)

R_kin,j=Ajexp

− E_j RuTi

, (20)

where R_{di f f,j}, R_kin,j, C_j, T_g, p_j, A_j and E_j represent respectively the diffusion rate coefficient, kinetic rate coefficient, mass diffu- sionrateconstant,gastemperature,partialpressureofthegasify- ing speciesinthegassurroundingtheparticleandArrhenius rate constantsforthecharreactionswithspeciesO₂,H₂OandCO₂ (cf., R7,R8andR9inTable1).C_j=5×10⁻¹²(s/K^0.75)[5].

2.2.2. Momentumequations

The kinematics of the ith particle simulated is governed by Newton’ssecondlaw,

dui/dt=

β

ρ

i

θ ⁽

^ug−ui

⁾

^+g

1−

ρ

g

ρ

i − 1

ρ

i

∇

^pg− 1

ρ

i

θ ∇ τ

, (21) whereu_i,

ρ

i,and

θ

^denoterespectivelythevelocityanddensityof the ith particle,and thesolid volume fractionin spatialposition x_i attime t.The right-hand side terms representthe sum of all forcesacting onthe ithparticle by thesurrounding gasand par- ticles.The forces considered include,from left to right, the drag, gravity, pressure gradient andthe interparticle stress.The virtual mass force,Basset force, andlift force such asthe Saffmanforce andMagnusforce,whichplayaminorroleinparticlemotion,are neglected[29].Withagivenu_i,thepositionoftheparticleiscom- putedbyintegrationoftheequation

dx_i/dt=u_i. (22)

InEq.(21),thesolidvolumefractionismodeledas

θ (

^xi,t

)

=

f

(

^mi,ui,xi,t

)(

^mi/

ρ

ⁱ

)

^dmidui, (23)

where f(^mi,u_i,x_i,t)^{is aparticle}distribution function,which de- scribesthestatisticaldistributionofmassandvelocityofparticles inspatialpositionx_iattimet.Namely, f(^mi,u_i,x_i,t)^dmidu_iisthe averagenumberofparticlesperunitvolumewithvelocitiesinthe intervals(^ui,u_i+du_i) ^{andmassin theinterval} (^mi,m_i+dm_i)^.In this model, it is assumed that the particle distribution function is independentof particletemperature andcomposition. Thisas- sumptionhasbeenwidelyusedintheliterature[14,31,45].

IntheMP-PICmodel f isobtainedfromtheLiouvilleequation, which isthe mathematical expression of conservationof particle numberspervolumemovingalongdynamictrajectoriesinthepar- ticlephasespace[15],

∂

^f

∂

^{t +}

∇

^·

(

^fui

)

⁺

∇

^u·

(

^fAi

)

^{= fG}

τ

^{−G f +}

fD−f

τ

^{D , (24)}

whereA_i=du_i/dtistheaccelerationoftheparticle, f_Gistheequi- librium isotropicparticle distribution function, f_D is the collision damping particle distribution function,

τ

G and

τ

D are relaxation times.

∇

^and

∇

^{uaredivergenceoperatorswithrespecttophysical}

spacexandvelocityu_i.InAppendixA,furtherdetailsaboutthese distributionfunctions,relaxationtimesandinterparticlestress are provided.

Thedragforcemodelwidelyusedfortheithindividualparticle isgivenbyKuetal.[5],Gidaspow[12],Yangetal.[40]

fi=Vi

β

θ ⁽

^ug−ui

⁾

^{, (25)}

whereV_iisthevolumeoftheithparticle.InEqs.(21)and(25),the dragforceparameter

β

^ismodeledas

β

⁼

150⁽¹⁻_α^α2^g)2μg

gd²_i +1.75⁽¹⁻_α^αg)ρg

gdi

|

^ug−ui

| α

^g<0.8

3

4C_d^(1−α_d^g)ρg

i

|

^ug−ui

| α

⁻g^2.65

α

^{g≥0.8 ,}

(26) whereC_disthedragcoefficient,modeledas[12]

C_d=

₂₄

Rei

(

¹+0.15Re⁰_i^.687

)

^Rei<1000

0.44 Re_i≥1000 , (27)

andtheparticleReynoldsnumberisdefinedas

Re_i=

α

^g

ρ

^gdi

|

^ug−ui

|

^/

μ

^{g. (28)}

2.2.3. Energyequation

Theparticletemperatureisobtainedfromtheenergyconserva- tionequationfortheithparticle,

m_iC_p,i^dT_dtⁱ=h_iAs_i

(

^Tg−Ti

)

^+εiAs4ⁱ

(

^G−4

σ

^Ti⁴

)

−h_vapor,im˙_vapor,i−h_devol,im˙_devol,i−₃

j=1hi,jm˙_{char,i j}, (29) whereC_p,i,

ε

i,

σ

^,hi=(^Nu

λ

g,s/d_i)^{,andGrepresenttheparticleheat} capacity, emissivity, Stefan-Boltzmann constant, interphase ther- maltransfercoefficient,andincidentradiation,respectively.h_vapor,i, h_devol,i,andh_i,jrepresentthelatentheat,theheatofpyrolysis,and heats ofchar reactionswithO₂,H₂OandCO₂,respectively. Nuis theNusseltnumber,whichismodeledusingRanz–Marshallcorre- lation,Nu=2+3/5Re¹_i^/2Pr^1/3 [44,46,47].PristhePrandtlnumber ofthesurroundinggas,and

λ

g,srepresentsthermalconductivityof thesurroundinggas.TheincidentradiationGisobtainedfromthe P-1radiationmodel.

2.3. Sourcetermsforparticle/gasinteraction

IntheCGMapproachafinitenumberofvirtualparticles(here- afterreferredasparcels)aresimulated.Assumethatthenumberof parcelsisNp.Theithparcelcontainsmultiplerealparticles; how- ever,all particleshavethe sameproperties,i.e., eachrealparticle in the ith parcelhas the samemass m_i,velocity u_i, temperature T_ianddiameterd_i.Thegoverningequationsfortheindividualreal particleintheithparcelhavebeenpresentedinSection2.2.

In thephysicalspacexattimet,thenumberofrealparticles per unit volume that pertain to the ith parcel is n_i. The source termsduetothegas/solidinteraction forthecontinuityequation, momentum equations, enthalpy equation, and the species trans- portequationsare,

Sm=−Np

i=1n_im˙_i, Su=−Np i=1n_if_i Sq=−Np

i=1niqi, SYk=−Np

i=1nim˙_k,i, (30)

where m˙_i=N

k=1m˙_k,i is givenin Eq.(14),f_i is givenin Eq.(25), q_i=m_iC_p,i^dT_dtⁱ isgiveninEq.(29),andm˙_k,iisduetopyrolysisreac- tion(Eq.(13))andcharreactions(R7,R8andR9,Table1).

2.4. Spatialredistributionofparcelsandsourceterms

InCGM-PCMapproach,thesolidvolumefractioninagivenpo- sitionandtimeiscalculatedusingEq.(23)withallparticlesinthe parcelstakenintoaccount,andthegasvolumefractionisthen

α

^g

(

^x,t

)

⁼¹⁻

θ (

^x,t

)

^{. (31)}

Physically,

θ

(^x,t)<1.However,since alargenumberofparti- clesareclusteredintoasingleparcel,whosevolumecouldexceed thevolumeofthelocalcell,thiscouldleadto

θ

(^x,t)>1or

α

g<0, whichisnon-physicalandnumericallyitcangiverise toinstabil- ity.Thisissueisparticularlytrueinthemostwidelyusedlow-cost

Fig. 1. Schematic illustration of PCM and DKM.

PCM.Inthepresentwork,weproposeanewmethod,theso-called distribution kernel method (DKM), which bears similarity to the DBMmethod,whilegivingtheadvantageoftheeasyimplementa- tionandlowcomputationalcost.AsshowninFig.1,theparcelsin a PCMcomputational cell (marked as cell-o)are clusters ofsand andbiomass particlesfromthesurroundingdomain(markedasa circularregion).Toavoidnumericalinstabilitycausedbylocallytoo manyparticles incell-o,in DKM theparticles andthe associated sourcetermsintheparcelsincell-oarere-distributedtothesur- rounding domain fromwhich the particles are clustered. The re- distributionalgorithmis constructedinsuch awaythatthesolid phase volume andsourceterms in theredistribution domain are conservedbeforeandafterdistribution.

A filtering kernel function g(^x,t)^{, which is defined based on} thedistanceofthesurroundingcellstocell-o,isemployedinthe presentwork.The integrationofthekernelfunction overtheen- tiregivenphysical spaceisunity. Similarstrategieswere used by Jesseetal.[48],Wangetal.[49],andSunetal.[22],aswellasin ourpreviouswork[24].Thesurroundingcellsofalocalcellcanbe locatedbyanewcellsearch algorithmbasedonagivendistance, asusedinthepresentwork.

Thetotal volume ofsolid phasein agivendomain ^isgiven by

Vs=

θ

^o

(

^x,t

)

^dV, (32)

wheresubscript “o” indicates thatthe quantity isbefore redistri- bution.Inthefollowing,subscript “r” willbe usedtodenote that the quantity isafter redistribution.Vs must stay thesame before andafterredistribution,

Vs=

θ

^o

(

^x,t

)

^{dV ≡}

θ

^r

(

^x,t

)

^dV=

g

(

^x,t

)

^{VsdV, (33)}

whichindicatesthat thesolid volumefractionafterredistribution is

θ

r

(

^x,t

)

=g

(

^x,t

)

^Vs=g

(

^x,t

)

θ

o

(

^x,t

)

dV, (34) Asimpleredistributionfunctiong(^x,t)^isemployed,

g

(

^x,t

)

=

(

¹−

|

^x−x0

|

dmax

)

², (35)

where x₀ is theposition of thecentroid of cell-o.dmax is a pre- scribed distance within which solid phase volume and source terms will be redistributed. The function g(^x,t) ^{maynot satisfy} Eq.(33). Bynormalizationofg(^x,t)^{,thefiltering kernelfunction} g(^x,t)^{canbeobtainedfromg}(^x,t)^,

g

(

^x,t

)

=g

(

^x,t

)

/

g

(

^x,t

)

^dV, (36)

Fig. 2. Cell searching strategy employed in the DKM model.

i.e.,

g

(

^x,t

)

^dV=1. (37)

g(^x,t) ^{can be used to} redistributethe source term for the mass conservationequation,

S_m,r

(

^x,t

)

=g

(

^x,t

)

^SM=g

(

^x,t

)

S_m,o

(

^x,t

)

dV, (38) whichisshowntosatisfymassconservation:

SM=

Sm,o

(

^x,t

)

^dV ≡

Sm,r

(

^x,t

)

^dV=

SMg

(

^x,t

)

dV, (39) Similarly, thesource termsforthe momentumequationsandthe energyequationandthespeciestransportequationscanberedis- tributedusingg(^x,t)^:

Su,r

(

^x,t

)

=g

(

^x,t

)

Su,o

(

^x,t

)

dV, (40) Sq,r

(

^x,t

)

=g

(

^x,t

)

Sq,o

(

^x,t

)

^dV, (41) SYk,r

(

^x,t

)

=g

(

^x,t

)

SYk,o

(

^x,t

)

dV. (42) 3. Numericalmethodsandcomputationcases

The MP-PIC and DKM model are implemented in an open- sourceCFDcode,OpenFOAMv6[50],basedoncoalChemistryFoam.

The MP-PIC modelisadoptedin coalCloudforthe discretephase while the solid/gasinteraction is taken intoaccount through the source termsinthegoverningequationsofthecontinuousphase.

Structured grid isused inthis study.The grid isgenerated using the“blockMesh” toolprovidedbyOpenFOAMpackage.

3.1. ImplementationofDKMandparallelcomputation

Inthe DKM,thesurroundingcellswithina certaindistanceto cell-o,d_max,areselectedforredistributionoperations.Thecellthat is partly located within the sphere of dmax is considered to be in the domain ofdistribution if the distancebetween the center of the givencell to the center of cell-o is less than or equal to dmax.Anefficientcellsearchalgorithmisneeded, since,thecom- putational cost wouldincrease hugely ifall cells arelooped dur- ing searching.Threesearch strategies,namelyshared-point-based, shared-edge-basedandshared-face-based,areusedinthepresent work. As an example,Fig. 2presentsa schematicdiagram ofthe searching procedure of the shared-face-based method. First, the neighboringcellsthatshareoneofthefacespertainingtothelocal cell areselected,cf.,theblackarrowinFig.2.Andthen,the blue

arrowsrepresentthesecondlayer,andthesearchingcontinuesun- tilthegivenmaximumdistancedmax,asshowninFig.2,isreached forallthelatestselectedcells.Aftereachsearch,theselectedcells thatareoutsideofthescopeofd_max areremoved.Notethat,ina staticmeshandafixeddmax,thesearchoftheneighboringcellsis onlyneededtobeperformedonce.

Since the searching strategy is cell-based without any spe- cificdirection,it worksperfectlyintheunstructuredmesh.How- ever,whenthesimulationdomainisdecomposedintoseveralsub- domains in a parallel computation, it is not straightforward to search theneighboringcellsacrosssub-domains.Thus,the whole simulationdomain isconsidered inthemaster processorandthe communicationbetween master andslave processors is achieved byapplyingthemessagepassinginterface(MPI).

As an example, asshown in Fig. 3,the domain of redistribu- tioninvolvesfoursub-domainsinfourslave processors.TheDKM procedureis doneasfollows.First,cell searchingisperformedin themasterprocessor. Thesolid volumefractionandsourceterms incell-oare transferredfromtheslave processor2to themaster processor. Then, the sourceterms of cell-owill be distributed to thecellswithinthesphereofdiameterdmax inthemasterproces- sor.Finally, thedistributedsourceterms arereturned tothe sub- domainsintheslaveprocessorsfromthemasterprocessor.

3.2. Numericalscheme

A finitevolume method is used forthe numericalsolution of the governing equations of the continuous phase. Central differ- ence scheme(CDS) is used forspatial directives and the Crank–

Nicolson scheme is used for temporal integration. PIMPLE algo- rithm, which combines the advantage of PISO (Pressure Implicit withSplittingofOperator)andSIMPLE(Semi-ImplicitMethodfor Pressure-Linked Equations)algorithms, isemployed to couplethe momentumequationsandthecontinuityequation.

3.3. Casesetup,initialandboundaryconditions

Twolab-scalefluidized bedreactors are employed to evaluate the MP-PIC model based on the DKM and PCM. The geometries of the two reactors are shown in Fig. 4. The first reactor, here- afterreferredtoasCase1,isabubblingfluidizedbed(BFB)reactor wheretheheat transferbetweenhotsandparticles andcold am- bientgaswasinvestigatedexperimentally[51].Thiscaseischosen hereto evaluatethe heat transfermodeland theperformance of DKMandPCM.Thesecondreactor,hereafterreferredtoasCase2, isbasedontheexperimentofbiomasscombustionandgasification ina lab-scalebubbling fluidized bedfurnace[45]. Inthe original experiment,thefurnacehasacylindricalgeometrywitha50mm internaldiameterand1200 mmheight. Inaprevious DEM simu- lation,thegeometrywassimplifiedtoacuboidwithanequivalent cross-sectionalarea [40].Thissimplified geometryischosen here toinvestigatethebiomasscombustionmodelandtheperformance ofDKM andPCM,which enableus tocompare directly withthe DEMresults.

3.3.1. Heattransferinlab-scalefluidizedbedreactor

Particle dynamicsandheattransfer inCase1are firstinvesti- gated.As showninFig.4,theBFBreactorhasadepthof15mm, widthof80mmandheightof250mm.ThefluidizationgasisN₂ with a temperature of 20 ^◦C. The gas is supplied from the bot- tomofthereactorwiththreedifferentsuperficialvelocities(U_sup).

Thereactorisinitially filledwithhotsandparticles of90^◦C.The sand particles havea uniform diameterof 1 mm. Table 2shows thephysicalpropertiesofthesandparticles.Sandtemperature(Ts) weremeasuredforfiveoperatingconditionswithacoldanodized aluminium background wall [51,52]. The total mass of the sand,

Fig. 3. Schematic diagram employed in the DKM parallelization.

Fig. 4. Schematic illustration of lab-scale bubbling fluidized bed reactors.

Table 2

Biomass and sand particle properties. ^∗The biomass mass of 0.002 kg is the total mass of biomass injected within 20 s. The biomass is injected to the reactor at a rate of 10 ⁻⁴kg/s.

Cases Particles d i[mm] ρi[kg/m ³] C p,i[J/kg/K] total mass ^∗[kg] T i[ ^◦C]

Case 1 [51] sand 1 2500 840 0.075 & 0.125 90

Case 2 [45] biomass 0.25 ∼0.35 750 1500 0.002 25

sand 0.25 ∼0.7 2300 840 1.25 750 ∼850