A CFD-based comparative analysis of drying in various single biomass particles

Applied Thermal Engineering 128 (2018) 1062–1073

Contents lists available at ScienceDirect

Applied Thermal Engineering journal homepage: www.elsevier.com/locate/apthermeng

Research Paper

A CFD-based comparative analysis of drying in various single biomass particles Hassan Khodaei a,c,⇑, Guan H. Yeoh b, Ferdinando Guzzomi a, Jacobo Porteiro d a

School of Engineering, Edith Cowan University, Joondalup, WA 6027, Australia School of Mechanical and Manufacturing Engineering, University of New South Wales, NSW, Sydney, Australia c Innovative Reduction Strategies Inc. (IRSI), Edmonton, AB, Canada d E.TS. Ingenieros Industriales universidade de Vigo, Vigo, Pontavedra, Spain b

h i g h l i g h t s
The geometrical shape of biomass and the external boundaries are investigated.
Reducing moisture content decreases drying time about 83%.
Non-uniform distribution of temperature increased drying time significantly.

a r t i c l e

i n f o

Article history: Received 2 June 2017 Revised 23 August 2017 Accepted 13 September 2017 Available online 20 September 2017 Keywords: Single biomass particles Moisture evaporation CFD modelling

a b s t r a c t The understanding of the behaviour of biomass particles is important, therefore, modelling different subprocesses of biomass thermal conversion is derived during these years. This paper addresses a comparative CFD based analysis of different drying models. Several sub-models are simulated to investigate the evaporation process of different geometries based on standard densified wood pellets. In order to predict the transient evolution (moisture to dry wood) of the wood particles, transport equations (energy and moisture evaporation) are solved considering the reaction heat loss, effective thermal conductivity, specific heat capacity, and radiative and convective heat transfer. These models are compared with the previous experimental and numerical works. The Heat Sink model demonstrates the closest agreement with the reported data based on a cylindrical shape of biomass particle. This model is further analyzed by increasing the moisture content and decreasing the surface area exposed to radiative and convective heat taking into consideration the particle density, effective thermal conductivity and specific heat capacity. The results show a remarkable decrease in drying time when the particle is fully exposed to external radiative and convective heat with the lowest moisture content. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction Biomass fuel is regarded as an important renewable energy resource. Owing to the importance of this resource, a better understanding of its thermal conversion behaviour is needed in order to improve process efficiency and reduce pollutant emissions. Coupling the heat and mass transfer with the chemical processes involved in both packed and fluidised bed requires the development of a sophisticated CFD model. The complexity of such models is due to intra-particle reaction processes such as drying, devolatilization, char burning and oxidisation, interconnection conductivity and mass transfer between particles and its surrounds ⇑ Corresponding author. E-mail address: [email protected] (H. Khodaei). https://doi.org/10.1016/j.applthermaleng.2017.09.070 1359-4311/Ó 2017 Elsevier Ltd. All rights reserved.

[1]. Convective heat exchange in packed-beds is not constant due to different counter-current air flow-rates and is calculated based on experimental correlation between the Nusselt and Sherwood numbers [2,3]. In addition, radiation models, geometry configurations, anisotropic shape of the fuel, phenomena like channelling and shrinkage and particles increase the level of complexity in CFD modelling [4]. This has resulted in focused attention on mathematical modelling and CFD analysis of single biomass particles. Biomass contains water in three forms: bound water, free water and water vapour. Bound water is found in plant cell walls. The wood particle is occupied by free water, water vapour is found in cell lumina and compared to the other forms is negligible at atmospheric temperature and humidity [5]. Most manufacturers establish the limit for moisture content in the fuels due to the influence of moisture in the design and selection of burners. In the wood pel-

H. Khodaei et al. / Applied Thermal Engineering 128 (2018) 1062–1073

1063

Nomenclature Av As cpm cp:p Ev h kv lpore LHmoist v qcon g R r dry Ss T t Dt

pre-exponential factor [s1] surface area [m2] moisture heat capacity [J/kg K] Particle heat capacity [J/kg K] activation energy [J/mol] convective heat transfer coefficient [W/m2 K] kinetic reaction rate [s1] pore diameter [m] evaporation enthalpy [J/kg] convective heat (gas phase) [W/m2] universal gas constant [J/mol K] evaporation rate [kg/m3s] depending source term temperature [K] time [s] current time step [s]

Greek letters

qmoisture ; qm moisture density [kg/m3] keff effective thermal conductivity [W/m K] krad radiative thermal conductivity [W/m K] rð5:67  108 Þ Stefan Boltzmann [W/m2 K4] x emissivity [–] e solid mass fraction [–]

let industry, a combination of low moisture content (8–10%) and high densified pellets (typically greater than 1000 kg/m3) offer fuels with higher bulk energy densities. Moisture content influences the combustion behaviour, particularly pyrolysis and the adiabatic temperature of combustion and the volume of emitted gases. Conversely, longer resistance times for drying leads to a later occurrence of other sub-processes of thermal conversion like pyrolysis and char combustion [6]. From an energetic stand point, Minkova et al. [7] demonstrated that the liberated steam (water vapour) has a strong effect on the yield and properties of pyrolysis products from biomass fuels. In indirect bio-char reactors, increasing moisture contents and non-uniformity of temperature reduce the quality of bio-char and leads to lower efficient syngas combustion and higher emissions [8]. Different experimental studies are carried out to study thermal degradation of single biomass particles. The temperature and density profiles, outflowing water and volatiles analysis during pyrolysis of large biomass particles were performed by Chan et al. [9]. Alves and Figueiredo [10] investigated pyrolysis of cylindrical, wet and dried woods. Tran and White [11] calculated the burning rate in different thick wood particles by determining mass loss, average heat release and charring rate in a heat-release calorimeter. Gronli [12] investigated the pyrolysis of different wet biomass particles taking into account temperature profiles and chemical composition in order to describe thermal degradation in single biomass particles during drying and pyrolysis. Lu et al. [13] studied drying, devolatilization and char oxidisation of different spherical and cylindrical biomass woods. Their experimental facility provided mass loss and particle surface and centre temperatures as a function of time. However, widely varying property (density, specific heat, thermal conductivity, porosity, permeability) among different woods, expensive equipment and the problem of installation (uncertainty in thermal conductivity, specific heat capacity and combustion progress) introduce uncertainties in the experimental studies. Therefore, modelling is useful not only from an

Subscripts conv convection cond conduction dry drying evp evaporation eff eeffective g gas m moisture p particle reac reaction rad radiation s solid HS Heat Sink AR Arrhenius Cy cylindrical geometry Cb cubic geometry Sp spherical geometry ⁄ mixed radiative temperature T = 800 K and natural convection is imposed on the half surface of the particles while the rest of the surfaces assume to be cooled by natural convection

economical aspect but also to provide more data than experimental results. The main advantage of single particle modelling is easy of implementation and assessment of essential physics. Saastamoinen and Aho [14] proposed a mathematical model of simultaneous drying and pyrolysis of single particle biomass pellets. They used temperature dependence functions based on experimental analysis and numerical approximations to express the phase change in evaporation and pyrolysis processes. They also determined that the Arrhenius model for larger particles is not accurate due to different wood varieties. The immediate outflow hypothesises used in the model was applied by Porteiro, Collazo et al. [15,16]. Based on this hypothesis, an outflow of gas species occurs in the reaction zone and moisture evaporated in the solid phase are incorporates into the gas phase domain. In fact, the time required for gases to travel through the outer layers of the particles is smaller than the simulation time step, therefore recondensation and moisture transfer in the liquid phase were not considered. This hypothesis is reinforced in the experimental study proposed by Gronli [12]. Porteiro et al. [17,15] investigated a one-dimensional model of a single wood particle incorporates sub-processes of thermal conversion. Mass loss and internal temperature were compared with the previous experimental studies. Yang et al. [18] applied a mathematical and experimental approaches to investigate the combustion characteristics of a single biomass particle for different sizes, taking into account volume shrinkage. The results demonstrated that assuming constant temperature for large particles appeared to be inadequate and more detail CFD modelling should be applied for those particles. Single particle studies have also been considered by several researchers [9,19– 22] to introduce intra-particle and interactions between single particles in the packed-beds. Peters [23] employed the Thiele modulus to determine whether particle combustion was dominated by reaction or diffusion and then classified the combustion regime based on the Damkohler number (indicating where the bed behaviour is homogeneous) and the Thiele modules. The use of drying num-

1064

H. Khodaei et al. / Applied Thermal Engineering 128 (2018) 1062–1073

ber (Dr), which is a qualitative measure of relative time for drying and devolatilisation, was proposed by Thunman et al. [24] in order to investigate separation of drying and devolatilisation processes in solid biomass fuels. In small particles (Dr  1), the process of drying and devolatilisation can be considered separately. Johansson et al. [25] studied the influence of intra-particle gradients on the internal rate of drying and devolatilisation by introducing a coupled single and packed-bed combustion model. The results showed a reasonable prediction of species and temperature profiles for thermally thin and thick particles. Thermally thick simulations of thermal conversion for different single particles was proposed by Mehrabian et al. [26] in order to describe the influence of particle interaction and movement within packed beds. The Heat Sink model was used for drying. Although the abovementioned studies considered drying in biomass combustion, there is a lack of comparative studies of drying for single biomass particles. 2. Model description The drying process is examined through three different approaches in the literature, namely the Equilibrium, Heat Sink and Arrhenius models [27]. In the Equilibrium model, liquid water is in equilibrium with the local water vapour (in the wood pores), therefore the partial pressure of water vapour is fixed by the saturation pressure. This assumption is typical for with lowtemperature drying models. This model is dependent on both heat and mass transfer within the wood pores [28–30]. In the Heat Sink model, evaporation is approximated by an energy balance that employs a given temperature for evaporation [15]. This approach is limited by heat transport. Whereas, the Arrhenius approach includes the reaction and the thermal time scales [31]. Bellais et al. [32] studied three different drying models under pyrolysis. The results of this study indicated that surface temperature and global drying rate in the first-order reaction rate (Arrhenius approach) is closer to the experimental results. Lu et al. [13] performed a comprehensive numerical simulation of different subprocesses of thermal conversion and demonstrated the influence of wood particle composition and temperature gradient in all combustion sub-processes. The drying model in this study was a combination of the equilibrium model for free water in the particles and Arrhenius model for bound water in the particles. In recent work, Gomez et al. [33] assumed a sub-grid scaled approach during the drying process. In this approach, it was assumed that the drying front occurs at an interface between the moist wood region and the dry wood layer (outer surface). The difference between the heat transferred from the dry wood to the moist wood zones was assumed to be the heat employed moisture evaporation. This methodology was applied to thermally thick particles. This work investigates CFD based comparative study of two different drying models: The Heat Sink and Arrhenius models. The focus of the study is on considering three different shapes of biomass fuels with the same mass and similar properties, while particle evolution during the drying process. The most reliable drying model is then recognized based on the comparison of the current simulation with the past experimental and numerical results. Considering the effect of geometry (shape and, external surface heat flux) and moisture content on drying time, thermal conductivity, specific heat capacity, density and solid particle temperature and evaporation rate are investigated. 2.1. Arrhenius model In the Arrhenius model, a first order kinetic reaction rate has been proposed to model the evaporation rate in single and packed-bed biomass particles. Based on this model, the water con-

version from liquid to gas state is assumed to be a first-order reaction rate: kv

Moisture ! Water Vapour

ð1Þ

K m ¼ Am expðEm =RT s Þ½1=s

ð2Þ

rdry ¼ kv qmoisture ½kg=m3 s

ð3Þ

Depending on the size and type of wood, different types of kinetic reaction rates were identified [9,20,31,34]. Chan et al. [9] proposed Av = 5.13  106 s1 and Ev = 88 kJ mol1, which results in a drying plateau at higher temperatures as expected for thick particles assuming P is 105 pa (1 atm) [9]. To provide a drying plateau around 100 °C, these drying rate parameters introduced by Chan et al. [9] modified by Bryden et al. [31] and are adopted in this study (Av = 5.13  1010 s1 and Ev = 88 kJmol1). They influence of moisture in combustion for large timber wood was investigated by Yuen et al. [35], used an Arrhenius expression proposed by Atreya et al. [36]. Similar to the previous investigations, a higher temperature plateau was observed. Peters et al. [37] used free and bound water kinetic reaction rate as heterogeneous reactions proposed by Chan et al. [9] and Krieger-Brockett et al. [38], however, it seems that the drying time increased dramatically compared to other drying models. Miltner et al. [39] proposed an additional temperature dependence term for the solid biomass temperature; however, this approach is not valid for solid temperatures higher than 475 K (202 °C). Bryden et al. [31] proposed a recondensation factor for the drying process to compensate for cooling during the process. This assumption is typically applied to large and thick biomass fuels in addition that, the evolution of pressure during combustion was considered. Since the objective of this study is not directly related to combustion and basically focus on fixed lower external temperatures (applicable in drying and pyrolysis technologies), the evolution of pressure and recondensation effects during thermal conversion are not included in these models. There is reasonably close agreement between the simulation results and previous experimental and numerical investigations as addressed in Section 7.1. 2.2. Heat Sink model In the Heat Sink model, the drying process approximates water evaporation based on an energy balance and a given evaporation temperature (saturation, boiling). The drying zone is assumed infinitely thin region. This drying zone acts as a Heat Sink and the drying rate is controlled by heat transfer. In Heat Sink model, drying is assumed to occur at a boiling temperature corresponding to the external pressure. Hence both the drying front and the particle surface are at external pressure and the water vapour cannot flow by convection out of the particle. In this approach, a threshold temperature is considered for evaporation; therefore, the evaporation process consumes heat above the evaporation temperature (373 K) without distinguishing bound and free water in the solid particles [18,20,29]. This model was widely used for modelling packed-beds when a single particle was resolved into the bed or defined as a spherical porous medium [1,20]. Peters et al. [20] proposed the density and thermal capacity of dry wood in the Heat Sink model, however, the analysis of this model highlighted instabilities in drying rate contours when the temperature reaches the evaporation temperature. Porteiro et al. [17] employed the thermal capacity of the particle to convert the excess temperature generated in the whole cell into the heat of evaporation, therefore, after the evaporation has taken place within a time-step, the whole cell should be cooled to the evaporation temperature. This produces a steep evaporation rate, but can generate numerical instability. In

1065

H. Khodaei et al. / Applied Thermal Engineering 128 (2018) 1062–1073

the Heat Sink model, drying and heating overlap in the outer surface of the particle. Gómez et al. [1] and Collazo et al. [16] assumed a parameter (s) which represented the heat fraction for drying. They used s = 0.5 based on their numerical studies. To avoid potential instability, Yang et al. [18] assumed X m where m is initial moisture content in the biomass fuel. The density and thermal capacity of moisture were used as indicated below:

8 < Xm ðT s T ev ap Þqm cpm rdry¼

:0

LHmoist :dt

; if

if T S P T ev p [ X m > 0

if T S < T ev p

ð4Þ

½kg=m3 s

3. Main assumptions The assumptions below are valid for both Arrhenius and Heat Sink models. 1. The energy and mass exchange to and from the single particle are through the boundary layers. 2. Gas and solid phases have the same temperature inside the particle when using equilibrium model. 3. The solid particle consists of moisture and dry wood. 4. The porosity of a wood particle is considered and the movement of solid phase during evaporation is ignored due to minor volume shrinkage during drying according to Johansson et al. [25]. 5. Uniform radiative and convective heat transfer exists between the particle external surface and surrounding environment. 6. A threshold temperature is assumed to start evaporation with the Heat Sink approach; however, this limitation is not considered in the Arrhenius approach. 7. An instantaneous outflow of gas species from the solid into the gas is assumed.

Table 2 Source terms for energy equation and mass loss. Sreac ¼ r dry :e:LHmoist s

(11)

[1]

v ¼ qconv ¼ h Scon conv As ðT g  T s Þ s g

(12)

[17,18]

Srad ¼ rxðT 4env  T 4s Þ s xloss ¼ xm;init  xm

(13)

[17,18]

(14)

[17]

the multiplication of solid mass fraction and drying rate, which is discussed in Section 2. The conversion of solid wood to gas is considered by removing the mass fraction of water vapour released from the physical domains during the drying process as mentioned in (Eq. (14)). 5. Fuel geometries and properties Three different fuel shapes are considered, namely cylindrical, spherical and cuboid. The fuel shape is selected according to the standard size of a cylindrical wood pellet [40]. Cuboid shape is chosen based on the mass and the length of standard cylindrical wood pellet while spherical matches only with the mass of the cylindrical wood pellet. The cuboid shape is selected due to its application in some home and industrial biomass burners or slow pyrolysis systems. The cylindrical shape is chosen due to the importance of cylindrical wood pellets. The spherical shape is considered to compare with the previous studies considering the spherical shapes of wood particles. Effective thermal conductivity is selected based on Johansson et al. [41] and specific heat capacity of the particles is identified based on specific heat capacity of wet wood proposed by Regland et al. [42]. Table 3 depicts the specification of fuels used in this study. 6. Simulation

4. Solid phase governing equations

6.1. Mesh and boundary conditions

In this study, user defined functions (UDFs) are employed to characterise the drying process in a research version of AnsysFluent 15. Two user defined scalars are identified to represent: solid temperature (T s ) and mass fraction of moisture (X m ) in both Heat Sink and Arrhenius models. Table 1 shows the transport equations (Eqs. (5 and 6)), the relation between the particle density, moisture and dry wood (Eqs. (7–9)) and the energy equation source terms (Eq. (10)). The left-hand side of the transport equations (Eqs. (5 and 6)) represents the unsteady term. The diffusive term of the energy equation (Eq. (5)) describes the effective thermal conductivity as a weighted sum of thermal conductivities of non-gas phase components (moisture and dry wood). The second term on the right-hand side of Eq. (5) represents the source terms. These source terms are: the energy of reaction due to drying process inside the particle (Eq. (11)), the radiative and convective energy exchanged through the particle’s surfaces respectively (see Table 2: Eqs. (12 and 13)). The right-hand side term in the moisture fraction (Eq. (5)) denotes

Twelve different simulations (six for Arrhenius and six for Heat Sink models) are examined. Appropriate meshes are generated, whilst considering the sensitivity of these meshes. More grid nodal points are concentrated near the walls to resolve the exposed boundary by external radiation as shown in Fig. 1. The mesh density independency is validated in different geometries. No significant change is observed in Sp, Sp⁄, Cb, Cb⁄ and Cy⁄, however, Cy (2D axisymmetric model) seems to be dependent on the mesh density. Table 4 shows the range of mesh independency tests in this geometry, moreover, different mesh densities (number of elements) versus solid wood temperature (Ts) and moisture fraction (Xmoist) for Cy are addressed in Fig. 2 The simulation results are almost similar in the mesh domain between 960,00 and 320,000 elements. The particles are assumed as standalone particles and Define Profile macros are used to define a custom boundary profile which covers heat transfer coefficients, heat flux, temperature and external emissivity. Radiative and convective heat transfers for laminar flow are identified on the external surfaces of the particles [18]. Notably, the result of this study could be used for modelling of biochar conversion in indirect fast pyrolysis techniques, therefore, assuming a stagnant (or laminar) flow around the particle and constant temperature (the desired temperature for fast pyrolysis is approximately 773 K) are valid as addressed in Eq. (18).

Table 1 Equation for solid phase variables. Energy equation Moisture fraction Total particle density

Energy equation source

@ð2qp Cp:pT s Þ ¼ rðkeff  rT S Þ þ @t 2qp @ðX m Þ ¼2 :r dry @t  p ¼ X m p;0 þ X drywood p;0  m ¼ X m p;0

Ss

(5)

[1]

(6)

[18]

q q q q q qdrywood ¼ X drywood qp;0

(7)

v þ Srad þ Scon Ss ¼ Sreac s s s

(10)

(8) (9) [1]

keff

  @T ¼ erad rðT 41  T 4x¼R=L Þ þ hðT 41  T 4x¼L=R Þ @x x¼L=R

ð18Þ

1066

H. Khodaei et al. / Applied Thermal Engineering 128 (2018) 1062–1073

Table 3 Fuel physical and thermal properties and dimensions. Fuel properties Moisture [wt.%] Dry wood [wt.%] Solid fraction [e] qp [kg/m3]

8 92 0.45 1240

lpore [m]

5  105

Fuel geometries

Fuel dimensions (mm)

Cylinder Sphere Cuboid

r = 4, h = 15 r = 5.63 b = 7, h = 15

W [–] 0.84 Thermal properties cpm [J/kg K] Cp.p [J/kg K] 103

Dn

½ð0:1031þ0:003867T s Þþ4:19xmoisture  ð1þxmoisture Þ

o

þ fð0:02355T S  1:32  xmoisture  6:191Þ  xmoisture g

4180 (15)

E

keff ¼ xm kmoisture þ xdrywood kdrywood þ krad [W/m K]

(16) (17)

krad ¼ 13:5rT 3s lpore =x [W/m K] LHmoist [J/kg]

2:25  106 4180 12

cpm [J/kg K] hconv [W/m2 K] wt%: percentage of weight

Fig. 1. Mesh configuration in different geometries: (a): 3D cylinder (Cy*); (b): 2D- Cuboid (Cb, Cb*); (c): 2D Sphere (Sp, Sp*); (d): 2D-Axsysemetric cylinder (Cy).

Table 4 Mesh independency results in 2D axisymmetric cylindrical biomass particle. Number of elements

Bias Number

Max aspect ratio

Axial

Radial

24000



40

96000



80

320000



158

Number of division

Mesh independency results

Radial

Axial

3:74  101

80

300

1:81  101

160

600

Independent

3:73  101

320

1000

Independent

Dependent

Note: mesh dependent range: 24  103 < Number of elements < 96  103

A mixed radiative temperature T = 800 K and natural convection are imposed on the half surface of particles while the rest of the surfaces are assumed to be cooled by natural convection or exposed by a similar level of external radiation (Table 5). A com-

mon value of convective heat transfer coefficient for natural convection, hconv = 12 [W/m2 K] is assumed based on Incropera [43]. The emissivity of the wall is assumed to be 0.95 based on Yang et al. [18]. To start calculation, we assumed the outer boundary

1067

H. Khodaei et al. / Applied Thermal Engineering 128 (2018) 1062–1073

Fig. 2. Number of elements (El) versus solid temperature and moisture fraction in Cy.

Table 5 Geometry and mesh identifications. Abbrv

Geometries

Number of Geometries

Model type 2Db

Sp Spc Cb Cbc Cy Cyc

Sphere Spherec Cuboid Cuboidc Cylinder Cylinderc

2 2 2 2 2 2

2D-axia

Element shape

Inflation

Number of elements

Quadrilateral Quadrilateral Quadrilateral Quadrilateral Quadrilateral Tetrahedral

U U  U U U

1090 1090 10,800 11,468 96,000 21,139

3D

U U U U U

Keys: a axi: Axisymmetric model. b 2D, 3D: two and three-dimensional geometries. c Geometry is divided into two parts in order to define different boundary conditions. Half of the particle surfaces are exposed to mixed 800 K radiative temperature and natural convection. Another half is exposed to natural convection and ambient temperature. For other geometries, all surfaces are exposed to a mixed 800 K radiative and natural convective.

of the domain is completely dried fast therefore Xmoist at the outer layer of wall is considered zero. In solution initialisation, the following values are assumed: Ts = 300 K and Xmoisture, initial = 0.08, 0.20, 030%. Furthermore, in all tests, the time step is set to 0.005 s and the number of time steps for calculation is assumed 80,000. A second order upwind schematic was selected for the user defined scalar Ts and Xmoist. Table 5 shows the specification of geometries presented in this study. Mixed thermal conditions (both convection and radiation) between the walls and surrounding area around the particles are incorporated into the UDFs. In addition, a diffusivity macro is used to evaluate the effective thermal conductivity of wood. Ten user defined memories are identified to track the variables during the evaporation process, including: solid temperature (Ts), moisture fraction (Xmoist), density of particle (qp Þ, moisture (qm Þ, specific heat capacity (C p;p ), effective thermal conductivity (keff ), mass loss (mloss ), drying rate (r dry ), Xloss and dT/dt. Fig. 3 illustrates a schematic of the solution algorithm. 7. Results and discussion 7.1. Experimental and numerical validation The discussed modelling methodology was developed for different types and sizes of wood particles. However, assessment of a

single model should be considered both numerically and experimentally. To validate drying models in the single particles, a separation between two sub-processes of thermal conversion (drying, devolatilization) could be considered [24]. In this regard, three set of validations were conducted to compare the behaviour of the models with the existing experimental data as well as the other numerical models. The first validation case is addressed in Fig. 4(a) which corresponds to the research study by Peters et al. [20]. Evaluation of these results indicates that the Heat Sink model is closer to the Heat Sink model proposed by Peters. The Heat Sink model in this study has been validated in three different studies: experiments on spherical particles of wet fir wood with 67% moisture content and wall temperature 743 K [44], the model with Arrhenius expression rates proposed by Chan et al. [9] (AR-2) and the Arrhenius model proposed by Peters [21]. It should be noted that Arrhenius model proposed in this study (the orange line) is closer to the experimental study [10] and numerical simulation proposed by Peters [20] than Arrhenius numerical model proposed by Peters [20]. It could be due to choosing the kinetic reaction rates (Av = 5.13  1010 s1 and Ev = 88 kJ mol1) for drying in the current investigation. In overall, the Arrhenius models have slower response than the Heat Sink model. The second validation case is addressed in Fig. 4(b), which corresponds to Bryden et al. [31] and experimental study of a 40 mm

1068

H. Khodaei et al. / Applied Thermal Engineering 128 (2018) 1062–1073

Fig. 3. Drying solution algorithm for Heat Sink and Arrhenius model.

[Current]

1.2

[19] Peters-AR [19] HS [Current] AR [Current]

[Current]

Peters-HS

m H2O /m H2O,0 [-]

1

[30] [30]

0.8

0.6

0.4

0.2

0 0

100

200

300

400

500

600

Time [s] Fig. 4. (a): 4 mm spherical fir wood Twall = 743 K with 67% moisture and (b): 40 mm red oak cube exposed to a 38 kW/m2 radiative heat flux on one surface respectively. The shaded area in (b) is the experimental domain proposed by Tran et al. [11]. The red dotted line is mass fraction of moisture related to x = 36 mm HS (red solid line) model. (HS*): is mass fraction in Heat Sink (current study). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

thick red oak sample with 8% moisture exposed on one surface to a 38 kW/m2 radiative heat flux by Tran and White [11]. The Heat Sink model in this study agrees well with the experimental study. Comparing to Bryden model, the current Arrhenius model appears to be closer to the numerical study by Bryden [31]. The third validation case concerns a cubic poplar wood with 2cm half thickness and 30% moisture exposed to 1500 K background temperature [45]. The evaporation for both Heat Sink and Arrhenius models compare well with the results. It is evident that the evaporation time in the Heat Sink model is less than the Arrhenius model modified by Bryden et al. [45]. The kinetic reaction proposed by Chan et al. [9] and Peters et al. [20] are also examined, however, the results were far from the Heat Sink and Arrhenius models adopted in this study. Due to considering all sub-

processes of combustion (drying, pyrolysis and gasification) in the Bryden model, a steeper profile in their model appears after completion of the evaporation compared to the current numerical study as illustrated in Fig. 5. 7.2. Modelling results Arrhenius and Heat Sink models are applied to six geometries considering 800 K radiative temperature on all and half of the external walls. The simulation results for the centre of the particles indicate that the quickest drying occurs in Cy (at 38 s) and the slowest is Sp⁄ (at 112 s) as illustrated in Fig. 6 (b). For particles Sp, Cb, Cy, the temperature profile during drying seems to plateau more when compared to Sp⁄, Cb⁄, Cy⁄ (Fig. 6(a)). Simulation

1069

500

1

Bryden-x=1.6mm [30]

450

0.9

Bryden-x=1.2mm [30]

400

0.8

AR-x=1.6mm

350

0.7

300

0.6

250

0.5

200

0.4

150

0.3

100

+

*

AR-x=1.2mm

mH2O/mH2O,0 [-]

Tempreature [ °C]

H. Khodaei et al. / Applied Thermal Engineering 128 (2018) 1062–1073

HS-x=1.6mm HS-x=1.2mm HS-x=1.6mm HS-x=1.2mm

0.2

50

0.1

0 0

20

40

0

*

60

80

100

Time [s] Fig. 5. Two-centimetre half thickness of a cubic poplar wood with 30% moisture exposed to 1500 K background temperature. The dotted lines are moisture fractions. (*): Zero moisture point (dried wood condition) in the half thickness of cubic poplar at t = 80 s for HS, x = 1.2 mm corresponds to Ts = 140 °C.

(a)

(b) HS-Cy

460

HS-Cy*

0.8

HS-Sp

0.7

HS-Cb*

400

HS-Sp* Arrhenius

380 360

M H2O /MH2O,0

Temperature [K]

420

HS-Cy HS-Cb HS-Cy* HS-Sp HS-Cb* HS-Sp* Arrhenius

0.9

HS-Cb

440

1

0.6 0.5 0.4 0.3

340

0.2

320

0.1 0

300 0

20

40

60

80

100

0

120

20

40

60

80

100

120

Time [s]

Time [s]

Fig. 6. solid temperature at r = 0 Sp and Sp*, r = 0, h = 7.5 mm Cy, Cy* and r = 0, h = 7.5 mm Cb, Cb* in Heat Sink (HS) and Arrhenius (AR) models; (b): moisture fraction of the mentioned geometries at the same positions.

results reveal that in all geometries, the Arrhenius model is slower than the Heat Sink model. The temperature profiles in Arrhenius models seem not to be plateau in cylindrical and cuboid compare to Heat Sink model. In fact, the threshold temperature in Heat Sink model (Eq. (4)) limits drying process to start at Tevp = 373 K (evaporation temperature), however, in Arrhenius model (Eq. (3)), evaporation starts slowly from the beginning of numerical calculation (Compare AR and HS in Fig. 6(b)) and there is no threshold temperature to limit the evaporation rate. The numerical results emphasis that the behaviour of Arrhenius and Heat Sink models are almost similar in spherical geometries which could be due to the uniformity of boundary condition and temperature distribution in these numerical domains. The drying and temperature profiles of Cb are close to Cy⁄, which clearly demonstrates the importance of geometry in the evaporation process. Table 6 shows the total drying time in different geometries. It is obvious that after 38 s the Cy is com-

pletely dried which is 16 and 26 s before Cb, Sp. The drying times in Cy⁄, Cb⁄, Sp⁄ are approximately 2.9, 2.7, 3.3 times greater than Cy, Cb, Sp. These results indicate the importance of uniform distribution of temperature around the particle in reducing evaporation time. To provide a better visual comparison between Heat Sink and Arrhenius models, temperature and evaporation rate contours are illustrated in Fig. 7. The dotted lines separate two half geometries to show HS and AR contours in a single geometry. The XY (Z = 0) plane is identified to represent the contours in Cy⁄ due to the 3dimentional simulations in Cy⁄ (Fig. 7(g)). The results at 31 s for Cy and Cy⁄ (HS, AR) are presented because of their quicker evaporation time (Fig. 7 (b, e)) and at 50 s for the other particles. The highest drying rate is observed in Cy for the Heat Sink model ( 17 kg/m3s) Fig. 7(e). It means that in the Heat Sink model, cylindrical geometries dried faster than the other geometries, whereas,

Table 6 Total drying time in different biomass particle. A mixed 800 K temperature and convective heat flux are fully exposed on the external surface and partially exposed on the half surface area of the particles. SArea [mm2] at T = 800 K

TExternal [K]

Cy*, Cb*, Sp*

Cy, Cb, Sp

Cy*, Cb*, Sp*

238.7 259 199.1

477.52 518 398.3

293 293 293

Particle drying Time [s]

800 (*) 800 (*) 800 (*)

Cy, Cb, Sp

Cy*, Cb*, Sp*

Cy, Cb, Sp

800 800 800

110 149 215 (Max)

38 (Min) 55 64

(*): Only half of the particle exposed on mixed T = 800 K and natural convection and the other part exposed by natural convection and T = 293 K

1070

H. Khodaei et al. / Applied Thermal Engineering 128 (2018) 1062–1073

HS-Sp

(a)

(b)

HS-Cy

HS-Cb

(c)

X (m)

AR-Cy

AR-Sp

HS-Sp

(d)

AR-Cb (e)

(f)

HS-Cy*

(h) HS-Cb*

(i)

Y (m)

(g)

HS-Cb

AR-Cb

AR-Cy

AR-Sp HS-Sp*

HS-Cy

0.004

AR-Sp*

AR-Cy*

AR-Cb*

Fig. 7. Values of solid temperature Ts [K] (a–c) and drying rate Rdry [kg/m3 s] (d–f) for Heat Sink model (HS) and Arrhenius model (AR) in Sp, Sp*, Cb, Cb* at t = 53 s and Cy and Cy* at t = 31 s. Below the dotted line represents AR and above the dotted line represents HS results.

in the Arrhenius model, the lowest drying rate occurs in Sp⁄, Cb⁄ ( 4.5 kg/m3s) respectively as shown in Fig. 7(d, i). In all geometries, the reaction thickness in Heat Sink model is thinner than Arrhenius model, however, the propagation of drying zone and the rate of evaporation in the Heat Sink model are faster than the Arrhenius models. Threshold temperature in HS model may contribute to intensification of drying evaporation rate compare to AR model. A minor change in the temperature contours appears because of short interval temperature range that is required for drying process (Fig. 7(a–c)), however, these subtle temperature changes are more visible in Fig. 7(a) particularly for Cy, Cb. It is happened due to limiting evaporation by threshold temperature in HS model rather than AR model.

Considering the abovementioned analysis, cylindrical geometry is chosen due to the fastest drying rate, in addition that, the Heat Sink model is considered for moisture analysis according to the experimental and numerical evaluations addressed in Section 7.1. 7.3. Moisture evaporation in a cylindrical biomass particle The following parameters are analysed to investigate the effect of moisture content on thermo-physical properties of a single cylindrical biomass particle. Specific heat capacity: To analyse moisture contents, 8%, 20% and 30% of moisture are assumed in a single cylindrical wood pellet. The moisture fraction is calculated based on Eq.6 and drying

1071

H. Khodaei et al. / Applied Thermal Engineering 128 (2018) 1062–1073

algorithm proposed in Fig. 3. The average moisture fraction (Xmoist) of the particles is used to describe the total moisture remaining during the drying process taking into accounts drying times and the level of radiation temperature exposed on the external surface. The average specific heat capacity (Cp.p) and effective thermal conductivity (keff) have been recorded by different user defined memories during the evaporation process. Drying time in the particle increased from 38 s for 8% moisture content to 58 s for 20% and 70 s for 30%. The average specific heat capacity is highly sensitive to higher amount of moisture contents, however, when the moisture content starts to decrease, Cp.p reduces significantly and the influence of wood temperature is dominated. While the term Xmoist decreases, the influence of Xmoist in CP.P decreases and average Cp.p shifts to dry wood Cp.p. The results indicate that after evaporation, the Cp.p rises again as temperature increases in the dry wood. In fact, the term Xmoist in Eq. (15) causes a fluctuation in Cp.p graph due to the moisture content of the wood. Fig. 8 shows the evolution of moisture contents (8, 20, 30%) and average specific heat capacity of the cylindrical biomass particles during evaporation process. Effective thermal conductivity: In this study, the term krad proposed by Bryden et al. [31], Yang et.al. [18] is also considered in UDFs. Nevertheless, due to the lower temperature range required for the drying process, no significant change was observed in the solid temperature and moisture fraction contours; therefore, the term krad (Eq. (17)) is ignored in the next tests. The average effective thermal conductivity in particle with 8% moisture decreases very fast due to lowest amount of moisture and higher

heating rate, however, longer time requires to reduce keff in the higher moisture fractions. Furthermore, the slope of keff is increasing from minimum moisture to higher moisture levels which means that effective thermal conductivity is highly sensitive to the level of moisture in biomass particles. It should be noted that, when the particle dried completely, the thermal conductivity shifts to thermal conductivity of dry wood regarding the mass weighted average method used in this study (Eq. (16)). Fig. 9 shows the evolution of moisture contents (8, 20, 30%) and average effective thermal conductivity of the cylindrical biomass particles during evaporation process. Wood density: The evolution of particle density during drying process is illustrated in Fig. 10. Notably, mixed thermal condition is considered in the particle surfaces in Cy, Cy⁄ according to Table 5. Increasing moisture content from 8% to 30% results in increasing drying time from 36 to 73 s (about twice) in fully exposed temperature surface. Likewise, the drying time in cylindrical particle which is exposed on half external surface (Cy⁄) increases from 104 to 220 s (more than twice). This result clearly shows that higher amount of moisture contents in biomass and non-uniform distribution of temperature in external surface of the biomass particles lead to increasing the drying time from 36 s in Cy with 8% moisture to 220 s in Cy⁄ with 30% moisture content (approximately 6 times greater than Cy). In fact, higher residence time in Cy⁄ (30%) requires much higher operating costs and may result in lower efficient thermal conversion.

0.35

2800

0.3

2600

0.25

2400

0.2

2200

0.15

2000

0.1

1800

0.05

1600

Xmoist 8% Xmoist 20%

Xmoist 30% cp.p at xmoist 8%

0

cp.p at Xmoist 30%

cp.p [Jkg-1K-1]

X moist [-]

cp.p at Xmoist 20%

Xmoist=0

1400

0

20

40

60

80

Time [s] Fig. 8. The evolution of moisture contents and average specific heat capacity of cylindrical biomass particles fully exposed by Text = 800 K and natural convective.

0.35

0.35

Xmoist 8% Xmoist 20%

0.3

0.30

Xmoist 30% λeff at Xmoist 8% λeff at Xmoist 20%

0.25

0.2 0.15

0.20

0.1

λeff [Wm-1K-1]

X moist [-]

0.25

λeff at Xmoist 30% Xmoist=0

0.15 0.05

0.10

0

0

20

40

60

80

Time [s] Fig. 9. The evolution of moisture contents and average effective thermal conductivity of cylindrical biomass particles fully exposed by Text = 800 K and natural convective.

1072

H. Khodaei et al. / Applied Thermal Engineering 128 (2018) 1062–1073


Significant change in thermochemical properties of biomass (Effective thermal conductivity and specific heat capacity) is observed when the moisture content is increased.
Non-uniform distribution of temperature in biomass external surfaces and higher moisture contents result in higher drying times which impose higher operating costs, emissions and create a non-efficient biomass thermal conversion.

1300 Cy-Xmoist 8% Cy*-Xmoist 8% Cy-Xmoist 20% Cy*-Xmoist 20% Cy-Xmoist 30% Cy*-Xmoist 30% X moist =0

1250 1200

ρParticle [kg/m3]

1150 1100 1050 1000

Acknowledgements

950 900 850 800 0

30

60

90

120

150

180

210

240

270

Time [s] Fig. 10. Cylindrical biomass particles with 8, 20, 30% moisture contents which are fully exposed on Text = 800 K and natural convective (Cy) and partially exposed by Text = 800 K and natural convective flux in half external surface and Text = 298 K and natural convective heat flux in the other part (Cy*).

The result of this simulation will be used in pre-treatment of biomass particles in an industrial biochar reactor. For example, using conveyor system and exposing wood particle by indirect heat transfer from the syngas chamber which is located above the hearth section causes a longer residence time and low quality of biochar (mostly torrefied wood) due to non-uniform radiative and convective heat flux exposed from the syngas chamber, however, in the new design, the company is planning to fabricate a rotary kiln auger or a conveyor with agitation system to provide a uniform distribution of temperature during drying and pyrolysis process and minimize resident time as well as operating costs. The result of this study could be used in wood pellet industries and wood pellet (and wood) and wood dryers. 8. Conclusion A CFD comparative study of the drying in a standard single biomass particle is performed through the consideration of two evaporation models. Three thermally thick single biomass particles are simulated. The models predict additional features of the drying process in biomass particles which are less amenable to be measured directly by experimental tests. Solid phase variables in both Heat Sink and Arrhenius models are established by introducing moisture and energy transport equations. The interaction between the solid particles and the environment is identified by external radiation exposed on the surfaces of the particle coupled with natural cooling and immediate outflow of mass loss from the surface of particles. Comparisons are made of the profiles of temperature, moisture fraction, and average mass loss, drying rate, effective thermal conductivity and specific heat capacity inside the particles. The results are compared with different experimental measurements and numerical studies of single thermally thick particles considering both Arrhenius and Heat Sink models. Reasonable agreements among the temperature profile and moisture fraction profile are attained. The outcomes are:
The Heat Sink model depicts quicker drying time and higher evaporation rate compare to the Arrhenius model in all cases.
Cylindrical shape of biomass which is fully exposed by radiative and convective heats shows the fastest evaporation rate and the minimum drying time rather than the other geometries.
Reducing moisture content from 30% to 8% decreases drying time about 83%.

Dr. Yasir M. Al- Abdeli of Edith Cowan University is gratefully acknowledged for supervising the project. Jacobo Porteiro gratefully acknowledges financial support from the Government of Spain through the project ENE2015-67439-R. Edith Cowan University International Postgraduate Research Scholarship (ECU-IPRS) is highly acknowledged. References [1] M. Gómez, J. Porteiro, D. Patiño, J.L. Míguez, CFD modelling of thermal conversion and packed bed compaction in biomass combustion, Fuel 117 (2014) 716–732. [2] N. Wakao, S. Kaguei, Heat and Mass Transfer in Packed Beds, Gordon & Breach, London, 1982, p. 364. [3] G. Brenn, H. Steiner, I. Kondor, C. Planchette, A. Mercuri, G. Scharrer, A model for the non-isothermal drying of particle-laden pastes, Int. J. Heat Mass Transf. 92 (2016) 558–571. [4] S. Hermansson, H. Thunman, CFD modelling of bed shrinkage and channelling in fixed-bed combustion, Combust. Flame 158 (5) (2011) 988–999. [5] R. Rowell, Handbook of Wood Chemistry and Wood Composites, vol. 2, USA. 98, 2005. [6] Y. Yang, V.N. Sharifi, J. Swithenbank, Effect of air flow rate and fuel moisture on the burning behaviours of biomass and simulated municipal solid wastes in packed beds, Fuel 83 (2004) 1553–1562. [7] V. Minkova, M. Razvigorova, E. Bjornbom, Rolando Zanzi, T. Budinova, N. Petrov, Effect of water vapour and biomass nature on the yield and quality of the pyrolysis products from biomass, Fuel Proc. Technol. 70 (1) (2001) 53–61. [8] M. Tripathi, J.N. Sahu, P. Ganesan, Effect of process parameters on production of biochar from biomass waste through pyrolysis: a review, Renew. Sustain. Energy Rev. 55 (2016) 467–481. [9] W. Chan, M. Kelbon, B. Krieger, Modelling and experimental verification of physical and chemical processes during pyrolysis of a large biomass particle, Fuel 64 (11) (1985) 1505–1513. [10] S. Alves, J.L. Figueiredo, A model for pyrolysis of wet wood, Chem. Eng. Sci. 44 (12) (1989) 2861–2869. [11] H. Tran, R.H. White, Burning rate of solid wood measured in a heat release rate calorimeter, Fire Mater. 16 (4) (1992) 197–206. [12] M. Grønli, A theoretical and experimental study of the thermal degradation of biomass, in Mechanical engineering. 1996, Norwegian University of Science and Technology. p. 339. [13] H. Lu, W. Robert, G. Peirce, B. Ripa, L. Baxter, Comprehensive study of biomass particle combustion, Energy Fuels 22 (4) (2008) 2826–2839. [14] J. Saastamoinen, M. Aho, The simultaneous drying and pyrolysis of single wood particles and American Flame Research Committee, vol. 9, 1984, p. 11 [15] J. Porteiro, E. Granada, J. Collazo, D. Patiño, J. Morán, A model for the combustion of large particles of densified wood, Energy Fuels 21 (2007) 3151– 3159. [16] J. Collazo, J. Porteiro, D. Patiño, E. Granada, Numerical modeling of the combustion of densified wood under fixed-bed conditions, Fuel 93 (2012) 149–159. [17] J. Porteiro, J.L. Míguez, E. Granada, J.C. Moran, Mathematical modelling of the combustion of a single wood particle, Fuel Proc. Technol. 87 (2) (2006) 169– 175. [18] Y. Yang, V.N. Sharifi, J. Swithenbank, L. Ma, L.I. Darvell, M. Jones, M. Pourkashanian, A. Williams, Combustion of a single particle of biomass, Energy Fuels 22 (2008) 306–316. [19] R. Mehrabian, A. Shiehnejadhesar, R. Scharler, I. Obernberger, Multi-physics modelling of packed-bed biomass combustion, Fuel 122 (2014) 164–178. [20] B. Peters, E. Schröder, C. Bruch, T. Nussbaumer, Measurements and particle resolved modelling of heat-up and drying of a packed bed, Biomass Bioenergy 23 (4) (2002) 291–306. [21] B. Peters, Thermal conversion of solid fuels, Computational Mechanics, vol. 15, 2003. [22] A.L. Ljung, T.S. Lundström, B.D. Marjavaara, Tano, K, Convective drying of an individual iron ore pellet–Analysis with CFD, Int. J. Heat Mass Transf. 54 (17) (2011) 3882–3890. [23] B. Peters, Classification of combustion regimes in a packed bed of particles based on the relevant time and length scales, Combust. Flame 116 (1) (1999) 297–301.

H. Khodaei et al. / Applied Thermal Engineering 128 (2018) 1062–1073 [24] H. Thunman, Davidsson, Kent, Leckner, Bo, Separation of drying and devolatilization during conversion of solid fuels, Combust. flame 137 (1) (2004) 242–250. [25] R. Johansson, H. Thunman, B. Leckner, Influence of intraparticle gradients in modeling of fixed bed combustion, Combust. Flame 149 (1) (2007) 49–62. [26] R. Mehrabian, S. Zahirovic, R. Scharler, I. Obernberger, S. Kleditzsch, S. Wirtz, V. Scherer, H. Lu, L.L. Baxter, A CFD model for thermal conversion of thermally thick biomass particles, Fuel Process. Technol. 95 (2012) 96–108. [27] H. Khodaei, Y.A. Al-Abdeli, F. Guzzomi, G.H. Yeoh, An overview of processes and considerations in the modelling of fixed-bed biomass combustion, Energy 88 (2015) 946–972. [28] J. Wurzenberger, Thermal conversion of biomass comprehensive reactor and particle modelling, AIChE J. 48 (2002) 2398–2411. [29] M. Bellais, Modelling of the pyrolysis of large wood particles, KTH, 2007, p. 95. [30] X. Zhang, Q. Chen, R. Bradford, V.N. Sharifi, J. Swithenbank, Experimental investigation and mathematical modelling of wood combustion in a moving grate boiler, Fuel Process. Technol. 91 (2010) 1491–1499. [31] K. Bryden, K.W. Ragland, C.J. Rutland, Modeling thermally thick pyrolysis of wood, Biomass Bioenergy 22 (1) (2002) 41–53. [32] M. Bellais K. Sjöström J. Svenson J.B.C. Pettersson A. Omrane F. Ossler M. Aldén Fast drying of large wood particles under pyrolysing conditions: experimental study and modelling Fuel 2007 0016–2361. [33] M. Gómez, Porteiro, J, Patiño, D, Míguez, JL, Fast-solving thermally thick model of biomass particles embedded in a CFD code for the simulation of fixed-bed burners, Energy Convers. Manage. 105 (2015) 30–44. [34] K. Bryden, Computational Modeling of Wood Combustion, in Agriculture, Mechanical engineering, University of Wisconsin-Madison, USA, 1998.

1073

[35] R. Yuen, G. de Vahl Davis, E. Leonardi, G.H. Yeoh, The influence of moisture on the combustion of wood, Numer. Heat Transfer Part A Appl. (3) (2000) 38257– 38280. [36] A. Atreya, Pyrolysis, Ignition and Fire Spread on Horizontal Surfaces of Wood, Harvard University, Cambridge, 1983. [37] B. Peters, Measurements and application of a discrete particle model (DPM) to simulate combustion of a packed bed of individual fuel particles, Combust. Flame 131 (1) (2002) 132–146. [38] B. Krieger-Brockett, David S. Glaister, in: Wood Devolat1lization-Sensitivity to Feed Properties and Process Variables, in Research in Thermochemical Biomass Conversion, Springer, 1988, pp. 127–142. [39] M. Miltner, A. Makaruk, M. Harasek, A. Fridel, CFD-modelling for the combustion of solid baled biomass, Parameters 1 (2006) 2. [40] J. Porteiro, E. Granada, J. Collazo, D. Patiño, J. Morán, Experimental analysis of ingnition front propagation of several biomass fuels in a fixed-bed combuster, Fuel 89 (2010) 26–35. [41] R. Johansson, H. Thunman, B. Leckner, Sensitivity analysis of a fixed bed combustion model, Energy Fuels 21 (2007) 1493–1503. [42] K. Ragland, Aerts, DJ, Baker, AJ, Properties of wood for combustion analysis, Bioresource Technol. 37 (2) (1991) 161–168. [43] F.P. Incropera, Adrienne S. Lavine, David P. DeWitt, Fundamentals of Heat and Mass Transfer, John Wiley & Sons Incorporated, 2011. [44] C. Heidenreich, D.K. Zhang, Measuring the temperature response of large wet coal particles during heating, Fuel 78 (8) (1999) 991–994. [45] K. Bryden, K.W. Ragland, Numerical modeling of a deep, fixed bed combustor, Energy Fuels (1996) 269–275.