# A Large Eddy Simulation Study on the Effect of Devolatilization Modelling and Char Combustion Mode Modelling on the Structure of a Large-Scale, Biomass and Coal Co-Fired Flame

Authors: Miriam Rabaçal; Mário Costa; Michele Vascellari; Christian Hasse; Martin Rieth; Andreas M. Kempf
Journal: Journal of Combustion
 (2018)

---

## Abstract

This work focuses on the impact of the devolatilization and char combustion mode modelling on the structure of a large-scale, biomass and coal co-fired flame using large eddy simulations. The coal modelling framework previously developed for the simulation of combustion in large-scale facilities is extended for biomass capabilities. An iterative procedure is used to obtain devolatilization kinetics of coal and biomass for the test-case specific fuels and heating conditions. This is achieved by calibrating the model constants of two empirical models: the single first-order model and the distributed activation energy model. The reference data for calibration are devolatilization yields obtained with predictive coal and biomass multistep kinetic mechanisms. The variation of both particle density and diameter during char combustion is governed by the conversion mode, which is modelled using two approaches: the power law using a constant parameter that assumes a constant mode during char combustion and a constant-free model that considers a variable mode during combustion. Three numerical cases are considered: single first-order reaction with constant char combustion mode, distributed activation energy with constant char combustion mode, and single first-order reaction with variable char combustion mode. The numerical predictions from the large eddy simulations are compared with experimental results of a high co-firing rate large-scale laboratory flame of coal and biomass. Furthermore, results from single particle conversion under idealised conditions, isolating the effects of turbulence, are presented to assist the interpretation of the predictions obtained with large eddy simulations. The effects of the devolatilization and conversion mode modelling on the flame lift-off, flame length, and spatial distribution and radial profiles of O2 and CO2 are presented and discussed. Both the devolatilization and conversion mode modelling have a significant effect on the conversion of particles under idealised conditions. The large eddy simulations results show that the devolatilization model has a strong impact on the flame structure, but not on the flame lift-off. On the other hand, for the tested numerical conditions, the char combustion mode model has a marginal impact on the predicted results.


---

## Body

## 1. Introduction

Co-firing biomass with coal can have a major role in accomplishing the European Union’s renewable targets for power production while contributing to energy security and grid stability [1]. Although biomass co-firing has been successfully demonstrated in over 150 installations worldwide, there are still many open questions concerning flame stability, combustion efficiency, and pollutant emissions [2].There are only a few examples in the literature of detailed experimental studies of large-scale laboratory co-firing under operating conditions that are relevant for power production. Relevant studies include those of Annamalai and coworkers. [3, 4], Damstedt et al. [5, 6], Lu et. al [7], and Casaca and Costa [8]. These studies have shown the positive impact of adding biomass on pollutant emissions, particularly the emissions of NOx. Among these, the work reported by Damstedt and coworkers [5, 6] provides the most complete set of in-flame experimental data for comparison with numerical predictions at a high co-firing ratio (≈50%). Numerical studies of pulverised biomass combustion or co-firing of coal and biomass under realistic conditions are scarce, which is not surprising given the shortage in detailed experimental measurements. Yin et al. [9] compared RANS predictions with the experimental data [5, 6] and obtained a reasonable agreement. They found that modelling intraparticle temperature gradients in large biomass particles had little influence on the predictions. Furthemore, a four-step gas reaction mechanism, combined with eddy disspation rate, outperformed a two-step mechanism, combined with eddy break-up. Other numerical studies are those of Ghenai and Janajreh [10], Ma et al. [11], Black et al. [12], and Blackreedy et al. [13], focused on purely numerical results of RANS simulations of industrial scale biomass combustion and co-firing biomass with coal. The work of Blackreedy et al. [13] focuses on how the devolatilization and the char combustion rate affect the particle burnout; however no comparisons with experimental data are provided. The authors emphasise the importance of the uncertainties in the absolute values of devolatilisation model constants under the heating conditions of the study, given that, for biomass combustion, devolatilisation is the process responsible for the major release of mass to gas phase. Note that the low co-firing ratio in this study is 3%, so that the characteristics of the flame were essentially the same as those of a base coal flame.The process of devolatilization has a significant impact on flame stabilisation and on flame length and consequently on particle burnout and combustion efficiency. A number of advanced coal and biomass conversion models have emerged in the past decade, such as the Chemical Percolation Devolatilization (CPD) model [14] and Bio-PoliMi [15]. The CPD model includes network modelling, coal structure characterization, depolymerization reactions, cross-linking reactions, and noncondensable gas, tar, and char formation. The Bio-PoliMi model has a “mechanist" approach based on conventional multistep devolatilization models of the three reference components of biomass (cellulose, lignin, and hemicellulose) and predicts the yields and lumped composition of the devolatilization products. The complexity of these models makes their use prohibitive in LES; therefore strategies are needed to bridge the advanced conversion models and the cost-effective empirical models that are suitable for CFD codes. Hashimoto et al. [16] used a tabulation strategy based on the particle temperature history to store the results of the devolatilization models for retrieval during the simulation. Vascellari et al. [17] used a number of heating rates from the CFD simulation to calibrate a simple kinetic expression for the volatile yield for later direct use. The latter method has been successfully used in the past for LES of large-scale of coal combustion [18] and Rieth et al. [19] showed that it yields similar averaged predictions as those obtained when directly coupling the CPD model to LES of a small-scale coal flame. The calibration methods have been applied to a single first-order devolatilization model (SFOR) but other empirical models are also good candidates, such as the distributed activation energy model (DAEM) that assumes multiple decomposition steps at a reduced computational cost [20, 21].In contrast to biomass char combustion, coal char combustion has been a subject of research for many years. Still, it can be assumed that the mechanisms of char combustion for coal are also applicable for chars originating from lignocellulosic fuels [22]. Heterogeneous char conversion depends on the char surface area, the surface accessibility, the carbon active sites and catalytic active sites created by indigenous or added inorganic matter, the particle temperature, and the local gaseous reactant concentration. These parameters influence the char combustion mode that in large-scale facilities typically is controlled by both kinetics and diffusion (regime II) or only by diffusion (regime III), depending on the operating conditions. Mitchell et al. [23] showed that the power law mode, based on a tunable model constant, is prone to failure in predicting the variations in size and apparent densities of char particles undergoing oxidation in a moderate zone II burning regime. More recently, Haugen and Mitchel [24] showed that the particle radius decreases with the ratio of the time derivative and the spatial derivative of the particle density at the surface of the particle. A model was proposed that can be used to describe the mode of conversion of reactive porous particles in a range of different applications such as entrained flow gasifiers, pulverized coal burners, and circulating fluidized bed combustors, without any tunable parameters.The main objective of this work is to evaluate the impact of the devolatilization model and the char combustion mode model on the structure of a large-scale, biomass and coal co-fired flame. The devolatization calibration procedure applied previously to PCC [18, 25] is extended to biomass, using SFOR and DAEM to model devolatilization, and the char combustion mode is modelled using two different approaches to describe the variation of density and diameter of the char particle. First, single particle conversion is analysed in order to assess the impact of the tested models on devolatilization and char combustion under idealised conditions and in the absence of turbulence. Subsequently, three LES test cases were simulated in order to assess the effect of devolatilization and conversion mode modelling on the particle conversion in the furnace, the flame structure, and the gas phase major species radial profile. The numerical predictions are compared to the experimental measurements of major gas species of the high co-firing ratio study performed by Damstedt et al. [5, 6] and the effects of the conversion models on the flame structure are discussed.


## 2. Experimental Conditions

Damstedt and coauthors measured the distribution of several gaseous species mole fractions for a set of coal and biomass flames [5, 6]. In this study, theco-firing 1.0 case is targeted. Figure 1 shows the burner and quarl geometry of the Brigham Young Univeristy (BYU) furnace and Table 1 lists the operating conditions. Table 2 lists the coal (high volatile bituminous Blind Canyon) and biomass (straw) properties, while Figure 2 shows the particle size distribution. Note that the particle size distribution of biomass is much wider than that of coal.Table 1
Operating conditions of the targeted test case [5].
Mass flow rate kg/h Center fuel - straw 15.1 Annular fuel - coal 7.5 Central air 9.0 Annular air 12.0 Secondary air 160.0  Swirl Number [-] 1.0Table 2
Properties of coal and straw [5].
Fuel properties Coal Straw Proximate Analysis (wt % ) Moisture (ar) 2.1 7.7 Volatiles (db) 40.6 79.5 Fixed carbon (db) 51.5 15.6 Ash (db) 7.89 4.91 Ultimate Analysis ( % wt, ar) Carbon 74.8 47.3 Hydrogen 5.08 5.68 Oxygen 10.1 41.6 Nitrogen 1.53 0.54 Sulfur 0.58 <0.01 HHV (kJ/kg, db) 29810 17750Figure 1
Burner and quarl geometry of the BYU furnace. Coal is fed through the central inlet, biomass through the annular middle inlet, and secondary swirled air through the outer annular inlet. [5].Figure 2
Particle size distribution represented by cumulative distribution function (CDF) and probability density function (PDF): (a) coal; (b) biomass.
(a) (b)Water-cooled probes were used to sample combustion products along the radius on several axial locations along the furnace. The gas samples were fed to an online Horiba gas analyser to measure major species mole fractions. The measurements of CO2 showed exceptional repeatability of 4.6%, whereas O2 showed a lower repeatability of 20.3%. The uncertainties of the experimental studies, concerning boundary conditions, measurement repeatability, and the experimental flame asymmetry, do not allow for an extensive validation. In this work, the global features of the experimentally studied flame provide valuable comparison data to access the effect of conversion models on the flame shape, namely, flame lift-off, flame break, and flame length.


## 3. Coal and Biomass Combustion Modelling

This work presents an extension of the modelling framework as used in the previous work of Rabaçal et al. [18], where an LES of large-scale laboratory coal flame was presented, by adding the DAEM devolatilization model and the Haugen and Mitchel [24] char conversion mode model. For brevity, the modelling details concerning particle motion, particle heat balance, volatile combustion, and gas-phase radiation are summarised here, whereas the DAEM and the Haugen and Mitchel [24] models are described in more detail.


### 3.1. Particle Motion

The particles’s trajectory is described in a Lagrangian framework. The particles are point sources assumed to be nondeformable and, given that the disperse phase is assumed to be dilute, particle-to-particle interactions are neglected. The only considered forces acting on the particle’s accelerationup are drag FD and gravity FG given by dup/dt=FD+FG. The gravity force is calculated as FG=mp(1-ρg/ρp)g, where mp is the particle mass, ρp is the particle density, ρg is the gas density, and g is the gravitational acceleration. The drag force is calculated as(1)FD=34ρgmpρpdpCDug-upug-up where dp is the particle diameter, ug is the velocity of the gas, and up is the velocity of the particle. The drag coefficient CD has been shown empirically to correlate with the particle Reynolds number Rep=ρgug-updp/μ, where μ is the viscosity of the gas. Given that the particle shape was not reported quantitatively in the experimental study, both coal and biomass particles are assumed spherical. For a spherical particle the correlation is defined by the drag correlation of Yuen and Chen [26]:(2)CD=24Rep1+Rep2/36Rep<10000.424Rep>1000 The particle size is smaller than the largest unresolved scales, which means that the effect of the unresolved scales on the particle trajectory has to be considered. This is done by an additional term added to the particles’ velocity equation as(3)dup=FDdt+FGdt+CokSGSτtdW The last term on the RHS of (3) represents the influence of the subgrid scales on the particle motion, where Co is the dispersion constant and kSGS is the unresolved kinetic energy of the gas phase, which can be obtained from the assumption of equilibrium of the subgrid scales by kSGS=2ΔCs2/3S~ijS~ij. The turbulent particle relaxation time τt is written as(4)τt=τpτpkSGS1/2Δ0.6 where τp is the particle relaxation time, given by(5)τp=43ρpρg1CDug-up Finally, dW represents the incrementalWiener process, which is sampled from a normal distribution with zero mean and standard deviation of dt [27]. The particle position xp evolves according to dxp=updt.


### 3.2. Particle Heat Balance

The timescale of heat conduction inside the particle is typically much higher than that of convection on the surface of the particle; therefore no internal temperature gradients are considered. In these conditions, the particle temperature can be obtained from the balance of heat transfer on the particle surface (convection and radiation) and internal consumption and generation of energy (drying and char combustion). The heat exchange between the particle and the surrounding gas due to convection through the surface of the particle is given by(6)Q˙conv=θNuPrcp,gcp,p1τmTg-Tp where cp,g and cp,p are the specific heat of the gas and the particle, respectively, Tg and Tp are the temperature of the gas and the particle, respectively, and τm=ρpdp2/18μ is the diffusion relaxation time. The Nusselt number is calculated using the Ranz-Marshall correlation, which for spheres is given as(7)Nu=2+0.552Rep1/2Pr1/3In addition, a blowing factor θ is included to account for the reduced heat transfer due to mass transfer between the particle and the gas phase [28]. During the drying phase, the particle will lose heat proportionally to the mass of evaporated water m˙p,H2O and the latent heat of evaporation Levap, calculated from the Watson relation, given by Q˙dry=-m˙p,H2OLevap. During the char combustion phase, heat is generated that contributes to particle heating as Q˙char=m˙p,charLHVchar. Assuming that the char is pure carbon, the heating value takes the value 32.76 MJ/kg. The radiative heat exchange between the particle and its surroundings is given by(8)Q˙rad=Apϵp4πIb,p-Gwhere the particle emissivity is denoted by ϵp, the particle black body intensity is denoted by Ib,p, and the total incident radiation is denoted by G.


### 3.3. Devolatilization

The direct coupling of detailed multistep devolatilization models to LES results in significant costs when simulating large-scale facilities. Empirical models offer cost effectiveness but require model constants which are specific to both the fuel and the heating conditions. In this work, the model constants of empirical models were calibrated to the specific fuel and heating conditions using predictions of detailed multistep devolatilization models for coal and biomass as reference data, following the same method proposed by Vascellari et al. [17], and applied for a large-scale coal flame by Rabaçal et al. [18].The phenomenological network model CPD developed by Grant et al. [14], for coal, and the mechanistic model Bio-PoliMi developed by Ranzi and coworkers [15], for biomass, were used to obtain the reference devolatilization rates and yields. Specific composition parameters are required as an input for both models that can be calculated from the proximate and the ultimate analysis of the fuels. In the case of CPD, the correlation developed by Genetti et al. [29], based on a database of 30 coals, is used to obtain the average coal chemical structural properties: the average molecular weight per side chain, the average molecular weight per aromatic cluster, the ratio of bridges to total attachments, and the total attachments per cluster. In the case of Bio-Polimi, the triangulation method proposed by Cuoci et al. [30] is used to obtain the composition of the biomass in terms of its main components: cellulose, hemicellulose, and lignin.Two empirical formulations have been considered to describe devolatilization in LES. The well established and cost-effective SFOR is used as a reference, where only one global decomposition reaction is considered and it is expressed as(9)dVdt=kV∗-V,withk=A0exp⁡-EaRTp where V is the mass fraction of volatiles released up to time t and V∗ is the total volatile content. The reaction constant k is given by the Arrhenius expression, where A0 is the preexponential factor, Ea is the apparent activation energy, R is the ideal gas constant, and Tp is the particle temperature.In the DAEM formulation, a large number of independent parallel reactions is considered. The contribution of theith reaction is described by a first-order reaction:(10)dVidt=kiVi∗-Vi,withki=A0exp-Ea,iRTp The amount of volatiles released in each reaction is obtained by integration of (10). Assuming that the preexponential factor remains constant and that the number of reactions is large enough to allow the distribution of energy to be expressed as Gaussian distribution with a mean activation energy E0 and a standard deviation σ, the total amount of volatiles released up to a certain time t can be described as(11)V∗-V=V∗σ2π1/2∫0∞exp-A0∫0texp⁡-ERTdt-E-Emean22σ2dE Equation (11) is discretised using the Modified Gauss-Hermite Quadrature method [20] using four quadrature nodes.Figure3 illustrates the calibration process. In the first LES run, SFOR is used considering an initial guess is considered for the model constants. In this work, the preexponential factors and activation energies used by Yin et al. [9] for coal and biomass were used as the initial guess. From the first LES, several particle heating rates from coal and biomass are extracted in order to obtain three representative heating profiles. In the next iteration, the devolatilization rates and yields are calculated with the predictive CPD and PoliMi models by imposing the heating profiles extracted from the previous iteration. The resulting devolatilization yields Yi,jpred, obtained for i=1,Ni heating profiles and j=1,Nj(i) time steps, are the reference data for the current iteration to calibrate the model constants of SFOR and DAEM. The devolatilization yields Yi,jemp are calculated with the empirical models considering a set of model parameters xk (A0 and Ea for SFOR, A0, Emean, and σ for DAEM). The calibration of the model parameters is achieved by minimising the following objective function OF with a genetic algorithm:(12)OFxk=1Ni∑i=1Ni∑j=1NjiYi,jpred-Yi,jempxk2NjimaxjYi,jpred-maxjYi,jemp More details on the calibration method can be found in the work of Vascellari et al. [17]. An LES is subsequently performed using the calibrated model constants obtained in the current iteration and new particle heating rates for the next iteration are obtained. In this study, only 2 iterations steps were required to achieve convergence of the particle heating profile, as in previous studies [17, 18].Figure 3
Flowchart of the calibration process. Adapted from [17].


### 3.4. Char Combustion

The instantaneous burning rate of an individual particle due to char combustion is modelled based on the assumption of a reacting particle surrounded by a chemically frozen boundary layer, the single film assumption. For the present test case, where 80% of the coal particles are under 130μm in diameter, this assumption is adequate. Biomass particles are larger, but given their significantly lower fixed carbon content, the impact of biomass char combustion on the gas phase is much less pronounced than that of coal. The combustion of lignocellulosic chars has been shown to be analogous to coal char combustion [22] and therefore both coal and biomass chars are assumed to be pure carbon. The effects of ashes (catalytic and thermal) are neglected in this work. For both coal and biomass, CO is assumed as the only product of char combustion and the char consumption rate was modelled using the Smith model [31] where both kinetic and diffusion rate limitations are considered:(13)m˙char=πdp2pX~O2RdiffRkinRdiff+Rkinwhere the product pX~O2 is the partial pressure of oxygen in the gas surrounding the particle. The bulk molecular mass diffusion rate is given:(14)Rdiff=CdiffTm0.75dpwhere Tm is the average temperature between the particle and the surrounding gas. The chemical rate depends on both the intrinsic chemical rate and pore diffusion as(15)Rkin=ηdp6ρpAskiwhere η is the ratio of the real combustion rate to the rate achievable if no pore diffusion resistance existed, As is the specific internal surface area, and ki is the intrinsic reactivity rate constant. The latter is calculated as ki=Aiexp⁡(-Ei/RTp), with the values of the intrinsic preexponential factor Ai and intrinsic activation energy Ei based on a line of best fit of several chars standardised at an oxygen partial pressure of 1 atm [31], specifically 305 g/cm2s(atm)1 and 42.8 Kcal/mol, respectively. The effectiveness factor can be related to the Thiele modulus η for a first-order irreversible reaction in a spherical porous particle:(16)η=3ϕ1tanh⁡ϕ-1ϕ The Thiele modulus is defined as(17)ϕ=dp2scharρpAspX~O2kiDeffpX~O2where schar=1.33 is the stoichiometric coefficient assuming that CO is the only product of char combustion. Mitchell et al. [23] have shown that the relation η-ϕ described by (16) can accurately describe the effect of pore diffusion resistance under zone I, moderate zone II, and intense zone II oxidising conditions. Additionally, Deff is the effective overall oxygen pore diffusion coefficient evaluated with the combined effects of bulk and Knudsen diffusion of oxygen through pores as described by the Bosanquet relation:(18)Deff=1DO2+1DKn-1 Following previous work [32], a tree pore model that describes the evolution of the specific area of the char with increasing char conversion is used, which is based on the fractional conversion of carbonaceous material xc:(19)AsAs,0=1-xcxcθc+1-xcDuring char combustion in large-scale applications, char combustion typically proceeds between burning zones II and II. In this work, two different approaches were used to describe density and diameter variation with particle conversion. In the first approach, the classical power law mode is used:(20)ρpρp,0=mCmC,0α(21)dpdp,char=mCmC,0βThe parameters in (20) and (21) are related by the expression α+3β=1 for spherical particles. For char burning under zone II, α takes the constant value of 1/3. Mitchell et al. [23] showed that the power law mode is prone to failure in predicting the variations in size and apparent densities of char particles undergoing oxidation in a moderate zone II burning regime. The second approach used in this work has been introduced very recently by Haugen and Mitchel [24], who derived a parameter-free model. In this model, a characteristic time τ is defined as the time taken for the mass in this outermost shell of the particle to be completely consumed (initial stages of char combustion). Based on the assumption that the reactivity at the outer shell of the particle and the effectiveness factor η are independent, the characteristic time corresponds to the one when conversion equals the effectiveness factor. When the characteristic time has passed, the assumption of a large Thiele modulus (ϕ>4) allows the derivation of the relationship α=η. Finally, the radius and density variation are given by(22)drpdt=0ift<τ,dmCdt1-η4πrp2ρpift≥τ.(23)dρCdt=dmCdt1Vpift<τ,dmCdtηVpift≥τ.In this work, the first approach using the power law and a constantα is named a constant conversion mode model, while the second approach is called a variable conversion mode model. In order to isolate the effects of the conversion mode and exclude any other source of variability in the char combustion rate, the char characteristics (porosity, tortuosity, and initial specific surface area) are kept the same for both fuels based on typical high-volatile bituminous coal values used in the previous study [18].


### 3.5. Gas Phase Combustion

The volatiles from each fuel are represented by a postulate substanceCaHbOc, which respects the proximate and ultimate analysis in a dry ash free (daf) and sulphur and nitrogen free (snf) basis, respectively, and assuming that the fixed carbon consistes of pure carbon, as described by Hashimoto et al. [16]. The molecular weight of the coal volatiles is assumed as 50 kg/kmol [17], whereas for the case of biomass it is assumed as 20 kg/kmol. The low heating value of the volatile matter is obtained from LHVvol=(LHVfuel-YFC,dafLHVFC)/YVM,daf, with LHVFC = 32.76 MJ/kg. Considering the global combustion reaction of the postulate volatiles with air, and assigning the low heating value of the fuel as the heat of reaction, the standard enthalpy of formation of the postulate volatiles can be calculated based on the enthalpy balance of the global combustion reaction.The homogeneous combustion is reduced to a two-step mechanism for both fuels. The reactionsR1 and R2 describe the oxidation of volatiles and carbon monoxide, respectively.R1CaHbOc+a2+b4-c2O2→aCO+b2H2OR2CO+12O2→CO2The turbulence-chemistry interaction on the subgrid level is described using the Eddy Break Up model [33]. The reaction rate is dependent on the inverse mixing time, which is calculated from the resolved strain rate tensor using 1/τ=S~=2S~ijS~ij, as proposed by Hu et al. [34].


### 3.6. Radiation

The general radiative transfer equation is solved using the discrete ordinates method [35–37], being discretized and solved for 24 directions using the S4 approximation. The participating medium’s spectral properties are modelled using a weighted-sum-of-gray-gases model [38]. The particle’s absorption coefficients κp and scattering coefficient σp in a given computational cell are defined according to Chui et al. [37]:(24)κp=ϵp∑i=1Nπdp,i24(25)σp=1-ϵp∑i=1Nπdp,i24 where ϵp depends on the char-burnout and the proportions of volatile content and ash found in the particle [37]. Equations (24) and (25) are valid for large, grey, and diffuse particles. Only particle scattering is considered as this is largely dominant over gas, soot, and fly-ash scattering. Furthermore, throughout this work only isotropic scattering is considered. For the experimental test cases used in this work, considering the measured and predicted temperatures and the smallest particle diameter (15 μm), it can be estimated that only negligible amounts of emitted radiative energy are on the wavelengths where the large particle assumption fails.


## 3.1. Particle Motion

The particles’s trajectory is described in a Lagrangian framework. The particles are point sources assumed to be nondeformable and, given that the disperse phase is assumed to be dilute, particle-to-particle interactions are neglected. The only considered forces acting on the particle’s accelerationup are drag FD and gravity FG given by dup/dt=FD+FG. The gravity force is calculated as FG=mp(1-ρg/ρp)g, where mp is the particle mass, ρp is the particle density, ρg is the gas density, and g is the gravitational acceleration. The drag force is calculated as(1)FD=34ρgmpρpdpCDug-upug-up where dp is the particle diameter, ug is the velocity of the gas, and up is the velocity of the particle. The drag coefficient CD has been shown empirically to correlate with the particle Reynolds number Rep=ρgug-updp/μ, where μ is the viscosity of the gas. Given that the particle shape was not reported quantitatively in the experimental study, both coal and biomass particles are assumed spherical. For a spherical particle the correlation is defined by the drag correlation of Yuen and Chen [26]:(2)CD=24Rep1+Rep2/36Rep<10000.424Rep>1000 The particle size is smaller than the largest unresolved scales, which means that the effect of the unresolved scales on the particle trajectory has to be considered. This is done by an additional term added to the particles’ velocity equation as(3)dup=FDdt+FGdt+CokSGSτtdW The last term on the RHS of (3) represents the influence of the subgrid scales on the particle motion, where Co is the dispersion constant and kSGS is the unresolved kinetic energy of the gas phase, which can be obtained from the assumption of equilibrium of the subgrid scales by kSGS=2ΔCs2/3S~ijS~ij. The turbulent particle relaxation time τt is written as(4)τt=τpτpkSGS1/2Δ0.6 where τp is the particle relaxation time, given by(5)τp=43ρpρg1CDug-up Finally, dW represents the incrementalWiener process, which is sampled from a normal distribution with zero mean and standard deviation of dt [27]. The particle position xp evolves according to dxp=updt.


## 3.2. Particle Heat Balance

The timescale of heat conduction inside the particle is typically much higher than that of convection on the surface of the particle; therefore no internal temperature gradients are considered. In these conditions, the particle temperature can be obtained from the balance of heat transfer on the particle surface (convection and radiation) and internal consumption and generation of energy (drying and char combustion). The heat exchange between the particle and the surrounding gas due to convection through the surface of the particle is given by(6)Q˙conv=θNuPrcp,gcp,p1τmTg-Tp where cp,g and cp,p are the specific heat of the gas and the particle, respectively, Tg and Tp are the temperature of the gas and the particle, respectively, and τm=ρpdp2/18μ is the diffusion relaxation time. The Nusselt number is calculated using the Ranz-Marshall correlation, which for spheres is given as(7)Nu=2+0.552Rep1/2Pr1/3In addition, a blowing factor θ is included to account for the reduced heat transfer due to mass transfer between the particle and the gas phase [28]. During the drying phase, the particle will lose heat proportionally to the mass of evaporated water m˙p,H2O and the latent heat of evaporation Levap, calculated from the Watson relation, given by Q˙dry=-m˙p,H2OLevap. During the char combustion phase, heat is generated that contributes to particle heating as Q˙char=m˙p,charLHVchar. Assuming that the char is pure carbon, the heating value takes the value 32.76 MJ/kg. The radiative heat exchange between the particle and its surroundings is given by(8)Q˙rad=Apϵp4πIb,p-Gwhere the particle emissivity is denoted by ϵp, the particle black body intensity is denoted by Ib,p, and the total incident radiation is denoted by G.


## 3.3. Devolatilization

The direct coupling of detailed multistep devolatilization models to LES results in significant costs when simulating large-scale facilities. Empirical models offer cost effectiveness but require model constants which are specific to both the fuel and the heating conditions. In this work, the model constants of empirical models were calibrated to the specific fuel and heating conditions using predictions of detailed multistep devolatilization models for coal and biomass as reference data, following the same method proposed by Vascellari et al. [17], and applied for a large-scale coal flame by Rabaçal et al. [18].The phenomenological network model CPD developed by Grant et al. [14], for coal, and the mechanistic model Bio-PoliMi developed by Ranzi and coworkers [15], for biomass, were used to obtain the reference devolatilization rates and yields. Specific composition parameters are required as an input for both models that can be calculated from the proximate and the ultimate analysis of the fuels. In the case of CPD, the correlation developed by Genetti et al. [29], based on a database of 30 coals, is used to obtain the average coal chemical structural properties: the average molecular weight per side chain, the average molecular weight per aromatic cluster, the ratio of bridges to total attachments, and the total attachments per cluster. In the case of Bio-Polimi, the triangulation method proposed by Cuoci et al. [30] is used to obtain the composition of the biomass in terms of its main components: cellulose, hemicellulose, and lignin.Two empirical formulations have been considered to describe devolatilization in LES. The well established and cost-effective SFOR is used as a reference, where only one global decomposition reaction is considered and it is expressed as(9)dVdt=kV∗-V,withk=A0exp⁡-EaRTp where V is the mass fraction of volatiles released up to time t and V∗ is the total volatile content. The reaction constant k is given by the Arrhenius expression, where A0 is the preexponential factor, Ea is the apparent activation energy, R is the ideal gas constant, and Tp is the particle temperature.In the DAEM formulation, a large number of independent parallel reactions is considered. The contribution of theith reaction is described by a first-order reaction:(10)dVidt=kiVi∗-Vi,withki=A0exp-Ea,iRTp The amount of volatiles released in each reaction is obtained by integration of (10). Assuming that the preexponential factor remains constant and that the number of reactions is large enough to allow the distribution of energy to be expressed as Gaussian distribution with a mean activation energy E0 and a standard deviation σ, the total amount of volatiles released up to a certain time t can be described as(11)V∗-V=V∗σ2π1/2∫0∞exp-A0∫0texp⁡-ERTdt-E-Emean22σ2dE Equation (11) is discretised using the Modified Gauss-Hermite Quadrature method [20] using four quadrature nodes.Figure3 illustrates the calibration process. In the first LES run, SFOR is used considering an initial guess is considered for the model constants. In this work, the preexponential factors and activation energies used by Yin et al. [9] for coal and biomass were used as the initial guess. From the first LES, several particle heating rates from coal and biomass are extracted in order to obtain three representative heating profiles. In the next iteration, the devolatilization rates and yields are calculated with the predictive CPD and PoliMi models by imposing the heating profiles extracted from the previous iteration. The resulting devolatilization yields Yi,jpred, obtained for i=1,Ni heating profiles and j=1,Nj(i) time steps, are the reference data for the current iteration to calibrate the model constants of SFOR and DAEM. The devolatilization yields Yi,jemp are calculated with the empirical models considering a set of model parameters xk (A0 and Ea for SFOR, A0, Emean, and σ for DAEM). The calibration of the model parameters is achieved by minimising the following objective function OF with a genetic algorithm:(12)OFxk=1Ni∑i=1Ni∑j=1NjiYi,jpred-Yi,jempxk2NjimaxjYi,jpred-maxjYi,jemp More details on the calibration method can be found in the work of Vascellari et al. [17]. An LES is subsequently performed using the calibrated model constants obtained in the current iteration and new particle heating rates for the next iteration are obtained. In this study, only 2 iterations steps were required to achieve convergence of the particle heating profile, as in previous studies [17, 18].Figure 3
Flowchart of the calibration process. Adapted from [17].


## 3.4. Char Combustion

The instantaneous burning rate of an individual particle due to char combustion is modelled based on the assumption of a reacting particle surrounded by a chemically frozen boundary layer, the single film assumption. For the present test case, where 80% of the coal particles are under 130μm in diameter, this assumption is adequate. Biomass particles are larger, but given their significantly lower fixed carbon content, the impact of biomass char combustion on the gas phase is much less pronounced than that of coal. The combustion of lignocellulosic chars has been shown to be analogous to coal char combustion [22] and therefore both coal and biomass chars are assumed to be pure carbon. The effects of ashes (catalytic and thermal) are neglected in this work. For both coal and biomass, CO is assumed as the only product of char combustion and the char consumption rate was modelled using the Smith model [31] where both kinetic and diffusion rate limitations are considered:(13)m˙char=πdp2pX~O2RdiffRkinRdiff+Rkinwhere the product pX~O2 is the partial pressure of oxygen in the gas surrounding the particle. The bulk molecular mass diffusion rate is given:(14)Rdiff=CdiffTm0.75dpwhere Tm is the average temperature between the particle and the surrounding gas. The chemical rate depends on both the intrinsic chemical rate and pore diffusion as(15)Rkin=ηdp6ρpAskiwhere η is the ratio of the real combustion rate to the rate achievable if no pore diffusion resistance existed, As is the specific internal surface area, and ki is the intrinsic reactivity rate constant. The latter is calculated as ki=Aiexp⁡(-Ei/RTp), with the values of the intrinsic preexponential factor Ai and intrinsic activation energy Ei based on a line of best fit of several chars standardised at an oxygen partial pressure of 1 atm [31], specifically 305 g/cm2s(atm)1 and 42.8 Kcal/mol, respectively. The effectiveness factor can be related to the Thiele modulus η for a first-order irreversible reaction in a spherical porous particle:(16)η=3ϕ1tanh⁡ϕ-1ϕ The Thiele modulus is defined as(17)ϕ=dp2scharρpAspX~O2kiDeffpX~O2where schar=1.33 is the stoichiometric coefficient assuming that CO is the only product of char combustion. Mitchell et al. [23] have shown that the relation η-ϕ described by (16) can accurately describe the effect of pore diffusion resistance under zone I, moderate zone II, and intense zone II oxidising conditions. Additionally, Deff is the effective overall oxygen pore diffusion coefficient evaluated with the combined effects of bulk and Knudsen diffusion of oxygen through pores as described by the Bosanquet relation:(18)Deff=1DO2+1DKn-1 Following previous work [32], a tree pore model that describes the evolution of the specific area of the char with increasing char conversion is used, which is based on the fractional conversion of carbonaceous material xc:(19)AsAs,0=1-xcxcθc+1-xcDuring char combustion in large-scale applications, char combustion typically proceeds between burning zones II and II. In this work, two different approaches were used to describe density and diameter variation with particle conversion. In the first approach, the classical power law mode is used:(20)ρpρp,0=mCmC,0α(21)dpdp,char=mCmC,0βThe parameters in (20) and (21) are related by the expression α+3β=1 for spherical particles. For char burning under zone II, α takes the constant value of 1/3. Mitchell et al. [23] showed that the power law mode is prone to failure in predicting the variations in size and apparent densities of char particles undergoing oxidation in a moderate zone II burning regime. The second approach used in this work has been introduced very recently by Haugen and Mitchel [24], who derived a parameter-free model. In this model, a characteristic time τ is defined as the time taken for the mass in this outermost shell of the particle to be completely consumed (initial stages of char combustion). Based on the assumption that the reactivity at the outer shell of the particle and the effectiveness factor η are independent, the characteristic time corresponds to the one when conversion equals the effectiveness factor. When the characteristic time has passed, the assumption of a large Thiele modulus (ϕ>4) allows the derivation of the relationship α=η. Finally, the radius and density variation are given by(22)drpdt=0ift<τ,dmCdt1-η4πrp2ρpift≥τ.(23)dρCdt=dmCdt1Vpift<τ,dmCdtηVpift≥τ.In this work, the first approach using the power law and a constantα is named a constant conversion mode model, while the second approach is called a variable conversion mode model. In order to isolate the effects of the conversion mode and exclude any other source of variability in the char combustion rate, the char characteristics (porosity, tortuosity, and initial specific surface area) are kept the same for both fuels based on typical high-volatile bituminous coal values used in the previous study [18].


## 3.5. Gas Phase Combustion

The volatiles from each fuel are represented by a postulate substanceCaHbOc, which respects the proximate and ultimate analysis in a dry ash free (daf) and sulphur and nitrogen free (snf) basis, respectively, and assuming that the fixed carbon consistes of pure carbon, as described by Hashimoto et al. [16]. The molecular weight of the coal volatiles is assumed as 50 kg/kmol [17], whereas for the case of biomass it is assumed as 20 kg/kmol. The low heating value of the volatile matter is obtained from LHVvol=(LHVfuel-YFC,dafLHVFC)/YVM,daf, with LHVFC = 32.76 MJ/kg. Considering the global combustion reaction of the postulate volatiles with air, and assigning the low heating value of the fuel as the heat of reaction, the standard enthalpy of formation of the postulate volatiles can be calculated based on the enthalpy balance of the global combustion reaction.The homogeneous combustion is reduced to a two-step mechanism for both fuels. The reactionsR1 and R2 describe the oxidation of volatiles and carbon monoxide, respectively.R1CaHbOc+a2+b4-c2O2→aCO+b2H2OR2CO+12O2→CO2The turbulence-chemistry interaction on the subgrid level is described using the Eddy Break Up model [33]. The reaction rate is dependent on the inverse mixing time, which is calculated from the resolved strain rate tensor using 1/τ=S~=2S~ijS~ij, as proposed by Hu et al. [34].


## 3.6. Radiation

The general radiative transfer equation is solved using the discrete ordinates method [35–37], being discretized and solved for 24 directions using the S4 approximation. The participating medium’s spectral properties are modelled using a weighted-sum-of-gray-gases model [38]. The particle’s absorption coefficients κp and scattering coefficient σp in a given computational cell are defined according to Chui et al. [37]:(24)κp=ϵp∑i=1Nπdp,i24(25)σp=1-ϵp∑i=1Nπdp,i24 where ϵp depends on the char-burnout and the proportions of volatile content and ash found in the particle [37]. Equations (24) and (25) are valid for large, grey, and diffuse particles. Only particle scattering is considered as this is largely dominant over gas, soot, and fly-ash scattering. Furthermore, throughout this work only isotropic scattering is considered. For the experimental test cases used in this work, considering the measured and predicted temperatures and the smallest particle diameter (15 μm), it can be estimated that only negligible amounts of emitted radiative energy are on the wavelengths where the large particle assumption fails.


## 4. Numerical Setup

The in-house code PsiPhi [18, 25, 39–43] was used to solve the implicitly filtered, low-Mach number Navier Stokes equations for a variable density flow. PsiPhi ensures continuity by a pressure-correction scheme using a projection method. A Cartesian, equidistant grid is used, providing good numerical accuracy, isotropic filtering, high parallel efficiency, and good vectorisation through avoiding slow, nonsequential memory access. The configuration allows decomposing the domain into blocks for maximum efficiency during MPI communication between CPUs.The convective fluxes are approximated with a 2nd-order central differencing scheme for momentum and a total variation diminishing scheme for scalars. A low storage, explicit Runge Kutta scheme is used for time-integration with a CFL number of 0.3 to ensure accuracy and stability. The particle motion is embedded into the Runge-Kutta procedure of the LES code and the coupling between particle and LES fields is done by trilinear interpolation schemes.In this work, a particle parceling strategy was applied to achieve a compromise between accuracy and computational cost. Given the biomass mass flow input and the biomass particle size distribution (c.f. Figure2), approximately ten particles are injected per time step. As for coal particles, approximately 200 particles are injected per time step. While the order of magnitude of the number of physical biomass particles injected per time step is affordable, in the case of physical coal particles that number is too large. Therefore ten physical coal particles are parcelled into one numerical particle. The particles were injected in the domain using a random distribution fitted to the particle size distribution and randomly positioned within the injection plane.The closure of the subgrid stress tensor is ensured through the eddy viscosity approach, where the turbulent viscosity is modelled using the Smagorinsky-Lilly model [44]. The model constant takes the value 0.173. Pseudo-turbulent inflow conditions were generated using an efficient implementation of Klein’s turbulence generator [45], for a turbulence level of 10% for the fluctuation magnitude. Inlets mean velocities were assigned using a 7th-power law profile matching the mass flow of each stream.Immersed boundaries describe the furnace walls and no additional wall modelling was attempted as the flame burns away from the walls and as the wall roughness is hardly known. A wall temperature profile was imposed to the immersed boundaries (D. Tree, personal communication). The furnace geometry was described by an equidistant Cartesian mesh with 2250x750x750 mm in size using a cell size ofΔ = 3.0 mm. This resolution has been shown to provide reasonable predictions of a coal flame in a furnace and burner of similar dimensions, operating under similar conditions [18]. The grid has a total of 40 million cells and 1024 cores were used in SuperMUC (LRZ, Munich, Germany). Statistics were collected after 5 seconds of initialisation over a period of 1 second. The numerical cost of each run was 100 thousand CPUh and 3 runs were necessary for the iterative optimisation process, corresponding to a total of 0.5 million CPUh. Three numerical setups were simulated, listed in Table 3, varying the devolatilization model and the char combustion mode model.Table 3
Numerical test case setup list.
Test case Devolatilization Conversion mode SC SFOR Constant DC DAEM Constant SV SFOR Variable


## 5. Results

### 5.1. Single Particle Conversion under Idealised Conditions

Figures4 and 5 show the total yield and rate curves of devolatilization for coal and biomass, respectively. The optimised constants are given in Table 4, where the Q factor refers to the ratio between the total volatile yield predicted by CPD and PoliMi under the heating conditions retrieved from the LES and the total volatile yield reported by the proximate analysis. Note that the latter is obtained following measuring standards where the heating rates and maximum temperatures are far from those in typical large-scale combustion (e.g., ASTM D3172 and ASTM E872). Note that coal and particles are exposed to different heating rates due to the different particle size distributions (see Figure 2). Both coal and biomass particles of similar size are subjected to heating rates in the order of 1×104 K/s, whereas larger biomass particles are additionally subjected to heating rates in the order of 1×103 K/s given their broader size distribution. Coal devolatilization presents several rate peaks whereas biomass presents one single peak at such high heating rates. In the case of coal, the DAEM curve presents better agreement with the CPD curve than that of SFOR. In the case of biomass, both calibrated curves present a good agreement with the PoliMi curve.Table 4
Optimised model constants of the empirical devolatilization models.
Coal Biomass  Coal Biomass SFOR   DAEM A0 1.95×105 1.19×108 A 0 12.29 8.73 E 7.775×103 13.418×103 E m e a n 192.47 125.09 Q 0.5558 0.8475 σ 26.75 4.73 Q 0.5598 0.8476Figure 4
Coal devolatilization. Top: total yield (daf); down: total rate (daf) normalised by the maximum value of the CPD rate. (a) From 300 to 1700 K in 0.06 s, (b) from 300 to 1600 K in 0.1 s, and (c) from 300 to 1500 K in 0.15 s.Figure 5
Biomass devolatilization. Top: total yield (daf); down: total rate (daf) normalised by the maximum value of the CPD rate. (a) From 300 to 1600 K in 0.04 s, (b) from 300 to 1500 K in 0.15 s, and (c) from 300 to 1400 K in 0.25 s.Figure6 shows the effect of the conversion mode modelling approach on the char conversion, normalised density, and normalised diameter variation during char conversion of a single particle. A char particle with 100 μm in diameter was immersed in an atmosphere with 6% O2 and at two temperatures, 1200 K and 1500K, corresponding to a moderate and an intensive zone II burning mode, respectively.Figure 6
Time evolution of char conversion conversion (xc), normalised diameter (dp∗=dp/dp,0), and normalised density (ρp∗=ρp/ρp,0) of a single particle char immersed in an atmosphere with 6%  O2 at 1200K. The subscripts C and V denote constant conversion mode and variable conversion mode, respectively.The conversion mode model affects the char diameter and density but has little influence on the burnout time. The trends obtained with the variable conversion mode are in good agreement with the ones obtained by Haugen and Mitchell [24] for similar conditions. At the lower temperature of 1200 K, both models predict a similar variation of the diameter, although a small deviation can be observed in the latter stages of conversion. The density variation curves present a more pronounced deviation throughout conversion. At the higher temperature of 1500 K, the diameter and density variation show clear deviations between conversion models. Under an intensive zone II burning mode (1500 K), where it is expected that the density variation is close to that of a shrinking core mode, the variable conversion mode yields reasonable results.


### 5.2. Particle Conversion in the Furnace Simulations

Figure7 shows the LES instantaneous devolatilization rate over time of a particle ensemble collected at 6 seconds of physical run time for the cases SC and DC, cf. Table 3. When SFOM is used for calibration (Figures 7(a) and 7(c)), corresponding to case SC, the devolatilization of each fuel type occurs mainly within a narrow location. When DAEM is used for calibration (Figures 7(b) and 7(d)), corresponding to case DC, devolatilization of both fuels spans through a longer axial distance and maximum rates decrease slightly. Since the local distribution of volatiles in the furnace affects the local gas phase mixture composition, these results provide a first indicator of the effect of the devolatilization model on the gaseous flame length, for the studied conditions.Figure 7
Devolatilization rate along the axial direction of the furnace: (a) SC coal, (b) DC coal, (c) SC biomass, and (d) DC biomass. Rates normalised by maximum rate obtained with SFOM of corresponding fuel. Particle ensemble collected after 6 seconds of physical runtime.Figure8 shows the LES instantaneous normalised particle diameter, density, temperature, and velocity magnitude during char conversion of a particle ensemble collected at 6 seconds of physical run time for the cases SC and SV (cf. Table 3). Particle temperature is normalised using the maximum temperature obtained in case SC and the analogous normalisation is applied to the velocity magnitude. Figures 8(c) and 8(d) show that the conversion mode model has an influence on the particles aerodynamic properties, as typified by the particle velocity variation during conversion. As a consequence, local char burning conditions will vary yielding the differences observed in the particle temperature during char combustion, as seen in Figures 8(a) and 8(b).Figure 8
Normalised diameter, density, velocity magnitude (up∗=up/up,max), and particle temperature (Tp∗=Tp/Tp,max) variation during char conversion as a function of particle temperature: (a) and (c) SC, (b) and (d) SV. Particle ensemble collected after 6 seconds of physical runtime.


### 5.3. Effects of the Particle Conversion Model on the Flame Structure

Figure9 shows the instantaneous radial mid-plane gas phase temperature distribution corresponding to 2.5 seconds of physical runtime, sufficient for the establishment of the gaseous flame. The temperature is used as an indicator of regions where intense heat release occurs and hence where combustion reactions are taking place. Figures 9(a) and 9(c) show that when using SFOR the reactions occur in a well defined envelope enclosing the internal recirculation zone and further downstream pockets of cold unreacted gas. On the other hand, Figure 9(b) shows that when using DAEM the high temperature regions are broader and spacial gradients are less intense downstream of the internal recirculation zone. Finally, Figure 9 shows that the flame lift-off height is marginally affected by the devolatilization model in the studied conditions, as typified by the early onset of reactions within the quarl.Figure 9
Instantaneous radial mid-plane gas phase temperature distribution corresponding to 2.5 seconds of physical runtime, sufficient for the establishment of the gaseous flame. Note that only the region of interest is shown, corresponding to a third of the diameter and half of the length. From left to right: SC, DC, and SV.
(a) (b) (c)Figures10 and 11 show experimental and numerical maps of CO2 and O2 mean concentration, respectively. The experimental maps show that the flame has a large volume, extending axially beyond 125 cm, and expanding outward radially for 15-20 cm. This large flame volume can be explained by several factors, namely, [5]: (1) given the low energy density in mass base of biomass, considerable mass flows of transport air need to be used to supply the furnace with a certain thermal input; therefore particles will be fed with a high momentum that allows the particles to penetrate through the internal recirculation zone; (2) large biomass particles, as compared to coal, typically show higher volatile yields and longer devolatilization times, spreading the reaction of the biomass particles over a larger volume; (3) the high volatile yields of biomass create local fuel-rich volumes in a section of the furnace where mixing is more limited due to the low gas velocities, which limits the homogenous gas-phase reactions. Another interesting aspect of the structure of flame is the “break” occurring between 20 and 35 cm. Damsteadt [5] attributes the occurrence of the “break” to the characteristics of the biomass used in the study. Straw contains two main structures: stalks and knees. The knees are much harder and less porous than the stalks. Upon grinding, the shape of the knees is near spherical, while the stalks are either flake-like or cylindrical. The knees proceed mostly unreacted through the internal recirculation zone into the lower regions of the reactor, where sufficient time in the hot reactor environment has elapsed to begin the combustion process. When comparing the different devolatilization models, specifically the cases SC and DC, it is possible to see that the flame break is only captured when DAEM is used. Furthermore, it is possible to observe in both figures the effect of large biomass particles that are marginally affected by the internal recirculation zone, penetrating through it and travelling downwards. The results indicate that SFOR overpredicts the devolatilization rates close to the burner yielding a faster burnout of coal and biomass particles. In parallel, gas-phase reactions are also more intense resulting in a larger radial expansion of the internal recirculation zone. In contrast, the case DC shows a radially compact internal recirculation zone where smaller particles are entrained and burn, and a radially compact secondary reaction zone after the flame break, owing to a more distributed volatile release as shown in Figure 7. Note that Figure 7 shows that the location in the furnace where the maximum devolatilization rate occurs varies depending on the choice of the devolatilization model. Essentially, it shows that with SFOR the predicted devolatilization rates occur very close to the burner exit and the devolatilization process is over at 400 mm. Experimentally, the flame break is observed at approximately 250 mm, downstream of which large particles go through devolatilization creating the second flame. This is captured by the DAEM in both Figure 7 and Figures 10 and 11. The O2 and CO2 maps obtained with the cases SC and SV show that under these conditions the average spatial distribution of the major species is not affected by the choice of the char combustion mode model, which is a result of high intensity reaction zone in the upper region of the furnace leading to early burnout under high temperatures.Figure 10
Experimental (adapted from [5]) and numerical maps of CO2 concentration. Dimensions on the the left are given in centimetres.Figure 11
Experimental (adapted from [5]) and numerical maps of O2 concentration. Dimensions on the the left are given in centimetres.Figure12 shows a comparison between the experimental measurements and the numerical predictions of O2 and CO2. The flame showed considerable asymmetry; therefore absolute comparisons between measurements and predictions should be regarded as qualitative. In order to minimise the asymmetry, points and corresponding error bars show the mean and standard deviation of measurements performed at opposing, radially symmetric locations. The predictions from the cases SC and DC show that the devolatilization model has a pronounced effect in the flame structure, whereas the char conversion mode model does not show significant differences. When SFOR is used, cases SC and SV, it is clear that the overall particle heating and reaction rates are overpredicted, given the underprediction of O2 and the overprediction of CO2. Furthermore, the particle burnout at the exit of the domain in those cases was close to 100%, but qualitative analysis of the experimentally collected char at the exit of the furnace showed that biomass, particularly the knees, showed limited burnout levels. When using DAEM, the predictions are closer to the experimental results and burnout levels at the exit of the computational domain are in the order of 60%.Figure 12
Comparison between measurements and predictions of O2 (left) and CO2 (right) for the modelling cases SCc, DCc, and SVc. Empty points and corresponding error bars show the mean and standard deviation of measurements performed at opposing, radially symmetric locations.


## 5.1. Single Particle Conversion under Idealised Conditions

Figures4 and 5 show the total yield and rate curves of devolatilization for coal and biomass, respectively. The optimised constants are given in Table 4, where the Q factor refers to the ratio between the total volatile yield predicted by CPD and PoliMi under the heating conditions retrieved from the LES and the total volatile yield reported by the proximate analysis. Note that the latter is obtained following measuring standards where the heating rates and maximum temperatures are far from those in typical large-scale combustion (e.g., ASTM D3172 and ASTM E872). Note that coal and particles are exposed to different heating rates due to the different particle size distributions (see Figure 2). Both coal and biomass particles of similar size are subjected to heating rates in the order of 1×104 K/s, whereas larger biomass particles are additionally subjected to heating rates in the order of 1×103 K/s given their broader size distribution. Coal devolatilization presents several rate peaks whereas biomass presents one single peak at such high heating rates. In the case of coal, the DAEM curve presents better agreement with the CPD curve than that of SFOR. In the case of biomass, both calibrated curves present a good agreement with the PoliMi curve.Table 4
Optimised model constants of the empirical devolatilization models.
Coal Biomass  Coal Biomass SFOR   DAEM A0 1.95×105 1.19×108 A 0 12.29 8.73 E 7.775×103 13.418×103 E m e a n 192.47 125.09 Q 0.5558 0.8475 σ 26.75 4.73 Q 0.5598 0.8476Figure 4
Coal devolatilization. Top: total yield (daf); down: total rate (daf) normalised by the maximum value of the CPD rate. (a) From 300 to 1700 K in 0.06 s, (b) from 300 to 1600 K in 0.1 s, and (c) from 300 to 1500 K in 0.15 s.Figure 5
Biomass devolatilization. Top: total yield (daf); down: total rate (daf) normalised by the maximum value of the CPD rate. (a) From 300 to 1600 K in 0.04 s, (b) from 300 to 1500 K in 0.15 s, and (c) from 300 to 1400 K in 0.25 s.Figure6 shows the effect of the conversion mode modelling approach on the char conversion, normalised density, and normalised diameter variation during char conversion of a single particle. A char particle with 100 μm in diameter was immersed in an atmosphere with 6% O2 and at two temperatures, 1200 K and 1500K, corresponding to a moderate and an intensive zone II burning mode, respectively.Figure 6
Time evolution of char conversion conversion (xc), normalised diameter (dp∗=dp/dp,0), and normalised density (ρp∗=ρp/ρp,0) of a single particle char immersed in an atmosphere with 6%  O2 at 1200K. The subscripts C and V denote constant conversion mode and variable conversion mode, respectively.The conversion mode model affects the char diameter and density but has little influence on the burnout time. The trends obtained with the variable conversion mode are in good agreement with the ones obtained by Haugen and Mitchell [24] for similar conditions. At the lower temperature of 1200 K, both models predict a similar variation of the diameter, although a small deviation can be observed in the latter stages of conversion. The density variation curves present a more pronounced deviation throughout conversion. At the higher temperature of 1500 K, the diameter and density variation show clear deviations between conversion models. Under an intensive zone II burning mode (1500 K), where it is expected that the density variation is close to that of a shrinking core mode, the variable conversion mode yields reasonable results.


## 5.2. Particle Conversion in the Furnace Simulations

Figure7 shows the LES instantaneous devolatilization rate over time of a particle ensemble collected at 6 seconds of physical run time for the cases SC and DC, cf. Table 3. When SFOM is used for calibration (Figures 7(a) and 7(c)), corresponding to case SC, the devolatilization of each fuel type occurs mainly within a narrow location. When DAEM is used for calibration (Figures 7(b) and 7(d)), corresponding to case DC, devolatilization of both fuels spans through a longer axial distance and maximum rates decrease slightly. Since the local distribution of volatiles in the furnace affects the local gas phase mixture composition, these results provide a first indicator of the effect of the devolatilization model on the gaseous flame length, for the studied conditions.Figure 7
Devolatilization rate along the axial direction of the furnace: (a) SC coal, (b) DC coal, (c) SC biomass, and (d) DC biomass. Rates normalised by maximum rate obtained with SFOM of corresponding fuel. Particle ensemble collected after 6 seconds of physical runtime.Figure8 shows the LES instantaneous normalised particle diameter, density, temperature, and velocity magnitude during char conversion of a particle ensemble collected at 6 seconds of physical run time for the cases SC and SV (cf. Table 3). Particle temperature is normalised using the maximum temperature obtained in case SC and the analogous normalisation is applied to the velocity magnitude. Figures 8(c) and 8(d) show that the conversion mode model has an influence on the particles aerodynamic properties, as typified by the particle velocity variation during conversion. As a consequence, local char burning conditions will vary yielding the differences observed in the particle temperature during char combustion, as seen in Figures 8(a) and 8(b).Figure 8
Normalised diameter, density, velocity magnitude (up∗=up/up,max), and particle temperature (Tp∗=Tp/Tp,max) variation during char conversion as a function of particle temperature: (a) and (c) SC, (b) and (d) SV. Particle ensemble collected after 6 seconds of physical runtime.


## 5.3. Effects of the Particle Conversion Model on the Flame Structure

Figure9 shows the instantaneous radial mid-plane gas phase temperature distribution corresponding to 2.5 seconds of physical runtime, sufficient for the establishment of the gaseous flame. The temperature is used as an indicator of regions where intense heat release occurs and hence where combustion reactions are taking place. Figures 9(a) and 9(c) show that when using SFOR the reactions occur in a well defined envelope enclosing the internal recirculation zone and further downstream pockets of cold unreacted gas. On the other hand, Figure 9(b) shows that when using DAEM the high temperature regions are broader and spacial gradients are less intense downstream of the internal recirculation zone. Finally, Figure 9 shows that the flame lift-off height is marginally affected by the devolatilization model in the studied conditions, as typified by the early onset of reactions within the quarl.Figure 9
Instantaneous radial mid-plane gas phase temperature distribution corresponding to 2.5 seconds of physical runtime, sufficient for the establishment of the gaseous flame. Note that only the region of interest is shown, corresponding to a third of the diameter and half of the length. From left to right: SC, DC, and SV.
(a) (b) (c)Figures10 and 11 show experimental and numerical maps of CO2 and O2 mean concentration, respectively. The experimental maps show that the flame has a large volume, extending axially beyond 125 cm, and expanding outward radially for 15-20 cm. This large flame volume can be explained by several factors, namely, [5]: (1) given the low energy density in mass base of biomass, considerable mass flows of transport air need to be used to supply the furnace with a certain thermal input; therefore particles will be fed with a high momentum that allows the particles to penetrate through the internal recirculation zone; (2) large biomass particles, as compared to coal, typically show higher volatile yields and longer devolatilization times, spreading the reaction of the biomass particles over a larger volume; (3) the high volatile yields of biomass create local fuel-rich volumes in a section of the furnace where mixing is more limited due to the low gas velocities, which limits the homogenous gas-phase reactions. Another interesting aspect of the structure of flame is the “break” occurring between 20 and 35 cm. Damsteadt [5] attributes the occurrence of the “break” to the characteristics of the biomass used in the study. Straw contains two main structures: stalks and knees. The knees are much harder and less porous than the stalks. Upon grinding, the shape of the knees is near spherical, while the stalks are either flake-like or cylindrical. The knees proceed mostly unreacted through the internal recirculation zone into the lower regions of the reactor, where sufficient time in the hot reactor environment has elapsed to begin the combustion process. When comparing the different devolatilization models, specifically the cases SC and DC, it is possible to see that the flame break is only captured when DAEM is used. Furthermore, it is possible to observe in both figures the effect of large biomass particles that are marginally affected by the internal recirculation zone, penetrating through it and travelling downwards. The results indicate that SFOR overpredicts the devolatilization rates close to the burner yielding a faster burnout of coal and biomass particles. In parallel, gas-phase reactions are also more intense resulting in a larger radial expansion of the internal recirculation zone. In contrast, the case DC shows a radially compact internal recirculation zone where smaller particles are entrained and burn, and a radially compact secondary reaction zone after the flame break, owing to a more distributed volatile release as shown in Figure 7. Note that Figure 7 shows that the location in the furnace where the maximum devolatilization rate occurs varies depending on the choice of the devolatilization model. Essentially, it shows that with SFOR the predicted devolatilization rates occur very close to the burner exit and the devolatilization process is over at 400 mm. Experimentally, the flame break is observed at approximately 250 mm, downstream of which large particles go through devolatilization creating the second flame. This is captured by the DAEM in both Figure 7 and Figures 10 and 11. The O2 and CO2 maps obtained with the cases SC and SV show that under these conditions the average spatial distribution of the major species is not affected by the choice of the char combustion mode model, which is a result of high intensity reaction zone in the upper region of the furnace leading to early burnout under high temperatures.Figure 10
Experimental (adapted from [5]) and numerical maps of CO2 concentration. Dimensions on the the left are given in centimetres.Figure 11
Experimental (adapted from [5]) and numerical maps of O2 concentration. Dimensions on the the left are given in centimetres.Figure12 shows a comparison between the experimental measurements and the numerical predictions of O2 and CO2. The flame showed considerable asymmetry; therefore absolute comparisons between measurements and predictions should be regarded as qualitative. In order to minimise the asymmetry, points and corresponding error bars show the mean and standard deviation of measurements performed at opposing, radially symmetric locations. The predictions from the cases SC and DC show that the devolatilization model has a pronounced effect in the flame structure, whereas the char conversion mode model does not show significant differences. When SFOR is used, cases SC and SV, it is clear that the overall particle heating and reaction rates are overpredicted, given the underprediction of O2 and the overprediction of CO2. Furthermore, the particle burnout at the exit of the domain in those cases was close to 100%, but qualitative analysis of the experimentally collected char at the exit of the furnace showed that biomass, particularly the knees, showed limited burnout levels. When using DAEM, the predictions are closer to the experimental results and burnout levels at the exit of the computational domain are in the order of 60%.Figure 12
Comparison between measurements and predictions of O2 (left) and CO2 (right) for the modelling cases SCc, DCc, and SVc. Empty points and corresponding error bars show the mean and standard deviation of measurements performed at opposing, radially symmetric locations.


## 6. Conclusions

A study on the effect of the devolatilization model and char combustion mode model on the structure of a large-scale laboratory coal and biomass flame is presented. Two devolatilization models were tested, single first-order reaction and distributed activation energy, and two char conversion modes were tested: power law mode with a model constants and a constant-free model. Single particle simulations under idealised conditions were performed to evaluate the impact of the different models in the absence of turbulence. Subsequently, large eddy simulations were performed and the results were compared with experimental results. Three numerical cases were considered: single first-order reaction with constant char combustion mode, distributed activation energy with constant char combustion mode, and single first-order reaction with variable char combustion mode. Simulations of single particle conversion under idealised conditions showed that both the devolatilization model and the char conversion mode model had a pronounced effect in the conversion of the particle when using typical furnace conditions, namely, heating rates, temperature, and oxygen concentration. Particle scatter data obtained from the large eddy simulations show that devolatilization tends to occur in a narrow region close to the burner region when using the single first-order reaction devolatilization model, while with the distributed activation energy devolatilization model devolatilization is retarded and occurs along the furnace axial distance. This effect has a profound effect on the spatial distribution of O2 and CO2, since a more intense and localised devolatilization stage leads to a more intense gas-phase combustion region, resulting in early burnout, strong radial expansion, and shorter flame when using the single first-order reaction devolatilization model. On the other hand, the distributed activation energy devolatilization model was capable of predicting the flame break observed experimentally, resulting in a more radially compact and longer flame. The authors of the experimental work attribute the flame break to particle dynamics, particularly different trajectories followed by the different particle categories; hence the results of this work motivate further investigations on the impact of the devolatilization model on particle dynamics. In what concerns the flame lift-off, in both cases the flame stabilises close to the burner within the quarl area. Finally, owing to the intense reaction zone close to the burner region, the different models used for the char combustion mode had a marginal impact on the spatial distribution of O2 and CO2.



---

*Source: 7036425-2018-11-01.xml*