Pyrolysis kinetics of a Wangqing oil shale using thermogravimetric analysis

A thermogravimetric analyzer was used to study the pyrolysis characteristics of an oil shale at different heating rates. Friedman, Flynn–wall–Ozawa (FWO), and distributed activation energy model (DAEM) methods were used to calculate the kinetic parameters of low‐temperature pyrolysis of oil shale. The results show that as the heating rate increases, (a) the initial pyrolysis temperature, peak temperature, and final pyrolysis temperature of oil shale shift toward the high‐temperature region, (b) the range of devolatilization temperature increases, and (c) the peak value of oil shale weight loss increases. The apparent activation energies calculated using the FWO and Friedman methods change with the conversion rate, and the average activation energies are 233.71 and 234.61 kJ/mol, respectively. Kinetic parameters obtained from the DAEM calculation of low‐temperature pyrolysis of oil shale have energy compensation effects. By combining F‐test and quantitative criterion test, it was found that the logistic DAEM model is suitable for solving oil shale kinetic parameters. The activation energy of low‐temperature pyrolysis of oil shale is 238.67 kJ/mol, and the reaction order n is 1.


| INTRODUCTION
In recent years, with the rapid economic development and upgrade of industrial structure, China's oil consumption has been increasing year by year, and the energy demand has also been increasing day by day. Conventional energy resources such as coal, oil, and natural gas cannot satisfy the needs of economic development. Therefore, new energy resources are needed. 1 Oil shale, also known as kerogen shale, is a finegrained sedimentary rock containing solid combustible organic matter in the mineral body. Oil shale is a high-mineral, saprolite-like coal, and solid fossil fuel with a low calorific value. 2 The native organic matter present in oil shale is insoluble in petroleum solvents, and a large amount of shale oil can be extracted during the pyrolysis. The amount of oil shale resources available in the world is 7.19937 × 10 12 t, equivalent to 4.7644 × 10 12 t of shale oil. 3 It is expected to become an alternative energy source for petroleum. China has 21.367 billion tons of proven oil shale resources, and China's shale oil production is increasing year by year. Shale oil production was ~700 000 tons in 2012, 750 000 tons in 2013, and reached 800 000 tons in 2014, ranking the first in the world in terms of oil production. Therefore, oil shale with abundant reserves is an alternative energy source for oil. 4,5 Nowadays, thermal analysis of fuels such as biomass, coal, oil shale, and oil sand is used as a method for | 913 HUA et Al.
analyzing the thermal kinetics samples during pyrolysis. [6][7][8][9][10][11][12][13] Thermogravimetric analysis can accurately determine the pyrolysis properties and kinetic parameters of fuels; therefore, this method is widely used for this purpose. Research on the pyrolysis characteristics and kinetics of oil shale is the basis of oil shale retorting and oil refining. This is important for the efficient conversion and utilization of oil shale and energy conservation and carbonization reduction. Oil shale is complex, diverse, and inhomogeneous. Studies on oil shale kinetics are the premise and basis to ensure the optimal transformation of oil shale. Solid-state reaction kinetics is the core of thermal kinetics analysis, and its main task is to determine the "three kinetic factors" of reaction, that is, apparent activation energy, frequency, and reaction mechanism. The performance of oil shale pyrolysis is mainly dependent on the properties of oil shale itself, heating rate, particle size, pyrolysis temperature, and other factors. To improve the yield and quality of shale oil, Researchers studied the effect of final pyrolysis temperature, 14 In this study, thermogravimetry and differential thermogravimetry (TG-DTG) method was used to study the effect of different heating rates on the pyrolysis characteristics of oil shale. Friedman and FWO methods were used to calculate the apparent activation energy of low-temperature pyrolysis of oil shale. Various probability density functions (PDFs) such as Gamma, Logistic, Lognormal, Gaussian, Rayleigh, and Weibull have been proposed, and the n th DAEM model is used to solve the kinetic parameters of oil sand. All these provide a theoretical basis for the development and utilization of oil shale.

| Experimental samples
A typical Chinese oil shale obtained from Wangqing mine in Jilin Province was selected for this study. The collected sample was crushed, ground, and screened to a particle size of up to 74 μm and labeled as WQ. The proximate and ultimate analyses of oil shale samples satisfied the national standards GB/T212-2011 and GB/T476-2001 of China, and the Fisher assay satisfied the national standard SH/T0508-92 of China. All the samples were stored in a dryer. The results of proximate and ultimate analyses of Wangqing oil shale are shown in Table 1. The results of proximate analysis indicate that Wangqing oil shale has the typical characteristics of a low content of fixed carbon (FC) and a high content of ash (A). The results of Fischer assay analysis of Wangqing oil shale are shown in Table 2. Wangqing oil shale has 5.2% mass percentages of shale oil.

| Thermogravimetric analysis
Thermogravimetric analysis of Wangqing oil shale was carried out using a Mettler-Toledo TGA/DSC1 thermogravimetric analyzer (Switzerland). During the pyrolysis, the mass of each sample was 10 ± 0.1 mg, and the heating rates were set at 5, 10, 20, 30, and 40°C/min. The test temperature was in the range of 50-900°C, and the flow rate of high-purity N 2 was 50 mL/min. Positive pressure conditions were established. A thermogravimetric analyzer continuously measures the weight of sample as a function of temperature, and the rate of weight loss can be obtained based on this measurement.

| Isoconversional method
Isoconversional methods are sometimes referred to as "modelless methods." Because they do not involve the reaction mechanism function, in a constant conversion rate, the reaction rate is only the function of temperature. 23 In the specified conversion with different heating rates, the apparent activation energy was calculated using the isoconversional method. The isoconversion principle lays the foundation for isoconversional methods, which are generally divided into two methods: differential and integral. 24 The two isoconversional methods are shown in Table 3. The apparent activation energy E and frequency factor A can be estimated from the least-squares linear fitting. Fifteen values were taken between the conversion rates α = 0.01-0.99. The Friedman method uses ln(dα/dt) as the ordinate and 1/1000 as the abscissa to obtain the slope E α /R by linear fitting. The FWO method uses ln(β) as the ordinate and 1/T as the abscissa to obtain the slope E α /R by linear fitting. The change in final apparent activation energy (E α ) with conversion rate is calculated from the slope E α /R. Isoconversional kinetic method are sometimes referred to as "modeless methods" because they do not involve reaction mechanism functions, but they are not freely assumed. The basic assumption is that the same reaction occurs in the same proportion independent of temperature at any given degree of reaction. 15 Apparent activation energies were calculated using the equal conversion rate method during the specified conversion at different heating rates. The principle of equal transformation lays the foundation for equal transformation method. In this paper, Friedman 17 and FWO method 18,19 were selected. Friedman method is a differential form, while FWO method is an integral form developed independently. The formulas of two equal conversion methods are shown in Table 1. Apparent activation energy E and frequency factor A can be estimated from the linear fitting of least-square method. Fifteen points were selected between conversion α = 0.01-0.99. The Friedman method used ln(dα/dt) as the ordinate and 1/T as the abscissa in the linear fitting to obtain the slope E α /R, while the FWO method used ln(β) as the ordinate and 1/T as the abscissa in the linear fitting to obtain the slope E α /R. The final apparent activation energy E α was calculated from slope E α /R with the change in conversion.

| n th DAEM
Least-squares method can be used to process DAEM calculations. Cai and Liu 28 reported that the effect of DAEM fitting for the reaction order (n ≠ 1) is much better than n = 1. The equation for distributed activation energy of n-order reaction obtained by assuming Different PDFs are shown in Table 4. The unknown kinetic parameters in DAEM were calculated using the most probable mechanism function. Nonlinear regression 29 involves the minimization of differences between calculated and experimental values, usually expressed as the sum of squared residuals: where e and c are the experimental and calculated values, respectively. Because the single use of RSS does not give the best fit results, the F-test was considered 30,31 : Ni, Nd, and P represent the total number of heating rates, number of experimental points selected, and number of kinetic parameters solved, respectively. Figure 1 shows the TG-DTG curves of oil shale samples at different heating rates. As shown in Figure 1A, oil shale has very similar pyrolysis trends, and the pyrolysis of oil shale can be roughly divided into the following three different stages: The first stage involves low-temperature weight loss from room temperature to ~200°C. The weight loss at this stage is mainly due to the removal of water, including the internal and interlayer water of clay minerals. This is also accompanied by the cracking of a certain amount of mineral such as sodium bicarbonate. The second stage involves the decomposition of organic matter at ~200-600°C, mainly

Methods Equations
Friedman 23 hydrocarbon decomposition accompanied by the escape of gas and oil vapor. The pyrolysis in this section is the main part of oil shale pyrolysis. On the other hand, the softening and intermolecular reorganization of kerogen occurring between 200 and 350°C are the physical changes before the cracking of kerogen into tar. The removal of volatiles mainly occurs in the range of 350-600°C. The macromolecular organics and some inorganic minerals in oil shale undergo chain cleavage, polycondensation, and decomposition reactions with increasing temperature. The third stage above 600°C is generally attributed to the cracking of carbonates such as calcite, dolomite, and iron-rich dolomite. 32 As the heating rate increases, the initial temperature of oil shale pyrolysis, peak temperature, and final pyrolysis temperature shift toward the high-temperature region. The devolatilization temperature range increases, and the peak value of oil shale weight loss increases. As the heating rate increases, more heat is supplied to a certain temperature point. Thus, the degree of decomposition of oil shale increases at this temperature point, and the peak value of weight loss increases. Figure 1B shows that as the heating rate increases, the decomposition temperature of oil shale shifts to the high-temperature region. This is mainly because when the heating rate is high, it is too late to transmit the heat flow to the inside of oil shale particles. The temperature difference between the surface and inside of particle is large, and the internal pyrolysis reaction occurs later, affecting the "thermal hysteresis." The higher the heating rate, the higher the temperature gradient and the promotion of movement of pyrolysis of organic matter to the high-temperature region, consistent with the literature. 33,34 DTG peaks of the second phase are significantly higher than those of the first and third phases, indicating that the second-phase reaction is more intense than the other phases. 35

| Isoconversional kinetics
During the low-temperature pyrolysis of oil shale, the conversion rate showed 15 values between α = 0.01-0.99. Figure 2A shows the kinetic parameters obtained from the Friedman method when the oil shale is heated at different heating rates of 5, 10, 20, 30, and 40°C/min. The relationship between reaction and reciprocal of temperature under different heating rates was analyzed for the low-temperature   (Af(α)). Thus, the apparent activation energy of oil shale pyrolysis was obtained; the average activation energy was 234.61 kJ/mol. When α is 0.1, the linear fit is 0.95 (<0.99), and the activation energy is also lower because of the physical change before the cracking of kerogen into tar. The linear fit between α = 0.02-0.99 works well, up to 0.99. Figure 2B shows that the activation energy of oil shale pyrolysis obtained by Friedman method shows instantaneity as the reaction progresses, and the activation energy is not constant as the conversion rate fluctuates. The activation energy of oil shale had almost no change at a conversion rate of 0.1-0.5 with a short-term surge between 0.5-0.68 and a decrease between 0.68-1.0. The apparent activation energy of oil shale volatilization precipitation stage ranged from 213.71 to 285.59 kJ/mol; a shoulder peak was observed in the entire activation energy with the change in conversion rate, consistent with the results of Klomp and Wright. 36 Figure 3A shows the kinetic analysis diagram of FWO equation when the relationship between heating rate and temperature is linearly fitted to oil shale pyrolysis using the least-squares method. This effectively avoids the reaction mechanism constant and directly calculates the activation energy. Using 15 data between α = 0.01-0.99 and 1/T as the abscissa and ln(β) as the ordinate, 15 straight lines were obtained by linear fitting at a heating rate of 5, 10, 20, 30, and 40°C/min. The slope of the obtained straight line is −1.0516 E/RT, and the intercept is ln(0.00484AE/RG(α)). The trend of activation energy of oil shale with conversion rate was obtained. The apparent activation energy of oil shale volatilization precipitation stage calculated using the FWO method was 192.55-279.88 kJ/mol, and the average activation energy was 233.71 kJ/mol, lower than that calculated using the Friedman method. Figure 3B shows the change in activation energy of lowtemperature pyrolysis of oil shale with conversion based on the FWO method. The trend of activation energy and correlation coefficient with conversion rate calculated using the FWO method was almost the same as that calculated using the Friedman method. The change in activation energy between α = 0.1-0.52 is relatively stable. When the conversion rate is between 0.5-1, shoulder peaks appeared during the change in activation energy. This is mainly because this stage is the hydrocarbon generation stage of oil shale pyrolysis formed by the superposition of two peaks. The correlation coefficient of oil shale pyrolysis calculated using the FWO method is higher, and the activation energy is slightly lower than the Friedman method, indicating that the FWO and Friedman methods can be used to calculate the kinetic parameters of oil shale pyrolysis.

| n th DAEM methods
To determine the kinetic parameters of oil shale volatilization precipitation stage, different PDFs were imported into the MATLAB software for damping the least-squares algorithm. 37 The number of iterations was set to 1500, and the maximum residual of the objective function is 1E−8. DAEM's calculations use a direct search algorithm in Matlab ® (R2014B version), which is coded by selecting 100 datapoints for each heating rate. Different DAEM models corresponding to different PDFs in oil shale pyrolysis are shown in Table 1. The kinetic parameters obtained from the isoconversion model range from 0 to 300 kJ/mol, including the maximum activation energy calculated using the Friedman method. Lobatto orthogonal numerical integration method was used to calculate the nonpositive integers of f (E) and guarantees different distributed fitting results (>99%). Using the average value of Eα (Friedman method) as E 0 , the value of unknown parameter was estimated using the PDF 2 corresponding to f(E). F-test was used to determine the rigor of calculation results. The fitting value r > 0.99 is considered as a good fit between the experimental and calculated results. Figure 4 shows a good fit between the experimental and calculation results of oil shale pyrolysis for DAEM models corresponding to different PDFs. The results of kinetic parameters calculated using different DAEM models are shown in Table 4. The calculated kinetic parameters are completely different. E (Logistic) > E (Gaussian) > E (Weibull) > E (Rayleigh) > E (Gamma) > E (Lognormal). The DAEM model was used to calculate the weight loss data of rice husks. The probability function of different densities fitted the result with r > 0.99; the standard residual ε was <0.48. Therefore, the Gamma, Logistic, Lognormal, Gaussian, Rayleigh, and Weibull DAEM models can be used to calculate the kinetic parameters of low-temperature pyrolysis of oil shale.
Although there are differences in different distributions, Figure 5A shows that the PDF of different DAEM models and their heights are almost unequal. The highest reaction order and height of f (E) and average reaction order of different DAEMs was n = 1.17. Figure 5B shows that the kinetic parameters calculated using the DAEMs corresponding to different PDFs have energy compensation effects, satisfying that lnA = 1.6092 × 10 −4 E − 4.33037, and the correlation coefficient r = 0.99368.
Regarding the fitting properties of different distributions, the correlation r (Table 4) was found to be very similar. By increasing the number of unknown parameters in the model, the Akaike information criteria (AIC) 38 index was used to compare different models with different parameters to determine the optimal distribution of oil shale devolatilization.