SYMTERM — program for modelling chemical processes in non-isothermal conditions

Computers and Chemistry 24 (2000) 527 – 532 www.elsevier.com/locate/compchem

SYMTERM — program for modelling chemical processes in non-isothermal conditions K.T. Wojciechowski *, A. Mal*ecki Institute of Inorganic Chemistry, Faculty of Materials Science and Ceramics, Uni6ersity of Mining and Metallurgy, Al. Mickiewicza 30, 30 -059 Krako´w, Poland Accepted 26 November 1999

Abstract The program SYMTERM allows to perform the numerical simulations of chemical processes which take place in the thermal analysis apparatus, considering two important factors: the exchange of heat between a sample and an environment, and the influence of heat evolved in the reaction on its kinetics. The program makes it possible to obtain kinetics profiles (DTA, TG, DTG or DSC curves) at various heating rates and permits to determine both the activation energy Eact and the reaction order of a single rate-limited step reaction with using DTA (DSC) and TG data simultaneously. © 2000 Elsevier Science Ltd. All rights reserved. Keywords: Kinetics; Computer simulation; Non-isothermal conditions; DTA; DSC; TG

1. Introduction The methods commonly used to investigate the reaction kinetics involve a series of experiments under isothermal conditions at different temperatures. However, these methods are laborious, and the convenience of obtaining data from non-isothermal methods, which enable a range of temperatures to be investigated relatively quickly, has aroused considerable interest. Nonisothermal methods are very useful in the investigations of rapid solid-state reactions, in which the beginning of the reaction cannot be controlled easily. Another advantage of these techniques is that temperature range can be covered continuously and that no regions are  5th International Conference Computers in Chemistry 99, Szklarska Poreba, Poland, 1–6 July 1999, Edited by W. Andrzej Sokalski and Morris Krauss. * Corresponding author. Tel.: + 48-12-6173442; fax: + 4812-6172493. E-mail address: [email protected] (K.T. Wojciechowski)

missed, what may occur in the series of isothermal experiments. Except for very simple cases, integration of the system of differential equations describing the non-isothermal chemical process is not straightforward. Hence, there are many assumptions, which have to be made in the mathematical model to simplify the studied problem. There are many of well-known, simple-to-use methods to determine the kinetic parameters of nonisothermal chemical reactions but, due to their limitations, these simplified methods fail, providing physically meaningless results. Therefore, in the last few years, many precise numerical models for selected nonisothermal techniques have been developed e.g.: TG (thermogravimetry) (Singh et al., 1996), DTA (differential thermal analysis) and DSC (differential scanning calorimetry) (Kurajica et al., 1996), TMDSC (temperature modulated differential scanning calorimetry) (Hohne, 1997; Kanari and Ozawa, 1997; Hohne and Shenogina, 1998; Cao, 1999), (Sbirrazzuoli et al., 1995; Sbirrazzuoli, 1996; Sbirrazzuoli et al., 1997), TPD (temperature-programmed desorption) (Yun-Hang et al.,

0097-8485/00/$ - see front matter © 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 0 9 7 - 8 4 8 5 ( 9 9 ) 0 0 0 8 8 - 1

528

K.T. Wojciechowski, A. Mal*ecki / Computers & Chemistry 24 (2000) 527–532

1996). However, only few programs described in the literature enable numerical evaluation of unknown parameters in the kinetic equation of the reaction (e.g. Kurajica et al., 1996). The SYMTERM-program is designed to examine the processes, which take place in the STA (simultaneous thermal analysis) apparatus, where the temperature of the sample and environment changes during the reaction progress. Two types of applications have been considered: the evaluation of kinetic parameters (activation energy Eact, reaction order a and pre-exponential factor k0,) and simulation of kinetic curves TG, DTG, DTA, DSC, which can be obtained from STA apparatus. There are two optimization tools implemented in the program: optimization method based on variations

of the location of DTA or DTG peaks caused by changes of heating rate, and the method using several kinetic profiles obtained at different heating rates. A few extra functions and features such as: an interactive interface, spreadsheet, statistical functions, graphic procedures (and many more), in which the SYMTERM has been equipped, make working with the reaction model and experimental data so much easier.

2. Mathematical model We have proposed the one-dimensional mathematical model to describe chemical and thermal processes which take place in the STA apparatus, giving two

Fig. 1. The flow diagram of the function which calculates the sum of residuals in Eq. (1).

K.T. Wojciechowski, A. Mal*ecki / Computers & Chemistry 24 (2000) 527–532

529

respectively). The coefficients Kr, Ks and Hs can be obtained from other measurements independently. The DTA signal is calculated from the equation: DTA(t) =Ts(t)− Tr(t)

(5)

and TG [%] signal is a linear function of x: TG(t) =

Fig. 2. Exemplary results of Eact and k0, for devitrification of gehlenite, obtained for different guess values of the estimation procedure.

signals: DTA, TG (or derivatives-DDTA and DTG, optional). The following assumptions have been made in the model: 1. the temperature of the environment (heater interior) is a linear function of time 2. the heat transfer from both the sample and the reference to the environment is by conduction only 3. the heat capacities of the samples in the holders are constant and independent of temperature over given temperature range 4. additionally, the assumption has been made that the reaction kinetic is described by the equation in which the right side can be split into two expressions: the Arrhenius term, and the function of x (amount of reactant) and independent of time and temperature parameters a, b…(e.g. reaction orders,…) Eq. (4). The above assumptions can be expressed by a system of differential equations: dTh =b dt

(1)

dTr =Kr(Th(t)−Tr(t)) dt

(2)

dTs dx =Ks(Th(t)−Ts(t))−Hs dt dt

(3)

dx − Eact = − k0 · exp · f(x,a,b..) dt RT

(4)

where: Th, Tr, Ts are the temperatures of the heater, respectively, the reference and the sample respectively, b-heating rate, Kr, Ks are the temperature coefficients of heat transfer for the reference and the sample, Hs is the temperature coefficient of the reaction (Hs =DH/ (mscps +mhcph), where: ms mh, cps cph-mass and thermal capacity coefficients of the sample and the holder,





Ms − Mp Mp · x(t) + 100 Ms Ms

(6)

where: Ms and Mp the molar masses of the substrate and the product, respectively. This model includes two important, but frequently neglected factors, which are: the heat exchange between the sample and the environment (heater interior) and the influence of heat evolved in the reaction on its kinetics.

3. Numerical procedures SYMTERM integrates the system of Eq. (1) – Eq. (4) numerically, by means of the fifth-order Gear’s method for the stiff systems of the equations (Gear, 1971). The Jacobians required for this procedure are calculated by a finite difference method. We have found Gear’s method to have more effectively addressed our problem in terms of time consumption at given level of accuracy and stability than any other tested procedures: modified fourth-order RK (Press et al., 1995), Bulirsch – Stoer’s (Holland and Liapis, 1983) or Adams – Bashforth – Moulton method (Gear, 1971). The Gear’s method may also be used to solve a system of coupled differential and algebraic equations, given the kinetic Eq. (4) in an integrated form. SYMTERM is equipped with optimization tool coupling Monte Carlo search and Simplex method. The Monte Carlo procedure performs a random search within the specified intervals. This is done by generating a number of random points, and selecting the set of points giving the lowest value of the sum of squares of the residuals. This set of points is used as a seed for the Simplex method. The polyhedrical Simplex is a zero-order method, which performs efficient search of minimum without the need of calculating the derivatives. Moreover, the method allows for applying, in our program, constrained optimization, with limitations imposed on the allowed values of optimized variables. In order to achieve the same weight for all the variables, which define the optimization hypersurface, the logarithm of pre-exponential factor log (k0) is being in preference to the k0 parameters itself, and the space of the parameters is scaled. Otherwise, since simplex does not work properly, optimization is much less efficient. Both of the optimization methods call the function calculating the sum of squares of the residuals:

K.T. Wojciechowski, A. Mal*ecki / Computers & Chemistry 24 (2000) 527–532

530 m

S= % (Xi −xi )2

(7)

i=1

where, m-number of points, Xi and xi the experimental and theoretical point, respectively.

4. Determination of the kinetic parameters SYMETERM can estimate kinetic parameters using two alternative methods. The basis of the first method rests on the well-known fact that the location of DTA, DSC or DTG peak depends on heating rate b, activation energy Eact and pre-exponential factor k0. In this method the peak extrema are evaluated for various

heating rates and compared with experimental data. The figure (Fig. 1.) shows the flow diagram of the function, which uses DTA data for optimization. If Hs parameter is B 0 then the reaction is exothermic, and a simple bisection procedure searches for maximum of DTA peak; otherwise, the minimum is searched. The minimized function returns the sum of squares of residuals of temperatures Ts. We have observed that the above method works properly for experimental data with the relatively small deviation errors. When the temperatures of the samples cannot be determined with accuracy DTB 3°C, then the optimization procedure returns significantly different pairs of parameters Eact and k0. The figure (Fig. 2.) presents the exemplary

Fig. 3. The flow diagram of the function which calculates the sum of residuals in Eq. (2).

K.T. Wojciechowski, A. Mal*ecki / Computers & Chemistry 24 (2000) 527–532 Table 1 The parameter values estimated for the Eq. (9) Parameter

Values

Reported values (kJ mol−1)

Eact

160–172 kJ mol−1

k0 a b

5·1010–2.7·1011 5·10−5–3·10−3 0.06–0.13

1779 3a 1889 4a 1869 11b – – –

a b

L’vov and Novichikhin, 1995 Wojciechowski and Malecki, 1999.

Fig. 4. Optimization results for kinetic data of the decomposition reaction of Cd(NO3)2.

531

trema of DTA peaks obtained take very close positions. That means that activation energy Eact and k0 can not be estimated separately with the satisfactory accuracy. This observation explains the significant differences in results obtained from analogous methods, which are based on the effect of the variations of maximum of the DTA peaks caused by changes of the heating rates F (Sharp, 1970). It is known that both parameters Eact and k0 also influence the shape of kinetic curves e.g. DTA, TG, and DTG. Hence, the use of the whole kinetic curves, obtained for different heating rates, should provide independent data, which can be used for more precise estimation of the parameters searched. The second method implements the procedure of fitting the whole theoretical kinetic profiles to experimental data (DTA, DSC and/or TG). The problem, which arises subsequently, is how to make exact calculation of the sample temperature Ts. The reaction kinetics in non-isothermal conditions depends directly on the temperature Ts, rather than on the time t elapsed from the beginning of the process. So, the residuals should be calculated for the same temperature values Ts instead for the time values t. It has been concluded, that due to the fact that numerous factors, such as: small variations of heating rate, the heat transport by convection and radiation, the changes of thermal capacity of the samples, etc., which have been neglected, the kinetic model proposed Eqs. (1) – (4) fails to describe the temperature lines with sufficient precision required for calculations. Unfortunately, even small variations of temperatures Ts may lead to substantial changes in the reaction rate. Hence, the method of adjusting time step h has been applied in the function to bring calculated values of Ts into close agreement with experimental data. The figure (Fig. 3) shows the flow chart for the function called by the optimization procedure. Given the arrays of temperatures Ts[n, m], amount of reactant x[n, m] and time t[n, m] this function returns the sum S of squares of residuals evaluated for m kinetic profiles. The sum S is calculated with consideration to the weight w[m] defined for each kinetic curve individually. The time step h is adjusted to achieve satisfactory of the temperatures matching.

5. An application example Fig. 5. Simulation results of DTA curve obtained for evaluated kinetic parameters.

results for the reaction of devitrification of gehlenite (Malecki et al., 1997). The analysis of the graph shows that the Eact and k0 are strongly correlated. If Eact and k0 increase in such a way that Eact –log (k0) is constant then for various heating rates b the appropriate ex-

The ability of SYMTERM to treat analytical problems was checked on the data obtained from the experiments. As an example of the application the reaction of thermal decomposition of cadmium nitrate Cd(NO3)2 is presented. The reaction mechanism was investigated by means of the TG, DTA and EGA techniques (Wojciechowski and Malecki, 1999). The kinetics of the process was studied using quasi-isother-

532

K.T. Wojciechowski, A. Mal*ecki / Computers & Chemistry 24 (2000) 527–532

mal and non-isothermal methods with various heating rates. The following reaction was proposed for a description of the last step of the decomposition process: Cd(NO3)2(s/l) “CdO(s) +2NO2(g) + 0.5O2(g)

(8)

The reaction can occur in solid or liquid state, depending on the conditions of the experiment. For the reaction in liquid state the general kinetic equation was considered: dx −Eact a = k0 exp( )x (1−x)b dt RT

(9)

For evaluation of unknown parameters Eact, k0, a and b the method of fitting the whole theoretical kinetic profiles to experimental data TG was applied. Up to three normalized TG curves obtained at different heating rates were simultaneously treated in the optimization process. The estimated parameters are collected in the Table 1, and the results of the best fitting are shown in the figure (Fig. 4). The analysis of the results indicates that the activation energy Eact is in good agreement with values reported elsewhere and the reaction was found to be of zero-order. This may suggest that the process of decomposition is limited by diffusion, or the reaction occurs on the surface only. Since the value of the activation energy Eact is too high for the diffusion processes occurring in liquid phase and is suitable rather for the processes limited by reaction; thus, the latter suggestion more likely seems to be the case. To check the correctness of the obtained results the simulation of DTA line was performed with the use of the estimated parameters (Fig. 5). The discrepancies in the lines of the calculated DTA and the experimental data, evident in the initial state of the reaction, are the results of the melting process (i.e. change of the thermal conductance Ks between sample and holder), which has not been included in the kinetic model.

6. Conclusions This article has described the new program, which allows the simulations of processes taking place in the STA apparatus. The program permits the evaluation of unknown parameters of the kinetic equations for the mechanisms studied and may be applied as a useful tool in the kinetic studies of the solid- or liquid-state reactions, having the kinetic data available from the nonisothermal measurements.

References Cao, J., 1999. Mathematical studies of modulated differential calorimetry-II: kinetic and non-kinetic components. Thermochim. Acta 329, 89. .

Gear, C.W., 1971. Numerical initial value problems in ordinary differential equations. Prentice-Hall, Englewood Cliffs, New Jersey. Holland, C.D., Liapis, A.I., 1983. Computer Methods for Solving Dynamic Separations Problems. McGraw-Hill Book Company, New York. Hohne, G.W.H., 1997. Evaluation of temperature-modulated differential scanning calorimetric measurements by phasesensitive rectification techniques. Thermochim. Acta 304305, 209. Hohne, G.W.H., Shenogina, N.B., 1998. Finite-element calculations on the behavior of a DSC in temperature-modulated mode. Thermochim. Acta 310, 47. Kanari, K., Ozawa, T., 1997. Simulation of temperature modulated DSC of temperature dependent heat capacity. Thermochim. Acta 304-305, 201. Kurajica, S., Bezjak, A., Tkalcec, E., 1996. Resolution of overlapping peaks and the determination of kinetic parameters for the crystallization of multicomponent system from DTA or DSC curves: I: non-isothermal kinetics. Thermochim. Acta 288, 123. L’vov, B.V., Novichikhin, A.V., 1995. Mechanism of thermal decomposition of anhydrous metal nitrates. Spec. Acta B 50, 1427. Malecki, A., Gajerski, R., Labus, S., Prochowska-Klisch, B., Oblakowski, J., 1997. Kinetic and mechanism of crystallization of gehlenite glass pure and doped with Co2 + , Eu2 + , Cr3 + and Th4 + . J. Non-Cryst. Sol. 212, 55. Press, W.H., Teukolsky S.A., Vetterling, W.T., Flannery B.P., 1995. Numerical recipes in C, Cambridge University Press, USA. Sharp, J.H., 1970. Reaction kinetics. In: Mackenzie, R.C. (Ed.), Differential Thermal Analysis, vol. 2. Academic Press, London. Sbirrazzuoli, N., Girault, Y., Elegant, L., 1995. Simulations for evaluation of kinetic methods in differential scanning calorimetry: part 1-application to single-peak methods: Freeman – Carroll, Ellerstein, Achar-13 Brindley – Sharp and multiple linear regression methods. Thermochim. Acta 260, 147. Sbirrazzuoli, N., 1996. Simulations for evaluation of kinetic methods in differential scanning calorimetry: part 2-effect of additional noise on single-peak methods. Thermochim. Acta 273, 169. Sbirrazzuoli, N., Girault, Y., Elegant, L., 1997. Simulations for evaluation of kinetic methods in differential scanning calorimetry: part 3-peak maximum evolution methods and isoconversional methods. Thermochim. Acta 293, 25. Singh, S.D., Mazumdar, P.S., Devi, W.G., 1996. A critical appraisal of the Dixit and Ray method for the analysis of non-isothermal kinetic data. Thermochim. Acta 284, 389. Wojciechowski, K.T., Malecki, A., 1999. Mechanism of thermal decomposition of cadmium nitrate Cd(NO3)2·4H2O. Thermochim. Acta, 331, 73. Yun-Hang, H., Hui-Lin, W., Khi-Rui, T., Chak-Tong, A., 1996. Computer simulation of derivative TPD. Thermochim. Acta 274, 289.