Analysis of Pyrolysis Kinetic Model for Processing of Thermogravimetric Analysis Data

Abstract

Pyrolysis has profound implications for coal as a raw material to make phase change material (PCM). It is necessary to derive a pyrolysis kinetic model for predicting the yield of volatiles and reaction performance during pyrolysis of coal, which is of significant importance for its thermal processing. The devolatilization of coal is characterized by thermogravimetric analysis (TGA) at different heating rates, and many kinetic models can be achieved by analyzing the TGA data. This work was aimed to find an appropriate model to describe the pyrolysis of coal and took Zhundong coal as an example. Four types of isoconversion kinetic methods, that is, Friedman, Flynn-Wall-Ozawa (FWO), Kissinger-Akahira-Sunose (KAS), Miura-Maki method, and different distributed activation energy models (DAEM) were employed here to fit TGA data for pyrolysis of Zhundong coal. The pre-exponential factors and activation energies obtained from different kinetic models were analyzed. An m-nth-DAEM was developed by considering that m classes of reactions took place with the same pre-exponential factor k0 but different distribution activation energy following logistic distribution or Gaussian distribution. The results showed that the FWO model was better for description of pyrolysis process of Zhundong coal, and the 2-nth-DAEM assuming Gaussian distribution of activation energy gave the best fit for the TGA data of Zhundong coal. The research provides a valuable reference to the development of thermal utilization technology of Zhundong coal.

1. Introduction

Energy storage technology can solve the contradiction between energy supply and demand in time and space, so it is an effective way to improve energy utilization [1]. Thermal energy storage is widely used in industrial and civilian applications, so it occupies an extremely important position in the field of energy storage technology. Phase change materials (PCM) will absorb or release a large amount of latent heat for energy storage when the state of matter changes. As a raw material for the preparation of phase change materials, coal has great energy storage density, and the output temperature and energy are relatively stable, so it has good application prospects [2]. Zhundong coalfield is the largest compositive basin in China, which has a large reserve of 390 Gt and is estimated to meet China’s next 100-year coal consumption requirement [3]. In addition, Zhundong coal is more environmentally friendly than other types of coal owing to its extremely low sulfur and ash contents, which has a great significance to coal industry [4, 5]. Therefore, the investigation on Zhundong coal has theoretical and practical significance for the utilization of coal resources of China.

Thermal pyrolysis is the first stage reaction in most of coal thermal conversion processes (combustion, gasification, and liquefaction) and it is also one of the most promising technologies for clean and effective utilization of low-rank coal. During the pyrolysis process, various non-condensible gas and liquid phase matter such as tar are generated, which will lead to the weightlessness of coal samples. The weight losses of the coal range from 10 to 55%, which is affected by pyrolysis temperature, coal size, pyrolysis time and other factors such as pressure, atmosphere, heating rate and mineral content, etc. [6, 7]. The heating rate and pyrolysis temperature have an important influence on the pyrolysis characteristics of coal. The heating rate affects the concentration gradient of the product inside the particle and the ratio of gas phase substance to tar in the volatiles. The pyrolysis temperature determines the pyrolysis of the macromolecule reaction rate in the coal structure and the total release of volatile matters. A thermogravimetric analyzer (TGA) is routinely employed by many researchers to characterize moisture and volatile contents of various coals. It enables the user to monitor the weight change of a certain quality of coal samples as a function of time and temperature, which endue this method to incorporate the advantages of simple operation and good repeatability [8]. From the TG curve, we can obtain the key parameters such as the temperature of initial point of pyrolysis volatilization, the maximum release rate of volatile matter, and its corresponding temperature.

It is necessary to obtain an accurate knowledge of the pyrolysis process, since all the models of coal thermal conversion processes involve pyrolysis, and the kinetic behavior of these products in pyrolysis is fundamental to the optimization of their use. However, the actual phenomena of pyrolysis are complex and the kinetic parameters of which are difficult to obtain. To date, a large number of kinetic models have been proposed in the literature [5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]. The models for processing TGA data and explaining the kinetic process of coal pyrolysis include first-order single reaction model (SRM), isoconversional or model-free methods, distributed activation energy model (DAEM). The SRM is the simplest, and it is applicable only over a narrow range of temperatures, thus limiting its applicability. Isoconversional or model-free methods are also used to process TGA data in which the kinetics parameters can be obtained without an explicit kinetics model [9]. Friedman method [10], Kissinger method [11], Flynn-Wall-Ozawa [12], and Miura-Maki method [13] are some of the outstanding methods in this classification. Several TGA curves at different heating rates for the same value of conversion are required for these models to determine activation energy for each conversion point. In recent years, nonlinear regressions (least square optimization) are increasingly used to process the pyrolysis process of coals. More often than not, the distributed activation energy model (DAEM), originally developed by Pitt [14], has become the distinguished representative of them. It assumes that the coal pyrolysis should include a large number of independent, parallel, and irreversible n-order reactions with different activation energies. This distribution of activation energy corresponding to these infinite reactions is often represented by a distribution function of activation energy $f\left(E\right)$, which is usually used as logistic distribution or Gaussian distribution. Based on the DAEM, a two-Gaussian DAEM (2-DAEM) [15] has been developed by dividing the pyrolysis process into two steps. According to these indexes, the distributed activation energy model will be better improved and applied.

In the present paper, TGA experiments with pyrolysis of Zhundong coal at different heating rate have been carried out. The data was first processed using various isoconversional methods and then processed using the distributed activation energy model (DAEM). An m-nth-DAEM was developed to predict the weight loss of Zhundong coal, considering that $m$ classes of reactions take place having the same pre-exponential factor $k_0$ and different distribution of activation energy.

2. Experiment

The experimental tests were with a TGA (Netzsch Model STA449F3). Before the pyrolysis experiments, coals were grounded and the finer particle fractions in a size range 0.075–0.15 (mm) were selected. Table 1 shows the results of the proximate and ultimate analyses of the coal samples. Small samples were used to ensure uniform heating and to avoid problems of transport phenomena through the sample bed in the crucible. The instrument was first purged for 30 min with nitrogen to remove any remains of air, then further heated from room temperature to 1000°C at a heating rate of 10 $\left(\mathrm{K}\bullet {\text{min}}^{-1}\right)$, 20 $\left(\mathrm{K}\bullet {\text{min}}^{-1}\right)$, 30 $\left(\mathrm{K}\bullet {\text{min}}^{-1}\right)$ and 50 $\left(\mathrm{K}\bullet {\text{min}}^{-1}\right)$, respectively. The tests were conducted in an ${\mathrm{N}}_2$ inert atmosphere, and the samples were maintained at the maximum temperature for 10 min until a constant mass was obtained at this temperature. All experiments at given heating rate were conducted two times and averaged to eliminate measure error. The experiments showed reproducibility within the range of ±3.7%.

Proximate analysis (%)				Ultimate analysis (%)
Mad	Aad	Vad	FCad	Cd	Hd	Od	Nd	St,d
9.96	15.46	22.49	52.09	70.55	3.49	7.64	0.62	0.53

Table 1.

Proximate and ultimate analyses of raw Zhundong coal samples.

Typical decomposition data of Zhundong coal obtained from the TGA experiments [9] are shown in Figure 1. The key characteristic parameters related to the main pyrolysis process are summarized in Table 2. The initial decomposition temperature is denoted as $T_0$, the corresponding temperature at the maximum point of the weight loss rate is denoted as $T_p$, the maximum weight loss rate at $T_p$ is denoted as $R\infty$, and $V_f$ is the weight loss corresponding to the final temperature. As can be seen from Figure 1, the mass fraction of releasing volatiles decreases with the increase of heating rate at the same temperature. The reason is that the increasing heating rate results in a thermal hysteresis, which makes the temperature inside the sample lower than the appearance, thereby leading to the value of $1-\mathrm{\alpha }$ shifts to lager temperature. As can be seen from Figure 1, the TGA curve first decreases slightly due to small molecular gases being removed as the weak bonds are decomposed. Then tar-based volatile components and the formation of char result in a significant decrease of the TGA curve. At the end of the curve, it is close to horizontal, where a cross-linking reaction occurs. The TGA curve contains the main chemical reactions in the coal pyrolysis process, so it is the basis for the kinetic model investigation on coal pyrolysis.

Figure 1.

TGA curves of Zhundong coal samples during pyrolysis at different heating rates taken from literature [9].

$\beta$ (K∙min^−1)	$T_0$ (K)	$T_p$ (K)	$R_{\infty }$ (%∙min^−1)	$V_f$ (%)
10	425.39	718.95	0.90	35.8
20	458.89	731.25	1.75	36.4
30	482.89	735.25	2.85	38.0
50	489.86	752.75	5.34	37.9

Table 2.

Pyrolysis characteristic values of coal samples.

3. Theoretical basis for thermal kinetics

For non-isothermal and heterogeneous systems, the thermal analysis kinetics is usually evolved from the theory of isothermal and homogeneous systems. The kinetics of reactions in isothermal and homogeneous systems for solids are generally described by Eq. (1)

$$\frac{\mathit{dV}/V^{\ast }}{\mathit{dt}}=\mathit{kf}\left(\frac{V}{V^{\ast }}\right)$$E1

where V is total mass of volatiles released at time ‘t’,$V^{\ast }$ is the total mass of volatiles originally available for the reaction, $f\left(\frac{V}{V^{\ast }}\right)$ is the reaction model mechanism function, and $k$ is the kinetic rate constant. The reaction rate constant varies with temperature and is assumed to follow the Arrhenius law [16].

$$k=k_0{e}^{\left(-\frac{E}{\mathit{RT}}\right)}$$E2

In the above formula, $k_0$ refers to pre-exponential factor, $E$ is the activation energy,$R$ is universal gas constant, and $T$ is the thermodynamic temperature. This equation is basically applicable to most elementary reactions and complex reactions, where $k_0$ and $E$ are two important parameters that need to be combined with experimental and theoretical models. If the TGA experiment is conducted at a linear heating rate of $\beta$ ($\mathrm{K}\bullet {\text{min}}^{-1}$), the temperature and time coordinates will be related to each other as:

$$T=T_0+\mathit{\beta t}$$E3

Thus,$\beta =\mathit{dT}/\mathit{dt}$. Substituting Eq. (2) and (3) into Eq. (1) gives,

$$\beta \frac{\mathit{dV}/d{V}^{\ast }}{\mathit{dT}}=k_0{e}^{\left(-\frac{E}{\mathit{RT}}\right)}f\left(\frac{V}{V^{\ast }}\right)$$E4

Assuming $\alpha =\frac{V}{V^{\ast }}$, then Eq. (4) becomes,

$$\beta \frac{\mathit{d\alpha }}{\mathit{dT}}=k_0{e}^{\left(-\frac{\mathrm{E}}{\mathit{RT}}\right)}f\left(\alpha \right)$$E5

For different reaction models, there are differences in the expression of $f\left(\alpha \right)$. One of the most commonly used reaction models is the reaction order model, which can be expressed as,

$$f\left(\alpha \right)={\left(1-\alpha \right)}^n$$E6

The core issue of the general reaction kinetics of coal pyrolysis is determining the three factors $k_0$, $f\left(\alpha \right)$, and $E$ [17], which is used as a key indicator to describe the process of devolatilization. In general, the integral form of the dynamical mechanism function $f\left(\alpha \right)$ is,

$$g\left(\alpha \right)={\int }_0^{\alpha }\frac{\mathit{d\alpha }}{f\left(\alpha \right)}$$E7

The above is the basic knowledge of thermal analysis dynamics. Next, two methods, that is, the isoconversional model and the distributed activation energy model (DAEM), will be used to align the TGA data of Zhundong coal pyrolysis.

4. Isoconversional methods

4.1. Isoconversional theories

The basic assumption of the isoconversional method is: under different conditions of temperature rise, the activation energy corresponding to the same conversion rate remains the same [18], and the pyrolysis is specified to be a first-order reaction, that is, n = 1. Thus, $g\left(\alpha \right)$ can be written as,

$$g\left(\alpha \right)=-\mathit{ln}\left(1-\alpha \right)$$E8

The activation energy and pre-exponential factor for each conversion point can be obtained by analyzing several TGA curves at different heating rates for the same value of conversion $\left(\alpha \right)$. Next, we will focus on the typical isoconversional method for discussion, including a temperature differential method and three temperature integration methods.

4.1.1. Friedman method

Friedman model is a typical representative of the temperature differential method. Taking the logarithm of the Eq. (5) on both sides of the equal sign to obtain,

$$\mathit{ln}\left(\beta \frac{\mathit{d\alpha }}{\mathit{dT}}\right)=\mathit{ln}\left(k_{0}f\left(\alpha \right)\right)-\frac{E}{\mathit{RT}}$$E9

The above formula is the general formula of the Friedman method. Thus, the plot of $\mathit{ln}\left(\beta \frac{\mathit{d\alpha }}{\mathit{dT}}\right)$ versus $1/T$ at the same conversion levels should be a straight line, the slope and intercept of which can be employed to specify the values of $E$ and $k_0$, respectively.

4.1.2. Flynn-Wall-Ozawa method

Combining Eq. (5), Eq. (8) can be rewritten as

$$g\left(\alpha \right)={\int }_0^{\alpha }\frac{\mathit{d\alpha }}{f\left(\alpha \right)}={\int }_{T_0}^T\frac{k_0}{\beta }{e}^{\frac{-E}{\mathit{RT}}}\mathit{dT}\approx {\int }_0^T\frac{k_0}{\beta }{e}^{\frac{-E}{\mathit{RT}}}\mathit{dT}=\frac{k_{0}E}{\mathit{\beta R}}p\left(u\right)$$ E10

where $u=E/\left(\mathit{RT}\right)$ and $p\left(u\right)$ can be simplified by using integral approximation,

$$p\left(u\right)={\int }_{\infty }^u\left(\frac{e^{-u}}{u^2}\right)\mathit{du}=\frac{e^{-u}}{u^2}\left|\genfrac{}{}{0pt}{}{u}{\infty }-{\int }_{\infty }^u{e}^{-u}d{u}^{-2}\right.=\frac{e^{-u}}{u^2}-{\int }_{\infty }^{u}2u^{-3}d{e}^{-u}$$

$$=\frac{e^{-u}}{u^2}-\frac{2e^{-u}}{u^3}+{\int }_{\infty }^{u}6u^{-4}d{e}^{-u}$$

$$=\frac{e^{-u}}{u^2}-\frac{2e^{-u}}{u^3}+\frac{6e^{-u}}{u^4}\left|\genfrac{}{}{0pt}{}{u}{\infty }-{\int }_{\infty }^u{e}^{-u}\right.d\frac{6}{u^4}=\cdots$$

$$=\frac{e^{-u}}{u^2}\left(1-\frac{2!}{u}+\frac{3!}{u^2}-\frac{4!}{u^3}+\cdots \right)$$ E11

Taking the first two items in parentheses on the right side of Eq. (11), we can get the approximate expression of $p\left(u\right)$ ,

$$p\left(u\right)=\frac{e^{-u}\left(u-2\right)}{u^3}$$ E12

Performing logarithmic operation on both sides of the Eq. (12),

$$\mathit{lnp}\left(u\right)=-u+\mathit{ln}\left(u-2\right)-3\mathit{lnu}$$ E13

For a given solid fuel, in the temperature range of the thermal reaction stage, the value of $u$ is generally ranging from 20 to 60 [19], which is, 20 ≤ u ≤ 60, then −1 ≤ (u − 40)/20 ≤ 1, take $s=\left(u-40\right)/20$ , thus it gives,

$$u=20s+40$$ E14

Substituting Eq. (14) into Eq. (13), it gives,

$$\mathit{lnp}\left(u\right)=-u-3\mathit{ln}40+\mathit{ln}\left(38\right)+\mathit{ln}\left(1+\frac{10s}{19}\right)-3\mathit{ln}\left(1+\frac{s}{2}\right)$$ E15

Using the Taylor series expansion for the two logarithmic expressions to the right of Eq.(15) and taking a first-order approximate expression of each logarithmic, the following expression can be obtained:

$$\mathit{lnp}\left(u\right)\approx -5.3308-1.0516u$$ E16

$$\mathrm{p}\left(u\right)=0.0048e^{-1.0516u}$$ E17

Substituting Eq. (17) into Eq. (10) yields,

$$g\left(\alpha \right)\approx \frac{k_{0}E}{\mathit{\beta R}}\left(0.0048e^{-1.0516\frac{E}{\mathit{RT}}}\right)$$ E18

Eq. (18) can be rewritten as,

$$\mathit{ln\beta }=\mathit{ln}\left(\frac{k_{0}E}{\mathit{Rg}\left(\alpha \right)}\right)-5.331-1.052\frac{E}{\mathit{RT}}$$ E19

Eq. (19) is the integral form of the FWO model. The plot of $\mathit{ln\beta }$ versus $1/T$ at the same $\alpha$ value for different heating rates gives the $E$ and $k_0$ value corresponding to the selected $\alpha$ .

4.1.3. Kissinger-Akahira-Sunose method

Taking the first item in parentheses on the right side of Eq. (11), the approximate expression of $p\left(u\right)$ can be rewritten as,

$$p\left(u\right)=e^{-u}/u^2$$ E20

Substituting Eq. (20) into Eq. (10) yields,

$$g\left(\alpha \right)=k_{0}E{e}^{-u}/\left(\mathit{\beta R}{u}^2\right)$$ E21

Performing Logarithmic operations on the two sides of Eq. (21) and then substituting $u=E/\left(\mathit{RT}\right)$ into the formula,

$$\mathit{ln}\left(\frac{\beta }{T^2}\right)=\mathit{ln}\left(\frac{k_{0}R}{\mathit{Eg}\left(\alpha \right)}\right)-\frac{E}{\mathit{RT}}$$ E22

Eq. (22) is the holonomic form of the KAS model. Using Eq. (22), we can estimate both $E$ and $k_0$ from the Arrhenius plot of $\mathit{ln}\left(\frac{\beta }{T^2}\right)$ versus 1/T at the selected $\alpha$ value and $\beta$ values.

4.1.4. Miura-Maki method

In the Miura-Maki integral method, the approximation of the temperature integral is the same as the KAS method. Rearranging Eq. (22) and inserting Eq. (8) gives,

$$1-\mathrm{\alpha }\cong \text{exp}\left[-\frac{{\mathrm{k}}_0\mathrm{R}{\mathrm{T}}^2}{\mathrm{\beta E}}\text{exp}\left(-\frac{\mathrm{E}}{\mathrm{RT}}\right)\right]=\varnothing \left(E_s,T\right)$$ E23

To derive the Miura-Maki integral method, the assumption of distributable activation energy model is used to approximate the right side of Eq. (23) as a step function at $E=E_s$ for a selected temperature. Furthermore, the activation energy $E_s$ can be chosen to meet $\varnothing \left(E_s,T\right)=0.58$ [12]. Substituting Eq. (8) into Eq. (22) yields

$$\mathit{ln}\left(\frac{\beta }{T^2}\right)=\mathit{ln}\left(\frac{k_{0}R}{E}\right)-\text{ln}\left[-\mathit{ln}\left(1-\alpha \right)\right]-\frac{E}{\mathit{RT}}$$ E24

Inserting the value of $\varnothing \left(E_s,T\right)$ into Eq. (24) gives

$$\mathit{ln}\left(\frac{\beta }{T^2}\right)=\mathit{ln}\left(\frac{k_{0}R}{E}\right)+0.6075-\frac{E}{\mathit{RT}}$$ E25

Eq. (25) is the final expression of Miura-Maki method. The linear expression used to calculate the activation energy is similar to the KAS method, but there are some differences in the calculation results of the pre-exponential factor.

Generally speaking, the core methods of the above four models for calculating activation energy and pre-exponential factor are the same. The procedure to estimate $f\left(E\right)$ and $k_0$ using the four methods is as follows:

1. Measure $\alpha$ versus $T$ relationships at least three different heating rates;

2. Calculate the values of $\mathit{ln}\left(\beta \frac{\mathit{d\alpha }}{\mathit{dT}}\right)$ , for Friedman model, $\mathit{ln\beta }$ , for FWO model, and $\mathit{ln}\left(\frac{\beta }{T^2}\right)$ , for KAS/Miura-Maki model, at selected $\alpha$ values from the $\alpha$ versus $T$ relationships obtained for different heating rates, respectively.

3. Plot $\mathit{ln}\left(\beta \frac{\mathit{d\alpha }}{\mathit{dT}}\right)$ versus $1/T$ , for Friedman model, $\mathit{ln\beta }$ versus $1/T$ , for FWO model,and $\mathit{ln}\left(\frac{\beta }{T^2}\right)$ versus $1/T$ , for KAS/Miura-Maki model, at the selected $\alpha$ and then determine the activation energies $E$ from the slopes and $k_0$ from the intercept, respectively.

The Friedman model is a representative of the differential method in the isoconversional method. The FWO model, KAS model, and Miura-Maki method are three typical integration methods. The integral equations all involve a general temperature integral, ${\int }_0^T\frac{k_0}{\beta }{e}^{\frac{-E}{\mathit{RT}}}\mathit{dT}$ , which cannot be solved exactly, and different approximation methods have been proposed for this integral [13]. The KAS method and the Miura-Maki method express the above integral expression by a step-by-step integration formula, as shown in Eq. (11), and adopt the first item of the integration formula, which produces a lower accuracy. The FWO model performs two terms of the step-integration method, which improves the prediction accuracy. Nevertheless, Taylor approximation for the logarithm gives some errors. Two approximations may cause the FWO method to have some deviation from the experimental value. The Friedman method avoids the problem of integral approximation, but the numerical differential calculation will amplify the experimental noise. Unlike the three methods, the Miura-Maki method uses not only the temperature integral approximation, the same as KAS method, but also the assumption in the DAEM, as discussed below. It is obvious that all the isoconversional methods will cause some inevitable errors, so it is necessary to use a method with a good linear correlation to describe the TGA data.

4.2. Processing thermogravimetric analysis data by isoconversional methods

In order to obtain the key parameters of the Friedman, FWO, KAS, and Miura-Maki models, the value of $\mathrm{\alpha }$ is chosen to ranges from 0.05 to 0.95 with an interval of 0.05 in this paper [10]. The activation energy and pre-exponential factor at selected $\alpha$ values can be obtained by the slope and intercept of the linear fit, respectively, as shown in Figures 4–6 in literature [9]. The activation energy and pre-exponential factor obtained with the isoconversional and Miura-Maki methods with respect to the conversion extent $\left(\alpha \right)$ have been shown in Figure 2. It can be noted that activation energies obtained from the four methods were almost the same, and they all increase continuously with increased conversion rates. This is because the refractory molecules in coal must require more energy to be released in the form of volatiles. Note also that as the reaction proceeds, the logarithmic value of the pre-exponential factor first increases and then levels off. This indicates that as the temperature rises, the frequency of molecular collisions in the coal rises and the reaction tends to be strenuous; when the reaction proceeds to a certain degree, that is, $\mathrm{E}=170\left(\mathrm{kJ}\bullet {\mathrm{mol}}^{-1}\right)$ , most of the molecules are in an activated state, and the reaction rate tends to be level, that is, pre-exponential factors fluctuate within an order of magnitude. The relationship between $\mathit{ln}{k}_0$ and $E$ for Zhundong coal are also shown in Figure 2. It shows that before the activation energy was 170 $\left(\mathrm{kJ}\bullet {\mathrm{mol}}^{-1}\right)$ , $\mathit{ln}{k}_0$ and $E$ showed a clear linear relationship, the so-called “kinetic compensation effect” [20], which would bring about some unavoidable errors to the solution of the activation energy and the pre-exponential factor of the reaction. When the activation energy is greater than 170 $\left(\mathrm{kJ}\bullet {\mathrm{mol}}^{-1}\right)$ , the value of the pre-exponential factor fluctuates above and below $\mathit{ln}{k}_0=23$ . The fixed value of $k_0$ obtained here is to be used as a constant value of the $k_0^{\prime }$ in the distributed activation energy model. The range of kinetic parameters obtained for the four models is shown in Table 3. Since the three integration methods use certain approximation methods in the derivation process, and the differential method will amplify the experimental noise, a method with better linear correlation is used here to discuss the distribution of activation energy in the Zhundong coal. The linear correlation coefficients of the four models are listed in Table 4. It can be seen that when the value of $\mathrm{\alpha }$ reaches 0.75 or more, the linear correlation coefficient square $R^2$ of the four methods shows a downward trend. For the FWO method, the linear correlation coefficient is kept above 0.92 when $\mathrm{\alpha }$ is between 0.15 and 0.95, and the value of $R^2$ is higher than the other three methods. It shows that the FWO method describes the experimental data more accurately, and the obtained kinetic parameters are more reliable over a wide range of temperatures. Figure 3 was obtained by differentiating the $\mathrm{\alpha }$ versus E (calculated by FWO method) relationship. It can be seen that the $f\left(E\right)$ shows a unimodal distribution, with the peak at 173 $\left(\mathrm{kJ}\bullet {\mathrm{mol}}^{-1}\right)$ and the activation energy mainly in the range of 160–180 $\left(\mathrm{kJ}\bullet {\mathrm{mol}}^{-1}\right)$ . The calculated activation energy distribution can be well fitted using the Gaussian function, where the correlation coefficient is 0.96. Comparison between the experimental and the simulated data by the FWO method (the best fit among the three methods) is shown in Figure 4. The error between the simulated values and the experimental results of $\mathrm{\alpha }$ is less than 12.3%, so FWO model achieves a good agreement with the experimental data over the whole range of conversion. But it needs more than three TGA curves to obtain the distributions of kinetic parameters, which will produce some inconvenience. Jain et al. [23] used Indian coal as an experimental sample and found that the results obtained by using the isoconversional methods at higher temperature had poor linear correlation. Therefore, the activation energy calculated by the isoconversional methods can only reflect the condition within a certain temperature range. However, the kinetic parameters obtained from isoconversional methods can provide inputs for selecting appropriate distributions for the DAEM.

Figure 2.

Comparison of activation energy and pre-exponential factor obtained with four different methods for Zhundong coal samples.

Method	$E_{\mathit{min}}(\mathrm{kJ}·{\mathrm{mol}}^{-1}$ )	$E_{\mathit{max}}(\mathrm{kJ}·{\mathrm{mol}}^{-1}$ )	$k_{\mathit{min}}$ ( ${\mathrm{s}}^{-1}$ )	$k_{\mathit{max}}$ ( ${\mathrm{s}}^{-1}$ )
FWO	55	237	120.16	1.58E + 10
KAS	48.5	231.6	9.46	1E + 10
Miura-Maki	48.5	231.6	100.48	3.82 E + 9
Friedman	69.7	247	528.75	5.07E + 9

Table 3.

Summary of the kinetic parameters obtained from different methods.

${\mathbf{R}}^2$ $\mathit{\alpha }$ methods	0.15	0.25	0.35	0.45	0.55	0.65	0.75	0.85	0.95
FWO	0.999	0.998	0.998	0.999	0.999	0.998	0.990	0.973	0.927
KAS/Miura-Maki	0.998	0.996	0.997	0.999	0.998	0.998	0.988	0.969	0.916
Friedman	0.999	0.996	0.997	0.998	0.997	0.993	0.973	0.974	0.881

Table 4.

Linearly dependent coefficients for different methods at various $\mathrm{\alpha }$ values.

Figure 3.

Pyrolysis activation energy distributions $f\left(E\right)$ based on FWO model.

Figure 4.

Comparison of calculation results from FWO method ( ${\mathrm{\alpha }}_{\mathit{cal}})$ and experimental data $\left({\mathrm{\alpha }}_{\mathit{exp}}\right)$ of conversion.

5. Distributed activation energy model

5.1. Distributed activation energy model theory

The DAEM assumes that the original materials of coal contain plenty of constituents, which are numbered i = 1, 2… $j$ . These components have different reaction activation energies [21]. Only when the reaction temperature reaches a certain value will they be released in the form of a volatile fraction. These reactions for diverse components are expressed as n-order reactions and are all irreversible and independent of the reactions in other components. The instantaneous volatilized mass fraction for the $i_{\mathit{th}}$ constituent is denoted by $V_i\left(t\right),$ and the total released mass fraction for the $i\mathrm{th}$ constituent is $V_i^{\ast }\left(t\right)$ . According to Eq. (4) and Eq. (6), the reaction rate of the mass fraction of releasing volatiles can be expressed as,

$$\beta \frac{d{V}_i/d{V}_i^{\ast }}{\mathit{dT}}=k_0{e}^{\left(-\frac{E_i}{\mathit{RT}}\right)}{\left(1-\frac{V_i}{V_i^{\ast }}\right)}^n$$ E26

To some solid state reactions, the pre-exponential factor is connected with the temperature through the following relationship,

$$k_0=k_0^{,}{T}^{\gamma }$$ E27

Combining Eq. (27), the integrating solution to Eq. (26) may be rewritten as,

$$\begin{array}{c}\frac{V_i}{V_i^{\ast }}=1-{\left[1-\left(1-n\right)\frac{k_0^{,}}{\beta }{\int }_0^T{T}^x{e}^{-E_i/\mathit{RT}}\mathit{dT}\right]}^{1/\left(1-n\right)},n\ne 1\\ \frac{V_i}{V_i^{\ast }}=1-\mathit{exp}\left[-\frac{k_0^{,}}{\beta }{\int }_0^T{T}^x{e}^{-E_i/\mathit{RT}}\mathit{dT}\right],n=1\end{array}$$ E28

Considering the continuous reactants have different activation energies, the DAEM devotes to describe the energy distribution of the infinite complicated set of reactions, that is $j\to \infty$ , as a function $f\left(E\right)$ . Thus, the amount of volatiles loss that has an activation energy between $E$ and $E+\mathit{dE}$ can be expressed as,

$$d{V}^{\prime }=V^{\ast }f\left(E\right)\mathit{dE}$$ E29

The integration of Eq. (29) brings about,

$${\int }_0^{\infty }f\left(E\right)\mathit{dE}=1$$ E30

The value of $\mathrm{d}{V}^{\prime }$ is actually the same concept as $V_i^{\ast }$ , which accounts for the total released mass fraction for the $i\mathrm{th}$ constituent which has an activation energy between $E$ and $E+\mathit{dE}$ . Let replace $V_i^{\ast }$ and $V_i$ with $\mathrm{d}{V}^{\prime }$ and $\mathit{dV},$ respectively. Introducing Eq. (29) to Eq. (28), Eq. (28) has following expression.

$$\begin{array}{c}\mathit{dV}=V^{\ast }\left\{1-{\left[1-\left(1-n\right)\frac{k_0^{,}}{\beta }{\int }_0^T{T}^x{e}^{-E/\mathit{RT}}\mathit{dT}\right]}^{1/\left(1-n\right)}\right\}f\left(E\right)\mathit{dE},n\ne 1\\ \mathit{dV}=V^{\ast }\left\{1-\mathit{exp}\left[-\frac{k_0^{,}}{\beta }{\int }_0^T{T}^x{e}^{-E/\mathit{RT}}\mathit{dT}\right]\right\}f\left(E\right)\mathit{dE},n=1\hfill \end{array}$$ E31

By integrating Eq. (31), the following expression can be obtained,

$$\begin{array}{c}1-\frac{{\int }_0^V\mathit{dV}}{V^{\ast }}=1-{\int }_0^{\infty }\left\{1-{\left[1-\left(1-n\right)\frac{k_0^{,}}{\beta }{\int }_0^T{T}^x{e}^{-E/\mathit{RT}}\mathit{dT}\right]}^{1/\left(1-n\right)}\right\}f\left(E\right)\mathit{dE},n\ne 1\hfill \\ 1-\frac{{\int }_0^V\mathit{dV}}{V^{\ast }}=1-{\int }_0^{\infty }\left[1-\mathit{exp}\left(-\frac{k_0^{,}}{\beta }{\int }_0^T{T}^x{e}^{-E/\mathit{RT}}\mathit{dT}\right)\right]f\left(E\right)\mathit{dE},n=1\hfill \end{array}$$ E32

Substituting Eq. (30) into Eq. (32), the equation yields,

$$\begin{array}{c}1-\alpha ={\int }_0^{\infty }{\left[1-\left(1-n\right)\frac{k_0^{,}}{\beta }{\int }_0^T{T}^x{e}^{-E/\mathit{RT}}\mathit{dT}\right]}^{1/\left(1-n\right)}f\left(E\right)\mathit{dE},n\ne 1\\ 1-\alpha ={\int }_0^{\infty }\mathit{exp}\left(-\frac{k_0^{,}}{\beta }{\int }_0^T{T}^x{e}^{-E/\mathit{RT}}\mathit{dT}\right)f\left(E\right)\mathit{dE},n=1\hfill \end{array}$$ E33

A double model [15] was developed assuming that the pyrolysis process occurs in two steps with different kinetic behaviors. The pyrolysis process is divided into two steps: the tar and light hydrocarbon gas formation during the primary pyrolysis and the char condensation, cross-linking reactions, and a further gas production during the secondary pyrolysis. Two sets of parallel reactions occur, sharing the same pre-exponential factor but not the same distributed activation energy. The distribution function of activation of energy $f\left(E\right)$ in 2-DAEM equation can be written as,

$$f\left(E\right)=\mathrm{\omega }{f}_1\left(E\right)+\left(1-\mathrm{\omega }\right)f_2\left(E\right)$$ E34

where $\omega$ is a parameter that weighs the different reaction classes, varying from 0 to 1. This parameter describes how many volatiles are released during the various reactions. $f_{\epsilon }\left(E\right)$ , where $\epsilon$ denotes the stage number, is a Gaussian or logistic function of the form with a mean activation energy $E_0$ and a standard deviation $\sigma$ ,

$$\begin{array}{c}f\left(E\right)=\frac{1}{\sigma \sqrt{2\pi }}{\mathrm{e}}^{-\left(E-E_0\right)/2{\sigma }_2}\\ f\left(E\right)=\frac{\pi }{\sqrt{3}\sigma }\frac{e^{\frac{-\pi \left(E-E_0\right)}{\sqrt{3}\sigma }}}{{\left[1+e^{\frac{-\pi \left(E-E_0\right)}{\sqrt{3}\sigma }}\right]}^2}\end{array}$$ E35

In the 1-DAEM case, there are four parameters to be estimated: the mean activation energy $E_0$ , the standard deviation $\sigma$ , the temperature exponent $x$ , and the chemical reaction order $n$ . In the 2-DAEM case, there are seven parameters to be estimated: two mean activation energies $E_{01}$ and $E_{02}$ , two standard deviations ${\sigma }_1$ and ${\sigma }_2$ , the temperature exponent $\gamma$ , the chemical reaction order, and $\omega$ .

In this work, the model-fitting exercise was carried out using MATLAB TM (R2014b). It is essential to use an appropriate method to compute the inner $\mathit{dT}$ integral of Eq. (33). An approximation proposed by Cai and Liu [21] was validated by literature [24], and was applied to estimate the general temperature integral,

$${\int }_0^T{T}^{\gamma }{e}^{\left(-\frac{E}{\mathit{RT}}\right)}\mathit{dT}=\frac{R{T}^{\gamma +2}}{E}{e}^{\left(-\frac{E}{\mathit{RT}}\right)}\frac{0.99954E+\left(0.58058-0.044967\gamma \right)\mathit{RT}}{E+\left(2.54+0.94057\gamma \right)\mathit{RT}}$$ E36

To deal with the outer $\mathit{dE}$ integral of Eq. (33), Simpson’s 1/3 rule [22] was employed. In the numerical integration of Eq. (33), $E_0+3\sigma$ value has been used for the upper limit of the outer $\mathit{dE}$ integral. Therefore, for given values of the parameters, the DAEM equation can be numerically solved.

The kinetic parameters of the DAEM are determined using an objective function. An effective criterion for achieving the objective function is to minimize the differences of squares. The functional form of the objective function used in this paper is,

$$\mathit{\text{OF}}={\sum }_x^l{\left({\alpha }_{x,\mathit{cal}}-{\alpha }_{x,\mathit{exp}}\right)}_i^2$$ E37

where ${\alpha }_{x,\mathit{cal}}$ and ${\alpha }_{x,\mathit{exp}}$ are the calculated and experimental values of the residual mass fraction, $l$ is data number, and the sub-index $x$ makes reference to data point $x$ . In this paper, the objective function ( $\mathit{\text{OF}}$ ) has been optimized using the pattern search method [23], which is a class of direct search methods. The Gaussian and logistic DAEMs show an interrelationship between activation energies and the pre-exponential factor. Thus, $E_0$ obtained for DAEM differed significantly over values of $k_0^{\prime }$ , and it was difficult to figure out their relationship. To remove this ambiguity, the $k_0^{\prime }$ was set to be a fixed value according to the results from isoconversional methods. The initial value of the mean activation energy $E_0$ and standard deviation $\sigma$ settled in pattern search method used the results from isoconversional methods.

5.2. Processing thermogravimetric analysis data by distributed activation energy models

In the distributed activation energy models (DAEMS) adopted here, the pre-exponential factor $k_0^{\prime }$ is shared by all the reactions. From the literature [20], it is known that there exists an interrelation between the pre-exponential factor $k_0^{\prime }$ and the mean activation energy $E_0$ , which means that several pairs of $k_0^{\prime }$ and $E_0$ values can fit equally well for a given set of experimental data. Thereby, the value of $k_0^{\prime }$ has to be fixed to remove the mutual influence between $k_0^{\prime }$ and $E_0$ . Nevertheless, when the value of $k_0^{\prime }$ is selected within two orders of magnitude, the observed correlation between $k_0^{\prime }$ and activation energy does not adversely influence the utility of kinetic models [24]. In this paper, the value of $k_0^{\prime }$ is set to be $1.0\mathrm{E}+10\left(S^{-}\right)$ , which is based on the results from the isoconversional methods discussed above. The initial guess value of $E_0\left(S^{-}\right)$ used in the pattern search method of DAEMs is also based on the results from isoconversional methods of FWO model, as shown in Figure 3, which is corresponding to the peak activation energy and can be set to be $173\left(\mathrm{kJ}\bullet {\mathrm{mol}}^{-1}\right)$ . The expected initial value of $\mathrm{\sigma }$ was set to be $20\left(\mathrm{kJ}\bullet {\mathrm{mol}}^{-1}\right)$ , which is in accordance with the peak width of the $f\left(E\right)$ obtained from isoconversional methods. The value of the objective equation ( $\mathit{\text{OF}}$ ) can be used as a benchmark for evaluating the prediction results obtained by different DAEMs.

The DAEMs were then applied to fit with the experimental data of $\beta =20\left(\mathrm{K}\bullet {\text{min}}^{-1}\right)$ . Figure 5 shows the comparison of the various 1-DAE models. The kinetic parameters of these models and the parameters achieved with these models are listed in Table 5. It can be noticed from Figure 5 that 1-DAE models with Gaussian/logistic distribution do not show good match with experimental values, which suggests that it is not appropriate to describe the pyrolysis process based on the single activation energy distribution for Zhundong coal. 2-DAEMs consider that pyrolysis occurs in two steps assuming that the first and the second stage share the same activation energy. The comparison of the various 2-DAE models with experimental results for Zhundong coal at a heating rate of $\beta =20\left(\mathrm{K}\bullet {\text{min}}^{-1}\right)$ are shown in Figure 6, and the kinetic parameters extracted from these models are listed in Table 6. Looking at the value of $\mathit{\text{OF}}$ listed in Tables 5 and 6, it can be concluded that 2-nth-DAEM assuming the Gaussian distribution gives best fit with the experimental data. From the values extracted with the 2-DAE models assuming the Gaussian distribution, as shown in Table 6, it can be seen that the mean activation energy for the primary pyrolysis (where breakage of weaker bonds takes place) is lower than the mean activation energy for the secondary pyrolysis (where breakage of tougher bonds takes place). However, the difference in activation energy between the two stages is not very large. This behavior is more visible with a view to the plots showed in Figure 7 where the $f\left(E\right)$ function versus activation energy for 2-nth-DAEM with Gaussian distribution is exhibited. It can be concluded that the amount of volatile production attributed to primary pyrolysis stage is a little more than the secondary stage for the Zhundong coal. The original 2-DAEM put forward by Caprariis et al. [15] considers coal pyrolysis as a first order reaction, that is, n = 1, and the pre-exponential factor is a constant. This work takes all these factors into account and will be more conducive to better predictions.

Figure 5.

Comparison of experimental data and simulated data with the 1-DAEM methods: (a) Gaussian distribution; (b) logistic distribution.

		$E_0\left(\mathrm{kJ}\bullet {\mathrm{mol}}^{-1}\right)$	$\sigma \left(\mathrm{kJ}\bullet {\mathrm{mol}}^{-1}\right)$	$\gamma$	n	OF
$n=1$	Gaussian	169.42	24.45	−0.444	1	0.0493
	Logistic	175.11	79.84	−0.400	1	0.2148
$n\ne 1$	Gaussian	172.96	18.73	0.127	5.766	0.0058
	Logistic	172.96	19.35	0.132	5.786	0.0060

Table 5.

Estimated values of kinetic parameters [1-DAEMs with Gaussian/logistic distribution functions].

Figure 6.

Comparison of experimental data and simulated data with the 2-DAEM methods: (a) Gaussian distribution; (b) logistic distribution.

Model		$E_{01}\left(\mathrm{kJ}\bullet {\mathrm{mol}}^{-1}\right)$	${\sigma }_1\left(\mathrm{kJ}\bullet {\mathrm{mol}}^{-1}\right)$	$E_{02}\left(\mathrm{kJ}\bullet {\mathrm{mol}}^{-1}\right)$	${\sigma }_2\left(\mathrm{kJ}\bullet {\mathrm{mol}}^{-1}\right)$	$\mathbf{\gamma }$	n	$\omega$	OF
$n=1$	Gaussian	150.00	20.00	173.00	36.85	−0.500	1	0.400	0.0492
	Logistic	151.02	41.89	171.08	25.06	−0.514	1	0.304	0.0160
$n\ne 1$	Gaussian	150.00	28.22	173.00	26.85	−0.502	2.872	0.646	0.0001
	Logistic	150.00	29.57	181.83	34.50	−0.513	2.826	0.765	0.0005

Table 6.

Estimated values of kinetic parameters [2-DAEMs with Gaussian/logistic distribution functions].

Figure 7.

Distribution activation energy curves as a function of activation energy for 2-nth-DAEM with Gaussian distribution functions.

The kinetic parameters extracted with 2-nth-DAEM with Gaussian distribution functions (best models to describe data at 20 $\left(\mathrm{K}\bullet {\text{min}}^{-1}\right)$ for Zhundong) were used to predict the experimental data obtained at different heating rates of $10\left(\mathrm{K}\bullet {\text{min}}^{-1}\right)$ , $30\left(\mathrm{K}\bullet {\text{min}}^{-1}\right),$ and $50\left(\mathrm{K}\bullet {\text{min}}^{-1}\right)$ , respectively. The results obtained with 2-nth-DAEM are shown in Figure 8. It can be seen that the 2-nth-DAEM with Gaussian distribution shows good match with other three heating rate curves.

Figure 8.

Comparison of predicted results with the experimental data at other heating rates with 2-nth-DAEM with Gaussian distribution functions: (a) $\mathrm{\beta }=10\left(\mathrm{K}\bullet {\text{min}}^{-1}\right)$ ; (b) $\mathrm{\beta }=30\left(\mathrm{K}\bullet {\text{min}}^{-1}\right)$ ; (c) $\mathrm{\beta }=50\left(\mathrm{K}\bullet {\text{min}}^{-1}\right)$ .

6. Conclusions

As a raw material for the preparation of phase change materials, coal is of great significance to energy utilization. In this work, Zhundong coal was taken as an example to investigate the pyrolysis process of it by means of thermogravimetric tests. In order to evolve appropriate kinetic parameters to characterize the pyrolysis process of coal, various methods and potential assumptions for processing TGA data were rigorously analyzed. The four isoconversional methods, that is, Friedman, FWO, KAS, and Miura-Maki method were used to calculate the activation energy and pre-exponential factor of Zhundong coal samples, respectively. It was found that the FWO model gave the best fit for the experimental data. During the pyrolysis process, the distribution of the activation energy of Zhundong coal, $f\left(E\right)$ , shows a unimodal distribution, with a peak at the activation energy of 173 $\left(\mathrm{kJ}\bullet {\mathrm{mol}}^{-1}\right)$ .In addition, Several DAEMs have been theoretically derived, considering the reaction order, the dependence of the pre-exponential factor on temperature, the steps mechanism for the pyrolysis, and the form of activation distribution, wherein the 2-nth-DAEM assuming Gaussians distribution of activation energy gave the best fit for the experimental data of Zhundong coal.

Acknowledgments

This work is currently supported by the National Natural Science Foundation of China through contract No.51606153 and China Postdoctoral Science Foundation (2017M3200).

E	activation energy (kJ·mol−1)
E0	mean activation energy (kJ·mol−1)
fE	distribution curve of the activation energy E
i	the number of constituents in original materials
k	rate coefficient s−1
k0	pre-exponential factor s−1
n	The reaction order
R	universal gas constant (kJ·mol−1·K−1)
R∞	maximum weight loss rate (%·min−1)
T	temperature (K)
T0	initial decomposition temperature (°C)
Tp	peak temperature (°C)
V	total mass of volatiles released at time ‘t’ (%)
Vf	the corresponding weight loss percentage when the pyrolysis reaches the final temperature (%)
V∗	total amount of volatile yields (%)
α	Conversion
β	heating rate (K∙min−1)
γ	Constant
ω	weighing factor
σ	standard deviation (kJ·mol−1)
Subscript
cal	calculated value
exp	experimental value