Article info

Volume         6    

Pages             92 - 98

DOI                 10.5027/jnrd.v6i0.10

Published     31/03/2017

Keywords     Biomass pyrolysis, Mathematical modeling, Approximation technique, Distribution function

Download (PDF, 957KB)


Alok Dhaundiyal a*and Suraj B. Singh b

a Department of Mechanical Engineering, Himgiri Zee University, Dehradun, U.K, India

b Department of Mathematics, Statistics and Computer Science, Govind Ballabh Pant University of Agriculture and Technology, Pantnagar, U.S. Nagar, U.K, India

* Corresponding author:



This article focuses on the influence of relevant parameters of biomass pyrolysis on the numerical solution of the isothermal nth-order distributed activation energy model (DAEM) using the Rayleigh distribution as the initial distribution function F(E) of the activation energies. In this study, the integral upper limit, the frequency factor, the reaction order and the scale parameters are investigated. This paper also derived the asymptotic approximation for the DAEM. The influence of these parameters is used to calculate the kinetic parameters of the isothermal nth-order DAEM with the help of thermo-analytical results of TGA/DTG analysis.


1. Introduction

Describing biomass decomposition mathematically is a complex task, as several decomposition reactions occur and their mechanisms are not familiar. Numerous mathematical approaches are considered to demonstrate the process of decomposition. One of them is the iso-conversion model, which postulates that kinetic parameters, such as the pre-exponential function and activation energy, are variable during the decomposition process and depend on the conversion rate. Another model is the lumped kinetic model that assumes large numbers of parallel decomposition nth-order reactions are taking place. These partial reactions contribute to the overall decomposition step. The most accurate and well defined model used for modeling the biomass pyrolysis is the distributed activation energy model (DAEM). This model assumes that several decomposition reactions with distributed activation energies are taking place simultaneously. The principle of the lumped kinetic model is very similar to that of the DAEM with the main difference being the number of expected decomposition reactions, though the flexibility of the DAEM makes it more versatile. Application of the DAEM is not only confined to biodegradable material and its main components, as it can also be used to describe other thermally decomposable materials, such as coal [1], [2], medical waste [3], and sewage sludge [4], [5].

The distributed activation energy model, or multiple reaction model (MRM) is applicable to either the total amount of volatiles released or the amount of an individual volatile constituent, such as carbon monoxide, or tar [6]. Numerical solutions of the DAEM involve several evaluations of double integral terms, posing significant numerical difficulties. To tackle the numerical complication, asymptotic expansion was adopted for accurate approximation to the integral term. The practical scope of asymptotic approximations falls within the area of numerical sciences, such as Computational Fluid Dynamics (C.F.D) modeling of coal-fired boilers, where it is very significant to evaluate the double integral quickly. The main aim of involving asymptotic expansion is to focus on the analytical insight from the solutions, since the asymptotic forms are explicit in form.

The purpose of this study is to find an accurate approximation of the DAEM with the help of the asymptotic technique and then assess the influence of relevant parameters of the biomass pyrolysis on the numerical solution. The nth-order DAEM equation is expressed by equation (1):

where F(E) is the initial distribution function of the activation energies, n is the reaction order, E is the activation energy, ko is the frequency factor, R is the universal gas constant and t is time.

Though in most of the studies, the distribution function is assumed to be symmetric in nature, Gaussian [7], it may sometimes be advantageous to select an asymmetric distribution for modeling the kinetics of biomass pyrolysis, such as the Rayleigh, Gamma, Erlang and Weibull distributions, as in practice, thermo-analytical data has an asymmetrical distribution. Data can either be positively skewed or negatively skewed. Therefore, F(E) was chosen to be the Rayleigh distribution in the present study. In addition, to date residual waste of cedrus deodara is not used for pyrolysis, though in the future government agencies may use it as a new alternative fuel.


then we have,

The relation between the scale parameter and the mean and variance is stated as:

2. Asymptotic methodology for isothermal pyrolysis

The double integral term and double exponential term in equation (1) are difficult to compute, as they needs lengthy computing resources. Therefore, the asymptotic technique is applied to solve the rapidly varying functions. In the previous approximation to the solution of equation (1), the drawback is the need to extrapolate to other heating regimes, thus complicating the process of obtaining the volatile distribution F(E) from the thermo-analytical data. Niksa and Lau derived the relationship between the DAEM and the single first-order reaction (SFOR) in their paper, based on keeping the activation energy fixed and hence defining an effective rate constant [8]. They also provided an analytical approximation to the DAEM for the exponentially and linearly varying profile of temperature. In order to demonstrate the approach, the double exponential term (DExp) has been considered as a piece-wise linear function. The unit-step function is incorporated to ensure an accurate approximation [6], [9]-[11]. In the subsequent subsections a systematic simplification of the double exponential term is discussed.

2.1. Approximation to the Double Exponential term (DExp)

Equation (1) comprises two terms. The first term, DExp depends on time through the temperature encountered by the sample, whereas the second term solely depends on the distribution of volatiles in the sample. The characteristics of DExp are discussed firstly and thereafter approximations are derived for the relevant problems.

The double exponential term can be expressed as:

Putting in the value of T(l), DExp becomes

Assumptions of typical values of parameters and functions on which it depends are considered to demonstrate the stepwise simplification of DExp. The frequency factor varies from k0~1010 - 1013 s-1; whereas the activation energies are in the range of 100-300 kJ/mol.

To keep the method simple, the temperature is assumed to follow an isothermal profile. According to equation (2), the double exponential term can then be rewritten as:

As the behavior of s(E) in the neighborhood of Es is of interest, it is expanded with the help of the Taylor series.


If the value of the activation energy E is less than the stationary central value Es, DExp approaches zero. On the other hand, for all E>Es the value of the function is approximately one. The variation of the function in the domain of [0, 1] in a range of E values within the step width of Ew of Es.

In the following section, we will discuss approximation of initial distribution function and one of the special cases of the distribution type.

2.2 Initial Distribution function (F(E))

The two different limits are discussed under this heading. The first is when the initial distribution function F(E) is relatively wide as compared with the width of DExp while in the second case, the width of F(E) is relatively narrow compared to the width of DExp. The significance of distribution type is seen in the fact that the shape of the whole integrand as stated in equation (1) varies with the width of the functions. The shape of the curve adheres to the distribution function F(E) during the initial phase of pyrolysis, but it is progressively chopped off from the left by the step-like DExp as time proceeds. In case of a narrow distribution, the total integrand will be more symmetrical than that of the wide distribution; however, the location of its maximum is not stationary. The scope of paper is limited to the wide distribution case, and therefore the limit applied to the isothermal problem is the wide distribution. The Rayleigh distribution was chosen as the initial distribution function, centered at E0 with standard deviation σ used to find the approximations of the DAEM.

Approximations of equation (1) can be taken as:



For constant ramp temperature T=T0,we have

2.2.1. Wide distribution case

In order to demonstrate the wide distribution type, the limit σyw <<1 is considered. After implementing this limit, DExp will vary from zero to one in the neighborhood of y= ys, and is thereafter approximated with the help of the step function H (y-ys ) as follows:

Equation (3) is modified in the form:


From equation (4), it can be seen that the first integral is multiplied by a function that is very small everywhere but near the point y=ys. Thus, the integrand is approximated with the help of a Taylor series expansion.

After the Taylor series expansion, the terms can be separately integrated to get the results for first order reaction (n=1)

Coefficients Mn (n=0,1,2,3) need to be evaluated once, as they are invariable to any parameters. The approximated values of the coefficients are:

The remaining coefficient values are equivalent to the following integral:

Furthermore, to carry out the simplification for the nth-order reaction, we have

These terms arose while solving equation (6) and can be integrated exactly in the same manner as performed with equation (4) to obtain the desired results:

The initial values of Pn, Qn, Rn coefficients are:

The values of the remaining coefficients are evaluated by the integrals as follows:

2.3 Application of Biomass

The TGA/DTG analysis of cedrus deodara was done in the presence of an inert atmosphere. The chemical characteristics of the sample were identified with the help of a CHNS (O) elemental analyzer (Euro-E3000). It is to be noted that the results of this paper was used in the process of simulating the experimental data. It can be seen that the nth-order Rayleigh DAEM shows good agreement with the thermo-analytical data. Table 1 shows the chemical composition of the cedrus deodara. Moreover, a Matlab algorithm was applied to solve equations (5) and (7) to minimize the root mean square error and obtain the simulated results for the nth-order Rayleigh DAEM.


3. Results and Discussions

In order to obtain the numerical solution of the nth-order Rayleigh DAEM, asymptotic expansion of the kinetics equations was performed. For the solution of equation (1), the upper limit of dE must be determined. The influence of the upper limit (E) on the numerical results of the isothermal nth-order DAEM equation is shown in Figure 1. At the initial stage of the pyrolysis reaction, the remaining mass fraction v must be close to 1. In Figure 1 however, it is seen that the remaining mass fraction v is less than 1 for 76.06 kJ/mol <E values. When E values less than 69 kJ/mol are used, the results are more accurate and closely proximate each other. Therefore, 66 kJ/mol can be used as the upper limit for the activation energy. The effect of the frequency factor (A) on the numerical solution is shown in Figure 2. With an increase in the value of A, the v curves shift to the left. The effect of the scale parameter (β) on the numerical solution is depicted in Figure 3. Increases in β values make the curves shift up and therefore the remaining mass fraction lines tend to be constant with respect to time in the case of isothermal pyrolysis. The effect of the reaction order (n) values on the numerical results is illustrated in Figure 4. It can be seen that the increase in n values causes the v curves to lead toward the right. The influence of varying these parameters on the final simulated results does not provide as good a fit as in the problem of non-isothermal pyrolysis. Comparative variation of simulated results, which is obtained with the help of the best suited unique values of the parameters, is shown in Figure 5.


4. Conclusions

It can be concluded from the results that varying the scale parameter merely changes the attribute of the remaining mass proportion curves. In the numerical solution of the isothermal nth-order DAEM using the Rayleigh distribution, 66 kJ/mol can be used as the upper limit for the activation energies. The results are helpful for determining the kinetic parameters of the isothermal nth-order Rayleigh DEAM from the TGA/DTG data on biomass pyrolysis. Besides knowing the influence the kinetic parameters have on the numerical solution, the behavior of the Rayleigh distribution towards the isothermal pyrolysis condition is not as flexible as it is with other distribution functions, such as the Gaussian distribution. The approximated results in the case of the Rayleigh distribution depend on other more decisive variables. In case of the Rayleigh distribution, asymptotic techniques do not provide relatively more accurate curve fitting than non-isothermal conditions [7], [12]. The variation in the pattern of the curves due to some relevant parameter of biomass pyrolysis is negligible. Mathematical rigidity of the Rayleigh distribution towards the isothermal pyrolysis of biomass is found in this study. However, the multivariate distribution function may provide more promising outcomes for the future scope of the Rayleigh distribution in modeling biomass pyrolysis.


5. References

[1] C. Quan, A. Li, and N. Gao, “Thermogravimetric analysis and kinetic study on large particles of printed circuit board wastes,” Waste Management, vol. 29, no. 8, pp. 2353-2360. 2009. Doi:

[2] J. H. Yan, H. M. Zhu, X. G. Jiang, Y. Chi, K. F. Cen, “Analysis of volatile species kinetics during typical medical waste materials pyrolysis using a distributed activation energy model,” Journal of Hazardous Materials, vol. 162 no. 2, pp. 646-651. 2009 Doi:

[3] H. M. Zhu, J. H. Yan, X. G. Jiang, Y. E. Lai, K. F. Cen, “Study on pyrolysis of typical waste material by using TG-FTIR analysis,” Journal of Hazardous Materials vol. 153, pp. 670-676. 2008. Doi:

[4] M. B. Folgueras, R. M. Díaz, J. Xiberta, I. Prieto, “Thermogravimetric analysis of the co-combustion of coal and sewage sludge,” Fuel, vol. 82, pp. 2051-2055. 2003. Doi:

[5] M. Otero, L. F. Calvo, M. V. Gil, A. I. García, A. Morán, “Co-combustion of different sewage sludge and coal: a non-isothermal thermogravimetric kinetic analysis,” Bioresource Technology, vol, 99, no. 14, pp. 6311-6319. 2008 Doi:

[6] J. B. Howard, “Fundamentals of coal pyrolysis and hydropyrolysis,” in Chemistry of Coal Utilization, Edited by M.A. Elliott: Wiley and Sons, 1981, ch. 12, pp. 359-386.

[7] A. Dhaundiyal and S. B. Singh, “Analysis of the non Isothermal Distribution Activation Energy Model for Biomass Pyrolysis by Fuzzy Gaussian Distribution,” Rural Sustainability Research, vol. 35, no. 330, 2016. Doi:

[8] S. Niksa and C. W. Lau, “Global rates of devolatilization for various coal types,” Combust. Flame, vol. 94, pp. 293-307. 1993.

[9] A. Vand, “A theory of the irreversible electrical resistance changes of metallic films evaporated in vacuum,” Proc. Phys. Soc. Lond. A, vol. 55, no.3, pp. 222-246. 1943.

[10] G. J. Pitt, “The kinetics of the evolution of volatile products from coal,”. Fuel vol. 1, pp. 267-274. 1962.

[11] E. M. Suuberg, “Approximate solution technique for nonisothermal, Gaussian distributed activation energy models,” Combust. Flame, vol. 50, pp. 243-245. 1983.

[12] A. Dhaundiyal and S. B. Singh “Distributed activation energy modelling for pyrolysis of forest waste using Gaussian distribution,” Proceedings of the Latvian Academy of Sciences Section B Natural Exact and Applied Sciences, vol. 70, no. 2, pp.64-70. 2016. Doi:

2351 View