Application of Monte Carlo Simulation for Energy Modeling | ACS

Evaluation of kinetic parameters and their effect on the biomass pyrolysis is ... of chir pine waste are obtained at different heating rates of 5, 10,...
0 downloads 0 Views 3MB Size
Download PDF
This is an open access article published under an ACS AuthorChoice License, which permits copying and redistribution of the article or any adaptations for non-commercial purposes.

Article Cite This: ACS Omega 2019, 4, 4984−4990

http://pubs.acs.org/journal/acsodf

Application of Monte Carlo Simulation for Energy Modeling Alok Dhaundiyal,*,†,‡ Suraj B. Singh,‡ Divine Atsu,†,§ and Rashmi Dhaundiyal∥ †

Doctoral School of Mechanical Engineering, Szent István University, Gödöllő 2100, Hungary Department of Mathematics, Statistics & Computer Science, Govind Ballabh Pant University of Agriculture and Technology, Pantnagar 263145, U.K., India § Department of Energy Systems Engineering, Koforidua Technical University, Koforidua KF-981, Ghana ∥ Department of Computer Science and Engineering, DIT University, Dehradun 248009, U.K., India ‡

Downloaded via 91.243.190.12 on March 14, 2019 at 01:13:02 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: This work is pivoted on the implementation of Monte Carlo simulation technique in the thermal decomposition of biomass. Evaluation of kinetic parameters and their effect on the biomass pyrolysis is investigated through proposed integral scheme. Activation energy modeling relies upon the pairing of random points generated through the MATLAB program. The numerical solution of the multireaction model is based on the ratio of numbers of pair points (xi, yi) and the points (i) generated. Thermoanalytical data of chir pine waste are obtained at different heating rates of 5, 10, and 15 °C min−1. The scheme involves the inverse formulation of kinetic model; thus, the kinetic parameters can be more precisely computed than the hit-and trail-methodologies. Rather than speculating the behavior through regression analysis, the experimental thermogravimetric plane provides more correct information. So, the conversion of biomass is converted into the area of curve bounded by the experimental data. For experimental purpose, nitrogen is employed as the purge gas to sterilize the furnace at the volumetric rate of 200 mL min−1. The relevant regime of temperature follows a non-isothermal profile for temperature range of 300−900 K.

■

stochastic scheme to solve a multireaction model. Cai and Lui16 implemented the basic numerical techniques to solve the model by approximating the Arrhenius equation, but the error component while allocating the steps is not perused in his studies. Later, pyrolysis proposed by Dhaundiyal et al.17−19 must be considered as a series of independent events and could be represented by the distribution function that imbibes the heterogeneous attribute of the marginal density functions. They used the Archimedean family of copula to solve the multireaction model. Previously, an attempt of replacing the basic form of the Gaussian distribution with the convex combination of probability functions is done by some researchers.20 But one must know the thermoanalytical data of biomass with respect to time or temperature continuously changes with time, as the activation energy is highly overlapped due to the large number of parallel reactions. Moreover, the objective of the analysis must not only be pivoted on the basic fitting of the thermoanalytical data to model but also comply with the intrinsic characteristic of biomass pyrolysis. The objective of this work is to develop a new method to collect useful information about the kinetic properties of the

INTRODUCTION Pyrolysis of biomass is one of thermal conversion techniques to extract the secondary fuels by the decomposition of biowaste in the absence of air. But the yield of the secondary fuels depends on the kinetic as well as the thermal mechanism of the process, albeit it is not necessary that thermally driven reactions may be supported by the kinetic conditions. These two aspects are largely guided by postulation or the well-defined models that help to evaluate the basic parameters required to trigger any chemical reactions. The quantitative aspect of the biomass pyrolysis modeling is categorized into the single-step and the multistep reactions.1−6 These models are comprehensively reviewed by Burnham et al.,7 who suggested that the discretize energy model could provide the reasonable kinetic properties. Moreover, the logistical distributions do not follow the same principle as the reactivity distribution.7 The Gaussian distribution model for the constant heating rate was initially proven to be influential step in the pyrolysis community7−9 despite its shortcoming in evaluating kinetic properties.7 Pitt10 and Hanbaba11 adopted the parallel reaction approach, which was later recognized by the work of Howard9 and other researchers.12−14 Dhaundiyal et al.15 reported that the Rayleigh distribution model provided a better insight into demarcate stages of thermal decomposition of biomass than the Gaussian model. However, there is a scope for qualitative development of © 2019 American Chemical Society

Received: December 23, 2018 Accepted: February 18, 2019 Published: March 7, 2019 4984

DOI: 10.1021/acsomega.8b03442 ACS Omega 2019, 4, 4984−4990

ACS Omega

Article

pyrolysis. The logical matrices are used to identify the target area of particle. Once the particle passes the given criteria, it is separated from the remaining particles. At end of computational cycle, the number of particles bombarded on the specified area are calculated by summing the arrays that qualify the prerequisite condition. The similar concept is used for all thermograms generated at different heating rates. Bombardment of 17 to 74 million particles are performed by using MATLAB coding. All the inbound particles are divided by total number of particles dropped on the plane of thermogravimetric curves. Thus, the effective fraction of successful particles is found to determine to the enclosed area.

thermal decomposition of biomass, as well as to implement similar approach to suit other thermal conversion processes. In order to make it versatile, it is necessary to develop the model that supports the multithermal profile, thus the non-isothermal profile is chosen. This study combines the asymptotic assumptions and the calculation of the actual values of kinetic parameters within the well-defined constraint. The proposed scheme works on the principle of the inverse mathematical problem. Thermogravimetric technique is used to measure the extensive property of biomass with respect to temperature. Moreover, it is assumed that the conversion function is the injective function of activation energy (E) (i.e., one-to-one mapping). Application of Biomass. The statistical modeling of DAEM used the thermoanalytical data of pine needles. The sample of pine-waste undergoes non-isothermal decomposition at the heating rate of 5, 10, and 15 °C min−1 in the presence of inert atmosphere of nitrogen (N2). Thermal evaluation of the pine sample is performed by thermogravimetric-differential thermal (Diamond TG/DTA, Perkin Elmer). The samples of similar weight, 10.58 mg, are kept inside an electrobalance for an hour duration. The temperature range of 300−982 K is used to decompose the pine-waste. To avoid the measurement error due to sudden increase in the pressure of inert gas, the horizontal differential type balance is chosen for the experimental task. Thermocouple type “R” is referred to measure the temperature inside the furnace. The volumetric flow of inert is fixed at 200 mL min−1. Indium and tin are considered as the reference materials. On the other hand, the chemical analysis is performed by CHNS-O analyzer (vario MACRO). Before bringing the analyzer into operation, the furnace is heated up to 1473 for 30 min. Once it reaches saturation temperature, the capsule form of samples is put inside the rotating turret. The flow of oxygen maintains the catalytic combustion, whereas helium gas is used as a carrier gas for product of combustion. Thus, the volatile gases with water (CO2, SO2, NO2, and H2O) can be absorbed at different reduction columns. These tubes are juxtaposed to combustion chamber and detection unit, where the gas mixture is separated into the constituent elements by purge or trap chromatography. Finally, each component is individually detected by thermal conductivity detector (TCD). The sequence of measurement is CO2, H2O, and SO2. However, nitrogen gets by-passed through all the three columns to be finally measured by TCD, which generates the pulse that is proportional to the concentration of elementary component of biomass. The bomb calorimeter at constant volume is implemented to measure the calorific value of pine needles. Elemental composition of pine waste is mentioned in Table 1.

■

RESULTS AND DISCUSSION Mathematical analysis of kinetic model is performed by using the Monte Carlo simulation scheme. The DAEM is transformed into the inverse form to obtain the kinetic parameters for thermal decomposition of pine waste. The estimated range of variation of activation energy is demarcated by the upper limit (EU) and the lower limit (EL). To carry out simulation process, random numbers are generated and scattered on the plane of thermogravimetric (TG) curves, which are experimentally obtained at different heating rates. The numbers of points are separated by the binary algorithm and the area under the target plane is estimated. The computed values of TG curves areas are provided at Table 2, whereas the graphical representation of the Brownian motion of random particles at different heating rates is illustrated in Figures 1, 2, and 3. The physical importance of energy distribution function is plotted in Figures 4 and 5. The estimated area of TG curve varies from 308 to 329 sq. unit. The variation of area happened, as the number of particles (inbound) bombarded on the target plane of TG curves vary relative to the constraints imposed on the total plane. From 14 to 74 million of particles are generated, but most of the particles are generated accordingly to the experimental matrix. The metric length of experimental matrix varies; therefore, the area under the curve also gets changed. It can be observed that some of particles are not able to pinpoint the target area properly. Some of them are sparsely located at the interface TG plots. The second reason of variation in area of curve is that the remaining mass fraction (1 − α) gets shifted to right as the heating rate changes, thus the area under curve becomes dynamic in nature. The experimental tests are performed to establish the trend among the Arrhenius parameters (E, A, and k) with the varying heating rate. Characteristic of energy distribution relative to the central value (Ew) is evaluated with the help of critical width (σc). The upper range of activation energy varies from E ∼ 51.54−55.07 kJ mol−1, whereas the lower range is almost constant, as the pyrolysis is divided into the three main regions (water evaporation, devolatilization, and char formation)38 and the activation energy remains constant during the initial stage of pyrolysis (dehydration of biomass). The statistical information of the energy distribution function derived from the Monte Carlo simulation is illustrated in Table 2. The width (σd) of the distribution energy function is comparatively narrow to the width (Ew) of the double exponential function (Dexp). The limiting condition of whether the energy distribution function is to be called narrow or wide distribution is tabulated in Table 3. It is graphically illustrated in Figure 4 that the metric space of energy distribution for thermal decomposition reaction is highly overlapped in the range of 30−40 kJ mol−1, whereas it is 32.5− 42.5 kJ mol−1 for the central value (Es). As the activation energy (E) goes higher than (Es), the rapidly varying function Dexp will

Table 1. Elemental Composition of Pine Needles C%

H%

N%

O%

S%

ash content (%)

53.7

6.48

0.56

36.13

0.10

2.92

higher heating value (MJ kg−1) 22.23

Computational Algorithm. The non-isothermal temperature profile is used to solve the DAEM. MATLAB 2015a is used for numerical purpose. The computational algorithm is based on the generated random number. Before generating the random points, the constraints are developed so that the generated particles do not violate the experimental conditions required for 4985

DOI: 10.1021/acsomega.8b03442 ACS Omega 2019, 4, 4984−4990

ACS Omega

Article

Table 2. Statistical Parameters Derived from Monte Carlo Simulation for DAEM E (kJ mol−1) −1

Es (kJ mol−1)

β (°C min )

area under TG curve (sq. unit)

EU

EL

Em

σ

5 10 15

329.3209 308.4438 325.3855

51.54 51.69 55.07

16.05 16.03 16.69

34.75 34.20 35.92

115.13 101.94 123.38

2

−1

Ew (kJ mol−1)

A (min )

Ems

σs2

3.9221 × 103 3.9549 × 103 5.1067 × 103

37.36 37.08 37.67

147.13 156 174.76

Emw

σw2

4.57 4.54 4.63

1.82 1.93 2.19

Figure 5. Normal distribution of energy spacing (Ew) obtained for the temperature integral.

Figure 1. Graphical representation of Brownian motion caused by the continuous bombardment on the thermogram obtained at 5 °C min−1.

Table 3. Limiting Value of Distribution Function Width,

1 λ

,

and Critical Value of the Standard Deviation (σc) To Determine the Nature of Energy Distribution β (°C min−1)

λ

5 10 15

5.24 5.73 5.22

σc =

RβtcE0 q(Rβtc + E0)

1.27 1.19 1.13

affect the conversion rate

Figure 2. Demonstration of the Monte Carlo simulation obtained at 10 °C min−1.

σd = 0.43 0.41 0.43

( ddαt ) of

1 λ

r2

tc (min)

q

0.95 0.97 0.96

42 21.38 17.41

20 17 17

chemical reactions. The

Figure 3. Graphical view of the numerical solution derived at 15 °C min−1.

numerical solution converges only when the activation energy remains near the central value (Es). The biomass decomposition reactions are sigmoidal reactions, which means the acceleratory phase of the chemical reaction is always accompanied by the deceleration stage of some kinetically stable products. It can be realized in two ways. If the activation energy goes beyond Es, the conversion (α) rate decreases and the reactants become more stable. Second, if the tangent on eq (A.13) decreases, the conversion rate will also decrease. The points of inflection on the TG curve is obtained at the activation energy range of 30−40 kJ mol−1. It implies all reactions beyond 40 kJ mol−1 that belong to the deceleration stage of biomass pyrolysis, whereas the reactions within the range of 16−30 kJ mol−1 are actively participating in the thermal-conversion process (Figure 6). The derived values of activation energies are tested by putting into

Figure 4. Physical significance of energies (E and Es) distribution obtained at different width and the mean values.

Figure 6. Variation of activation energy (E) as the function of conversion (α). 4986

DOI: 10.1021/acsomega.8b03442 ACS Omega 2019, 4, 4984−4990

ACS Omega

Article

Figure 7. Comparison between experimental data and the predicted solution obtained from Monte Carlo simulation (a, c, e: the experimental curves obtained at 5, 10, and 15 °C min−1, respectively; b, d, f: the residual plots obtained at 5, 10, and 15 °C min−1, respectively).

Table 4. Validation of Results with the Other Statistical Method of Computing Kinetic Models E (kJ mol−1) Monte Carlo simulation (51.54) ≤ E ≤ 55.07)

Copula based

18

41.5 ≤ E ≤ 84

model-free methods39−41,43 67.63 ≤ E ≤ 89.59 (cedrus deodar) 49.5 ≤ E ≤ 109.5 (Botryococcus braunii)40 153.92 ≤ E ≤ 157.27 (Poplar wood)41 39

model-fitting method42

the Rayleigh model17

43.7 ≤ E ≤ 197.31

39.4 ≤ E < 84.04

■

CONCLUSIONS The parametric information about the kinetics of thermal conversion process is retrieved through the implemention of the Monte Carlo simulation. The upper bound of the activation energy varies from 51 to 55 kJ mol−1. The frequency factor changes from 3.921 × 103 to 5.10 × 103 min−1. It is observed that the heating rate changes the characteristic of distribution function. Moreover, the predicted solution is also affected by relative variation of width of double exponential and distribution function. The temperature integral or Dexp is relatively wide to the distribution function of activation energy. The critical value of width is determined by comparing the energy spacing width Ew with qσc . The Monte Carlo simulation process provided the

the given DAEM and it is observed that the coefficient of regression (r2) varies from 0.96 to 0.97 (Table 3). Comparison between the predicted solution and the experimental curves is shown in Figure 7. It is also clear from the mathematical analysis that the increasing heating increases the energy spacing among chemical reactions (Figure 5); therefore, the rate of conversion increases for the similar temperature range. To ensure whether the solution obtained from the Monte Carlo resembles the pyrolysis phenomenon, the residual plots are drawn against each heating rate (Figure 7b−f). It is observed that the there is no symmetricity in the plot, hence the obtained solution represents the true value for pyrolysis of pine needles. Validation of results (Table 4) obtained for activation energy is done by comparing the results with different regression models. It is found that the activation energy required for thermal conversion of Cedrus deodara,39 which is very similar to pine needles, at different heating rates falls within the range of activation energies estimated by Monte Carlo scheme. Likewise, the values of E obtained from the Rayleigh model16 and the Archimedean families of Copula18 for pine needles also comply with the computed solution of Monte Carlo. Ali et al.40 derived the kinetic parameters for the components of the green Algae, and the range of variation of activation energies is within permissible limit of numerical solution obtained from the inverse solution of DAEM. However, the obtained kinetic parameters for different kind of biomass have different structural behavior; moreover, adopted methodologies propagate the computational error. Model-free methods estimate kinetic parameters at different conversion points, so the uncertainty in the estimated values adds up with the final solution. But the proposed scheme does not involve uncertainty, as it is a holistic approach to determine the kinetic parameters. It directly correlates the experimental data with kinetic model. Only way of ingression of uncertainties is the stochastic behavior of Brownian motion of the bombarded particles, or asymptotic approximations of kinetic model.

E0

solution that is similar to the results obtained for the biomass of similar characteristic. But approximation of DAEM and the computation of area under the TG plane is very stochastic, so the exact solution cannot be obtained, but the stochastic component can be refined to some extent.

■

MATERIALS AND METHODS Characteristic of Density Function p(E). To solve eq 10, it is very important to demarcate the characteristic of density function p(E) from the double exponential term of DAEM. Therefore, two different limits of interest are defined on the basis of the relative width of density function and Dexp. The shape of the total integrand changes with time, depending on the limit applied to DAEM. If the width of the density function is relatively wide to Ew, the total integrand is initially represented by the density function f(E). But as time proceeds, it is progressively surpassed by the step-like function, Dexp. The location of the maximum of the total integrand get shifted with time and the resultant shape shows skewness. On the other hand, merely the amplitude of the total integrand is affected if the density function is relatively narrow to Dexp. Thus, the final shape of integrand remains more symmetrical than that of wide 4987

DOI: 10.1021/acsomega.8b03442 ACS Omega 2019, 4, 4984−4990

ACS Omega

Article

approach may not truly be able to describe the actual process, or it is also possible that there is no analytical solution available for the given problem. Hence, simulation technique can be implemented to figure out level of uncertainties among different models. The decision-making ability of a model makes it more accurate, as the observer can correlated the actual experimental data with the proposed model.37 In this study, the basic concept of Monte Carlo simulation is used to correlate the enclosed area of TG curves with the DAEM model. It is assumed that the conversion of biomass at particular conversion has only one value, so the fuzzy characteristic in determining the kinetic parameters is refuted and the function is absolutely injective. In eq 10, there is an integral within the given range of temperature and the solution of DAEM is derived with the help of Monte Carlo algorithm. The required eqs 17 and 18 are used to collect the kinetic properties of the model. The term ϵ, which denotes the error component, comes into existence while solving DAEM using the stochastic approach. The DAEM has no analytical solution; therefore, Laplace method of integration is used to approximate the function Q(A, E, T). ÄÅ T ÉÑ Å f ÑÑ −1Å Å Q (A , E , T ) = Φ ÅÅÅ (1 − α) dT ÑÑÑÑ ÅÅÇ T0 ÑÑÖ (7)

distribution limit, albeit the location of its maximum (ze) will not be static either. But the most crucial part of the analysis is to determine whether the distribution is narrow or wide in nature. Therefore, it is necessary to determine the critical point that determines the nature of distribution function (σc). If the estimated value of the standard deviation or width (σ) is less than σc, density function will follow a narrow distribution (α ≪ RT), else it will be considered the case of a wide distribution (yw√λ ≪ 1). Therefore, it is necessary to evaluate the limiting condition of width, which can be applied to the numerical solution. An important note is included to assist in the decision-making phase whether the distribution is wide or narrow. σc =

RβtcE0 q(Rβtc + E0)

(1)

where tc denotes the time at which Dexp passes through the center of distribution function (z = 1). tc =

tmin + tmax 2

(2)

∫

tmin is the time when Dexp is at the verge to encounter the distribution function, whereas tmax is the time when Dexp has surpassed the distribution function. The time interval can be obtained by solving the following eqs 3 and 4 qσ zs + 2z w = 1 − c E0 (3) zs − z w = 1 +

qσc E0

N ji Area of rectangle zyz zz) = Q (A, E , T ) + ϵ Φ−1(∑ Mijjjj zz j Ji i=1 k {

Here • Mi are the inbound particles enclosed by the thermogram • Ji are the total number of particles bombarded on the thermogravimetric plane. Distributed Activation Energy Model. The notion of introducing distribution in reactivity was initially borne out by the work of Constable in 1925.21 He noticed that the frequency factor (A) and the activation energy (E) were highly interrelated to each other for surface-catalyzed reactions. On the basis of his findings, he proposed that the reactant surface had various sites of attack, which led to the distribution of activation energies controlling the reaction.21 In 1943, a reactivity distribution model was developed for irreversible electrical resistance changes upon annealing of evaporated metal films.22 Later, the same concept was introduced for biomass pyrolysis.23−26 The form of distribution function and its applicability are still questionable, as different distribution forms may not provide similar information about the kinetic properties of biomass pyrolysis. According to another perspective of modeling, the discrete energy distribution is also found to be effective for knowing the kinetic properties of petroleum products. It was initially introduced by the workers of Institut Francais du Petrole, but the numerical solution was first published by Ungerer and Pelet.27 The tenet of the regression-based analysis proposed for determining the kinetic mechanism of any thermal conversion process follows eq 1, where the study variable is connected to the explanatory variables via the characteristic variables (k1, k2, k3, ..., kn). The stochastic nature of the regression analysis is reflected by ϵ, so if this term gets minimized, the solution will become more precise. However, a model cannot be robust and precise at the same time if the experimental conditions are overlooked and the premise is merely based up pseudo-assumptions.

(4)

E − qσc Tf W (A(T0 + βtmin))[3 + W (A(T0 + βtmin))] = 0 1 + W (A(T0 + βtmin)) R (5)

E + qσc Tf W (A(T0 + βtmax ))W (A(T0 + βtmax )) = 0 1 + W (A(T0 + βtmax )) R

(8)

(6)

Equations 5 and 6 are the required expression for determining tmin and tmax. Monte Carlo Algorithm. Monte Carlo algorithm is a computational scheme to retrieve an unknown parameter.32 It constitutes a fundamental tool for calculating and simulating a set of maps that are independent and identically distributed with the same distribution as the set of all maps.33 The Monte Carlo provides a robust platform to evaluate absorbed dose and other quantities in the field of medical physics. It can also be used for simulation of radiation detectors.34 Guo et al.35 simulated the evolution of bonds in pyrolysis by combining the Boltzmann− Monte−Carlo algorithm with the percolation theory, which provides a new scheme to simulate coal pyrolysis. In pyrolysis problem of coal, the Boltzmann distribution is implemented to represent the portion of the bonds dissociated, whereas the Monte Carlo algorithm delineates the bond formation by using coupling or collision of radical fragment.36 Basis of Monte Carlo analysis relies on a random number generation. A random number generator provides a random value for each variable within the constrain limit set by the model and develops a probability distribution for all possible outcomes. Generally, performance assessment of the given model is done by deterministic values of parameters. But sometimes deterministic

y = f (k1 , k 2 , k 3 , ... kn , x1 , x 2 , x3 ... xn ) + ϵ 4988

(9)

DOI: 10.1021/acsomega.8b03442 ACS Omega 2019, 4, 4984−4990

ACS Omega

Article

A parallel approach to reactivity distributions emerges from the isoconversional methods of Friedman28 and Ozawa.29 These methods were mainly conceived to determine a single activation energy and reaction order, but later they had been used to establish a correlate between the activation energy (E) and conversion (α). The concept of isoconversional was initially developed by using the differential method.28 Ozawa and Flynn and Wall independently developed the integral method.29,30 In 1995, Miurra24,31 suggested a proposal that overlapping of parallel reactions constituted the reactivity distribution, so instead of the actual fraction reaction, the effective fraction reaction would rather be used in the isoconversional methods. There are two hypotheses proposed to determine the characteristic of activation energy. One of them considered that the activation energy at a particular instant of conversion is basically the average value. Although Miura24,31 denies the fact that reaction profiles can be described by a delta function. Thus, the given level of conversion is described by overlapping of parallel reactions of different activation energies. Later, he demarcated the activation energies by their significant contribution to the conversion at particular instance of time. The conversion function can be expressed by the simple differential equation describing the disappearance and appearance of reactants and products.

(1 − α)i , t =

Differentiating eq 15 with respect to time, we get ÅÄÅ d(1 − α) ÑÉÑ d(1 − α)i , t Å ÑÑ = ∑ kiÅÅÅ ÑÑ ÅÅÇ ÑÑÖ dt d t i i,t

∫

k(T (t )) = A e(−E / RT )

E is the activation energy, R denotes the universal gas constant, T represents the temperature, and A denotes the pre-exponential or frequency factor. Substituting the value of k(T(t)) in eq 9, we have ÉÑ ∞ Ñ i t 1−α= [expjjj −A e(−E / RT ) dt )ÑÑÑÑp(E) dE ÑÑÖ (18) 0 0 k

∫

(

(10)

■

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acsomega.8b03442. Appendix equations (PDF)

■

• k(T(t)) is the rate constant or Arrhenius equation depending on the temperature, T, and time, t.

*E-mail: [email protected]. Tel: +36-70-2924045. ORCID

m0 is the initial mass of the biomass. mt is the mass of biomass at time, t. mr is the residual mass of biomass. h(P) is the pressure dependence function represented by P h(P) = 1 − P .

Alok Dhaundiyal: 0000-0002-3390-0860 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by the Stipendium Hungaricum Programme. The author(s) express their sincere gratitude to Govind Ballabh Pant University of Agriculture and Technology, Pantnagar, Uttarakhand, India.

equilbrium

• P denotes the partial pressure of gaseous product. • Pequilbrium represents the vapor pressure at equilibrium. • f(α) and g(α) are the differential and integral functional forms of conversion, respectively.

■

The foundation of the distributed activation energy model is deduced from eq 1, which is for a single reaction. To generalize to an indefinite number of parallel reaction network, the reaction must be defined for the ith reaction. In the beginning, the conversion is represented by eq 4. A(1) → B(0) + C(0) att = 0

(12)

A(1 − αi) → B(αi) + C(αi) at t = ti

(13)

∑ ki(1 − α)i i

AUTHOR INFORMATION

Corresponding Author

) denotes the conversion of biomass.

1 = (1 − α)0 =

ASSOCIATED CONTENT

S Supporting Information *

Here

• • • •

∫

Equation 10 is the required form of distributed activation energy model (DAEM) for the first reaction order. ÑÉ t The term [exp ∫ −A e(−E / RT ) dt )ÑÑÑÑ refers to the temperature 0 ÑÖ integral or the double exponential term (Dexp).

(11)

m0 − m t m0 − m r

∫

where

d(α) = k(T ) f (α) h(P) (rate of appearance of product) dt

(

(16)

Now, replace the discrete summation with the continuous mathematical form for i = ∞ ÉÑ ∞ Ñ i t 1−α= [expjjj −k(T (t )) dt )ÑÑÑÑp(E) dE (17) 0 ÑÖÑ k 0

and

• α=

(15)

i

d(1 − α) = −k(T (t ))f (1 − α) dt h(P)(rate of disappearance of reactant)

∑ ki(1 − α)i ,t

REFERENCES

(1) Capart, R.; Khezami, L.; Burnham, A. K. Assessment of Various Kinetic Models for the Pyrolysis of a Microgranular Cellulose. Thermochim. Acta 2004, 417, 79−89. (2) Conesa, J. A.; Caballero, J.; Marcilla, A.; Font, R. Analysis of Different Kinetic Models in the Dynamic Pyrolysis of Cellulose. Thermochim. Acta 1995, 254, 175−192. (3) Conesa, J. A.; Marcilla, A.; Caballero, J. A.; Font, R. Comments on the Validity and Utility of the Different Methods for Kinetic Analysis of Thermogravimetric Data. J. Anal. Appl. Pyrolysis 2001, 58−59, 617− 633. (4) Hanaya, M.; Osawa, I.; Pysiak, J. J. Kinetic Equations for Thermal Dissociation Processes. J. Therm. Anal. Calorim. 2004, 76, 521−528.

(14) 4989

DOI: 10.1021/acsomega.8b03442 ACS Omega 2019, 4, 4984−4990

ACS Omega

Article

(5) Yaroshenko, A. P. Theoretical Model and Experimental Study of Pore Growth during Thermal Expansion of Graphite Intercalation Compounds. J. Therm. Anal. Calorim. 2005, 79, 515−519. (6) Criado, J. M.; Pérez-Maqueda, L. A. Sample Controlled Thermal Analysis and Kinetics. J. Therm. Anal. Calorim. 2005, 80, 27−33. (7) Burnham, A. K. Introduction to Chemical Kinetics. In Global Chemical Kinetics of Fossil Fuels; Springer International Publishing: Cham, 2017; pp 25−74. (8) Burnham, A. K.; Braun, R. L.; Gregg, H. R.; Samoun, A. M. Comparison of Methods for Measuring Kerogen Pyrolysis Rates and Fitting Kinetic Parameters. Energy Fuels 1987, 1, 452−458. (9) Anthony, D. B.; Howard, J. B. Coal Devolatilization and Hydrogastification. AIChE J. 1976, 22, 625−656. (10) Pitt, G. J. The Kinetics of the Evolution of Volatile Products from Coal. Fuel 1962, 41, 267−274. (11) Hanbaba, P. Reaktionkinetische Untersuchungen Sur Kohlenwasserstoffenbindung Aus Steinkohlen Bie Niedregen Aufheizgeschwindigkeiten, University of Aachen, 1967. (12) Galgano, A.; Blasi, C. Di Modeling Wood Degradation by the Unreacted-Core-Shrinking Approximation. Ind. Eng. Chem. Res. 2003, 42, 2101−2111. (13) Ferdous, D.; Dalai, A. K.; Bej, S. K.; Thring, R. W. Pyrolysis of Lignins: Experimental and Kinetics Studies. Energy Fuels 2002, 16, 1405−1412. (14) Burnham, A. K.; Braun, R. L. Global Kinetic Analysis of Complex Materials. Energy Fuels 1999, 13, 1−22. (15) Dhaundiyal, A.; B. Singh, S. Mathematical Insight to NonIsothermal Pyrolysis of Pine Needles for Different Probability Distribution Functions. Biofuels 2018, 9, 647−658. (16) Cai, J. M.; Liu, R. H. Parametric Study of the Nonisothermal n Th-Order Distributed Activation Energy Model Involved the Weibull Distribution for Biomass Pyrolysis. J. Therm. Anal. Calorim. 2007, 89, 971−975. (17) Dhaundiyal, A.; Singh, S. B. Approximations to the NonIsothermal Distributed Activation Energy Model for Biomass Pyrolysis Using the Rayleigh Distribution. Acta Technol. Agric. 2017, 20, 78−84. (18) Dhaundiyal, A.; Singh, S. B.; Hanon, M. M. Study of Distributed Activation Energy Model Using Bivariate Distribution Function, f ( E 1, E 2 ). Therm. Sci. Eng. Prog. 2018, 5, 388−404. (19) Dhaundiyal, A.; Singh, S. B.; Hanon, M. M. Application of Archimedean Copula in the Non-Isothermal n Th Order Distributed Activation Energy Model. Biofuels 2018, 8, 1−12. (20) de Caprariis, B.; De Filippis, P.; Herce, C.; Verdone, N. DoubleGaussian Distributed Activation Energy Model for Coal Devolatilization. Energy Fuels 2012, 26, 6153−6159. (21) Zhang, J.; Chen, T.; Wu, J.; Wu, J. Multi-Gaussian-DAEMReaction Model for Thermal Decompositions of Cellulose, Hemicellulose and Lignin: Comparison of N2 and CO2 Atmosphere. Bioresour. Technol. 2014, 166, 87−95. (22) Constable, F. H. The Mechanism of Catalytic Decomposition. Proc. R. Soc. A 1925, 108, 355−378. (23) Vand, V. A Theory of the Irreversible Electrical Resistance Changes of Metallic Films Evaporated in Vacuum. Proc. Phys. Soc. 1943, 55, 222−246. (24) Miura, K. A New and Simple Method to Estimate f(E) and K0(E) in the Distributed Activation Energy Model from Three Sets of Experimental Data. Energy Fuels 1995, 9, 302−307. (25) Cai, J.; Li, T.; Liu, R. A Critical Study of the Miura−Maki Integral Method for the Estimation of the Kinetic Parameters of the Distributed Activation Energy Model. Bioresour. Technol. 2011, 102, 3894−3899. (26) Lakshmanan, C. C.; Bennett, M. L.; White, N. Implications of Multiplicity in Kinetic Parameters to Petroleum Exploration: Distributed Activation Energy Models. Energy Fuels 1991, 5, 110−117. (27) Ungerer, P.; Pelet, R. Extrapolation of the Kinetics of Oil and Gas Formation from Laboratory Experiments to Sedimentary Basins. Nature 1987, 327, 52−54. (28) Friedman, H. L. Kinetics of Thermal Degradation of CharForming Plastics from Thermogravimetry. Application to a Phenolic Plastic. J. Polym. Sci., Part C: Polym. Symp. 1964, 6, 183−195.

(29) Ozawa, T. A New Method of Analyzing Thermogravimetric Data. Bull. Chem. Soc. Jpn. 1965, 38, 1881−1886. (30) Flynn, J. H.; Wall, L. A. A Quick, Direct Method for the Determination of Activation Energy from Thermogravimetric Data. J. Polym. Sci. Part B Polym. Lett. 1966, 4, 323−328. (31) Miura, K.; Maki, T. A Simple Method for Estimating f (E) and k 0 ( E ) in the Distributed Activation Energy Model. Energy Fuels 1998, 12, 864−869. (32) Fishman, G. S. Monte Carlo, Concepts, Algorithms, and Applications; Springer Science & Business Media, 1995. (33) Cho, W. K. T.; Liu, Y. Y. Sampling from Complicated and Unknown Distributions: Monte Carlo and Markov Chain Monte Carlo Methods for Redistricting. Phys. A 2018, 506, 170−178. (34) Polo, I. O.; Santos, W. de S.; Caldas, L. V. E. Determination of Correction Factors in Beta Radiation Beams Using Monte Carlo Method. Appl. Radiat. Isot. 2018, 140, 50−54. (35) Guo, X.; Liu, Z.; Xiao, Y.; Xu, X.; Xue, X.; Liu, Q. The BoltzmannMonte-Carlo-Percolation (BMCP) Model on Pyrolysis of Coal: The Volatiles’ Reactions. Fuel 2018, 230, 18−26. (36) Guo, X.; Liu, Z.; Liu, Q.; Shi, L. Modeling of Kraft Lignin Pyrolysis Based on Bond Dissociation and Fragments Coupling. Fuel Process. Technol. 2015, 135, 133−149. ́ (37) Yoriyaz, H. Método de Monte Carlo: Principios e Aplicações Em ́ Fisica Médica. Rev. Bras. Fis. Med. 2009, 3, 141−149. (38) Dhaundiyal, A.; Singh, S. B.; Hanon, M. M.; Rawat, R. Determination of Kinetic Parameters for the Thermal Decomposition of Parthenium Hysterophorus. Environ. Clim. Technol. 2018, 22, 5−21. (39) Dhaundiyal, A.; Hanon, M. M. Calculation of Kinetic Parameters of the Thermal Decomposition of Residual Waste of Coniferous Species: Cedrus Deodara. Acta Technol. Agric. 2018, 21, 75−81. (40) Ali, I.; Naqvi, S. R.; Bahadar, A. Kinetic Analysis of Botryococcus braunii Pyrolysis Using Model-Free and Model Fitting Methods. Fuel 2018, 214, 369−380. (41) Slopiecka, K.; Bartocci, P.; Fantozzi, F. Thermogravimetric Analysis and Kinetic Study of Poplar Wood Pyrolysis. Appl. Energy 2012, 97, 491−497. (42) Font, R.; Conesa, J. A.; Moltó, J.; Muñoz, M. Kinetics of Pyrolysis and Combustion of Pine Needles and Cones. J. Anal. Appl. Pyrolysis 2009, 85, 276−286. (43) Dhaundiyal, A.; Mohammad, A. T.; Laszlo, T. Thermo-Kinetics of Forest Waste Using Model-Free Methods; Universitas Scientiarums, 2019; Vol. 24, pp 1−31.

■

NOTE ADDED AFTER ISSUE PUBLICATION This paper was originally published on March 7, 2019. Due to a production error, the paper was published before all requested corrections were incorporated. The corrected version was reposted on March 8, 2019.

4990

DOI: 10.1021/acsomega.8b03442 ACS Omega 2019, 4, 4984−4990