A numerical study on biomass fast pyrolysis process: A comparison between full lumped modeling and hybrid modeling combined with CFD

A numerical study on biomass fast pyrolysis process: A comparison between full lumped modeling and hybrid modeling combined with CFD

Computers and Chemical Engineering 82 (2015) 202–215 Contents lists available at ScienceDirect Computers and Chemical Engineering journal homepage: ...

4MB Sizes 2 Downloads 75 Views

Computers and Chemical Engineering 82 (2015) 202–215

Contents lists available at ScienceDirect

Computers and Chemical Engineering journal homepage: www.elsevier.com/locate/compchemeng

A numerical study on biomass fast pyrolysis process: A comparison between full lumped modeling and hybrid modeling combined with CFD Yu Ri Lee, Hang Seok Choi ∗ , Hoon Chae Park, Ji Eun Lee Department of Environmental Engineering, Yonsei University, Wonju 220-710, Republic of Korea

a r t i c l e

i n f o

Article history: Received 19 September 2014 Received in revised form 21 June 2015 Accepted 9 July 2015 Available online 17 July 2015 Keywords: Biomass Fast pyrolysis CFD Process analysis Hybrid model

a b s t r a c t This study focuses on process modeling and simulation of a biomass fast pyrolysis system. For the simulation of the biomass fast pyrolysis process, two types of simulation models were developed: lumped model and hybrid model. Employing the above models, the effect of reaction temperature on reaction rate and final product yields were analyzed. It was found that the hybrid model exhibited a peculiar characteristic of displaying multiphase reacting flow occurring in fluidized bed. This behavior of hybrid model could have attributed for the difference in product yields of the models (hybrid and lumped). For the yields of the tar and NCG, hybrid model prediction was very consistent with the experimental results than the lumped model. However, for char yield, the results of both the lumped model and the hybrid model were close to that of experimental results. © 2015 Elsevier Ltd. All rights reserved.

1. Introduction Due to fast depletion of fossil fuels and their effects on global warming, there is a growing concern for energy and environmental issues. A clean and sustainable renewable energy is being highlighted as a substitute for fossil fuels. Among the renewable energy sources, biomass is more promising and shows great potential to replace fossil fuels as it is a sustainable and carbon neutral fuel (Francois et al., 2013; Long et al., 2013; McKendry, 2002). Fast pyrolysis is one of the existing thermochemical conversion methods for transforming biomass into useful products. This method can generate higher oil yield compared to that of other conversion methods, such as gasification or carbonization (Meier and Faix, 1999). The biocrude oil produced through fast pyrolysis of biomass can be easily transported and stored compared to solid or gaseous fuels. In fact, biocrude oil is considered to be an interesting prospective to replace petroleum crude. In fast pyrolysis, the yield of the biocrude oil depends on various factors such as reaction temperature, heating rate, condensation method, residence time, biomass type, reactor type, etc. Therefore, many experimental studies have been conducted to determine the

∗ Corresponding author at: Department of Environmental Engineering, Yonsei University, 1 Yonseidae-gil, Wonju 220-710, Gangwon-do, Republic of Korea. Tel.: +82 33 760 2485; fax: +82 33 760 2571. E-mail address: [email protected] (H.S. Choi). http://dx.doi.org/10.1016/j.compchemeng.2015.07.007 0098-1354/© 2015 Elsevier Ltd. All rights reserved.

optimum conditions to obtain oil of higher yield and good quality (Bok et al., 2012, 2013; Choi et al., 2012; Ellens and Brown, 2012; Kim et al., 2014; Salehi et al., 2011; Wilson et al., 2013). For numerical studies, basic studies using computational fluid dynamics (hereinafter abbreviated as CFD) were carried out for the fast pyrolysis of biomass (Boateng and Mtui, 2012; Bruchmuller et al., 2012; Choi et al., 2009; Choi, 2011; Mellin et al., 2013; Papadikis et al., 2009a,b,c; Xiong et al., 2013a,b, 2014; Xiong and Kong, 2014; Xue et al., 2011, 2012). Although lumped model study on mass analysis of fast pyrolysis is common, investigations on process analysis (using lumped model) is rare (Wright et al., 2010). The process analysis is essential for the scale-up study and further optimal design for its commercialization. Therefore, in this study, computational models were developed for a laboratoryscale fast pyrolysis system. Most of the computational models are derived using mathematical functions of temperature, pressure, species mass fraction, etc. Lumped models are models (hereinafter abbreviated as LM) represented by point functions. LM is normally adopted for basic process simulation of physical or chemical processes involved in various industries. However, LM has some limitations viz.: inaccurate prediction, cannot fully resolve complex three-dimensional (3D) physical phenomena occurring in each element of equipment. Therefore, in this study, CFD simulation was applied for essential equipment, such as the fast pyrolysis reactor. Then, the data from the CFD was linked with other process LSs to increase the calculation accuracy which can be referred as a hybrid model.

Y.R. Lee et al. / Computers and Chemical Engineering 82 (2015) 202–215

203

Fig. 1. Schematic diagram of the fast pyrolysis system.

There have been several previous studies to couple process model with CFD model for improving accuracy of the simulation. Firstly, through CFD model simulation of essential equipment, an empirical model is developed and applied to system level model. Alberton et al. (2009) employed CFD to simulate complex catalytic reformer and developed an empirical metamodel to predict the effectiveness factor of catalyst (Alberton et al., 2009). In this study, the effectiveness factor is a function of effective diffusivity, catalyst thermal conductivity, temperature and composition of gases. This empirical metamodel can resolve the complex geometry of the chosen catalyst and give better computational fidelity when it is applied to the full system level simulation. For a boiler of coal power plant, heat flux is calculated using CFD and 3D heat flux data is regressed to 1D function of boiler height (Edge et al., 2013). Then, this heat flux model is integrated into system level simulation of the coal power plant. Also, Fei et al. (2015) developed reduced order model for the heat transfer coefficients of boiler in an oxy-firing coal power plant considering O2 concentration and coal feeding rate. The reduced model for heat transfer coefficients has ±4% error in comparison with full CFD simulation. This model is applied to retrofit design of the power plant. Through these developed models, the calculation cost is reduced (compared with the full CFD results) however, the reduced model value has interpolation errors. For another approach, compartment model can be used for system level simulation by gPROMS. The bath cooling suspension crystallization vessel was simulated by CFD and its compartment model was developed for combined simulation using gPROMS (Kougoulos et al., 2005). However, to design the compartment structure high resolution simulation is required and special care should be taken during compartment of the flow field of the object. In addition, the real time communication model between the CFD and process model was developed to analysis rotary kiln incinerator and absorption storage tank by Marias (2003) and Moto et al. (2004). The above studies indicate that the developed model can enhance prediction accuracy. However, developed model study involves high computational cost due to different size of time-step and spatial nodes. And, for process optimization which needs much more case studies, the computational cost becomes huge. In this study, to combine process simulation with CFD, simple CFD database is composed of the interpolation technique, which also increases the calculation speed compared to the integrated simulation of the CFD and process analysis. The results of the hybrid simulation are compared with the experimental results and also with the simulation results of LM. The differences between the two simulations are examined, and the possibility of the combined

model for its actual application is evaluated for the biomass fast pyrolysis system. 2. Numerical method 2.1. Overview of the fast pyrolysis process Fig. 1 shows a schematic diagram of the fast pyrolysis system that is used for the thermochemical conversion of biomass into liquid biocrude oil. The fast pyrolysis system is composed of a dry biomass feeder, bubbling fluidized bed reactor, heaters, cyclone, three-stage condensers and electrostatic precipitator. First, the dry biomass is fed into the reactor. Then, the temperature of the biomass feedstock is rapidly increased by hot fluidizing sand and preheated nitrogen. Here, nitrogen is used for the fluidization of the solid mixture, including sand and biomass particles. Then, the biomass is decomposed into three main products, including noncondensable gases (hereinafter abbreviated as NCG), tar (volatiles) and char. After the fast pyrolysis of the biomass, the products are taken to the cyclone to eliminate char in order to increase the quality of the biocrude oil. Then, the gases and volatiles from the cyclone are issued into the three-stage condensers to collect the biocrude oil. Finally, the remaining gases are passed to the electrostatic precipitator to capture tiny biocrude oil droplets before moving into the gas sink. For reference, gPROMS software (gPROMS, 2005) was used for the process simulation of the biomass fast pyrolysis system. However, all of the equipment was modeled by the following mathematical models instead of using pre-installed models. 2.2. Computational procedure of lumped or hybrid model simulation Fig. 2 shows the flowchart of the calculation procedure for the fast pyrolysis system. For the full lumped model simulation, first, the input data are read. Then, the initial mass fractions and enthalpies are calculated using the input data. The governing equations for the fast pyrolysis reactor are calculated. The reaction rates of the reactant and product, and the mass fractions of the chemical species are computed. Then, the mathematical formulas for the cyclone, condenser, and electrostatic precipitator are calculated. The final yield of products can then be obtained. In this study, two types of process models were developed: LM and hybrid model. When the hybrid simulation is turned on, the governing equations of the LM for the fast pyrolysis reactor are not solved. Instead, the input data, including the reaction temperature and inlet gas velocity, are considered for the CFD model. Then, the searching

204

Y.R. Lee et al. / Computers and Chemical Engineering 82 (2015) 202–215

Fig. 2. Flowchart of the calculation procedure for the fast pyrolysis system.

subroutine determines appropriate values for each yield for tar, NCG and char. For reference, the yield data map was generated using the CFD calculation results and installed in the program prior to the process simulation. Finally, the interpolation operation was carried out to obtain the product yields based on the reaction temperature and inlet gas velocity. Then, these data were taken for process simulation, and the governing equations except for the fast pyrolysis reactor were solved. The implicit method was used for the temporal discretization of the governing equations. The time integration of the governing equations was continued until the convergence criteria was reached, as shown in Fig. 2. Fig. 3 represents the yield data map of the biocrude oil for various inlet gas velocities for various reaction temperatures. In general, fast pyrolysis is intended to generate more biocrude oil. Therefore, the prediction of biocrude oil yield is very important for the study. As stated above, the inlet gas velocity is an important factor in determining the biocrude oil yield as the yield and quality of biocrude oil is dependent on residence time and flow regime of the fluidized bed which are closely connected with inlet gas velocity. The yield of biocrude oil is reduced when the residence time is longer than few seconds as it leads to secondary tar cracking reactions. Therefore, the optimum inlet gas velocity must be set for the maximum biocrude oil yield. In addition, the reaction rate of the secondary reaction is much faster for a higher reaction temperature than for a lower temperature. Therefore, the biocrude oil yield significantly decreases as the reaction temperature increases, especially for longer residence time. In Fig. 3, the biocrude oil yield increases until 773 K and then decreases as the reaction temperature increases when the inlet gas velocity is 18.6 cm/s. This result

agrees well with the experimental data (Bok et al., 2012, 2013; Choi et al., 2012; Ellens and Brown, 2012; Kim et al., 2014; Salehi et al., 2011; Wilson et al., 2013). For reference, CFD can resolve the complex multiphase reacting flow fields by changing pyrolysis variables however only volume-averaged temperature can be considered for the lumped model.

Fig. 3. Yield data map of biocrude oil for various inlet gas velocities with respect to temperature.

Y.R. Lee et al. / Computers and Chemical Engineering 82 (2015) 202–215

For the hybrid simulation, the real-time communication between CFD and process analysis as discussed in introduction requires higher computational cost which is originated from different time-step size and spatial node. Also, for system optimization it needs many simulation cases. Hence, when the hybrid simulation is adopted for the optimization this needs huge computational cost. To overcome this, in the present study, before carrying out main hybrid simulation, the CFD library is made as follows. Firstly, most important flow regimes and flow variable were selected and CFD calculations were performed. Then, using interpolation technique CFD library was developed as discussed in Fig. 3. After that, using CFD library, the process analysis was carried out following the algorithm in Fig. 2. Hence, the calculation cost of hybrid simulation is the same as that of LM simulation. Total CUP time for the construction of the CFD map is around 67.8 h with 6 machines (total 8 computational nodes with 1.4 GHz). 2.3. Process analysis modeling For the system level, the following model assumptions were made in this study. • System variables such as temperature, pressure, velocity, gas and solid compositions, etc. are represented by volume-averaged value. • All gas phase chemical components in the system are treated as ideal gas. • For simplified representation, pyrolysis products such as tar, noncondensable gases and char were only considered. • Barring condensers (heat exchangers) all the concerned equipment were assumed to be adiabatic. • The tar condensation efficiency of the condensers was assumed to be 90%. • The collection efficiency of liquid droplet in the electrostatic precipitator was taken from the experimental data (see Section 2.3.3). 2.3.1. Modeling of the fast pyrolysis reactor In the present study, the lignocellulosic biomass was selected as the feedstock. The lignocellulosic biomass mainly consists of three major components such as cellulose, hemi-cellulose and lignin. The characteristics of biomass are influenced by these components. There are several kinds of reaction mechanism for the biomass pyrolysis viz.: three-step (Agrawal, 1988), BroidoShafizadeh (Bradbury and Skai, 1979), multi-step (Calonaci et al., 2010) and two-stage, semi global mechanism (Chan et al., 1985). Among them, two-stage, semi-global mechanism is widely used for the biomass pyrolysis modeling using CFD (Blasi, 1996, 1998; Papadikis et al., 2009b) for its simplicity and moderate accuracy. Hence, in the present study, two-stage, semi-global mechanism is adopted considering computational cost and accuracy. In this mechanism, tar contains all vapors from the biomass fast pyrolysis including water. This mechanism can be applied to all kinds of biomass; hence in this study the product yields were compared with several kinds of biomass. However, it should be noted that the full validation of two-stage, semi-global mechanism is not carried out till now. Also, the multi-step mechanism can be used for the modeling, considering the specific contents of three major components of the biomass. Then, light, moderate, and heavy molecules coming from the pyrolysis can be represented in the simulation. However, for some biomass, the amount of lignin content becomes negative value and this is unrealistic. Hence, the modification of the multi-step mechanism is required and furthermore high computational cost is required. Fig. 4 shows the two-stage semi-global reaction mechanism for the fast pyrolysis of the biomass. Primarily, the biomass is

205

Fig. 4. Two-stage semi-global mechanism for the fast pyrolysis of biomass (Chan et al., 1985).

decomposed into NCG1, tar and char1. Then, the tar that is generated by the primary reaction is further decomposed into NCG2 and char2 by the secondary tar-cracking reactions. Following the reaction mechanism, first, the LM of the fast pyrolysis reactor was generated. The LM was adopted for basic process simulation, although complex physical phenomena occur in the bubbling fluidized bed reactor (Bok et al., 2013). In addition, the CFD model was generated for the fast pyrolysis reactor using the same reaction mechanism. The CFD model is discussed in detail in Section 2.4. This CFD model was combined with other LMs, and the hybrid process simulation was carried out. Eqs. (1)–(5) represent the material balance and fast pyrolysis reaction rates for the LM. ˙ biomass = m

n 

˙ p+m ˙ residue m

for p = 1, 2, ..., n,

(1)

p=1

dmj dt

=

L 

vji qi for j = 1, 2, ..., m,

(2)

i=1

ji = (ji − ji ), qi = ki

n 

(3)

v

(4)

mj ji ,

j=1



ki = Ai exp −

Ei Ru T

 .

(5)

where subscripts ‘p’ and ‘residue’ indicate the fast pyrolysis prod˙ indicates the ucts and unreacted biomass residue, respectively. m mass flow rate, and n is 5 in this study. In addition, subscripts ‘j’ and ‘i’ indicate the chemical species and reaction number, respectively. Here, m = 6, and L = 5. vji and ji are stoichiometric coefficients of the fast pyrolysis reaction for the jth chemical species in the ith reaction. mj is the mass of chemical species j. ki is the reaction rate constant for the ith reaction and is expressed by the Arrhenius equation (see Eq. (5)). Ai and Ei are the pre-exponential factor and activation energy for the ith reaction rate constant, respectively. Table 1 shows the pre-exponential factors and activation energies for each reaction rate constant. For the calculation of Eqs. (1)–(5), the physicochemical properties for each material refer to CHEMKIN data (http://creckmodeling.chem.polimi.it/). For the discretization of Eq. (2), the implicit method was applied for the time integration. In this study, the time integration was carried out for over 3.76 s for the LM of the fast pyrolysis reactor for all the reaction

Table 1 Kinetics data for the two-stage semi-global mechanism (Papadikis et al., 2009c). Pre-exponential factor A1 A2 A3 A4 A5

−1

1.3 × 10 s 2.0 × 108 s−1 1.08 × 107 s−1 2.6 × 106 s−1 1.0 × 106 s−1 8

Activation energy E1 E2 E3 E4 E5

140 kJ/mol 133 kJ/mol 121 kJ/mol 108 kJ/mol 108 kJ/mol

206

Y.R. Lee et al. / Computers and Chemical Engineering 82 (2015) 202–215

temperatures, corresponding to the gas residence time in the fast pyrolyzer. In addition, the heat balance equation for the fast pyrolysis reactor is as follows.



o Hf,biomass

+ T



Tin,biomass

Cp,biomass dT + YN2



o YNCG Hf,NCG



Tin,N

To

Tout

+ YNCG To

Cp,NCG dT + YN2

2



Cp,N2 dT + Qheater =

To

Cp,N2 dT

UA , Cmin

(7)

1 . (1/htube ) + (1/hshell )

– Tube side Nu =

hdtube 0.668Re · Pr(d/L) . = 3.66 + 2/3 k 1 + 0.4[Re · Pr(d/L)]

(9)

where h is the convection heat transfer coefficient, dtube is the tube diameter, and L is the length of the condenser. k represents the fluid thermal conductivity. Re and Pr indicate the Reynolds and Prandtl numbers, respectively. – Shell side Nu =

=

hdshell k

⎧ 1.04Re0.4 Pr 0.36  0.25 ⎪ ⎪ ⎪ ⎪ ⎨ 0.71Re0.5 Pr 0.36  0.25 ⎪ ⎪ ⎪ ⎪ ⎩

for 1 ≤ Re < 500 for 500 ≤ Re < 1000

,

5

0.35ς 0.2 Re0.63 Pr 0.36  0.25

for 1000 ≤ Re < 2 × 10

0.031ς 0.2 Re0.8 Pr 0.36  0.25

for 2 × 105 ≤ Re < 2 × 106

=

Pr XT ,ς = . Pr␻ XL

(10)

where dshell is the shell diameter, and the subscript ‘␻’ indicates the wall. XL and XT represent the longitudinal tube spacing and transverse tube spacing, respectively. Finally, to determine the gas and cooling water temperatures, the following equation was employed: ˛condenser =

Ch (Th,i − Th,o ) Cc (Tc,i − Tc,o ) = . Cmin (Th,i − Tc,i ) Cmin (Th,i − Tc,i )

(11)

where subscripts ‘h’ and ‘c’ indicate the hot and cold fluids, respectively. If the heat capacity of the hot fluid was greater than

Cp,tar dT + To

(6)

Tout

+ Ychar

Cp,char dT. To

2.3.3. Miscellaneous The biomass feeder and the nitrogen vessel were modeled as material sources. Therefore, a constant biomass of 0.00028 kg/s and constant volume flow rate of N2 were supplied to the fast pyrolysis reactor, respectively. In addition, the heaters were considered as the heat source, and the constant heat flux condition was applied to the reactor as well as in the N2 vessel. The elimination efficiency of the char in cyclone was calculated employing Lapple method as defined in Eq. (12) (Theodore and Paola, 1980). For the electrostatic precipitator, the collection efficiency of the biocrude oil was calculated using Deutsch equation as defined in Eq. (15) (Deutsch, 1922) cyclone =

1 1 + (dp50 /d)



(8)

where U represents the overall heat transfer coefficient and was calculated using two convective heat transfer coefficients for the tube and shell in Eq. (8). A is the wetted area, and Cmin is the minimum heat capacity between the cold and hot fluids. To calculate the convective heat transfer coefficients of Eq. (8), following formulas were used:

Tout

+ Ytar

 o + Ychar Hf,char

2.3.2. Modeling of the condenser In this study, a shell- and tube-type heat exchanger was considered. For the modeling of the shell and tube heat exchanger, NTU (number of transfer units) method was used as follows (Bejan and Kraus, 2003; Incropera et al., 2006):

U=

 o Ytar Hf,tar

Tout

Here, the subscripts ‘tar’ and ‘NCG’ represent the tar and NCG, respectively. The subscripts ‘f’, ‘in’ and ‘out’ represent the formation enthalpy, input and output, respectively. In addition, ‘o ’ indicates the standard condition. Qheater is the external heat that is supplied from the heater to the reactor, and T is the temperature. Yi is the species mass fraction, and Cp,i is the specific heat at constant pressure for each species, Cp,i is a function of temperature.

NTU =

that of the cold fluid, Cmin becomes the heat capacity of the cold fluid. Otherwise, Cmin is the heat capacity of the hot fluid.

dp50 (␮m) = Ne =

2

(12)

,

9W 2(p − g )Ne V

1/2 × 106 ,

(13)



1 Lc L + . H b 2

(14)

where d indicates the particle diameter, and dp50 indicates the cut diameter that is defined in Eq. (13). In Eq. (13),  and W indicate the gas viscosity and cyclone inlet width, respectively. p is the particle density, and g is the gas density. In addition, Ne represents the number of gas revolutions within the cyclone, which is defined in Eq. (14). V is the gas inlet velocity in the cyclone. H indicates the cyclone inlet height. In Eq. (14), Lb and Lc indicate the cyclone cylinder height and cyclone cone height, respectively. ESP = 1 − exp



−w ·

A Q



.

(15)

where ESP means the electrostatic precipitator. w represents the migration velocity, and A represents the area of the collecting surface. Q represents the gas volume flow rate through the electrostatic precipitator. In this study, the experimental result of Jedrusik et al. (2001) was applied to calculate w. 2.4. CFD modeling For the CFD modeling, both the gas and solid phase were assumed to be continuum. All the solid particles have the same characteristics such as physicochemical properties and effective diameter. The assumptions and approximation to derive the CFD governing equations can be found in MFIX documentation theory guide (Syamlal et al., 1993). For the hybrid simulation, the fast pyrolyzer was modeled by CFD method and linear interpolation technique was applied to develop CFD library as discussed in Section 2.2. Eqs. (16)–(23) show the governing equations for the gas–solid multiphase reacting flow (Wachem et al., 2001) (Table 2). Eqs. (16) and (20) are the continuity equations. Subscripts ‘g’ and ‘sj’ indicate the gas and solid phase. v represents the velocity vector, and ε indicates the volume fraction for the gas phase or solid phase. In addition, R indicates the reaction rate for the gas and solid phase

Y.R. Lee et al. / Computers and Chemical Engineering 82 (2015) 202–215

207

Table 2 Governing equations. Fluid phase Continuity equation

∂ (εg g ) + ∇ · (εg g g ) = ∂t

n 



(16)

˛=1

Momentum equation

Energy equation

∂ (εg g g ) + ∇ · (εg g g g ) = −εg ∇ Pg + ∇ ·  g + ∂t εg g Cpg (

∂Tg + g · ∇ Tg ) = −∇ qg + ∂t

n 

n 

Fgsj (␯sj − g ) + εg g g −

j=1

n 

RMgsj ( gsj sj + gsj g )

(17)

j=1

gsj (Tsj − Tg ) − Hrg

(18)

j=i

Species equation Solid phase Continuity equation

∂ (εg g Yg˛ ) + ∇ · (εg g Yg˛ g ) = ∇ · Dg˛ ∇ Yg˛ + Rg˛ ∂t

  ∂ (εsj sj ) + ∇ · εsj sj sj = ∂t

n 

Rsjˇ

(20)

ˇ=1

Momentum equation

(19)

∂ (εsj sj sj ) + ∇ · (εsj sj sj sj ) = −εsj ∇ Pg + ∇ · S sj − Fgsj (sj − g ) + ∂t

 n

Fsksj (sk − sj ) + εsj sj g + RMgsj ( gsj sj + gsj g )

k=1



n 

RMsksj ( sksj sk + sksj sk )

(21)

k=1

Energy equation Species equation

∂Tsj + sj · ∇ Tsj ) = −∇ qsj − gsj (Tsj − Tg ) − Hrsj ∂t ∂ (εsj sj Ysjˇ ) + ∇ · (εsj sj Ysjˇ v sj ) = ∇ · Dsjˇ ∇ Ysjˇ + Rsjˇ ∂t

εsj sj Cpsj (

that are produced by the chemical reaction as discussed in Section 2.3.1. Eqs. (17) and (21) present the momentum equations. g is the stress tensor for the gas phase, and Ssj is the jth stress tensor for the solid phase. Fgsi is the interphase drag force coefficient between the gas and solid phases, and Fsksj is the interphase drag force coefficient between the solid phases. In addition, MRml indicates the mass transfer term from the mth phase to the lth phase. If MRml is  = 1 − . Eqs. (18) not equal to zero, ml is equal to one, and ml ml and (22) represent the energy equations. T represents the temperature for each phase. q is the heat flux from the conduction of the gas phase or solid phase.  gsi is the coefficient of the heat transfer between the gas and solid phases. In addition, Hrg and Hrsj are the heat of reaction for the gas and solid phases, respectively. Eqs. (19) and (23) represent the species equations. Y represents the mass fraction for the chemical species, and D stands for the coefficient of the molecular diffusion for each chemical species. In this study, a bubbling fluidized bed was used for the fast pyrolysis of biomass. The bubbling fluidized bed reactor is the most important equipment in the entire fast pyrolysis system. Therefore, the bubbling fluidized bed reactor was modeled and simulated by CFD for the hybrid process simulation. The CFD calculation was performed using MFIX code (Syamlal et al., 1993). Then, the CFD data were linked to other process models during the hybrid process simulation to increase the fidelity of the prediction. To combine the process models with the CFD model, the simple CFD database set (see Fig. 3) was developed as it can increase the total calculation speed. For the CFD simulation of the multiphase flow in the reactor, an Eulerian-Granular approach was applied for the simplicity of the calculation compared to other particle methods considering both the continuum gas and discrete solid phases. The calculation conditions of the computational domain are given in Fig. 5 and Table 3, respectively. To validate the present CFD method, a bubbling fluidized bed was calculated, and the results were compared with the experimental results (Hulme et al., 2005) as like earlier studies (Choi, 2011, 2012, 2015). The hydrodynamics and heat transfer behavior which are very important for the simulation of mixing and chemical reaction of gas–solid multiphase reacting flow are analyzed in the bubbling fluidized bed. In this study, the

(22) (23)

bubbling fluidized bed that was used in the experiments (Hulme et al., 2005) was considered for the fast pyrolysis simulation. Considering inlet temperature and corresponding flow rate, the maximum inlet Reynolds number of the pyrolyzer was found to be around 300 and the gas flow was observed to be laminar. In bubbling fluidized bed (dense bed), the gas phase turbulence is mostly damped out because of the inertia of the dense particles (Ding and Gidaspow, 1990; Enwald et al., 1997; Mckeen and Pugsley, 2003; Plain et al., 2001), Hence, the governing equation for the gas phase can be approximated as laminar equation. However, for the solid phase, granular temperature (intensity of solid fluctuation)

Fig. 5. Computational domain.

208

Y.R. Lee et al. / Computers and Chemical Engineering 82 (2015) 202–215

Table 3 Calculation conditions. Computational domain (2D) Length (x) Height (y) Grid allocation (x,y) Boundary conditions Gas inlet Biomass inlet Outlet Wall

10 cm 70 cm 20 × 210 Uinlet = 18.6 cm/s, Tinlet = 673–873 K ˙ inlet = 0.00028 kg/s, Tinlet = 300 K m Neumann No-slip, Twall = 753 K

was considered and the effect of solid turbulence was included by solving the governing equation of granular temperature which was derived from gas kinetic theory. However, for the simulation of dilute gas–solid flow such as circulating fluidized bed flow the gas turbulence becomes more important and the effect should be included in the governing equation. Fig. 6 shows the results for grid independency test. In Fig. 6, the mass flow rate of tar at the center location of the bed was examined based on the number of grids. The mass flow rate of tar was saturated over grid number of 4200. Hence, considering computational cost, the minimum number of grid was selected for the present CFD simulation as indicated by the red arrow in Fig. 6. 3. Results and discussion 3.1. Process analysis with the lumped model

flowrate (kg/s) M Mass fl t off ttar (k /)

In this study, the data from the LM process analysis and hybrid process analysis were compared. Fig. 7 presents the reaction rates for the fast pyrolysis products with respect to reaction temperature from 673 to 873 K, and the reaction rates were obtained from the LM calculation. The reaction rate of the tar initially increases until the magnitude becomes maximum then, decreases for a longer reaction time because of the secondary tar-cracking reactions (Papadikis et al., 2009c). At high reaction temperatures the reaction rate of tar reaches maximum very quickly. The maximum reaction rate was attained at approximately 0.5 s for 873 K. However, for the lower reaction temperature, a longer reaction time was required to reach the maximum reaction rate owing to its lower reaction constant. In addition, the peak value increases with the increasing reaction temperature. The maximum reaction rate at 873 K was 1.01 e−4 kg/s, and the maximum reaction rate at 673 K was 6.35 e−5 kg/s. Thus, the maximum reaction rate at 873 K was 59.1% higher than that at 673 K. However, the reaction rates of NCG1 and char1 continuously increase until they reach the saturated value for all of the reaction temperatures. The reaction time that is required to reach the maximum value decreases with increasing reaction temperature. 5.0e-4

2.5e-4

0.0 0

60000

12000

18000

Number of ggrid Fig. 6. Grid independency test. (For interpretation of the references to color in this figure text, the reader is referred to the web version of this article.)

This means that the reaction rate is getting faster with increasing reaction temperature. In addition, the maximum reaction rate of char1 decreases as the reaction temperature increases, indicating that the total amount of char1 decreases as the reaction temperature increases. This result matches well with the experimental fast pyrolysis results (Bok et al., 2012, 2013; Choi et al., 2012; Ellens and Brown, 2012; Kim et al., 2014; Salehi et al., 2011; Wilson et al., 2013). For the secondary reactions of NCG2 and char2, the curves of the reaction rates exhibited a similar pattern to those of the primary reaction. However, the maximum rate of char2 increases with increasing reaction temperature. The secondary reaction rates are functions of the temperature as well as the mass fraction of the tar that is produced by the primary reaction, indicating that intensive reaction regions for the primary reaction and secondary reaction are locally different in the bubbling fluidized bed reactor, and resulting in a significant difference in the final product yields between the full LM and hybrid simulations. For the LM simulation, as previously mentioned, the integration time is was 3.76 seconds, which corresponds to the gas residence time for the fast pyrolyzer, as indicated by the red arrow in Fig. 7. 3.2. CFD simulation for the fast pyrolysis reactor Fig. 8 presents the instantaneous gas volume fraction with respect to reaction temperature. It was noticed that tiny bubbles were generated at the bottom of the bed and grew along the height of the bed. This was due to their coalescence as well as the decreased static pressure following positive streamwise direction. Therefore, the bubble size increases along the top of the bed (Choi, 2011; Plain et al., 2001; Taghipour et al., 2005; Zimmermann and Taghipour, 2005). Fig. 9 represents the instantaneous solid particle velocity with respect to reaction temperature. In this figure, the subscript ‘s1’ indicates the sand, namely the fluidized media. Figs. 8 and 9 present the CFD results at the same instant. In Fig. 9, it was interesting to note that the magnitude of the solid velocity was higher near the top of the bed where larger bubbles are located. As the bubble rises and grows, the bubble rising velocity increases; consequently, the solid velocity near the bubbles increases as well. In addition, the solid velocity near the walls has greater negative values (Hamzehei et al., 2010). These bubbles and the associated solid flow play an important role in the mixing of solid particles, heat transfer and final fast pyrolysis reaction. Fig. 10 presents the CFD results for the instantaneous contours of reaction rates with respect to reaction temperature. The biomass feeder location ranges from y = 20 cm to y = 22 cm on the left sidewall of the pyrolyzer. Considering Fig. 8, the solid volume fraction was greater near the walls; therefore, the solid biomass volume fraction was greater at those locations. Accordingly, the reaction rate of the tar shows a higher magnitude near the wall as the tar is generated from biomass feedstock. The heat transfer rate is more increased between solid–solid interfaces than that of solid–gas interface (Choi, 2011, 2012, 2015). Hence, the higher heat transfer rate between solid–solid interfaces increases the emulsion temperature and finally enhances the primary degradation rate of the biomass. In addition, the tar reaction rate was found to be lower near the gas bubble and very low in the freeboard due to the absence of the reactant. For the secondary reaction of the NCG2, the contour of the reaction rate was very different from that of the primary reaction. The secondary reaction occurs in the bubbles as well as the freeboard because the tar, which is the reactant for the secondary cracking reaction, is in the gas phase. The reaction rates of secondary products are a function of the tar mass fraction as well as the temperature. Therefore, the reaction rate of the NCG2 has a higher value where the tar is actively produced (Choi, 2011). To investigate these phenomena, the time-averaged quantities were analyzed as below.

Y.R. Lee et al. / Computers and Chemical Engineering 82 (2015) 202–215

209

Fig. 7. Reaction rates for fast pyrolysis products with respect to reaction temperature. (For interpretation of the references to color in this figure text, the reader is referred to the web version of this article.)

Fig. 11a shows the distribution of the time-averaged gas volume fraction at 823 K for three different streamwise locations. The gas volume fraction near the wall was observed to be lower than the gas volume fraction in the middle of the reactor, indicating that the solid volume fraction is higher near the wall (Taghipour et al., 2005). The gas volume fraction increases as the streamwise

position increases because of the increasing bubble size. Finally, the gas volume fraction dramatically increases at the freeboard (approximately y = 50 cm). Fig. 11b shows the streamwise solid particle velocity at 823 K for three different streamwise locations. In this figure, subscript ‘s1’ indicates the sand. The solid particle velocity in the middle of the reactor was positive for y = 20 cm and

210

Y.R. Lee et al. / Computers and Chemical Engineering 82 (2015) 202–215

Fig. 8. Instantaneous gas volume fraction with respect to reaction temperature.

Fig. 9. Instantaneous solid velocity vector with respect to reaction temperature.

Y.R. Lee et al. / Computers and Chemical Engineering 82 (2015) 202–215

Fig. 10. Contours of the instantaneous reaction rates with respect to reaction temperature.

211

212

Y.R. Lee et al. / Computers and Chemical Engineering 82 (2015) 202–215

Fig. 11. Distributions of the time-averaged gas volume fraction and solid velocity at 823 K.

y = 40 cm. In contrast, the solid particle velocity at the wall was negative for y = 20 cm and y = 40 cm. These results correspond to the previous instantaneous results. The negative velocities can be explained by the following physical phenomena. In bubbling fluidized bed, particles go upward at the middle of the bed and fall down near the walls increasing the particle volume fraction. These physical phenomena have been reported in many experiments (Pain et al., 2001). Also, the solid particle velocity was found to be negative at the top of the bed (y = 50 cm), indicating that the solid particle falls down at the freeboard. Fig. 12 presents the distribution of the time-averaged reaction rates of the tar and NCG2 for three different streamwise locations. The reaction temperature was noticed to be 823 K. In Fig. 12a, the reaction rate on the left side was higher than the reaction rate on the right side due to the larger biomass mass fraction on the left side because the biomass is inserted at the left inlet. The reaction rate of the tar is a function of the biomass mass fraction as well as the temperature. As a result, the reaction rate at the left side

Fig. 12. Distribution of time-averaged reaction rates at 823 K.

of y = 20 cm was highest in Fig. 12a. Overall, the time-averaged tar reaction rates near the wall were higher than those of the middle for all of the y positions. This result can be explained as follows. Considering Fig. 10, the solid volume fraction was higher, and the negative solid flow was active near the wall; therefore, solid mixing and consequent heat transfer could be more vigorous than the middle of the reactor, leading to the intensive primary reaction near the wall. With increasing y position, the tar reaction rate decreases. Fig. 12b presents the time-averaged NCG2 secondary reaction rate. Here, the tar is a reactant for the secondary reaction that produces the NCG2 and char2. Along the outlet of the fast pyrolysis reactor, the overall magnitude of the time-averaged NCG2 reaction rate increases. Near the wall, the NCG2 reaction rate was found to be higher than that of the middle of the reactor. The NCG2 reaction rate was noticed to be higher at the freeboard (y = 50 cm), while the tar reaction rate was observed to be lower. In addition, the distribution of the time-averaged tar reaction rate was still higher near the left side than that of the right side, however the distribution of the timeaveraged NCG2 reaction rate indicates the nearly symmetry for all

Yield of biocrude oil Maximum yield of biocrude oil (wt/wt)

Y.R. Lee et al. / Computers and Chemical Engineering 82 (2015) 202–215

213

1.2 1.0 0.8 0.6 0.4 0.2 0.0 673

723

773

823

873

Reaction temperature(K) Hybrid CFD, LM , Upper and lower boundaries of experimental data Fig. 14. Yield of biocrude oil with respect to reaction temperature (- (Beis et al., 2002), ⊕, , , , 䊐, (Blasi, 2009), 夽, × (Bok, 2013), 䊎 (Onay and Kockar, 2004)).

reaction rate at y = 40 cm at 873 K is lower than those of other temperature conditions. Fig. 13b presents the time- and line-averaged NCG2 secondary reaction rate with the changing reaction temperature. The NCG2 reaction rate steadily increases as the reaction temperature increases. The overall magnitude of the NCG2 reaction rate is greater toward the freeboard. It is interesting to see that these phenomena affect the final yields of the products however, it cannot be resolved by the LM simulation. Therefore, there should be differences in the final yields between the full LM and hybrid CFD simulations. This possibility is discussed in using the figures below. 3.3. Comparison of the final yield between the lumped and hybrid models

of the streamwise locations. This result occurs because the timeaveraged NCG2 reaction rate is a function of the tar mass fraction, and the tar is distributed more evenly than solid biomass over the reactor due to vigorous gas bubble motion. This complex gas–solid flow and consequent reaction phenomena greatly influence on the final yield of the biocrude oil. Fig. 13 presents the time- and line-averaged reaction rates of the tar and NCG2 with respect to reaction temperature. In Fig. 13a, the reaction rates increase until 823 K as the temperature increases for all of the streamwise locations because the reaction rate depends on the reaction rate constant, which is exponentially proportional to the reaction temperature. At greater than 823 K, the reaction rates decrease. However, the reaction rate at y = 20 cm still increases at temperatures greater than 823 K. The biomass was actively converted into tar around y = 20 cm, and the tar reaction rate at y = 20 cm was found to be maximum at 873 K. Therefore, the biomass concentration becomes rapidly lower at y = 20 cm when compared to other temperatures. Thus, the tar

1.2 Yield of NCG Maximum yield of NCG (wt/wt)

Fig. 13. Time- and line-averaged reaction rates with respect to reaction temperature.

Figs. 14–16 present the yields of fast pyrolysis products with respect to reaction temperature for the LM, hybrid analysis and experimental data (Beis et al., 2002; Blasi, 2009; Bok, 2013; Onay and Kockar, 2004). The biocrude oil yield was calculated by taking into account the biocrude oil yield collected at the condensers and electrostatic precipitator. The NCG and char were calculated from the gas and char sinks as sum of all the products from both the primary and secondary reactions. In addition, experimental results indicate the upper and lower limits for the various types of biomasses and pyrolyzers and the consequent reaction conditions. For the LM, the reaction rates of fast pyrolysis products change with the reaction time (see Fig. 7). In addition, the gas residence time is an important factor when determining the biocrude oil

1.0 0.8 0.6 0.4 0.2 0.0 673

723

773

823

873

Reaction temperature(K) Hybrid CFD, LM , Upper and lower boundaries of experimental data Fig. 15. Yield of NCG with respect to reaction temperature (- (Beis et al., 2002), ⊕, , , , 䊐, (Blasi, 2009), 夽, × (Bok, 2013), 䊎 (Onay and Kockar, 2004)).

214

Y.R. Lee et al. / Computers and Chemical Engineering 82 (2015) 202–215

Yield of char Maximum yield of char (wt/wt)

1.2 1.0 0.8 0.6 0.4 0.2 0.0 673

723

773

823

873

Reaction temperature(K) Hybrid CFD , LM , Upper and lower boundaries of experimental data Fig. 16. Yield of char with respect to reaction temperature (- (Beis et al., 2002), ⊕, , , , 䊐, (Blasi, 2009), 夽, × (Bok, 2013), 䊎 (Onay and Kockar, 2004)).

yield because of the secondary tar-cracking reaction. Therefore, in this process analysis, the total reaction time is selected according to the gas residence time of the CFD. In Fig. 14, the biocrude oil yield first increases as the temperature increases until the yields peak. Then, the biocrude oil yield decreases as the temperature further increases because of active secondary tar-cracking reaction (Papadikis et al., 2009c). For the experimental data for analysis of pyrolysis gas with gas chromatography/mass spectrometry, the high molecular weight chemical components in the biocrude oil such as levoglucosan or syringol are decreased with the high reaction temperature (Branca et al., 2003; Demirbas, 2007). The calculation result of the hybrid model was found to be much closer to the experimental result than that of the LM. In particular, the hybrid result was closer to the upper limit of the experimental result, however the LM result indicates a smaller biocrude oil yield except for the temperature range from 753 to 803 K. In Fig. 15, the NCG yield increases as the temperature increases for both the LM and hybrid model results. The LM and hybrid model results present a comparable trend with the experimental data, nevertheless at lower temperatures, the LM underestimates the NCG yield. In Fig. 16, the char yields of the LM and hybrid analyses present a similar trend as the experimental yields. The char yields of the hybrid model and LM decrease as the temperature increases, and these results correspond well to the experiment data. Based on these overall results, the hybrid model agrees well with the real phenomena than that of the LM, which presents a smaller value compared to that of the hybrid model and experimental results. 4. Conclusions In this study, the process analysis was carried out for the fast pyrolysis of the biomass. Two models, the LM and hybrid model combined with the CFD, were studied and their results were compared. The results of comparing the LM and hybrid model simulations are as follows: (1) For the reaction rate of tar, the LM presents increasing and decreasing patterns as the reaction time progresses. The other reaction rates, such as the NCG1, NCG2, char1, and char2, increase and then gets saturated. With the increasing reaction temperature, the maximum reaction rate increases in all of the cases except for char1. (2) From the CFD results, for the primary reaction, the instantaneous higher reaction rates are located at the region of the instantaneous higher solid volume fraction. However, the secondary reaction occurs at the region of the higher gas volume fraction, such as the bubbles and freeboard.

(3) The time-averaged reaction rate of the tar is much higher in the wall region than that in the middle of the reactor because the solid volume fraction and magnitude of the solid velocity are higher at this location. In contrast, the distribution of the time-averaged NCG2 reaction rate shows y-axis symmetry. Because the reactant tar for the secondary reaction is the gas phase, it can be more freely convected to the entire reactor than the solid biomass. The time- and line-averaged tar reaction rate increases with the increasing reaction temperature for most of the streamwise locations except for the reaction temperature at 873 K. In addition, the time- and line-averaged NCG2 reaction rate increases with the increasing reaction temperature. (4) Finally, there is a significant difference in the magnitude of the product yield between the LM and hybrid model simulations, even though the product yield patterns of the two models are similar to the experimental data. For the LM, the magnitude of the biocrude oil yield was far below the lower boundary of the experimental data except for the temperature range from 753 to 803 K. In contrast, the results of the hybrid model for biocrude oil yield are located within the lower and upper boundaries of the experimental data. Thus, through the combination of the CFD and process models, it can be concluded that the calculation accuracy of the process simulation can be greatly improved. Acknowledgment This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea Government (MSIP) (No. 2014R1A2A2A03003812). References Agrawal RK. Kinetics of reaction involved in pyrolysis of cellulose. I. The three reaction model. Can J Chem Eng 1988;66:403–12. Alberton AL, Schwaab M, Fontes CE, Bittencourt RC, Pinto JC. Hybrid modeling of methane reformers. 1. A metamodel for the effectiveness factor of a catalyst pellet with complex geometry. Ind Eng Chem Res 2009;48:9369–75. Beis SH, Onay O, Kockar OM. Fixed-bed pyrolysis of safflower seed: influence of pyrolysis parameters on product yields and compositions. Renew Energ 2002;26:21–32. Bejan A, Kraus AD. Heat transfer handbook. 1st ed. New Jersey: John Wiley & Sons, Inc.; 2003. Blasi CD. Heat momentum and mass transport through a shrinking biomass particle exposed to thermal radiation. Chem Eng Sci 1996;51:1121–32. Blasi CD. Comparison of semi-global mechanisms for primary pyrolysis of lignocellulosic fuels. J Anal Appl Pyrol 1998;47:43–64. Blasi CD. Combustion and gasification rates of lignocellulosic chars. Prog Energ Combust 2009;35:121–40. Boateng AA, Mtui PL. CFD modeling of space-time evolution of fast pyrolysis products in a bench-scale fluidized-bed reactor. Appl Therm Eng 2012;33–34:190–8. Bok JP, Choi HS, Choi YS, Park HC, Kim SJ. Fast pyrolysis of coffee grounds: characteristics of product yields and biocrude oil quality. Energy 2012;47:17–24. Bok JP, Choi HS, Choi JW, Choi YS. Fast pyrolysis of Miscanthus sinensis in fluidized bed reactors: characteristics of product yields and biocrude oil quality. Energy 2013;60:44–52. Bradbury AGW, Skai Y, Shafizadeh F. A kinetic model for pyrolysis of cellulose. J Appl Polym Sci 1979;23:3271–80. Branca C, Giudicianni P, Blasi CD. GC/MS characterization of liquids generated from low temperature pyroysis of wood. Ind Eng Chem Res 2003;42:3190–202. Bruchmuller J, Wachem BGM, Gu S, Luo KH, Browm RC. Modeling the thermochemical degradation of biomass inside a fast pyrolysis fluidized bed reactor. AIChE J 2012;58:3030–42. Calonaci M, Grana R, Hemings E, Bozzano G, Dente M, Ranzi E. Comprehensive kinetic modeling study of bio-oil formation from fast pyrolysis of biomass. Energy Fuel 2010;24:5727–34. Chan WCR, Kelbon M, Krieger BB. Modelling and experimental verification of physical and chemical processes during pyrolysis of a large biomass particle. Fuel 1985;64:1505–13. Choi HS. The fast pyrolysis characteristics of lignocellulosic biomass in a bubbling fluidized bed reactor. J Comput Fluid Eng 2011;16:94–101. Choi HS. The characteristics of gas–solid flow and wall heat transfer in a fluidized bed reactor. Heat Mass Transf 2012;48:1513–24. Choi HS. Wall heat transfer of a small blunt body immersed in a fluidized bed. Numer Heat Transf A-Appl 2015;68:288–311. Choi HS, Choi YS, Kim SJ. Numerical study of fast pyrolysis of woody biomass in a gravity-driven reactor. Environ Prog Sustain Energy 2009;28:418–26.

Y.R. Lee et al. / Computers and Chemical Engineering 82 (2015) 202–215 Choi HS, Choi YS, Park HC. Fast pyrolysis characteristics of lignocellulosic biomass with varying reaction conditions. Renew Energ 2012;42:131–5. Demirbas A. The influence of temperature on the yields of compounds existing in bio-oils obtained from biomass samples via pyrolysis. Fuel Process Technol 2007;88:591–7. Deutsch W. Bewegung und Ladung der Elektrizitatstrager im Zylinderkondensator. Ann Phys 1922;68:335–44. Ding J, Gidaspow D. A bubbling fluidization model using kinetic theory of granular flow. AIChE J 1990;36:523–38. Edge PJ, Heggs PJ, Pourkashanian M, Stephenson PL. Integrated fluid dynamicsprocess modelling of a coal-fired power plant with carbon capture. Appl Therm Eng 2013;60:242–50. Ellens CJ, Brown RC. Optimization of a free-fall reactor for the production of fast pyrolysis bio-oil. Bioresour Technol 2012;103:374–80. Enwald H, Peirano E, Almstedt AE. Eulerian two-phase theory applied to fluidization. Int J Multiph Flow 1997;22:21–66. Fei Y, Black S, Szuhanszki J, Ma L, Ingham DB, Stanger PJ, et al. Evaluation of the potential of retrofitting a coal power plant to oxy-firing using CFD and process co-simulation. Fuel Process Technol 2015;131:45–58. Francois J, Abdelouahed L, Mauviel G, Patisson F, Mirgaux O, Rogaume C, et al. Detailed process modeling of a wood gasification combined heat and power plant. Biomass Bioenerg 2013;51:68–82. gPROMS. Introductory user guide. Process System Enterprise Ltd; 2005. Hamzehei M, Rahimzadeh H, Ahmadi G. Computational and experimental study of heat transfer and hydrodynamics in a 2D gas–solid fluidized bed reactor. Ind Eng Chem Res 2010;49:5110–21. Hulme I, Clavelle E, Lee L, Kantzas A. CFD modeling and validation of bubble properties for a bubbling fluidized bed. Ind Eng Chem Res 2005;44:4254–66. Incropera FP, Dewitt DP, Bergman TL, Lavine AS. Introduction to heat transfer. 5th ed. Toronto: John Wiley & Sons, Inc.; 2006. Jedrusik M, Gajewski JB, Swierczok AJ. Effect of the particle diameter and corona electrode geometry on the particle migration velocity in electrostatic precipitators. J Electrost 2001;51–52:245–51. Kim JW, Lee HW, Lee IG, Jeon JK, Ryu CK, Park SH, et al. Influence of reaction conditions on bio-oil production from pyrolysis of construction waste wood. Renew Energ 2014;65:41–8. Kougoulos E, Jones AG, Kaczmar MW. CFD modelling of mixing and heat transfer in batch cooling crystallizers: aiding the development of a hybrid predictive compartmental model. Chem Eng Res Des 2005;83:30–9. Long H, Li X, Wang H, Jia J. Biomass resources and their bioenergy potential estimation: a review. Renew Sust Energ Rev 2013;26:344–52. Marias F. A model of a rotary kiln incinerator including processes occurring within the solid and the gaseous phases. Comput Chem Eng 2003;27:813–25. Mckeen T, Pugsley T. Simulation and experimental validation of a freely bubbling bed of FCC catalyst. Powder Technol 2003;129:139–52. McKendry P. Energy production from biomass (part 1): overview of biomass. Bioresour Technol 2002;83:37–46. Meier D, Faix O. State of the art of applied fast pyrolysis of lignocellulosic materials – a review. Bioresour Technol 1999;68:71–7. Mellin P, Zhang Q, Kantarelis E, Yang W. An Euler–Euler approach to modeling biomass fast pyrolysis in fluidized-bed reactors – focusing on the gas phase. Appl Therm Eng 2013;58:344–53.

215

Moto JPB, Esteves IAAC, Abadi MR. Dynamic modelling of an adsorption storage tank using a hybrid approach combining computational fluid dynamics and process simulation. Comput Chem Eng 2004;28:2421–31. Onay O, Kockar OM. Fixed-bed pyrolysis of rapeseed (Brassica napus L.). Biomass Bioenerg 2004;26:289–99. Pain CC, Mansooradeh S, Oliveira CRE. A study of bubbling and slugging fluidized beds using the two-fluid granular temperature model. Int J Multiph Flow 2001;27:527–51. Papadikis K, Gu S, Bridgwater AV. CFD modelling of the fast pyrolysis of biomass in fluidised bed reactors: modelling the impact of biomass shrinkage. J Chem Eng 2009a;149:417–27. Papadikis K, Gu S, Bridgwater AV. CFD modelling of the fast pyrolysis of biomass in fluidised bed reactors. Part B: Heat, momentum and mass transport in bubbling fluidised beds. Chem Eng Sci 2009b;64:1036–45. Papadikis K, Gu S, Bridgwater AV, Gerhauser H. Application of CFD to model fast pyrolysis of biomass. Fuel Process Technol 2009c;90:504–12. Plain CC, Mansoorzadeh S, Oliveira CRE. A study of bubbling and slugging fluidised beds using the two-fluid granular temperature model. Int J Multiph Flow 2001;27:257–551. Salehi E, Abedi J, Harding T. Bio-oil from sawdust: effect of operating parameters on the yield and quality of pyrolysis products. Energy Fuel 2011;25:4145–54. Syamlal M, Rogers W, O’Btien TJ. MFIX documentation theory guide. U.S. Department of Energy. Office of Fossil Energy. Morgantown Energy Technology Center; 1993. Taghipour F, Ellis N, Wong C. Experimental and computational study of gas–solid fluidized bed hydrodynamics. Chem Eng Sci 2005;60:6857–67. Theodore L, Paola VD. Predicting cyclone efficiency. J Air Pollut Control Assoc 1980;30:1132–3. Wachem BGM, Schouten JC, Bleek CM, Krishna R, Sinclair JL. Comparative analysis of CFD models of dense gas–solid systems. AIChE J 2001;47:1035–51. Wilson DM, Dalluge DL, Rover M, Heaton EA, Brown RC. Crop management impacts biofuel quality: influence of switchgrass harvest time on yield, nitrogen and ash of fast pyrolysis products. Bioenerg Res 2013;6:103–13. Wright MM, Daugaard DE, Satrio JA, Brown RC. Techno-economic analysis of biomass fast pyrolysis to transportation fuels. Fuel 2010;89:S2–10. Xiong Q, Kong SC. Modeling effects of interphase transport coefficients on biomass pyrolysis in fluidized beds. Powder Technol 2014;262:96–105. Xiong Q, Kong SC, Passalacqua A. Development of a generalized numerical framework for simulating biomass fast pyrolysis in fluidized-bed reactors. Chem Eng Sci 2013a;99:305–13. Xiong Q, Aramideh S, Kong SC. Modeling effects of operating conditions on biomass fast pyrolysis in bubbling fluidized bed reactors. Energy Fuel 2013b;27: 5948–56. Xiong Q, Aramideh S, Kong SC. Assessment of devolatilization schemes in predicting product yields of biomass fast pyrolysis. Environ Prog Sustain Energy 2014;33:756–61. Xue Q, Heindel TJ, Fox RO. A CFD model for biomass fast pyrolysis in fluidized-bed reactors. Chem Eng Sci 2011;66:2440–52. Xue Q, Dalluge D, Heindel TJ, Fox RO, Brown RC. Experimental validation and CFD modeling study of biomass fast pyrolysis in fluidized-bed reactors. Fuel 2012;97:757–69. Zimmermann S, Taghipour F. CFD modeling of the hydrodynamics and reaction kinetics of FCC fluidized-bed reactors. Ind Eng Chem Res 2005;44:9818–27.