Eulerian CFD modelling for biomass combustion. Transient simulation of an underfeed pellet boiler

Energy Conversion and Management 101 (2015) 666–680

Contents lists available at ScienceDirect

Energy Conversion and Management journal homepage: www.elsevier.com/locate/enconman

Eulerian CFD modelling for biomass combustion. Transient simulation of an underfeed pellet boiler M.A. Gómez a,⇑, J. Porteiro b, D. Patiño b, J.L. Míguez b a b

Defense University Center, Spanish Naval Academy, Spain Industrial Engineering School, University of Vigo, Lagoas-Marcosende s/n, 36310 Vigo, Spain

a r t i c l e

i n f o

Article history: Received 22 February 2015 Accepted 1 June 2015 Available online 18 June 2015 Keywords: Biomass Combustion Boiler CFD

a b s t r a c t This paper describes a transient model for biomass combustion in a fixed bed boiler. This method implements several submodels which address a variety of conversion processes and interactions between solid and gas phases. A set of Eulerian variables will be defined representing the solid phase state as well as the governing equation model for both the evolution and the thermal conversion of the bed. The solid phase components including the ash content is divided into solid and volatile elements. The fuel feeding is modelled by an advective flux term in the transport equations of the solid phase variables. The advective flux term includes a treatment that prevents the numerical diffusion beyond the bed’s surface. This method also includes heat and mass transfer models between phases, particle and bed shrinkage, porous media and gas reactions. Several experimental tests have been simulated to contrast model behaviour. The primary variable profiles of the solid phase in the bed and gas phase in the furnace have been analysed. The results show reasonably good predictions for the exchanged heat and the flue gas concentrations. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction The study of biomass combustion is currently receiving increased interest due to the effects of greenhouse gas emissions and the lack of fossil fuel resources in certain countries. Understanding biomass behaviour in the combustion process is a key factor in designing and improving combustion systems such as boilers, burners and furnaces. CFD techniques have become efficient tools in the search for a better understanding of the biomass thermal conversion phenomena. However, CFD commercial codes still require extensive development regarding the combustion and/or gasification of solid fuels, particularly those that involve the thermal conversion of solid biomass. Numerous works have proposed a wide variety of strategies to implement several submodels into commercial CFD codes to couple the thermal conversion of the solid phase with the CFD simulation of the remaining processes. A recent review of the numerical modelling of combustion systems was presented by Karim and Naser [1]. Most published CFD models of biomass combustion systems are based on compact domains, where the combustion of the packed bed occurs through either a zero- or one-dimensional approach. Porteiro et al. [2] used a zero dimensional bed which took into account the ⇑ Corresponding author. Tel.: +34 986 812604; fax: +34 986 812 201. E-mail address: [email protected] (M.A. Gómez). http://dx.doi.org/10.1016/j.enconman.2015.06.003 0196-8904/Ó 2015 Elsevier Ltd. All rights reserved.

individual contribution of each particle of the bed in a transient simulation. The mass, energy and species were exchanged by the bed and the furnace through a boundary surface. In [3,4] a similar approach was used in a steady and computationally fast model that implemented the zero-dimensional bed inside the CFD domain. Other works use a one-dimensional bed model coupled with the CFD domain by a boundary surface to simulate the biomass advancing along the grate where it is consumed. This method was applied by other authors [5–7] who modelled the one-dimensional bed as a ‘‘walking column’’ of biomass that advances to the next position in the grate while the thermal conversion of biomass takes place. Zhou et al. [8] extended the Van der Lans [5] model to a more detailed account of the properties in order to achieve more accurate predictions. Yang et al. [9] used the finite volume method in a grid with the shape of the bed and resolved the transport equations of several variables which established the state of the biomass at each point of interaction with the gas phase in the furnace through their boundary surface. Shin and Choi [10] simulated an experimental packed bed of waste by calculating the solid phase variables inside the CFD domain and calculating the interaction between the solid and gas phases in each computational cell. Cooper and Hallett [11] also meshed the char bed to perform a one-dimensional simulation in which the char was fed from above and the particles advanced downwards. Thunman and Leckner [12] meshed the solid phase of a bed and

M.A. Gómez et al. / Energy Conversion and Management 101 (2015) 666–680

667

Nomenclature Ai Afi Ain Av b Cp CD deq Dii Ei F fv g G1 h hs k K kg Km LH LVH Mi _ 00f m P Pt R Ri,kin Re Rd St Sh Sc S T t th V v vs

pre-exponential factor (s1) area of the face i (m2) area of the inlet (m2) area–volume ratio (m1) blowing parameter (–) specific heat (J kg1 K1) drag coefficient (–) equivalent diameter (m) diffusivity of the specie i in the mixture (m2 s1) activation energy (J mol1) fouling conductivity factor (–) porosity (–) gravity acceleration (m s2) mass flux of the gas current through the cell (kg m2 s1) convection coefficient (W m2 K1) solid phase enthalpy (J kg1) thermal conductivity (W m1 K1) char reaction constants (m s1) geometric constant () mass transfer constant (m s1) latent heat (J kg1) lower heating value (J kg1) molecular weight (kg kmol1) mass flux of the released gases (kg m2 s1) pressure (Pa) power of the boiler (W) ideal gas constant (J K1 mol1) kinetics equation of the reaction i (reaction dependant) Reynolds number (–) reduction factor (–) Stanton number (–) Sherwood number (–) Schmidt number (–) source term (W m3) temperature (K) time (s) thickness (m) volume (m3) gas velocity (m s1) solid phase velocity vector (m s1)

calculated the effect of elutriation and shrinkage of biomass during thermal conversion. Peters and Bruch [13] used a discrete phase model (DPM) in which they applied the thermal conversion of a single particle and used them in the simulation of a one-dimensional packed bed that took into account the mechanics of the contact between particles. Yang et al. [14–17] employed a higher level of discretisation to model a two-dimensional packed bed of thermally thick particles. Rather than developing an experimental setup to study the entire combustion process at the industrial scale, some of the above mentioned researchers have constructed facilities that were impractical from the energy-production point of view but more convenient and suitable for the study of certain subprocesses which would be otherwise difficult to isolate or control [8–13]. Furthermore, most of these works focused on the prediction of species and temperature fields in a steady state; however, only a few cases have tackled the transient evolution of the combustion system. Collazo et al. [18] and Gómez et al. [19] simulated a three-dimensional burner in which the solid phase was modelled through a set of Eulerian variables in which transport equations represented the evolution of the

Greek symbols e solid fraction (–) g boiler performance (–) l viscosity (kg m1 s1) q density (kg m3) s fraction of heat received by the particle employed in drying (–) w cell filling parameter (–) u char oxidation parameter (–) x_ 000i generation or consumption rates of the wood components (kg m3 s1) Subscripts bl blowing p particle s solid g gas eff effective c consumption G generation in inlet moist moisture wood dry wood char char ash ash evap evaporation glob global C carbon dep deposit vol volatile ash,v volatile ash cr critical Superscripts ox char oxidation g,1 gasification reaction with CO2 g,2 gasification reaction with H2O s specific ash,v volatile ash surf surface vap vapour

biomass in the bed. Along with other models of heat transfer, this study employed a compaction model that shrinks the layers of the bed that are weakened by the combustion process. Thunman et al. [20] proposed an Eulerian approach for calculating thermally thick beds that used a subgrid scale inside the particles and reduced the number of grid points to the number of stages in biomass thermal conversion. Mehrabian et al. [21,22] applied the moving-layer methodology proposed by Thunman et al. [20] at the particle scale to represent the thermal conversion of each layer by employing a dense discrete phase model (DPM) to simulate the biomass particles and their interaction with the gas phase. Ström and Thunman [23] used the same layer model to implement a three-dimensional single-particle model. This particle approach is appropriate when addressing a thermally thick difficulty at a reasonably low computational cost. Chaney et al. [24] has shown an overview of small-scale biomass pellet boiler simulations; different approaches were used in this review to introduce the effect of bed combustion in the CFD framework. Several authors [2,5–7,25] used the top of the bed as a gas inlet to the furnace where the by-products of the solid fuel

668

M.A. Gómez et al. / Energy Conversion and Management 101 (2015) 666–680

Fig. 1. Scheme of the commercial boiler used in the experiments.

Table 1 Composition and properties of the fuels. Proximate analysis

Fuel 1

Fuel 2

Fuel properties

Fuel 1

Fuel 2

Moisture [wt%, a.r.]

6.8

8.5

CH1.57O0.74

CH1.55O0.6

1200

1200

8.8

8.8

74.7

Equivalent formulation (dbaf) Density (q) [kg m3] Equivalent diameter (deq) [mm] Sphericity

Ash [wt%, a.r.]

0.3

0.6

Fixed Carbon [wt%, a.r.]

24.4

16.2

Volatile matter [wt%, a.r.] LHV [MJ/kg, a.r.]

68.5

0.85

0.85

16.17

16.56

Solid fraction (e)

0.56

0.57

it is generally insufficient to validate the bed models. On the other hand, some other experimental works showed detailed data taken from the bed zone and thus are useful to validate the bed submodels. Dieguez-Alonso et al. [31] characterised the species released during pyrolysis through Laser-Induced Fluorescence. Porteiro et al. [32,33] experimented with several solid fuels in a biomass burner to analyse the combustion properties of each fuel and studied the combustion fronts and ignition rates in different stoichiometric conditions. Girgis and Hallett [34] measured gas and tar concentrations, temperatures and bed properties to model devolatilisation kinetics. Our previous approach [19] successfully establish a comprehensive bed model into a commercial CFD code (Ansys-Fluent) to investigate the transient combustion of biomass. The model was validated with the data obtained in an experimental combustor. This paper extends the models presented in [19] by simulating the fuel transport that occurs inside the bed and model the feeding of fuel. This work applies a transient three-dimensional model to the simulation of an up-to-date biomass domestic boiler under different working conditions to reach its operational steady state. Experimental tests are used to contrast the model predictions. The solid phase is modelled using Eulerian variables that are defined and calculated through transport equations that model the thermal conversion of the solid particles. The treatment of the solid and gas phases as well as the heat transfer and bed compaction models follow the formulation used in our previous work [19]. However, in that work the model was validated against a batch-fed pilot facility, hence only bed consumption was considered. In the work presented here, a feeding model is introduced and is coupled with bed consumption and compaction and adapted to commercial boiler geometry and operational conditions. The fuel-feeding model accounts for the fuel movement inside the bed from the fuel inlet to the consumption region. This work also accounts for the ash content of the fuel and tracks ash transport, consumption by partial evaporation and accumulation within the boiler furnace as an initial step towards the study of slagging and fouling processes. Although local ash concentrations and temperatures were obtained as results of the model, this feature is still under development and requires further experimental and numerical probing.

2. Experimental boiler and fuels combustion were introduced. These products introduce the mass, energy and species that were previously calculated by the bed submodels. This methodology is commonly employed to couple a 0-D or 1-D model to a CFD gaseous domain. For instance, Eskilsson et al. [26] used experimental measurements to select the composition and properties of the gas leaving the bed. However, several authors simulated the packed bed region as a porous zone and introduced the combustion products as volumetric sources [7,27,28]. A further advancement was the introduction of these mass, energy and species sources through user-defined functions (UDFs) that allowed them to be dependent on the embedded models, thereby accounting for local conditions when calculating the behaviour of each specific point of the packed bed. Some recent numerical works [29,30] have used CFD codes to describe different phenomena as biomass co-combustion, radiative heat transfer or swirl flow, which broadens the knowledge of combustion modelling to other processes. However, despite the effort made during the last decades, there is still a lack of comprehensive empirical data necessary to analyse all the mechanisms involved and validate the submodels developed in order to account for them. A limited number of studies have included measurements from biomass boilers [4,6,7], however, these data is collected far from the reaction zone and hence,

The experimental boiler used in this study is a commercially available 60 kW pellet boiler (KWB Multifire) that is schematically represented in Fig. 1. The fuel is introduced into the bed from the bottom by a screw. The fuel is forced to ascend towards a circular grate, where it spreads while air is introduced through the grate by several small holes arranged in a circular pattern. Secondary air is then introduced by a ring located over the fuel bed that injects the air towards the centre with a tangential component to generate swirl and intensify the turbulence of the flame. A metal dome is located over the secondary air ring to prevent the flame from reaching the water walls of the heat exchanger. The heat exchanger, which is composed of fifteen tubes, is placed over the furnace. Inside each tube, several spiral praises are installed to force the gases to spin, thus improving the convective thermal transfer. The furnace walls are refrigerated, with the exception of the frontal wall, where a door is placed. Several tests were performed in the boiler with different operating conditions and two different fuels. The properties of the two fuels used are shown in Table 1. The operating conditions of the boiler depend on several parameters that strongly influence boiler performance, including the fuel feeding rate, the air–fuel excess ratio, the distribution of the primary and secondary air

669

M.A. Gómez et al. / Energy Conversion and Management 101 (2015) 666–680 Table 2 Operating conditions for the different experimental tests. Test

Fuel

Fuel feeding rate (g s1)

Primary air flow rate (g s1)

Secondary air flow rate (g s1)

Air infiltration (g s1)

E1 E2 E3 E4

F1 F1 F2 F2

4.02 3.69 2.06 2.06

10.4 (30%) 4.2 (15%) 4.7 (13%) 7.3 (22%)

14.5 (41%) 14.3(50%) 9.6 (27%) 8.0 (24%)

10.2 (29%) 9.9 (35%) 20.6 (60%) 17.6 (54%)

fluxes and the water temperature. Air infiltration is unavoidable in conjunction with these parameters and is measured and regulated by controlling the chamber pressure through a combination of air supply and induction fans (this is further explained in Section 4.1). Four different experimental tests are used to contrast the simulation of the boiler with the proposed models. The operating conditions of the boiler in these tests are shown in Table 2. The water temperature of the boiler was 70 °C for all of the tests. 3. Model description The present paper is an extension of the author’s previous work [19] in which the equations that model the solid phase are adapted to simulating a boiler. As the thermal behaviour of the solid phase and its interaction with the gas phase were previously modelled, the variables transport, the addition of fuel feeding and the inclusion of ash movement, evaporation and deposition in the bed are the major improvements in the present model. This chapter describes the primary features and equations, with special attention paid to parameters that affect the fuel feeding and ash modelling; the compaction, heat and mass transfer models have already been detailed [19]. 3.1. Assumptions The model is based on the following hierarchical assumptions, which are sorted based on the anticipated impact of the hypothesis on the model results. This order should dictate the aspects of the model to be addressed in future works and its validity of application.  Thermally thin spherically equivalent particles are assumed to exist in the bed with the following corollaries: – The product gases generated inside the particles are instantaneously released into the gas phase, leading to an immediate outflow condition. – The gases emitted by the solids are introduced into the gas phase at the reference temperature (condition imposed by the CFD code), and the temperature of the removed solid phase is compensated by an energy source.  The bed is modelled as a porous dispersed medium with local volume-averaged properties within each computational cell. One of the limitations of this volume-averaged approach is that particle size should be much smaller than the cell size so that each cell can be considered a porous media made of small particles. Due to the geometric complexity, the mesh used in this work contains cells smaller than particles, which is a infringement of the theoretical basis of volume-averaging for porous media. Therefore, each cell value must be considered as the representative values of a volume averaged specific region of the bed and not of the particles inside the cell.  The density and solid fraction of the solid particles varies during the drying, devolatilisation and char reactions, while the solid fraction varies during the devolatilisation and char reactions. Therefore, particle shrinkage is only produced in the last two stages.  Fuel feeding is modelled as a continuous flux of the solid phase.

 Local collapses in the bed cause compaction that is driven by the critical porosity being locally achieved.  Drying is assumed to be thermally controlled and to occur at a specific temperature.  The ash content is divided in volatile and solid ash. The volatile ash content is consumed on the particle surface through a partial pressure balance.  The solid ash content is considered to be non-reactive and is dispersed for the entire particle volume. The dragged solid ashes are instantaneously removed from the bed. 3.2. Solid phase CFD commercial codes are prepared to calculate the gas phase combustion. The modelling of the solid phase combustion requires the implementation of several submodels that calculate the evolution of the major variables involved in the thermal conversion of the particles. Authors in previous works [18,19] proposed a set of variables that modelled the presence of biomass solid particles in the bed region and their interaction with the gas phase during the combustion process. In this work, eight scalar variables were defined to represent the major characteristics of biomass conversion. These variables are as follows:  Enthalpy of the solid (h).  Local solid fraction (e).  Volume of the particle, as represented by the third power of its diameter (dp).  Local moisture density (qmoist).  Local dry wood density (qwood).  Local char density (qchar).  Local volatile ash density (qash,v).  Local ash solid density (qash,s). Table 3 shows the transport equations of these scalars (Eqs. (1)– (8)) and the relationships between the wood component densities and the total density (Eq. (9)). In the scalar transport equations (Eqs. (1)–(8)), the first term on the left-hand side denotes the transient term that represents the temporal variation of the variables. The second term is the advective term that accounts for the transport of the biomass flux from the feeding process. Because there is no mass diffusion between the particles, the diffusive term is only present in the enthalpy equation, thus accounting for the equivalent conductivity between the solids. The term Ss in Eq. (1) and the entire right-hand side of Eqs. (2)–(7) represent the generation or consumption of each variable. Eqs. (2) and (3) illustrate a significant difference with the model presented in our previous work [19] regarding the source term. Table 3 Equations of the solid phase variables. Solid phase enthalpy Solid fraction

@ðeqp hs Þ @t

Third power of particle diameter Moisture density

3 @deq

Dry wood density Char density Volatile ash density Solid ash density Total particle density

@e @t

þ rðe  qp  Rd  v s  hs Þ ¼ rðks;eff  rT s Þ þ Shs

þ rðRd  v s  eÞ ¼



_ x

000 wood

qswood þ





_ 000  _ 000 G;char c;char s char _ 000  _ 000 _ 000 G;char c;char wood s s wood char

x

x

q

(1) (2)

e



3

(3)

@ðeqmoist Þ þ rð @t @ðeqwood Þ þ rð @t @ðeqchar Þ þ rð @t @ðeqash;v Þ þ rð @t

e  Rd  v s  qmoist Þ ¼ x_ 000moist e e  Rd  v s  qwood Þ ¼ x_ 000wood e e  Rd  v s  qchar Þ ¼ ðx_ 000G;char  x_ 000c;char Þe

(4)

e  Rd  v s  qash;v Þ ¼ x_ 000ash;v e

(7)

@ðeqash;s Þ @t

(8)

@t

3

þ rðe  Rd  v s  deq Þ ¼

x

q

þ

x

x

q

þ rðe  Rd  v s  qash;s Þ ¼ 0 qmoist + qwood + qchar + qash,v + qash,s = qp

deq

(5) (6)

(9)

670

M.A. Gómez et al. / Energy Conversion and Management 101 (2015) 666–680

Table 4 Equations of the solid consumption rates. q C

Drying rate

(10)

x_ 000moist ¼ s LHpmoistp @T@ts ; T s P T ev ap

  x_ 000wood ¼ qwood i¼1 Ai exp  RTEi s   E3 x_ 000G;char ¼ qwood A3 exp  RT s P3

Devolatilisation rate Char generation rate

(11) (12)

g;1 g;2 x_ 000c;char ¼ K ox glob Av ½O2 M C þ K glob Av ½CO2 M C þ K glob Av ½H2 OM C

Char consumption rate

(13)

Table 5 Kinetics of the solid heterogeneous reactions [39]. Pyrolysis reactions

Kinetics

Drywood ? Gas Drywood ? Tar Drywood ? Char

A1 = 111  109 (s1), E1 = 177  103 (J mol1) A2 = 9.28  109 (s1), E2 = 149  103 (J mol1) A3 = 30.5  109 (s1), E3 = 125  103 (J mol1)

Heterogeneous char reactions

Kinetics

(R.1)

C + uO2 ? 2(1  u)CO + (2u  1)CO2

(R.2)

C + CO2 ? 2CO

(R.3)

C + H2O ? CO + H2

This source term is now a function of the relationship between the wood consumption and char generation and the consumption with their reference densities. We consider the reference density of a component (wood, char or ash) as the density of a particle completely created from that component; the values empirically obtained for each fuel were used as constants throughout the model. This formulation calls for shrinkage of the particles during both the devolatilisation and char consumption stages, and the resulting fractions of wood and char fulfil the ratios of wood and char densities in both the particles and their respective specific densities. The advective term of the scalar transport equation represents the energy carried by the movement of the solid phase particles and is related to the velocity vector vs in each cell. The calculation procedures for this vector and the parameter Rd (which is involved in the advective flux) are explained in Section 3.4. _ 000 The terms x i that are involved in the source terms of Eqs. (2)– (7) represent the generation or consumption rates of the wood _ 000 components. Eqs. (10)–(13) show the terms x i without the ash consumption, which is explained below. The devolatilisation process is affected by several factors such as solid temperatures, heating rates or thicknesses of the particle reaction layers [35,36]. The dry wood consumption is modelled as a conversion into gas, tar and char [37]. Therefore, the devolatilisation rate considers three different Arrhenius rates that correspond to each conversion process (Eqs. (11) and (12) in Table 4). The char consumption rate (Eq. (13)) is modelled as the addition of the direct char oxidation, (R.1) [20], and two gasification reactions, (R.2) and (R.3) [38]. The kinetics of the heterogeneous

  K ox ¼ 1:715  T s  exp  9000 Ts   4 K g;1 ¼ 3:42  T s  exp  1:5610 Ts   4 K g;2 ¼ 5:7114  T s  exp  1:5610 Ts

g;1 reactions are shown in Table 5. The parameters u, K ox and g ; Kg

K g;2 g control the char oxidation CO/CO2 ratio (Eq. (14)) and the thermal and diffusion of oxygen, carbon dioxide and water vapour, respectively (Eq. (15)).

  2 þ 4:3exp 3390 Ts   u¼  2 1 þ 4:3exp 3390 Ts K ox glob ¼

1 K ox

1 1 ; K g;1 glob ¼ 1 þ K1ox g;1 þ K m

ð14Þ

1 g;1 Km

; K g;2 glob ¼

1 K

g;2

1 þ

1 g;2 Km

ð15Þ

3.3. Ash This work also examines the evolution of ash in the bed; the consumption rates are described in Eqs. (16) and (17). The ash content is divided into volatile and solid mineral ash, as proposed by several authors [40–43]. A partial pressure balance is used to model the evaporation of the volatile content; KCl was chosen as a representative species of the volatile content due to the reasonably high potassium (K) content in the analysis of ashes in the commercial pellet employed for this work. The volatile ash consumption is modelled as shown in Eq. (16), which considers both the mass transference from the particle surface and the volatile ash concentration [44]. Eq. (17) shows the concentration of volatile ash on the particle surface that functions as the vapour pressure equilibrium for the volatile ash [44]. The mass transfer coefficient

Table 6 Combustion reactions and kinetics. Homogeneous reactions

Kinetics

(R.4)

C6 H6 þ 92 O2 ! 6CO þ 3H2 O

(R.5)

CH4 þ 32 O2 ! CO þ 2H2 O

(R.6)

H2 þ 12 O2 ! H2 O

(R.7)

CO þ 12 O2 ! CO2

(R.8)

H2O + CO ? CO2 + H2

(R.9)

CO2 + H2 ? H2O + CO

  8 ½C6 H6 0:1 ½O2 1:85 R1;kin ¼ 1:3496  109 exp  1:25610 RT   11 2108 R2;kin ¼ 5:012  10 exp  RT ½CH4 0:7 ½O2 0:8   7 ½H2 ½O2  R3;kin ¼ 9:87  108 exp  3:110 RT   12 1:702108 R4;kin ¼ 2:239  10 exp  ½CO½O2 0:25 ½H2 O0:5 RT   1:255107 R5;kin ¼ 2:780 exp  ½H2 O½CO RT   7 ½CO2 ½H2  R3;kin ¼ 93690 exp  4:65910 RT

M.A. Gómez et al. / Energy Conversion and Management 101 (2015) 666–680

Fig. 2. Representation of the cell filling approach.

(Eq. (18)) was calculated from the Sherwood dimensionless number with a Whitaker correlation [45] (Eq. (19)). We assume that the solid ash content is dispersed in the biomass particles that are released when the particle is consumed. Part of that solid ash is elutriated by the gas flow while the balance remains in the bed. The solid ash is considered to be non-reactive in this work. Therefore, the only mechanism for removing the ash from the bed is the aerodynamic drag, apart from the possibility of rich ash regions of the bed reaching the boundaries of the grate and falling to the bottom of the furnace.



v gas x_ 000ash;v ¼ Av  kash;  qsurf m ash;v  qash;v v ap qsurf ash;v ¼ P ash;v ash;v

km

¼



ð16Þ

M ash;v R  Ts

ð17Þ

Sh  Diash;v deq

ð18Þ

 0:25   l1 Sh ¼ 2 þ 0:4  Re0:5 þ 0:4  Re2=3  Sc0:4 

ls

ð19Þ

When the biomass is consumed, the solid ash remains as a residue. A portion of the solid ash content of a bed region is maintained inside particles that have yet to be consumed, and the remaining content is thus detached from the particles. The solid ash content is affected by the level of bed compaction and any bed advective movements that cause displacements in the grate. However, the ash is removed from the bed not because of the combustion process but through aerodynamic drag or by reaching the borders of the grate. To address the elutriation of the solid ashes, a function was created to modify the solid ash density variable if the gas vertical velocity was sufficient for dragging the ash particles. Because the particles are assumed to be spherical, a relationship between the gas velocity and the critical diameter must be established (Eq. (20)). Thus, the solid phase content is removed for all cells whose particle diameter is smaller than the critical diameter. At this point of the combustion usually most of the solid phase content is ash since the gas velocity in the bed usually reach low values and hence when the particle reaches its critical size, it is almost completely depleted.

dcr ¼

3C D qg v 2 4qp g

ð20Þ

3.4. Gas phase Gas phase modelling has been extensively used in combustion simulations for a large variety of systems [46,47]. Commercial

671

CFD codes are typically equipped with a set of models that solve the variables involved in gas dynamics. These variables include pressure, momentum, energy, turbulence, radiation transport, chemical species, soot formation and combustion. In this work, the models used for these variables were the modified discrete ordinate modelling for radiation previously presented by the authors [48], realisable k–e modelling with enhanced wall treatment for turbulence, Finite-rate/Eddy-dissipation (FR-ED) modelling for species reaction and Moss-Brookes modelling [49] for soot nucleation and reaction. The flow of gas through the porous zone was modelled with a dispersed porous media approach. The effect of the bed on the gas is considered through a momentum source [3,18,19] and a physical velocity condition that accounts for the reduction of the flow-section caused by the presence of the solid phase. As in previous works, the thermal conversion of the solid phase is modelled using consumption and generation sources of solid phase variables. These sources include the gaseous emissions of H2O from drying, the volatiles from devolatilisation and the CO, CO2 and H2O released in the heterogeneous char reactions. The emitted volatile species are primarily composed of light hydrocarbons (represented by CH4), tars (represented by C6H6), CO, CO2, H2 and H2O [50]. The method of obtaining the volatile species fractions is a variation of those proposed by Thunman et al. [50] and Neves et al. [51]. This method searches the combination of species fractions that minimises the error of the C, H, O and energy balance. The combustion scheme consists in a set of chemical reactions where the hydrocarbons from the wood devolatilisation are oxidised in two steps (reactions (R.4)–(R.7)) [2,3,18,19] with a two-way reaction of water and carbon monoxide ((R.8) and (R.9)) [52,53]. The reactions and kinetics are shown in Table 6. As previously mentioned, the modelling of the heat and mass exchange between the phases in this work is based on [19]; the new approach involves the blowing of the gases, leaving the particles that reduce the effect of the heat convective exchange and mass transport between the particles and gases. This effect is calculated by reducing the convection coefficient as proposed by Ouedraogo et al. [54]. The convection coefficient that considers blowing follows Eq. (21), which corrects the convection coefficient without blowing (h) using a shape factor and a blowing parameter (b). This blowing parameter is a function of the mass fluxes of the outgoing gases from the particle, the free stream gases and the Stanton number (Eq. (22)).

hbl ¼ h

b¼

kg  b expðbÞ  1

_ 00f m G1 St

ð21Þ

ð22Þ

3.5. Advective motion generated by the feeding process Because the simulated boiler in this paper is fed from the bottom through a feeding screw that pushes the fuel upwards, the entire bed is subjected to movement. This movement affects the solids and must therefore be accounted for in the scalar Eqs. (1)– (8) by the second left-hand term. This term is the divergence of the flux in the scalar quantity. The parameter Rd, as explained below, is introduced to avoid the diffusion of the solid phase beyond the bed surface. The simulations in this work were successfully computed using first and second order upwind schemes for the spatial discretisation of the solid phase variables. Although higher order discretisation schemes such the QUICK scheme or the Third-Order MUSCL

672

M.A. Gómez et al. / Energy Conversion and Management 101 (2015) 666–680

Fig. 3. Representation of the changes in the filling parameter (w) in a sequence of the fuel feeding.

Fig. 4. Simulation of the bed advance without and with the cell-filling model.

Scheme can be used if more accuracy is desired, these present a higher numerical instability and convergence is usually difficult to reach, especially when using unstructured meshes. For simulations using multiphase models, upwind schemes are generally unsuitable for interface tracking because of their overly diffusive nature [49]. Central differencing schemes while generally able to retain the sharpness of the interface are unbounded and often give unphysical results. In order to overcome these deficiencies, the CFD

code ANSYS FLUENT offers a modified version of the High Resolution Interface Capturing (HRIC) scheme for multiphase simulations. The modified HRIC scheme is a composite NVD scheme that consists of a non-linear blend of upwind and downwind differencing [55]. However, in this work the solid phase was added by introducing a porous region in the single-phase model in which the modified HRIC is not available. The approach of the porous region has proved to be very efficient to model biomass beds,

M.A. Gómez et al. / Energy Conversion and Management 101 (2015) 666–680

673

Fig. 5. Representation of the mesh and boundary conditions.

and among others, we are exploring that solution. One of its drawbacks is that the implementation has to deal with some intrinsic hypothesis of single-phase formulations, such as the one we are discussing. The mesh used in this work is an unstructured one in which, in general, the flow direction does not match the alignment of cells. In addition, the time step size needed to make the simulation feasible in computation time is higher than 1 s, which does not favour convergence. Several attempts to calculate the bed movements were performed to choose the most accurate discretisation scheme. Only first and second order upwind schemes obtained convergence. Therefore, the second order upwind was chosen to compute the solid phase variables. To confront the problem of interface diffusion, a simple cell filling model was introduced, which also uses an upwind scheme to prevent bed surface distortion. The velocity vector was calculated in previous simulations, where the fuel was modelled as a non-reactive flow that represented the solid phase located in the bed zone. This flow is modelled with the density of the packed bed and high viscosity values. Once the simulation of the fuel motion converged, the coordinates of the velocity vectors were stored as static memories. These were later used in the complete boiler simulation to be applied in the solid phase scalar equations. Mehrabian et al. [27] performed a simulation with a similar process in which an Euler-granular model was used to simulate the fluid that represented the bed motion. However, that method was discarded due to difficulties in convergence. In a fuel motion simulation, the velocity of the fuel flux is imposed on the fuel inlet by considering the boiler power, the boiler performance, the inlet surface and the properties of the virgin fuel (Eq. (23)).

v in ¼

Pt LHV  q0s  e0s  Ain  g

ð23Þ

When the velocity field vectors are defined, the bed feeding commences by applying a flux into the transport equations. When doing this, the bed advances in a diffusive way, blurring

its surface because the solid flux is not constrained to the bed by gravity; therefore, the bed will slowly spread onto the entire gas volume of the boiler. To address this problem, a parameter w (cell filling) is introduced into the model that represents the fraction of the volume filled by the bed (solid and voids). This parameter is applied to each cell of the boiler, as represented in Fig. 2. Bed cells are assigned with w = 1, while gas cells receive w = 0 and the bed’s surface 0 < w < 1. In this way, the non-physical advance of the bed is avoided because the feeding flux through a face is not allowed unless the emitting cell is totally filled. The working of the bed advance model is schematically represented in Fig. 3. The filling (w) is also modelled by a transport equation (Eq. (24)) with transient and advective terms. The value of w is obtained by integrating this equation in each cell volume (Eq. (25)). These two equations (as well as Eqs. (1)–(8)) are affected by the parameter R. This parameter represents the reduction of the feeding flux that occurs in a cell and prevents the diffusion of the bed beyond its surface. Therefore, when w < 1 (surface or gas cells), Rd is null. The Rd also reduces the feeding flux when some of the upwind cells are not filled (case t1 in Fig. 3). In addition, the Rd parameter reduces the flux that occurs when a partially filled (w < 1) cell in a previous instant (t  Dt) becomes full during a time step. In this case, the outgoing flux must be reduced with respect to the incoming flux to compensate for the volume storage in the cell during the time step (1  wtDt). The three tasks of Rd are compiled in Eq. (26). The parameter Rd affects all of the solid phase variables (Eqs. (1)–(8), (24) and (25)) where the face value Ri is required. Because this follows an upwind formulation, the face value of Rd is the upwind cell value (Eq. (27)). The effect of the cell-filling model is shown in Fig. 4, which shows a simulated porous bed that is fed from the bottom. In the first case, the cell-filling model is not applied, thus showing the non-physical advance of the bed (only the solid fraction is represented, while all of the solid phase variables are affected). In the second case, the cell-filling model allows the bed surface to advance more realistically.

674

M.A. Gómez et al. / Energy Conversion and Management 101 (2015) 666–680

Fig. 6. Gas temperature, velocity contours and velocity vectors of the primary and secondary air injections.

@w ¼ rðRd  v s Þ @t w ¼ wtDt þ

Rdcell

ð24Þ

N faces Dt X Rdi  v s;i  Afi V cell i

8 PN ð1wtDt ÞV cell faces > < i ðRdi v i Afi Þ Dt PNfaces ¼ v i Afi i > : 0

ð25Þ 9 > if w ¼ 1 = > ; if w < 1

Rdface ¼ Rdcell ðupwindÞ

ð26Þ

ð27Þ

4. Simulation 4.1. Boiler mesh and simulation conditions In order to reach the mesh independency, several simulations in identical conditions were performed and the critical zones (bed region and immediately above) were refined until the differences between the simulations were negligible. The definitive boiler model is discretised in a mesh of approximately 1  106 polyhedral elements, with a special refinement in the zones where stronger gradients are expected, such as the bed region, the near wall region and the areas surrounding the air entrances. The number of

elements employed is far shy from what might be expected in a steady simulation, but the higher computational cost of transient simulations necessitated that a compromise be found between element size and cost. The average length of the 3D elements is 2.8 mm in the bed region and its proximities and 8.1 mm in the rest of the boiler. The average orthogonal quality is approximately 0.9. A representation of the mesh is shown in Fig. 5; the major zones and boundary conditions are also indicated. Three different zones act as mass airflow inlets; the primary and secondary air inlets and another surface located in the boiler floor that represents the infiltrations due to the negative pressure existing inside the boiler. Though enormous efforts have been made to avoid the entrance of air from the surroundings and completely control the fluxes entering the experimental boiler, experience has shown that this is almost impossible to achieve; every type of commercial boiler tested in our facility experienced infiltration from either the ash removal system, the fuel-feeding system or boiler body itself. Therefore, it has become common practice in our simulations to derive the infiltration air from the experimental balance between the airflow rate measured at the metres and the O2 balance at the chimney. Table 2 shows the mass flow rates for primary, secondary and infiltration inlets observed in each test and used as boundary conditions in the simulations. A reduction in size and complexity of the grid is performed to allow the transient simulation of the entire boiler. The water-walls of the heat exchanger and the furnace are modelled

M.A. Gómez et al. / Energy Conversion and Management 101 (2015) 666–680

Fig. 7. Gas temperature field in the middle plane of the boiler for the four simulated working regimens.

as a metal wall with an outer convective coefficient of 2500 (W m2 K1) and a water temperature of 70 °C. The chimney, which is surrounded by still air, is modelled with an outer convection coefficient of 12 (W m2 K1). These are common values of convective heat transfer for flowing water and still air, respectively [45]. The remaining external walls of the boiler are reasonably well insulated and are modelled as adiabatic. The emissivity of the internal walls is assumed to be 0.8 due to the expected level of

675

fouling. The average wall material is steel with a thickness of 5 mm. However, the thermal conductivity of the walls has been modified to account for any deposits that reduce the effective conductivity (Eq. (28)). The correlation proposed by Richards et al. [56] was used in this work to calculate the deposit thermal conductivity as a function of the gas and ash conductivities and the porosity of the deposit (Eqs. (29) and (30)). The deposit thickness was measured in the experimental boiler with values in the range of 1 mm. Because the porosity of the deposits is very difficult to measure, a value of 0.7 produced reasonable results and was therefore used for all of the simulations. This porosity (0.7) produces a conductivity of around 0.1 (W m1 K1) for the ash deposit considering an ash composed mostly of silica and KCl. With greater values of porosity the conductivity approaches to that of the gas. The minimum porosity measured by Richards et al. [56] were around 0.4 that would increase the conductivity to approximately 0.45 (W m1 K1). A layer of cells is selected in the bed bottom to maintain fixed initial values of all the solid phase variables because this layer is the fuel feeding origin. The specific density (Eqs. (2) and (3)) of the wood and char components are 1200 kg m3 and 400 kg m3, respectively, and were estimated as typical values based on our own measurements. The solid ash specific density is different for each fuel, which is calculated as a closure value to complete the addition of the component fractions. The moisture and the volatile ash are considered to be dispersed in the particle and their consumption adheres to a shrinkage density approach; therefore, they do not have values of specific density. Because the primary goal of this work is performing a transient simulation from ignition to pseudo-steady performance at each operational point, each simulation begins with an ignition source applied to the bed of raw pellets. This ignition is modelled as a heat source of 1000 W spread on a small area of the bed surface. Once the ignition source is terminated, the bed burns for 300 s without activating the feeding model to reach a stable state of combustion. The fuel-feeding model is subsequently activated, and the boiler operates for at least 2 h of simulation time until the monitored variables show either steadiness or a steady oscillation. The time step is set to 5 s for the entire simulation. The aforementioned simulation sequence attempts to mimic the experimental sequence employed in each test.

keff ¼

thwall þ thdep thwall kwall

þ

thdep kdep

Fig. 8. Solid phase temperature field in the middle plane of the fuel bed for the four simulated working regimens.

ð28Þ

676

M.A. Gómez et al. / Energy Conversion and Management 101 (2015) 666–680

Fig. 9. Contours of the gas temperature and the main species molar fractions in the boiler furnace.

5. Results and discussion 5.1. Results of the model

Fig. 10. Contours of the KCl molar fractions in the boiler furnace.

kdep ¼ ð1  FÞ  kash þ F  kgas 0

F¼

2n B 1 @1   2 1 1þf n

v dep

ð29Þ 1

C n A

ð30Þ

This chapter analyses the major variables of the model to understand the behaviour of the model and boiler. Other relevant data such as the heat fluxes, temperatures and concentrations in the flue gases are compared with experimental records. Fig. 6 shows contours and vectors that help in visualising the general behaviour of the boiler. These results are extracted from the final simulation time step when a pseudo-steady state has been reached. Fig. 6 also shows a general view of the boiler geometry with the gas temperature results in the central planes. The contour of the solid fraction of the fuel bed that determines the boundary of the bed is also visible. The details at the left of the image show the velocity contours inside a heat exchanger tube, the secondary air injection vectors (coloured according to velocity) and the primary air injection vectors that pass through the fuel bed in the grate. Figs. 7 and 8 show the temperature fields of the gas and solid phases for the four simulations in which a steady state was achieved. The differences in the air distribution can be observed in Fig. 7, where the overall furnace temperature was higher in E1 and E2 in which the boiler power is approximately the double than the other two cases. In simulation E2, the flame was more voluminous due to the air distribution; the primary air was insufficient for a complete combustion of the gases released from the bed, and the combustion was thus displaced to the secondary air injection. In cases of lower power (E3 and E4), the flame was concentrated in

M.A. Gómez et al. / Energy Conversion and Management 101 (2015) 666–680

677

Fig. 11. Contours of the solid phase main variables.

Fig. 12. Contours of the solid ash density after 1800 s and 10,800 s of simulation.

the centre of the furnace. This effect may have been produced by a lower flux of released volatiles and a higher value of air infiltrations, which increased the oxygen concentration in the furnace and depleted the reactive species before they escaped from the centre of the furnace. Fig. 8 shows a similar pattern in the profile of the four cases. This was caused by a primary air injection that was located in a crown of orifices around the grate. Therefore, these reactions occurred along the sides of the bed while the centre was not receiving oxygen. The temperature of the solid phase was high in test E1 due to the concentration of reactions caused by a higher primary airflow. This effect was also present when comparing cases E3 and E4, where E4 reached a higher temperature. There are visible differences in the bed size that were caused by different feeding rates (released power) and air distribution. The feeding

rate of E1 and E2 was double the other cases, which favoured a fuller bed. However, a lower primary air flux stimulated bed growth because of lower bed kinetics and slower particle combustion. The air distribution favoured the bed size in cases E2 and E3, where the primary airflow was lower. From this point onwards, all figures are derived from the simulation of test E1. Figs. 9 and 10 show the gas temperature and main species mass fraction fields in the boiler furnace (Fig. 9) and the KCl concentration (Fig. 10), while Fig. 11 shows the solid phase variables in the bed region delimited by the cells whose solid fraction has non-zero values. The first picture in Fig. 9 shows the flame rising from the grate in the area where the primary air was injected while the temperature was lower in the centre of the bed due to the lack of oxygen. The effect of the dome was also visible, which

678

M.A. Gómez et al. / Energy Conversion and Management 101 (2015) 666–680

Fig. 13. Schematic of the experimental setup.

prevented the flame from reaching the cold walls of the heat exchanger; this avoided a high concentration of unburned species in the emissions due to the flame-quenching mechanisms that are typical in other commercial boilers. A higher concentration of hydrocarbons (C6H6 and CH4) was witnessed just above the bed due to devolatilisation and a lack of oxygen in the lower part of the flame. The flame shape was almost completely governed by the amount of oxygen and turbulence available in each region. The O2 concentration had a reverse bias to the combustion, as did CO2. CO was more concentrated just over the bed. The CO was released in devolatilisation and the heterogeneous char reactions caused by the substoichiometric work conditions of the primary air supplied to the bed. Volatile hydrocarbon combustion was not completely developed in this area, which supplied more CO to the flame where the reaction was completed; O2 was introduced by the secondary airflow and by infiltrations that flowed from the bottom of the bed through natural draft. Fig. 10 shows the concentration of KCl in the boiler furnace. The KCl acts as a precursor specie in the formation of ash deposits in the walls. Therefore, the prediction of KCl concentrations in the gas phase can be useful in the deposition models.

Fig. 11 shows that the central region of the bed was not burned, as had been previously shown in the thermal maps of Fig. 8. This was primarily caused by the lack of oxygen in this zone. Otherwise, in the outer region of the bed, the biomass was at an advanced state of degradation because the moisture, dry wood and volatile ash were non-existent and the unique components were char and a solid ash deposit. The solid phase temperature was approximately 1500–1600 K in this area due to the heterogeneous char reactions and the combustion of the gas with the primary airflow. The solid fraction contour shows channelling due to the faster char consumption caused by the jets of primary air. The dry wood and volatile ash contours show concentrations of these components in the bed surface. This was caused by bed compaction, which not only affects the bed surface but also spreads the particles to the sides. Although a pseudo-steady state was reached in the boiler with regards to heat transfer and gas emissions, the bed continuously evolved due to the presence of solid ash. The solid ash was not reactive and could only be removed by the gas stream in areas of higher velocity. Therefore, the solid ash was continuously displaced by the advective movement of the bed towards the exterior area of the grate. Fig. 12 shows the solid ash density contours in two different instances of the simulation (at 1800 and 10,800 s). The concentrations of solid ash are visible in the sides of the grate. The density of the solid ash in this area was 30 kg m3, which is the specific density of the solid ash. The fraction of ash in the particles was 1, i.e., the region was full of ash. This effect has been commonly observed in actual boilers; when stopped after several hours of operation, a ring of solid ash has formed a ring around the exterior of the grate. In their new models, the boiler manufacturer has actually introduced an external mechanism to scrape the ash off the exterior of the grate. 5.2. Experimental contrast The boiler is installed in an experimental facility to measure the main parameters. Fig. 13 represents the experimental setup. Two flowmeters are installed in the primary and secondary air inlets and another one is located in the gas chimney to measure both air injections and air infiltrations. A gas analyser is located in the boiler chimney that reports measurements of CO and CO2 volumetric concentrations considering a regular 10% oxygen concentration as well as the flue gases temperature. The water circuit is also

Useful heat [kW]

CO2 emissions [%]

80 60 40 20 0

11 10.5 10 9.5 9 8.5

E1

E2

E3

E4

E1

E2

E3

E4

CFD

59.6

55

28.5

30.4

CFD

10.2

10.3

10.1

10.1

EXP

60.1

55.3

29.9

31.1

EXP

10.7

10.8

9.7

9.4

Flue gas temperature [ºC] 225 200 175 150 125 100

CO emissions [ppm] 800 600 400 200 0

E1

E2

E3

E4

E1

E2

E3

E4

CFD

194

166

152

154

CFD

1.5

1

269

298

EXP

195

182

196

180

EXP

6

4.9

227

591

Fig. 14. Experimental results vs. model prediction (time averaged values at steady state).

679

M.A. Gómez et al. / Energy Conversion and Management 101 (2015) 666–680

12

10.5

10

10

8

9.5

6

9

4

8.5

2

8

CO emissions (ppm)

CO2 emissions (%)

CO and CO2 emissions 11

0 5

25

45

65

85

105 125 145 165 185 205 225 245 265 285

me (s) CO2 simulaon

CO2 experimental

CO simulaon

CO experimental

Fig. 15. Unsteady CFD predictions and experimental values of CO and CO2 emissions.

equipped with a flowmeter and thermocouples to calculate the heat exchanged in the boiler. The heat is dissipated to the air in a sink. All the data are collected in a LabVIEW program and they are stored with a frequency of approximately 3 s. The fuel mass flux was previously measured by weighing the mass of fuel displaced by the feeding screw in each movement. The experimental measurements of outlet gas temperatures and the composition and heat transfer to the water are examined to evaluate the model. Fig. 14 shows the numerical results of the four simulated tests and their respective experimental data. These data represent the time-averaged values once the steady state had been achieved. The experimental values were averaged for approximately 2 h. The predicted data were collected for approximately half an hour since the simulation reached the steady state. No significant variations were seen in longer simulation times. The predictions of the heat exchanged to the water are reasonably close to the experiments, showing a maximum discrepancy of 1.4 kW; this represents 4.7% of the thermal output of the boiler during the test and 2.3% of its nominal power. The predictions of the outlet temperatures are reasonably accurate for tests E1 and E2, which represent full load tests. However, during tests E3 and E4, the discrepancies are noticeable, which coincides with the cases in which the heat exchange errors were greatest. As only a few quantities are measured and there is a lack of measurements in the bed and flame zones, it is not possible to discuss possible reasons for discrepancies. The predictions of CO2 emissions were reasonably concurrent, and the CO emissions were in the same magnitude range as the experiments, which can be considered a good result when considering the simplicity of the homogeneous gas reactions. Because the simulations were performed in a transient regime, the evolution of the different variables can also be analysed. Fig. 15 shows the temporal evolution of the CO and CO2 emissions and their comparison with the experimental data during the E1 test. It should be noted that the comparison is not straightforward, as it is well known that the dynamic response of gas analysers has a significant effect on instantaneous readings. In general, they can introduce a delay and filter effect on the actual value of the emissions at the chimney. This has not been addressed in this work because it is only shown for comparison purposes. In addition, the feeding system was not synchronised between the simulations and the experiments, which impedes any possible coincidences between the phases of the oscillations. In any event, the predictions of the CO2 emissions showed a similar trend in terms of

oscillation amplitude. However, the predictions of CO emissions were more regular than the measurements, which were very low in both cases. 6. Conclusions This work simulated a commercial biomass boiler by implementing CFD classical techniques for the gas phase and a comprehensive model to account for the main solid variables and their interaction with the gas phase. Several simulations were performed in a commercial boiler to study the behaviour of the model and to compare its results with those obtained from experimental tests under similar conditions. The boiler was simulated with a fully embedded transient bed model composed of several sub-models, some of which had been previously presented by the authors. Several new features were incorporated to enable the transient modelling of a complete boiler. The major additions were the introduction of an advective motion for the fuel in the bed, which was necessary to represent the fuel feeding, the simulation of the solid and volatile ash content of the fuel, particle shrinkage during the devolatilisation and char combustion stages and the blowing effect of the outgoing gases in the convective heat and mass exchange between the reacting particles and the surrounding gases. The transient simulation of the boiler allowed a detailed visualising of the evolution of the variables for both solid and gas phases. Specific variables such as the heat transferred to the water and the temperature/composition of the flue gases were used to compare the simulation with the experimental results. These results showed reasonably good predictions for the main emissions and the heat exchanged. The temperatures of the flue gases were concurrent in high load cases but showed significant discrepancies for partial loads, which may be due to the uncertainty in modelling the thermal effect of fouling walls. The dynamic pattern of the boiler emissions was analysed through the results of the dynamic model, although further research is required on these methods to experimentally characterise those instantaneous emissions. Acknowledgments The authors wish to acknowledge financial support from Project EN2012-36545 of the Ministry of Economy of Spain. The authors

680

M.A. Gómez et al. / Energy Conversion and Management 101 (2015) 666–680

also acknowledge the valuable suggestions by J.L. Toste de Azevedo. This study is dedicated to his memory.

References [1] Karim MR, Naser J. Progress in numerical modelling of packed bed biomass combustion. In: 19th Australasian fluid mechanics conference, Melbourne, Australia; 2014. [2] Porteiro J, Granada E, Collazo J, Patiño D, Morán JC. Numerical modeling of a biomass pellet domestic boiler. Energy Fuels 2009;23:1067–75. [3] Collazo J, Porteiro J, Míguez JL, Granada E, Gómez MA. Numerical simulation of a small-scale biomass boiler. Energy Convers Manage 2012;64:87–96. [4] Gómez MA, Comesaña R, Álvarez Feijoo MA, Eguía P. Simulation of the effect of water temperature on domestic biomass boiler performance. Energies 2012;5:1044–61. [5] van der Lans RP, Pedersen LT, Jensen A, Glarborg P, Dam-Johansen K. Modelling and experiments of straw combustion in a grate furnace. Biomass Bioenergy 2000;19:199–208. [6] Kær SK. Numerical modelling of a straw-fired grate boiler. Fuel 2004;83:1183–90. [7] Yin C, Rosendahl L, Kær SK, Clausen S, Hvid SL, Hiller T. Mathematical modeling and experimental study of biomass combustion in a thermal 108 MW gratefired boiler. Energy Fuels 2008;22:1380–90. [8] Zhou H, Jensen AD, Glarborg P, Jensen PA, Kavaliauskas A. Numerical modeling of straw combustion in a fixed bed. Fuel 2005;84:389–403. [9] Yang YB, Goh YR, Zakaria R, Nasserzadeh V, Swithenbank J. Mathematical modelling of MSW incineration on a travelling bed. Waste Manage 2002;22:369–80. [10] Shin D, Choi S. The combustion of simulated waste particles in a fixed bed. Combust Flame 2000;121:167–80. [11] Cooper J, Hallett WLH. A numerical model for packed bed combustion of char particles. Chem Eng Sci 2000;55:4451–60. [12] Thunman H, Leckner B. Co-current and counter-current fixed bed combustion of biofuel—a comparison. Fuel 2003;82:275–83. [13] Peters B, Bruch C. Drying and pyrolysis of wood particles: experiments and simulation. J Anal Appl Pyrol 2003;70:233–50. [14] Yang YB, Ryu C, Khor A, Yates NE, Sharifi VN, Swithenbank J. Fuel size effect on pinewood combustion in a packed bed. Fuel 2005;84:2026–38. [15] Yang YB, Sharifi VN, Swithenbank J. Numerical simulation of the burning characteristics of thermally-thick biomass fuels in packed beds. Trans IChemE B Process Saf Environ Protect 2005;83:1–11. [16] Yang YB, Sharifi VN, Swithenbank. Sub-stoichiometric conversion of biomass and solid wastes to energy in packed beds. AIChE J 2006;52:809–17. [17] Yang YB, Newman R, Sharifi V, Swithenbank J, Ariss J. Mathematical modelling of straw combustion in a 38 MWe power plant furnace and effect of operating conditions. Fuel 2007;86:129–42. [18] Collazo J, Porteiro J, Patiño D, Granada E. Numerical modeling of the combustion of densified wood under fixed-bed conditions. Fuel 2012;93:149–59. [19] Gómez MA, Porteiro J, Patiño D, Míguez JL. CFD modelling of thermal conversion and packed bed compaction in biomass combustion. Fuel 2014;117:716–32. [20] Thunman H, Leckner B, Niklasson F, Johnsson F. Combustion of wood particles – a particle model for Eulerian calculations. Combust Flame 2002;129:30–46. [21] Mehrabian R, Scharler R, Obernberger I. Effects of pyrolysis conditions on the heating rate in biomass particles and applicability of TGA kinetic parameters in particle thermal conversion modelling. Fuel 2012;93:567–75. [22] Mehrabian R, Zahirovic S, Scharler R, Obernberger I, Kleditzsch S, Wirtz S, et al. A CFD model for thermal conversion of thermally thick biomass particles. Fuel Process Technol 2012;95:96–108. [23] Ström H, Thunman H. CFD simulations of biofuel bed conversion: a submodel for the drying and devolatilization of thermally thick wood particles. Combust Flame 2013;160:417–31. [24] Chaney J, Liu H, Li J. An overview of CFD modelling of small-scale fixed-bed biomass pellet boilers with preliminary results from a simplified approach. Energy Convers Manage 2012;63:149–56. [25] Zhang X, Chen Q, Bradford R, Sharifi V, Swithenbank J. Experimental investigation and mathematical modelling of wood combustion in a moving grate boiler. Fuel Process Technol 2010;91:1491–9.

[26] Eskilsson D, Ronnback M, Samuelsson J, Tullin C. Optimisation of efficiency and emissions in pellet burners. Biomass Bioenergy 2004;27:541–6. [27] Mehrabian R, Stangl S, Scharler R, Obernberger I, Weissinger A. CFD simulation of biomass grate furnaces with a comprehensive 3D packed bed model. In: Proceedings of the 25th German flame day, Karlsruhe, Germany, September; 2011. [28] Miltner M, Makaruk A, Harasek M, Friedl A. CFD-modelling for the combustion of solid baled biomass. In: Fifth international conference on CFD in the process industries, Melbourne, Australia; 2006. [29] Bhuiyan AA, Naser J. Computational modelling of co-firing of biomass with coal under oxy-fuel condition in a small scale furnace. Fuel 2015;143:455–66. [30] Bhuiyan AA, Naser J. Numerical modelling of oxy fuel combustion, the effect of radiative and convective heat transfer and burnout. Fuel 2015;139:268–84. [31] Dieguez-Alonso A, Anca-Couce A, Zobel N. On-line tar characterization from pyrolysis of wood particles in a technical-scale fixed-bed reactor by applying Laser-Induced Fluorescence (LIF). J Anal Appl Pyrol 2013;102:33–46. [32] Porteiro J, Patiño D, Collazo J, Granada E, Moran J, Miguez JL. Experimental analysis of the ignition front propagation of several biomass fuels in a fixedbed combustor. Fuel 2010;89:26–35. [33] Porteiro J, Patiño D, Miguez JL, Granada E, Moran J, Collazo J. Study of the reaction front thickness in a counter-current fixed-bed combustor of a pelletised biomass. Combust Flame 2012;159:1296–302. [34] Girgis E, Hallett WLH. Wood combustion in an overfeed packed bed, including detailed measurements within the bed. Energy Fuels 2010;2(4):1584–91. [35] Di Blasi C. Modeling wood gasification in a countercurrent fixed-bed reactor. AIChE J 2004;50:2306–19. [36] Zobel N, Anca-Couce A. Slow pyrolysis of wood particles: characterization of volatiles by laser-induced fluorescence. Proc Combust Inst 2013;34:2355–62. [37] Wagenaar BM, Prins W, van Swaaij WPM. Flash pyrolysis kinetics of pine wood. Fuel Process Technol 1993;36:291–8. [38] Bryden KM, Ragland KW. Modeling of a deep, fixed bed combustor. Energy Fuels 1996;10:269–75. [39] Hagge M, Bryden K. Modeling the impact of shrinkage on the pyrolysis of dry biomass. Chem Eng Sci 2002;57:2811–23. [40] Kær SK. Straw combustion on slow-moving grates—a comparison of model predictions with experimental data. Biomass Bioenergy 2005;28:307–20. [41] Christensen KA, Livbjerg H. Aerosol Sci Technol 1996;25:185–99. [42] Baxter L, Miles Jr TR, Milne T, Bryers RW. Fuel Process Technol 1998;54:47–78. [43] Jensen PA, Stenholm M, Hald P. Deposition investigation in straw-fired boilers. Energy Fuel 1997;11:1048–55. [44] Kær SK. Numerical investigation of deposit formation in straw-fired boilers. PhD thesis, Institute of Energy Technology, Aalborg University, Denmark; 2001. ISBN: 87-89179-39-0. [45] Incropera FP, De Witt DP. Fundamentals of heat transfer. Wiley; 1981. [46] Collazo J, Porteiro J, Patiño D, Míguez JL, Granada E, Morán JC. Simulation and experimental validation of a methanol burner. Fuel 2009;88:326–34. [47] Yin C, Rosendahl L, Clausen S, Hvid SL. Characterizing and modeling of an 88 MW grate-fired boiler burning wheat straw: experience and lessons. Energy 2012;4:473–82. [48] Gómez MA, Patiño D, Comesaña R, Porteiro J, Álvarez Feijoo MA, Míguez JL. CFD simulation of a solar radiation absorber. Int J Heat Mass Transf 2013;57:231–40. [49] Ansys Fluent. 12.0 Theory guide; 2009. [50] Thunman H, Niklasson F, Johnsson F, Leckner B. Composition of volatile gases and thermochemical properties of wood for modeling of fixed or fluidized beds. Energy Fuel 2001;15:1488–97. [51] Neves D, Thunman H, Matos A, Tarelho L, Gómez-Barea A. Characterization and prediction of biomass pyrolysis products. Prog Energy Combust 2011;37:611–30. [52] Andersen J, Jensen PA, Meyer KE, Hvid SL, Glarborg P. Experimental and numerical investigation of gas-phase freeboard combustion. Part 1: main combustion process. Energy Fuels 2009;23:5773–82. [53] Jones WP, Lindstedt RP. Global reaction schemes for hydrocarbon combustion. Combust Flame 1988;73:233–49. [54] Ouedraogo A, Mulligan JC, Cleland J. A quasi-steady shrinkage core analysis of wood combustion. Combust Flame 1998;114:1–12. [55] Muzaferija S, Peric M, Sames P, Schellin T. A two-fluid Navier-Stokes solver to simulate water entry. In: Proc 22nd symposium on naval hydrodynamics, Washington DC; 1998. [56] Richards G, Slater P, Harb J. Simulation of ash deposit growth in a pulverized coal-fired pilot scale reactor. Energy Fuels 1993;7:774–81.