Assessment of isoconversional methods and peak functions for the kinetic analysis of thermogravimetric data and its application to degradation processes of organic phase change materials

Vol.:(0123456789) Journal of Thermal Analysis and Calorimetry (2024) 149:13879–13899 https://doi.org/10.1007/s10973-024-13494-w Assessment ofisoconversional methods andpeak functions forthekinetic analysis ofthermogravimetric data andits application todegradation processes oforganic phase change materials RocíoBayón1 · RedlichGarcía‑Rojas1· EstherRojas1· MargaritaM.Rodríguez‑García1 Received: 17 October 2023 / Accepted: 9 July 2024 / Published online: 3 October 2024 © The Author(s) 2024 Abstract In this work, theoretical kinetic curves of both singleand multi-step reaction mechanisms were simulated by using different sets of kinetic parameters. Various isoconversional methods were applied for the kinetic analysis of these curves so that the corresponding activation energy vs. conversion degree curves were obtained and then compared with the energy values used in the simulations. For single-step reaction mechanisms Friedman method resulted to be the most accurate while for multi-step reaction mechanisms, Kissinger–Akahira–Sunose and Coats–Redfern methods led to the most accurate estimation of the activation energy. On the other hand, conversion rate curves of different single-step reaction mechanisms were fitted with two kinds of peak functions (normalized Fraser–Suzuki and generalized logistic) so that the relationships between the parameters of these functions and the kinetic parameters used in the simulations were obtained. These relationships were then used in the mathematical deconvolution analysis of conversion rate curves simulated for multi-step reaction mechanisms. In general, the curves resulting from deconvolution fitted quite well the simulated conversion rate curves and the analysis of the resulting single-step reaction curves with Kissinger method led the kinetic parameters close to the ones used in the simulations. Finally, a similar kinetic analysis was applied to experimental thermogravimetric measurements taken both under N2 and air for two phase change materials (PCMs) based on polyethylene glycol, PEG6000 and PEG12000. Activation energy values obtained with isoconversional methods for the measurements under N2, varied from 40kJ mol−1 at low conversions up to 150kJ mol−1 at high conversions, whereas for the measurements under air the energy values remained almost constant in the range of 50–75kJ mol−1. The lower activation energies obtained for the measurements under air are clearly associated with the polymer combustion. The experimental conversion rate curves were deconvoluted with the most appropriate peak functions so that the possible single-step reaction mechanisms occurring in these PCMs were separated and further analyzed with Kissinger method. The activation energies obtained with this method were in good agreement with the values resulting from the isoconversional methods. Keywords Kinetic analysis· Isoconversional methods· Fraser–Suzuki function· Generalized logistic function· Fitting· Mathematical deconvolution analysis· TG measurements· PCM Abbreviations A Frequency factor or pre-exponential factor in s−1 E Activation energy in kJ mol−1 CR Coats–Redfern dTG-T curve Conversion rate and also derivative of the thermogravimetric curve dTA Differential thermal analysis FWO Flynn–Wall–Ozawa GLOG Generalized logistic probability density function KAS Kissinger–Akahira–Sunose MDA Mathematical deconvolution analysis NFS Normalized Fraser–Suzuki function PCM Phase change material PEG Polyethylene glycol TG curve Thermogravimetric curve α Conversion degree m ini− m m ini * Rocío Bayón rocio.ba[email protected] 1 Thermal Energy Storage Unit, CIEMAT-PSA, Av. Complutense 40, 28040Madrid, Spain 13880 R.Bayón et al. α-T curve Conversion degree variation with temperature β Heating rate used in the thermogravimetric measurement in K min−1 Introduction Phase change materials (PCMs) with melting temperatures (Tm) between 30 and 160°C are of particular interest for thermal energy storage applications in both lowand midtemperature ranges. In terms of practical implementation, one of the most critical issues when choosing a PCM for a certain application is to be sure that the material keeps its performance along the whole service life of the storage system. In most references found in the literature, PCM long-term performance is assessed through melting/freezing cycles or just by analyzing thermogravimetric (TG) curves (i.e., mass loss over temperature) measured at only one heating rate. However, if the material suffers some kind of degradation after melting, cycles may lead to misleading results because such degradation will be hindered during the freezing period. On the other hand, the information obtained from TG curves in terms of temperature limit for material stability strongly depends on the heating rate used in the experiment. Hence, these tests are not sufficient for validating either the stability nor the successful life performance of a PCM. Therefore, for assessing the long-term stability of PCM, the starting point is to carry outakinetic analysis of the TG measurements in order to determine the possible degradation processes [1, 2]. The majority of kinetic methods used in the area of thermal analysis consider the reaction rate equation [3]: which is a function of only two variables: temperature, T, and conversion degree, 𝛼 = m ini −m mini . The dependence of the process rate on temperature is represented by the rate constant, k(T) which is typically parameterized through the Arrhenius equation in which A is the so-called frequency or pre-exponential factor, E is the activation energy and R is the molar gas constant. The dependence on the conversion degree is represented by f(α), which is a function whose formulation depends on the mathematical model describing the reaction mechanism. For the particular case of dynamic TG measurements taken at constant heating rate 𝛽 = dT dt , Eq.1 becomes: (1) d 𝛼 dt =k(T)f(𝛼)=Ae ( -E RT ) f(𝛼 ) (2) d 𝛼 dT =Af (𝛼)e -E/RT 𝛽 This equation represents the conversion rate, usually named as dTG-T curve, which is also measured by the TG apparatus. By integrating Eq.2, the variation of the conversion degree with temperature (α-T curve) can be obtained: In this equation, g(α) is a function calculated from f(α) so that it also depends on the model used to describe the reaction mechanism [3]. ∫T T 0 e-E/RTd T is the temperature integral [4], which cannot be calculated explicitly so that either numerical methods or analytical approximations have been proposed by different authors for its evaluation [5]. In our case, we have taken the Coats and Redfern approximation for calculating the temperature integral [6]. The aim of any kinetic analysis is to obtain the kinetic triplet: E, A and f(α). If the reaction mechanism f(α) is unknown, the activation energy, E, can be estimated in a first approach by using one or more model-free isoconversional methods that are already described in the literature. These methods can be grouped in differential, if they are based in Eq.2, or integral, if they are based in Eq.3 [3]. Among the differential, we have the method of Friedman [7], and among the integral, we have the methods of Coats and Redfern [6], Ozawa [8], Flynn–Wall [9], Akahira and Sunose [10], Kissinger [11], and Miura and Maki [12]. These are what we can call the traditional isoconversional methods; however, there are other methods more sophisticated that are based on numerical integration like the one developed by Vyazovkin [13, 14]. If the reaction takes place in a single-step process, E should not vary significantly with α, and hence, A and f(α) could be obtained by applying some of the methods already stablished by the Kinetics Committee of the International Confederation for Thermal Analysis and Calorimetry (ICTAC) [3, 15]. Examples of these methods are Kissinger [11], Coats–Redfern with discriminating approach [6] or the method based on the so-called master plots [3]. Conversely, for multi-step processes, the E calculated with the isoconversional methods is expected to change with α due to the different relative contributions of each single-step to the overall reaction rate. The occurrence of multi-step processes is also detected if dTG-T curve does not present a unique maximum but it presents various peaks, which may be more or less overlapping. In such case, each reaction step should be separated and analyzed independently [15]. There are several ways to separate overlapping peaks of single-step reactions. The most common one consists in performing a mathematical deconvolution analysis (MDA) by using asymmetric peak functions [15]. However, some authors have proposed alternative methods based in (3) 𝛼 ∫ 0 d𝛼 f(𝛼)=g(𝛼)=A 𝛽 T ∫ T 0 eE/RTdT 13881 Assessment ofisoconversional methods andpeak functions forthekinetic analysis of… nonlinear regressions for estimating kinetic parameters of complex multi-step processes [16]. After single-step processes are separated, the kinetic triplet for each reaction could in principle be obtained by applying the different methods mentioned above. One of the most widely used functions for MDA of dTG-T curves is the one proposed by Fraser and Suzuki (FS) [17, 18]. In fact, there are many examples in the literature that use FS function in the kinetic analysis of thermogravimetric measurements for different kinds of materials undergoing multi-step reactions [19–22]. Some authors have used other peak functions like Gaussian, Lorentzian or Weibull for fitting dTG-T curves [19], but in the end the best fitting was always obtained with the FS curve. There is, however, a peak function that, to our knowledge, has never been used for MDA of dTG-T curves but could represent and advantage over the FS function due to its less restrictive mathematical formulation. This function is the generalized logistic probability density function (GLOG) [23] also known as Richards curve. In this context, this work had three main objectives: The first one was to assess which isoconversional method is the most adequate for having a preliminary estimation of the activation energy, E, of multi-step reaction process; the second one was to determine which peak function better represents the different reaction mechanisms, i.e., which one fits dTG-T curves more accurately and, the third one, was to find a correlation between the kinetic triplet and the parameters of the fitting functions. For that purpose, theoretical kinetic curves of conversion degree (α-T) and conversion rate (dTG-T) were constructed for different reaction mechanisms, f(α), known Arrhenius parameters, E and lnA, and various heating rates, β. Also similar sets of curves were simulated for multi-step reaction mechanisms with two parallel single-step independent reactions. Both single-step and multi-step theoretical curves were analyzed by different isoconversional methods in order to obtain the variation of activation energy with the conversion degree (E-α) and determine which method is the most adequate for analyzing real experimental data from TG measurements. Subsequently, the simulated dTG-T curves of single-step reactions were fitted with two peak functions: the normalized Fraser–Suzuki (NFS) and the generalized logistic (GLOG). With this analysis, it would be possible to determine which function is the most appropriate for representing each reaction mechanism and also the relationship between both the function parameters and the kinetic parameters used in the simulations. Additionally, mathematical deconvolution analysis (MDA) was applied to dTG-T curves simulated for multi-step reaction mechanisms so that the single-step reaction curves were separated. The kinetic parameters E and lnA of these single-step reactions were obtained and then compared with the corresponding values used in the simulations. Finally, the procedure of kinetic analysis developed for the theoretical kinetic curves was applied to the particular case of polyethylene glycol (PEG) with two molecular weights: 6000 and 12000. This PCM was chosen because up to now the kinetics of thermal degradation of polyethylene glycol have not been obtained. In fact, all references found in the literature where its thermal stability is studied establish the temperature interval in which degradation occurs and the temperature at which degradation rate is maximum by analyzing TG measurements taken at only one heating rate [24–28]. Therefore, in order to have a first insight into PEG thermal degradation kinetics, TG measurements were taken for both PEG3000 and PEG 12000 under N2 and air atmospheres at heating rates from 2 to 20K min−1. These measurements were analyzed by the different isoconversional methods applied in this work so that E-α curves were obtained. Also, mathematical deconvolution analysis was applied to dTG-T curves by using the most appropriate peak functions so that single-step reactions were separated and the corresponding activation energies calculated. Finally, the activation energy values obtained by both procedures were compared and discussed. Construction ofthetheoretical kinetic curves According to the previous literature, the mechanisms of the chemical reactions can be described by specific mathematical expressions for the function f(α) [3]. Šesták–Berggren general equation puts together most of the kinetic mechanisms by using a unique empirical equation with three parameters: m, n and p. f(α) function associated with different reaction mechanisms can be obtained by giving specific values to the Šesták–Berggren parameters. In Table1, the mathematical expressions of both f(α) and g(α) are recorded for the reaction mechanisms whose kinetics have been theoretically simulated in this work. Using Eqs.2 and 3, sets of theoretical kinetic curves α-T and dTG-T were simulated for all the mechanisms displayed in Table1 at the heating rates, β, recorded in Table2. The values of E and lnA used in the simulations are also displayed in Table2, and they were chosen by taking into account the values considered by other authors for constructing similar theoretical curves to carry out kinetic analysis [19, 29–32]. Moreover, according to our previous experience (4) f(𝛼)=𝛼m(1−𝛼)n[−ln(1−𝛼)]p 13882 R.Bayón et al. in kinetic analysis of TG measurements of various PCMs [2, 33], the values of E and lnA of Table2 selected for the theoretical simulations are in good agreement with the ones experimentally obtained. The shape of simulated α-T and dTG-T curves depended on the reaction mechanism and hence on the associated function f(α) being the differences more clearly observed in dTG-T curves (i.e., peak curves). In this sense, mechanisms F1, A2, A3 and A4 led to quite symmetric curves, mechanisms R2, R3, F2 and D3 led to curves with higher asymmetry, whereas mechanisms F0, D1 and P2 led to spikelike curves. In Fig.1, simulated dTG-T curves for the heating rates recorded in Table2 are displayed for the reaction mechanisms A2 (a), R2 (b) and D1 (c), each one representing an example of curve shape. The values of E and lnA used in the simulations are displayed in each graph. The curves simulated for the other reaction mechanisms included in Table1 are shown in Fig.S1 of the supplementary material. It is important to mention here that not all combinations of E and lnA were used for simulating all mechanisms since the effect of these parameters on dTG-T curves is only the temperature interval they cover. This is clearly shown in Fig.1 (d), where curves simulated with β = 2K min−1 and different values of E and lnA are displayed. In Fig.S1 of supplementary material, further examples of this effect can be found for various reaction mechanisms. Following a procedure similar to the one used by Granado [29] and Sbirrazzouli [32], kinetic curves (α-T and dTG-T) were simulated for multi-step reaction mechanisms with two parallel independent reactions whose overall reaction rate can be expressed as: Combinations of two reaction mechanisms with different kinetic parameters were simulated including one of the cases already studied by Sbirrazzouli [32] (CASE 7). The combination of mechanisms with their corresponding kinetic parameters is given in Table3. In Fig.2, both (1-α)-T (a) and dTG-T (b) curves calculated at 2K min−1 are displayed for each combination CASE. In Fig.2a, we can see how (1-α)-T curves of CASES 1, 2 and 6 clearly show one step at the α value that corresponds to the contribution fraction of each single-step reaction. However, (1-α)-T curves of CASES 3, 4 and 5 do not show any clear step, which means that both mechanisms are highly overlapping. This behavior is also observed in dTG-T curves of Fig.2b; while CASES 1, 2 and 6 show two peaks clearly separated and quite independent, the curves of CASES 3, 4, 5 and 7 present overlapping peaks. (5) d 𝛼 dt =x1A1e ( −E1 RT ) f1 ( 𝛼1 ) +x2A2e ( −E2 RT ) f2 ( 𝛼2 ) with 𝛼 = 𝛼 1+ 𝛼 2and x 1+ x 2=1 Table 1 Reaction mechanisms and their corresponding f(α) and g(α) expressions used for simulating the theoretical kinetic curves α-T and dTG-T [3] Model ID Reaction mechanism f( 𝛼 ) g( 𝛼 ) F0 Zero order 1 α F1 (Mampel) First order (1−𝛼)1 −ln(1−𝛼) F2 Second order (1−𝛼)2 (1−𝛼)−1−1 D1 One-dimensional diffusion 1 2𝛼 𝛼2 D3 Three-dimensional diffusion [ 3(1−𝛼) 2∕3] ∕ [ 2 ( 1−(1−𝛼) 1∕3)] [ 1−(1−𝛼)1∕3 ]2 Random nucleation and instantaneous growth of nuclei A1≡F1 Avrami–Erofeev (n = 1) (1−𝛼)1 −ln(1−𝛼) A2 Avrami–Erofeev (n = 2) 2(1−𝛼)[−ln(1−𝛼)]1∕2 [−ln(1−𝛼)]1∕2 A3 Avrami–Erofeev (n = 3) 3(1−𝛼)[−ln(1−𝛼)]2∕3 [−ln(1−𝛼)]1∕3 A4 Avrami–Erofeev (n = 4) 4(1−𝛼)[−ln(1−𝛼)]3∕4 [−ln(1−𝛼)]1∕4 R2 Phase boundary controlled (contracting area or cylinder) 2(1−𝛼)1∕2 1−(1−𝛼)1∕2 R3 Phase boundary controlled (contracting volume or sphere) 3(1−𝛼)2∕3 1−(1−𝛼)1∕3 P2 Power law 2𝛼1∕2 𝛼1∕2 Table 2 Kinetic parameters used in the simulation of α-T and dTG-T curves E/kJ mol−1 lnA/s−1 β/K min−1 60 10, 12 2, 5, 10, 20 80 15, 17, 18, 19 2, 5, 10, 20 90 18, 19 2, 5, 10, 20 100 20, 25 2, 5, 10, 20 120 25, 30 2, 5, 10, 20 13883 Assessment ofisoconversional methods andpeak functions forthekinetic analysis of… As example, Fig.3 shows the evolution of dTG-T curves with heating rate, β, for CASE 1, with slight peak overlapping (a) and CASE 3 (b) with strong peak overlapping. Model‑free isoconversional methods fortheanalysis ofsingle‑ andmulti‑step reaction mechanisms Model-free isoconversional methods allow the activation energy to be estimated as a function of conversion without choosing any reaction mechanism. The basic assumption of these methods is that the reaction rate at constant conversion degree, α, depends only on temperature so that constant E values can be expected [3]. However, if multi-step processes occur, E varies with α due to the different relative contributions of each single-step to the overall reaction rate and also to the fact that those reactions can take place either in parallel or in consecutive form. For this reason, the results from the isoconversional methods may not be fully accurate, but still can serve as a preliminary estimation and, of course, for 300 0.00 0.01 0.02 d  /dT/K–1 d  /dT/K–1 d  /dT/K–1 d  /dT/K–1 0.03 0.04 0.00 0.01 0.02 (d) (b) (a) (c) 0.03 0.04 0.06 0.05 0.00 375 400 425 450 475 500 525 550 375 400 425 450 475 500 0.01 0.02 0.03 0.04 0.05 0.00 0.01 0.02 0.03 0.04 0.05 325 D1 mechanism A2 mechanism R2 mechanism F1 mechanism E = 120 kJ mol–1 E = 90 kJ mol–1 -InA = 25s–1 E = 100 kJ mol–1 -InA = 25s–1 E = 80 kJ mol–1 -InA = 15s–1 E = 80 kJ mol–1 -InA = 18s–1  = 2 Kmin–1 InA = 30s–1 E = 80 kJ mol–1 InA = 15 s–1 E = 120 kJ mol–1 InA = 25 s–1 350 375 T/K T/K T/K T/K 400 425 450 300 350 400 450 500 550 600  = 2 K min–1  = 5 K min–1  = 10 K min–1  = 20 K min–1  = 2 K min–1  = 5 K min–1  = 10 K min–1  = 20 K min–1  = 2 K min–1  = 5 K min–1  = 10 K min–1  = 20 K min–1 Fig. 1 Examples of dTG-T curves simulated for different reaction mechanisms f(α): A2 (a), R2 (b), D1 (c) and effect of the kinetic parameters E and lnA in dTG-T curve of 2Kmin.−1 for the mechanism F1 (d) Table 3 Kinetic parameters of single-step reactions 2and 2used for simulating the multi-step reaction mechanismsof Eq.5 Mechanism E1, 2 /kJ mol−1 lnA1, 2/s−1 x1, 2 CASE 1 F1 80 18 0.5 A2 100 20 0.5 CASE 2 F1 80 17 0.8 A2 100 20 0.2 CASE 3 F1 80 17 0.5 A2 90 18 0.5 CASE 4 F2 60 15 0.5 R3 120 30 0.5 CASE 5 F2 60 15 0.2 R3 120 30 0.8 CASE 6 D3 80 18 0.5 P2 120 25 0.5 CASE 7 [32]F1 80 19 0.5 F1 90 19 0.5 13884 R.Bayón et al. checking whether a reaction is composed by one or various single-step processes [15]. The traditional isoconversional methods are based on linear plots with 1⁄T in the abscissa axis when the same conversion degree occurs. This is done for TG measurements taken at different heating velocities, and the corresponding activation energy is usually obtained from the slope of the linear plot. Table4 summarizes the mathematical expressions of the isoconversional methods applied in this work for analyzing the simulated kinetic curves [6, 7, 9, 10]. These four methods were implemented in a self-developed MATLAB© code for automatic calculation. The code allows obtaining E-α curves for each of the four methods recorded in Table4 and also the intercepts of both Friedman 1.0 0.07 0.06 0.05 0.04 0.03 d α /dT/K–1 0.02 0.01 0.00 0.9 0.8 0.7 0.6 0.5 CASE 1: F1+A2 CASE 2: F1+A2 CASE 3: F1+A2 CASE 4: F2+R3 CASE 5: F2+R3 CASE 6: D3+P2 CASE 7: F1+F1 CASE 1: F1+A2 CASE 2: F1+A2 CASE 3: F1+A2 CASE 4: F2+R3 CASE 5: F2+R3 CASE 6: D3+P2 CASE 7: F1+F1 0.4 0.3 0.2 0.1 0.0 250300 350 T/K 1α 400450 500250 300350 T/K 400450 500 β = 2 K min–1 β = 2 K min–1 (a) (b) Fig. 2 (1-α)-T (a) and dTG-T (b) curves simulated for CASES 1–7 that combine two single-step reaction mechanisms (see Table3 for kinetic parameters details) 0.03 0.02 0.01 F1 A2 F1 A2 β = 2 K min–1 β = 5 K min–1 β = 10 K min–1 β = 20 K min–1 β = 2 K min–1 β = 5 K min–1 β = 10 K min–1 β = 20 K min–1 0.00 300 CASE 1: F1+A2 CASE 3: F1+A2 325350 375400 425 T/K d α /dT/K–1 d α /dT/K–1 450475 500525 550300 325350 375400 425 T/K 450475 500525 550 0.03 0.02 0.01 0.00 (a) (b) Fig. 3 Evolution of dTG-T curves with heating rate, β, for simulated for CASE 1 with slight peak overlapping (a) and CASE 3 (b) with strong peak overlapping (see Table3 for kinetic parameters details) Table 4 Isoconversional methods used in this work for the kinetic analysis of simulated kinetic curves [3] Method Kind Linear plot equation Friedman (FR) Differential ln[( d𝛼 dt )𝛼] =ln [ 𝛽 ( d𝛼 dT )𝛼] =ln [ f(𝛼)A𝛼 ] −E RT 𝛼 Coats–Redfern (CR) Integral ln 𝛽 T2=ln [ AR Eg ( 𝛼 )( 1−2RT E )] −E RT Kissinger–Akahira–Sunose (KAS) Integral ln 𝛽 T 2 =Const. − E RT Flynn–Wall–Ozawa (FWO) Integral dlog𝛽 d1∕T ≅ 0.457 R E 13885 Assessment ofisoconversional methods andpeak functions forthekinetic analysis of… and Coats–Redfern ones. It is important to mention here that in this work we have only used the traditional isoconversional methods because according to the results obtained by Luciano and Svoboda [29], some of these methods lead to similar E-a curves as the numerical method developed by Vyazovkin [13, 14]. For assessing the most accurate isoconversional methods of Table4, kinetic curves (α-T and dTG-T) constructed for all the reaction mechanisms included in Table1 with some of the parameters of Table2 were analyzed. In Fig.4, the variation of E with α obtained from the different isoconversional methods is displayed as example for various reaction mechanisms: F1 (a), F1 (b), F2 (c), F2 (d), D1 (e), A2 (f), R3 (g) and P2 (h). The corresponding E and lnA values used in the simulation of the kinetic curves are included in each graph. The E-α curves for the mechanisms F0, D3, A3, A4 and R2 mechanisms are given in Fig.S2 of supplementary material. As we can see, except Flynn–Wall–Ozawa method, which clearly leads to the less accurate values, the other isoconversional methods lead to E values very close to the used in the simulation of the kinetic curves. In Fig.5, the mean relative error of E has been represented for all the reaction mechanism simulated with different activation energies. Clearly the method leading to the highest error is the Flynn–Wall–Ozawa; however, it is interesting to note that this error increases as activation energy decreases. On the other hand it is also clear that the most accurate results are obtained with the Friedman’s, which uses dα/dT values (i.e., dTG-T curves) in the linear plots for estimating the activation energy (see Table4). For this method, it is also possible to calculate lnA from the intercept of the linear fit performed for each α value. In Fig.6, lnA calculated from the intercept of Friedman method has been plotted for various reaction mechanisms with different Arrhenius parameters (E and lnA). As we can see, lnA values calculated from the intercept are coincident with the corresponding ones used in the simulations of the kinetic curves. From these results, we can conclude that, although all isoconversional methods lead to activation energy values within an error range quite low (< 7%), Friedman method is clearly the most accurate, at least for the case of single-step reaction mechanisms. Actually, being a differential isoconversional method that does not make use of any approximation, it should potentially be more accurate than the integral methods. However, as discussed by Vyazovkin etal. [3] the differential methods should not be considered as being necessary more accurate and precise than the integral ones. In their work, Luciano and Svoboda [29] also mention the relative errors of both KAS and FWO methods. For the first one, they obtained 0,15% while for the second one the error was about 0,8% and hence higher. In our case, relative errors are clearly above the values obtained by Luciano and Svoboda, but this may be due to the fact that our curves were simulated at a much lower points/curve density (100 instead of 10,000) that the E values used by these authors are in the range of 110–300kJ mol−1 but also that they only simulated F1 and A2 mechanisms. In fact, our calculations for mechanism F1 simulated with E = 120kJ mol−1 showed that relative errors of KAS and FWO methods were 0,228 and 1,06%, respectively, so that very close to the values obtained by Luciano and Svoboda. Isoconversional analysis was also applied to the multistep reaction mechanisms of Table3 so that E-α curves were obtained for each method used in this work. The results for all simulated cases (1–7) are displayed in the graphs of Fig.7a–g with the E values used in the single-step mechanisms included for comparison. As we can see, for CASES 1 and 2 in which the two mechanisms are clearly differentiated (i.e., peaks in dTG-T curves are slightly overlapping), the transition between the two E values is a well-defined step. This step occurs at the same α value where the step was observed in (1-α)-T curves (see Fig.2 a), which also corresponds to the contribution of each reaction mechanism to the overall one: CASE 1 with x1 = x2 = 0,5 and CASE 2 with x1 = 0.8 and x2 = 0.2. As occurred in the analysis of the single-step mechanisms, the curves from FWO method deviate from the E values used in the simulations more than the curves obtained with the other methods, especially when the energy values of the singlestep mechanisms are close (CASES 1, 2, 3 and 7). However, it leads to most accurate activation energy values when those energy values are not very close (CASES 4, 5 and 6). As for FR method, it is interesting to note that it usually produces strong discontinuities (CASES 1, 2 and 6) in the transition interval of E and, in many cases, strong deviations from the activation energy values used in the simulations (CASES 3, 4 and 7). In contrast, KAS and CR methods clearly lead to the most accurate values of activation energies for all simulated cases of multi-step reaction mechanisms, with a smooth and continuous variation of E in the transition interval. It is important to mention that the energy value used for the mechanism R3 (120kJ mol−1) is only attained in CASE 5, in which the contribution of F2 mechanism is x1 = 0.2 but it is not attained in CASE 4 even for the highest values of conversion when the contribution of F2 mechanism is x1 = 0.5. This must be due to the fact that both mechanisms are highly overlapping, and hence, the one with the lowest energy prevents the one with the highest energy to prevail. Similar results were obtained by Luciano and Svoboda [29] when they compared the results of FR, KAS and FWO methods applied to kinetic curves simulated for multi-step reaction mechanism with different extents of peak overlapping. From these results, we can conclude that for the case of multi-step reactions, the traditional isoconversional methods are able to give only an estimation of the range of the activation energies associated with the different single-step 13886 R.Bayón et al. 110 108 106 104 102 100 98 96 94 92 90 0.0 132 130 128 126 124 122 120 118 116 114 112 110 108 0.2 E = 100 kJ mol–1 F1 mechanism F2 mechanism F2 mechanism E/kJ mol–1 InA = 25 s–1 0.4 α 0.60.8 1.0 FR KAS CR FWO 0.00.2 E = 120 kJ mol–1 E/kJ mol–1 InA = 25 s–1 0.4 α 0.60.8 1.0 FR KAS CR FWO 132 130 128 126 124 122 120 118 116 114 112 110 108 D1 mechanism A2 mechanism 0.00.2 E = 120 kJ mol–1 E/kJ mol–1 InA = 30 s–1 0.4 α 0.60.8 1.0 FR KAS CR FWO 0.0 88 86 84 82 80 78 76 74 72 0.2 E = 80 kJ mol–1 E/kJ mol–1 InA = 15 s–1 0.4 α 0.60.8 1.0 FR KAS CR FWO R3 mechanism 0.0 88 86 84 82 80 78 76 74 72 0.2 E = 80 kJ mol–1 E/kJ mol–1 InA = 18 s–1 0.4 α 0.60.8 1.0 FR KAS CR FWO P2 mechanism 0.0 66 65 64 63 62 61 60 59 58 57 56 55 54 0.2 E = 60 kJ mol–1 E/kJ mol–1 InA = 12 s–1 0.4 α 0.60.8 1.0 FR KAS CR FWO 0.0 66 65 64 63 62 61 60 59 58 57 56 55 54 0.2 E = 60 kJ mol–1 E/kJ mol–1 InA= 12 s–1 0.4 α 0.60.8 1.0 FR KAS CR FWO 66 65 64 63 62 61 60 59 58 57 56 55 54 0.00.2 E = 60 kJ mol–1 F1 mechanism E/kJ mol–1 InA = 10 s–1 0.4 α 0.60.8 1.0 FR KAS CR FWO (a) (b) (c) (d) (e) (f) (g) (h) 13887 Assessment ofisoconversional methods andpeak functions forthekinetic analysis of… mechanisms. Moreover, although Friedman differential method resulted to be the most accurate when analyzing single-step reaction cases, we have seen that the integral methods of Coats–Redfern and Kissinger–Akahira–Sunose lead to better estimations when multi-step reactions are taking place. Mathematical deconvolution analysis ofdTG‑T curves Definition ofthepeak functions As discussed above, the use of isoconversional methods in the kinetic analysis of TG measurements is expected to lead to reliable results when single-step reactions occur. However, the most usual situation in real cases is that different processes of decomposition overlap each other so that E varies with α. In general, the overlapping processes are better observed in dTG-T curves (see Fig.2b) and one of the procedures to separate themis the so-calledmathematical deconvolution analysis (MDA). In this kind of analysis, two or more peak curves are used for fitting the experimental dTG-T curves so that single-step reaction mechanisms can be separated and subsequently analyzed. According to the previous experience of different authors [19–22], Fraser–Suzuki (FS) function is very appropriate for deconvoluting curves into single peaks. Most of them have used FS function for deconvoluting either dα/dt-t curves [19, 22] or dα/dT-T curves [20], and in some cases also for deconvoluting heat flow curves [21]. In this work, FS function has been be applied to dTG-T curves so that the general mathematical expression is: The four parameters of the function are: • h: amplitude (peak height) • s: shape parameter • p: position • w: half height width The main advantage of using dα/dT-T curves in the deconvolution is that their integral along the temperature range should be equal to 1 because, usually, the final (6) d 𝛼 d T=hexp { −ln2 s2 [ ln(1+2sT−p w) ] 2 } conversion degree of a reaction is α = 1, and this applies for both singleand multi-step mechanisms. Therefore, if a dα/dT-T curve is deconvoluted in various curves with a certain contribution, this implies that the integral of each single curve must be equal to 1. Hence, FS curves used in the deconvolution process should meet this requirement, i.e., have a normalized area. According to the authors [17], the area under the Fraser–Suzuki curve is calculated with the following expression, which depends on h, w and s: If this area has to be equal to 1, then the parameter h can be expressed in terms of the other two: Therefore, the expression for the normalized Fraser–Suzuki (NFS) function has only 3 parameters: s, p and w. It is important to note that, in the field of real numbers, the argument of a natural logarithm cannot be zero or negative, so this has to be taken into account in the computing procedure used in dα/dT-T curve deconvolution. For overcoming this issue, in this work we wanted to check another asymmetric peak function for fitting deconvoluting dTG-T curves. The selected alternative curve is the generalized logistic probability density function (GLOG) [23] whose mathematical expression is: Being a statistic function, it is already normalized so that the area under the curve is equal to 1. GLOG function has three parameters as well: peak position, Tc; B and Q, while peak height is calculated from the first derivative at Tc as: It is important to highlight that both NFS a GLOG functions have three parameters, and hence, many combinations could lead to good fitting in a certain MDA procedure. This means that stablishing the range for the initialization parameters and their most appropriate (7) area =hw 2 exp ( s 2 4ln2)( 𝜋 ln 2)1∕2 (8) h =2 wexp ( −s2 4ln2 )( ln2 𝜋 )1∕2 (9) d 𝛼 d T= [ 2 wexp ( −s2 4ln2 )( ln2 𝜋 ) 1∕2 ] exp { −ln2 s2 [ ln(1+2sT−p w) ] 2 } (10) d 𝛼 d T= B eB(T-Tc) [ 1+Qe-B(T-Tc) ] (1 Q+1 ) (11) d 𝛼 d T |||| Tc =B (1+Q) (1 Q+1 ) Fig. 4 Results of FR, KSAS, CR and FWO isoconversional methods for simulated kinetic curves of different reaction mechanisms and Arrhenius parameters. F1 (a), F1 (b), F2 (c), F2 (d), D1 (e), A2 (f), R3 (g) and P2 (h) ◂ 13894 R.Bayón et al. 600 70 60 50 40 30 20 550 500 E80-InA15 F1 mechanism (a) (b) F1 mechanism E80-InA18 E80-InA20 E90-InA20 E90-InA25 E100-InA25 E80-InA15 E80-InA18 E80-InA20 E90-InA20 E90-InA25 E100-InA25 450 p/K w/K 400 350 300 0510 15 β /K min –1 20 25 0510 15 β /K min –1 20 25 Fig. 9 Variation of NFS parameters (p and w) with β, E and lnA for a reaction mechanism F1 600 7 6 5 4 3 2 550 500 F0-E100-InA25 F0-E100-InA20 F0-E120-InA25 (a) (b) D1-E120-InA25 D1-E120-InA30 D1-E100-InA25 F0-E100-InA25 F0-E100-InA20 F0-E120-InA25 D1-E120-InA25 D1-E120-InA30 D1-E100-InA25 450 400 350 0510 β /K min–1 Tc/K B/K–1 15 20 25 0510 β /K min–1 15 20 25 Fig. 10 Variation of GLOG parameters: Tc (left) and B (right) with β and Arrhenius parameters, for reaction mechanisms F0 and D1 13895 Assessment ofisoconversional methods andpeak functions forthekinetic analysis of… CASE 1: 0.5F1+0.5A2 CASE 3: 0.5F1+0.5A2 CASE 4: 0.5F2+0.5R3 CASE 6: 0.5D3+0.5P2 0.025 0.030 0.020 Simulated curve GLOG Dec 0.5F1 GLOG Dec 0.5A2 β = 2 K min–1 β = 2 K min–1 β = 2 K min–1 β = 2 K min–1 GLOG Dec Total Simulated curve GLOG Dec 0.5D3 GLOG Dec 0.5P2 GLOG Dec Total Simulated curve NFS Dec 0.5D3 NFS Dec 0.5P2 NFS Dec Total Simulated curve NFS Dec 0.5F2 NFS Dec 0.5R3 NFS Dec Total 0.015 0.010 0.03 0.02 0.01 0.00 0.005 0.000 0.025 0.030 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0.00 0.020 0.015 0.010 0.005 0.000 300 275 300 325 350 375 400 425 325 350 375 400 T/K T/K 275 300 325 350 375 400 425 450 475 T/K d α /dT/K–1 d α /dT/K–1 d α /dT/K–1 d α /dT/K–1 425 450 475 500 350 375 400 425 450 475 T/K (a) (b) (c) (d) Fig. 11 Deconvolution of dTG-T curves for some simulated multi-step reaction mechanisms: a CASE 1; b CASE 3; c CASE 4; and d CASE 6 (see Table3) Table 6 Kinetic parameters obtained with Kissinger method for the multi-step reaction mechanisms simulated and then deconvoluted with both NFS and GLOG functions E and lnA from Kissinger method Simulated NFS GLOG Mechanism x E/kJ mol−1 lnA/s−1 E/kJ mol−1 lnA/s−1 E/kJ mol−1 lnA/s−1 CASE 1 F1 0.5 80 18 79.71 (0.36%) 17.89 (10.4%) 79.90 (0.13%) 17.97 (2.9%) A2 0.5 100 20 99.73 (0.27%) 20.02 (22.14%) 99.63 (0.37%) 19.94 (5.8%) CASE 2 F1 0.8 80 17 78.50 (1.87%) 79.76 (0.30%) A2 0.2 100 20 100.39 (0.39%) 99.45 (0.55%) CASE 3 F1 0.5 80 17 79.58 (0.53)% 79.07 (1.16%) A2 0.5 90 18 89.54 (0.51%) 89.89 (0.12%) CASE 4 F2 0.5 60 15 59.57 (0.72%) 14.81 (17.3%) 57.52 (4.13%) R3 0.5 120 30 120.30 (0.25%) 32.76 (1479%) 119.96 (0.03%) CASE 6 D3 0.5 80 18 78.47 17.43 (43.4%) 79.11 17.63 (30.9%) P2 0.5 120 25 119.46 24.81 (17.3%) 119.24 24.74 (22.9%) CASE 7 F1 0.5 80 19 79.97 (0.04%) 79.04 (1.2%) F1 0.5 90 19 89.93 (0.08%) 89.30 (0.08%) 13896 R.Bayón et al. 0.020 0.030 0.025 0.020 0.015 0.010 0.005 0.000 PEG 12000 under N2 (a) (b) PEG 12000 under air 0.015 0.010 0.005 0.000 400 500 600 β = 2 K min–1 β = 5 K min–1 β = 10 K min–1 β = 20 K min–1 β = 2 K min–1 β = 5 K min–1 β = 10 K min–1 β = 20 K min–1 T/K d α /dT/K–1 d α /dT/K–1 700 800 400 500 600 T/K 700 800 Fig. 12 dTG-T curves of PEG 12000 obtained at different heating rates under air (a) and N2 (b) atmospheres 200 125 100 75 50 25 0 FR KAS-CR FWO FR KAS-CR FWO 150 100 E/kJ mol–1 E/kJ mol–1 50 0 0.00.2 PEG 3000 in N2PEG 3000 in air 0.40.6 α 0.81.0 0.00.2 0.40.6 α 0.81.0 200 125 100 75 50 25 0 FR KAS-CR FWO FR KAS-CR FWO 150 100 E/kJ mol–1 E/kJ mol–1 50 0 0.00.2 PEG 12000 in N2PEG 12000 in air 0.40.6 α 0.81.0 0.00.2 0.40.6 α 0.81.0 (a) (b) (c) (d) Fig. 13 E vs. α curves of PEG3000 and PEG12000 for TG measurements taken under N2 and air obtained from the isoconversional methods used in this work 13897 Assessment ofisoconversional methods andpeak functions forthekinetic analysis of… combined kinetic analysis and master plots, which for the moment are beyond the scope of this paper. Conclusions In this work, theoretical kinetic curves α-T and dTG-T were simulated for both single-step and multi-step reaction mechanisms by using different sets of kinetic triplets (f(α), E and lnA). These curves were analyzed with various isoconversional methods so that E-α curves were obtained and compared with the energy values used in the simulations. For the single-step reaction mechanisms, all methods lead to activation energies within an error range below 7% in relation to the values used in the simulations, being Friedman method the most accurate. For the multi-step reaction mechanisms, the obtained E-α curves were only a rough estimation of the activation energies used in the simulations and, in this case, Kissinger–Akahira–Sunose and Coats–Redfern methods were the ones leading to the most accurate results. Simulated dTG-T curves of single-step mechanisms were fitted with two kinds of peak functions (normalized Fraser–Suzuki—NFS and generalized logistic—GLOG) in order to determine which one better represented the different reaction mechanisms. Although both functions led to regression values above 0,99, the NFS proved to better fit dTG-T curves of the majority of reaction mechanisms analyzed. From the fitting, the relationship between the kinetic triplet used for simulating dTG-T curves and the parameters of NFS and GLOG functions was obtained as well. Both functions had only one parameters that depended on the reaction mechanism used in the simulation (s for NFS; Q for GLOG) while the other two depended on both the Arrhenius parameters, E and lnA, and the heating rate, β (p and w for NFS; Tc and B for GLOG). These relationships 0.030 0.015 0.010 0.005 0.000 0.020 0.015 0.010 0.005 0.000 0.025 0.020 PEG 12000 under air β = 2 K min–1 PEG 12000 under air β = 20 K min–1 PEG 12000 under N2 β = 20 K min–1 PEG 12000 under N2 β = 5 K min–1 0.015 Experimental data NFS1 NFS2 NFS3 NFS total Experimental data NFS1 NFS2 NFS3 NFS total Experimental data NFS1 NFS2 NFS3 NFS total Experimental data NFS1 NFS2 NFS total 0.010 0.005 0.000 400 500 600 T/K d α /dT/K–1 d α /dT/K–1 d α /dT/K–1 d α /dT/K–1 700 400 500 600 T/K 700 800 400 500 600 T/K 700 800 400 500 600 T/K 700 800 0.020 0.015 0.010 0.005 0.000 (a) (b) (c) (d) Fig. 14 Examples of dTG curves deconvolutions using NFS functions for PEG12000 at different heating rates and TG measurements under air and N2 atmospheres Table 7 Estimated E for deconvoluted peaks obtained with NFS functions E from Kissinger method/kJ mol−1 Atmosphere Curve 1 Curve 2 Curve 3 PEG 3000 N264 196 – Air 50 40 93 PEG 12000 N247 162 – Air 57 37 93 13898 R.Bayón et al. will be strongly helpful when experimental dTG-T curvesof multi-step reactions have to be analyzed by applying mathematical deconvolution analysis (MDA). In this sense, the values of s and Q parameters associated with the different reaction mechanisms should be used as boundary conditions. Moreover, to the author knowledge this is the first time both normalized Fraser–Suzuki and generalized logistic functions have been used for fitting dTG-T curves of different reaction mechanisms. MDA was applied to dTG-T curves of multi-step reaction mechanisms by using both NFS and GLOG functions and taking into account the relationships between the kinetic triplet and the function parameters obtained in this work. With this procedure, single-step reaction peaks were separated and their corresponding kinetic parameters calculated. In general, the curves resulting from deconvolution fitted quite well the simulated curves and the analysis of the single-step peaks with Kissinger method leads kinetic triplets quite close to the ones used in the simulations mainly for single-step mechanism with dTG-T curves with high symmetry. Compared with similar studies found in the literature, which are mostly focused in reaction mechanism F1, our results extend the validity of using either NFS or GLOG functions for deconvoluting dTG curves composed by many other reaction mechanisms. However, it must be taken into account that deconvolution of experimental dTG-T curves may be complicated if single-step reaction mechanisms are strongly overlapped and their contribution to the overall reaction varies with the conversion. A similar procedure of kinetic analysis was applied to thermogravimetric measurements taken under both air and N2 atmospheres for two PCMs: PEG3000 and PEG12000. In first stage E-α curves were obtained for the different isoconversional methods used in this work. The curves from the measurements under N2 showed E values going from around 40kJ mol−1 at low conversions up to 150kJ mol−1 at high conversions, whereas the curves from TG measurements under air showed more constant E values of about 50–75kJ mol−1. The lower activation energies obtained for the measurements under air aremost probably an indication of the polymer combustion. Finally, experimental dTG-T curves were deconvoluted with the most appropriate peak functions in order to have a preliminary estimation of the possible single-step reaction mechanisms occurring in these PCMs. For both PEG, fairly good deconvolutions were achieved with either 3 or 2 NFS functions and the resulting single-step reaction curves were analyzed with the Kissinger method. The activation energies obtained from Kissinger method were in good agreement with the E-α curves calculated with the isoconversional methods. As for the other kinetic parameters associated with the single-step mechanisms, they would require a further kinetic analysis, which is beyond the scope of this paper. Supplementary Information The online version contains supplementary material available at https:// doi. or g/ 10. 1007/ s1097302413494-w. Acknowledgements This work has been supported by Comunidad de Madrid and European Structural Funds through ACES2030 Project (S2018/EM-4319), the European Union’s Horizon H2020 Research and Innovation Programme through StoRIES Project (GA Nº 101036910) and SFERA III Project (GA Nº 823802) and also by the STES4Dmat project, Grant TED2021-131061B-C33 funded by MCIN/AEI/ https:// doi. org/ 10. 13039/ 50110 00110 33 and by the “European Union NextGenerationEU/PRTR.” Author contributions Rocío Bayón and Redlich García-Rojas were involved in conceptualization, methodology, formal analysis and investigation; Rocío Bayón took part in writing-original draft preparation; Rocío Bayón, Redlich García-Rojas, Esther Rojas and MargaritaM. Rodríguez-García participated in writing-reviewing and editing; and Rocío Bayón, Esther Rojas and Margarita M.Rodríguez-García were responsible for funding acquisition and resources. Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. References 1. Bayón R, Rojas E. Development of a new methodology for validating thermal storage media: application to phase change materials. Int J Energy Res. 2019. https:// doi. org/ 10. 1002/ er. 4589. 2. Bayón R, Bonanos A, Rojas E. Assessing the long-term stability of fatty acids for latent heat storage by studying their thermal degradation kinetics. Proceedings Eurosun. 2020; https:// doi. org/ 10. 18086/ euros un. 2020. 07. 10 3. Vyazovkin S, Burnham AK, Criado JM, Pérez-Maqueda LA, Popescu C, Sbirrazzouli N. ICTAC kinetics committee recommendations for performing kinetic computations on thermal analysis data. Thermochim Acta. 2011. https:// doi. org/ 10. 1016/j. tca. 2011. 03. 034. 4. Doyle CD. Kinetic analysis of thermogravimetric data. J Appl Polym Sci. 1961. https:// doi. org/ 10. 1002/ app. 1961. 07005 1506. 5. Órfão J. Review and evaluation of the approximations to the temperature integral. AICHE J. 2007. https:// doi. org/ 10. 1002/ aic. 11296. 6. Coats W, Redfern JP. Kinetic parameters from thermogravimetric data. Nature. 1964. https:// doi. org/ 10. 1038/ 20106 8a0. 13899 Assessment ofisoconversional methods andpeak functions forthekinetic analysis of… 7. Friedman HI. Kinetics of thermal degradation of char-forming plastics from thermogravimetry. Application to a phenolic plastic. J Polym Sci Part C. 1964;6:183–95. https:// doi. org/ 10. 1002/ polc. 50700 60121. 8. Ozawa T. A new method for analyzing thermogravimetric data. Bull Chem Soc Japan. 1965. https:// doi. org/ 10. 1246/ bcsj. 38. 1881. 9. Flynn JH, Wall LA. A quick, direct method for the determination of activation energy from thermogravimetric data. Polymer letters. 1966. https:// doi. org/ 10. 1002/ pol. 1966. 11004 0504. 10. Akahira T, Sunose T. Method of determining activation deterioration constant of electrical insulating materials. Res Report Chiba Inst Technol (Sci Technol). 1971;16:22–31. 11. Kissinger HE. Reaction kinetics in differential thermal analysis. Anal Chem. 1957. https:// doi. org/ 10. 1021/ ac601 31a045. 12. Miura K, Maki T. A simple method for estimating f(E) and k0(E) in the distributed activation energy model. Energy Fuels. 1998. https:// doi. org/ 10. 1021/ ef970 212q. 13. Vyazovkin S. Evaluation of the activation energy of thermally stimulated solid state reactions under an arbitrary variation of the temperature. J Comput Chem. 1997. https:// doi. org/ 10. 1002/ (SICI) 1096987X(199702) 18:3% 3c393:: AIDJCC9% 3e3.0. CO;2-P. 14. Vyazovkin S. Modification of the integral isoconversional method to account for variation in the activation energy. J Comput Chem. 2001;22(2):178–83. 15. Vyazovkin S, Burnham AK, Favergeon L, Koga N, Moukhina E, Pérez-Maqueda LA, Sbirrazzouli N. ICTAC kinetics committee recommendations for analysis of multi-step kinetics. Thermochim Acta. 2020. https:// doi. org/ 10. 1016/j. tca. 2020. 178597. 16. Pomerantsev AL. Kinetic analysis of non-isothermal solid-state reactions: multi-stage modelling without assumptions in the reaction mechanism. Phys Chem Chem Phys. 2017. https:// doi. org/ 10. 1039/ C6CP0 7529K. 17. Fraser RDB, Suzuki E. Resolution of overlapping bands: functions for simulating band shapes. Anal Chem. 1969. https:// doi. org/ 10. 1021/ ac602 70a007. 18. Rusch PF, Lelieur JP. Analytical moments of skewed gaussian distribution functions. Anal Chem. 1973. https:// doi. org/ 10. 1021/ ac603 30a060. 19. Perejón A, Sánchez-Jiménez PE, Criado JM, Pérez-Maqueda LA. Kinetic Analysis of complex solid-state reactions. A new deconvolution procedure. J Phys Chem. 2011;115(8):1780–91. https:// doi. org/ 10. 1021/ jp110 895z. 20. Cheng Z, Wu W, Ji P, Zhou X, Liu R, Cai J. Applicability of Fraser-Suzuki function in kinetic analysis of DAEM processes and lingnocellulosic biomass pyrolysis process. J Therm Anal Calorim. 2015. https:// doi. org/ 10. 1007/ s109730144215-3. 21. Svoboda R, Málek L. Applicatility of Fraser-Suzuki function in kinetic analysis of complex crystallization processes. J Therm Anal Calorim. 2013. https:// doi. org/ 10. 1007/ s109730122445-9. 22. Stankovic B, Jovanovic J, Adnadjevic B. Application of the Suzuki-Fraser function in modelling the non-isothermal dihydroxylation kinetics of fullerol. React Kinet Mech Cat. 2018. https:// doi. org/ 10. 1007/ s111440181380-6. 23. Richards FJ. A Flexible growth function for empirical use. J Experimental Botany. 1959; https:// www. jstor. org/ stable/ 23686 557 24. Sheng M, Sheng Y, Wu H, Liu Z, Li Y, Xiao Y, Lu X, Qu J. Bio-based poly (lactic acid) shaped wood-plastic phase change composites for thermal energy storage featuring favorable reprocessability and mechanical properties. Sol Energy Mater and Sol Cells. 2023. https:// doi. org/ 10. 1016/j. solmat. 2023. 112186. 25. Chen X, Guo X, Lin X, etal. pH-responsive wood-based phase change material for thermal energy storage building material application. J Mater Sci. 2022. https:// doi. org/ 10. 1007/ s1085302207474-4. 26. Yan D, Zhao S, Ge C, Gao J, Gu C, Fan Y. PBT/adipic acid modified PEG solid-solid phase change composites. J Energy Storage. 2022. https:// doi. org/ 10. 1016/j. est. 2022. 104753. 27. Wang Z, Zhang X, Jia S, etal. Influences of dynamic impregnating on morphologies and thermal properties of polyethylene glycol-based composite as shape-stabilized PCMs. J Therm Anal Calorim. 2017. https:// doi. org/ 10. 1007/ s109730165958-9. 28. Liu Z, Zhang Y, Hu K, Xiao Y, Wang J, Zhou C, Lei J. Preparation and properties of polyethylene glycol based semi-interpenetrating polymer network as novel form-stable phase change materials for thermal energy storage. Energy Build. 2016. https:// doi. org/ 10. 1016/j. enbui ld. 2016. 06. 009. 29. Luciano G, Svoboda R. Activation energy determination in case of independent complex 2 kinetic processes. Processes. 2019;7(10):738. https:// doi. org/ 10. 3390/ pr710 0738. 30. Muravyev NV, Pivkina AN, Koga N. Critical appraisal of kinetic calculation methods applied to overlapping multistep reactions. Molecules. 2019. https:// doi. org/ 10. 3390/ molec ules2 41222 98. 31. Granado L, Sbirrazzouli N. Isoconversional computations for nonisothermal kinetic predictions. Thermochim Acta. 2021. https:// doi. org/ 10. 1016/j. tca. 2020. 178859. 32. Sbirrazzouli N. Model-free isothermal and nonisothermal predictios using advanced isoconversional methods. Thermochim Acta. 2021. https:// doi. org/ 10. 1016/j. tca. 2020. 178855. 33. Bayón R, García RJ, Quant L, Rojas E. Study of thermal degradation of adipic acid as PCM under stress conditions: a kinetic analysis. Eurosun Proceedings. 2022; https:// doi. org/ 10. 18086/ euros un. 2022. 13. 02 34. Pérez-Maqueda LA, Criado JM, Sánchez-Jiménez PE. Combined kinetic analysis of solid-state reactions: a powerful tool for the simultaneous determination of kinetic parameters and the kinetic model without previous assumptions on the reaction mechanism. J Phys Chem A. 2006. https:// doi. org/ 10. 1021/ jp064 792g. 35. Málek J. The kinetic analysis of non-isothermal data. Thermochim Acta. 1992. https:// doi. org/ 10. 1016/ 00406031(92) 85118-F. Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.