Modeling for Pyrolysis of Solid Biomass

In comparison to coal, biomass is characterized by a higher content of volatile matter. It is a renewable source of energy which has many advantages from an ecological point of view. Understanding the physical phenomena of pyrolysis and representing them with a mathematical model is the primary step in the design of pyrolysis reactors. In the present study, an existing mathematical model is used to describe the pyrolysis of a single solid particle of biomass. It couples the heat transfer equations with the chemical kinetics equations. A finite difference method is used for solving the heat transfer equation and the two-step pyrolysis kinetics equations. The model equation is solved for a slab particle of equivalent dimension 0.001 m and temperature ranging from 300 to 923 K. An original numerical model for the pyrolysis of wood chips is proposed and relevant equations solved using original program realized in MATLAB. To check the validity of the numerical results, experimental results of pyrolysis of woody biomass in laboratory facility was used. The samples were heated over a range of temperature from 300 to 923 K with three different heating rates of 21, 32 and 55 K/min, and the weight loss was measured. The simulation results as well as the results obtained from thermal decomposition process indicate that the temperature peaks at maximum weight loss rate change with the increase in heating rate. The experimental results showed that the simulation results are in good agreement and can be successfully used to understand the degradation mechanism of solid reaction.


Introduction
Biomass, especially wood, is one of the first sources of energy used by humans. In the developing countries it still represents a major source of energy, while the developed countries have seen a renewed interest in using biomass since the nineteen-seventies.
In Serbia, energy demand is primarily satisfied from fossil fuels. At the same time the available domestic primary energy sources cannot satisfy the country's needs and therefore 30-40 % of the annual energy needs is imported (5.31 Mtoe in 2010. of oil) [1]. On the other hand, the total renewable energy potential in Serbia is more than 4.3 Mtoe per year [2].
Unfortunately, the production of energy from renewable sources is present only in some small plants in Serbia. Biomass is not currently used for electricity generation, and its utilization is limited to some new facilities installed in food and processing industry.
In spite of that, the production of energy from renewable sources shows a significant potential for Serbia at this moment. The most important source of biomass in Serbia is in agriculture and forestry (2.7 Mtoe) [2]. Geographically, the northern part of the territory of Serbia is mostly plain agricultural area, while the southern part is a mountainous region rich in forests. The biomass from wood processing industries consists of residues and waste from cutting, grinding and planning. Industrial wood residues represent a significant energy potential in Serbia. The biomass from wood processing industries can be used as fuel in boilers and for bio-fuel production.
It is clear that biomass cannot totally replace fossil fuels, but it may be a partial answer to the problems of CO 2 emissions and oil dependency. It appears that biomass, especially wood, is a much cleaner fuel than coal. Wood has very low sulfur content. There is usually no need for desulfurization treatment of the flue-gas in wood combustion. The combustion temperature is lower due to a lower HHV, which reduces the fuel and thermal NO x formation. However, a denitrification installation might still be necessary.
On the other hand, energy density for coal is 3 to 5 times greater than for wood. Hence, for wood to be cost competitive, it is important to limit fuel transportation and storage needs. For that reason, wood power plants are usually located in forested areas.
According to the present situation of the market in Serbia (investment limitations, ways of collecting, storage and transportations of biomass, technology and prices), it is estimated that utilities with biomass pyrolysis would be the most efficient.
Many kinetic models for wood pyrolysis have been reported in the literature. A good review is given by Di Blasi [3]. The simplest models were based on a single decomposition reaction, and they do not allow to predict the influence of pyrolysis conditions on the amount of products [4]. Other models assume some parallel reactions to predict the production kinetics of gas, tar and char [5,6].
Extensive reviews based on a number of papers done by Di Blasi [3] and Guedes et al. [7] show that the pyrolysis of biomass involves a complex series of reactions that are influenced by many factors, such as heating rate, temperature, pressure, residence time, moisture, composition of biomass material and size of particles. Most of these observations which are made from reported publications have primarily focused on woody material and in a few cases on agricultural residues. Pyrolysis in wood is typically initiated at 473 K and lasts till 723-773 K depending on the species of wood.
Miljković et al. [8] developed an original numerical model for the surface combustion of straw in a moving bed and relevant equations solved using original programme realised in C++ to simulate heterogeneous characteristics of the burning process. The model provides information concerning temperature front propagation, reaction front rate and remaining gas species composition in bed.
Pešenjanski et al. [9] used a one-step global reaction to describe the thermal degradation of wheat straw samples. The influence of different factors was investigated, such as particle size, humidity levels and the heating rate in the kinetics of devolatilization.
Peters et al. [10] developed a novel kinetic reaction model for biomass pyrolysis process and implemented in Aspen Plus. The model is based on the three main building blocks of lignocellulosic biomass, cellulose, hemicellulose and lignin. The kinetic reaction model is found to be suited for predicting pyrolysis yields and product composition for any lignocellulosic biomass.
Xianjun et al. [11] established the biomass pyrolysis model and calculated a series of the biomass pyrolysis parameters.
Ranzi et al. [12] discussed in this paper a comprehensive and unifying mathematical model that describes the chemistry of fast biomass pyrolysis. Emphasis is given to the multicomponent, multiphase, multiscale nature of this problem, together with the several simplifications for both the gas and solid phase kinetic mechanisms.
Slopiecka et al. [13] investigated the kinetic study of the pyrolysis process of poplar wood using a thermogravimetric analyzer. The results obtained from thermal decomposition process indicate that there are three main stages such as dehydration, active and passive pyrolysis. In the DTG thermograms the temperature peaks at maximum weight loss rate changed with increasing heating rate.
Gašparovič et al. [14] examined the pyrolysis of wood and main wood compounds by thermogravimetry and revealed that the decomposition process of wood depends on the composition and concentration of the main components.
Grønli et al. [15] compared the thermogravimetric curves of several hardwoods and softwoods. A comparison, between both types of wood, shows that the decomposition of softwood starts at lower temperatures.
Grieco and Baldi [16] compared experimental results of wood biomass pyrolysis with suggested kinetic model, focusing on mass loss, gas and bio-oil production.
Miljković [17] presented an experimental facility which offers an effective way to simulate grate combustion in a real facility.
A different mechanism for the pyrolysis of single particle wood and the corresponding kinetic parameters has been reported in literature [3,18,19]. In the present work, a twostep kinetic model proposed by Babu and Chaurasia [18] is modified for a slab particle. The present study examines how a wood solid generates fuel volatiles in response to a prescribed time-varying heat flux. Kinetic modeling is done to find the thermal behavior of biomass. The model enables the prediction of both char and volatile production yields. The proposed model is simulated and the results are compared with the experimental results. Non-isothermal pyrolysis data are analyzed. This procedure is applied for different values of pyrolysis heating rate. By using the model proposed in this study, it is possible to predict the pyrolysis rate over a wide range of heating rates. In this work, a wide range of temperatures (300-923 K) and relevant particle dimension (1 mm) are considered.
The aim of this research is to investigate the pyrolysis behavior of wood, to check the influence of experimental conditions on kinetic parameters and to provide kinetic information for evaluation of the processing of these materials.

Mathematical modeling 2.1 Theoretical background
Biomass is a renewable source of energy and has many advantages from an ecological point of view. Thermochemical conversion (pyrolysis, gasification, combustion) represents the most common commercial utilization of biomass. Pyrolytic degradation is recognized as an effective method for the production of high yields of char or liquid / gaseous fuels. Pyrolysis is the first stage of both combustion and gasification. So, pyrolysis is not only an independent conversion technology, but also part of the gasification and combustion processes. Pyrolysis consists of a thermal degradation of the initial solid fuel (in the absence of oxygen/air) into solid (charcoal), liquid (tar and other organics) and gaseous products. Depending on the type of pyrolysis, different species can be obtained.
For the purposes of the current study, biomass pyrolysis models may be divided into two primary categories, micro-and macro-particle models. Micro-particle pyrolysis involves the thermal decomposition of virgin matter with sample sizes sufficiently small such that diffusion effects become negligible and the intensity of pyrolysis is controlled with kinetics. Thus, micro-particles are desirable in experiments focusing on identification of kinetic schemes. Critical particle size estimates for kinetic control are generally ~0.1~1 mm and decrease with increasing pyrolysis temperatures [20]. Particles larger than the critical limit are characterized by relatively large diffusion effects which can strongly affect the pyrolysis evolution due to internal and external temperature gradients, thermal inertia due to heat capacity effects as well as temperature variations resulting from endothermic (or exothermic) reactions. In this article, a micro-particle mathematical model of pyrolysis is developed i.e. a kinetic mathematical model.
Understanding the physical phenomena of pyrolysis and representing them with a mathematical model is essential in the design of pyrolysis reactors and is also necessary for the better understanding of the behavior of engines fueled by pyrolysis products.
When a solid particle of biomass is heated in an inert atmosphere the following phenomena occur. Heat is first transferred to the particle surface by radiation and/or convection and then to the inside of the particle. The temperature inside the particle increases, causing the removal of moisture that is present in the biomass particle, and after that main pyrolysis reaction takes place.
Biomass particles commonly have more irregular shapes and much larger sizes than pulverized coal. Larger particle sizes establish the potential for large internal temperature and composition gradients that complicate combustion models. Slab represents a characteristic case because wood chips can be considered as a bunch of small slabs and therefore are especially appropriate for modeling biomass particles.

Kinetic model
Pyrolysis is a thermo-chemical decomposition process, resulting in the production of a huge number of chemical compounds [21]. For engineering applications, reaction products are often lumped into three groups: gases, liquid and char, or simply into volatiles and char. They result from both primary decomposition of the solid fuel and secondary reactions of volatile condensable organic products into low-molecular weight gases and char. Different factors affect the pyrolysis rate and the composition of the product classes. Temperature, pressure and heating rate are chief operating parameters. Biomass properties (chemical composition, ash content and composition particle size and shape, density, moisture content, etc.) also play an important role. The pyrolysis products are formed from both primary decomposition of the solid biomass as well as secondary reactions of volatile condensable organic products into low-molecular weight gases, secondary tar and char [3].
The pyrolysis reaction can be described by means of the Scheme 1 [18]: reaction III

Scheme 1 Simplified reaction scheme
This simplified reaction scheme is generally accepted nowadays and the kinetic parameters presented by Babu and Chaurasia [18] are frequently quoted and used in simulations. However, a scheme can be criticized as an over-simplification of extremely complex chemical and physical phenomena. On the one hand, the corresponding mechanisms are too complex for a mathematical modeling. On the other hand, as Velten wrote: "The best model is the simplest model that still serves its purpose, that is, which is still complex enough to help solve the original problem" [22].
When biomass is heated, it decomposes into volatiles, gases and char. The products of the first reaction (volatiles and gases) further reacts with char produced by second reaction to produce volatiles, gases and char of different composition. Therefore, the primary pyrolysis products participate in secondary interactions causing a modified final product.

Heat transfer model
During the pyrolysis process, the pores of the solid are enlarged and the solid particle merely becomes more porous, the biomass converts into gases, volatiles and char.
Inside the pyrolysis particle, heat is transmitted by the following mechanisms: • conduction inside the solid particle, • convection inside the particle pores and • convection and radiation from the surface of the particle / pellet.
For simplicity, it is assumed that heat is transmitted inside the solid by conduction only. The heat transfer coefficient represents the overall effect of the above mechanisms.
Consider a slab of thickness 2L. Assume the length of the two other dimensions of the slab to be so large that the heat transfer takes place in the one direction only. This approach was previously adopted by many researchers [23].
As ρ and T are function of t, Eq. (7) can be simplified as: where: . .
The thermal properties of the solid (i.e. specific heat, density and thermal conductivity) are presumed to vary continuously from their values for the virgin wood to their values for the char.
The initial and boundary conditions for Eq. (8) are: • Initial conditions: • Boundary conditions: In Eq. (11), the external heat transfer is considered to occur by a combination of convective and radiative mechanisms.
The governing equations form a set of nonlinear partial differential equations which can be solved numerically.
The common explicit forward finite difference method is suitable in this case. This method, however, requires simultaneous solution for the unknowns at each new time step.

Numerical solution
Mathematically speaking, the system of nonlinear differential conservation equations, with defined initial and boundary conditions, does not have an analytical solution, so numerical solution becomes the main tool in the study of these kinds of problems. Therefore, a finite difference method is used to solve the foregoing governing equations. In this study the particle is treated as a continuous porous medium and numerical calculations are carried out by dividing the particle into many cells, as shown in Fig. 1.
Inside each cell the parameters of concern (e.g. temperature, percentage of gases, carbon) are assumed to be uniform. In discretizing the governing equation by a finite difference approach, the partial differential Eqs. (1), (6) and (8) are transferred into a system of ordinary equations. Due to the arrangement of particles in a packed bed, they form a void space -a porous media. Assuming a constant volume, pyrolysis process of the particles occurs, which results in the decrease of the particle density and increase of the particle porosity. The equations are discretized in a finite difference formulation and subsequently solved. To solve this system, an iterative scheme is employed. The number of grid points was chosen to provide good resolution without excessively delaying the process. In the presented model, it is assumed that during the pyrolysis process, the cell size remains constant and only its porosity changes. In order to numerically solve mathematical model, initial conditions for differential Eqs. (1) -(6) and initial and boundary conditions for differential Eq. (8), are required to complete the system of equations.

Initial conditions
At the start of the simulation process, all gradients are set to zero. The temperatures of solid particle and gas phase in the pyrolysis chamber are 300 K, while particle density is 650 kg / m 3 and gas and char density of the particle are 0 kg / m 3 .

Boundary conditions
The gas temperature around the particle is assumed to change from 300 K to 923 K and the average heating rates are: case I-21 K / min, case II-32 K / min and case III-55 K / min. Due to intensive heat transfer between the gas and solid phase, the wood particle is heated by convection and radiation from the gas phase. Thus, the energy on the particle surface can be balanced according to Eq. (11).
Relevant equations are solved using original program realized in MATLAB and program algorithm is given in Appendix A.

Model results
The non-isothermal pyrolysis data were analyzed to determine the reaction kinetics for wood particle at three heating rates. Consider a piece of wood, of thickness 1 mm, initially at ambient temperature. When it is exposed to heat flux, the initial temperature changes from 300 to 923 K, and the average heating rates are: case I-21 K / min, case II-32 K / min and case III-55 K / min.
The set of equations is solved numerically subjected to the initial and boundary conditions and the results are presented. Given the symmetry of the problem, only half a particle is simulated. Therefore, while the particle axes of symmetry (x = 0) is subjected to conditions of zero gradients on the variables, the external surface (x = L) is exposed to heating. The external heat-transfer coefficient at the particle surface is representative of the conditions in the reactors. This transfer is made by convective transport and radiation.
The simulations have been made with non-shrinking slab-shaped particles to simulate the effects of reaction temperature. Similar approach could be used to simulate Heat transfer by convection and radiation (i=0,j=0) (i=1,j=0) (i=n/2,j=0) (i=0,j=1) Fig. 1 Cross section grid diagram of the particle the effects of particle size, external heat transfer coefficient etc. The main result of the mathematical investigation is the mass loss curve for the wood sample as a function of reaction time and heating rate. The effect of heating rate is shown in Figs. 2 and 3. Heating rate affects mass loss curve position, as well as maximum decomposition rate. Fig. 2 shows the loss of mass as a function of reaction time at different heating rates for the particle thickness of 1 mm. It is observed that, as the heating rates increase, the pyrolysis process is faster, but when the process starts it ends very fast and the slope of the curve is similar in every case. Similar results were reported by other authors [24]. As the heating rates increase, the mass loss curve is just shifted towards lower temperature. According to the model, the duration of the pyrolysis process was 1450, 1050 and 680 s with heating rates values of 21, 32 and 55 K / min, respectively. As can be seen, char yields are lower at higher heating rates and the amount of char at the end of the process is about 15 % in every case. The differences in char yields between three cases were relatively small but noticeable. Fig. 3 shows the loss of mass with temperature at different heating rates (21 K / min, 32 K / min and 55 K / min) for the particle thickness of 1 mm. It is observed that the pyrolysis process starts at the same temperature (about 350 °C), and proceeds rapidly with increasing temperature. But, when the heating rate increases, the final temperature of pyrolysis process also increases. Higher heating rate has a short reaction time and the temperature needed for the sample to decompose is also higher.
Figs. 4, 5 and 6 show the temperature profile as a function of particle distance from the axis at various times of progression of pyrolysis (5, 10, 15, 20 and 25 min). It is observed that as the pyrolysis time increases, the temperature increases at all positions. The rate of increase is higher at the end of the process compared to that at the beginning. In the initial stages of pyrolysis, the heating  Case II: heating rate 32 K / min process is very fast, and as the time progresses, the heating process decreases. This can be explained by the fact that the heat transfer by convection and radiation from the gas phase to the wall surface at the initial stages of pyrolysis is very high. On the other hand, when the particle temperature is close to gas temperature, heat transfer to the wall surface is much smaller. As expected, the results show that as the heating rate increases, the time required for temperature rise decreases.
It is also observed that the increase in temperature at various positions at different times of progression of the pyrolysis is not significant, for the particle dimension of 1 mm, i.e. the temperature profile is flat, which is obviously due to the reasons mentioned above [20].
The above results shown in Figs. 4-6 and the justifiable and logical explanation given there have a lot of practical importance and physical significance in industrial pyrolysis applications. The results obtained consolidate the fact that it is possible to get the same extent of conversion of biomass and with less pyrolysis time under controlled conditions by increasing the heating rate. Fig. 7 shows the temperature profile as a function of time at the center (i.e. x = 0) at different heating rates (21 K / min, 32 K / min and 55 K / min) for a slab particle of thickness 1 mm. As can be seen, the temperature linearly increases with time and when it reaches about 723 K it begins to rise faster. That is the temperature when rapid degradation begins, Fig. 3. But, when the heating rate increases, the time to reach the critical temperature decreases. Higher heating rate has a short heating time.
Mass loss of wood sample as a function of reaction time and position at different heating rates are given in Appendix B.
The results obtained through model investigation of pyrolysis process of wood particles are tested by comparing them to the results obtained by experimental investigation.

Experimental analysis
To compare mathematical and experimental results, woody biomass was used as experimental sample. The sample was a mixture of six different kinds of woody biomass. Proximate and ultimate analysis of woody biomass is shown in Table 1. The laboratory facility is designed and constructed in the Institute for Energy and Process Engineering and Environmental Engineering at the Faculty of Technical Sciences in Novi Sad [25].   During the experimental investigation of the pyrolysis process, initial sample mass of the wood chips was 10 g. The mixture samples were placed into the biomass sample container and then into the reactor. The heating process was performed by using electrical heaters and after achieving the temperature of 923 K, that temperature was maintained within a narrow range around 923 K until the sample mass in the reactor was stabilized. Gaseous pyrolysis products were emitted to the atmosphere through the cooler.
The process of biomass mixture pyrolysis was carried out with controlling the temperature and mass of the sample. The samples were subjected to the temperature range of 300-923 K, and the average heating rates were 21, 32 and 55 K / min. All the experiments were conducted at atmospheric pressure.
The loss of mass profiles as a function of reaction time for quantitative comparison of the present model with experimental results is demonstrated in Figs. 8, 9 and 10. As can be seen, the mass change profiles obtained in the present study are in agreement with the experimental data. The difference (between the experimental and model results) is slightly more prominent with decreasing heating rate, but model results are in much better agreement with the experimental data in Case III.
According to the experimental data, the char yield is lower at higher heating rates. The higher char yields at lower heating rates are the consequence of slow reactions that occur in the reactor. Duration of the pyrolytic process of woody biomass mixture was 14, 18 and 30 min with heating rates values of 55, 32 and 21 K / min, respectively. After that period, mass ratio m/mo stabilizes. Mass ratio of char and biomass mixture sample was in the range 0.15-0.23 % at pyrolytic temperature around 923 K, in every case.
A simple model, such as the one developed in the present study with restrictive assumptions, combined with the thermal properties variation with temperature, can describe the overall progress of a set of processes of great complexity such as pyrolysis.
With some modifications the model developed here can be utilized to predict the temperature and concentration profiles for different types of biomass for a wide range of particle dimensions and temperatures.

Conclusions
In practice, pyrolysis of wood very often occurs in the form of wood chips (or some other form of wood processing by-products). That is why most of the kinetic models published in the literature were based on experimental tests on very small samples with a weight of several mg.
In this work, a mathematical model for single-particle processes has been set and numerically solved to predict   In particular, the solutions compared the effects of heating rate. Fig. 2 shows that with the lower heating rate the same size of biomass sample takes more time to pyrolyze than in the case with the higher values of heating rate. It clearly shows that lower heating rate pyrolysis produces more char and higher heating rate of pyrolysis produces more amounts of volatiles and gases and smaller amount of char.
The theory has been compared to the experiment. Because wood is a heterogeneous material, the mass loss curves show the existence of several zones associated with the devolatilization of the different components, as shown in Figs. 8, 9 and 10. Following this experimental observation, two-step reaction schemes of primary wood degradation have been proposed to describe the evolution of the reaction zones. The pyrolysis rate has been simulated by two parallel primary reactions and a third secondary reaction between the volatile and gaseous products and char.
Unlike the experimental results, the mathematical results show that once it reached the temperature of about 300-350 °C, pyrolysis process starts very fast and does not depend much on the temperature change.
This simplified model (where thermal degradation depends only on temperature, but not on the heating rate) could be further improved in order to obtain better results.
Apart from this difference which relates to the slope of the curve, the simulated results obtained from the model developed in the present study are in agreement with the experimental data (e.g. the temperature at the beginning and the end of the process, the amount of the char).
The main conclusion resulting from this study concerns the possibility of modeling the pyrolysis of a single biomass particle by coupling the heat transfer equation with the pyrolysis chemical kinetics equation.
As the model is simple and accurate, it is very useful in the design of industrial pyrolysis units. C B concentration of B, kg / m 3 C C 1 concentration of C 1 , kg / m 3 C C 2 concentration of C 2 , kg / m 3 C G 1 concentration of G 1 , kg / m 3 C G 2 concentration of G 2 , kg / m 3 c p specific heat, J / kg K h convective heat transfer coefficient, W / m 2 K k thermal conductivity, W / m K t time, s ρ 0 density, kg / m 3 σ Stefan-Boltzmann constant k 1 , k 2 , k 3 rate constants, s −1 A 1 , A 2 , A 3 constants, s −1 L 1 , L 2 constants, K 2 D 1 , D 2 constants, K n 1 , n 2 , n 3 orders of reactions T 0 initial temperature, K T temperature, K T f final temperature, K ΔH heat of reaction, kJ / kg ε emissivity coeficient