An alternative approach to kinetic analysis of temperature-programmed reaction data.

To date, kinetic computations have been carried out efficiently for a great variety of physico-chemical processes including crystallization, melting and solid-solid transitions. However, appropriate methods for the kinetic analysis of chemical reactions, especially multi-staged reactions, are currently lacking. Here we report on an alternative way of treating temperature-programmed reaction data using the reduction of iron(iii) oxide as an example. The main principle in the suggested approach is to take into account every stage of the studied process, resulting in a system of kinetic differential equations. Kinetic parameters (activation energy and preexponential factors) are optimized for each of the stages, and cubic splines are used to approximate the conversion functions that reflect changes in reaction-specific surface area throughout the process. The applicability of the suggested method has been tested on temperature-programmed reduction (TPR) data for iron(iii) oxide samples produced from the original Fe2O3 powder by annealing it at 600, 700 and 800 °C. Results of kinetic analysis obtained at different temperature regimes demonstrate the good stability and performance of the method. Peculiarities of iron(iii) oxide reduction have been revealed, depending on the stage and heating rate. The influence of material morphology on the reduction kinetics has been assessed by comparing preexponential factors corresponding to the first reduction stage. This approach allows a comparison of the structural characteristics of the materials based on the kinetic analysis of the TPR data. Using optimized conversion functions, the initial particle size distribution has been reproduced. Theoretically found particle size distribution was found to correlate well with the experimental distribution obtained via laser diffraction. This journal is © The Royal Society of Chemistry.

Introduction

Experimental techniques for thermal analysis are widely used in the characterization of solids, enabling one to trace processes occurring during the thermal treatment of the material. Investigations of thermic decomposition, phase transformations and heterogeneous chemical reactions are conducted most effectively by non-isothermal methods.[1-3] Information about the dynamics of such processes allows the targeted selection of conditions for the fabrication of more effective functional materials. Methods of temperature-programmed reactions are convenient instruments in a researcher's toolbox for collecting such information.[4-6] The main advantages of this group of methods include the simplicity and accuracy of experimental apparatus and procedure: the gas mixture, consisting of reactive and inert gases, flows through the investigated sample under heating at a constant rate. The concentration of the reactive gas is measured during the experiment at the outlet. Most frequently used variations of this group of methods are temperature-programmed reduction and oxidation (TPR and TPO), deriving from the importance of the reduction or oxidation step in the catalyst production cycle. A remarkable feature of the TPR method is the possibility of comparing the dispersion of oxides prepared by different methods.[7] In some cases, it is possible to determine the quantity of the oxide phase in the investigated sample.[8] According to the location of the peak maximum in the TPR curve, one can compare different samples, revealing strong metal–support interactions and the influence of the support on the catalyst operation.[9-12]

Despite the wide scope of TPR applications in investigational practice, analysis of experimental results is carried out far more superficially. As usual, conclusions are based on such facile assessments as analysis of peak shape alterations or changes in the location of the peak maxima of the TPR curve.[13-15] However, more detailed information about the underlying processes of temperature-programmed reactions can be obtained by a quantitative description of the TPR processes conducted via kinetic analysis.[16]

There are two main groups of method used to process non-isothermal experimental data: isoconversional and model fitting methods. These methods are commonly based on eqn (1):where α is the conversion degree, β is the heating rate (kelvin per second), T is the absolute temperature in kelvin that changes via T = T0 + βτ, τ is time in seconds; A is a preexponential factor (s−1), Ea is activation energy and R is the gas constant. Thus, the temperature function of the rate constant is described by Arrhenius law and the conversion function, f(α), which is frequently associated with the mechanism of heterogeneous reaction or kinetic model, defining the dependence of the reaction rate on the degree of conversion.

The basic assumption of isoconversional (model-free) methods is that the reaction rate depends only on temperature at the same degree of conversion. In order to implement this assumption, it is necessary to have experimental data recorded at several heating rates. Processing non-isothermal experimental data is carried out by means of various model-free methods.[17] In the method suggested by Friedman,[18] the logarithmic representation of rate eqn (1) is used:

Then, for a chosen α, the plot ln(βdα/dT) vs. T−1 obtained from the TPR curves recorded at different heating rates is built up. If it is a straight line, its slope allows evaluation of the activation energy.

The Flynn–Wall–Osawa method,[19] based on the integral form of eqn (1), includes Doyle's approximation for temperature integral:[20]

The procedure for evaluating activation energies is similar to the procedure for the Friedman method.

If the process is characterized by one rate-limiting step or several separated in time, then it is possible to obtain a unique set of kinetic parameters. However, in the case of a blend of oxides, in which reduction is not a separate process, or oxide reacting through multiple intermediate stages (e.g. reduction of iron oxide), kinetic analysis performed by isoconversional methods gives us variable (apparent) Ea, which depends on α. Although it is useful to use isoconversional kinetics to reveal the process mechanism, evaluation of the apparent activation energy and conversion function is complicated due to the absence of sound theory behind chemical solid–gas reactions. Even though there are a lot of mathematical expressions for f(α), each of which refers to a certain type of topochemical mechanism, it is difficult to choose between them due to the fact there can be several mechanisms acting when the wide temperature region is studied. Also, several kinetic models can fit well and be equally accurate but correspond to different activation energies and preexponential factors. Further, in the case of the multistage process mentioned above, each stage can be characterized by an individual kinetic triplet (Ea, A, f(α)), which cannot be addressed by eqn (1).

In the case of a multistage process, it is advisable to implement methods of model fitting which imply a numerical solution to the system of differential equations, each of which corresponds to one of the stages.[21] Thus, the TPR curve is simulated while the kinetic parameters are adjusted to fit the experimental curve. However, because of the correlation between kinetic parameters, results of model fitting strongly depend on the initial approximation.[22,23] Consequently, if no additional assumptions about the process are made, it is impossible to achieve the unique set of kinetic parameters via model fitting methods.

Here we present an attempt to overcome the aforementioned disadvantages of existing methods of kinetic analysis. The aim of this study was to create a universal approach for the kinetic analysis of non-isothermal solid–gas reaction experimental data, based on minimal model assumptions. The key point of the suggested method is to reproduce functions of “active” specific surface area S(α) (where α is conversion degree) from the set of TPR experimental curves. The suggested method allows separating temperature (Arrhenius) dependence and temperature independent (specific surface area changes) parts of the TPR curve and, therefore, obtaining the possibility of determining stable (mathematically) activation energies of separate stages in a multistage TPR process. TPR experimental data of iron(iii) oxide were selected to test the suggested method.

Experimental

Materials and methods

Iron(iii) oxide (99.98%) powder was purchased from Sigma-Aldrich. To prove the influence of sample morphology on the results of the kinetic analysis, a comparison of several samples with different morphologies and surface areas was conducted. Thus, the original Fe2O3 sample was annealed at 600, 700 and 800 °C for 3 hours.

The Brunauer–Emmett–Teller (BET) surface area of the iron oxides was determined on a Autosorb IQ (Quantachrome, USA) device. Scanning electron microscope (SEM) images were taken using a Carl Zeiss Crossbeam 1540xb (Germany) instrument. Particle size distributions were obtained using Fritsch particle sizer Analysette 22 (Germany). TPR measurements were carried out on an automated chemisorption analyzer ChemBET Pulsar TPR/TPD (Quantachrome Instr., USA). Sample powders (35–40 mg) were loaded into the quartz sample cell. To separate water forming during oxide reduction, a liquid nitrogen cold trap was used. Prior to the experiment, samples were annealed at 350 °C in a nitrogen flow for 30 minutes to degas and remove moisture. TPR curves were recorded under different temperature programs, at a flow rate of 50 mL min−1, under flow of 6% H2 + N2 and ambient pressure. Purity of gases was 99.995 vol%.

Theory

Model description

The main disadvantages of existing methods of non-isothermal kinetic analysis include low versatility with respect to complex systems, whose reaction kinetics is not fully described by eqn (1) and low reproducibility due to the strong influence of the initial approximation on the final results. In the present work, we suggest a universal approach to the kinetic analysis of heterogeneous reaction data, obtained in non-isothermal conditions. The main principle is based on the reproduction of relative specific surface functions S(α) from the set of experimental TPR curves. Obtained kinetic parameters are used to characterize investigated materials.

In TPR experiments, small amounts of solid material and high gas flow rates are usually used. In connection to this, the following assumptions are made: (1) because reaction by-products (in the case of TPR–water) are removed from the reaction zone with gas flow, reverse reaction can be neglected; (2) curve broadening caused by longitudinal diffusion is not to be taken into account; (3) the elementary stage of the reduction reaction is first order with respect to hydrogen.

Let's consider the kinetic model of the reduction of bivalent metal with hydrogen that corresponds to the scheme MeO + H2 = Me + H2O. Hydrogen consumption over time, dτ, is defined by eqn (4):where k(T) is the rate constant of the reduction reaction at temperature T; S(τ), s(τ) are absolute and specific “active” surface areas of the metal oxide at time τ; n(τ) is the molar amount of the substance at time τ; Pg is the hydrogen partial pressure in the system. Conversion degree α is connected with molar content via the expression:

Substituting (5) into (4) and taking into account that α is a single-valued function of τ we get the following equation:

Comparing eqn (6) and (1), and assuming k(T) = A* exp(−Ea/RT), we find that:

Thus, the conversion function in that case is connected with the function of specific surface area change; i.e. if f(α) is known, then the s(α) can be evaluated and vice versa. At the present time, there are plenty of mathematical forms to express f(α).[24] The choice of a certain type of conversion function is based either on a priori knowledge about the geometry of the interfacial reacting surface or just on the lowest residual dispersion between the theoretical and experimental curves. However, in both cases the chosen conversion function far from always reflects the real surface geometry of the investigated system: a priori considerations can be wrong; function choice based on minimal residual dispersion also can be wrong if none of the test functions describe the surface geometry. In connection with this, we suggest another approach to the choice of optimal conversion function. The idea is simple: instead of determination of what known function exactly corresponds to the investigated system, we evaluate the optimal function of the surface area change, corresponding to the investigated TPR curves. Further, the realization of the suggested approach will be given taking into account experimental parameters (gas flow rate, hydrogen concentration in the gaseous mixture, heating rate etc.).

Calculations

Consider the kinetics of an n-staged reduction of some oxide (Ox0) by hydrogen to metal state (Me) in the process of TPR:where k is the rate constant of the ith stage of reduction; is the apparent preexponential factor (preexponential factor multiplied by the specific surface area (S0,) of the oxide at zero conversion); E is the activation energy of the ith stage.

Within the small time interval that corresponds to the time needed to fill the tube with the gaseous mixture, when the temperature change and convective transfer can be neglected, the process is described by the following system of differential equations:where α is the degree of conversion of ith oxide; s(α) is the function representing the dependence of the specific relative surface area of the ith oxide from its degree of conversion; V is the reactor volume, containing reduced oxide. In order to solve this system of equations within the certain time interval, one needs to know n pairs of values of , E and n functions of s(α). Let all the s(α) functions belong to one class and differ by only some vector of values λ(y, …, y) of length m. Thus, if we know n pairs of values , E and n*m values of y, modeling of the TPR curve in the process of Ox0 reduction can be carried out in the following way.

Let's introduce the modeling step, which corresponds to the time interval needed to fill the reactor's volume with the gaseous mixture. For a tubular reactor we have:where ωV is the gas flow rate; LMe is the length of the tube segment with sample; Rp is the tube radii; mMeO and ρMeO are the sample mass and density, respectively.

We propose that within the time interval Δτ, the temperature does not change and equals:

After Δτ has passed, the temperature changes abruptly up to the value T = T(τ) = T(τ + Δτ). On every time interval, the system of differential equations is solved numerically, and the change in hydrogen partial pressure and the composition of the reactive mixture are estimated. During the next-time step, the reaction mixture is filled with a new portion of gaseous mixture and the process of modeling is repeated. Because Δτ is small, such a scheme allows the avoidance of a numerical solution of a differential equation of convective transfer and, at the same time, provides sufficient accuracy of the modeling (taken from preliminary numerical experiments). As a result, the model TPR is built:

Because the time step is small, the obtained model curve can be considered as continuous. To solve the inverse problem; i.e. to identify the set of parameters A*,E,λ0,…,λ that allow such model curves to fit the experimental curves (with various heating rates) as best as possible, one has to minimize the following function:where is the hydrogen partial pressure at the jth point of the experimental curve with ith heating rate. In order to solve this problem, one needs to know the certain type of function s(α, λ). All that is known about it is that s(α = 0, λ) ≡ 1.0. Because any other information about this function is absent (other than it should be continuous, single valued and limited in its domain (0.0; 1.0)), in present work we suggest it is sought in the class of cubic splines.[25] At the present time, cubic splines are successfully used for the interpolation and approximation of various functions and experimental data. Cubic spline with m knots is determined by a set of m + 2 parameters; i.e. m values in spline knots and two parameters, defining the type of the cubic spline. Usually, these two additional parameters are values of the second derivative of the spline in the first and in the last knot (if one takes this value to equal 0 it will result in a natural cubic spline). However, in the present work, these two parameters were used to build the cubic spline with the minimum norm of the first derivative, as in ref. 26. Spline with a minimum norm of the first derivative differs from the ordinary one by the absence of oscillations of values in between the knots, that is in agreement with the most general assumptions about the character of dependence of specific surface area from conversion. Finally, in our work to describe s(α, λ), the cubic splines with minimum norm of the first derivative have been used. Also, λ = (y1, …, y) for the spline with m knots; i.e. m−1 value in spline knots, because the value in the first knot (y0) is always 1. To minimize (14) it is better to use methods of nonlinear minimization of zero order that do not require derivative evaluation, because analytical expression for the object function is unknown (the system of differential eqn (10) is solved numerically) and using numerical approximation of derivatives drastically reduces the efficiency of methods similar to Newton's. After a series of numerical experiments in the present work to minimize (14), we used the Covariance Matrix Adaptation Evolution Strategy (CMA-ES) method.[27] This method is of zero order, possesses good convergence and is undemanding in terms of the choice of initial approximation. Preliminary calculations showed that using the object function mentioned above for the minimization allowed no stable solutions. Such behavior is typical for ill-posed problems. In order to solve such problems, regularization methods are usually used.[28] It is known that the formation of the new interface surface requires some energy; namely ΔES = σΔS, where σ is surface tension and ΔS is surface change. Because spontaneous processes progress toward an energy minimum, we presume it is logical to use the value of a general increase of the interfacial surface in the process of oxide reduction as a stabilizer. Therefore, for some oxides:

Let's introduce new function φ(α) as the following:

The following function is used as a stabilizer in the present work:

Finally, the object function (14) is transformed into:where q is a regularization parameter.

After we obtained an optimal set of kinetic parameters and splines of relative specific surface area we tried to evaluate the initial particle size distribution, ω, based on the assumption that initial particle size distribution defines explicitly how relative specific surface area changes throughout the process of TPR. We take into account only the initial iron oxide, Fe2O3, because its reduction is characterized by the process of contraction of the particles during TPR, while the other oxides proceed via the growth and contraction stages that cannot be evaluated and separated precisely. In order to obtain an initial particle size distribution using the optimal s(α) spline and n(α) dependence, we conducted an iterative procedure until the particle size distribution on the nth step reproduces the optimal s(α) spline. The following iterative formula was used:where and ω are the particle size distribution values of the i size interval on the next and current step, respectively; n and Δn are the total mole amount and its change; is the normalization factor; r and Δr are the size and width of the size interval. The full derivation can be found in the ESI.†

In present work, we conducted research on the applicability of the suggested method to describe the kinetics of a heterogeneous chemical reaction on the example of TPR of iron(iii) oxide. We provide the results of calculations and compare them with the corresponding experimental data.

Results and discussion

The reduction process of iron(iii) oxide has been studied in depth[29-34] due to the wide scope of applications of iron and steel in modern technologies. Thus, iron-based catalysts are prospective because of their low price and low methane selectivity in the Fischer–Tropsch process.[12,35,36]

TPR of iron(iii) oxide is very complicated in terms of kinetic analysis because of the multistage character of the reduction (Fe2O3 → Fe3O4 → FeO → Fe).[37] Every stage of this reduction process is characterized by a separate set of kinetic parameters that makes model-free methods, based on eqn (1), inapplicable to the kinetic analysis of such a process. This is also proved by the wide range of activation energies (18–246 kJ mol−1)[32] caused not only by the various compositions and morphologies of the investigated sample, but also by inaccurate kinetic analysis. In our work we used the Friedman method to obtain dependence Ea from α in order to compare it with the corresponding values evaluated via the suggested method. The obtained curves (Fig. 1) clearly demonstrate that the reduction passes through three consecutive steps as there are three regions on the Friedman plots, where Ea values remain more or less constant. The disadvantage of the Friedman method in this case is that Ea cannot be evaluated precisely for each stage due its variation with respect to the studied material or the degree of conversion. In particular, the activation energy of the first reduction step Fe2O3–Fe3O4, proceeding while the conversion degree has not reached 0.11, lies in the range from 7 to 115 kJ mol−1, therefore making evaluation of this parameter inaccurate. Such a wide confidence interval occurs because the reduction of each oxide contributes to the dα/dt term. Because reduction reactions proceed with different rates, the contribution of each reduction reaction to the dα/dt term varies when the heating rate changes, thus widening the interval of the activation energy values. This complication arises from eqn (1), which can be utilized properly only in the case of a single-stage processes. Also, it can be seen that the activation energies of the other reduction stages oscillate from one curve to another, thus suggesting that the single-rate eqn (1) cannot be applied to a multistage process.

Fig. 1

Dependence of activation energy of iron(iii) oxide reduction from conversion degree, calculated via Friedman method. (A) Original Fe2O3 powder; (B) Fe2O3 annealed at 600 °C; (C) Fe2O3 annealed at 700 °C; (D) Fe2O3 annealed at 800 °C.

In our study we used samples of different morphologies to assess the variation in the kinetic analysis results with respect to the studied material. To ensure that all samples have identical chemical composition, we annealed iron(iii) oxide at different temperatures to obtain series of samples with different morphologies.

SEM images of the initial iron oxide powder are shown in Fig. 2. The original powder consists of spherical particles aggregated from flakes of Fe2O3. This sample possesses a highly-developed surface structure, proved by the highest BET surface area among all studied materials (25 m2 g−1; Table 1).

Fig. 2

SEM images of the original Fe2O3 powder.

Optimized values of activation energy and preexponential factor for the reduction of iron oxide(iii) samples and BET surface area values

Sample	E Fe2O3/kJ mol−1	E Fe3O4/kJ mol−1	E FeO/kJ mol−1	A Fe2O3/s−1	S spec./m2 g−1
Fe2O3 original powder	114.1	98.8	72.9	304.9	25.0
Fe2O3 annealed at 600 °C	120.0	111.7	72.2	264.7	4.9
Fe2O3 annealed at 700 °C	111.8	115.9	73.6	44.8	4.4
Fe2O3 annealed at 800 °C	104.3	113.9	83.9	7.9	3.0

Annealing of the iron oxide powder drastically changes the surface morphology of the material. SEM images of annealed samples are given in Fig. 3. Annealed oxides have a labyrinth structure and become more amorphous from the surface with growth of annealing temperature. Increasing annealing temperature leads to particle growth and to a decrease of the inner particle pores that were in original powder (Fig. 2 and 3). At 800 °C, particles tend to sinter into larger grains and the number of surface pores became minimal (Fig. 3C). Specific surface area of the annealed samples decreases with increasing annealing temperature (4.9, 4.4 and 3.0 m2 g−1; Table 1). Although values of the specific surface area of annealed samples are close to each other, such small differences in morphology affects the TPR spectrum and can be revealed by the analysis of optimized kinetic parameters, as discussed below.

Fig. 3

SEM images of the Fe2O3 samples, annealed at various temperature. (A) 600 °C; (B) 700 °C; (C) 800 °C.

TPR curves of iron(iii) oxide samples, recorded at several heating rates (3, 6, 9 and 12 °C min−1), are presented in Fig. 4A. Using data obtained at several temperature regimes improves the reliability of the results of the kinetic analysis according to the International Congress on Thermal Analysis and Calorimetry (ICTAC) project recommendations.[38] Curve fitting was carried out using a three-stage scheme of reduction, with the following kinetic equations:

Fig. 4

TPR curves of the original Fe2O3 powder (A) recorded under various heating rates, and TPR curves of iron(iii) oxide samples (B) recorded under 12 °C min−1 heating rate.

The shape of the reduction curve changes with heating rate, with process shifts to higher temperatures (Fig. 4A). The first reduction peak shifts by 50 °C from 3 °C min−1 curve to 12 °C min−1. There are also variations in TPR curves if we consider them in the row of studied materials (Fig. 4B). The first maximum in the TPR of iron(iii) oxides also shifts to higher temperatures in the progression from original powder to the most annealed. This clearly demonstrates the influence of the surface morphology of the investigated samples on the TPR spectra and, as a result, on the kinetic analysis. However, all TPR curves are approximated well, and the optimization procedure resulted in close values of kinetic parameters demonstrating high reliability of the suggested method (Table 1). It is noteworthy that the preexponential factors found for the first step of reduction correlate well with the specific surface areas of the materials. If we recall that the preexponential factor is a multiplication of the specific surface area at zero conversion and true preexponential factor (9) at a given temperature, then it should correlate to some extent with the specific surface area. Although there is a loose correlation between them, because we deal with the “active” specific surface area that is not the one measured by nitrogen, preexponential factors form the same dependence as the specific surface areas (A1 > A2 > A3 > A4 as well as S1 > S2 > S3 > S4). Thus, we can conclude that the suggested method takes into account the difference in material morphology and is capable of providing reliable kinetic results.

Besides kinetic parameters, we also obtained s(α) curves as presented in Fig. 5. Change in relative specific surface area with degree of conversion for Fe2O3 oxide is a monotone increasing function (Fig. 5A) that can be explained by the contraction particle model. When particles of hematite react with hydrogen, their size decreases and at the same time the specific surface area of the particles increases due to the fact that the ratio between the surface atoms and the bulk atoms grows. The shape of the s(α) curves of the intermediate oxides Fe3O4 and FeO (Fig. 5B and C) is caused by the same factors. The first particles of Fe3O4 formed from the reduced hematite began to grow; thus, their specific surface area decreases. Then, when the reduction rate of Fe3O4 exceeds the growth of Fe3O4, the relative specific surface area begins to grow. It is for this reason that we observe a minimum on the s(α) curve. The same hypothesis we be applied for s(α) dependence of FeO.

Fig. 5

Optimized relative specific surface area vs. conversion degree curves of oxide forms obtained for 12 °C min−1 heating rate. (A) Fe2O3; (B) Fe3O4; (C) FeO.

From s(α) of Fe2O3, we reproduce the initial particle size distribution using the iterative procedure described in the Calculations section. We evaluated particle size distribution and compared it with the one measured experimentally using a laser diffraction method (Fig. 6). The results have shown that the proposed method allows accurate reproduction of the particle sizes, in particular, the mode in particle size distribution based on very simple theoretical assumptions. Deviations of the model from the experimental values come from the fact that the procedure of nonlinear minimization used in the work is very sensitive to the initial values of the optimized parameters. Good consistency between the model and the theory proves that the suggested method for kinetic analysis of TPR satisfies all the requirements and can be used in wider research practice.

Fig. 6

Particle size distribution bar charts of the studied samples obtained via iterative procedure and experimentally measured. (A) Original Fe2O3 powder; (B) Fe2O3 annealed at 600 °C; (C) Fe2O3 annealed at 700 °C; (D) Fe2O3 annealed at 800 °C.

Conclusions

A new approach to the kinetic analysis of the kinetics of heterogeneous reactions in terms of a TPR has been suggested. The method is based on the optimization of activation energy and the preexponential factor separately from the relative surface area function and evaluating this function as cubic spline, which allows non-isothermal data to be described with controlled precision depending on the number of spline knots. Analysis of the optimized relative surface area function enables one to qualitatively assess the influence of temperature and heating rate on the kinetics of the heterogeneous chemical processes.

The presented approach was implemented in the kinetic analysis of the TPR of iron(iii) oxide samples that possess different surface morphology under different heating rates (3–12 °C min−1). Close values of kinetic parameters were obtained for the studied materials by the presented method, proving its reliability. Correlation between the found preexponential factors and the specific surface areas of the iron(iii) oxide samples demonstrates the capability of the method to probe morphology variation among the series of studied materials (from the original Fe2O3 sample to the sample annealed at 800 °C: A* – 304.9 > 264.7 > 44.8 > 7.9; Sspec. – 25.0 > 4.9 > 4.4 > 3.0). Using optimized s(α) dependence, the initial particle size distribution of the sample was obtained and compared with the experimental distribution obtained via laser diffraction. Deviations between model and experimental values can be attributed to the strong dependence of the results from the initial values of the optimized parameters.

Conflicts of interest

There are no conflicts to declare.