A granular-biomass high temperature pyrolysis model based on the Darcy flow

Abstract

We established a model for the chemical reaction kinetics of biomass pyrolysis via the high-temperature thermal cracking of liquid products. We divided the condensable volatiles into two groups, based on the characteristics of the liquid prdoducts., tar and biomass oil. The effects of temperature, residence time, particle size, velocity, pressure, and other parameters on biomass pyrolysis and high-temperature tar cracking were investigated numerically, and the results were compared with experimental data. The simulation results showed a large endothermic pyrolysis reaction effect on temperature and the reaction process. The pyrolysis reaction zone had a constant temperature period in several layers near the center of large biomass particles. A purely physical heating process was observed before and after this period, according to the temperature index curve.

Keywords

biomass pyrolysis / high temperature pyrolysis model / condensable volatile cracking / Darcy flow

Introduction

Energy shortages and emissions from energy application processesimpose severe limits on the development of human society and threaten the survival of humanity ( Chen et al., 2007; Yang and Chen, 2011). These problems have drawn the long-term attention of the international community. China is one of the largest energy consumers and Green House Gas (GHG) emitters in the world; and the government is taking multiple measures to solve the problem ( Chen and Chen, 2006a, b). Biomass is a clean and renewable energy, especially in rural areas ( Yang et al., 2011; Chen and Chen, 2012a, b; Zhang et al., 2012), and pyrolysis is an important way to use biomass. Anaerobic digestion and biogas production are promising ways to achieve energy and environmental benefits at both local and global levels ( Ji et al., 2009; Chen et al., 2010, 2012c). Pyrolysis is a thermo chemical transformation process in anoxic or anaerobic conditions in which the thermal degradation of solid particles results in a large number of chemical products, consisting of non-condensable gas (light gas), condensable volatiles (biomass oil and tar), and charcoal ( Jiang et al., 2007; Di Blasi, 2008; Ju and Chen, 2011). The reactive products of biomass pyrolysis originate from the secondary cracking of the initial condensable volatile pyrolysis products, during which the condensable organic products are decomposed further into low molecular weight gases and charcoal ( Chen and Chen, 2011). Such procedures are of particular importance in coping with the increasing pressure of problems related to global energy scarcity and climate change ( Qu et al., 2009a, b).

A large number of factors influence the output, reaction rate, product components, and response characteristics of biomass pyrolysis. Temperature, pressure, and the pyrolysis rate are the primary operating parameters, but biomass parameters (such as chemical components and ash content, and particle size, shape, and density) are also critical factors ( Chen et al., 2006b; Dai et al., 2012). The composition of the non-condensable gases includes carbon dioxide, carbon monoxide, methane, and low molecular weight hydrocarbons ( Piskorz et al., 1998; Su et al., 2012). The composition of the condensable volatiles only depends on the pyrolysis temperature, residence time, and charcoal layer ( Babu and Chaurasia, 2004; Hubacek et al., 2012). Therefore, it is essential to understand the influence of the temperature, transfer processes, and chemical reactions on the output and composition of biomass pyrolysis products ( Feng et al., 2003; Yang et al., 2009). At high temperatures, heat is delivered to the biomass particle surface by convection and radiation, and heat is transferred from the surface to the interior of the particle by conduction ( Di Blasi, 1996; Chen and Chen, 2012c). The interior temperature of the particle gradually rises as heat is delivered from the surface, resulting in a chemical reaction with volatile products and a solid carbon residue ( Di Blasi, 2000; Chen et al., 2012a, b). The heat and volatile chemical fluxes influence heat transfer. The biomass pyrolysis process can be divided into three regions: the first region is the outer part of the solid particle, where the pyrolysis reaction is complete; the second is the interlayer, where pyrolysis is occurring; and the third is where pyrolysis has not yet occurred ( Jiang et al., 2008, 2009; Chen et al., 2012c; Liu et al., 2011b). Condensable volatiles produced by secondary thermal cracking change the composition of the pyrolysis products. Diffusion effects caused by the temperature difference between the exterior and interior of the particle influence the pyrolysis product distribution ( Zeng et al., 2010; Yang et al., 2012). In fact, the temperature is also affected by heat absorption and release during the pyrolysis reaction itself ( Jiang and Chen, 2011; Liu et al., 2012). It is, therefore, difficult to model pyrolysis because of the complexity of the processes, and because information on the thermochemical equilibrium, heat and mass transfer, and the chemical reaction kinetics model is required.

Chemical reaction kinetics model

Diffusion effects caused by the temperature difference between the exterior and interior of the particle influence the pyrolysis product distribution to a great extent ( Di Blasi, 1993a, b; Zhang et al., 2011). The pyrolysis temperature is also affected by heat absorption and release during the pyrolysis process ( Koufopanos et al., 1989, 1991). It is difficult for large particle models to simulate the pyrolysis process entirely because of its complexity ( Kansa et al., 1977; Lu et al., 2012; Song et al., 2012). However, in industrial practice, reducing the size of biomass particles is uneconomic. It is, therefore, increasingly important for pyrolysis processes to be simulated for large particles.

Particle characteristic changes and physical and chemical changes during the pyrolysis process must be considered when modeling large particle pyrolysis ( Grønli and Melaaen, 2000; Chen et al., 2006a; Wei et al., 2009). Test measurements indicate that the production of non-condensable gases and condensable volatiles by pyrolysis is an endothermic process. The secondary cracking process is an exothermic process. Many physical parameters, including the conductivity factor, specific heat capacity, density, porosity, blackness, average molecular weight of products, viscosity, and diffusion coefficient, that are required for a pyrolysis model for large particles can be found in the literature ( Piskorz et al., 1998; Liu et al., 2011a, b). Since different assumptions are used in different pyrolysis models, and chemical reaction dynamics are simplified, a suitable model for predicting biomass pyrolysis and high-temperature pyrolysis over a wide temperature range has not yet been developed ( Melaaen, 1996; Chen and Chen, 2009).

Non-condensable gases, condensable volatiles, and solids are generated by pyrolysis. It is essential to control the production of condensable volatiles in a successful biomass gasification system. A tar measurement conference was held in Brussels in 1998, during which a number of experts agreed that tar should be defined as organic chemicals with higher molecular weights than benzene ( Maniatis and Beenackers, 2000; Su et al., 2009, 2011). For the purposes of this study, biomass oil was defined as pyrolysis-generated hydrocarbon or oxy-hydrocarbon compounds with lower molecular weights than benzene , the main components of which are acids, furans, alcohols, phenols, and paraffins. These low molecular weight hydrocarbons can be used as fuel in gas turbines and engines ( Li and Suzuki, 2009).These compounds are mainly generated at temperatures around 500°C. During this study, liquid generated by biomass pyrolysis was divided into the categories of tar and biomass oil according to these definitions ( Hagge and Bryden, 2002). The biomass pyrolysis process for making oil can be simulated using this means of classification. At high temperatures, the relationship between tar content and pyrolysis temperature can also be found, the details of which are shown in Fig. 1.

There are two steps in achieving biomass pyrolysis. The first step is the initial pyrolysis reaction (k₁, k₂), which describes how the biomass is transformed into gas, coke, and biomass oil. The second reaction, cracking, is assumed a first order reaction with Arrhenius temperature characteristics. This process can lead to the biomass’s being aggregated into tar. The biomass can also be cracked into non-condensable gas at high temperatures, and tar can be further cracked into non-condensable gas and coke.

Based on the pyrolysis reaction kinetics scheme (k₁–k₅) shown in Fig. 1, the relationship between the generation and consumption of each component can be expressed as follows:

(1) {\displaystyle \frac{\partial {\rho }_B}{\partial t}}=-k_{\mathrm{1}}{\rho }_B-k_{\mathrm{2}}{\rho }_B,

(2) {\displaystyle \frac{\partial {\rho }_C}{\partial t}}=k_{\mathrm{1}}{\rho }_B+ϵk_{\mathrm{5}}{\rho }_T,

(3) {\displaystyle \frac{\partial ϵ{\rho }_{BO}}{\partial t}}=k_{\mathrm{2}}{\rho }_B-ϵk_{\mathrm{3}}{\rho }_{BO}-ϵk_{\mathrm{4}}{\rho }_{BO},

(4) {\displaystyle \frac{\partial ϵ{\rho }_T}{\partial t}}=ϵk_{\mathrm{3}}{\rho }_{BO}-ϵk_{\mathrm{5}}{\rho }_T,

(5) {\displaystyle \frac{\partial ϵ{\rho }_G}{\partial t}}=k_{\mathrm{1}}{\rho }_B+ϵk_{\mathrm{4}}{\rho }_{BO}+ϵk_{\mathrm{5}}{\rho }_T,

where

(6) k_i=A_i\exp \left[\left(-E_i/RcT\right)\right],i=\mathrm{1},\mathrm{2},\cdots ,\mathrm{5},

(7) {\rho }_B=M_B/V,

and

(8) {\rho }_C=M_C/V.

ƒM_B is the biomass weight, M_C is the charcoal weight, V is the volume of all of the particles, ρ_B is the biomass density, and ρ_C is the charcoal density. M_BO is the biomass oil weight, M_T is the tar weight, M_G is the non-condensable gas weight, V_g is the volume of both the gas and volatiles produced by the particle, ρ_BO is the mass concentration of biomass oil, ρ_T is the mass concentration of tar, and ρ_G is the mass concentration of non-condensable gases.

(9) {\rho }_{BO}=M_{BO}/V_g=M_{BO}/\left(ϵV\right),

(10) {\rho }_T=M_T/V_g=M_T/\left(ϵV\right),

(11) {\rho }_G=M_G/V_g=M_G/\left(ϵV\right),

where ϵ is the porosity of the medium,

(12) ϵ=V_g/V.

Mathematical model

Two steps are involved in achieving pyrolysis, the initial pyrolysis reaction, and the secondary thermal cracking reaction. Both pyrolysis and combustion processes take place in a gasification process. Heat is delivered to the biomass particle surface by convection and radiation, and then transferred to the interior. The increasing temperature results in pre-pyrolysis and pyrolysis, and the chemical reactions and phase changes lead to heat changes. The temperature is, therefore, a function of the conversion time. The volatile and gaseous products flow through the particle pores and participate in the heat transfer process. The pyrolysis reaction rate also depends on the temperature of the particle environment. In our study, because tar condensation was not considered, further assumptions were as follows:

(a) Every substance in the biomass particle maintains the same volume during pyrolysis. Thermal expansion and contraction are neglected, i.e., the particle volume is considered constant.

(b) Gases and solids are in local thermal equilibrium, that is, the temperature of the gas is the same as that of the solid.

(c) The volatile matter follows the Darcy flow, without diffusion flow during pyrolysis. That is, the effect of molecular diffusion on energy change is not considered.

(d) Particles are well distributed and either spherical or cylindrical.

(e) The heat flux is one-dimensional.

(f) Drying is a physical process, which will not have a material impact on the pyrolysis reaction. Therefore, it is assumed that the material does not contain moisture.

Mass conservation equation

(13) h{\displaystyle \frac{\partial \rho }{\partial t}}+{\displaystyle \frac{h_g}{r}}{\displaystyle \frac{\partial }{\partial r}}\left(r{\rho }_{g}u\right)=h_g{\omega }_g-h_s{\omega }_g={\omega }_g\left(h_g-h_s\right)={\omega }_g\mathrm{\Delta }h.

The detailed derivation for this equation can be found in the literature ( Qi, 2010; Qi et al., 2011).

The Darcy law

The momentum conservation equation can be replaced by the Darcy law when we predict the porous medium gas speed. The Darcy law is described by

(14) u=-{\displaystyle \frac{\kappa }{\mu }}{\displaystyle \frac{\partial P}{\partial r}}.

We used the ideal gas equation to calculate the gas pressure:

(15) P={\rho }_g{R}_{g}T/M_g,

where M_g is the average molecular weight of the gaseous phase.

The gas density is for the whole gaseous phase, calculated as

(16) {\rho }_g={\displaystyle \sum _{i=\mathrm{1}}^{n_g}{\rho }_i}.

Energy conservation equation

(17) \left[{\rho }_C{c}_C+{\rho }_B{c}_B+ϵ{\rho }_g{c}_g\right]{\displaystyle \frac{\partial T}{\partial t}}=\left(-{\rho }_g{c}_{g}u+k_{eff}{\displaystyle \frac{p}{r}}\right){\displaystyle \frac{\partial T}{\partial r}}+{\displaystyle \frac{\partial }{\partial r}}\left(k_{eff}{\displaystyle \frac{\partial T}{\partial r}}\right)+q{\omega }_g.

The detailed derivation of this equation can be found in the literature ( Qi, 2010; Qi et al., 2011).

Initial and boundary conditions

The initial conditions for Eq. (17) were

(18) t=\mathrm{0}\mathrm{;}T=T_{\mathrm{0}},{\rho }_B={\rho }_{B\mathrm{0}},{\rho }_C={\rho }_T={\rho }_G={\rho }_{BO}=\mathrm{0},u=\mathrm{0.}

At the initial conditions, the solid is in a static environment at atmospheric pressure.

The boundary conditions for Eq. (17) were

(19) \begin{array}{l}t>\mathrm{0},r=\mathrm{0}\hfill \\ {\left({\displaystyle \frac{\partial T}{\partial r}}\right)}_{r=\mathrm{0}}=\mathrm{0};u=\mathrm{0}\hfill \end{array},

(20) \begin{array}{l}t>\mathrm{0},r=R\hfill \\ -k_{eff}{\displaystyle \frac{\partial T}{\partial r}}={\alpha }_{\mathrm{0}}\sigma (T_f^{\mathrm{4}}-T^{\mathrm{4}})+h(T_f-T);P=P_{\mathrm{0}}\hfill \end{array}.

Convective and radiative transfers at the surface are considered in Eq. (20).

Dimensionless forms of the energy conservation equation

The dimensionless form of Eq. (17) is

(21) {\displaystyle \frac{\partial \vartheta }{\partial \tau }}=\left(-\mathrm{\Omega }+{\displaystyle \frac{p}{x}}\right){\displaystyle \frac{\partial \vartheta }{\partial x}}+\mathrm{\Pi }{\displaystyle \frac{{\partial }^{\mathrm{2}}\vartheta }{\partial x^{\mathrm{2}}}}+\mathrm{\Theta }.

The dimensionless forms of the initial conditions are

(22) \tau =\mathrm{0}\mathrm{;}\theta \left(x,\mathrm{0}\right)=\mathrm{0},{\rho }_B={\rho }_{B\mathrm{0}},{\rho }_C={\rho }_T={\rho }_G=\mathrm{0},u=\mathrm{0.}

The dimensionless forms of the boundary conditions are

(23) \tau >\mathrm{0}\mathrm{;}x=\mathrm{0},{\displaystyle \frac{\partial \theta }{\partial x}}=\mathrm{0},

(24) \tau >\mathrm{0}\mathrm{;}x=\mathrm{1},{\displaystyle \frac{\partial \theta }{\partial \mathrm{X}}}=-\theta \mathrm{\Phi }.

The dimensionless parameters are

(25) \vartheta ={\displaystyle \frac{T-T_f}{T_{\mathrm{0}}-T_f}},

(26) x={\displaystyle \frac{r}{R}},

(27) \tau ={\displaystyle \frac{\alpha t}{R^{\mathrm{2}}}},

(28) \mathrm{\Omega }={\displaystyle \frac{R{\rho }_g{c}_{g}u}{\alpha \rho c_p}},

(29) \mathrm{\Pi }={\displaystyle \frac{k_{eff}}{\alpha \rho c_p}},

(30) \mathrm{\Theta }={\displaystyle \frac{q{\omega }_g{R}^{\mathrm{2}}}{\left(T_{\mathrm{0}}-T_f\right)\alpha \rho c_p}},

(31) \mathrm{\Phi }={\displaystyle \frac{R}{k_{eff}}}\left[h+{\alpha }_{\mathrm{0}}\sigma \left(T^{\mathrm{3}}+T^{\mathrm{2}}{T}_f+T_f^{\mathrm{2}}T+T_f^{\mathrm{3}}\right)\right].

In accordance with L Hospital’s rule:

(32) \underset{x\to \mathrm{0}}{\lim }{\displaystyle \frac{\mathrm{1}}{x}}{\displaystyle \frac{\partial \vartheta }{\partial x}}={\displaystyle \frac{{\partial }^{\mathrm{2}}\vartheta }{\partial x^{\mathrm{2}}}}.

Thus, at x = 0, Eq. (21) is converted to

(33) {\displaystyle \frac{\partial \vartheta }{\partial \tau }}=-\mathrm{\Omega }{\displaystyle \frac{\partial \vartheta }{\partial x}}+\left(p+\mathrm{\Pi }\right){\displaystyle \frac{{\partial }^{\mathrm{2}}\vartheta }{\partial x^{\mathrm{2}}}}+\mathrm{\Theta }.

Numerical solution of the pyrolysis model for granular biomass

A simulated granular biomass particle is shown in Fig. 2, and the initial density, temperature, and pressure for the particle are given. The boundary conditions for pyrolysis of the particle can be seen in Eqs. (18)–(20).

The energy equation, initial and boundary dimensions, and calculation process are shown in Eqs. (34)–(43).

Discrete convection calculations were performed using first-order backward differences, and diffusion calculations were performed using the central difference method. The finite difference scheme equations are as follows:

(34) {\displaystyle \frac{{\theta }_{\mathrm{1}}^{n+\mathrm{1}}-{\theta }_{\mathrm{1}}^n}{\mathrm{\Delta }\tau }}=\mathrm{2}\left(p+\mathrm{\Pi }\right){\displaystyle \frac{{\theta }_{\mathrm{2}}^{n+\mathrm{1}}-{\theta }_{\mathrm{1}}^{n+\mathrm{1}}}{{\left(\mathrm{\Delta }X\right)}^{\mathrm{2}}}}+\mathrm{\Theta },

(35) {\displaystyle \frac{{\theta }_i^{n+\mathrm{1}}-{\theta }_i^n}{\mathrm{\Delta }\tau }}=-\left(\mathrm{\Omega }+{\displaystyle \frac{p}{x_i}}\right){\displaystyle \frac{{\theta }_{i+\mathrm{1}}^{n+\mathrm{1}}-{\theta }_{i-\mathrm{1}}^{n+\mathrm{1}}}{\mathrm{2}\mathrm{\Delta }X}}+\mathrm{\Pi }{\displaystyle \frac{{\theta }_{i+\mathrm{1}}^{n+\mathrm{1}}-\mathrm{2}{\theta }_i^{n+\mathrm{1}}+{\theta }_{i-\mathrm{1}}^{n+\mathrm{1}}}{{\left(\mathrm{\Delta }X\right)}^{\mathrm{2}}}}+\mathrm{\Theta },

(36) {\displaystyle \frac{{\theta }_M^{n+\mathrm{1}}-{\theta }_M^n}{\mathrm{\Delta }\tau }}=-\left(\mathrm{\Omega }+p\right)\mathrm{\Phi }{\theta }_M^{n+\mathrm{1}}+\mathrm{2}\mathrm{\Pi }{\displaystyle \frac{\left(\mathrm{\Delta }X\mathrm{\Phi }-\mathrm{1}\right){\theta }_M^{n+\mathrm{1}}+{\theta }_{M-\mathrm{1}}^{n+\mathrm{1}}}{{\left(\mathrm{\Delta }X\right)}^{\mathrm{2}}}}+\mathrm{\Theta }.

To facilitate the use of the numerical method, Eqs. (34)–(36) were completed as follows:

(37) \left({\displaystyle \frac{\mathrm{2}\mathrm{\Pi }\mathrm{\Delta }\tau }{{\left(\mathrm{\Delta }X\right)}^{\mathrm{2}}}}+\mathrm{1}\right){\theta }_{\mathrm{1}}^{n+\mathrm{1}}+\left(-{\displaystyle \frac{\mathrm{2}\left(p+\mathrm{\Pi }\right)\mathrm{\Delta }\tau }{{\left(\mathrm{\Delta }X\right)}^{\mathrm{2}}}}\right){\theta }_{\mathrm{2}}^{n+\mathrm{1}}={\theta }_{\mathrm{1}}^n+\mathrm{\Theta }\mathrm{\Delta }\tau ,

(38) \left(-\left(\mathrm{\Omega }+{\displaystyle \frac{p}{x_i}}\right){\displaystyle \frac{\mathrm{\Delta }\tau }{\mathrm{2}\mathrm{\Delta }X}}-{\displaystyle \frac{\mathrm{\Pi }\mathrm{\Delta }\tau }{{\left(\mathrm{\Delta }X\right)}^{\mathrm{2}}}}\right){\theta }_{i-\mathrm{1}}^{n+\mathrm{1}}+\left({\displaystyle \frac{\mathrm{2}\mathrm{\Pi }\mathrm{\Delta }\tau }{{\left(\mathrm{\Delta }X\right)}^{\mathrm{2}}}}+\mathrm{1}\right){\theta }_i^{n+\mathrm{1}} +\left(\left(\mathrm{\Omega }+{\displaystyle \frac{p}{x_i}}\right){\displaystyle \frac{\mathrm{\Delta }\tau }{\mathrm{2}\mathrm{\Delta }X}}-{\displaystyle \frac{\mathrm{\Pi }\mathrm{\Delta }\tau }{{\left(\mathrm{\Delta }X\right)}^{\mathrm{2}}}}\right){\theta }_{i+\mathrm{1}}^{n+\mathrm{1}}={\theta }_i^n+\mathrm{\Theta }\mathrm{\Delta }\tau ,

(39) \left(-{\displaystyle \frac{\mathrm{2}\mathrm{\Pi }\mathrm{\Delta }\tau }{{\left(\mathrm{\Delta }X\right)}^{\mathrm{2}}}}\right){\theta }_{M-\mathrm{1}}^{n+\mathrm{1}}+\left({\displaystyle \frac{\mathrm{2}\mathrm{\Pi }\mathrm{\Delta }\tau }{{\left(\mathrm{\Delta }X\right)}^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{2}\mathrm{\Pi }\mathrm{\Phi }\mathrm{\Delta }\tau }{\mathrm{\Delta }X}}+\mathrm{\Omega }\mathrm{\Phi }\mathrm{\Delta }\tau \right){\theta }_M^{n+\mathrm{1}}={\theta }_M^n+\mathrm{\Theta }\mathrm{\Delta }\tau .

Eq. (37) can be reorganized into the following form,

(40) b_{\mathrm{1}}{\theta }_{\mathrm{1}}^{n+\mathrm{1}}+c_{\mathrm{1}}{\theta }_{\mathrm{2}}^{n+\mathrm{1}}={\theta }_{\mathrm{1}}^n+\mathrm{\Theta }\mathrm{\Delta }\tau ,

where, b_{\mathrm{1}}=\left({\displaystyle \frac{\mathrm{2}\mathrm{\Pi }\mathrm{\Delta }\tau }{{\left(\mathrm{\Delta }X\right)}^{\mathrm{2}}}}+\mathrm{1}\right),c_{\mathrm{1}}=\left(-{\displaystyle \frac{\mathrm{2}\left(p+\mathrm{\Pi }\right)\mathrm{\Delta }\tau }{{\left(\mathrm{\Delta }X\right)}^{\mathrm{2}}}}\right) .

Eq. (38) can be reorganized into the following form

(41) a_m{\theta }_{i-\mathrm{1}}^{n+\mathrm{1}}+b_m{\theta }_i^{n+\mathrm{1}}+c_m{\theta }_{i+\mathrm{1}}^{n+\mathrm{1}}={\theta }_i^n+\mathrm{\Theta }\mathrm{\Delta }\tau ,

where

\begin{array}{l}{a}_m=\left(-\left(\mathrm{\Omega }+{\displaystyle \frac{p}{x_i}}\right){\displaystyle \frac{\mathrm{\Delta }\tau }{\mathrm{2}\mathrm{\Delta }X}}-{\displaystyle \frac{\mathrm{\Pi }\mathrm{\Delta }\tau }{{\left(\mathrm{\Delta }X\right)}^{\mathrm{2}}}}\right),b_m=\left({\displaystyle \frac{\mathrm{2}\mathrm{\Pi }\mathrm{\Delta }\tau }{{\left(\mathrm{\Delta }X\right)}^{\mathrm{2}}}}+\mathrm{1}\right),\hfill \\ c_m=\left(\left(\mathrm{\Omega }+{\displaystyle \frac{p}{x_i}}\right){\displaystyle \frac{\mathrm{\Delta }\tau }{\mathrm{2}\mathrm{\Delta }X}}-{\displaystyle \frac{\mathrm{\Pi }\mathrm{\Delta }\tau }{{\left(\mathrm{\Delta }X\right)}^{\mathrm{2}}}}\right).\hfill \end{array}

Eq. (39) can be reorganized into the following form

(42) a_M{\theta }_{M-\mathrm{1}}^{n+\mathrm{1}}+b_M{\theta }_M^{n+\mathrm{1}}={\theta }_M^n+\mathrm{\Theta }\mathrm{\Delta }\tau ,

where a_M=\left(-{\displaystyle \frac{\mathrm{2}\mathrm{\Pi }\mathrm{\Delta }\tau }{{\left(\mathrm{\Delta }X\right)}^{\mathrm{2}}}}\right){\theta }_{M-\mathrm{1}}^{n+\mathrm{1}},b_M=\left({\displaystyle \frac{\mathrm{2}\mathrm{\Pi }\mathrm{\Delta }\tau }{{\left(\mathrm{\Delta }X\right)}^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{2}\mathrm{\Pi }\mathrm{\Phi }\mathrm{\Delta }\tau }{\mathrm{\Delta }X}}+\mathrm{\Omega }\mathrm{\Phi }\mathrm{\Delta }\tau \right). ƒEqs. (40)–(42) show that the discrete energy equation can be expressed as Eq. (43)

(43) \left\{\begin{array}{lllll}{b}_1\hfill & c_1\hfill & \hfill & \hfill & \hfill \\ a_2\hfill & b_2\hfill & c_2\hfill & \hfill & \hfill \\ \hfill & \cdots \hfill & \cdots \hfill & \cdots \hfill & \hfill \\ \hfill & \hfill & \cdots \hfill & \cdots \hfill & \cdots \hfill \\ \hfill & \hfill & a_{M-1}\hfill & b_{M-1}\hfill & c_{M-1}\hfill \\ \hfill & \hfill & \hfill & a_M\hfill & b_M\hfill \end{array}\right\}\vartheta =\left\{\begin{array}{c}{\vartheta }_1^n+\mathrm{\Theta }\mathrm{\Delta }\tau \hfill \\ {\vartheta }_2^n+\mathrm{\Theta }\mathrm{\Delta }\tau \hfill \\ ⋮\hfill \\ ⋮\hfill \\ {\vartheta }_{M-1}^n+\mathrm{\Theta }\mathrm{\Delta }\tau \hfill \\ {\vartheta }_M^n+\mathrm{\Theta }\mathrm{\Delta }\tau \hfill \end{array}\right\}.

The matrix in Eq. (43) contains M linear algebraic equations and M unknown quantities. The discrete equations were solved using the tridiagonal matrix algorithm (TDMA). Δτ=0.01, M = 201, and Δx=1/M were chosen for the simulation. Based on the initial conditions, the temperature distribution was obtained from numerical calculations using the heat transfer model; and then the temperature distribution was put into the chemical kinetic equation from Eqs. (1)–(5), which were then solved using a fixed step Runge-Kutta method. Because the control equation is nonlinear, the unknown quantities had to be assessed by iteration.

Experimental verification

For the discrete equations ( Qi, 2010), we used the Fortran programming language for the numerical calculations necessary to obtain simulated data. The simulated data were compared with experimental data published by Pyle and Zaror (1984), who investigated the pyrolysis process for cylindrical wood pellets. The distribution of wood pellet mass with temperature, and the time spent at a constant temperature were included in the investigation. Pyle and Zaror used pine wood in their tests, with diameters of 0.6 cm, 1.5 cm, and 2.2 cm, and a length of 6–9 cm so that heat conduction in the axial direction could be disregarded. The mass and temperature changes in the furnace during the biomass pyrolysis process were recorded in the temperature range of 623–780 K. The kinetic parameters used in the calculations can be found in Tables 1 and 2, in which E₃–E₅ are activation energies, and the pre-exponential factors A₃–A₅ are derived from the linear regression between the calculated values and experimental data.

Comparison of the computational and experimental results for the mass conversion yield in biomass pyrolysis

The mass conversion yield of the biomass particle as a function of time is shown in Fig. 3 (the diameter was 6 mm, the initial temperature was 303 K, and the final temperature was 643 K). The results show that the simulated and experimental data are in agreement, apart from between 1 and 2 min. In general, the simulation curve is in good agreement with the test values.

Temperature distribution for biomass pyrolysis

The temperature profile as a function of time for biomass pyrolysis using a particle diameter of 15 mm is shown in Fig. 4. Both the simulated and experimental data show the center of the particle to have an initial temperature of 303 K, and to reach a final temperature of 660 K. The simulated data are very consistent with the experimental data at both low temperatures and at the final temperature. There is, however, a little difference between the simulated and experimental data, which may have been caused by the sharp increase in temperature caused by the fast exothermic pyrolysis reaction. The main mode of heat transfer to the center of the particle at 303 K is conduction, and the biomass starts to pyrolyse as the temperature increases to 450 K. There is a delay when the endothermic pyrolysis reaction occurs, and the exothermic reaction appears when secondary cracking of the condensable volatiles occurs, at 550 K.

Numerical simulation results

Internal temperature distribution during pyrolysis of biomass particles

Radial temperature profile in biomass particles over time

The biomass temperature profile by radial distance, as a function of pyrolysis time (5 s, 10 s, 20 s, 30 s, 40 s, and 50 s), with an initial temperature of 303 K, a final temperature of 1,473 K, and a particle radius of 4 mm, is shown in Fig. 5. The temperature at each radial distance clearly increases with pyrolysis time. The temperature at the surface increases much faster than the temperature at the center of the particle. In the initial stages of pyrolysis (5 s), the temperature distribution curve around the particle surface is steeper than the other three curves shown because of the higher heat transfer resistance near the surface. The temperature at the center stays almost constant at 5 s and 10 s. However, there is a temperature rise at 20 s, and at 50 s the temperature at the center and at the surface are almost the same, which shows that the pyrolysis reaction has nearly finished. The six temperature profiles shown in Fig. 5 demonstrate that the slope is gentler with increasing temperature. That is, the temperature difference between the center and the surface becomes smaller, and the heat transfer rate slows down with time.

Radial distance temperature profile as a function of conversion time

The internal temperature profiles as a function of conversion time for particles with different diameters are shown in Figs. 6 and 7. The initial temperature was 303 K, the environmental temperature was 1,473 K, and the particle diameters were 1 mm (Fig. 6) and 32 mm (Fig. 7). The initial temperature of the small particles is well distributed, and the temperature difference between the interior and exterior is less than 500 K. However, the temperature difference between the interior and exterior of the large particles reaches 1,000 K during the initial stages of pyrolysis. There is an obvious difference in the temperature distribution between large and small particles, because large particles have a thicker geometric structure with a low coefficient of heat conductivity.

The pyrolysis reaction zone appears as a period at constant temperature in several layers near the center of the large biomass particles. Before and after that period, the pure physical heating processes were observed, following the temperature index curve, which is embodied in the endothermic effect on the temperature field and the reaction process. However, there is little endothermic effect on the temperature field for large particles. Therefore, the chemical reaction is a speed control mechanism for only the chemical reaction in small particles; but it is a speed control mechanism for both internal heat transfer and the chemical reaction in large particles. The upper section of the temperature curve describes the temperature variation for char particles. Because of low porosity and radiative heat transfer, large particles have stronger thermal resistance and internal temperature gradients, which will have a significant influence on the secondary cracking of condensable volatiles.

Pressure as a function of radial distance for biomass pyrolysis

Pressure as a function of radial distance at different times is shown in Fig. 8 for particles of diameter 4 mm, an initial temperature of 303 K, and an environmental temperature of 1,473 K. A large number of factors affect the internal pressure of the particle, and the main cause of overpressure is the production of gas phase products. The main reason for a decrease in pressure is pressure-driven gas phase loss through pores. The pressure distribution maximum appears adjacent to the reaction boundary with the unreacted core because of the low intrinsic penetration rate. The maximum appears because of high gas phase product concentrations despite the low temperature. The low intrinsic penetration rate, low porosity, and almost zero convection diffusion lead to this pressure distribution maximum. In addition, the low intrinsic penetration rate combined with the pressure gradient results in a low calculated gas velocity. The low mass flow rate results in increasing flow resistance in the pores, so the pressure increases with time. The intrinsic penetration rate for char is six orders of magnitude higher than the intrinsic penetration rate for the biomass raw material. When the gas phase products flow to the surface, resistance declines as soon as a small amount of char is generated. At that time, the gas phase pyrolysis products accumulated in the particle interior flow out quickly, and the internal pressure of the particle decreases rapidly. The porosity increases with the quantity of gas phase pyrolysis products. The internal pressure of the particle gradually decreases to the ambient pressure when a high intrinsic penetration rate is achieved.

Condensable volatile flux as a function of conversion time

The flux of condensable volatiles as a function of conversion time is shown in Fig. 9 for a particle diameter of 4 mm, an initial temperature of 303 K, and an environmental temperature of 1,473 K. It can be seen that the mass flux of condensable volatiles rapidly decreases to 1.2 kg/(m²·s), increases slightly to 0.8 kg/(m²·s), and then gradually decreases to zero with time. Two peaks can be seen in Fig. 9. The first peak appears as biomass oil is generated during the first pyrolysis process. The second peak appears when the condensable volatiles are cracked and polymerized to noncondensable gas and char during the second pyrolysis process. More biomass oil than tar is generated, and, therefore, the first peak is significantly higher than the second peak.

Conclusions

1) A high-temperature biomass pyrolysis model was established. Based on the characteristics of the liquid products, the condensable volatiles were divided into two types: tar and biomass oil. The effects of temperature, residence time, particle size, velocity, pressure, and other parameters on biomass pyrolysis and high-temperature tar cracking were investigated numerically. We found that the model was effective for describing the biomass and tar pyrolysis processes by comparing it with experimental data.

2) The endothermic process affects the temperature and the reaction process. There is a small endothermic effect on temperature for small particles. Therefore, the chemical reaction is a speed control mechanism for only the chemical reaction in small particles. However, it is a speed control mechanism for both the internal heat transfer and the chemical reaction in large particles.

3) A pressure distribution maximum appears adjacent to the boundary between the reaction and the unreacted core. The intrinsic penetration rate for char is six orders of magnitude higher than that of the raw biomass materials. Gas phase pyrolysis products accumulated in the particle interior flow out quickly, and the internal pressure of the particle decreases rapidly. The internal pressure of the particle gradually decreases to the ambient pressure when the intrinsic penetration rate becomes high.

4) There are two peaks in the pyrolysis gas mass fraction curve. The first peak appears as biomass oil is cracked into non-condensable gases and tar, and the second peak appears as tar is cracked into non-condensable gases and coke.

References

[1] Babu B V, Chaurasia A S (2004). Heat transfer and kinetics in the pyrolysis of shrinking biomass particle. Chem Eng Sci, 59(10): 1999–2012

[2] Chen B, Chen G Q (2006a). Exergy analysis for resource conversion of the Chinese society 1993 under the material product system. Energy, 31(8–9): 1115–1150

[3] Chen B, Chen G Q (2006b). Ecological footprint accounting based on emergy—A case study of the Chinese society. Ecol Modell, 198(1–2): 101–114

[4] Chen B, Chen G Q, Yang Z F (2006a). Exergy-based resource accounting for China. Ecol Modell, 196(3–4): 313–328

[5] Chen B, Chen G Q, Yang Z F, Jiang M M (2007). Ecological footprint accounting for energy and resource in China. Energy Policy, 35(3): 1599–1609

[6] Chen B, He G X, Qi J, Su M R, Zhou S Y, Jiang M M (2012a). Greenhouse gas inventory of a typical high-end industrial park in China. ScientificWorldJournal, doi: 10.1155/2013/717054

[7] Chen B, He G X, Yang J, Zhang J R, Su M R, Qi J (2012b). Evaluating ecological and economic benefits of a low-carbon industrial park based on millennium ecosystem assessment framework. ScientificWorldJournal, 2012, 909317: 1–5

[8] Chen G Q, Chen B (2009). Extended exergy analysis of the Chinese society. Energy, 34(9): 1127–1144

[9] Chen G Q, Jiang M M, Chen B, Yang Z F, Lin C (2006b). Emergy analysis of Chinese agriculture. Agric Ecosyst Environ, 115(1–4): 161–173

[10] Chen S Q, Chen B (2011). Assessing inter-city ecological and economic relations: an emergy-based conceptual model. Front Earth Sci, 5(1): 97–102

[11] Chen S Q, Chen B (2012a). Defining indirect uncertainty in system-based risk management. Ecol Inform, 10: 10–16

[12] Chen S Q, Chen B (2012b). Network environ perspective for urban metabolism and carbon emissions: a case study of Vienna, Austria. Environ Sci Technol, 46(8): 4498–4506

[13] Chen S Q, Chen B (2012c). Sustainability and future alternatives of biogas-linked agrosystem (BLAS) in China: an emergy analysis. Renew Sustain Energy Rev, 16(6): 3948–3959

[14] Chen S Q, Chen B, Song D (2012c). Life-cycle energy production and emissions mitigation by comprehensive biogas–digestate utilization. Bioresour Technol, 114: 357–364

[15] Chen Z M, Chen G Q, Zhou J B, Jiang M M, Chen B (2010). Ecological input-output modeling for embodied resources and emissions in Chinese economy. Commun Nonlinear Sci Numer Simul, 15(7): 1942–1965

[16] Dai J, Fath B D, Chen B (2012). Constructing a network of the Social-economic Consumption System of China using Extended Exergy Analysis. Renew Sustain Energy Rev, 16(7): 4796–4808

[17] Di Blasi C (1993a). Modeling and simulation of combustion processes of charring and non-charring solid fuels. Pror Energy Combust Sci, 19(1): 71–104

[18] Di Blasi C (1993b). Analysis of convection and secondary reaction effects within porous solid fuels undergoing pyrolysis. Combust Sci Technol, 90(5): 315–340

[19] Di Blasi C (1996). Heat, momentum and mass transport through a shrinking biomass particle exposed to thermal radiation. Chem Eng Sci, 51(7): 1121–1132

[20] Di Blasi C (2000). Modelling the fast pyrolysis of cellulosic particles in fluid-bed reactors. Chem Eng Sci, 55(24): 5999–6013

[21] Di Blasi C (2008). Modeling chemical and physical processes of wood and biomass pyrolysis. Pror Energy Combust Sci, 34(1): 47–90

[22] Feng J K, Shen Y T, Yang R C (2003). Principle and Calculation of the Boiler (3rd ed). Beijing: Science Press, 691

[23] Grønli M G, Melaaen M C (2000). Mathematical model for wood pyrolysis—Comparison of experimental measurements with model predictions. Energy Fuels, 14(4): 791–800

[24] Hagge M J, Bryden K M (2002). Modeling the impact of shrinkage on the pyrolysis of dry biomass. Chem Eng Sci, 57(14): 2811–2823

[25] Hubacek K, Feng K S, Chen B (2012). Changing lifestyles towards a low carbon economy: an IPAT analysis for China. Energies, 5(12): 22–31

[26] Ji X, Chen G Q, Chen B, Jiang M M (2009). Exergy-based assessment for waste gas emissions from Chinese transportation. Energy Policy, 37(6): 2231–2240

[27] Jiang M M, Chen B (2011). Integrated urban ecosystem evaluation and modeling based on embodied cosmic exergy. Ecol Modell, 222(13): 2149–2165

[28] Jiang M M, Chen B, Zhou J B, Tao F R, Li Z, Yang Z F, Chen G Q (2007). Emergy account for biomass resource exploitation by agriculture in China. Energy Policy, 35(9): 4704–4719

[29] Jiang M M, Zhou J B, Chen B, Chen G Q (2008). Emergy-based ecological account for the Chinese economy in 2004. Commun Nonlinear Sci Numer Simul, 13(10): 2337–2356

[30] Jiang M M, Zhou J B, Chen B, Yang Z F, Ji X, Zhang L X, Chen G Q (2009). Ecological evaluation of Beijing economy based on emergy indices. Commun Nonlinear Sci Numer Simul, 14(5): 2482–2494

[31] Ju L P, Chen B (2011). Embodied energy and emergy evaluation of a typical biodiesel production chain in China. Ecol Modell, 222(14): 2385–2392

[32] Kansa E J, Perlee H E, Chaiken R F (1977). Mathematical model of wood pyrolysis including internal forced convection. Combust Flame, 29(3): 311–324

[33] Koufopanos C A, Maschio G, Lucchesi A (1989). Kinetic modelling of the pyrolysis of biomass and biomass components. Can J Chem Eng, 67(2): 75–84

[34] Koufopanos C A, Papayannakos N, Maschio G, Lucchesi A (1991). Modelling of the pyrolysis of biomass particles. Studies on kinetics, thermal and heat transfer effects. Can J Chem Eng, 69(4): 907–915

[35] Li C, Suzuki K (2009). Tar property, analysis, reforming, mechanism and model for biomass gasification-an overview. Renew Sustain Energy Rev, 13(3): 594–604

[36] Liu G Y, Yang Z F, Chen B (2011a). Extended exergy analysis of urban socioeconomic system: a case study of Beijing, 1996–2006. International Journal of Exergy, 168–191

[37] Liu G Y, Yang Z F, Chen B, Ulgiati S (2011b). Monitoring trends of urban development and environmental impact of Beijing, 1999–2006. Sci Total Environ, 409(18): 3295–3308

[38] Liu G Y, Yang Z F, Su M R, Chen B (2012). The structure, evolution and sustainability of urban socio-economic system. Ecol Inform, 10: 2–9

[39] Lu Y, Su M R, Liu G Y, Chen B, Zhou S Y, Jiang M M (2012). Ecological network analysis for a low-carbon and high-tech industrial park. ScientificWorldJournal, 2012, 305474, 1–9

[40] Maniatis K, Beenackers A A C M (2000). Tar protocols. IEA bioenergy gasification task. Biomass and Bioenergy, 18(1): 1–4

[41] Melaaen M C (1996). Numerical analysis of heat and mass transfer in drying and pyrolysis of porous media. Numerical Heat Transfer, Part A: Applications, 29(4): 331–355

[42] Piskorz J, Majerski P, Radlein D, Scott D S, Bridgwater A V (1998). Fast pyrolysis of sweet sorghum and sweet sorghum bagasse. J Anal Appl Pyrolysis, 46(1): 15–29

[43] Pyle D L, Zaror C A (1984). Heat transfer and kinetics in the low temperature pyrolysis of solids. Chem Eng Sci, 39(1): 147–158

[44] Qi G L (2010). Experimental research and numerical simulation of biomass pyrolysis and tar thermal cracking. Dissertation for Ph.D degree. Harbin Institute of Technology, 76–85

[45] Qi G L, Dong P, Zhang Y (2011). Construction of large biomass high temperature pyrolysis model and numerical simulation. Acta Energiae Solaris Sinica, 32(7): 1058–1063

[46] Qu Y H, Lin C, Zhou W, Li Y, Chen B, Chen G Q (2009a). Effects of CO₂ concentration and moisture content of sugar-free media on the tissue-cultured plantlets in large growth chamber. Commun Nonlinear Sci Numer Simul, 14(1): 322–330

[47] Qu Y H, Wei X M, Hou Y F, Chen B, Chen G Q, Lin C (2009b). Analysis for an environmental friendly seedling breeding system. Commun Nonlinear Sci Numer Simul, 14(4): 1766–1772

[48] Song D, Su M R, Yang J, Chen B (2012). Greenhouse gas emission accounting and management of low-carbon community. ScientificWorldJournal, 2012, 6137211–6

[49] Su M, Liang C, Chen B, Chen S, Yang Z (2012). Low-carbon Development patterns: observations of typical Chinese cities. Energies, 5(2): 291–304

[50] Su M R, Yang Z F, Chen B (2009). Set pair analysis for urban ecosystem health assessment. Commun Nonlinear Sci Numer Simul, 14(4): 1773–1780

[51] Su M R, Yang Z F, Chen B (2011). An emergy-based analysis of urban ecosystem health characteristics for Beijing city. International Journal of Exergy, 192–209

[52] Wei X M, Chen B, Qu Y H, Lin C, Chen G Q (2009). Emergy analysis for ‘Four in One’ peach production system in Beijing. Commun Nonlinear Sci Numer Simul, 14(3): 946–958

[53] Yang J, Chen B (2011). Using LMDI method to analyze the change of industrial CO₂ emission from energy use in Chongqing. Front Earth Sci, 103–115

[54] Yang J, Chen B, Qi J, Zhou S Y, Jiang M M (2012). Life-cycle-based multicriteria sustainability evaluation of industrial parks: a case study in China. ScientificWorldJournal, 2012917830

[55] Yang J, Chen W C, Chen B (2011). Impacts of biogas projects on agro-ecosystem in rural areas — A case study of Gongcheng. Front Earth Sci, 5(3): 317–322

[56] Yang Q, Chen B, Ji X, He Y F, Chen G Q (2009). Exergetic evaluation of corn-ethanol production in China. Commun Nonlinear Sci Numer Simul, 14(5): 2450–2461

[57] Zeng L, Luo Z L, Chen B, Yang Z F, Li Z, Lin W X, Chen G Q (2010). Numerical analysis of a lock-release oil slick. Commun Nonlinear Sci Numer Simul, 15(8): 2222–2230

[58] Zhang L X, Feng Y Y, Chen B (2011). Alternative scenarios for the development of a low-carbon city: a case study of Beijing, China. Energies, 4(12): 2295–2310

[59] Zhang L X, Song B, Chen B (2012). Emergy-based analysis of four farming systems: insight into agricultural diversification in rural China. J Clean Prod, 28: 33–44

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg