A numerical study on the combustion of a resolved carbon particle

ContentslistsavailableatScienceDirect

Combustion and Flame

journalhomepage:www.elsevier.com/locate/combustflame

A numerical study on the combustion of a resolved carbon particle

Ewa Karchniwy

^a,b,∗,∗

, Nils Erland L. Haugen

^c

, Adam Klimanek

^a

aDepartment of Thermal Technology, Silesian University of Technology, Konarskiego 22, Gliwice 44-100, Poland

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

cSINTEF Energi A.S., Sem Saelands vei 11, Trondheim 7034, Norway

a rt i c l e i nf o

Article history:

Received 16 August 2021 Revised 13 November 2021 Accepted 15 November 2021

Keywords:

Resolved particle Char conversion Overset grid

Mutliphase reactive flows Solid fuel combustion

a b s t r a c t

Combustionofasingle,resolvedcarbonparticleisstudiedusinganovelnumericalapproachthatmakes useofanoversetgrid.ThemodelisimplementedintotheframeworkofacompressibleDirectNumerical Simulation(DNS)code.Amethodtoartificiallyreducethespeedofsoundispresented.ForMachnum- berslowerthan∼0.1thismethodmaydramaticallyimprovenumericalefficiencywithoutaffectingany physicalaspects exceptfortheacoustics.Theability ofthemodel tosimulatesolidfuelcombustionis demonstratedandallpartsofthemodelarevalidatedagainstexperimentalandnumericaldata.Asensi- tivityofthecarbonconversionratetoselectedparameters(diffusioncoefficientsandhomogeneousand heterogeneouskinetics)isinvestigated.Astrongdependenceonthe oxygendiffusivityisobservedand explained.

© 2021TheAuthors.PublishedbyElsevierInc.onbehalfofTheCombustionInstitute.

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

1. Introduction

Solid fuels are among the most important energy sources worldwide. On one hand, some countries, like e.g. China, India or Poland, are still vastly dependent on coal [1]. On the other hand, the contributionto the energy productionfrom solid fuels in the form of biomass and refuse-derived fuel is increasing ev- eryyear[2].Duetoitsstrongeffectonglobalwarming,emission of carbon dioxide from solid fuels conversion is a serious envi- ronmental problem. This, in connection with the global increase in energy demand [3], necessitates developmentof low-emission andefficientsolidfuel-basedtechnologies.Such technologiescan- notbedesignedwithoutathoroughknowledgeaboutfuelproper- tiesandunderstandingoftheunderlyingfuelconversionphenom- ena.Thisunderstandingiscurrentlyprovided byexperimentsand by numerical simulations. Experimental investigation of solid fu- elscombustionisdifficultbecauseofcomplexphysicalandchemi- calprocessesoccurringatdifferentscales.Asaconsequence,infor- mation provided by experiments maynot be complete. Adeeper insight can be gained through detailed numericalsimulations, in which all flow scales are resolved on a numerical grid.It should bestressed,however,thatbothresearchmethodsarecomplemen- taryandequallyimportant.

∗Corresponding author at: Department of Thermal Technology, Silesian Univer- sity of Technology, Konarskiego 22, Gliwice 44-100, Poland.

E-mail address: ewa.karchniwy@sintef.no (E. Karchniwy).

InDirectNumericalSimulation(DNS)studiesonsolidfuelscon- version inturbulentsystems, particlesare commonlyrepresented aspoint sources.This approachhas previouslybeen employed to study differentaspects of pulverized coal combustion, for exam- ple in jet flames [4–6] and mixing layers [7–9]. The approxima- tionofpointparticlesisapplicableonlytoverysmallparticles,i.e.

toparticleswithdiameterssmallerthanKolmogorovlengthscales ofturbulence[10].Also,insuchsimulations,interactionsbetween thefluidandparticlesmustbemodeledusingclosureexpressions.

These expressions can be supplied by simulations in which the particlesurfaceanditsboundarylayerareresolvedonthenumeri- calmesh.Eventhoughsuchresolvedsimulationsaretypicallylim- ited to one or a few particles, thisapproach has a great poten- tialto providean understanding of thesolid fuel conversionand gas-particleinteractionsataveryfundamentallevel.Theresolved particleapproachhasrecentlybeenemployedinseveralnumerical investigationsofcoalorcarbonconversion.Devolatilizationandig- nitionstagesoftheresolvedpulverized coalparticlewereconsid- eredbyVascellari etal.[11],whosestudieswereextendedbyTu- fanoetal.[12]toaccountfordifferentatmospheresandamoreac- curatedescriptionofthevolatileyieldandcomposition.Thesame researchgroupfurtherbroadenedthefocusoftheirstudiesonre- solvedcoalparticlesbyconsideringparticlearrays[13],higherpar- ticleReynoldsnumbersandeffectsofturbulence[14].Anumberof publications neglect the devolatilization and investigate resolved char particlecombustion andgasification insteady state. Forex- ample,Kestel etal.[15]studied theimpact ofsteamcontentand

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

0010-2180/© 2021 The Authors. 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/ )

Reynoldsnumberonthecharoxidationinair,whiletheeffectsof the ambientgas temperature,gas velocity andoxygenmass frac- tionsinO₂/CO₂atmospherewereconsideredbyRichteretal.[16]. Asimilar analysiswasalsoperformedbySafronov etal.[17]who indicated differencesincombustionbehaviorbetweenmicro-and milimeter-sizedparticles.Theconversionofacollectionofresolved carbonparticleswasalsoinvestigatedinasimilarwaybySchulze et al.[18]. Furthermore,thesteadystate approach wasemployed inafewstudies [19–21]that attemptedtoresolveporousparticle andunderstandintrinsicreactivity.Itwasshownthat bothporos- ityandporestructurecanaffectcharconversion.

As demonstrated by the above-mentioned examples, a great deal of understanding can be reached with the steady state as- sumption. However, all transientphenomenaandprocesses(igni- tion,volatilesburnout,progressofcharconversion,combustionin non-laminar flow) requireunsteady approach. The first transient simulations of resolved particle combustion in a non-quiescent (two-dimensional) flow were performed by Lee et al. [22] using the spectral element method. Recently, Farazi etal. [23] used an unsteadyapproachandadetailedchemicalmechanism,andinves- tigated charparticle combustioninairandoxy-fuel atmospheres.

Thecombustioncharacteristicsinthesetwoatmosphereswereex- plored, as well as interactions between kinetics andmass trans- fer. This work was further extended to particle arraysby Sayadi et al. [24]. Anotherstudy on the resolved particle conversion in whichthegoverningequationsweresolvedintheirunsteadyform was done by Luo et al.[6]. In their work, an immersed bound- arymethodandasimplesemi-globalmechanismwereutilized.Fi- nally,Tufanoetal.[25]performedthemostcompletestudyupto date,inwhichallstagesofthecoalparticleconversionareconsid- ered, i.e. heating, drying, ignition, volatiles combustion andchar particle conversion. Moreover, in addition to detailed chemistry, their numerical model accounts for complex features of particle interior, such as time evolution of porosity and tortuosity. Most recently, Nguyen etal. [26] performedunsteady particle-resolved simulations toinvestigatethe evolutionofchar particlemorphol- ogy.Based ontheirresults,improvedexpressions forthemode of burningandtheRandomPoreModelwereproposed.

In the existing literature on resolved particle conversion,very differentlevelsofnumericalmodelcomplexityarepresented.The currenttrendseemsto betowardsmoreandmoredetailedmod- elsandmodelsthatareabletocapturetransienteffects.However, highaccuracy isachievedattheexpenseofefficiency. Theobjec- tive ofthisworkistoproposeanovelnumericalapproachforre- solvedcharparticlecombustionmodeling.Contrarytothepresent trendintheliterature,we aimforthemodeltobeassimpleand efficientaspossible,whilestillpreservinghighaccuracyandbeing able topredictunsteady phenomena.Thisisaccomplishedby us- ing structured, overset gridsand byintroducing carefullyverified assumptionsandsimplifications.

2. Governingequationsandnumericalmethods

An open-source, compressible solver called the Pencil Code [27] is used to perform the simulations presented in this work.

The PencilCode uses a 6thorder finite difference scheme anda 3rdorderRunge–Kuttaschemeforspatialandtemporaldiscretiza- tion, respectively. One ofthe main features of the numericalap- proach employedin thisstudyisthe oversetgrid.The particleis surroundedby acylindricalbody-fittedgrid(lateralsoreferredto as‘ogrid’),whichspansthespacebetweenr=r_ptor=3r_p=r_ogrid, where r is a radial coordinate andrp is the particle radius. The restofthecomputationaldomainisresolvedontheCartesiangrid.

Such an approachallowsone tousevery highresolutioncloseto the particle,whichisnecessary toresolveits boundarylayer and the surrounding flame. Further away from the particle, the grid

is much coarser, making the computational effort relatively low.

The solutionis interpolatedbetween theogrid andthe Cartesian grid using a 4thorder, explicit Lagrangian interpolation method, whichhasbeenshowntobeanoptimalchoiceinconnectionwith a6thorderfinitedifferencescheme[28,29].Inordertoavoidspu- riousoscillations,Padé filtering[30,31]isappliedonthecylindrical gridto density,temperatureandvelocityfields.The detailsabout the implementation of the overset grid and performance of this methodcanbefoundin[32,33].

2.1. Fluidequations

The continuity and momentum equations are solved in their non-conservative,compressibleform:

∂ρ ∂

^{t +}

∇

·

( ρ

^u

)

=0, (1)

ρ ∂ ∂

^{ut +}

ρ

^u·

∇

^u=−

∇

^p+

∇

^·

τ

^{+f, (2)}

where

ρ

^{andparethedensityandpressure,}respectively,andthe bold symbols representthe velocity (u) andvolumetric force (f) vectors.Thestresstensor,

τ

^,isgivenby

τ

=

μ ( ∇

^u+

( ∇

^u

)

^T

)

−2

3

μ ( ∇

·u

) τ

, (3)

where

μ

^{stands for the dynamic viscosity and}

τ

^{is the identity}

matrix.Themassfractionofchemicalspeciesk,givenbyY_k,obeys thefollowingtransportequation

ρ ∂ ∂

^Ytk+

ρ

^u·

∇

^Yk=−

∇

^·Jk+

ω

˙k, (4) inwhichthediffusiveflux, Jk,issimplifiedby usingtheassump- tionofFickiandiffusion,suchthat

Jk=−

ρ

^Dk

∇

^Yk, (5)

whereD_kisthediffusioncoefficientofspecieskand

ω

˙krepresents thegasphasereactionrateofthesamespecies.

Byneglectingviscousheating,theenergyequationisexpressed intermsoftemperatureas[34]

ρ ∂ ∂

^{Tt +}

ρ

^u·

∇

^T

=

k

( ω

^˙k−

∇

·Jk

)

TR

c_vMk−hk

c_v

−

ρ

^TR

c_vM

∇

·u−

∇

·q

c_v , (6) whereTrepresentsthetemperature,c_vistheheatcapacityatcon- stantvolume,R isthe universalgas constantandM isthemolar mass for the mixture, 1/M=

kY_k/M_k. The heat flux, , is com- putedas

=

k

h_kJk−

λ∇

^{T, (7)}

where

λ

^{representsthermal} conductivity andh_k=^hs,k+h⁰_f,k is theabsoluteenthalpyofspeciesk,whichisthesumofitssensible enthalpy, ^hs,k, and its heat of formation, h⁰_f,k. Finally, to relate densitywithpressure,theidealgasequationofstateisused, p=

ρ

^RT

M . (8)

2.2. Chemicalmechanismandboundaryconditions

Asimplified chemical mechanism that consistsoftwo surface reactionsandonereversiblegasphasereactionisemployed:

2C+O2→2CO (R1)

Table 1

Kinetic parameters. Here, [ a ] denotes concentration of species a , k iis given by Eq. (9) and r irep- resents the rate-of-progress variable. Note that for surface reactions units of r iare mol / cm ²/ s , while for gas phase reactions it is mol / cm ³/ s .

reaction B i E i[ kcal / mol ] r i source

R1 1 . 97 ×10 ⁹cm / s 47.3 k 1[ O 2] [36]

R2 1 . 29 ×10 ⁷cm /s 45.6 k 2[ CO 2] [36]

R3 (forward) 3 . 98 ×10 ¹⁴(^cmmol³)^3/4/ s 40.7 k 3,f[ CO ][ H 2O ] ^1/2[ O 2] ^1/4 [37]

R3 (reverse) 5 ×10 ⁸1 / s 40.7 k 3,r[ CO 2] [37]

C+CO₂→2CO (R2)

CO+0.5O2↔CO2 (R3)

It should be noted that the gasification reaction through H₂O is not considered inthe presentstudy, even though watervapor is present inthe atmosphere. Thereason thisreactionwas omitted wastheverylowconcentrationofH₂O(Y_H

2O=8×10⁻⁴atthein- let),whichhasbeenshownbyKesteletal.[15]tohaveessentially no effecton the conversionrate. To studycasescharacterized by highercontent ofwatervapor, theadditionalgasificationreaction andthewater-gasshiftreactionshouldbeincludedinthemecha- nism.TheArrheniusexpressionforreactionireads

k_i=B_iexp

(

^−Ei/RT

)

^{. (9)}

Theempiricalkineticparameters:pre-exponentialfactorB_i,activa- tion energyE_iandreactionorders arelistedinTable1.The reac- tiontermforthegasphasereactioninEq.(4)iscomputedas

ω

˙k=M_k

nr,gas

i=1

( ν

ki−

ν

ki

)

^ri, (10)

where

ν

_ki ^and

ν

_ki ^{are the}stoichiometric coefficientsofgasphase species k in reaction i on the reactant and product side, respec- tively,while n_r,gas isthe numberofgasphase reactions,andr_i is the rate-of-progressvariable (adoptingterminologyfromCh. 4in [35]),asgiveninTable1.

Sincetheparticleinteriorisnotincludedinthecurrentframe- work,itisassumedthatallcontributions tothereactionratedue to internal reactions are accounted for through the apparent ki- neticparameters,andthatthetemperaturegradientinsidethepar- ticleissmallenoughtobeneglected.Also,theparticleisassumed to be entirely madeof carbon andthe modeldoes not incorpo- rateparticle shrinkageduringits conversion.In reality,theparti- clesizeanddensityareslowlychangingascombustionprogresses [38].However,thetypicaltimeofoursimulationsismuchshorter than the burnout time ofthe particle such that the reduction of theparticlediametercanbeconsiderednegligible.

As statedabove, theinterior oftheparticle isnot included in the computational mesh. The interaction between the solid and the surrounding gas is therefore incorporated through the parti- cleboundaryconditions.Wewillnowcontinuebydescribingthese boundary conditions. The species balance at the cylinder surface canbeexpressedas[6]:

ρ

^Dk

∂

^Yk

∂

^{r +m˙cYk+m˙k=0, (11)}

where

˙ m_k=M_k

n_r,heter

i=1

( ν

ki−

ν

ki

)

^ri, (12)

istheproductionrateofspecieskduetoheterogeneousreactions, and n_r,heter is the number of heterogeneous reactions. The char conversionrateisgivenby

m˙c=−MC

(

^2k1[O2]+k2[CO2]

)

=−

(

^m˙O2+m˙CO2+m˙CO

)

=−

ns,gas

k=1

˙

m_k, (13)

wherethefinalsummationisoverallgas-phasespecies.Adetailed deductionof Eq. (11) can be found inAppendix A. It should be notedthat both m˙c andm˙_k depend onthe speciesconcentration on the surface, whichmakes it necessaryto solve Eq. (11) in an iterativemanner. Anotherpossibility isto use speciesproduction rates fromthe previous time step, this can howeverlead to nu- merical instabilities and non-physical results.Here, we employ a simpleiterativealgorithmtosimultaneouslyfindsolutionsforY_O2 andY_CO

2 atthesurface,whiletheremainingspeciesaresolvedfor directly.

Mass conservation at the particle surface requires that (see AppendixA)

k

( ρ

^Yku+Jk

)

·rˆ=

k

m˙k=−m˙c, (14)

where rˆ is the vector normal to the particle surface. From the aboveequation,andsince

kJk·rˆ=0,theboundaryconditionfor velocitybecomes:

ur=−m˙c/

ρ

, (15)

where ur is the outward velocity in the radial direction, corre- spondingtothesocalledStefanflow.

Dirichletboundaryconditionisemployed forthetemperature.

The intention behind the Dirichlet boundary condition for tem- perature isto validate the code against theexperimental data of Makino et al.[39], where the temperature was maintained con- stant. The last variable that needs to be defined at the cylinder surfaceisdensity,which issolved fordirectlyfromthe transport equation anddoestherefore not requireanyspecialtreatment at theboundary.

2.3. Transportproperties

Insimulations ofreactingflows,itiscommonpracticetocom- pute transport coefficients, such as

μ

k, D_k and thermal diffusiv- ityD_th,based onthe kinetic theory ofgases, asdescribede.g. in [34].Thisapproach,whileaccurate,significantlyincreasescompu- tationalcost.Thisisespeciallythecaseforspeciesdiffusioncoeffi- cientsforwhichbinarydiffusioncoefficientsneedtobeevaluated first. In order to maximize computational efficiency, a simplified approachisemployedinthiswork.Atthesametime,careistaken nottocompromisetheaccuracyoftheresults.

Thekinetic viscosityisrelatedtotemperaturethroughSuther- land’slaw

ν

⁼

ρ (

^CT1T+3/C22

)

⁽¹⁶⁾

with constants C1=1.52·10⁻⁶ kg/m/s/K^1/2 and C2=110 K. The above expression is fully applicable to single-component gases.

However, if a mixture is dominated by components with similar

Table 2

Polynomial coefficients for heat capacity in the temperature range 10 0 0 K < T < 50 0 0 K .

species CO CO 2 H 2O N 2 O 2

a 1 3.025 4.454 2.672 2.927 3.698

a 2 1 . 443 ·10 ⁻³ 3 . 140 ·10 ⁻³ 3 . 056 ·10 ⁻³ 1 . 488 ·10 ⁻³ 6 . 135 ·10 ⁻⁴ a 3 −5 . 631 ·10 ⁻⁷ −1 . 278 ·10 ⁻⁶ −8 . 730 ·10 ⁻⁷ −5 . 685 ·10 ⁻⁷ −1 . 259 ·10 ⁻⁷ a 4 1 . 019 ·10 ⁻¹⁰ 2 . 394 ·10 ⁻¹⁰ 1 . 201 ·10 ⁻¹⁰ 1 . 010 ·10 ⁻¹⁰ 1 . 775 ·10 ⁻¹¹ a 5 −6 . 911 ·10 ⁻¹⁵ −1 . 669 ·10 ⁻¹⁴ −6 . 392 ·10 ⁻¹⁵ −6 . 753 ·10 ⁻¹⁵ −1 . 136 ·10 ⁻¹⁵

Fig. 1. Kinetic viscosity as obtained using Sutherland’s law ( Eq. (16) ) and multi- component approach for the mixture consisting of Y N2= 0 . 7292 , Y O2= 0 . 05 , Y H2O= 0 . 0 0 08 , Y CO= 0 . 02 and Y CO2= 0 . 2 .

properties(asisthecasehere),Eq.(16)isreducedtoadecentap- proximation.Furthermore,constantsC₁ andC₂ wereselectedsuch thatforawiderangeoftemperaturesandcompositionsthekinetic viscosity resultingfromEq.(16)is in agood agreement withthe kinetic viscositydeterminedusingthe multi-componentapproach (i.e.basedonkinetictheory).InFig.1,thesetwomethodsarecom- paredfora typicalcomposition encountered inthecurrentwork.

For other compositions that are likely to occur, a deviation from the kinetic theory remains below 7% for the temperature range presentedinFig.1.

The main assumption allowing us to compute the remaining transport coefficientsisthatthetransport coefficientsare propor- tionaltoeachother,i.e.

ν

=PrDth=PrLekDk, (17)

with the constants of proportionality beingthe Prandtl (Pr) and Lewis (Le_k) numbers. Such an assumption of constant Prandtl and/or Lewis numbers has successfully been applied in re- cent studies on resolved particle devolatilizationandcombustion [12,23].Typically,Pr=0.7andLe_k=1forallspeciesareassumed.

This was shownto havea negligible impact on the devolatiliza- tionstagewhencomparedwiththecomplexmulti-componentap- proach [12]. However, in some conditions, the combustion rate mightbeaffectedbydiffusioncoefficients,aswillbedemonstrated in the next section. Therefore, a more careful approach is em- ployed,asdescribedbelow.

Theheatcapacityatconstantpressureisgivenby

cp=

k

Y_kc_p,k= R M

k

Y_k 5

i=1

a_iTⁱ⁻¹, (18)

where the polynomial coefficients a_i are takenfrom Gordon and Mcbride[40]andarelistedinTable2fortherelevanttemperature range.Theheatcapacityatconstantvolumeisrelatedtotheheat capacityatconstant pressure throughthegas constant,such that

Fig. 2. Thermal conductivity as obtained using Eq. (17) with Pr = 0 . 9 and multi- component approach for the mixture consisting of Y N2= 0 . 7292 , Y O2= 0 . 05 , Y H2O= 0 . 0 0 08 , Y CO= 0 . 02 and Y CO2= 0 . 2 .

Table 3

Selected Lewis numbers.

species CO CO 2 H 2O N 2 O 2

Le k 0.78 1.01 0.58 0.7 0.78

cp−c_v=R/M. (19)

Usingtheheat capacitygivenby Eq.(18)andthethermaldif- fusivitygivenbyEq.(17),thethermalconductivity,definedas

λ

=cp

ρ

^Dth, (20)

isshowninFig.2asafunctionoftemperatureforthesamemix- ture asused in Fig. 1.In Fig. 2, the thermalconductivity asob- tainedusingthemulti-componentapproachisalsopresented.The bestagreement betweenthesetwo functionsforawide rangeof mixturesisachievedbysettingthePrandtlnumberequalto0.9.

Foreachspecies,Le_kischosensuchthattheresultingdiffusion coefficientdoesnotdifferbymorethanaround10% fromthedif- fusioncoefficient computedbasedon themulti-componentdiffu- sionapproach.Thiswasverifiedforthefullrangeofcompositions andtemperaturesthat arelikelytoappearinthecasesweexam- ine.Figure 3presentsa comparisonbetweenthe diffusioncoeffi- cientsasafunctionoftemperatureascomputedfromEq.(17)and asobtainedusingthemulti-componentdiffusion.Themagnitudes oftheLewisnumbersleadingtotheseresultsarelistedinTable3. Agoodagreementbetweenthetwoapproachesisachievedforall transportcoefficients(

ν

^,

λ

^andDk), whichjustifiestheuseofthe simplifiedapproachforthetransportcoefficients.

In order to quantify the efficiencygain obtained by simplify- ingtheformulationofthetransportcoefficients,aone-dimensional flame was simulated for two cases (details regarding the one- dimensional flame simulations are given in the next section).

In the first case, transport properties were computed according

Fig. 3. Diffusion coefficients as obtained using Eq. (17) (referred to as ‘simplified’) and multi-component approach for the mixture consisting of Y N2= 0 . 7292 , Y O2= 0 . 05 , Y H2O = 0 . 0 0 08 , Y CO = 0 . 02 and Y CO2= 0 . 2 .

to Eqs. (16), (17) and (20), while in the second case, a multi- component approach wasemployed. A comparisonof the execu- tiontimeofsubroutinesresponsibleforcomputingtransportprop- ertiesrevealedthat7.5timeslesscomputationaltimewasrequired for thecasein whichthe simplifiedapproach wasused.Further- more, since these subroutines are computationally the most ex- pensive(i.e.theirexecutiontakesalargefractionofthesimulation time),thiscorresponds toa reductioninthe totalexecutiontime byafactorof3.4.Itshouldalsobenotedthattheefficiencygainis dependentonthenumberofspeciespresentinthesimulation.The reasonforthisisthat one additionalnestedloop overall species must beexecuted andasignificantly largernumberofoperations have tobe performedtocompute transport coefficientsbasedon thekinetictheory.Inourcase,thefactorof7.5wasachievedfor5 species.

2.4. Speedofsoundreduction

Numerical stability of the simulations requires several condi- tions to befulfilled.First ofall, arequirementdueto convection, oftencalledtheCFLconditionlimitsthemaximumtimestepto:

^t≤ C

^x

max

(

^cs+u

)

^{, (21)}

whereCisaconstantthat dependsona numericalscheme(typi- callyC≈1)and

cs=

γ

^{RT/M (22)}

isthespeedofsoundand

γ

=cp/c_v.Forreactingflows,thelength of the time step and the grid spacing is most often limited by chemicalscales.However,itturnsoutthatinthecaseofflowsthat arebothreactingandcompressible,theresolutionrequirementdue to the ratio between viscosityand the speed of soundmight be more restrictive. Fortheparticular numericalapproach employed inthePencilCode,ithasbeenshown[41]thatthegridspacingis constrainedby

^x<

βν

cs , (23)

where

β

∼50. It follows fromEq. (23)that larger grid spacing, and hence lessmesh points, maybe used ifthe speed of sound is reduced. A goodrule of thumb is that, aslong as we are not interestedinthermo-acoustics,theresultsareindependentofthe Mach number,Ma=u/cs,forall Machnumbersbelow 0.1.In our

Table 4

Initial conditions for 1D carbon monoxide flame.

reactant side product side

Y O2 0.165 0.0

Y CO 0.29 0.0

Y CO2 0.0 0.455

Y H2O 0.0008 0.0008 Y N2 0.544 0.544

T [K] 298 2000

case,theMachnumberistypicallyoftheorderof10⁻³.Thespeed ofsound can therefore be reducedby up to two orders ofmag- nitudewhilestillmaintainingMach-independentresults.Sincethe timestepisoftenlimitedbytheCFLcondition, whichistypically thecaseforlowertemperatures,areductionofthespeedofsound wouldalsoallowustouselargertimesteps.

Inthe previous paragraph we showedthat areduction in the speedofsoundcouldbe verybeneficialfortheCPU consumption ofoursimulations,andthattheeffectsuchareductionhasonthe resultsshould be negligible if theMach number iskept below a certainvalue.Thequestionnowishowthespeedofsoundcanbe changedwithout affectinganyother aspect ofthe results.Thisis donebydividingthegasconstantbyafactor

α

^2,suchthat

R→R/

α

^{2, (24)}

whichimpliesthat(22)

cs→cs/

α

^{. (25)}

Thegasconstantischangedconsistentlyforallequations,withthe exceptionofEq. (9)inwhichtheoriginalmagnitudeofRmustbe usedinorderforthereactionratenottobeaffected.Itshouldbe notedthatthereductionofRmeansthatcp,c_v and

λ

^arealsore-

ducedbythesamefactorof

α

^{2,ascanbe seenfromEqs.(19)to}

(20). However, this has no effect on the energy equation as all thesereductionscanceloutineverytermofEq. (6).Theonlyterm that is affected is the pressure gradient term in the momentum equation,since

∇

^p∼c²_s,whichisasintended.

We will now validate the assumption that a reduction in the speed of sound does not affect the main results, except for the acoustic waves, as long as the Mach number is below 0.1. This is done by simulating reacting flows of a one-dimensional car- bonmonoxideflamewiththreedifferentvaluesofcs.Inthebase case, the speed ofsoundwas kept unchanged, whichresulted in Ma≈0.001, in the other cases the speed of sound wasreduced by factors of 10 and 50,which led to Ma≈0.01 and Ma≈0.05, respectively. The initial conditions for these cases are given in Table4, whilethe one-step mechanismgivenin Section 2.2gov- ernstheflame.

Theresultingtemperatureandspeciesmassfractionprofilesat steadystatearepresentedinFig.4,fromwhichitcanbeseenthat theresultsarenot affectedbythespeed ofsoundreduction.Fur- thermore,forallthreecases,thesameflamespeed,S_L=14cm/s,is obtained.Havingverifiedthat thespeedofsoundcanbereduced withoutaffecting the results,thistactic isemployed forall cases discussedinthenextsection,whichresultedinamajorreduction ofCPUpowerconsumption,inparticularforthosecaseswherethe time-step wasnot limitedby chemical reactions. Itis alsoworth mentioning that the efficiency gain resulting from the speed of soundreductionisverycase-dependent.Thiscanbeillustratedby subsequentlyreducingthespatialresolutionofthe1D flamesim- ulation inwhich the speed of soundwas reducedby a factorof 10(correspondingtothegreenlineinFig.4).Despitethefactthat themaximum gridsize, asdefinedbyEq. (23),isinverselypro- portionaltothespeedofsound,itwaspossibletoreducetheres- olutiononlybyafactorof∼3duetothefactthat,forstabilityrea-

Fig. 4. Comparison of temperature and species profiles across the flame obtained before and after the speed of sound reduction. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 5. Schematic representation of the analyzed case (not drawn to scale).

sons,acertainnumberofgridpointsarerequiredacrosstheflame front.Itcanbethereforeconcluded,thatthespeedofsoundreduc- tion allows one to eliminate the grid size/time step requirement due to the speed of sound in low Mach number flows, butthe efficiency gain associated withthis cannot be quantified in gen- eralbasis sinceitdependsonother case-specifictimeandlength scales.

2.5. Numericalset-up

Theset-upforallsimulatedcasescorrespondtotheexperimen- tal set-up of Makino et al. [39] and can be summarized as fol- lows. Acylindrical particle of5mm in diameteris placed inthe middleofa10cm×8cmcomputationaldomain.Thefluid,which has a composition that is typical for air (Y_N

2=0.77, Y_O

2=0.23, Y_H2O=0.0008)entersthedomainthroughonesidewithavelocity

of1m/sinthey−direction.Periodicboundaryconditionsarespec- ifiedinthetwocross-flowdirections.Initially,thetemperaturein- sidethedomainiseverywhereequalto1280K.Theinitialspecies distribution on the ogrid is such that the oxygen mass fraction decreases exponentiallyfromY_O2=0.23atr=r_ogrid toY_O2=0at theparticlesurface(r=rp),whilecarbondioxideisintroducedin placeofoxygen, i.e.Y_CO

2(^r)=Y_O

2(^rogrid)−Y_O

2(^r)^{.The initial com-} positionon the Cartesian grid isthe sameas thecomposition at the inlet. Such initial conditions do not reflect the experimental set-upandwereselectedpurelytoimprovestabilityofsimulations duringtheinitialstage.

For most cases, a grid resolution of 720 x 896 (x x y direc- tions) grid points on the Cartesian grid and 208 x 432 (r x

θ

directions) on the ogrid was sufficient to accurately resolve all flow features. It should be notedthat the ogridis stretched in a non-linearmannerintheradialdirection. Fortheresolutiongiven

abovethisresultedin^rmin=8.3·10⁻⁴cmattheparticlesurface and^rmax=6.8·10⁻³cmattheouteredgeofthecylindricalgrid.

Aschematicrepresentationofthenumericalgridtogetherwithini- tial condition is presented in Fig. 5. If the particle temperature is relatively low (Tp1800K) the maximumtime-step is limited to ∼10⁻⁷ sby convection,whileforhigherparticletemperatures thetime-step needstobereducedto∼10⁻⁸ sduetotheshorter chemicaltimescales.

3. Resultsanddiscussion

3.1. Implementationofchemistrymodule-validation

Various aspects of the Pencil Code have been validated and tested a number of times and the results have been published in a large numberof papers available in theopen literature. See [27] for an overview of some relevant papers. In this work, we have, however, implemented several new methods and approxi- mations tospeed upthe calculations, suchas: simplifiedcalcula- tionoftransportdata,simplifiedglobalreactionmechanisms,het- erogeneous reactions attheparticlesurfacewiththeoversetgrid method,andvariablespeedofsound.Inordertovalidate thecur- rent numerical model beyond the more specific validations pre- sentedintheprevioussection,theexperimentalset-upofMakino etal.[39]isreproducednumerically.IntheexperimentofMakino et al., combustionof agraphite rodwasstudied atdifferent sur- facetemperatures,fordifferentairvelocitiesandtemperatures.An importantfeatureoftheexperimentisthattheheatlossfromthe graphitesurfaceduetoradiationisbalanced byelectricalheating, such thataconstantparticlesurfacetemperatureismaintainedat all times. As a result, a quasi-steadystate is achievedfor a rela- tively largefractionoftheparticleconversiontime.Inthecurrent work,thecasecharacterizedbyan airtemperatureof1280Kand avelocity of102.5cm/sisanalyzedforarangeofparticlesurface temperatures. Thisparticularselection ofexperimentalconditions wasmotivatedbythefactthatthesamecasewasstudiednumer- ically byLuoetal.[6],whodemonstratedthat agoodagreement withtheexperimentalresultscanbe obtainedusingthechemical mechanismgivenbyreactions(R1)–(R3).Despitethefactthat Luo etal.alsousedthePencilCode,therearetwomaindifferencesbe- tween theirapproachandtheapproachusedinthecurrentwork:

(1)Luoetal.usedkinetictheorytocomputetransportcoefficients, and (2)their particlewas resolved ona Cartesian gridusing im- mersedboundaryconditionsfortheparticlesurface.

Figure6presentsthecarbonconversionrateobtainedwiththe current numerical approach (green squares) in addition to what was foundexperimentally by Makino etal.[39] (red circles) and numericallyby Luoetal.[6](bluecircles).Infact, whatisshown is the conversion rate inthe forward stagnation point. Addition- ally,kinetic(solidblueline)anddiffusion(dottedblackline)limits foroxidationarealsoincludedinFig.6.Thefirstlimitcorresponds tothecaseofinfinitelyfastdiffusion(Y_O

2,sur f ace=Y_O

2,∞), whilethe latter to the reaction ratebeing controlled by diffusion (∼T^1/2).

It canbeseenthat uptoTp=1200K,thecarbonconversionrate is governed by kinetics,while around T_p=1600 Kthe slope cor- responding to the diffusion limit is achieved. There is one more limiting slopeincluded inFig. 6,which iscalled‘flame diffusion’

limit. This limit arises due to the fact that at around T_p=1700 K the flame begins to detach fromthe particle surface. The rea- son forthisdetachmentisthelargeCOproductionatthesurface anditssubsequenttransportbymeansoftheStefanflowanddif- fusion. The resultis that mostofthe O₂ is consumedin the gas phase atthe position ofthe flamethat is formedaway fromthe surface. As a consequence, mostlyCO2 can diffuseto the surface andthecarbonconversionisduetotheBoudouardreaction(R2). Fromtheperspectiveoftheoxidationreaction,theoxygendiffuses

Fig. 6. Comparison of carbon conversion rates as a function of particle surface tem- perature. The results for Luo et al. [6] are reproduced from their Fig. 8 . (For inter- pretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 7. Effective flame radius in a function of temperature.

Fig. 8. Comparison of temperature profiles along the centerline with and without gas-to-gas radiation model.

now towards the flame surface, not the particle surface. This ef- fective surfacegrows proportionally to T^g,where theexponent g canbefoundbyafittingprocedure.Thiswasdone inFig.7,from whichitcanbeseenthatthe‘effectiveradius’scalesasT^0.78.Here, the effectiveradius wascomputed asthe average radial distance fromtheparticlecentertothe flame,whereitwasassumedthat theflamelocationcorresponds tothegridpointinwhichthegas phase reaction rateis thehighest. Thecarbon conversion ratein

Fig. 9. Contributions to CO production from gasification and oxidation.

thediffusionlimitisproportionaltotheproductofthemasstrans- fercoefficient(k_i)andtheeffectivesurface:

˙

mc∼d²_{p,e f f}k_i, (26)

where d_{p,e f f} isthe effectivediameter ofthe flamesurface. Since the masstransfercoefficient scalesask_i∼D_i/d_{p,e f f} andD_i∼T^1/2 (see Eqs. (16)and(17)),the conversionratedependenceon tem- peraturebecomes:

˙

mc∼d²_{p,e f f}k_i∼d_{p,e f f}D_i∼T^1/2T^0.78=T^1.28. (27) ThisistheflamediffusionlimitseeninFig.6,whichisreachedfor thehighestofthestudiedparticlesurfacetemperatures.

Compared to the experimental results, slightly too low con- version is obtained for most temperatures. On the other hand, verysimilarmagnitudesofconversionrateswereobtainedbyLuo et al. [6], which indicates that the difference is most probably caused by the reaction kinetics. It is in fact common that there isnoagreementonthereactionkineticsandquiteoftenanumber ofmechanismsare suggested,resultingindifferentreactionrates.

The influenceof heterogeneouskineticshas alreadybeeninvesti- gated by Nikrityuk etal.[42],who revealedthat a factorof2 or even3differenceinthecarbonconsumptionratecanbeexpected betweendifferentsetsofkinetic parametersthatarefound inthe literature.Asetofkineticparametersforsurfacereactionswasalso proposedbyMakinoetal.[39]basedontheirexperimentalresults andtheconversionratesresultingfromtheseparametersarepre- sentedinFig.6(cyantriangles).Itcanbeseenthatthisyieldeda significantlyhighercarbonconversionrateathighsurfacetemper- atures, but didnot lead to noticeable difference for T_p≤1600 K.

Thiscould be expectedasthe gasificationreactionismuch faster in Makino’s mechanism, while there is only a tiny difference in the oxidationrateswhencomparedwiththemechanismgivenin Table1.

Anotherexperimentalfeaturethatisnotcapturedproperlywith the currentapproachisa sudden decreaseofthe conversionrate forsurfacetemperaturesaround1700–1800K.Thisdecreaseisalso presentintheresultsshowninFig.6inLuoetal.[6](althoughthe resultsintheirFigs.6and8seemtobeinconsistentregardingthis feature). The main difference between their and the presentnu- merical approach ishow thetransport coefficientsare computed.

Inthatrespect,ourapproachismuchsimplerand,potentially,less accurate. Therefore, a further validation is essential. Such a vali- dation wasperformedusingtheANSYSFluentsoftware,inwhich thesamecaseswerereproducedandtheresultingcarbonconver- sionratesareshownasblackx-signsinFig.6[addcontourshere oraplotshowingTcomparisonalongcenterline].TheFluentsim- ulations wereperformedwiththediffusioncoefficientscalculated

fromkinetic theory,as wasalso done by Luo etal. [6]. In addi- tion,incompressibleandsteadystate flowwasassumed. Bothas- sumptionsarevalidsincetheMachnumberislowandthechange in particle radius is very slow. As can be seen in Fig. 6, almost the sameconversion rates were obtained usingthe complex for- mulationforthetransport coefficientsinANSYSFluentasforthe simplified formulation usedin the Pencil Code. In particular, the conversionrates inbothcases aremonotonically increasing func- tions,i.e.noreductionoftheconversionratewasobservedaround T_p= 1700–1800 K. This verifies that the simplified approach for thetransportisnotresponsibleforthisqualitativediscrepancybe- tweentheexperimentalresultsandournumericalresults,andal- lowsustogainconfidenceinthepredictionsofourapproach.

Thecasethatwassetup inANSYSFluentwasalsousedtoes- timatetheinfluenceofgasphaseradiation,whichwasomittedin the energy equation in the Pencil Code. While it is not uncom- mon to omitgas-to-gas radiationin simulations ofconversion of resolvedcharparticles,somestudiessuggestthatitseffectisnon- negligible.Forexample,asignificantreductionofthecharparticle surfacetemperature dueto gas phase radiationwas observed by Richteretal.[16],especiallyforcaseswithhighambienttempera- ture.Ontheother hand,Tufanoetal.[12]showedthattheeffect ofgas-to-gasradiationonignitionisratherweak.Inourstudy,the gasphase radiationwasaccountedforthrough theDiscreteOrdi- nates model,andits influencecan be seeninFig. 8,which com- paresthetemperaturedistributionalong thecenterline ofthere- actingparticleforthecaseswithandwithoutradiation.Thecase withTp=2000Kisshownheresincetheeffectofradiationisthe highest for cases with highparticle temperature. It can be seen thattheeffectonthegasphasetemperaturefieldiscertainlynon- negligibleintheregionbehindtheparticle.Nevertheless,thecon- versionrateremainedunchangedduetotheexperiment-imitating assumption of constant temperature at the particle surface and virtually no influence of radiation on species concentrations. It shouldbe noted,however,thatbasedon theresultspresentedby Luoetal.[6],itisexpectedthat theparticlesurfacetemperature is unlikely to change by more than a few percent for the cases studiedin thepresentpaper, evenifheat transferat theparticle surface(chemical heatrelease, conduction,convection andradia- tion) wasaccountedforthrough theparticle boundarycondition.

Fig. 10. Upper: conversion rates for different diffusion coefficients, lower: contributions from gasification/oxidation to the CO production rate, T p= 1700 K.

Itstill remainstobe understoodwhyconversionrateobtained withthe PencilCode (andANSYSFluent) doesnot followthe ex- perimental trendwhen it comes tothe dipin carbonconversion ratearoundTp=1700K.Thereexistseveralphysicalexplanations of thistrend in the literature, e.g.:it is attributedto the change ofthe effectivereactionzone thickness[43],it islinked withthe change ofmolecular structureof graphite[44,45],itis causedby thermal rearrangement of surface-covering sites, fromhighly re- active atlow temperaturestolessreactiveathighertemperatures [46,47].Makinoetal.[39]arguethatthepresenceofthedipstems fromthefactthatthedominantsurfacereactionshiftsfromoxida- tion to gasification around T_p=1700 K.The reason forthis shift isthatatlowtemperaturestheoxygenisusedtooxidizethecar- bon directly at the surface, while athigh temperatures the oxy- genis usedtooxidizeCOina CO-flamesurroundingthe particle, while thecarbonconversionproceedsthrough gasificationofCO₂ that diffusetothe surfacefromtheCOflame.Thischangeinthe dominantmechanismforCOproductionatthesurfaceiscorrectly predicted by the Pencil Code, ascan be seen in Fig. 9.However, theshiftisgradualanddoesnotresultinthenon-monotonicityof m˙_C(^Tp) ^{assuggested by Makinoetal.[39].Another plausibleex-} planation forthedipinm˙_c is thatsincethe shapeoftheconver- sionfunctiondependsonthegasphasekinetics,asshownin[48], the kinetic parameters we use might not yield the right behav- ior. While all theabove explanations are probable, itis alsopos- siblethattheresultsareaffectedbythemeasurementmethod.In theexperiment,thesurfacetemperatureoftherodwasmeasured using two-colorpyrometer[49].Thesemeasurementsare usedto controltheinternalheatingthatisrequiredtomaintainaconstant temperatureofthegraphiterod. Thismethodisindirect,it might therefore be difficult to precisely measure the surface tempera- turewithout theresultsbeingaffectedby thesurrounding flame.

Atrelativelylowsurfacetemperatures,theflameremainsattached to the surface, sothe difference between theflame and thesur- facetemperatureissmall.However,attemperaturesatwhichthe dropin the conversion rateis observed,the flame startsdetach- ing from the rod surface. As such, the flame temperature might besignificantly higher,givingafalseimpression ofhighersurface temperature.Sincetheexperimentattemptstomaintainaconstant surfacetemperature, itis likelythat therodwascooled to lower temperaturethanintended,whichresultedinasuddendecreaseof theconversionrate.Theseare,however,onlyconjectures,andthe reason for the qualitative inconsistency between the experiment andour resultsmight be a combinationof severalof the above- mentionedfactors.

3.2. Sensitivityanalysis

Inordertobetterunderstandwhichparametersthatcontrolthe carbonconversionrate,we havedone aseriesofparameterstud- ies.The firststudyinvestigatesthe effectofspeciesdiffusivity.In thisrespectwevariedthediffusivitiesofO₂,CO₂andCOfromhalf oftheir originalvalue upto twicethe originalvalue, andinvesti- gated how thisinfluenced thesolid (carbon) conversionrate. For thisinvestigation,we concentrateonthesituationwherethepar- ticletemperatureis1700K.

From theupperpanel ofFig.10 we seethat thesolid conver- sionratehasastrongdependenceondiffusivityofO₂.Thisisex- pectedsince higherdiffusivityofO₂ willyield ahigher transport rateofO₂tothesolid,whichwillthenbeabletoconvert(oxidize) moresolid.FromthelowerpanelofFig.10weseethatincreasing theoxygendiffusivityresultsinan increaseofboth theoxidation andgasificationratesofthesolid.Atfirstglance,itmaylooksur- prisingthateventhegasificationrateincreaseswithincreasedO₂

Fig. 11. Oxygen and carbon dioxide mass fractions and temperature profiles in the particle stagnation region as obtained for different diffusion coefficients.

diffusivity,butthereasonissimplythatatthesurfacetemperature of1700Kthatwefocusonhere, weexperienceahighersurface- fraction of CO₂, resulting fromoxidation of COvery close tothe surface.

LetusnowmoveontotheeffectofCOdiffusivity.Weseefrom theupperpanelofFig.10thatthesolidconversionrateisweakly increasing with increasing diffusivity of CO. This effect is, how- ever,morecomplicatedthanthatofO₂ diffusivity,ascanbeseen fromthelowerpanelofFig.10,whichshowsthatsolidconversion dueto oxidationincreaseswithCOdiffusivity, whiletheopposite is true forgasification. Toelucidate this behaviour inFig. 11,we show CO2 and O2 concentrationsalong they-axis infront ofthe solid.Thedashedverticallineinthefigurecorrespondstothesolid surface.FromtheleftpanelweseethatlowerCOdiffusivityyields higherconcentrationofCO₂atthesurface,whichexplainswhythe gasificationratedecreases withincreasing COdiffusivity.Therea- son for the increased CO₂ concentration at the surface is that a lower CO diffusivitymovesthe flame closerto thesurface. Since the CO₂ concentration is highest close to where it is produced, which is in the CO flame, this means that the concentration of CO₂atthesurfaceisalsohigher.StudyingthegradientsofO₂very closetothesurface(rightpanel)weseethatthecasewithhigher COdiffusivityhasasteepergradientofO₂veryclosetothesurface.

For a given O₂ diffusivity, a steeper O₂ gradient results in more transportofO2 tothesurface,and,hence,moresolidoxidation.

Finally, when increasing the diffusivity of CO₂, we see from Fig. 10 that the solid conversion rateis actually reduced. Thisis despitethefactthatthesolidoxidationrateisindependentofthe diffusivityofCO₂(see thelowerpanel ofFig.10).Thequestionis thereforewhythesolidgasificationrateisreducedwhentheCO₂ diffusivityisincreased.Theanswertothatquestionisthatforthe currentcase,whichhasasolid temperatureof1700K,theCO₂ is alwaysproducedclosetothesolidsurfaceduetotheCOflamenot beingsignificantly lifted.Consequently,an increasedCO₂ diffusiv- itywilltendtotransportCO₂awayfromthesurface,loweringthe surfaceconcentration,and,bythat,reducingthegasificationrate.

Another parameter that can influence the carbon conversion rateischemicalkinetics,both ofsurfaceandgasphase reactions.

Inthefollowingwewillproceedbystudyingthesensitivityofthe carbon conversion rate to the chemical reactivity. The reactivity isvariedby changingthepre-exponentialfactor.First,thesurface reaction rate is varied. This is done separately for the oxidation (denoted by R1) and gasification (denoted by R2) reactions. The effect of this variation on the conversion rate can be seen in Fig.12fortwodifferentsurfacetemperatures:1200and1800K.

Forthehighertemperature,theconversionrateisalmostunin- fluencedbychangesintheoxidationrate,whichisduetothefact that at such high temperatures the reaction is controlled almost purelybydiffusion.ThisisconfirmedinthelowerpanelofFig.12, whichshowsthattheoxidationrate(R1)variationshavenoeffect neither onthe contribution fromoxidation, nor onthe contribu- tionfromgasification.Atthesamesurfacetemperature,variations inthe gasificationrate(R2) haveonly aweak effecton the solid conversionrate.However, thereasonforthisisquite different,as inthiscasebothcontributionsfromgasificationandoxidationare significantly affected, as can be observed in the lower panel of Fig. 12. These two contributions are affected insuch a way that theincreaseinthecarbonconversionrateduetothehighergasifi- cationrateisalmostexactlybalanced bythedecreaseinthesolid conversionrateduetothefasteroxidation.

ForTp=1200 K,thecarbonconversionrateisdirectlypropor- tionalto the changeof the oxidationrate (R1),butdoes not de- pend on the gasification rate(R2). This is expectedsince atthis temperaturethe surfacereaction rates are controlled by kinetics, butthecontributiontothesolid conversionratefromgasification isaround two ordersofmagnitude smallerthan thecontribution fromoxidation.

The effectof thegas phase kinetics isshown inFig. 13,from whichit isclearthat thesolid conversionrateis notsensitive to the gas phase reactionrate variations, aslong as the surfacere- actions are controlled by kinetics, i.e. forT_p=1200 K. At higher particlesurfacetemperatures,thesolidconversionbecomesfaster upon decreasing the gas phase reaction rate. This is consistent

Fig. 12. upper: conversion rates as obtained for modified surface reaction rates, lower: contributions from gasification/oxidation to the CO production rate. R 1 and R 2 denote oxidation and gasification, respectively, and indicate which reaction has been modified, while ox. and gas. denote contribution from oxidation and gasification to the CO production rate.

with theoretical predictions of Libby and Blake [50] and Makino [51] who showedthat the solid conversionrateis highestinthe limit of the gas phase reaction rate approaching zero (so called

’frozen mode’),andlowest inthelimit ofvery fasthomogeneous reaction rate. This tendency can be linked to the fact that the higherthe gasphasereactionrate, themoreoxygenis consumed inside theCO-flame beforereachingtheparticlesurface,thus,the

contributiontothecombustionratefromoxidationdecreases(see thelower panelofFig.13). Furthermore,the flamecharacteristics are also directly linked to the gas phase reaction rate. In partic- ular, when the rateis increased, the flame becomesthinner and candetachfromtheparticlesurfaceorshiftfurtherfromthesur- faceifit wasalreadydetached.Thissituationcan be observedin Fig.14,which presentscontours oftheCOflameforthe casesin

Fig. 14. CO flame contours, T p = 1800 K, left: k 3,f→ 0 . 5 k 3,f, right: k 3,f→ 2 k 3,f.

Fig. 13. upper: conversion rates as obtained for modified homogeneous reaction rate, lower: contributions from gasification/oxidation to the CO production rate.

Fig. 15. Oxygen and carbon dioxide mass fraction profiles along the centerline behind the particle as obtained for different gas phase reaction rates.