Energy Fuels 2010, 24, 2133–2145 Published on Web 02/22/2010
: DOI:10.1021/ef901497b
Numerical Comparison of the Drag Models of Granular Flows Applied to the Fast Pyrolysis of Biomass K. Papadikis,† S. Gu,*,† A. Fivga,‡ and A. V. Bridgwater‡ †
School of Engineering Sciences, University of Southampton, Highfield, Southampton, SO17 1BJ, United Kingdom, and ‡ School of Engineering and Applied Science, Aston University, Aston Triangle, Birmingham B4 7ET, United Kingdom Received December 8, 2009. Revised Manuscript Received February 9, 2010
The paper presents a comparison between the different drag models for granular flows developed in the literature and the effect of each one of them on the fast pyrolysis of wood. The process takes place on an 100 g/h lab scale bubbling fluidized bed reactor located at Aston University. FLUENT 6.3 is used as the modeling framework of the fluidized bed hydrodynamics, while the fast pyrolysis of the discrete wood particles is incorporated as an external user defined function (UDF) hooked to FLUENT’s main code structure. Three different drag models for granular flows are compared, namely the Gidaspow, Syamlal O’Brien, and Wen-Yu, already incorporated in FLUENT’s main code, and their impact on particle trajectory, heat transfer, degradation rate, product yields, and char residence time is quantified. The Eulerian approach is used to model the bubbling behavior of the sand, which is treated as a continuum. Biomass reaction kinetics is modeled according to the literature using a two-stage, semiglobal model that takes into account secondary reactions.
reported in the literature6,9-14 to account for the momentum exchange between the gas and solid phases and have been applied to dispersed two-phase flow simulations with great success. Most of these models have been studied and compared with experimental data in several cases, such as the studies of Kafui et al.,4 Li and Kuipers,5 and Taghipour et al.15 The work of Kafui et al.,4 concerned the simulation of pressure drop-gas superficial velocity using two models (pressure gradient force model, fluid density based buoyancy model) in an EulerianLagrangian framework, using the continuous single-function correlation of Di Felice13 to account for the drag forces in dense and dilute regimes. The study showed significant differences in the pressure drop-superficial gas velocity profiles in the fixed bed regime and corresponding significant differences in the prediction of minimum fluidization velocity. The work of Li and Kuipers5 examined the occurrence of heterogeneous flow structures in gas-particle flows and the impact of nonlinear drag force for both ideal particles (elastic collision, without interparticle friction) and nonideal particles (inelastic collision, with interparticle friction). Several drag models were compared,10,13,14,16 and the study showed that heterogeneous flow structures exist in systems with both nonideal and ideal particle-particle interaction. For the ideal particle interaction it was shown that heterogeneous flow structures are caused purely by the nonlinearity of the gas drag and that the stronger the dependence of drag to voidage, the more heterogeneous flow structures develop. The weak dependence of voidage to drag causes more homogeneous flow
1. Introduction Fluidized beds have been the center of modeling attention for many years, and various simulation approaches have been developed, ranging from Lattice-Boltzmann methods (LBM)1 to discrete particle (DPMs, Eulerian-Lagrangian)2-5 and two-fluid models (TFMs, Eulerian-Eulerian).6,7 More recently, discrete bubble models (DBMs)8 have been developed and been applied to fluidized beds, increasing the potential of industrial scale simulations beyond the capabilities of TFMs. The Eulerian formulation of the granular medium, using the kinetic theory of granular flows, has made the realization of fluidized bed simulations, less computationally intensive. The particulate phase is treated as a continuum with an effective viscosity, thus the method is also called a two-fluid approach. The calculation of the momentum exchange coefficient of gas-solid systems is very different in comparison of a single sphere surrounded by a fluid. When a single particle moves in a dispersed two-phase flow, the drag force is affected by the surrounding particles. Several correlations have been *To whom correspondence should be addressed. Telephone: 023 8059 8520. Fax: 023 8059 3230. E-mail:
[email protected]. (1) Ladd, A. J. C.; Veberg, R. J. Stat. Phys. 2001, 104, 1191. (2) Feng, Y. Q.; Xu, B. H.; Zhang, S. J.; Yu, A. B.; Zulli, P. AIChE. J. 2004, 50 (8), 1713–1728. (3) Hoomans, B. P. B.; Kuipers, J. A. M.; Mohd Salleh, M. A.; Stein, M.; Seville, J. P. K. Powder Technol. 2001, 116 (2-3), 166–177. (4) Kafui, K. D.; Thornton, C.; Adams, M. J. Chem. Eng. Sci. 2002, 57, 2395–2410. (5) Li, J.; Kuipers, J. A. M. Chem. Eng. Sci. 2003, 58, 711–718. (6) Gidaspow, D. Multiphase Flow and Fluidization: Continuum and Kinetic Theory Descriptions; Academic Press: New York, 1994. (7) Enwald, H.; Peirano, E.; Almstedt, A. E.; Leckner, B. Chem. Eng. Sci. 1999, 54, 311–328. (8) Bokkers, G. A.; Laverman, J. A.; Van Sint Annaland, M.; Kuipers, J. A. M. Chem. Eng. Sci. 2006, 61, 5590–5602. (9) Syamlal, M.; OBrien, T. J. AIChE. Symp. Ser. 1989, 85, 22–31. (10) Wen, C.-Y.; Yu, Y. H. Chem. Eng. Prog. Symp. Ser. 1966, 62, 100–111. r 2010 American Chemical Society
(11) McKeen, T.; Pugsley, T. Powder Technol. 2003, 129, 139–152. (12) Gibilaro, L. G.; Di Felice, R.; Waldram, S. P.; Foscolo, P. U. Chem. Eng. Sci. 1985, 40, 1817–1823. (13) Di Felice, R. Int. J. Multiphase Flow 1994, 20, 153–159. (14) Happel, J. AIChE. J. 1958, 4, 197. (15) Taghipour, F.; Ellis, N.; Wong, C. Chem. Eng. Sci. 2005, 60, 6857–6867. (16) Hill, R. The effects of fluid inertia on flow in porous media. Ph.D. Dissertation; Cornell University: 2001.
2133
pubs.acs.org/EF
Energy Fuels 2010, 24, 2133–2145
: DOI:10.1021/ef901497b
Papadikis et al.
structures. The simulations were also carried out in an Eulerian-Lagrangian framework. The work of Taghipour et al.15 concerned the Eulerian-Eulerian simulation of a fluidized bed comparing the effect of the drag correlations of Gidaspow,6 Syamlal-O’Brien,9 and Wen and Yu10 at the pressure drop and sand bed expansion. All studies showed good qualitative gas-solid flow patterns. Fluidized beds are the most widely used type of reactor for fast pyrolysis, as they offer a number of advantages, such as high heat transfer rates and good temperature control. Several models have been developed and applied to biomass pyrolysis.17-30 Most of the models apply numerical analysis to determine heat, mass, and momentum transport effects in a single biomass particle varying the pyrolysis conditions. Usually, the influence of parameters such as, size, shape, moisture, reaction mechanisms, heat transfer rates, and particle shrinkage is the main objective of study. The numerical investigations of di Blasi.17-22 are typical examples of single particle models that study the heat, mass, and momentum transport through biomass particles by varying the pyrolysis conditions, feedstock properties (implementing different chemical kinetics mechanisms), particle sizes, effect of shrinkage, etc. The models provide very useful information about different conditions of pyrolysis and product yields and can be an excellent guide for every researcher on the field. Saastamoinen,23 and Saastamoinen and Richard24 studied the effect of drying to devolatilization stating that the surface temperature of the particles after drying can exceed those of the initiation of devolatilization, meaning that drying and pyrolysis may overlap. The models of Babu and Chaurasia25-30 treat the pyrolysis process in a similar way as the one of di Blasi and focus more on the heat transfer effects and estimation of optimum parameters for biomass pyrolysis. They also study the effect of the enthalpy of reaction on product yields and vary the shrinkage parameters of biomass close to the complete degradation of the particles. They found that their approach gave better validation with experimental data compared to those of di Blasi. The main research of the authors has been focused on the incorporation of the single particle models presented in the
literature into a unified computational model that simulates the complete process of fast pyrolysis inside a bubbling fluidized bed.31-37 In this way, several important parameters, such as heat transfer coefficients from the bubbling bed, char and pyrolysis vapor residence times, biomass particle trajectories, and effect of particle size and shape, can be monitored and quantified, something that is impossible by single particle models alone. This unified CFD model can greatly aid the design and optimization of pyrolysis equipment by specifying the operation limits of the fluidized bed reactor regarding the type, the size, and the shape of biomass particles, as well as the high heat transfer regimes for more efficient biomass degradation and consequently efficient char elutriation at the freeboard. The scope of the current paper is to quantify the effect of the various drag models for granular flows that have been presented in the literature, to the fast pyrolysis of biomass. The different flow pattern and bubble distribution that will be created by these models will have a certain impact on the heat transfer coefficient, the momentum transport from the bed to the biomass particle, the final product yields, and more importantly the char and vapor residence times in the reactor. The reaction kinetics is based on a two-stage, semiglobal mechanism, with kinetic constants suitable for wood pyrolysis according to Chan et al.,38 for the primary stage and Liden et al.,39 and di Blasi17 for the secondary reactions. This scheme has been chosen for this study because it can predict the correct behavior of wood pyrolysis including the dependence of the product yields on temperature.18,22,38 The simulation results are compared with experimental data, whenever it is possible, that comes from the pyrolysis of spruce in the same reactor. 2. Model Description 2.1. Experiment. The biomass fast pyrolysis experiment was carried out in Aston University’s fluidized bed reactor. The reactor has a capacity of 100 g/h, is constructed of stainless steel, and is 40 mm in diameter and 260 mm high. The fluidizing gas was nitrogen, whereas sand was used as a fluidizing and heat transfer medium with an average particle diameter of 440 μm. The flow rate of nitrogen used was sufficient to provide 3 times the minimum fluidizing velocity (MFV) to the bed. The residence time of the pyrolysis vapors in the reactor was approximately 0.93s at a reaction temperature of 538 C. The main apparatus of the fast pyrolysis unit consists of the feeder, the reactor, and the liquid product collection system and is shown in Figure 1. 2.2. Numerical Model. The dimensions of the 100 g/h fast pyrolysis reactor are illustrated in Figure 2. Nitrogen flows
(17) Di Blasi, C. Combust. Sci. Technol. 1993a, 90, 315–339. (18) Di Blasi, C. Chem. Eng. Sci. 1996, 51, 1121–1132. (19) Di Blasi, C. Chem. Eng. Sci. 2000, 55, 5999–6013. (20) Di Blasi, C. Biomass Bioenerg. 1994, 7, 87–98. (21) Di Blasi, C. Ind. Eng. Chem. Res. 1996b, 35, 37–46. (22) Di Blasi, C. J. Anal. Appl. Pyrol. 1998, 47, 43–64. (23) Saastamoinen, J. Model for drying and pyrolysis in an updraft gasifier. In Advances in Thermochemical Biomass Conversion; Bridgwater, A. V. Ed.; Blackie: London, 1993; pp 186-200. (24) Saastamoinen, J.; Richard, J.-R. Combust. Flame 1996, 106, 288– 300. (25) Babu, B. V.; Chaurasia, A. S. Modelling and Simulation of Pyrolysis: Influence of Particle Size and Temperature. Proceedings of the International Conference on Multimedia and Design, 2002a, 4, Mumbai, India; pp 103-128. (26) Babu, B. V.; Chaurasia, A. S. Modelling and Simulation of Pyrolysis: Effect of Convective Heat Transfer and Orders of Reactions. Proceedings of International Symposium and 55th Annual Session of IIChE (CHEMCON-2002), 2002b, OU, Hyderabad, India; pp 105-106. (27) Babu, B. V.; Chaurasia, A. S. Energ. Convers. Manage. 2003a, 44, 2251–2275. (28) Babu, B. V.; Chaurasia, A. S. Energ. Convers. Manage. 2003b, 44, 2135–2158. (29) Babu, B. V.; Chaurasia, A. S. Modelling and Simulation of Pyrolysis of Biomass: Effect of Heat of Reaction. Proceedings of International Symposium on Process Systems Engineering and Control (ISPSEC 03) for Productivity Enhancement Through Design and Optimisation 2003c, IIT-Bombay, Mumbai, India; pp 181-186. (30) Babu, B. V.; Chaurasia, A. S. Energ. Convers. Manage. 2004c, 45, 1297–1327.
(31) Papadikis, K.; Bridgwater, A. V.; Gu, S. Chem. Eng. Sci. 2008, 63 (16), 4218-4227. (32) Papadikis, K.; Gu, S.; Bridgwater, A. V. Chem. Eng. Sci. 2009, 64 (5), 1036-1045. (33) Papadikis, K.; Gerhauser, H.; Bridgwater, A. V.; Gu, S. Biomass Bioenerg. 2009, 33 (1), 97-107. (34) Papadikis, K.; Gu, S.; Bridgwater, A. V.; Gerhauser, H. Fuel Process. Technol. 2009, 90 (4), 504-512. (35) Papadikis, K.; Gu, S.; Bridgwater, A. V. Chem. Eng. J. 2009, 149 (1-3), 417–427. (36) Papadikis, K.; Gu, S.; Bridgwater, A. V. Biomass Bioenerg. 2010, 34 (1), 21–29. (37) Papadikis, K.; Gu, S.; Bridgwater, A. V. Fuel Process. Technol. 2010, 91 (1), 68–79. (38) Chan, W. R.; Kelbon, M.; Krieger, B. B. Fuel 1985, 64, 1505– 1513. (39) Liden, A. G.; Berruti, F.; Scott, D. S. Chem. Eng. Commun. 1988, 65, 207–221.
2134
Energy Fuels 2010, 24, 2133–2145
: DOI:10.1021/ef901497b
Papadikis et al.
size of the particles used in the experiments. Bigger rigs and commercial plants use much larger particles in the range of 2-5 mm. The different bubble distribution that will occur due to the different granular flow drag models will have an impact on the heat transfer conditions and consequently will affect the degradation rate of the particles and subsequently the final product yields and the residence time of char and vapors. The numerical model will try to identify any significant differences between the models and compare the results with reallife experiments. The simulations last until the particles are entrained from the reactor independently of the simulation time needed to achieve that. 3. Mathematical Model 3.1. Multiphase Flow Governing Equations. The simulations of the bubbling behavior of the fluidized bed were performed by solving the equations of motion of a multifluid system. An Eulerian model for the mass and momentum for the gas (nitrogen) and fluid phases, was applied, while the kinetic theory of granular flow was applied for the conservation of the solid’s fluctuation energy. The Eulerian model is already incorporated in the main code of FLUENT, and its governing equations are expressed in the following form. Mass Conservation. Eulerian-Eulerian continuum modeling is the most commonly used approach for fluidized bed simulations. The accumulation of mass in each phase is balanced by the convective mass fluxes. The phases are able to interpenetrate, and the sum of all volume fractions in each computational cell is unity. Gas Phase.
Figure 1. Experimental apparatus. Adapted from Coulson M, EC Project Bioenergy Chains Final Report, Aston University, ENK6-2001-00524.
∂ðεg Fg Þ þ r 3 ðεg Fg υg Þ ¼ 0 ∂t
ð1Þ
∂ðεs Fs Þ þ r 3 ðεs Fs υs Þ ¼ 0 ∂t
ð2Þ
Solid Phase.
Momentum Conservation. Newton’s second law of motion states that the change in momentum equals the sum of forces on the domain. In gas-solid fluidized beds the sum of forces consists of the viscous force r 3 τCs , the solids pressure force rps, the body force εsFsg, the static pressure force εs 3 rp and the interphase force Kgs(ug - us) for the coupling of gas and solid momentum equations by drag forces. Gas Phase.
Figure 2. Fluidized bed reactor geometry used in the simulation.
through a porous plate at the bottom of the reactor at a velocity of U0 = 0.6 m/s and temperature of 303 K. The superficial velocity is approximately 3 times greater than the minimum fluidizing velocity Umf of the reactor, which is typically around 0.2 m/s using a sand bed with average particle diameter of 440 μm.40 The sand bed has been preheated to 815 K as it was also performed in the experimental procedure. The biomass particle of density of 700 kg/m3 is injected at the center of the sand bed, which has been previously fluidized for 1 s. Momentum is transferred from the bubbling bed to the biomass particle as well as from the formed bubbles inside the bed. The studied biomass particle was chosen to be 350 μm in diameter, which is more-or-less the
∂ðεg Fg υg Þ þ r 3 ðεg Fg υg Xυg Þ ∂t ¼
¼ -εg 3 rp þ r 3 τg þ εg Fg g þ Kgs ðug -us Þ
ð3Þ
Solid Phase. ∂ðεs Fs υs Þ ¼ þ r 3 ðεs Fs vs Xυs Þ ¼ -εs 3 rp -rps þ r 3 τs ∂t þ εs Fs g þ Kgs ðug -us Þ ð4Þ where the solid-phase stress tensor is given by, ¼
τs ¼
(40) Geldart, D. Powder Technol. 1973, 7, 285.
2135
εs μs ðrus þ ruTs Þ þ εs
¼ 2 λs - μs r 3 us I s 3
ð5Þ
Energy Fuels 2010, 24, 2133–2145
: DOI:10.1021/ef901497b
Papadikis et al.
motion. This phenomenon is quantified by the term called frictional viscosity
The Gidaspow interphase exchange coefficient is, 3 εs εg Fg jus -ug j -2:65 Kgs ¼ Cd εg for εg > 0:8 ds 4 Kgs ¼ 150
ε2s μg εs Fg jus -ug j þ 1:75 for εg e 0:8 2 ds εg ds
ð6Þ
μs, fr ¼
24 ½1 þ 0:15ðεg Res Þ0:687 εg Res
ð8Þ
and ds Fg jus -ug j Res ¼ μg
The Syamlal-O’Brien kinetic viscosity is pffiffiffiffiffiffiffiffiffi εs Fs ds Θs π 2 μs, kin ¼ 1 þ εs g0, ss ð1 þ ess Þð3ess -1Þ 6ð3 -ess Þ 5
ð9Þ
The Syamlal-O’Brien interphase exchange coefficient is, ! 3εs εg Fg Res jus -ug j Cd ð10Þ Kgs ¼ 2 4 ur, s ds ur, s
ð22Þ The solids pressure ps, which represents the normal force due to particle interactions is given by ps ¼ εs Fs Θs þ 2Fs ð1 þ ess Þε2s g0, ss Θs
where
!
Cd ¼
4:8 0:63 þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Res =ur, s
and ur, s ¼ 0:5ðA -0:06Res qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ ð0:06Res Þ2 þ 0:12Res ð2B -AÞ þ A2 Þ
ð12Þ
with ð13Þ
A ¼ εg4:14 , B ¼ ε2:65 g , for εg > 0:85
ð14Þ
¼
where τCs is defined in eq 5. The diffusion coefficient of granular temperature kΘs according to eq 6 is given by:
The Wen-Yu interphase exchange coefficient is,
kΘs ¼ ð15Þ
where Cd ¼
24 ½1 þ 0:15ðεg Res Þ0:687 εg Res
ð25Þ
ð26Þ
the collision dissipation energy as γΘs ¼
The tangential forces due to particle interactions are summarized in the term called solids shear viscosity, and it is defined as
where the collision viscosity of the solids μs,col is rffiffiffiffiffiffi 4 Θs μs, col ¼ εs Fs ds g0, ss ð1 þ ess Þ 5 π
pffiffiffiffiffiffiffiffiffi 2 150Fs ds Θs π 6 1þ εs g0, ss ð1þ ess Þ 384ð1 þ ess Þg0, ss 5 rffiffiffiffiffiffi Θs 2 þ 2Fs ds εs g0, ss ð1 þ ess Þ π
The radial distribution function g0,ss is defined as 2 !1=3 3 -1 ε s 5 g0, ss ¼ 41 εs, max
ð16Þ
The bulk viscosity λs is a measure of the resistance of a fluid to compression, which is described with the help of the kinetic theory of granular flows rffiffiffiffiffiffi 4 Θs ð17Þ λs ¼ εs Fs ds g0, ss ð1 þ ess Þ 3 π
μs ¼ μs, col þ μs, kin þ μs, fr
¼
ð -ps I s þ τs Þ : r 3 us þ r 3 ðkΘs 3 r 3 Θs Þ -γΘs þ φgs ð24Þ
or
εs εg Fg jus -ug j -2:65 3 Kgs ¼ Cd εg ds 4
ð23Þ
Fluctuation Energy Conservation of Solid Particles. The solid phase models discussed above are based on two crucial properties, namely the radial distribution function g0,ss and granular temperature Θs. The radial distribution function is a measure for the probability of interparticle contact. The granular temperature represents the energy associated with the fluctuating velocity of particles. 3 ∂ ðεs Fs Θs Þ þ r 3 ðεs Fs us Θs Þ ¼ 2 ∂t
ð11Þ
A ¼ εg4:14 , B ¼ 0:8εg1:28 , for εg e 0:85
ð20Þ
which contains the angle of internal friction φ. This expression of frictional viscosity is valid if the volume fraction of solids is constant, which is reasonable in the dense packed bed state of particles. The Gidaspow6 kinetic viscosity is pffiffiffiffiffiffiffiffiffi 2 10Fs ds Θs π 4 1þ εs g 0, ss ð1þess Þ ð21Þ μs, kin ¼ 96εs g0, ss ð1 þ ess Þ 5
ð7Þ
where the drag coefficient is given by Cd ¼
ps sinðφÞ pffiffiffiffiffiffiffi 2 I2D
12ð1 -e2ss Þg0, ss pffiffiffi Fs ε2s Θ3=2 s ds π
ð27Þ
and the transfer of kinetic energy φgs as
ð18Þ
jgs ¼ -3Kgs Θs
ð28Þ
ð19Þ (41) Ranz, W. E.; Marshall, W. R. Chem. Eng. Prog. 1952a, 48, 141– 146. (42) Ranz, W. E.; Marshall, W. R. Chem. Eng. Prog. 1952a, 48, 173– 180.
At the packed state of the bed, the stresses are dominated by interparticle friction rather than collisions and fluctuating 2136
Energy Fuels 2010, 24, 2133–2145
: DOI:10.1021/ef901497b
Papadikis et al.
An analytical discussion of the solid-phase properties can be found in ref 43. 3 2. Forces on Discrete Particles. The coupling between the continuous and discrete phases has been developed in a UDF to take into account the bubbling behavior of the bed. For an analytical discussion of this section the reader is referred to the previous work done by the authors in this aspect.31 Assuming a spherical droplet with material density of Fd inside a fluid, the rate of change of its velocity can be expressed as44 ! dud f Fcon þ Fυm ¼ ðucon -ud Þ þ g 1 ð29Þ τu dt Fd
Figure 3. Two-stage, semiglobal model.49
and the drag force between solid and liquid is computed for a modified particle volume fraction εp εd εp ¼ ð37Þ εs þ εd and an effective continuum viscosity μeff,con εp -1:55 μeff , con ¼ 1 εdm
where f is the drag factor and τu is the velocity response time Fd dd 18μcon 2
τu ¼
ð30Þ
If the volume fraction of the space among the solid particles, if they were closely packed, is larger than the liquid fraction
There are several correlations for the drag factor f in the literature.45-47 The one used in this study is the correlation of Putnam47
εs > εs
ð31Þ
f ¼ 0:0183Rer for 1000 e Rer < 3 105
ð32Þ
εdg ¼ εd ð1 -εs =εs Þ
3.3. Reaction Kinetics. The reaction kinetics of biomass pyrolysis is modeled using a two-stage, semiglobal model.49 The mechanism is illustrated in Figure 3. The mechanism utilizes the Arrhenius equation which is defined as ð42Þ Ki ¼ Ai expð -Ei =RTÞ
According to Kolev,48 if bubble three-phase flow (i.e., solid particles in bubbly flow) is defined, two subcases are distinguished. If the volume fraction of the space among the solid particles, if they were closely packed, is smaller than the liquid fraction (in this case the Eulerian sand)
The values of the kinetic parameters were obtained by Chan et al.38 for the primary pyrolysis products, while the fourth and fifth reactions were obtained from Liden et al.39 and Di Blasi,17 respectively. The choice of a different set of parameters used in the simulation was mainly done to predict as correctly as possible the fast pyrolysis of wood. The kinetic parameters of the three-step model of Chan et al.38 were used for the primary products because, when coupled with secondary reaction kinetics, they can predict, at least qualitatively, the correct behavior of wood pyrolysis.22 However, the secondary reaction kinetic constants were taken from the work of Liden et al.39 and di Blasi17 because these studies included secondary reaction effects that depend on the concentration of tar vapors from flash pyrolysis in the porous matrix of the solid, something that is more suitable for the current study. The elemental analysis of the pine wood used in the experiments of Chan et al. can be found in refs 38 and 50. The model associates, via a UDF, the reaction kinetics mechanism with the discrete biomass particle injected in the fluidized bed. The particle’s properties change according to the reaction mechanism due to the phase transition phenomena. Momentum and heat transfer on the particle are
ð34Þ
where
εs ¼
1 -εdm εd εdm
ð40Þ
are surrounded by gas and the drag force can be calculated between one single solid particle and gas as for a mixture εdg εp ¼ ð41Þ εg þ εdg
The second term on the right-hand side of the equation represents the gravity and boyancy force, while the third term represents the unsteady force of virtual mass, which is expressed as F Vd ducon dud ð33Þ Fυm ¼ con 2 dt dt
ð39Þ
then only
Rer ð2=3Þ for Rer < 1000 f ¼ 1þ 6
εs < εs
ð38Þ
ð35Þ
then the theoretical possibility exists that the particles are carried only by the liquid. The hypothesis is supported if we consider the ratio of the free setting velocity in gas and liquid sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Fd -Fg Fs wdg ¼ .1 ð36Þ wds Fd -Fs Fg Due to great differences between gas and liquid densities, the particles sink much faster in gas than in a liquid. Therefore, the drag force between gas and solid particles is zero, (43) Boemer, A.; Qi, H.; Renz, U. Int. J. Multiphas. Flow 1997, 23 (5), 927–944. (44) Crowe, C. T.; Sommerfeld, M.; Tsuji, Y. Multiphase Flows with Droplets and Particles; CRC Press LLC: 1998. (45) Schiller, L.; Naumann, A. Ver. Deut. Ing. 1933, 77, 318. (46) Clift, R.; Gauvin, W. H. Proc. Chemeca 1970, 1, 14. (47) Putnam, A. ARS JNl. 1961, 31, 1467. (48) Kolev, N. I. Multiphase Flow Dynamics 2. Thermal and Mechanical Interactions, 2nd ed.; Springer: 2005.
(49) Shafizadeh, F.; Chin, P. P. S. ACS Symp. Ser. 1977, 43, 57-81. (50) Chan, W. R.; Kelbon, M.; Krieger, B. B. Ind. Eng. Chem. Res. 1988, 27, 2261–2275.
2137
Energy Fuels 2010, 24, 2133–2145
: DOI:10.1021/ef901497b
Papadikis et al. Table 1. Simulation Parameters
property
value
biomass density, Fw biomass particle diameter, dp biomass specific heat capacity, Cpw char specific heat capacity, Cpc biomass thermal conductivity, kw char thermal conductivity, kc superficial velocity, U0 gas density, Fg gas viscosity, μg gas specific heat capacity, Cp,g gas thermal conductivity, kg solids particle density, Fs sand specific heat capacity, Cp,s sand thermal conductivity, ks mean solids particle diameter, ds restitution coefficient, ess initial solids packing, εs static bed height bed width heat of reaction
700 kg/m 350 μm 1500 J/(kg K) 1100 J/kg K 0.105 W/(m K) 0.071 W/(m K) 0.6 m/s 1.138 kg/m3 1.663 10-5 kg/(m s) 1040.67 J/(kg K) 0.0242 W/(m K) 2500 kg/m3 835 J/(kg K) 0.35 W/(m K) 440 μm 0.9 0.6 0.08 m 0.04 m ΔH = -255 kJ/kg
Nu ¼ 2 þ 0:9Rep, mf 0:62 ðdp =db Þ0:2
Rep, mf ¼
ð45Þ
The effective thermal conductivity keff and effective specific heat capacity Cp,eff are given by
Cp eff ¼ Cp c þ jCp w -Cp c jψw
ð47Þ
The heat transfer coefficient is evaluated from the wellknown Ranz-Marshall41,42 correlation, when the particle is carried only by the fluidizing gas hdp ¼ 2:0 þ 0:6Rep 1=2 Pr1=3 kg
ð50Þ
For the implementation of the model, certain parameters have been quantified and assumptions have been made in order to provide, as much as possible, an insight to the fast pyrolysis process in bubbling beds. (1) The particles used in the simulation were assumed to be totally spherical, whereas the particles used in experiments can be found on all sorts of shapes. The actual sphericity of the particles greatly differs from 1. This would have an impact on the drag and virtual mass forces and consequently on the trajectory of the particle inside the reactor. (2) The model assumes a plug flow profile at the inlet of the reactor. (3) The geometry of the reactor has been discretized using a structured grid. The average side length of the computational cells is about 1 mm, resulting in a total number of 10 440 cells. (4) The time-step used for the simulation was on the order of 10-5 s. The reason for this small time advancements was that the computational algorithms in the UDF had to converge simultaneously with the algorithms of FLUENT. The explicit central difference scheme used for the discretization of the heat diffusion equation demands very small time-steps due to the small radial discretization resulting from particles of 350 μm in diameter. For further analysis of the CFD code structure together with the flow-chart and advantages and drawbacks of the model, the reader is referred to the previously published work of the authors in this aspect.32 The simulation parameters are shown in Table 1.
r ¼0
ð46Þ
Fg Umf dp μg
4. Model Parameters
r ¼R
keff ¼ kc þ jkw -kc jψw
ð49Þ
where
The boundary condition at the surface of the particle is given by ∂T hðT¥ -Ts Þ ¼ -keff ð44Þ ∂r and at the center of the particle ∂T ¼0 ∂r
wood fixed wood char wood char ≈3Umf nitrogen (303 K) nitrogen (303 K) nitrogen (303 K) nitrogen (303 K) sand fixed fixed uniform distribution value in literature fixed value fixed value fixed value Koufopanos et al.52
coefficient according to the findings of Collier et al.,51
calculated according to the new condition in the UDF, and the variables associated with it are updated in each time-step. Intraparticle secondary reactions due to the catalytic effect of char are taken into account, resulting in secondary vapor cracking. 3.4. Heat Transfer. The heat conduction along the radius of the particle is calculated by solving the heat diffusion equation for an isotropic particle. The particle is discretized in the UDF for the complete pyrolysis model and the reader is referred to the work previously done by the authors in this aspect.32 ∂ 1 ∂ ∂F 2 ∂T ð FCp eff TÞ ¼ 2 keff r þ ð -ΔHÞ ð43Þ ∂t r ∂r ∂r ∂t
Nu ¼
comment
3
ð48Þ
However, since the diameter of the biomass particle is smaller than the diameter of the sand particles, a modified Nusselt number is used to calculate the heat transfer
5. Results 5.1. Fluidized Bed Hydrodynamics and Particle Trajectories. Figures 4-6 illustrate the fluidized bed hydrodynamics
(51) Collier, A. P.; Hayhurst, A. N.; Richardson, J. L.; Scott, S. A. Chem. Eng. Sci. 2009, 59, 4613–4620.
(52) Koufopanos, C. A.; Papayannakos, N.; Maschio, G.; Lucchesi, A. Can. J. Chem. Eng. 1991, 69, 907–915.
2138
Energy Fuels 2010, 24, 2133–2145
: DOI:10.1021/ef901497b
Papadikis et al.
Figure 4. Volume fraction of sand with instantaneous position of the biomass particle, Gidaspow drag model.
Figure 6. Volume fraction of sand with instantaneous position of the biomass particle, Wen-Yu drag model.
mainly caused by the tendency of the particles to move toward the wall of the reactor as it is shown in Figure 7 where gas velocities tend to be much lower than at the central axis of the reactor. This is valid for all the drag models compared in this study, however the bubble distribution created from the Gidaspow drag model moved the particle to the right-hand side of the reactor in contrast to the Syamlal-O’Brien and Wen-Yu drag models, which moved the particle on the left-hand side. The 2D nature of the simulation excludes the mixing of the biomass particles in the z-direction of the reactor, however the same behavior was also observed in the previous papers by the authors,35,37 where 3D geometry was considered. The inclusion of the pyrolysis vapors mass source in the simulation greatly discourage the use of a 3D geometry since the simulation times approach the order of months even with the use of 16 processors. Another phenomenon that is of great importance is the transverse mixing of the particles. Despite the dominance of the motion in the y-direction of the reactor, we can easily notice that the particles have moved across the whole width of the bed (x-direction) before their entrainment. This comes in great agreement with the experimental study of Hoomans et al.,3 where they showed the validity of this phenomenon using the positron emission particle tracking (PEPT) technique, and the numerical study of Xu and Yu.53 Despite the fact that these studies were performed using 2D geometries, the same phenomenon was also observed in the previous numerical investigations performed by the authors31,35,37 where a 3D geometry was used. 5.2. Heat Transfer. Figure 8 shows the temperature rise of the particles for the different drag models. The temperature rise in the Gidaspow and Syamlal-O’Brien model appears to be relatively similar, however a higher temperature rate is observed for the Wen-Yu drag model, represented by a steeper surface temperature profile at the first 0.1 s after the
Figure 5. Volume fraction of sand with instantaneous position of the biomass particle, Syamlal-O’Brien drag model.
together with the instantaneous position of the biomass particle indicated by the black spot. The snapshots were taken for 1 s time intervals, from the time of injection at 1 s until close to the entrainment time of approximately 4 s. Figure 7 illustrates the complete trajectory of the particles for the three different drag models. The residence time of the particles was finally determined as 2.92, 3.3, and 2.6 s for the Gidaspow, Syamlal-O’Brien, and Wen-Yu drag models, respectively. One can observe that the residence time of the particles is somewhat higher, though not unreasonable, compared to what the fast pyrolysis process suggests, i.e., < 2 s. This is
(53) Xu, B. H.; Yu, A. B. Chem. Eng. Sci. 1997, 52, 2785–2809.
2139
Energy Fuels 2010, 24, 2133–2145
: DOI:10.1021/ef901497b
Papadikis et al.
Figure 7. Particle trajectories inside the fluidized bed reactor.
transfer between a fluidized bed of particles larger than the injected particles. In our case the bed consists of particles of 440 μm in diameter, whereas the biomass particles are 350 μm in diameter. Collier et al.51 showed that when a particle finds itself among fluidized particles of larger diameter, the heat transfer coefficient appears to have a constant value depending on the ratio of the diameter of the particles and the minimum fluidization velocity of the bed. This is because the conductive effect of the larger particles to the smaller ones is negligible and convective effects are dominant. This approach is used in the current study. Therefore, whenever the heat transfer coefficient appears to have this constant value, it means that the particle finds itself in a densely packed zone of sand particles. The other values of the heat transfer coefficient represent the motion of the particle inside a bubble, or the splash zone, or the freeboard of the reactor. This gives a good indication of the region that the particle is instantly located. The high heat transfer values, close to 700-750 W/(m2 K), that are observed in all of the three diagrams clearly indicate the motion of the particle close to the wall of the reactor, where the wall heat transfer effects become more noticeable. In the current simulation only the convective effects from the wall were taken into account and the collision of the particles with the wall were neglected, since the actual contact time of the particles with the wall is extremely small. However, Figure 9 in conjunction with Figure 7 clearly show the wall effect on the rise of the heat transfer coefficient. This becomes even more noticeable when the particle flies toward the outlet of the reactor and the heat transfer coefficient drops significantly in all cases, only to rise again in the actual outlet of the reactor where the x-velocity of the gases increase significantly due to the small diameter of the outlet.
Figure 8. Surface and center temperatures of the particles.
particle injection. In the case of the Wen-Yu drag model, the particle reaches the splash zone of the reactor much faster than in the Gidaspow and Syamlal-O’Brien drag models. In this region, high velocity gradients are dominant due to the eruption of bubbles, and the slip velocity between the fluidizing gas and the particles is high, resulting in a very high convective heat transfer coefficient. This can be easily seen in Figure 9, where the heat transfer coefficient becomes very large for the case of the Wen-Yu model almost instantly. In the other two cases a certain amount of time (≈0.3 s) is needed for the particles to experience the same high heat transfer rates. The constant value of the heat transfer coefficients, which is close to 320 W/(m2 K), comes from the correlation of Collier et al. in eq 49, which is appropriate for the heat 2140
Energy Fuels 2010, 24, 2133–2145
: DOI:10.1021/ef901497b
Papadikis et al.
Figure 9. Heat transfer coefficient.
appears to have more intense y-velocity fluctuations, and this is because, as it was discussed earlier, it reaches the splash zone of the reactor almost 0.2 s faster than the other particles. This high y-velocity fluctuation gave rise to the very high heat transfer coefficients that were observed earlier and consequently the more rapid density drop of the particle. Comparing Figures 8, 9, and 11. one can conclude that the unstable splash zone of the reactor is the most desirable region for fast pyrolysis. In other words, the violent velocity gradients that occur due to bubble eruption and sand recirculation have a higher impact on the particles that float on the top of the bed, in contrast with the particles that are immersed deeply inside it. Also, the high velocity gradients give rise to the virtual mass force that can greatly contribute to the violent ejection of the particles from the sand bed to the freeboard of the reactor, thus making the elutriation phenomenon easier, when the density drop of the particles has reached the desirable level. It has to be noted that for the current simulation the particles were considered to be totally spherical, something that differs from reality. Also, shrinking and attrition effects were assumed not to play an important role and therefore were neglected in the simulation. For the complete analysis of the momentum transport model, the reader is referred to the previous study of the authors,31 and for the modeling of the impact of shrinkage to ref 35. 5.4. Product Yields. Figure 12 shows the final product yields (vapors, gases, char) percentage by weight of original wood. We can observe that the different drag models actually have an impact on the final product yield of condensable volatiles (vapors), char, and gases, while there is still a small percentage of wood that has not yet reacted. The final product yields for the three cases are given in Table 2, and the experimental data is given in Table 3. In general terms, we can observe that there is a similarity between the final vapors and the gases in the experiment and the simulation. However there is a significant difference on
Figure 10. Density drop of the particles.
The effect of the heating rates can be seen in Figure 10, where the density drop due to the reaction mechanism is illustrated. The Gidaspow and the Syamlal-O’Brien models appear to have the same effect on the density drop of the particle, whereas the Wen-Yu model has a more rapid effect due to the higher heat transfer rates that appear during the simulation. From Figure 10 one can also observe the onset of pyrolysis, which comes close to 0.2 s, where, according to Figure 8, the temperature is close to 750 K. The density drop of the particles has a major effect on their trajectory inside the reactor since their motion is heavily dependent upon their density and consequently their mass that is involved on the particle motion equations given in Section 3.2. 5.3. Particle Dynamics. Figure 11 shows the instantaneous velocity components of the particles inside the reactor until their entrainment time. All the drag models appear to have a similar effect on the average velocity magnitude of the particles with some exceptions regarding their instantaneous positions. For example, the particle on the Wen-Yu case 2141
Energy Fuels 2010, 24, 2133–2145
: DOI:10.1021/ef901497b
Papadikis et al.
Figure 11. Particle velocities inside the fluidized bed reactor.
Figure 12. Product yields. Table 2. Product Yields from the Simulation
Gidaspow Syamlal-O’Brien Wen-Yu
Table 3. Experimental data and product yields
wood
vapors
char
gases
pyrolysis temperature C
3.5% 2.0% 6.4%
62% 63% 60%
20.5% 20.8% 20.0%
14.0% 14.2% 13.6%
≈525 ≈525 ≈528
the final char content, and differences of up to 7.5% can be observed. This is mainly due to the reaction kinetics mechanism and not to the fluidized bed operation. Due to the lack of analytic kinetic data regarding spruce fast pyrolysis, the selection of kinetic data that would as closely as possible 2142
feedstock moisture content (wt %) ash content (wt %, dry basis) pyrolysis temperature vapor residence time
spruce 6.96 0.32 538 C
