Chemical Engineering Research and Design 191 (2023) 414–425

Available online at www.sciencedirect.com

Chemical Engineering Research and Design

journal homepage: www.elsevier.com/locate/cherd

A distributed activation energy model for cellulose

pyrolysis in a fluidized bed reactor
]]
]]]]]]
]]

⁎ ⁎⁎
Samreen Hameed a,c, Abhishek Sharma b,d, , Vishnu Pareek c,
a
Chemical, Polymer & Composite Materials Engineering Department, University of Engineering & Technology, New
Campus, Lahore, Pakistan
b
Department of Chemical Engineering, Manipal University Jaipur, Rajasthan, India
c
Western Australian School of Mines: Minerals, Energy and Chemical Engineering, Curtin University, WA,
Australia
d
Chemical & Environmental Engineering, School of Engineering, RMIT University, Melbourne, Victoria, Australia

a r t i c l e i n f o a b s t r a c t

Article history: Cellulose pyrolysis is used to give solid char, condensable vapors and non-condensable
Received 8 September 2022 gases. This is a complex process and to model this using CFD simulations, Euler-Euler
Received in revised form 10 January approach with multi-fluid is applied. The objective of the present study is to develop a CFD
2023 model to combine the reaction kinetics with the hydrodynamics to study the cellulose
Accepted 26 January 2023 pyrolysis in the fluidized bed reactor. The distributed activation energy model to compare
Available online 29 January 2023 the product yield with and without distribution of activation energy has been in­
corporated. Simulations of cellulose pyrolysis in a fluidized bed reactor using two different
Keywords: kinetic schemes have been conducted. Model has been validated using global kinetic
Modeling scheme and then optimized kinetic parameters are introduced with a distribution of ac­
Kinetic tivation energy. The predicted values of activation energies, frequency factor and stan­
Fluidized bed reactor dard deviation impact the overall pyrolysis product yield. The yield of tar and char
Pyrolysis increases, however the fraction of gases decreases significantly because of higher acti­
Cellulose vation energy. Tar yield is 82.1 wt% in the presence of DAEM, while the gas yield is 8 wt%
DAEM and char yield comes out to be 9.9 wt%.
© 2023 Published by Elsevier Ltd on behalf of Institution of Chemical Engineers.

developed a multi-fluid model for different biomass and cellu­

1. Introduction
lose pyrolysis using CFD simulations and predicted the key
features of the fast pyrolysis in a fluidized bed reactor (Xue et al.,
To model the whole process of pyrolysis in fluidized bed sys­
2012, 2011). Recently, a model was developed by Jalalifar et al.,
tems, a multiphase fluid model must be integrated with reaction
(2017) using Eulerian-Eulerian approach for pure cellulose pyr­
kinetics using computational fluid dynamics (CFD). Several re­
olysis in a bubbling fluidized bed reactor. The model (Jalalifar
search articles (Bharadwaj et al., 2004; Czajka et al., 2016; Jalalifar
et al., 2018) was later used for biomass pyrolysis in a bubbling
et al., 2018, 2017; Sharma et al., 2015; Soria-Verdugo et al., 2015;
fluidized bed reactor to investigate the impact of different op­
Xue et al., 2012, 2011) are available in the literature discussing
erating parameters on the product yield.
the coupling of biomass pyrolysis models with CFD.Xue et al.

Abbreviations: BFB, bubbling fluidized bed; CFD, computational fluid dynamics; DAEM, distributed activation energy model; EE,
Eulerian-Eulerian; FBR, fluidized bed reactor; TGA, thermo-gravimetric analysis
⁎
Corresponding author at: Department of Chemical Engineering, Manipal University Jaipur, Rajasthan, India.
⁎⁎
Corresponding author.
E-mail addresses: abhishek.sharma@jaipur.manipal.edu (A. Sharma), v.pareek@curtin.edu.au (V. Pareek).
https://doi.org/10.1016/j.cherd.2023.01.048
0263-8762/© 2023 Published by Elsevier Ltd on behalf of Institution of Chemical Engineers.
 Chemical Engineering Research and Design 191 (2023) 414–425 415

biomass pyrolysis and the fluctuations in the tar yield with

Nomenclature
hydrodynamics of bubbling bed. The model was further ex­
A Pre-exponential factor (s−1). tended and coupled with DAEM to study the influence of
Cp Specific heat (J/kg-K). distribution of activation energy (Xiong et al., 2016b). The
D Diffusivity (m2/s). researchers focused on the comparison of two different dis­
E Activation Energy (J/kgmol). tribution schemes (Gaussian and Logistic) and their effect on
h Convective heat transfer coefficient. the exit tar yield. Also, the effect of variation in activation
M Molecular weight (kg/kmol). energy was accounted for by using coefficient of variation,
P Particle pressure (atm). skewness and kurtosis. DAEM has widely been applied to
R General gas constant. better understand the complexity of biomass pyrolysis and
T Temperature (K). predict the reaction kinetics for different lignocellulosic
U Gas inlet velocity (m/s). biomass materials (Cai et al., 2011; Sonobe and
Umf Minimum fluidization velocity (m/s). Worasuwannarak, 2008; White et al., 2011) Morphological
υ Vapor velocity (m/s). changes in the cellulose structure during the pyrolysis pro­
V Particle volume (m3). cess have revealed that the process involves several reac­
tions which have different activation energies. Hence, using
Greek Letters DAEM for cellulose pyrolysis indeed gives accurate product
α Volume fraction. yield when compared with single activation energy value.
k Thermal conductivity (W/m-C). The objective of the present work is to develop an in­
ρ Density (kg/m3). tegrated model of the cellulose pyrolysis in a fluidized bed
μ Viscosity (kg/m-s). reactor encompassing both CFD and DAEM. Broido-
θ Granular Temperature (m/s)2. Shafizadeh scheme for cellulose decomposition has been
ϕ Specularity coefficient. selected and the kinetics obtained from TGA (Hameed, 2019;
Hameed et al., 2022). Distributed activation energy model is
Subscript used to specify the activation energies and the effect of this
Cell Cellulose. distributed activation energy on the pyrolysis product yield
A.Cell Active cellulose. (yield of tar, char and gases) has been determined. This study
C Char. is divided into three steps: first step includes the kinetic
G Gases. model validation in CFD using Briodo-Shafiadeh scheme and
T Tar. the kinetic parameters available in the literature (Di Blasi,
g Gas phase. 1994). The second step is to include the reaction kinetics
v Vapor phase. predicted and optimized by TGA experiments and DAEM in
s Solid phase. multiphase model. The third step then focuses on the cou­
pling of DAEM to multiphase model using Gaussian dis­
tribution. Pyrolysis product yields determined from the CFD
Since cellulose is the major component of lignocellulose
simulations will be compared with the experimental data
biomass, there are several kinetic schemes reported in the
available in the literature and then simulations will be ex­
literature which are used to understand the cellulose pyr­
tended for the case when DAEM has been integrated
olysis but not all of them have been used for the numerical
with CFD.
simulation of the pyrolysis. Broido-Shafizadeh scheme is the
most commonly used for cellulose pyrolysis modeling. This
2. Mathematical modeling description
scheme only includes a single value of activation energy for
each reaction. Since cellulose pyrolysis involves multiple
Multiphase reactive CFD models are formulated by coupling
reactions with different kinetic energies, it is essential to
the kinetics of devolatilization of cellulose particles with the
include the range of activation energies. In addition, as re­
hydrodynamics of fluidizing bed in the reactor. The selection
ported by Burnham et al., (2015) first-order reaction con­
of a kinetic model and the multi-fluid model used for cellu­
stants derived from single heating rate do not always give
lose pyrolysis have been discussed below in detail.
true values of frequency factor, A and activation energy, E
and specially in case of cellulose, the product yield is affected
2.1. Kinetic model
because of the high values of A.
The distributed activation energy model (DAEM) is a
The kinetic scheme used in this work is shown in Fig. 1. This
model having multiple reactions. This is used to understand
is the Broido-Shafizadeh scheme and was first proposed by
biomass decomposition by considering large number of in­
Shafizadeh and Bradbury (Bradbury et al., 1979; Shafizadeh,
dependent, parallel, first order or nth-order reactions with
1982; Shafizadeh and Bradbury, 1979) as a set of three reac­
distributed activation energies.(Cai et al., 2014) DAEM has
tions and was later modified to include the tar decomposi­
gained researcher’s attention in recent years and has been
tion. According to this scheme, cellulose (Cell) is converted to
implemented to better understand the thermal decomposi­
some intermediate named as “active cellulose” (A.Cell) in the
tion of biomass (Cai et al., 2014; Czajka et al., 2016; Kirtania
first reaction. The intermediate (A.Cell) then immediately
and Bhattacharya, 2015; Soria-Verdugo et al., 2015). In the
decomposes to other products tar (T), char (C) and gases (G)
previous studies DAEM kinetics implemented to biomass
following two parallel reaction schemes. Tar (T) produced as
considered only single particle kinetics (Cai et al., 2014) ig­
a result of second reaction undergoes decomposition to
noring the effect of turbulent mixing and transport phe­
produce gases (G). Cellulose, active cellulose and char are
nomena. Xiong et al., (2016a) simulated a laboratory-scale
included in the solid phase (s1) while tar vapors and the
bubbling fluidized bed reactor to predict the product yield of
gases are considered as the gas phase (g) along with the inert
416 Chemical Engineering Research and Design 191 (2023) 414–425

position reaction as specified by Xiong et al. (Xiong et al.,

2016b) can be described as:
[dm(E)] E
= Aexp [dm(E)]
t RT (7)
Fig. 1 – Cellulose reaction scheme (Bradbury et al., 1979).
However,
dm(E) = f(E)mdE (8)
gas nitrogen (N2). Sand particles (s2) make the secondary
Upon integration the equation will become:
solid phase, this phase does not react and act as an inert
fluidizing media in the reactor. m E
= Am exp f(E)dE
The rate of generation and consumption of the solid t 0 RT (9)
phase (s1) and gas phase (g) species can be defined as:
This distribution function f (E), is an important parameter
Cellulose, R s1,Cell = k1 s1 s1 X s1,Cell (1) in the DAEM which may be specified by Weibull, Logistic or
Gaussian distribution function. In this work Gaussian dis­
Active Cellulose, R s1,A.Cell tribution function has been used which can be represented
= k1 s1 s1 X s1,Cell (k2 + k3) s1 s1 Xs1,A.Cell (2) by the following equation:

1 (E E0 )2
Tar Vapors, Rg,T = k2 s1 s1 Xs1,A.Cell k4 g g Xg (3) f(E) = exp 2
2 2 (10)

Gas, Rg,G Where E0 is the mean activation energy and σ is the standard
deviation. Kinetic parameters, frequency factor, A and acti­
= k3 s1 s1 Xs1,A.Cell (1 Y) + k 4 g g
vation energy, E were predicted using model fitting approach
g g
Xg + [ s1 s1 Xs1,A.Cell (k2 + k3)] +[ s1 s1 Xs1,A.Cell k3Y] in which the objective function was minimized by using
b c
fmincon in MATLAB as discussed in detail in the work of
(4) Hameed et al. (Hameed et al., 2022). These optimized kinetic
Char, R s1,C = parameters are given below in Table 3:
s1 s1 Xs1,A.Cell k3Y (5)

“Y” represents the ratio of char formed in the third reac­ 2.3. Multi-fluid model
tion, which is 0.35 (Xue et al., 2012).
All the reactions in this scheme are considered first-order The multiphase model used in this study describes the gas
and the reaction kinetics are determined according to and solid phases as inter-penetrating continua in Eulerian-
Arrhenius law, as: Eulerian framework. During the pyrolysis in a fluidized bed
k = Ae E/RT (6) reactor there exist three different phases: gas phase as pri­
mary phase and the two solid phases as the secondary
The data of kinetic parameters and the heat of reaction is
phases. Conservation of mass, momentum, energy and the
presented in Table 1:
other basic laws is specified using the models described
Thermo-physical properties of species involved in the
below (Fluent, 2009). In addition to basic conservation laws,
pyrolysis of cellulose are presented in Table 2.
some other correlations necessary for the implementation of
multi-fluid model are also discussed in the following sec­
2.2. Distributed activation energy model
tions:

Distributed activation energy model (DAEM) is comprised of a

2.3.1. Gas phase
large number of independent parallel reactions with dif­
For the gas-mixture phase the conservation of mass, mo­
ferent activation energies, represented by a continuous dis­
mentum and energy equations are formulated. The con­
tribution function (Hu et al., 2016). There are certain
tinuity equation can be defined as:
assumptions which are important for the derivation of DAEM
as reported in the literature (Zhu et al., 2021): (1) thermal ( g g) + .( g g vg) = Rg
decomposition of biomass includes several independent, t (11)
parallel and first order reactions (2) each reaction has its Rg represents the source term for the consumption/gen­
unique activation energy but share the same frequency eration of the gas phase species in the reactor. This source
factor (3) a continuous distribution can be used to describe all term defines the interphase mass transfer taking place at the
reactions. In DAEM the kinetics are not described by the solid-fluid interphase due to the chemical reactions and is
conventional Arrhenius law but are specified using the the sum of the reaction rates of all the species in the gas
probability density function. The rate of biomass decom­ phase, which can be defined as:

Table 1 – Kinetic Parameters and Heat of Reaction Data.

Reaction Pre-exponential Activation Reference Heat of reaction Reference
constant factor, A (sec−1) energy, E (J/kmol) (kJ/kg)

k1 2.8 × 1019 2.42 × 108 (Di Blasi, 1994) 0 (Koufopanos et al., 1991)
k2 3.28 × 1014 1.965 × 108 (Di Blasi, 1994) 255 (Koufopanos et al., 1991)
k3 1.3 × 1010 1.505 × 108 (Di Blasi, 1994) -20 (Koufopanos et al., 1991)
k4 4.25 × 106 1.080 × 108 (Di Blasi, 1994) -42 (Curtis and Miller, 1988)
 Chemical Engineering Research and Design 191 (2023) 414–425 417

Table 2 – Thermo-physical properties of species (Jalalifar et al., 2017; Lathouwers and Bellan, 2001).
Species Density, ρ Particle Molecular Weight, Heat Capacity, Dynamic Thermal Conductivity, k
(kg/m3) Diameter, ds (m) M.W (kg/kmol) Cp (J/kg.K) Viscosity, μ (J/kg.K)
(kg/ms)

Cellulose 400 5.0 × 10−5 342.297 2300 - 0.18

Active 400 5.0 × 10−5 342.297 2300 - 0.18
Cellulose
Tar - - 100 950 1.3 × 10−5 1.32
Char 2330 5.0 × 10−5 12.01 1100 - 0.0878
Gases - - 23 4545.09 1.32 × 10−5 0.0582
Nitrogen - - 28.04 1040.7 3.49 × 10−5 0.0242
Sand 2649 5.2 × 10−4 60.08 860 0.2

chemical reactions is included using source term Hrg and

Table 3 – Optimized kinetic parameters for pure
the inter-phase heat transfer between the gas and the solid
cellulose including standard deviation (Hameed, 2019).
phase is given by Hsg . These terms can be defined as:
Reaction Parameters Optimized Values
Hg = Cp,g dTg (16)
Reaction 1 Activation energy, E1 (kJ/mol) 232.75
Frequency factor, A1 (s−1) 8.52 × 1017
Standard deviation, σ1 (kJ/mol 2.418 Hsg = hsg (Tsn Tg ) (17)
Reaction 2 Activation energy, E2 (kJ/mol) 196.5
Frequency factor, A2 (s−1) .28 × 1014
Standard deviation, σ2 (kJ/mol) 0 2.3.2. Solid phase
Reaction 3 Activation energy, E3 (kJ/mol) 150.5 For each sth solid phase, in a system of n solid phases the
Frequency factor, A3 (s−1) 1.3 × 1010 equations of continuity, momentum and energy are for­
Standard deviation, σ3 (kJ/mol) 0
mulated and discussed in this section. Continuity equation
Reaction 4 Activation energy, E4 (kJ/mol) 255.97
for solid phase has been defined as:
Frequency factor, A4 (s−1) 5.51 × 1015
Standard deviation, σ4 (kJ/mol) 33.99
( sn sn) + .( sn sn vsn) = R sn
t (18)

sn , sn and vsn are the bulk density, volume fraction and

Rg = Rg,T + Rg,G (12) velocity of the nth solid phase respectively. The source term
for solid phase Rsn is the sum of the reaction rates of all the
For ith specie in the gas phase, the conservation equation
species in the solid phase including cellulose, active cellulose
is given as:
and char.
( g g Xg,i) + .( g g vg Xg,i) = .( g g Dg,i Xg,i) + Rg,i Rsn = R s1,Cell + R s1,A.cell + Rs1,C (19)
t (13)
For ith specie in the gas phase, the conservation equation
Where Rg,i represents the reaction rate term for the ith specie
is given as:
and Xg,i and Dg,i are the mass fraction and diffusion coeffi­
cient respectively.
( sn sn Xsn,i) + .( sn sn vsnX sn,i)
The momentum balance for the gas phase can be defined t
using Navier-Stokes equation as: = .( sn sn Dsn,i Xsn,i) + R sn,i (20)
( g g vg ) Momentum equation for sth solid phase is:
+ ( g g vg vg)
t
n
( sn sn vsn)
= P+ . + g + ( Rgs + vsg msg vgsmgs) + ( sn sn vsn vsn)
g g g g s=1 t
(14) = sn p psn + . sn + sn sn
N
In the first term on the right side of the equation, P re­ g + l =1 (Rln + vsg msg vgsmgs) (21)
presents the pressure shared by all the phases. The second
In the momentum transport equation, psn represents of
term, g is the stress-strain tensor and the acceleration due
nth solid phase. sn is the stress-strain tensor for nth solid
to gravity is represented by g . Rgs is the interaction force
phase and the acceleration due to gravity is represented by
between the gas phase and the solid phases due to mo­
mentum transfer happening because of the drag force be­ g . Rln is the interaction force between lth gas phase or sth
tween the different phases and the last two terms are solid phase and nth solid phases including momentum
representing the mass transfer between the phases due to transfer due to the drag force between the gas and the solid
thermal degradation of particles. phase. The last two terms are representing the mass transfer
For gas phase the energy conservation equation is between the phases same as explained for gas phase in
given as: Eq. 11.
The energy conservation equation for nth solid phase
( g g Hg ) n
will be:
+ .( g g vg Hg ) = . qg + Hrg + Hsg
t s=1 (15) n
( sn sn Hsn )
+ .( sn sn vsnHsn) = . qsn + Hrsn + Hgs
Where, Hg is the specific heat of the gas phase and qg is the t s=1

conductive heat flux. Enthalpy change occurring due to (22)

418 Chemical Engineering Research and Design 191 (2023) 414–425

Where, Hsn and qsn are the specific heat and conductive heat
flux of the nth solid phase. Hrsn is the source term which
incorporates the enthalpy change due to chemical reactions
in the solid phase. However, the inter-phase heat transfer
between the gas and solid phases is given by Hgs. Since there
is no radiative heat transfer in the system hence the terms
including radiative heat transfer have been neglected. These
terms can be defined as:

Hsn = Cp,sn dTg (23)

Hgs = hgs (Tg Tsn) (24)

According to kinetic theory of granular flow (KTGF) kinetic

energy of fluctuating particles needs to be defined for a
multi-phase system. This kinetic energy is usually specified
using granular temperature of the nth solid phase, θsn.
Granular temperature is proportional to the mean square of
the fluctuating velocity of solid particles. The algebraic form
of this fluctuating energy is given as:

3 ( sn sn sn )
= (psn I + sn): vsn + sn + ln
2 t (25)

This equation is a combination of different forms of en­

ergy. These terms may be defined as:
(ps I + s): vs = Energy produced due to solid stress
tensor.
s = Energy dissipation due to collisions for nth solid
phase.
ls = Energy exchange between the lth gas or sth solid
phase and nth solid phase.
Other than basic conservation laws, the constitutive re­
lations for gas-solid interaction in a multiphase system are
required which have been described in detail in Appendix A.
This includes the stress-strain tensor relations for gas and
Fig. 2 – Fluidized bed reactor for cellulose pyrolysis.
solid phase, solids pressure and drag models. Fourier’s law
has been used to calculate the conductive heat flux and for
convective heat transfer coefficient Nusselt number correla­
Table 4 – Reactor dimensions and computational
tion by Gunn (Gunn, 1978) has been used.
domain data.
Parameter Value
3. CFD simulation strategy
Reactor Height, H 0.3429 m
The governing equations of momentum, pressure and vo­ Reactor Diameter, D 0.0381 m
lume fractions were solved using phase coupled SIMPLE al­ Bed Height, h 0.055 m
Biomass inlet Diameter, d 0.0073 m
gorithm in ANSYS FLUENT version 17.2. The coupled
Grid Cells in horizontal direction 10
algorithm solves the momentum and pressure-based con­
Grid Cells in vertical direction 91
tinuity equations together and pressure correction equation
coefficients are calculated based on the momentum equa­
tions. According to this scheme the coupled per phase velo­
cities are solved in a segregated manner (Fluent, 2009). izontal and vertical direction. The 2-D reactor is shown in Fig. 2
However, the algorithm is coupled with iterative time ad­ with the reactor configuration data presented in Table 4.
vancement scheme in which all the equations are solved Sand being used as an inert fluidizing material in the bed
through iterations until the convergence occurs. was patched to the height of 0.055 m with the volume frac­
tion of 0.59 and with an initial temperature of 773 K. The inlet
3.1. Solution procedure with initial and boundary velocity of the bed material was zero. Cellulose was entered
conditions to the reactor at a flow rate of 2.77 × 10−5 kg/sec through the
side inlet in a continuous manner at a temperature of 300 K.
The configuration used for this part of work is same as the one The gas phase entered through the gas inlet at the bottom of
used in experiments and simulations of Xue et al., (2012, 2011). the reactor with the superficial gas velocity of Ug = 3.5 Umf
The computational domain for this system is a 2-dimensional and a temperature of 773 K. At the outlet of the reactor
(2-D) reactor which has been divided into fixed number of pressure boundary condition was specified. The walls were
control volumes through specific number of grid cells in hor­ kept isothermal with a constant temperature of 773 K. The
 Chemical Engineering Research and Design 191 (2023) 414–425 419

Table 5 – Pyrolysis product yield (wt%) comparison between experiment and simulations for pure cellulose.
Data Bio-oil (wt%) Char (wt%) Gases (wt%)

Experimental data (Xue et al., 2012) 82.1 2.2 12.4

Current Simulations 81.1 3.1 15.8

value of co-efficient of restitution was taken as 0.97 and the than the total volume as stated by Xue et al., (2012). The si­
angle of internal friction was 55. For boundary conditions of mulated results showed an increase in amount of gases
walls, no slip was specified for the gas phase while for the that is the actual amount of gases leaving the fluidized bed
solid phase partial slip condition was used. reactor.
Energy conservation and species transport models were
used to account for the pyrolysis reactions. The reaction
terms have been integrated using a stiff ODE (Ordinary
Differential Equations) solver. For time discretization, second 4.2. Implementation of DAEM to CFD
order implicit method was used with the time step of 10−4
sec. QUICK algorithm has been used to calculate the volume Kinetic parameters used in the current model were predicted
fraction of each phase. Conservation equations for mo­ using thermogravimetric analysis and DAEM (Hameed, 2019).
mentum, energy and the species transport equations have The predicted kinetic parameters with the mean activation
been solved using second order (upwind) method. All simu­ energy and the standard deviation for pure cellulose are
lations were run for the simulation time of 150 s when summarized previously in Table 3. Kinetic parameters for
pseudo-steady state was achieved in the system. reaction 1 and 4 were specified with standard deviations of
activation energy; however, for reactions 2 and 3 only single
values of activation energy were used. These kinetic para­
4. Results and discussion
meters have been used to understand the cellulose pyrolysis
behavior following the same reaction scheme as described in
4.1. Model validation
Figure1, in a fluidized bed reactor using CFD simulations.
The analytical solution of the DAEM is difficult therefore
In this study, the CFD model was validated using the ex­
researchers have focused on the numerical solution of this
perimental data of Xue et al. (2012), and the same model has
model. Braun and Burnham used E0 3 andE0 + 3 as the
been used to calculate the pyrolysis product yield with dis­
integral lower and upper limits respectively and solved the
tributed activation energy model which will be discussed in
model using 96 intervals for activation energy (Braun,
the next section in detail. In this part of the research, all
Burnham, 1987). Xiong et al., have previously carried out the
operating parameters and material properties have been
integration of DAEM with CFD, using different intervals such
kept same as of experimental data. The product yield (η) was
as: (E0 3 , E0 + 3 ), (E0 4 , E0 + 4 ), (E0 5 , E0 + 5 ) and
calculated using the equation:
(E0 6 , E0 + 6 ). They concluded that if the integration in­
t tervals and their spaces are chosen optimally, then this
0 outlet
( X)dAdt + reactor
( X)dV
(t) = coupling does not increase the computational power prohi­
Mcell_feed (26)
bitively (Xiong et al., n.d.). Keeping in view the outcome of
The comparison of the product yield obtained from si­ previous reported literature, the integration space chosen for
mulations using the literature kinetic data with the experi­ this case was (E0 3 , E0 + 3 ) with the 60 equal intervals.
mental data is presented in Table 5. Mean activation energies have been calculated and pre-ex­
The comparison of the present simulation results with the ponential factor has also been modified using the Gaussian
experimental data shows that the predicted bio-oil yield was distribution function to incorporate the probability function.
comparable to the experimental data reported in literature. This has led to 124 chemical reaction equations. Reaction
The char and non-condensable gases yield, however, was terms along with the transport equations have been solved in
slightly over-predicted by the simulation when compared ANSYS FLUENT to predict the pyrolysis product yield. The
with the experimental results. Nevertheless, it is reasonable yield of bio-oil, char and gases is presented in Table 6.
to conclude that the product yields predicted by the model A comparison of product yield from experimental data
were in close agreement with the experimental data, thus and simulation results using Broido-Shafizadeh scheme,
validating the CFD model. The mass closure from the ex­ with and without DAEM is shown in Table 6. It is clear that
perimental data is around 96.7% which is expected due to the with DAEM the predicted bio-oil yield increased significantly.
error in calculating the mass percentage of gases. The cal­ This observation is in agreement with those reported in the
culated volume of non-condensable gases came out to be less literature (Xiong et al., 2016b) However, the yield of other

Table 6 – Pyrolysis product yield for pure cellulose with and without using DAEM.
Data Bio-oil (wt%) Char (wt%) Gases (wt%)

Experimental data (Xue et al., 2012) 82.1 2.2 12.4

Without DAEM 77.8 3.4 18.9
With DAEM 82.1 9.9 8.0
420 Chemical Engineering Research and Design 191 (2023) 414–425

velop near the biomass inlet due to decomposition of cellu­

lose to volatiles. These bubbles combine with those produced
in the bed due to fluidizing gas (N2) and form big bubbles
which move upward and then leave the surface of the bed
(Fig. 3(a)). This process continues and promotes the mixing
behavior of the particles. Similar behavior can be observed
while looking into the profiles of the sand phase (Fig. 3(b)).
Time-averaged volume fractions of the three phases in
the fluidized bed reactor have also been plotted in Fig. 4. This
figure shows that there exists some segregation in the bed at
this superficial gas velocity, which may affect the rate of heat
and mass transfer and ultimately the rate of decomposition
of cellulose. Hence, to let the pyrolysis happen at faster rates,
it is recommended to use an optimized value of gas velocity.
Fig. 5 shows instantaneous temperature profiles of gas
phase (gas mixture of N2, tar and gases), sand phase and the
cellulose phase (mixture of cellulose, active cellulose
and char).
Since the operating temperature for reactor had been set
to 773 K, there was a small variation in the temperatures of
gas and sand phase. Sand was packed at the temperature of
773 K, and the gas phase also entered at the same tempera­
ture. The small change in temperature of both phases
through the bed (772–773 K) could be due to production of
volatiles n the fluidized bed as a result of cellulose decom­
position. However, for cellulose there was a significant
change in temperature through the fluidized because it was
heated from its initial temperature of 300–300–773 K.
However, the heating was rapid due to both conductive and
convective heat transfers from the other two phases re­
sulting in cellulose to stay at a temperature of 728 K through
the most part of the reactor.
Time-averaged temperature profiles for the three phases
gas phase, cellulose phase and the sand phase as a function
of the reactor height are plotted in Fig. 6. It is clear that the
gas phase temperature was uniformly distributed in the re­
cator near its intial value of 773 K. The temperature of sand
also showed little variation through the bed. For the cellulose
Fig. 3 – Instantaneous volume fraction of (a) gas phase (b) phase, the temperature profile did show a different behavior.
sand phase in a bubbling fluidized bed reactor with time The temperature of cellulose decreased at the inlet of the
interval of 0.1 s between each profile. feed where fresh cellulose entered in the reactor at a tem­
perature of 300 K. However, thereafter the temperature of
pyrolysis product species was also affected due to standard cellulose increased rapidly due to the heat transfer between
deviation of activation energies. For example, the char yield cellulose and heated gas.
was over-predicted and that for gases was under -predicted. The time-averaged vertical velocity vectors for all three
These results show that the model including distributed phases are shown in Fig. 7, which shows that for gas phase
activation energy predicted product yield more analogous to the maximum value of gas velocity was more than its initial
experimental data than those using a single value of activa­ value due to production of volatiles (tar and gases). The gas-
tion energy. This is because cellulose pyrolysis is a combina­ mixture flows in the upward direction in the center of the
tion of multiple reactions, using an averaged single value of column while solids particles recirculate in the downward
activation energy may result in erroneous results if the dif­ direction along the reactor walls. The cellulose particles are
ference between activation energies of reaction is significant. mixed with the sand and gas phase near the feed inlet and
A small change in the value of activation energy may affect then start to decompose. Solid particles flow down with the
the rate of reaction exponentially because in CFD simulations minimum averaged velocity. Sand particles keep circulating
rate of reaction is determined using Arrhenius law. Using a in the upward direction with the maximum velocity and flow
distribution of activation energies includes more number of down with the minimum velocity. The recirculation of par­
reactions and helps in predicting better product yield. The ticles promotes the mixing among the phases resulting in
prominent decrease in the gas yield may be attributed to the increased heat and mass transfer and rate of reaction.
higher sigma value for the gas formation reaction. Plotted in Fig. 8 are the molar concentrations of gas phase
Fig. 3(a) and (b) show the instantaneous volume fraction and tar as a function of reactor height. The concentration of
of gas and sand phase respectively. The bubbles start to de­ both species was zero at the bottom of the reactor and before
 Chemical Engineering Research and Design 191 (2023) 414–425 421

Fig. 4 – Time-averaged volume fraction profiles of three phases in the fluidized bed reactor.

starting to increase near the feed inlet where cellulose en­ bed region. As the tar particles moved along the bed height, it
tered in the reactor and started to pyrolyze. The rate of tar started to decompose to form gases, and its concentration
formation was higher than the gas formation reaction and decreased progressively. Since both tar and gases are elu­
hence the concentration of tar increased more rapidly in the triated out of the reactor along with the gas phase, their

Fig. 5 – Instantaneous temperature profiles for gas, sand and cellulose phase in the fluidized bed reactor.
422 Chemical Engineering Research and Design 191 (2023) 414–425

Fig. 6 – Time-averaged temperature profiles of three phases Fig. 8 – Molar concentration of tar and gases.
in fluidized bed reactor.

along the bed height was consistent in both cases. Fraction of

concentrations become constant after a certain height in the gases decreased significantly with DAEM because of the higher
reactor. activation energy for this reaction and inclusion of standard
Fig. 9 compares the mass fraction of tar and gases leaving deviation of activation energy. Fraction of char increased in the
the reactor and fraction of char accumulating inside the reactor bed region and then became nearly zero in the freeboard re­
with and without DAEM kinetics. The results presented in the gion. With DAEM, charring was favored because of high acti­
graph below validated the numerical values shown in Table 6. vation energy of reaction 2 as compared to reaction 3. Including
The fraction of tar was maximum near the feed inlet and then the distribution of activation energies has significant impact on
started to decrease as a result of tar cracking reaction. The yield pyrolysis product yield.
of tar increased with DAEM, however trend for tar formation

Fig. 7 – Time-averaged velocity vectors of the three phases (a) gas (b) cellulose (c) sand.
 Chemical Engineering Research and Design 191 (2023) 414–425 423

Fig. 9 – Mass fraction of pyrolysis product species (a) tar, (b) gas and (c) char.

5. Conclusions yield with DAEM were in better agreement with experimental

data in comparison to single activation energy model.
A multi-phase CFD model has been developed for cellulose
pyrolysis in a fluidized bed reactor. Firstly, the model was
Declaration of Competing Interest
validated by making a comparison between pyrolysis pro­
duct yields from experiments. The yield of bio-oil, char and
The authors declare that they have no known competing fi­
gases was in good agreement with the experimental data.
nancial interests or personal relationships that could have
Once the model was validated, it was integrated with the
appeared to influence the work reported in this paper.
distribution activation energy model. Kinetic parameters in­
cluding the activation energy and standard deviation were
collected from the data reported in author’s previous work. Acknowledgement
Using activation energy with the distribution function had a
profound effect on the tar yield. The yield of char and gases The author acknowledges Curtin University, University of
were also being affected significantly. The predicted product Engineering and Technology and Faculty Development
424 Chemical Engineering Research and Design 191 (2023) 414–425

Program of Higher Education Commission of Pakistan for the In the particles collision term there are two important
support to her PhD studies. parameters: ess is known as restitution coefficient and g0,ss
is the radial distribution function.
Appendix A For multiphase systems having N number of phases, the
same equation has been modified as:
Granular temperature model has been specified as phase
N dls3
property, neglecting the convective and diffusive terms, pl = l l l + l=1
2 (1 + els ) g0, ls l s l l
dl3 (38)
hence energy dissipation due to collisions and energy ex­
change will be given as: The radial distribution is calculated according to Syamlal-
2 O′brien (Syamlal et al., 1993):
12(1 ess ) g0, ss 2 3/2
s = s s s N
ds (27) 1 3( k = 1 k / dk )
g0, kl = + 2
ds dl
(1 s) (1 s ) (ds + dl ) (39)
ls = 3Kls s (28)

Rls = Kls ( vl vs ) (29) References

Kls = Momentum exchange co-efficient between the lth-
Bharadwaj, A., Baxter, L.L., Robinson, A.L., 2004. Effects of in­
fluid and sth-solid phase. traparticle heat and mass transfer on biomass devolatiliza­
Syamlal-O′brien symmetric model (Syamlal, 1987) has tion: Experimental results and model predictions. Energy
been used to specify the drag between solid particles. The Fuels 18, 1021–1031. https://doi.org/10.1021/ef0340357
model is of form: Bradbury, A.G.W., Sakai, Y., Shafizadeh, F., 1979. Kinetic model
2 /8))
for pyrolysis of cellulose. J. Appl. Polym. Sci. 23, 3271–3280.
3(1 + els )(( /2) + Cfr, ls ( (d
s s l l l
+ ds )2 g0, ls https://doi.org/10.1002/app.1979.070231112
Kls = 3 3
| vl vs|
2 ( l dl + d
s s ) Braun, R.L., Burnham, A.K., 1987. Analysis of chemical reaction
(30) kinetics using a distribution of activation energies and sim­
pler models. Energy Fuels 1, 153–161.
Stress-strain tensor for solid phase Burnham, A.K., Zhou, X., Broadbelt, L.J., 2015. Critical review of
T 2 the global chemical kinetics of cellulose thermal decomposi­
s = s µs [ vs + vs ] + s s µ . vs I tion. Energy Fuels 29, 2906–2918. https://doi.org/10.1021/acs.
3 s (31)
energyfuels.5b00350
Solid shear viscosity is obtained by adding the kinetic, Cai, J., Li, T., Liu, R., 2011. A critical study of the Miura – Maki
collisional and frictional viscosities, such as: integral method for the estimation of the kinetic parameters
of the distributed activation energy model. Bioresour.
µs = µs, kin + µs, col + µs, fr (32) Technol. 102, 3894–3899. https://doi.org/10.1016/j.biortech.
2010.11.110
Syamlal-O′brien model (Syamlal et al., 1993) has been
Cai, J., Wu, W., Liu, R., 2014. An overview of distributed activation
used to define these viscosities which results in the expres­ energy model and its application in the pyrolysis of lig­
sions as: nocellulosic biomass. Renew. Sustain. Energy Rev. 36, 236–246.
https://doi.org/10.1016/j.rser.2014.04.052
s ds s s 2
µs, kin = 1+ (1 + ess )(3ess 1) s g0, ss Curtis, L.J., Miller, D.J., 1988. Transport model with radiative heat
6(3 ess ) 5 (33) transfer for rapid cellulose pyrolysis. Ind. Eng. Chem. Res. 27,
1775–1783. https://doi.org/10.1021/ie00082a007
2
4 s Czajka, K., Kisiela, A., Moroń, W., Ferens, W., Rybak, W., 2016.
µs, col = s ds s g0, ss (1 + ess ) s
5 (34) Pyrolysis of solid fuels: thermochemical behaviour, kinetics
and compensation effect. Fuel Process. Technol. 142, 42–53.
Schaffer’s relation (Schaeffer, 1987) has been used to in­ https://doi.org/10.1016/j.fuproc.2015.09.027
corporate frictional stress for the solid phase near the Di Blasi, C., 1994. Numerical simulation of cellulose pyrolysis.
packing limit: Biomass Bioenergy 7, 87–98.
Fluent, A., 2009. ANSYS Fluent 12.0 user’s guide. Ansys Inc., pp.
ps sin
µs, fr = 1–2498. https://doi.org/10.1016/0140-3664(87)90311-2
2 I2D (35) Gunn, D.J., 1978. Transfer of heat or mass to particles in fixed and
fluidized beds. Int. J. Heat. Mass Transf. 21, 467–476.
To specify frictional pressure the model based on kinetic
Hameed, S., 2019. Kinetic Analysis and Modelling of Cellulose
theory of granular flow has been used. Pyrolysis.
Additional to shear viscosity, there is solid bulk viscosity Hameed, S., Wagh, A.S., Sharma, A., Pareek, V., Yu, Y., Joshi, J.B.,
which includes the resistance of the granular particles to 2022. Kinetic modelling of pyrolysis of cellulose using CPD
compression and expansion. The expression is given by Lun model: effect of salt. J. Therm. Anal. Calorim. 147, 9763–9777.
et al. (Lun et al., 1984): https://doi.org/10.1007/s10973-021-11192-5
Hu, M., Chen, Z., Wang, S., Guo, D., Ma, C., Zhou, Y., Chen, J.,
1/2
4 s Laghari, M., Fazal, S., Xiao, B., Zhang, B., Ma, S., 2016.
s = s ds s g0, ss (1 + ess )
3 (36) Thermogravimetric kinetics of lignocellulosic biomass slow
pyrolysis using distributed activation energy model, Fraser-
The pressure for the sth solid phase has been defined Suzuki deconvolution, and iso-conversional method. Energy
using the model given by Lun et al. (Lun et al., 1984) which is Convers. Manag. 118, 1–11. https://doi.org/10.1016/j.
the sum of the kinetic term and the pressure because of enconman.2016.03.058
collisions among particles. Jalalifar, S., Ghiji, M., Abbassi, R., Garaniya, V., 2017. Numerical
modelling of a fast pyrolysis process in a bubbling fluidized
2
ps = s s s + 2 s (1 + ess ) s g0, ss s (37) bed reactor.
 Chemical Engineering Research and Design 191 (2023) 414–425 425

Jalalifar, S., Abbassi, R., Garaniya, V., Hawboldt, K., 2018. Soria-Verdugo, A., Goos, E., García-Hernando, N., 2015. Effect of
Parametric analysis of pyrolysis process on the product yields the number of TGA curves employed on the biomass pyrolysis
in a bubbling fl uidized bed reactor. Fuel 234, 616–625. kinetics results obtained using the distributed activation en­
Kirtania, K., Bhattacharya, S., 2015. Coupling of a distributed ac­ ergy model. Fuel Process. Technol. 134, 360–371. https://doi.
tivation energy model with particle simulation for entrained org/10.1016/j.fuproc.2015.02.018
flow pyrolysis of biomass. Fuel Process. Technol. 137, 131–138. Syamlal, M., 1987. The particle-particle drag term in a multi­
https://doi.org/10.1016/j.fuproc.2015.04.014 particle model of fluidization. https://doi.org/MC/21353–2373,
Koufopanos, C.A., Papayannakos, N., Maschio, G., Lucchesi, A., NTIS/DE87006500.
1991. Modelling of the pyrolysis of biomass particles. Studies Syamlal, M., Rogers, W., O`Brien, T.J., 1993. MFIX documentation
on kinetics, thermal and heat transfer effects. Can. J. Chem. theory guide. Other Inf. PBD Dec 1993 1004, Medium: ED; Size:
Eng. 69, 907–915. 49 p. https://doi.org/METC-9411004, NTIS/DE9400087.
Lathouwers, D., Bellan, J., 2001. Yield optimization and scaling of White, J.E., Catallo, W.J., Legendre, B.L., 2011. Biomass pyrolysis
fluidized beds for tar production from biomass. Energy Fuels kinetics: a comparative critical review with relevant agri­
15, 1247–1262. https://doi.org/10.1021/ef010053h cultural residue case studies. J. Anal. Appl. Pyrolysis 91, 1–33.
Lun, C.K.K., Savage, S.B., Jeffrey, D.J., Chepurniy, N., 1984. Kinetic https://doi.org/10.1016/j.jaap.2011.01.004
theories for granular flow: inelastic particles in Couette flow and Xiong, Q., Xu, F., Ramirez, E., Pannala, S., Daw, C.S., 2016a.
slightly inelastic particles in a general flowfield. J. Fluid Mech. Modeling the impact of bubbling bed hydrodynamics on tar
140, 223–256. https://doi.org/10.1017/S0022112084000586 yield and its fluctuations during biomass fast pyrolysis. Fuel
Schaeffer, D.G., 1987. Instability in the evolution-equations de­ 164, 11–17. https://doi.org/10.1016/j.fuel.2015.09.074
scribing incompressible granular flow. J. Differ. Equ. 66, 19–50. Xiong, Q., Zhang, J., Xu, F., Wiggins, G., Stuart Daw, C., 2016b.
Shafizadeh, F., 1982. Introduction to pyrolysis of biomass. J. Anal. Coupling DAEM and CFD for simulating biomass fast pyrolysis
Appl. Pyrolysis 3, 283–305. https://doi.org/10.1016/0165- in fluidized beds. J. Anal. Appl. Pyrolysis 117, 176–181. https://
2370(82)80017-X doi.org/10.1016/j.jaap.2015.11.015
Shafizadeh, F., Bradbury, A.G.W., 1979. Thermal degradation of Xue, Q., Heindel, T.J., Fox, R.O., 2011. A CFD model for biomass
cellulose in air and nitrogen at low temperatures. J. Appl. fast pyrolysis in fluidized-bed reactors. Chem. Eng. Sci. 66,
Polym. Sci. 23, 1431–1442. https://doi.org/10.1002/app.1979. 2440–2452. https://doi.org/10.1016/j.ces.2011.03.010
070230513 Xue, Q., Dalluge, D., Heindel, T.J., Fox, R.O., Brown, R.C., 2012.
Sharma, A., Wang, S., Pareek, V., Yang, H., Zhang, D., 2015. Multi- Experimental validation and CFD modeling study of biomass
fluid reactive modeling of fluidized bed pyrolysis process. fast pyrolysis in fluidized-bed reactors. Fuel 97, 757–769.
Chem. Eng. Sci. 123, 311–321. https://doi.org/10.1016/j.ces. https://doi.org/10.1016/j.fuel.2012.02.065
2014.11.019 Zhu, H., Dong, Z., Yu, X., Cunningham, G., Umashanker, J., Zhang,
Sonobe, T., Worasuwannarak, N., 2008. Kinetic analyses of X., Bridgwater, A.V., Cai, J., 2021. A predictive PBM-DEAM
biomass pyrolysis using the distributed activation energy model for lignocellulosic biomass pyrolysis. J. Anal. Appl.
model. Fuel 87, 414–421. https://doi.org/10.1016/j.fuel.2007. Pyrolysis 157. https://doi.org/10.1016/j.jaap.2021.105231
05.004