Numerical Calculation and 3-D Imaging of the Arrhenius Temperature Integral

: The Arrhenius temperature integral is typically used in non-isothermal kinetic analysis, which is widely applied in gas–solid reactions in separation processes. In previous studies, researchers provided various methods to solve the temperature integral, but the error usually became signiﬁcant when the value of x ( x = E a / RT ) was too large or too small. In this paper, we present a new series method and design a computer program to calculate the temperature integral. According to the precise calculation of the temperature integral, we ﬁrst reveal the relationship among the integral, the temperature, and the activation energy, and we ﬁnd an interesting phenomenon in which the 3-D image of the temperature integral is of self-similarity according to fractal theory. The work is useful for mechanism and theoretical studies of non-isothermal kinetics.


Introduction
Non-isothermal chemical kinetics are very important for gas-solid reactions in separation processes.As a result, they have been extensively introduced in disciplines, such as chemical engineering [1,2], mineral technology [3], metallurgical engineering [4], materials science [5], biomass energy [6,7], etc.The reaction rate of a solid-involved reaction is now generally accepted to be dα dt = k f (α) where α denotes the conversion ratio of the concerned substance at time t, f (α) is the mechanism function, and k is the rate constant expressed by the Arrhenius equation: where A, T, and E a are the pre-exponential factor, temperature, and activation energy of the reaction, respectively.For a linear heating process, the heating rate β is usually a function of time, Combining all of the above equations, re-arranging them, and integrating Equation (1), we obtain where T 0 represents the initial reaction temperature and is set at 298 K in this study, and the initial reaction extent α 0 usually equals zero.
The Arrhenius temperature integral is also known as the integral of Boltzmann factor Ψ(T), which is derived from Equation (4), where another variable x could be introduced as Then, the temperature integral could be calculated as x 2 dx (7) where the upper limit of Equation ( 7) is usually considered to be ∞ regarding T 0 in the denominator, which is equal to 0. Then, the expression of the temperature integral as shown in Equation ( 7) is extensively regarded as the most concise and the best one for calculation.
The solution of the temperature integral is usually the basic demand for a nonisothermal kinetic study [8].However, it cannot be analytically solved.With the solution of the temperature integral, constant activation energy and the corresponding kinetic mechanism could be determined in model-based (single-step) reactions, and variable activation energies could be calculated in model-free iso-conversional (multi-step) reactions.
The model-based non-isothermal analysis was dominant for a long time because of its ability to determine the kinetic triplet (pre-exponential factor, activation energy, and reaction model) [9,10].However, these methods suffer from several problems, among which is their inability to uniquely determine the reaction model.This issue led to the decline of these methods in favor over the latest two decades of iso-conversional methods that evaluate kinetics without model assumptions [11].Regardless of the kind of method, the evaluation of the temperature integral is an inevitable issue.As a result, there is a bottleneck in kinetic analysis, namely the calculation of the temperature integral.
Flynn [25] and Órfão [17] provided excellent reviews of the temperature integral.According to Flynn's categorization, these solutions of the temperature integral are classified into three categories: series solutions, complex approximations [29,30], and simple approximations.It is given accurately for x greater than 2 by the rational fraction expression.Truncating the number of terms in the fraction expression, we can obtain one, two, three, and four "rational approximations".Other series solutions, such as Schlömilch expansion, Bernouli number expansion, and Asymptotic expansion, are all available for the calculation limitation of x in different ranges [13,25].
In model-based non-isothermal analysis, simple and complex approximations were preferred and truncated from those series solutions because the truncated solutions could give the temperature integral an analytical expression and make the whole non-isothermal kinetics analytically solvable.For example, truncating the two-term asymptotic approximation, we can obtain the famous linear plot of ln(g(α)/T 2 ) against 1/T, resulting in a slope of E a /R.In the other words, it is exactly the famous Coats-Redfern method [12].In addition to the truncations, some modified approximations were also proposed to improve the calculation precision [8,[17][18][19][20][26][27][28].
The numerical calculation of the temperature integral is usually carried out by series methods, and model-free iso-conversional non-isothermal kinetics are usually required.Although the precision of the numerical calculation is better today, there are still backward calculations because of the limitations of the series method itself.Numerical calculation based on the quadrature method was occasionally mentioned in previous papers [24,28].However, the residual error usually increases with the rise in the upper limit temperature of the integration, or it increases with the decrease in x [25].As a disquieting situation, some published papers provided the wrong quadrature results as standard integral values, such as in Ref. [24].
There were hundreds of papers focused on the solution of the temperature integral once it was suggested that a relative error of 10% for the temperature integral could be acceptable for many calculations.However, more accurate values are still necessary for scientific purposes, such as distinguishing between reaction mechanisms [25].Furthermore, a more accurate solution possible today with computer technology.All in all, accurate values and numerical calculation methods of the temperature integral are always attractive in non-isothermal kinetic analysis.

Theory
So far, most of the solutions for the temperature integral have been based on the concise Equation ( 7).However, the upper limit of the integration may be far from ∞, resulting in the corresponding inaccurate mathematical equations of series methods.
At the same time, all approximations of the temperature integral take a detour; errors and troubles cannot be eradicated.The only way to obtain the exact numerical solution of the integral is to develop an appropriate series solution, which easily reaches convergence.Here, we provide a new series solution based on Taylor's expansion: where ∆T is the difference between T and T 0 , T = T 0 + ∆T, and n denotes the order of Taylor expansion.Then, the exact derivative value at temperature T 0 could be calculated using Equations ( 7)- (9).
The n-th-order derivative expression shown as Equation ( 9) is established according to Leibniz formula [31].With the help of the general form of the n-th-order derivative equation, we designed a computer program to calculate the integral of the Boltzmann factor with Equation (6).The numerical value can converge within a sufficiently small residual error when the derivative order n is large enough, and ∆T is small enough.Computing the results shows that the value usually converges with reasonably little error when n is greater than 4 and ∆T is no greater than 50 K.Consequently, a multi-step integral method is applied to compute the integral value when the difference between the upper and lower limit temperature is larger than 50 K.The corresponding error analysis of the present calculation method is also provided below.

Error Analysis
Table 1 shows the converging evolution of the integral value with increasing order of Taylor expansion.Most of the values converge within a relative error of 1‰ when the Taylor order is larger than 4. A relative error within 10 −8 can be achieved when the Taylor order increases to 10.To reach high accuracy, a 3D image (Figure 1a) of the integral function is established based on 200 × 200 calculation nodes, and each node is computed with a Taylor order of 20.As shown in Figure 1a-d, the integral values have differences in orders of magnitude from each other with different temperatures and activation energies.To control computing error, the relative residual was defined as follows: where Ψ(n, T, E a ) denotes the integral of Boltzmann factor with upper limit temperature of T, activation energy of E a , and Taylor order of n.The relative residual δ denotes a tolerance when it is presupposed as a maximum limitation that should not be exceeded.Figure 1a-d demonstrates the minimum required Taylor order for calculating the temperature integral under error tolerances of 10 −4 , 10 −8 , 10 −12 , and 10 −16 , respectively.According to Figure 1a,b, the relative residual errors in Figure 1a are no greater than 10 −4 under the domain of definition, and most of the values have a residual error of 10 −8 .The smaller that the minimum required Taylor order is, the more easily that the iteration reaches tolerance.In previous studies, it is simply concluded that the errors become significant when the value of x (x = E a /RT) is too small or too large [25].However, according to these figures, a different trend is discovered.Generally, the trend that is proved by these four figures is that the minimum required Taylor order rises with a decreasing upper limit or an increase in the activation energy when the upper limit temperature is greater than about 350 K, and the minimum required Taylor order rises with increasing the upper limit or increase in the activation energy when the temperature is between 298 and 350 K.
Furthermore, the relationships among T, E a , and Taylor order n are more complicated than monotonic, and some singular points also exist in the plots.However, the maximum required Taylor order is no greater than 45 under a tolerance of 10 −16 , which proves that it is easier to approach convergence than the former methods [25].Another series solution is established to perform a comparison.Using the Newton-Leibniz quadrature method as a usual idea to resolve it, we obtain: where However, we find that the Newton-Leibniz quadrature method was difficult to converge, as shown in Table 2.Moreover, it is not realistic to adopt this method in kinetic analysis because it always requires a massive amount of time to reach high precision.As shown in Table 1, the results converge within 8 accurate digits after the decimal point when the Taylor order n is only larger than 15.Coats and Redfern [12] exp(−x) Gorbachev (1st degree) [33] exp(−x) x 1 x+2

Results and Discussion
The dependence of the temperature integral on E a and T has always been a mysterious relationship.Here, we show the true value of the integral in Figure 2a.This image is drawn based on 200 × 200 calculation nodes, and each node is computed with the Taylor order of 20 to ensure the numerical value is accurate enough.It can be seen that the integral value declines by an order of magnitude with the growth of the activation energy, proving that the integral value is extremely sensitive to the activation energy.Therefore, the approximation method often provides an inaccurate value and fails to distinguish the kinetic mechanisms [26].The integral value tends to be equivalent to the value of ∆T when E a infinitely approaches 0 kJ/mol, while the integral value declines by an order of magnitude and tends to be 0 when E a is greater than about 75 kJ/mol.
There is also an interesting phenomenon in that the profile of the 3D surface is relatively similar to others if the independent variables (temperature, activation energy) are arbitrarily truncated, such as in images (b), (c), and (d).Image (a) is drawn based on a 200 × 200 node network, and each node is computed with the Taylor order of 20 to ensure the numerical value is accurate enough.Independent variables E a and T vary from 2 to 200,002 J/mol and 300 to 1200 K, respectively.In the aforementioned node network (200 × 200), mid nodes (E a varies from 50,002 to 149,002 J/mol, and T varies from 529.5 to 975 K) are chosen to form a new node network (100 × 100); then, image (b) can be drawn.The node network of image (c) is 50 × 50, and its values of E a and T change from 20,002 to 69,002 J/mol and from 300 to 520.5 K, respectively.Regarding image (d), its node network is 49 × 63, and corresponding ranges of E a and T are 52,002 to 100,002 J/mol and 651 to 930 K, respectively.As we can see, the truncated images b, c, and d are similar to the original image a.This phenomenon indicates the self-similarity property of the figure according to fractal theory [38].

Results and Discussion
The dependence of the temperature integral on Ea and T has always been a mysterious relationship.Here, we show the true value of the integral in Figure 2a.This image is drawn based on 200 × 200 calculation nodes, and each node is computed with the Taylor order of 20 to ensure the numerical value is accurate enough.It can be seen that the integral value declines by an order of magnitude with the growth of the activation energy, proving that the integral value is extremely sensitive to the activation energy.Therefore, the approximation method often provides an inaccurate value and fails to distinguish the kinetic mechanisms [26].The integral value tends to be equivalent to the value of ΔT when Ea infinitely approaches 0 kJ/mol, while the integral value declines by an order of magnitude and tends to be 0 when Ea is greater than about 75 kJ/mol.
There is also an interesting phenomenon in that the profile of the 3D surface is relatively similar to others if the independent variables (temperature, activation energy

Conclusions
Because of the demand for accurate solutions related to non-isothermal kinetics, previous researchers have provided various methods to solve the temperature integral, but the error usually becomes significant when the value of x (x = Ea/RT) is too small or too large.A numerical method is proposed to obtain an accurate value of the temperature integral.Multi-step and 20th Taylor expansion are applied to ensure the high accuracy of results calculated by the new method.
The results show that, when the upper limit temperature is greater than about 350 K, the error rises with the decreasing upper limit or the increasing activation energy; when the temperature is between 298 and 350 K, the error rises with the increasing upper limit or the increasing activation energy.
The temperature integral function is very sensitive to the activation energy according to the error analysis and could be well displayed by the 3-D image.The reason why the solution is always divergent when the value of x (x = Ea/RT) is too small or too large, as widely encountered in previous study, can be well interpreted by the contour image of the minimum Taylor order required when using the present method within the iterative tolerances.More interesting, the 3-D image of the temperature integral demonstrates selfsimilarity according to fractal theory.

Resource Codes
Here we share our resource codes for scientific purposes.Anyone who directly uses/translates these codes to produce new papers must cite this paper.All codes below are written in Visual Basic language.
These codes are easily transplanted into a project because they are written in functions.To use these codes, one should first define the environmental temperature (T0, default as 298 K), the presupposed Taylor expansion order (s), the Arrhenius energy (Ea),

Conclusions
Because of the demand for accurate solutions related to non-isothermal kinetics, previous researchers have provided various methods to solve the temperature integral, but the error usually becomes significant when the value of x (x = E a /RT) is too small or too large.A numerical method is proposed to obtain an accurate value of the temperature integral.Multi-step and 20th Taylor expansion are applied to ensure the high accuracy of results calculated by the new method.
The results show that, when the upper limit temperature is greater than about 350 K, the error rises with the decreasing upper limit or the increasing activation energy; when the temperature is between 298 and 350 K, the error rises with the increasing upper limit or the increasing activation energy.
The temperature integral function is very sensitive to the activation energy according to the error analysis and could be well displayed by the 3-D image.The reason why the solution is always divergent when the value of x (x = E a /RT) is too small or too large, as widely encountered in previous study, can be well interpreted by the contour image of the minimum Taylor order required when using the present method within the iterative tolerances.More interesting, the 3-D image of the temperature integral demonstrates self-similarity according to fractal theory.

Resource Codes
Here we share our resource codes for scientific purposes.Anyone who directly uses/translates these codes to produce new papers must cite this paper.All codes below are written in Visual Basic language.
These codes are easily transplanted into a project because they are written in functions.To use these codes, one should first define the environmental temperature (T0, default as 298 K), the presupposed Taylor expansion order (s), the Arrhenius energy (Ea), and the upper limit temperature (T1) of the integral.If the precision worsens, it only requires a smaller iteration temperature (default as 50 K in the code) and a larger Taylor expansion order.Usually, an iteration temperature of 50 K and a Taylor expansion order of 10 are good enough to reach a rational precision.

Standard Integral Values
To be compared with the other values of temperature integrals, the standard integral values calculated by the present method are reported in Table 5
) are arbitrarily truncated, such as in images (b), (c), and (d).Image (a) is drawn based on a 200 × 200 node network, and each node is computed with the Taylor order of 20 to ensure the numerical value is accurate enough.Independent variables Ea and T vary from 2 to 200,002 J/mol and 300 to 1200 K, respectively.In the aforementioned node network (200 × 200), mid nodes (Ea varies from 50,002 to 149,002 J/mol, and T varies from 529.5 to 975 K) are chosen to form a new node network (100 × 100); then, image (b) can be drawn.The node network of image (c) is 50 × 50, and its values of Ea and T change from 20,002 to 69,002 J/mol and from 300 to 520.5 K, respectively.Regarding image (d), its node network is 49 × 63, and corresponding ranges of Ea and T are 52,002 to 100,002 J/mol and 651 to 930 K, respectively.As we can see, the truncated images b, c, and d are similar to the original image a.This phenomenon indicates the self-similarity property of the figure according to fractal theory [38].

Figure 2 .
Figure 2. Three-dimensional image of the integral function of Boltzmann factor.

Figure 2 .
Figure 2. Three-dimensional image of the integral function of Boltzmann factor.

Table 1 .
The integral values of Boltzmann factor with varying Taylor order n.

Table 2 .
Results of the Newton-Leibniz quadrature method.

Table 3
[29]ides specific expressions of the commonly used approximate solutions for the temperature integral.The values of temperature integral calculated by MATLAB were used to check the accuracy of other approximation methods.Table4reports the temperature integral values calculated by different methods.Except for the present solution, the other values are calculated based on a published article[29], in which the original values are equivalent to Ψ/[T•exp(−E a /(RT))].In Table4, the temperature integral values of the present method and Órfão III's approximation are fully consistent with MATLAB's values.Therefore, the new series solution can provide sufficient, accurate values.In contrast, other approximations provide unreliable values when the value of x tends toward minority.

Table 3 .
Expressions of the approximation for temperature integral.