Rapid change of particle velocity due to volatile gas release during biomass devolatilization

(1)

ContentslistsavailableatScienceDirect

Combustion and Flame

journalhomepage:www.elsevier.com/locate/combustflame

Rapid change of particle velocity due to volatile gas release during biomass devolatilization

Ángel David García Llamas

^a^,^∗

, Ning Guo

^b

, Tian Li

^b^,^c

, Rikard Gebart

^a

, Kentaro Umeki

^a

aEnergy Engineering, Division of Energy Science, Luleå University of Technology, Luleå, Sweden

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

cRISE Fire Research, Tiller 7092, Norway

a rt i c l e i nf o

Article history:

Received 14 July 2021 Revised 17 November 2021 Accepted 18 November 2021 Available online 21 December 2021 Keywords:

Biomass devolatilization TR-PTV

In-situ measurements Rocket effect

Non-isothermal modelling

a b s t r a c t

Ourearlierstudyshowedsignificantdifferencesinaverageparticlevelocitybetweensimulationandex- perimentalresultsfordevolatilizingbiomassparticlesinanidealisedentrainedflowreactor[N.Guoet al.,Fuel,2020].Thisindicatesthatthesimulationsdonotaccuratelydescribethephysicochemicaltrans- formationsandfluiddynamicprocessesduringdevolatilization.Thisarticleinvestigates thereasonsfor thesediscrepanciesusingtime-resolvedanalysesoftheexperimentaldataandcomplementarymodelling work. The experimentswere conductedin adowndraft drop-tubefurnace with opticalaccess, which usesafuel-richflat flame(CH₄–O₂–CO₂)toheat theparticles. Gasflow was characterizedusingpar- ticle imagevelocimetry, equilibriumcalculationsand thermocouplemeasurements. High-speedimages ofdevolatilizingNorwayspruce(PiceaAbies)particleswerecapturedandanalysedusingtime-resolved particletracking velocimetrymethods.The datawereusedto estimatethebalance offorcesand fuel conversion.Thrustand“rocket-like” motionswerefrequentlyobserved,followedbyquickentrainmentin thegasflow.Rocketingparticleswere,onaverage,smaller,moresphericalandconvertedfasterthantheir non-rocketingcounterparts.Thesedifferencesinconversionbehaviourcouldbecapturedbyaparticle- sizedependent,0-D devolatilizationmodel, correctedfornon-isothermaleffects.The resultsfrom this investigationcanprovideabasisforfuturemodellingandsimulationworkrelevantforpulverizedfiring technologies.

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

1. Introduction

ThebenefitsofbiomassasaCO₂-neutralenergysource[1]have ledtoarenewedinterestinindustrialapplicationsduringthelast decades to mitigate the global warming problem [2,3]. Biomass fuels are compatible with existing large-scale energy conversion technologies, suchaspulverized suspensionfiring. Suspension fir- ingisalsorelevantinbiofuelproductiontechnologies,suchasen- trainedflowgasification,andtechnologiesforCO₂ emissionreduc- tion,suchasoxy-fuelcombustion[4,5].However,theuniqueprop- erties of biomass, e.g. much higher reactivity, non-spherical particle morphology,anddifferent ashcomposition than fossilfuels, create aneed forfurther investigation beforethey can be imple- mentedonaglobalscale[6].

During suspension ﬁring, particles undergoa rapidconversion that canbe separatedintothefollowing threestages:drying,de-

∗Corresponding author at: Lulea University of Technology: Lulea Tekniska Uni- versitet, Sweden.

E-mail address: [email protected] (Á.D.G. Llamas).

volatilization, and char gasification/combustion [7]. The conversion behaviour depends on the heating rate, peak temperature, residencetime athigh temperatureand the local gas concentra- tionaround theparticle.Duringdevolatilization, theparticlesun- dergo significantmorphologicaltransformations andreleasemore than 70% of their initial mass in the form of vapour and gases [8]. This is an intensely heat-driven process, which is promoted by highheating rates[9].However, theapparent ratecan be re- stricted by blowing [10], evaporative cooling [11], endothermic- ityofreactions[12]andinternalconvective flowofvolatiles[13], whichtakesplacepreferentially inthedirectionofthepores[14], usually aligned withthe longest dimension. In addition, biomass particles tend to be elongated rather than spherical, which af- fects heat transfer [14]. All these heat transfer resistances make theparticlesnon-isothermalduringthedevolatilizationstagefora widerangeoffuelsizefractionsunderindustriallyrealisticcondi- tionsforsuspensionfiring.However,computationalfluiddynamics (CFD)simulationsgenerallyconsiderparticlestobethermallythin, such as below 100μm of equivalentspherical diameter[15] and especially at very high heating rates [11]. Apparent devolatiliza-

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

(2)

Nomenclature

A: area,m²

AR: aspectratio,max./min.Dimension,- B: blowingcoeﬃcient,-

C: speciﬁcheat,J•kg⁻¹•K⁻¹ C_d: dragcoeﬃcient,- d: diameter,m e: unitvector

Er: energyemittedbyradiation,J F: force,N

Lv: latentheatofvaporization,J•kg⁻¹ m: mass,kg

Nu: Nusseltnumber,- Pr: Prandtlnumber,- q: yield,kg•s⁻¹ Re: Reynoldsnumber,- T: temperature,K t: time,s

v

: velocity,m•s⁻¹ V: volume,m³ Greekletters

α

^: Absorptivity

^H^: endothermicheatofreaction,J•kg⁻¹

ε

^: Emissivity,-

μ: dynamicviscosity,Pa•s

ν

^: ^kinematic^viscosity,^m²•s⁻¹

ρ

^: ^density,^kg•m⁻³

σ

SB: Boltzmannconstant,5.6703•10⁻⁸,W•m⁻²•K⁻⁴ ^: ^Solid^angle,^sr

Subscripts

B: Bassetforce,N D: dragforce,N eff: effective

eq: equivalentforaspherewiththesamevolume

f: ﬁlm

g: gas

L: liftforce,N p: particle P: pressure,Pa r: ratio sf: Stefan T: thrustforce,N VM: virtualmass,kg vol: volatiles

tion kineticsmodels forhighheatingratesmustaccount fornon- isothermal particles,such asinthe modeldevelopedbyJohansen et al. [9,15,16], which simpliﬁes the non-isothermal problem by providing apparent devolatilizationkinetics parameters fordiffer- entsizefractions.

Theadvectionofdevolatilizationproductsfromtheparticlealso affectsthegasvelocityﬁeldandviscosityaroundtheparticle,potentially inﬂuencingviscous andpressure forces[17].Momentum exchange can also cause thrust if blowing is directional. In fact, suddenaccelerationhasbeenobservedduringthefastdevolatiliza- tionofbiomassduetodirectionalgasejection[18].Thishasbeen related to heterogeneous heating and preferential gas advection through the anisotropic pore structures of the particle [19]. Fol- lowingthislineofthought,Elfasakhanyetal.[18]andourprevious work[20]modelledthisphenomenonasathrustforcecausedbya heterogeneousreleaseofvolatiles,inaphenomenon referredhere

as“rocketing”.An alternative explanation,based onexperimental observationsofcelluloseparticlesundergoingreactiveboiling[21], couldbethepresenceofanintermediatemoltenphasethatforms a bubble witha high internal pressure that suddenly bursts and releases the enclosed pyrolysis products. Under this assumption, the devolatilization model developed by Montoya et al. [22] in- cludes bubble formation, coalescence, and rising in the molten phasetoexplaintheseburstsofgas.Furtherevidenceforthislat- termechanismissupportedbytheinspectionofparticlesthathad undergone devolatilization. Fig. 1 shows SEM images of Norway Spruce particles and their char, obtained in a drop tube furnace under pyrolysis conditions at 1200°C. Fig. 1b depicts that particles formed spherical, hollow char structures (cenospheres) with distinctholesonthesurfaceduetomelting[23,24].Thepresence oftheseholes inthe molten charstructures is commonforvari- ousbiomassspecies[24,25]andindicatestheviolentbubblingand boilingprocesses duringthe meltingprocess. Itispossible thatif meltingoccurred, volatilegasescouldhaveescapedthroughthese holes,causingthrust.Additionally,experimentsperformedbyRiaza etal.[26]showedthatelongatedparticlestendtoheatuphetero- geneously,moreintenselyintheedges.Undersuchcircumstances, localised melting could closethe pores andproduce microexplo- sionsunderveryhighheatingrates(>10⁵K/s).

Despite the pile of evidence, volatile-driven momentum ex- changebetweenparticlesandbulkgasﬂowisusuallydisregarded insimulationmodelssinceitisassumedtooccur homogeneously in all directions, therefore cancelling out thrust forces. In addition, not much has been investigated experimentally about the relevance, mechanisms and implications of the “rocketing” phe- nomenonunderindustriallyrealisticconditions,usingin-situmea- surements. Disregarding these forces can potentially lead to in- accurate estimations of particle residence time, which is cru- cial when modelling industrial burners. Furthermore, the non- isothermal behaviour, blowing effects on heat transfer, unsteady forces, morphological and density changes, etc., are often disre- garded inCFD simulations toreduce computational time andre- ducemodelcomplexitywithoutsigniﬁcantexperimentalevidence.

Thereisalackofexperimentalmeasurementsofthe“rocketing”

phenomenon forstreams of biomass particles underdevolatiliza- tion, especially with respect to its relevance, predictability, and relationship to the heterogeneous blowing. In this work, we de- scribethe“rocketing” phenomenonduringbiomassdevolatilization andinvestigateits relevanceby estimatingits frequencywithin a stream of devolatilizingparticles. We provide a simple statistical predictive model based on particle size and shape, valid for the studiedexperimentalconditions.Inaddition,weestimatethemass andmagnitudeof theforces on a“rocketing” particle during devolatilization with the aid of existing models which capture the complexheatandmasstransfer effectsandthethrust forcefrom experimentaldata.

2. Methodology

A schematic of the experimental apparatus can be seen in Fig. 2. The reactor is a drop-tube setup with optical access. A supporting flame supplies the heat and reaction environment to the devolatilizingparticles. The combustion products ofthe supporting flame were chosen to achieve a similar composition to the one found in the near-burner zone of entrained-flow gasi- fiersandoxy-fuelburners.Thesupportingflatflamewasproduced bythe fuel-richcombustionofCH₄/CO₂/O₂,whoseproducts were mainlyH₂/CO₂/CO/H₂Oandfree ofoxygen.Compositionsofpost- combustiongasweremeasuredusinggaschromatographyandde- tailedinTableS2inthesupplementarymaterial.Thebiomasspar- ticleswereinjected froma centraltubewitha streamofCO₂ gas at a feeding rate of approximately 10 g•h⁻¹. The feedstock used

(3)

Fig. 1. SEM images of Norwegian Spruce char particles from high temperature (1200 °C) pyrolysis experiments in a drop tube furnace, (a) typical particles before heating, highlighted in brown, (b) cenospheres formed during heating, highlighted in blue. Notice the holes in the cenospheres. Adopted from [26] with permission from ACS.

Fig. 2. Schematic representation of the reactor setup, all measurements are in mm.

Reprinted from [28] , with permission from Elsevier.

wasNorwegianspruceparticles (PiceaAbies) producedby aham- mermillfollowedbysievingwithasievesizeof200–250μm.The resulting particles hada highaspect ratio (AR=3.9±2.9, defined as the ratioof longestto shortest diameter). Fuelproperties and reactionconditionsforthesupportingflameandcarriergasescan be seeninTablesS1andS2inthesupplementarymaterial.Imag- ing was performed by two high-speed cameras, which collected the scatteredlight of a pulsed-laser sheet that shone across the streamofparticlesandthereactor’saxis.Imagesweresampledat 800Hzwithan exposuretime of625μs during3.75s, forafield ofviewof20×74mm.Thecylindricallensforthelaseropticswas placedapproximately1.5mfromthereactor’saxis,wellbelowthe flatflameburner.Thisarrangementwasmadeinordertotakead- vantageofthedispersionangleofthelaser,soastoilluminatethe particles enteringthereactor(SeeSection 1ofthesupplementary material for further schematics).The minimum spatialresolution oftheimagingsystemwas53.9μm.DynamicStudio6.8fromDan- tecDynamicsandMatlabwereusedtocollectandpost-processthe imagesinordertoobtaintime-resolvedmeasurementsofvelocity, positionanddimensionsoftheparticles.Thereaderisreferredto [27,28] formoredetails abouttheexperimental conditions,setup descriptionandmethodologyusedforimageprocessing.

Gasvelocitywithoutparticleswasmeasuredusingparticleim- age velocimetry(PIV),seeding thecarriergas ﬂowwith titanium dioxide particles. The same arrangement and software used for

particletrackingvelocimetry(PTV)werealsousedforthePIVmea- surements. In this study, the particle slip velocity is defined as thedifferencebetweenthevelocityofthegasflowattheparticle positionmeasured without particles andthe instantaneous particlevelocity.Thisapproximationseemsadequatesince thevolume fractionis below10⁻⁴, andthereforethe flow can be considered dilute,with negligible effectsof the particleson the gas velocity field[29].Additionally,measurementsofgasflowvelocitywithout particlesindicatethatthegasflowislaminar,andtimeoscillations canbedisregarded(seeSection 4ofthesupplementarymaterial).

Therearemethodsto performsimultaneousmeasurementsofgas andparticlevelocity,withacombinationofPIVandPTV,suchas the one describedby Khalitovand Longmire[30].However, such methodswerenotconsideredinthisstudybecausethesimultane- ousmeasurementofgasandparticlevelocitywouldhaverequired additionalseeding particles,possibly obstructingtheperformance ofthecurrentmethodologyforvolatilecloud edgedetection.Ad- ditionally,gasseedingwouldhaveaffectedtheradiativeproperties of the gas due to the incandescence of the seeding particles, as well as acting as a heat ballast. Statistics for gas measurements withoutparticles are providedin Section4 ofthesupplementary material.

A ﬂow scheme with the methodology for data analysis and modelling forthis work is presented in Fig. 3. In this chart, the centralsequenceofdatatreatmentcorrespondstoretrievinginfor- mationfromTR-PTVdata.Thisinformationwasusedtodeducea time-averaged,statisticalregressionmodeltopredicttheprobabil- ityof “rocketing” (sequence tothe right). Additionally,modelling wascarriedout using theexperimental data ofselected particles (sequencetotheleft).Furtherexplanationofthemodelsusedcan befoundinsubsequentsections.

It should be noted that this experimental studycannot track particlerotationsinceitisnotbasedonvolumetricimaging.With thissetup,onlytheprojectedareaoftheparticlecanbeaccounted for.Rotationproducesoscillatoryeffectsintheminimumandmax- imumdimensionsoftheprojectedarea. Ifthenumberofsamples ishighenough, andthe samplingrateishigherthan theparticle rotation, the average measurements of minimum and maximum dimensionsshouldcorrespond tothereal ones.Thisprocedure is similartotheoneappliedbycommerciallyavailabletools,suchas shadowgraphparticle-sizeanalysers.

However,time-resolvedmeasurements can bebiasedby parti- clerotation.Fortunately,theseoscillationscanbecompensatedby ﬁlteringtocapturethegeneraltrends.Forthetime-resolvedstudy includedin thiswork,particles showingsmalloscillations dueto rotationwere selected.Othermoresophisticated techniques,such asmachinelearningcanbeusedtopredict3Drotationfromplanar measurements.

(4)

Fig. 3. Flow chart with the process for data analysis performed in this work.

2.1. Estimationofparticlemassduringdevolatilizationfrom experimentalmeasurementsandestimatedforces

Theparticlemasswasestimatedfromexperimentaldatabased on momentumconservationduringdevolatilization. Theobjective is to provide an estimate of the mass loss during “rocketing”, alongwiththethrustforcewhichcouldexplainthisphenomenon.

Eq. (1)expresses thetransient motionofa devolatilizingparticle immersedinagasﬂow[29]:

F_B+F_D+F_L+F_P+F_V_M+mpgj+F_T=mp

d

v

p

dt (1)

Theunitvector j representsthedirectionofgravity.Theorigin ofcoordinatescanbeseeninFig.4,andparticlemovementisas- sumedtooccurwithintheplanedefinedbythexandy-axis,without out-of-plane movements. Thiswaspossible sincetherewasa sufficient numberof particles not exhibiting out of plane movements (around80% ofthe incoming particles). More information about out-of-plane movements can be found inSection 6 ofthe supplementary material. Aschematicrepresentationofthe domi- nantforcesontheparticle(drag,weight,inertiaandthrust)isde- pictedinFig.4.Theforcesareestimatedfromavailablemodelsin theliterature (Table1),usingthespatialfieldofgaspropertiesat eachparticlepositionandtime-resolvedparticleproperties(velocity,diameter,etc.).Thrustforceisdependentonthemasslossand requires the estimation of the Stefan flow from the particle surface. Fig. 4 also shows a schematicrepresentation of a cloud of volatile products beingexpelled fromthe particle witha heterogeneous Stefanflow,inthiscasewithahighervelocityoftheex- pelledgasesintheleewarddirection.ThelocalStefanflowvelocity of devolatilizationproductsatthe particlesurfaceisidentified as

v

_{s f}. Dueto the heterogeneityof theStefan ﬂowﬁeld, its integral across asurfaceenclosing theparticleis non-zeroandequivalent to a resultant velocity. Thisconsequent velocity has beennamed thereafter “effectivevelocity”,

v

_{e f f},asit is commonlyreferred to inpropulsiontheory:

v

e f f=

v

s f

r² ·dS=

v

s f·d

⁽²⁾

Thiseffectivevelocityisassumedtobethecauseforthethrust forceF_T:

FT=dmp

dt

v

e f f (3)

Propagationofuncertaintyinthecalculationofthesolutionfor m(t) can be minimized by projecting all terms of Eq. (1) in the directionofthe acceleration.Thisis because thehorizontalcom- ponentof particlevelocity isusually closetozero and,therefore, itsdifferentiationcarrieshighexperimental uncertainty.Including thedeﬁnitionofthrust forcefromEq.(3)intoEq.(1),andmulti- plyingbytheunitvectorparalleltotheacceleration,theresulting equationbecomes:

FB+FD+FL+FP+FVM+mpgj+dmp

dt

v

e f f

·ea=mp

d

v

p

dt ·ea

(4) wheree_a= ^d_dt^v^p/

|

^d_dt^v^p

|

îs^theûnit^vector ⁱⁿ^the^directionôf^theâc-

celeration.Eq.(4)isaﬁrst-orderordinarydifferentialequationfor theparticlemass,whichcanbesolvednumericallyunderthecon- dition that all the forces andother unknowns can be estimated fromtheexperimentaldata.Theparticlemasswascalculatednu- merically,usingavariable-step,variable-order(VSVO)solver,suit- ableforstiff ODEs,usingtheMatlabfunction“ode15s” forversion R2021a[31].Theinitialdensityoftheparticlewasassumedtobe 440kg/m³ [32].Tsiolkovski’s“Idealrocket” equation[33]wasused tocheckwhetherthemasslossduring“rocketing” canexplainthe observedchangeofmomentumthatisunaccountedforbytherest oftheforces(seeSection2.1.2).

2.1.1. EstimationoftheStefanvelocityandtheeffectivevelocityfrom recordedimages

TheStefan ﬂowﬁeld emanatingfromthesurfaceofthe parti- clecanbeestimatedfromtheexpansionofthecloudofincandes- centmatterthatsurroundstheparticleunderpyrolysisconditions.

The camerasensor wasableto capturethe lightfromthe incan- descentsootycloudintheabsenceoflaserillumination,indicating thatitresultedfromthecombinationofintegratedemissionalong theopticalpathandcross-sectionalscatteringfromthelasersheet.

Thisestimationcanonlybeaccurateaslongastheparticlemoves

(5)

Fig. 4. Representation of the main forces on the particle.

Table 1

Summary of models used to calculate forces.

Equation Description References

F D= ¹₂ρgC dC _d,rA p|vslip|vslip

With:

C d= ²⁴_Re(1 + 0 . 15 R e ⁰^.⁶⁸⁷) And:

C d,r= ₆_π^c^¯_d_eq c ¯= ⁶^π^d^eq√

AR²−1 ln(AR+√

AR²−1)

Drag force

Drag coeﬃcient according to Schiller-Naumann correlation

Drag ratio for a prolate spheroid with random orientation: Clift correction.

[29]

[32]

F B= _d⁹_eqαd

_ρ_g_μ_g

π

∫ t 0

dvslip t−dtτsdt With:

τs= ^ρ₁₈^p^d_μ²^eq_g

Basset force [33]

F V M= ^M₂^f(^DDt^v^g−^d_dt^v^p)

M fis the mass of ﬂuid displaced by the particle, calculated using the equivalent diameter, d eqand the local gas density ρg

Virtual mass force [33]

F P = −V eq·∇^p

V eqis the displaced volume of gas by the particle, calculated using the equivalent diameter, d eq.

Pressure-gradient force [34]

F L= FS= 1 . 61 μgd eq|vslip|(

R e _G,x·i + R e _G,y·j ) With:

R e G,x= ^d²^eq^·⁽^d^v^dx^g,x⁺^d

vg,y dx ) νg

R e _G,y= ^d²^eq^·⁽^d

vg,x dy+^d^vdy^g,y)

νg

Saffman lift force [35]

withina2D plane,turbulenceislowandvolatilematterisheated up enough tobecomeincandescent. Consequently,thismethodis onlyapplicabletothoseparticlesthatdonotmoveoutofthethin lasersheetusedforparticledetection.Themethodconsistsofap- plyinganedgeﬁltertorawimagestodetecttheedgeofthecloud ofvolatiles.Afterwards,theStefanﬂowvelocity canbecalculated from theexpansion of the edge,relative to the particle centroid, betweenconsecutiveframes.Aschemeofsuchaprocessisshown inFig.5.

IftheStefanﬂowvelocityisassumedtoemanateradiallyfrom theparticlecentroid,itslocalvaluecanbeexpressedasafunction

of

θ

^,^whichîs^theângle^formed^with^the^verticalâxis:

v

s f

( θ )

⁼

^L^t 1

θ

^e^θ ⁽⁵⁾

UsingEq.(2),

v

s f(

θ

)^can^be^used^to^calculate

v

e f f,asinEq.(6):

v

e f f=

v

s f

( θ )

^·^eθd

θ

⁽⁶⁾

Themethodologyherepresentedisonlyapplicableunderlami- nar conditionsforsmall particleReynolds numbers.Thesecondi- tions differ signiﬁcantly from those found under realistic indus-

(6)

Fig. 5. Estimation of Stefan ﬂow velocity from edge detection of the incandescent cloud of volatiles around the particle. Green crosses identify the particle centroid. These images and edges have been obtained from real experimental data, with t = 3 ms between t1 and t2, using a Sobel ﬁlter for edge detection. Contrast has been adjusted for easier visualization.

trialpulverizedburners,suchasentrained-ﬂowgasiﬁcation.How- ever,theapplicabilityforindustrially-realisticconditionswouldbe more plausiblewith a similar methodologyto the one proposed, involving stereoscopic measurements on an electrodynamic ther- mogravimetricanalyser.Thiswouldenableveryhighheatingrate, controlled atmosphereandsimultaneous visualisation, such asin theworkbyE.Bar-Zivetal.[34]andthemorerecentbyBiaginiet al.[35],incombinationwithpowerfultoolsforthemeasurement ofvolatilegasvelocity outofthe particle,suchasinterferometry, asit wasperformedby Lycksamet al.[36] orLIF(Laser Induced Fluorescence).

2.1.2. Theidealrocketequation

To check whether the velocity changes experienced by a particlecan beexplained bythemass estimatedasdescribedabove, theequation forthemovement ofanidealrocketwithoutdragis used for comparison.This equation wasdescribed by Tsiolkovski [33],andrelatesthechangeofvelocityoftheidealrocketwiththe

“effectivevelocity” ofthepropelledgasesandthechangeofmass duringthisprocess:

v

p=

v

e f f·ln

m0

m_f

(15)

where

v

p isthe increase invelocity of therocket (in thiscase, the particle),

v

e f f is theeffective velocity of the gases propelled out of therocket (volatiles) and ^m_m⁰

f isthe ratioof initialto ﬁnal massoftherocket.

2.2. Estimationofparticleconversionfromexperimentaldataand non-isothermaldevolatilizationmodels

Conversion of a biomass particle during devolatilization was estimated, considering the effect of changes in particle size and shapeaswellasgastemperatureonheattransfer.Non-isothermal heating and blowing effects were also taken into account. The aim was to ﬁnd out which model reproduces the estimated masslossobtainedfrommomentumconservationmoreaccurately (Section 2.1), andto providean estimate of the heatingrateand product composition during conversion of a “rocketing” particle.

The spatialﬁeld of gas properties and time-resolved particle di- mensionswereusedtoestimate thetime-dependentheattransfer ratebyconvectionandradiationfromavailablemodels,usingthe energy conservation equation. The obtained particle temperature is thenusedto updatethekinetic parametersofa 0Dconversion model,takingintoaccountnon-isothermalheating.

Particle temperatureis calculated foreach experimental point fromtheenergyconservationequation:

mpCp

dTp

dt =−hAp

(

^Tp−Tg

)

−Er+q_vol

^Hv^ol (16) Inthisequation,particlemassandparticletemperaturecanbe obtainedfromthe simultaneoussolution ofthechemical kinetics ofTable 2andEq. (16). Thermochemical gaspropertieswere obtained from interpolation from temperature measurements from theNISTChemicalKineticsDatabase[37].Theheatofdevolatiliza- tion,^Hvol,wasobtainedfrom[10].Initial densityoftheparticle wasassumedas440kg/m³ [32]andtheinitialtemperatureofthe particlewasassumedtobe300K.

External convective heat transport is usually modelled using heat transfer correlations for spheres using a Nusselt number correlation (Nu=hD/

λ

^). ^One ^commoncorrelation is that of Ranz- Marshall[38]:

Nu=2+0.6Re¹²Pr¹³ (17) wherethermochemicalpropertiesatfilmconditionareusedtocal- culateReynoldsandPrandtlnumbers.However,thepresenceofan outflowofdevolatilizationproductsfromtheparticlesurfacemust be takeninto account witha correction to the estimationof the film temperature, using the 1/3 rule, as suggested by Yuen and Chen[17]:

Tf=Tp+Tg−Tp

3 (18)

Toaccount forblowingeffects, severalauthorssuggest a Nus- seltnumbercorrectionbasedonablowingcoeﬃcient(alsocalled mass transfer ratio orSpalding heat transfer number) [10,38,39].

TheNusseltnumbercorrectionfromtheheattransfernumbercan be obtained from the model developed by Renksizbulut & Yuen, correctingtheRanz-MarshallcorrelationfromEq.(17):

Nur= Nu

(

¹⁺^B

)

⁰^.⁷ ⁽¹⁹⁾

whereBistheblowingcoeﬃcient,calculatedas:

B=Cp

(

^T^p−Tg

)

L_v (20)

where Lv is the latent heat of vaporization of a liquid droplet, which can be approximatedby using Trouton’srule forthe boil- ingtemperatureofthemoltenphase:

L_v

T_boiling ≈85−88 J

Kmol (21)

Inthiswork,the boilingtemperatureofthemolten phasehas beentakenasthat ofLevoglucosan(384°C)[40]. Gasradiationis

(7)

Table 2

Devolatilization kinetics parameters.

Parameter Unit Constant References

Small size fractions ( < 112 μm): High-temperature kinetics, one-step reaction mechanism

A 1 s ⁻¹ 8.56 ×10 ¹⁰ [16]

Ea 1 KJ •mol ⁻¹ 171.8

Medium size fractions (112–616 μm): High-temperature kinetics, one-step reaction mechanism

A 2 s ⁻¹ 3.99 ×10 ⁹ [16]

Ea 2 KJ •mol ⁻¹ 162.3

Large size fractions (616–2000 μm): High-temperature kinetics, one-step reaction mechanism

A 3 s ⁻¹ 2.62 ×10 ⁶ [16]

Ea 3 KJ •mol ⁻¹ 118.7

Kinetic parameters for the two-step reaction mechanism, valid for low heating rates.

A V s ⁻¹ 1.11 ×10 ¹¹ [41]

Ea V KJ •mol ⁻¹ 177

A T s ⁻¹ 9.28 ×10 ⁹ [41]

Ea T KJ •mol ⁻¹ 149

A V,2 s ⁻¹ 4.28 ×10 ⁶ [42]

Ea V,2 KJ •mol ⁻¹ 108

Fig. 6. Reaction models (a) One-step reaction mechanism, used at high temperatures and heating rates (b) Two-step reaction mechanism, based on the Broido-Shaﬁzadeh scheme, valid for low heating rates and lower temperatures.

especially important in environments withhigh partial pressures of radiating gases(CO₂ and H₂O) and can be obtained fromthe HottelcorrelationsformixturesofCO₂andH₂O,correctedformu- tual radiation[41].Back-radiationfromincandescent particlesand sootycloudscanbedisregardedunderverydiluteﬂows,although they can contribute to maintaining particle temperature. Particle emissivity wasassumed as0.3[42]. Finally,the energyexchange byradiationbetweengasandparticlescanbeobtainedfrom:

Er=−

σ

sb

ε

gT_g⁴−

ε

p

α

gT_p⁴

(22)

whereEristheenergytransferredtotheparticlebyradiation.

The estimated temperature fromthe energy balance for each experimental pointisusedtodeterminethekineticcoeﬃcientsof a devolatilizationmodelin ordertocalculate theconversion.The biomass particles during the devolatilization stage under indus- trially realistic conditionsforsuspension ﬁring isnon-isothermal.

Therefore, itis necessaryto considera non-isothermal correction andtheuseofkineticparametersoptimizedforhightemperatures andheating.Thekineticparametersforazero-dimensional,single ﬁrst-order reaction(0D SFOR)model(Fig.6a)by Johanssen etal.

[16] were used in this study.Kinetic parameters were optimized with the results of an experimentally validated model for single particles,wherelocalconversionissolvedoverparticleradius,and intermediate species are considered using the reaction model of Fig. 6b. Kinetic parameters are provided for three different size fractions:small(below100μm),medium(100to600μm),andbig (above600μm)particles,asshowninTable2.

Particletemperatureandconversionproductsarecalculatednu- merically using the Runge-Kutta method. Char gasiﬁcation reactions have not been included to avoid making further assumptions sincetheaimofthearticleistostudydevolatilisationreac- tions.Thisassumptionseemssafe,sincegasiﬁcationreactionshave a much higher characteristic time than pyrolysis for the studied cases: _t^t^pyrolisis

gasi f ication ≈10³−10⁴ forparticlesizesfrom 50to 1000μm

Fig. 7. Sequence of images of a particle (black) exhibiting the jet effect caused by the sudden release of volatile matter (grey area). Arrows represent particle velocity and its length correlates to velocity magnitude, indicated with the reference vector in the upper left image.

[43]. Oxidation reactions are very unlikely since the atmosphere aroundtheparticlesisheavily reducingduetothefuel-richcom- bustion products from the ﬂat ﬂame. No assumptions on homo- geneous chemistry have been madesince only devolatilisation is considered,anditwouldexceedtheassumptionsforthisarticle.

3. Resultsanddiscussions

3.1. Introductiontothe“Rocketing” phenomenon

Figure 7 illustrates the aforementioned “rocketing” phe- nomenonwitha sequenceofimagesofa biomassparticleunder- goingdevolatilization.Theparticlevelocityvectorsthataresuper- imposedontheimageshavebeendeterminedwithPTVandindi- catethe velocity magnitude. The raw images havebeen inverted to enhance contrast, andthus, all radiatingmatter in the visible

(8)

Fig. 8. Particle trajectories for: (a) all time-resolved particles (b) “non-rocketing”

particles (c) “rocketing” particles.

spectraappears darkerthanthebackground.Thetimestep ontop oftheimagesindicatestheresidencetimeoftheparticlefromits entrance tothe reactor.Theparticle exhibitedlateralmotionina similarmannertothatofarocket,apparentlycausedbyafastre- lease ofvolatilematter. This phenomenon began after23.0ms of residence time forthisparticle,when a smallcloud of incandescent matter emerged behindthe particle.The appearance of this cloudwasaccompaniedbyasuddendeflectionoftheparticletra- jectoryintheoppositedirectiontothereleaseofvolatiles.Judging bythewedgeshapeofthevolatilecloudneartheparticle,there- leaseofvolatilesseemedtocomefromanarrowgapatthesurface of the particle. Thiseffect, which we refer to as“rocketing”, has beenobservedrepeatedly,forasignificantfractionoftheparticles, in all the experiments performed forthis work. The direction of the deflection was random, sometimes directing the particle up- stream. Theresidencetime atwhichthe phenomenon tookplace variedfromparticletoparticle.

3.2. Thecollectivebehaviourof“rocketing” versus“non-rocketing”

particles

3.2.1. Categorizationofbehaviourandparametricanalysis

The “rocketing” phenomenonwaseasily recognizablefrom di- rectobservationofthetime-resolvedparticletrajectoriesandcom- parison with the recorded experimental images, as it was seen in Section 3.1. However, not all the detected particles remained within the thicknessof thelaser sheet, andthe trajectorieswere rarely complete. Therefore, velocity and properties of“rocketing”

and“non-rocketing” particles havebeenextractedfromastatisti- cally signiﬁcant number of particles, which remained within the laser sheet throughout the ﬁeld of view. Then, these properties have been averaged at each residence time for each category.

Fig. 8 presents the trajectories forthe particles which remained within the thickness of the laser sheet along the ﬁeld of view.

Trajectories in red and blue correspond to “rocketing and “non- rocketing” particles,respectively.Withintheonesexhibiting“rocketing”, it is possible to see large motion deﬂections, presumably duetotheejectionofanarrowjetofgasfromtheparticlesurface that givesrise to a netthrust. Itcan be alsonoted that there is migrationtowardsnegativeradiiforalltrajectories.Thiscanbeex- plainedbythenon-axysimmetryofthegasﬂow,causedbyamis- alignment in the carriergas injection line.Flow inhomogeneities causingliftcanbedisregarded,asitisfurtherdiscussedinsection

3.3.2 of this manuscript. Further discussion on this topiccan be foundinSection4ofthesupplementarymaterial.

Figure 9 depictstime-averaged properties for “rocketing” and

“non-rocketing” particles, namely: acceleration (Fig. 9a), slip velocity(Fig.9b), particlevelocity (Fig.9c),andparticledimensions (Fig.9dtof),includingvolume,minimumdiameterandaspectra- tio.Continuouslinesrepresentmeanvalues,andshadedareasin- dicate standarddeviation aroundthe mean.Graphswith theraw datausedforthesegraphscanbefoundinthesupplementaryma- terial(FigureS1).

Figure9arepresentsaverageaccelerationversusresidencetime andrevealsintense ﬂuctuationsinthe accelerationof“rocketing”

particles. Meanwhile, the average acceleration of “non-rocketing”

particleschangedsmoothlywithasmallstandarddeviation.These ﬂuctuationsarecausedbythesuddenaccelerationsaccompanying the “rocketing” phenomenon. The average accelerationof all tra- jectoriesattheentrancetothereactorwasslightlylowerthanthe accelerationofgravityandwithhigherdispersionforthe“rocketing” particles. Thiscould be attributed tothe smaller sizeof the

“rocketing” particles,thereforepresentinglessgravimetricforce,as itisdiscussedinfurthersections.

Figure9brepresentsaverageslipvelocityasafunctionofresi- dencetime.Averageslipvelocityattheentrancetothereactorwas verysimilarandhadananalogousstandarddeviationforbothcat- egories ofparticles. As residencetime increased, the averageslip velocity for“rocketing” particles tended towards zero,indicating that“rocketing” particleswereentrainedinthegasﬂowfasterthan the“non-rocketing” ones.Thedifferencesinaverageslipvelocities between“rocketing” and“non-rocketing” particles startedappear- ingafterapproximately20msofresidencetime.

Figure 9crepresentsaverage particlevelocity versus residence time. Here,it can be seen that the“rocketing” particlestravelled ata signiﬁcantlylower velocity, compared tothe “non-rocketing”

ones,alreadyfromtheentranceofthereactor.Thisinterestingre- sultindicatesthat,fromatime-averagedperspective,thetendency towards particle entrainment in the gas ﬂow is a more relevant factorto representthemotionof the“rocketing” particles, rather thanthe shortbutintensevelocity ﬂuctuationsthathelp identify it.

RegardingmorphologypresentedinFig.9d–f,particlesthatex- hibited “rocketing” motions were on averageslightly smallerand consistently less elongated when they entered the reactor. Dur- ing conversion, intense shrinking wasnoticeable for both “rocketing” and “non-rocketing” particles. Additionally, later entrainment of “rocketing” particles in the gas ﬂow was accompanied by spheroidization.Thisis indicated bya decrease inparticleas- pect ratio and a simultaneous increase of the minimum diameter. It is unknown whether“non-rocketing” particles also tended tospheroidize,sincethey didnotremain longenoughwithin the ﬁeldofview.

As summarized above, already atthe entrance to the reactor, the “rocketing” particles were smaller, rounder and travelled on average slowerthan the“non-rocketing” counterparts. Thisresult hasbeenusedtodevelopapredictivemodelofthe“rocketing” effect basedonthe originalsize andshapeofthefeedstock, which canbefoundinsubsequentsections.Moreover,asconversionpro- ceeded,“rocketing” particlesgot entrainedinthegasflowasthey shrankandspheroidized.Incontrast,for“non-rocketing” particles, gasflow entrainment didnot happen within the field ofview of thecamera.Therefore,fromatime-averagedperspective,the“rocketing” phenomenon affectsthe response time of theparticles in a fluid,allowingthemtogetentrainedfasterinthegasflow.

3.2.2. Frequencyof“rocketing”

ResultsfromSection3.2.1.indicatethatparticlesexhibitingthe

“rocketing” effecthad a smaller minimum diameter anda lower

(9)

Fig. 9. Time-averaged properties of “rocketing” and “non-rocketing” trajectories.

Fig. 10. (a) Scatter plot of with all particles detected in the experiment, colorscale based on minimum diameter (b) Projection of a Loess ﬁt of the points deﬁned by velocity, distance from burner outlet and minimum diameter. Overlayed on top of this graph: data from Fig. 9 c.

aspectratiothanthenon-rocketingones.Theresultsfromtheav- eraged propertiesofthesetrajectoriescanonly representqualita- tively thedifferences betweeneach behaviour.To drawquantita- tiveinformationonthedifferencesbetweenbehaviours,itisneces- sarytoanalysealargernumberofparticles,usingatime-averaged approach.However,not allparticlescouldbe followedfortheen- tire ﬁeld-of-view dueto: (1) out-of-plane motions(see Section 6 of the supplementary material), (2) hindered particlerecognition in regions withhighlight intensity, (3)particles moving closeto eachother causingfailureintrackreconstructionand(4)thesud- den motions of“rocketing” particles preventing trackreconstruc- tion.Instead,allsamplesdetected fromtheexperimenthavebeen takenintoaccount,includingallrepeatedsamplesfromeachpar- ticle.

Figure 10shows(a) thescatter plotwithall samplesdetected during the experiment, with the colour scale indicating particle size, and (b) the 2D projection of the surface ﬁt for the z-axis.

Overlayedontopofgraph(b),theaveragesoftime-resolvedparti- clevelocityversus distancefromtheburneroutletfor“rocketing”

and “non-rocketing” behaviourscan be found. To allow compari- sonwithqualitativetime-resolveddatafromsection3.2.1,velocity

contoursforaveragedtime-resolveddatafromFig.9havebeenin- tegratedtoberepresentedagainstdistancefromtheburneroutlet.

The scatter plot in Fig. 10a shows two distinct velocity behaviours diverging from 20mm from the burner outlet. Figure 10bindicatesthat theupperbranchcorrespondstothebehaviour of “non-rocketing” particles, while the lower branch corresponds to the “rocketing” ones. A signiﬁcant amount of samples were detected during the experiment exhibiting the “rocketing” phenomenon.Therefore,theparticlesbelongingtothesesampleswere responsible for the discrepancies with our previous simulation work[27] duetotheir tendencytowardsgettingentrained inthe gasﬂow.FurtherinformationcanbefoundinSection7ofthesup- plementarymaterial.

The clear velocity branching caused by “rocketing” particles, seen after 20mm from the burner outlet in the scatter plot of Fig.10a,andthehighnumberofdetectedsamplesthroughoutthe experiment, allows a statistically signiﬁcant quantiﬁcation of the frequencyofthetwobehaviours.Thecategorizationhasbeendone byhistogramsofnormalizedparticlevelocityfrom0to5mmfrom theburner outletandfrom40to 45mm fromtheburner outlet.

Allvelocitiesusedforthesehistogramsweretranslatedtotheini-

(10)

Fig. 11. Histograms and tri-normal ﬁts of particle velocity at (a) 0–5 mm from the burner outlet and (b) 40–45 mm from the burner outlet. Dashed lines indicate average velocity from time-resolved measurements at the beginning of each range of distances from the burner outlet. Np stands for number of particles.

tial value of the range.To avoid the inclusion of repeated measurements ofthesameparticle trajectoryin thefrequencyanaly- sis, multiplehistogramswereobtainedandaveraged forthe min- imumdisplacementpossible,giventhevelocitiesobserved.Fig.11 representsthehistogramsofparticlevelocity for5mm belowthe burner outlet and40to 45mmbelow theburner outlet. The ve- locityobtainedfromthelimitednumberofthecompletetrajectory (section3.2.1)isalsorepresentedwithdashedlines.Thehistogram of velocitiesat40–45mmfromtheburner outletpresenteda bi- modaldistribution.Toestimatethefrequencyofthe“rocketing” ef- fectwiththeconsiderationoftheoverlapbetweenthetwomodes, thehistogramwasdeconvolutedtoatrinomialdistribution,forcing themodestobetheclosesttotheaveragetime-resolvedvelocities whilemaximizingtheR².Theareasunderthesecurveswereused to obtain the fractions of “rocketing”,“non-rocketing”, andparti- cles with unclear behaviour due to overlap.The trinomial distribution wasfit tothehistogramdatausingthe“fit” function with a fittingequation in Matlab. An iterativeseekwas performedso as two Gaussians hadmeans asclose aspossible to the average time-resolved velocitieswhile maximisingtheR².This wasmade bychangingtheoptionsfortheMatlabfittingfunction(upperand lower limitsofthefunction “fit”).Noother thresholdswere used.

The resultingGaussianswere classiﬁedinto“rocketing”,“unclear”

and“non-rocketing” basedontheirdistancestotheaveragetime- resolvedvelocities.

AsitcanbeseenfromthedatapresentedinFig.11,ataround 40–45mm from the burner outlet, the percentage of “rocketing”

particles wasat least37%.Many oftheseparticles also exhibited anintensedeceleration,asitcanbeattestedbycomparingtheav- erage time-resolved velocity. The bimodal distribution from0 to 5mm fromtheburner outletdoesnot allowenough accuracyfor quantifyingthefrequencyofthe“rocketing” phenomena,giventhe number of unclear particles. However, the ﬁtted distributions of

“rocketing” and“non-rocketing” particlesaresuﬃcientlyseparated fromeachothertoassignprobabilitiesof“rocketing” basedonpar- ticlevelocity.Probabilityplotscanbefoundinthesupplementary material.

The probability of“rocketing”,obtainedfromFig.11a,wasex- pressedasafunctionofminimumdiameterandaspectratioforall theparticlesusedforthehistograminFig.11a.Thesametransla- tiontocoincidewiththeinitialvalue oftherangewasperformed in the same way as withvelocities for the previous histograms.

Thefunctionwasexpressedasthesumoftwologisticregressions.

Figure 12representstheresultsofthisestimation,indicatingthat for a spherical particle with 200μm of diameter, the probability

ofrocketingisapproximatelyp≈20+15=35%.Moreinformationon thesurfaceﬁtcanbeseeninthesupplementarymaterial.

3.2.3. Discrepanciesbetweentheestimatedandmeasuredvelocityat theentrancetothereactor

Thecauseforthedifferentvelocitiesattheentrancetothere- actorfor“rocketing” and“non-rocketing” particles isintriguing.It must be examined whether it can be explained solely based on size and shape differencesor if other phenomena such as addi- tionalforcescouldbeinvolved.Thisstudymustconsiderthesharp temperaturegradientattheentrancetothereactor,whichwillaf- fect dragforces in its immediate surroundings. For this, the velocity of the particles from the feeding mechanism, through the conveyingtubeuntilthe entrancetothe reactor(y=0),hasbeen estimatedfromthegeometricaldata,carriergasvelocityandtem- peratureproﬁle,assumingonlydragandweightastheonlyforces actingontheparticles.

Figure13arepresentsthedragforcetoweightratiobeforeand afterenteringthereactorfromthefeedingtube.Itindicatesadras- ticreduction ofthe dragforce comparedto theweight upon the injection from the feeding line to the burner. It was caused by the changein carrier gasproperties while heatingup. Therefore, weightwasmostprobablythemaincontributingforcetomomen- tum at the entrance to the reactor, which explains why the average accelerationof the particles at the entrance to the reactor wasclosetogravity.Fig.13bshowsthecomparisonbetweenesti- matedparticlevelocityandmeasuredoneattheoriginofcoordi- nates.Calculationofthebalancebetweendragandweightfailsto predictthevelocity of“rocketing” particlesattheentranceofthe reactorbasedontheirsizeandshape.

Potential reasons for this discrepancy are either: (1) drag force isgreatly underestimatedor(2)an additionalforce orphe- nomenonispresent.Giventheintensechangesintemperatureand gas properties,this discrepancycould be explained by the effect of devolatilization products, such as changes in the gas properties orthe momentum exchange betweenvolatiles, gas andpar- ticles during conversion. Scenario (1) is not very plausible since mostmodels correcting forStefan flow (gas emanatingfrom the surfaceoftheparticle) predictalubricatinglayeraroundthepar- ticle,causingadecreaseinthedragcoefficientand,therefore,the dragforce [44].Only one modelpredicts an increase in thedrag force due toStefan flow, but it isexpected to occur undercom- bustion [45].Forscenario (2),additionalforcessuch asthose de- scribedinthetheorysectioncouldberelatedtothisphenomenon.

Subsequent sectionswill investigate theestimated forces ontwo

(11)

Fig. 12. Stochastic model for “rocketing” probability as a function of (a) Minimum diameter and (b) aspect ratio.

Fig. 13. (a) Estimated drag-weight ratio before and after entering the reactor (b) Estimated versus measured particle velocity at the entrance to the reactor. For the measured values, errorbars indicate standard deviation. For the estimation, errorbars indicate the minimum and maximum values obtained taking into account the dispersion of the raw data.

particles exhibiting “rocketing” and “non-rocketing” phenomena and try to distinguish which force could be related to such an event.

3.3. Behaviourofindividual“rocketing” particles

3.3.1. Particlemotionofa“rocketing” anda“non-rocketing” particle Twoparticles havebeenchosen asrepresentative examples of the“rocketing” and “non-rocketing” behaviours.Fig.14containsa series of snapshots of these two particles at different residence times. Inthese images,particles always appear whitedue tothe scattered light from the laser. Therefore their pixel intensity is not related to temperature. As devolatilization proceeds, volatile matter is released, which eventually becomes incandescent and appears in the images as diffused grey areas. The pixel intensity of thesecloudsis a combinationof Miescattering andradi- ation inthe visible spectrum. Underneath each image, it is indicated whethertheparticles were accelerating, “rocketing”,orde- celerating. Their behaviour is identiﬁed as “rocketing” and “non- rocketing”. At their entrance intothereactor, the“non-rocketing”

particle was substantially more elongated and bigger in volume thanthe“rocketing” one.Forthe“rocketing” particle,thecloud of incandescent volatiles appeared ataround 15ms (notincluded in the set of images). This cloud followed the particle as it moved downstream, stretching vertically. After 40ms, theparticle began toescapethecloudinthedirectionofgravity,“rocketing” violently

out of itat around 54ms. The particlecan be seen escaping the cloud ofvolatilesfromthe bottomin theimage at63.8ms.After rocketing,theparticlemoveddownstream withoutavisiblecloud ofvolatilegasessurroundingitandappearedslightlyswollencom- paredtothesizeduringrocketing.Forthe“non-rocketing” particle, thecloudofincandescentvolatilesalsoappearedataround15ms.

However,itwaslessintense(ﬁgurenotincludedinthissetofim- ages),eitherfromvolatilesbeingreleasedinlesseramountsorby being at a colder temperature. The cloud of the “non-rocketing”

particle achieved its maximum size at around 40ms. Eventually, similartothe“rocketing” particle,the“non-rocketing” onealsoes- caped the cloud of volatiles and did not show more signs of a cloudaroundit.

Figure 15 shows the particle velocity and effective velocity against residence time for these two particles. Coloured back- groundsindicatedeceleration,accelerationand“rocketing” stages, and average acceleration during these stages is characterised by textonthegraph.Fortheparticlevelocityplots,apiecewisepoly- nomialﬁtwith95%conﬁdenceintervalsisaddedtothegraphfor easier interpretation of the results. Additionally, the gas velocity measured without particles is included after being converted to theparticle frameofreference. Forthe effectivevelocity, thereis a lack of experimental pointsbefore 15ms ofresidence time for the“rocketing” particle.Thesepointshavebeenobtainedbylinear extrapolationofsubsequentdataandare indicatedinthe plotby dashedlines.

(12)

Fig. 14. Snapshots of the particles studied during the different stages.

Fig. 15. Different stages of particle motion for a rocketing and a non-rocketing particle. The particle velocity is positive for motion in the direction of gravity.

The“rocketing” particledeceleratedimmediatelyafterentering thereactor.Giventhatgasandparticleareinco-ﬂow,itwouldbe expectedthat thegasaidedparticlemotionattheentrance.Such deceleration could also be attributedto thrusting dueto drying, startingneartheburneroutlet.Theeffectivevelocityofthesegases wouldnotbedetectedbythemethodpreviouslyexplainedsinceit requiresanincandescentcloudofvolatiles.

Afterwards,theparticlecontinuedalmostatfreefall.Apossible explanationforthisisthepresenceoftheStefanﬂowfromthepar- ticle surfacedue todevolatilization. Thiscould havecaused drag reductionby creatingathickenedlubricatinglayer,andiftheSte- fanﬂowisheterogeneous,thrust/dragcompensation.Bothpossibil- ities are plausiblesincethereis asmalleffectivevelocity present duringthisstage,ascanbeseeninFig.15.Eventually,“rocketing”

was observedwitha sudden increase in themagnitude ofaccel- eration, explained by thethrust causedby the releaseof volatile gases.

Lateron,the“rocketing” particledeceleratedsuddenlyandfol- lowed the stream of gas. By contrast, the “non-rocketing” parti- cledeceleratedonlyslightlytowardstheendofits trajectory.The

“Rocketing” particle wasaccompaniedby a noticeableincrease in the effective velocity, which caused a net thrust force. Interest- ingly,inthe caseofthisspeciﬁc particle,thrust initially opposed theparticle’s motion,duetovolatiles beingreleasedin thesame direction as particle velocity, temporarily slowing it down. After this, the direction of the effective velocity quickly changed, and the particlewas propelleddiagonally. This indicates that forsuf- ﬁcientlyhigheffectivevelocities,theparticlecouldhavebeenpro- pelledupstream.Particlerotationcanexplainthegeneraltendency ofparticles forbeingpropelled diagonally. Bycontrast,the effective velocity observed in the “non-rocketing” particleis substan- tiallylower.

3.3.2. Estimatedmasslossduring“rocketing” fromestimatedforces Figure 16presents theestimatedparticle massduring conver- sionforthesame “rocketing” particlefromSection3.3.1. Asensi- tivityanalysistothemodelappliedinthissection,andanuncer- taintyanalysistoallderivedparameterscanbefoundinSection3 ofthesupplementary material.Thisresultwasobtainedusingthe methodology describedinSection 2.1. Note that thevertical axis

(13)

Fig. 16. Comparison of mass obtained from equilibrium of forces for a “rocket- ing” particle. Overlayed to this plot, the estimated mass loss required for the ideal rocket equation (blue) and the ﬁnal calculated mass from experimental measurements without considering thrust force.

is on a logarithmicscale. Forcomparison,a graph withthe esti- matedmasslossduringrocketing,calculatedusingtheidealrocket equation fromSection2.1.4, hasbeenoverlayedwithbluedashed lines.Thesolutiontotheidealrocketequationusedasinitialmass is givenbythe estimationrightbefore“rocketing”.Thisis shown with blue dashed lines in Fig. 16. In addition, a direct solution of Eq. (4) without thrust force has been obtained after rocketing (represented withred dashed lines). Thiswas done to check whetherparticle motionafterrocketingcould beexplained with- outtheneedforathrustforce.

Figure 16 indicates that the mass of the particle during the rocketingphenomenonisverysmallcomparedtoitsoriginalvalue attheentrance.Therefore,the“rocketing” phenomenon took place forthisparticleatanadvancedstageofconversionwherethemass oftheparticlehadlostmorethan90% ofitsoriginal value.How- ever, during “rocketing”, the particle loses 80% of the remaining masswhilebeingquicklypropelledaway.Fig.16alsoindicatesthat the idealrocket assumption isinagreement withthe increase of momentumexperiencedbytheparticleduetotheamountofmass

releasedandthatnomorethrustneededtobeassumedafterrock- etingfortheestimationoftheparticlemass.

Figure 17 shows the forces on the particle during conversion that were obtained from the solution for the mass, using the methodologypresentedin Section2.1.The initialdeceleration re- quiresaninitialthrustforceonthe“rocketing” particle.Thisforce wouldbepresenteven atvery smalleffectivevelocities,anditis almost not present for the “non-rocketing” particle. Thrust force alsoexplainsthesuddenaccelerationduringthe“rocketing” effect.

Other forces, such asSaffman and Basset force are relevant only when the slip velocity becomes zero. Although these forces are irrelevant duringmostofthe residencetime,they might become importanttoavoidintegrationerrorsfromthesolutionofEq.(4), duetotheinstantaneouszerovalueofthedragforce atthezero- crossing. However,they cannot bethe only causeforthe sudden lateraldisplacementsorthelateraldeviationofthetrajectories.

The rocketing particle began thrusting towards its windward side whenslipvelocity wasalreadylow,decelerating theparticle tothepointofalmostzerodragbeforeeventuallybeingpropelled diagonally.Oncetheparticlebeginstogain velocityduetoitsin- creasedmomentum,dragincreasesagain,slowingdowntheparti- cle.

3.3.3. Estimatedparticletemperature,compositionanddensity during“rocketing” fordifferentdevolatilizationmodels

Figure18showsthecalculatedmassandtemperatureusingthe methodologyfromSection2.2.forthesame“rocketing” and“non- rocketing” particlesstudied inSection 3.3.1 and 3.3.2.A sensitiv- ityanalysisofthe modelused inthissection andan uncertainty analysis to all derived parameters can be found in Section 3 of thesupplementary material.Apparentdevolatilizationkineticsfor smallparticlesagreesbestwiththe“rocketing” particle,whilecon- version ofthe“non-rocketing” one isbetter approximatedby the kineticsmodelformediumparticlesize.Particletemperatureesti- mationsusingthekineticsmodelforsmallparticlesshowamuch fasterrisethanothermodelsandfollowgastemperature.

Figure19representstheestimatedyieldofproductsandparti- cledensityforthe“rocketing” and“non-rocketing” particles.These were obtainedusing thekinetic models that agree best withthe mass in Fig. 18. Modelling results of product yield and density shown in Fig. 19 indicate that most of the mass was lost at a moreorlessconstantdensityinitially.Thisresultisinaccordance withHolmgrenetal.[46].Duringthesubsequent“rocketing” stage, the“rocketing” particlebecameverydenseandturnedﬂuffy(with verylowdensity)afterthrusting.Rocketingcouldbepossibledue tosomeunconvertedmaterialorgastrappedintheparticle’score dueto thenon-isothermal heatingorbybubbleformationcaused

Fig. 17. Estimated forces on the particles. Background colours indicate the regions of acceleration, deceleration and rocketing.

(14)

Fig. 18. From left to right and top to bottom: (a,b) Estimated mass from equilibrium of forces and thermochemical models for a “rocketing” and “non-rocketing” particle (c,d) Estimated particle and gas temperature from thermochemical models and experimental measurements, respectively.

Fig. 19. From left to right and top to bottom: (a,b) Estimated yield of devolatilization products, using the models best ﬁtting to particle conversion (c,d) Estimated particle density from kinetic models and balance of forces.

by ametaplasticstage.However, thelatterishardtopredictwith the available kinetic parameters for metaplast formation, which have been measured at much slower heating rates. The product compositionshowninFig.19doesnotindicateanyyieldofmeta- plast for the “rocketing” particle. However, the modellingresults for the “non-rocketing” particle indicated a signiﬁcant formation ofmetaplast.

3. Conclusions

Thisstudyprovidesexperimentalevidenceofthesuddenaccel- erationoffuelparticles,referredtoas“rocketing”,duringbiomass devolatilization.ThefeedstockstudiedwasNorwegianSprucepar- ticlesundersuspensionfiringconditions,withanatmospheresim- ilar to that encountered under entrained flow gasification. The