Introduction
There are various models which are proposed to describe the pyrolysis of biomass, such as single-reaction and multi-reaction models [1-6]. The most authentic and accurate model considered in this study is the Distributed Activation energy model (DAEM)[7-10]. While the aim of this study is mainly concentrated around the parametric values relevant to the loose biomass, the DAEM also applies to the pyrolysis of other conventional sources including residual oils, resin chars[11], and kerogen[12]. Calculations of the solution to this model may require iterative loops of double integral and rapidly varying functions which in turn create significant numerical complications. To tackle this problem, some asymptotic methods have been adopted to make an accurate approximation to the integral and therefore curtail the computation time. Another main source of numerical difficulty is due to the double exponential term (DExp) in the DAEM, which is investigated for isothermal pyrolysis. The domain of DExp acts over a narrow range of activation energies which changes as time progresses. In this study, the key to our approach is to know the significance of the relative width of the DExp term as compared with the width of the initial distribution.
In this study, the numerical solutions have been obtained by using the asymptotic expansions. The derived results are used to determine the kinetic parameters of the kinetic model. For predicting the realistic results, the parameters affecting the behavior solution must be estimated. The effect of these parameters on the single-reaction model is reviewed in literature[11]. The DAEM is very flexible, and it can successfully describe the pyrolysis of different types of biomass. The main aim of this study is to focus on the relevant parameters which affect the kinetics of pyrolysis. The DAEM is also applicable to the pyrolysis of other conventional sources of energy as Coal, residual oils, resin chars[12], and kerogen[13]. This analytical method is not only used for thermal decomposition of plant or animal biomass[14,15] but also for other materials such as medical wastes[16] waste car tyres[17], printed circuit board wastes[18], or sewage sludge[19,20].
The purview of this paper is to use asymptotic techniques to make an accurate approximation to the double as well as DExp terms and then predict the behavior of the isothermal nth order DAEM by involving the effect of some relevant parameters on the numerical solution.
Materials and methods
Distributed Activation Energy Modeling:
This model includes a reaction time scale, which gained acceptance since it is the significant part of biomass devolatilization[21,22\. The complication related to DAEM is that the function f (E) and k0 (E) are highly correlated, hence it is very difficult to evaluate both the functions accurately. Therefore, it has been assumed that all the frequency factors k0i to have the same value k0, thus it makes analysis easy. The uncertainty of reactant distribution is highly emphasized. The nonisothermal nth DAEM is expressed by equation (1)
where X is conversion, n is the order of reaction and f(E) is the Gamma distribution function of activation energies.
The value of conversion is found with the help of TGA analysis of cedrus deodara.
where, m r is the residual mass, m0 is the mass of the sample at beginning of decomposition, and mt is the mass of sample at given time.
Although in most of the study, the symmetrical function is assumed (Gaussian), yet it would be advantageous to choose an asymmetric distribution for modelling the kinetics of biomass pyrolysis, such as the Gamma distribution over a symmetrical one[23]. In addition to that the chosen function is mathematically flexible and can be expressed as:
where λ is the scale parameter expressed in kJ/mol and n is dimensionless positive shape parameter. The mean and the variance of distribution are given by equation (4) and equation (5) respectively:
Approximation methodology
Equation (1) comprises of two terms. The first term is the DExp which varies with time through temperature history experienced experimentally. The second part is the initial distribution function of activation energies that is independent of time. Primarily, the temperature dependant part is considered and thereafter one can derive approximations to solve the DExp.
The isothermal temperature profile only has been involved in this study. In the nonisothermal, linear ramping temperature, the amount of released volatile does not change until the critical time is attained when the two parts of integrand overlap significantly, whereas the variation of volatile release begins to change appreciably at very early times. In case of ramping temperature, the mean value of DExp changes its location with time in a similar manner to the isothermal temperature case, with logarithmic term replaced by a LambertW function[24,25]. However, the step width of DExp, Ew (= E0 y w), is narrow at early times, which in turn leads to a major difference in the appearance of the remaining mass fraction curves versus time, between the ramping and constant temperature cases.
Systematic approximations to distributed activation energy model (DAEM)
For approximation of double exponential term, the primal step is to assume the typical values of dependent parameters and rapidly varying functions. The frequency factors (k0) are in range of k0 ~ 1010 - 1013 s-1, whereas the activation lies in domain of 100 - 300 kJ/mol[24[. The double exponential term is expressed as:
To demonstrate the simplification method, the constant temperature history is assumed to be as follows:
Here, l is any instant of time t.
After implementing the isothermal condition, DExp becomes
Assume typical values, , while tk0~1010. The large size of both these parameters make function, DExp varies rapidly with E. The equation (8) can be further illustrated as:
Where Es - RT0 ln(tk0) and Ew - RT0
When the value of E is much less than Es, the function approaches zero. Whereas for E much greater than Es, DExp tends to one. The DExp varies in domain of (0, 1) in a range of E with step width of Ew. The distribution can either be wide or narrow distribution[24]. It depends on the width of DExp as compared to the width of the initial distribution function f (E). In the total integrand of equation (1), DExp is multiplied by the initial distribution function f (E). The shape of the total integrand relies upon the applied limit. When the initial distribution is relatively wide compared to Ew, the total integrand is initially following the shape of the initial distribution function, but as time proceeds, it is progressively truncated from the left by the step function, DExp. The location of maximum can move significantly and the shape becomes quite skewed. Conversely, the relatively narrow width of initial as compare to DExp makes the shape of the total integrand similar to the initial distribution, with amplitude that is progressively reduced by DExp as time proceeds. The total integrand remains more symmetrical than that of wide distribution. The location of maximum of total integrand does move in the time-dependant manner.
But the scope of this paper is confined itself to the wide distribution case, wherein the width of f (E) is wider than that of DExp.
Using equation (1) and (3), the remaining mass fraction can be expressed as:
where Es and Ew are functions of t as mentioned earlier.
Energy is now rescaled as , so that the problem becomes
let, . For the practical purpose α << 1.
Time is also rescaled as T = k0t
For linear ramping temperature T = θt,
The wide distribution case
It is assumed that the initial distribution is wider than DExp. According to this, DExp jumps from zero to one near y = ys in such a way that it has been approximated by the step function[21, 22, 26, 27]. In order to adopt the wide distribution method, the limit is considered. The Heaviside function is written as
Equation (13) can be written in the form:
where is the upper incomplete Gamma function.
The second term in equation (13) is complementary function, and hence easily computed. In fact, in many previous simplification (the step-function approximations) include just this term and neglect the remaining terms. The first integral term is negligibly small everywhere except in a neighbourhood of size yw about the point y = ys . Therefore, it can be expanded with the help of Taylor expansion about y = ys .
Suppose
Putting , d y = yw dx in equation (17), we have
From equations (17) and (18), we have
or
We know that
where is the lower cumulative distribution.
Coefficient Lt is independent of any parameters and the some few values are evaluated as:
The remaining integral terms are estimated by the expression
The equation (20) is the required expression for the first order reaction.
In the same manner, the approximations can be obtained for nth order reactions by invoking equation (1).
The equation (1) for nth order reactions is written as:
Using the wide distribution limit, the equation (23) can be expressed as:
Or
The coefficients Pn, Mn and Nn values are evaluated as:
The other integral terms can be evaluated as:
Application of loose biomass and computation methodology
The sample of Cedrus deodar underwent isothermal pyrolysis. Elemental composition is evaluated with the help of CHNS (O) analyser (Flash-EA 1112 series). TGA/DTA (Exstar 6300) analysis has been performed in the presence of inert atmosphere of Nitrogen. An alumina crucible is used to hold the sample. The purge flow rate is fixed to be 200 mL/min. Dulong’s formula for estimating the calorific value of a solid fuel is used to evaluate the high calorific value of sample (28). Thermocouple ‘R’ type is used to measure the furnace temperature.
It is to be noted that the results of this analysis are used for the prediction of nth order DAEM using the Gamma distribution. Computation of equations (20) and (25) are done with the help of MATLAB algorithm. Accurately approximated results are obtained by minimizing the root mean square error between experimental and predicted value. Iterative loops are used for lengthy computation of DAEM equations. Fig. 6. demonstrates that the nth order Gamma DAEM provides the good fit with experimental data.
Result and discussions
The numerical solution of equation (1) is carried out with the help of asymptotic expansion. The parametric influence of upper limit (E∞) is illustrated in Fig. 1. At the initial stage of the pyrolysis, the remaining mass proportion (1-X) must be close to one. While it has been observed in Figure 1 that the remaining mass fraction is less than one for E∞< 8.33 kJmol-1 values. When more than 11.33 kJ is applied, the results are more accurate and closely overlapped with each other. Therefore, 13.8 kJ mol-1 can be adopted as an upper limit of the ‘dE’ integral. The behaviour of (1-X) curves with variation of frequency factor (A) on the numerical results is depicted in Fig. 2. According to these curves, increment in A shifts the curves to the left direction.
The effect of scale parameters on the numerical solution is shown in Fig. 3. Where it is visible that the remaining mass proportion curves shifted up for the lower value of scale parameter (λ). Therefore, the initial value of (1-X) becomes more than one, which provides inaccuracy of the predictive algorithm for the values of λ less than equal to 35.85 kJ mol-1. The effect of shape parameter (ŋ) on the numerical solution is depicted in Fig. 4. For the values of ŋ less than 13.36, the remaining mass curves shifted down and become constant with time.
Thus, it is concluded that the converse rate of biomass becomes constant and hence, the remaining mass fraction curves show the asymptotic behavior with time as the shape parameter decreases. The effect of the reaction order (n) values on the numerical results is illustrated in Fig. 5. from which it is seen that increase in reaction order (n) causes (1-X) curves to shift up.
Unlike Gaussian[29], Weibull[30] and Rayleigh[31] distribution functions, Gamma distribution relatively converges at very low value of activation energies. It also implies that sensitivity of Gamma distribution while modeling the biomass pyrolysis is quite high. The sudden variation in relevant parameters of biomass pyrolysis led to the drastic variation in the curvature of the remaining mass fraction. Therefore, the variation of activation energies in between different reactions is negligibly small.
Conclusions
Numerical solution of Gamma distribution is derived by using asymptotic approximations to double exponential term. It has been found that the upper limit of ‘dE’ provided a good curve fitting at 13.8 kJmol-1; whereas the remaining mass fraction curves did not predict an accurate result for all the value of E∞ ≤ 11 kJmol-1. Furthermore, the effect of frequency factor, reaction order, and shape and scale parameters of Gamma distribution merely influenced attribute of (1-X ) curves. However, the remaining mass fraction curves did not simulate the mechanism of the reaction occurring at the high-temperature regime. If we replace the initial distribution function from some other multivariate function which bifurcates the behavior of biomass pyrolysis at lower as well as higher temperature regimes, it may provide a much more promising outcome than that of univariant function.