Advanced Isoconversional Kinetic Analysis for the Elucidation of Complex Reaction Mechanisms: A New Method for the Identification of Rate-Limiting Steps

Two complex cure mechanisms were simulated. Isoconversional kinetic analysis was applied to the resulting data. The study highlighted correlations between the reaction rate, activation energy dependency, rate constants for the chemically controlled part of the reaction and the diffusion-controlled part, activation energy and pre-exponential factors of the individual steps and change in rate-limiting steps. It was shown how some parameters computed using Friedman’s method can help to identify change in the rate-limiting steps of the overall polymerization mechanism as measured by thermoanalytical techniques. It was concluded that the assumption of the validity of a single-step equation when restricted to a given α value holds for complex reactions. The method is not limited to chemical reactions, but can be applied to any complex chemical or physical transformation.


Introduction
The control of polymerization reactions of thermosetting material is of major importance to reach the optimal properties of the cured polymer or composite. These complex reactions often involve several chemical and diffusion steps that make the elucidation of the reaction mechanism very difficult. The use of classical analytical techniques is often limited by high change in the viscosity of the reaction medium during polymerization. Thus, Differential Scanning Calorimetry (DSC) and rheometry have been widely used to monitor complex cure kinetics that depend on both time, temperature and heating (cooling) rate. The final degree of crosslinking, which is correlated to the extent of conversion (α), is linked to the temperature program used for curing. In addition, there is a direct link between the final properties of the polymer and (i) the molecular weight, (ii) the extent of conversion as evaluated by DSC and (iii) the glass transition temperature (T g ). As an illustration of this, the Di Benedetto equation shows the link between the glass transition temperature and the extent of conversion [1][2][3].
In order to gain more insight into the elucidation of the reaction mechanism of complex reactions, the objective of the work is to show how isoconversional kinetic analysis can provide information on the rate-limiting steps involved in complex chemical reactions or physical transformations by an analysis of the activation energy dependency, called the E α dependency. In the case of single-step processes, the activation energy computed with an isoconversional method leads to a constant value with the extent of conversion (α). Nevertheless, this case is not frequently observed. In the case of complex reactions or complex physical transformations, analysis of the E α dependency and its variations indicate the presence of a complex mechanism and may give important insight into change in the rate-limiting steps. For this purpose, a complex chemical polymerization was selected and simulated using data obtained

Theoretical Part
Isoconversional methods are amongst the more reliable kinetic methods for the treatment of thermoanalytical data, see for example [5][6][7]. The Kinetics Committee of the International Confederation for Thermal Analysis and Calorimetry (ICTAC) has recommended the use of multiple temperature programs for the evaluation of reliable kinetic parameters [6]. The main advantages of isoconversional methods are that they afford an evaluation of the effective activation energy, E α , without assuming any particular form of the reaction model, f (α) or g(α), and that a change in the E α variation, called the E α dependency, can generally be associated with a change in the reaction mechanism or in the rate-limiting step of the overall reaction rate, as measured with thermoanalytical techniques.
Polymerizations are frequently accompanied by a significant amount of heat released; thus, cure kinetics can be easily monitored by DSC. It is generally assumed that the heat flow measured by calorimetry is proportional to the process rate [5][6][7]. Thus, the extent of conversion at time t, α t , is computed according to Equation (1), as follows: where t i represents the first integration bound of the DSC signal and t f is the last integration bound selected when the reaction is finished. (dQ/dt) is the heat flow measured by DSC at time t, Q tot is the total heat released (or absorbed) by the reaction and Q t is the current heat change.
The general form of the basic rate equation is usually written as [5]: where T is the temperature, f (α) is the differential form of the mathematical function that describes the reaction model that represents the reaction mechanism, E is the activation energy and A is the pre-exponential factor. The advanced non-linear isoconversional method (NLN) [8][9][10][11] used in this study is presented in Equations (3) and (4) and was derived from Equation (2): where E α is the effective activation energy. The E α value is determined as the value that minimizes the function Φ(E α ). This non-linear kinetic method (referred as NLN) allows one to handle a set of n experiments carried out under different arbitrary temperature programs T i (t) and uses a numerical integration of the integral with respect to the time. For each i-th temperature program, the time t α,i and temperature T α,i related to selected values of α are determined by an accurate interpolation using a Lagrangian algorithm [11,12]. Numerical integration is performed using trapezoidal rule. Several possibilities are proposed for the initial estimate E 0 of E α in the non-linear procedure. The method developed by Sbirrazzuoli and implemented in his internally generated software can treat any kind of isothermal or non-isothermal data from DSC, calorimetry (C80, for example), Thermogravimetric Analysis (TGA), Dynamic Mechanical Analysis (DMA), or rheometry [9,[11][12][13][14][15][16]. This software was used in this study to compute a value of E α for each value of α between 0.02 and 0.98 with a step of 0.02. This advanced non-linear isoconversional method (NLN) was applied in this study. Another isoconversional method can be derived by the linearization of Equation (2) and is known as Friedman's method [12,17]: Application of this method requires the knowledge of the reaction rate (dα/dt) α,i and of the temperature T α,i corresponding to a given extent of conversion, for the i temperature programs used. The advantages of differential methods such as Friedman's method (referred as FR) are that they use no approximations and can be applied to any temperature program. As for NLN, the interpolation is made using a Lagrangian algorithm. This does not hold for usual integral methods, but is also the case for the non-linear advanced isoconversional method previously described. Nevertheless, simulations have shown that differential isoconversional methods can sometimes reveal numerical instability [12]; therefore, before using Friedman's method it was checked that the obtained results were consistent with those obtained with the NLN method.
Equation (5) shows that the intercept of the Friedman's plot led to the determination of the term [A α f (α)]. This term represents the product between the pre-exponential factor A α and the mathematical function f (α) that describes the reaction mechanism. Once E α and [A α f (α)] have been evaluated it is possible to compute the reaction rate (dα/dt) for each value of α using Equation (6): The terms (dα/dt) α , [A α f (α)] and exp(−E α /RT) of Equation (6) were evaluated for each α. If the reaction rate increased and the term exp(−E α /RT) decreased (i.e., E α increased), then it was concluded that an increase of the term [A α f (α)] compensated for the decrease of the exponential term. This corresponds to a change in the pre-exponential factor and/or to the reaction mechanism. Thus, the aim of this work is to show how the analysis of these variations in association with the E α dependency can be used to identify changes in the reaction mechanism. Figure 1 show the variation of the extent of conversion (α) with temperature (T) for three heating rates and the corresponding reaction rate (dα/dt). Equation (2) is the equation of a single-step process that does not apply to complex polymerization, which is a multi-step process. When applying an isoconversional method, the computations are performed for a constant value of α. Thus, Equation (2) is transformed into Equation (6). In this case, the hypothesis of a single-step process is only applied for each constant α value used for the computation, which corresponds, for non-isothermal data, to a very narrow temperature range. Therefore, it can be assumed that the validity of a single-step equation for a given α value generally holds, even for complex reactions. In addition, the Arrhenius equation is only applied to a narrow temperature region related to this α value.
Molecules 2019, 24, x FOR PEER REVIEW 4 of 16 non-isothermal data, to a very narrow temperature range. Therefore, it can be assumed that the validity of a single-step equation for a given α value generally holds, even for complex reactions. In addition, the Arrhenius equation is only applied to a narrow temperature region related to this α value.

Figure 1.
Non-isothermal data. Example of variation of the extent of conversion (α) with temperature (T). Inset: corresponding variation of the reaction rate (dα/dt) with temperature. Isoconversional methods are based on the assumption of the hypothesis of a single-step process only for each α value and the Arrhenius equation applies to a narrow temperature region related to this α value.
Isoconversional methods require the performance of a series of experiments at different temperature programs and yield the values of effective activation energy Eα as a function of the extent of conversion α. A significant variation of Eα with α indicates that the process is kinetically complex and the Eα dependencies evaluated by an isoconversional method allow for meaningful mechanistic and kinetic analyses and for understanding multi-step processes, as well as for reliable kinetic predictions. Model-free isoconversional kinetic methods are a powerful tool to gain information on the reaction complexity through the Eα dependency determination. A change in the slope of the Eα dependency may generally be associated with a change in the rate-limiting step of the overall reaction mechanism. Isoconversional methods are based on the isoconversional principle that states that the reaction rate is only a function of temperature for a given constant value of the extent of conversion. Thus, the Eα dependency can be evaluated without any assumption of the reaction mechanism, as illustrated by Equation (7): In this equation, the derivative of the term containing f(α) is zero because each computation is performed for a constant value of α (isoconversional methods).

Data Simulation
Two sets of simulated data were generated and analyzed for non-isothermal conditions using the heating rates of 1, 2 and 4 Kmin −1 and for isothermal conditions at four temperatures of 120, 140, 160 and 180 °C. Usually, it is recommended to use three to five temperature programs [6]. In this work only three heating rates were used for the Eα dependency computations because using three or  Isoconversional methods require the performance of a series of experiments at different temperature programs and yield the values of effective activation energy E α as a function of the extent of conversion α. A significant variation of E α with α indicates that the process is kinetically complex and the E α dependencies evaluated by an isoconversional method allow for meaningful mechanistic and kinetic analyses and for understanding multi-step processes, as well as for reliable kinetic predictions. Model-free isoconversional kinetic methods are a powerful tool to gain information on the reaction complexity through the E α dependency determination. A change in the slope of the E α dependency may generally be associated with a change in the rate-limiting step of the overall reaction mechanism. Isoconversional methods are based on the isoconversional principle that states that the reaction rate is only a function of temperature for a given constant value of the extent of conversion. Thus, the E α dependency can be evaluated without any assumption of the reaction mechanism, as illustrated by Equation (7): In this equation, the derivative of the term containing f (α) is zero because each computation is performed for a constant value of α (isoconversional methods).

Data Simulation
Two sets of simulated data were generated and analyzed for non-isothermal conditions using the heating rates of 1, 2 and 4 Kmin −1 and for isothermal conditions at four temperatures of 120, 140, 160 and 180 • C. Usually, it is recommended to use three to five temperature programs [6]. In this work only three heating rates were used for the E α dependency computations because using three or four heating rates led to the same values as the computation was performed using simulated data.
In the first set (set 1), a complex reaction involving an autocatalytic step, a first-order step and a diffusion process was simulated according to the following equations [18]: and using the following parameters [4]: Here k D , k C and k ef respectively represent the specific rate constant for diffusion, the rate constant for the chemically controlled reaction and the effective rate constant. D 0 represents the pre-exponential factor of the diffusion-controlled reaction, A 1 the pre-exponential factor for the non-catalyzed reaction and A 2 the pre-exponential factor for the catalyzed reaction. K is a constant accounting for the effect of the chemical reaction on the change in diffusivity. E D , E 1 and E 2 respectively represent the activation energy of the diffusion-controlled reaction, the activation energy of the non-catalyzed reaction and the activation energy of the catalyzed reaction. m and n are kinetic exponents [19,20]. A second set of data (set 2) involving a first-order step and a diffusion process was simulated according to the same procedure, wherein Equation (9) was replaced by Equation (13): and using the following parameters: Figure 2 shows the variation of the reaction rate and of the extent of conversion with temperature for three heating rates (non-isothermal data). The curves shifted to higher temperatures when the heating rate was increased. Figure 3 shows the variation of the reaction rate and of the extent of conversion with time for four temperatures (isothermal data). It can be observed that the maximum of the reaction rate was not obtained for α = 0, as would be the case for a reaction order mechanism. Another characteristic feature of the autocatalytic model is that the isothermal α-t curves presented typical sigmoidal shapes, as shown in Figure 3.

Dependence of the Effective Activation Energy and of the Pre-Exponential Factor
The dependence of the effective activation energy (Eα) with the extent of conversion (α) is presented in Figure 4 (left axis). Note that FR and NLN methods gave similar results in this case. The complexity of the mechanism is perfectly reflected by the important Eα dependence observed. The first value obtained was Eα = 63.4 kJ·mol −1 for α = 0.02 with the NLN method. This value is close to the value of E1 = 70.0 kJ·mol −1 used in the simulation for the activation energy of the uncatalyzed reaction. For α = 0.001 this value would be Eα = 69.0 kJ·mol −1 which is very close to E1. The lowest activation energy value was Eα = 5.1 kJ·mol −1 for α = 0.98 (NLN method), which is very close to the activation energy of diffusion (ED = 4.4 kJ·mol −1 ). Figure 4 (right axis) give the results obtained for the dependence of the logarithm of the pre-exponential factor (lnAα) with the extent of conversion (α). This dependence was computed using the compensation parameters method described by Sbirrazzuoli [11]. For this computation the

Dependence of the Effective Activation Energy and of the Pre-Exponential Factor
The dependence of the effective activation energy (Eα) with the extent of conversion (α) is presented in Figure 4 (left axis). Note that FR and NLN methods gave similar results in this case. The complexity of the mechanism is perfectly reflected by the important Eα dependence observed. The first value obtained was Eα = 63.4 kJ·mol −1 for α = 0.02 with the NLN method. This value is close to the value of E1 = 70.0 kJ·mol −1 used in the simulation for the activation energy of the uncatalyzed reaction. For α = 0.001 this value would be Eα = 69.0 kJ·mol −1 which is very close to E1. The lowest activation energy value was Eα = 5.1 kJ·mol −1 for α = 0.98 (NLN method), which is very close to the activation energy of diffusion (ED = 4.4 kJ·mol −1 ). Figure 4 (right axis) give the results obtained for the dependence of the logarithm of the pre-exponential factor (lnAα) with the extent of conversion (α). This dependence was computed using the compensation parameters method described by Sbirrazzuoli [11]. For this computation the

Dependence of the Effective Activation Energy and of the Pre-Exponential Factor
The dependence of the effective activation energy (E α ) with the extent of conversion (α) is presented in Figure 4 (left axis). Note that FR and NLN methods gave similar results in this case. The complexity of the mechanism is perfectly reflected by the important E α dependence observed. The first value obtained was E α = 63.4 kJ·mol −1 for α = 0.02 with the NLN method. This value is close to the value of E 1 = 70.0 kJ·mol −1 used in the simulation for the activation energy of the uncatalyzed reaction. For α = 0.001 this value would be E α = 69.0 kJ·mol −1 which is very close to E 1 . The lowest activation energy value was E α = 5.1 kJ·mol −1 for α = 0.98 (NLN method), which is very close to the activation energy of diffusion (E D = 4.4 kJ·mol −1 ).
(three-dimensional diffusion), R3 (contracting sphere) and D2 (two-dimensional diffusion) of ref. [6]. For α = 0.02, ln(Aα/s -1 ) was found to be 9.16, which is very close to the value used in the simulation, i.e., ln(A1/s -1 ) = 9.94. The lowest value for ln(Aα/s -1 ) was −4.82, which can be associated with ln(D0/s -1 ) = 0.36. The dependence of the effective activation energy (Eα) with temperature (T) is presented in Figure 5 in association with the corresponding activation energies of the different steps, i.e., E1, E2 and ED.     Figure 4 (right axis) give the results obtained for the dependence of the logarithm of the pre-exponential factor (lnA α ) with the extent of conversion (α). This dependence was computed using the compensation parameters method described by Sbirrazzuoli [11]. For this computation the model-fitting method of Achar-Brindley-Sharp was used in the interval 0.10 < α < 0.40 to evaluate the relationship between E and lnA. A good correlation (r 2 = 0.99903) between these two parameters was obtained using the models F1 (Mampel, first order), A2 (Avrami-Erofeev), D3 (three-dimensional diffusion), R3 (contracting sphere) and D2 (two-dimensional diffusion) of ref. [6]. For α = 0.02, ln(A α /s −1 ) was found to be 9.16, which is very close to the value used in the simulation, i.e., ln(A 1 /s −1 ) = 9.94. The lowest value for ln(A α /s −1 ) was −4.82, which can be associated with ln(D 0 /s −1 ) = 0.36. The dependence of the effective activation energy (E α ) with temperature (T) is presented in Figure 5 in association with the corresponding activation energies of the different steps, i.e., E 1 , E 2 and E D .

Variation of the Reaction Rate with the Extent of Conversion
The terms [A α f (α)] and exp[−E α /(RT α )] can be evaluated using Equation (6), then a value of the reaction rate can be computed for each extent of conversion as the product of [A α f (α)] by exp[−E α /(RT α )]. Figure 6 shows the relation between the dependence of the effective activation energy (E α ) with the extent of conversion (α) and the reaction rate (dα/dt). The maximum of the reaction rate (dα/dt) was located in the range of 0.46 < α < 0.48. Figure 6 shows that this corresponds to a change in the slope of the E α dependence.  Figure 6 shows the relation between the dependence of the effective activation energy (Eα) with the extent of conversion (α) and the reaction rate (dα/dt). The maximum of the reaction rate (dα/dt) was located in the range of 0.46 < α < 0.48. Figure 6 shows that this corresponds to a change in the slope of the Eα dependence.  (6)).
The comparison of the terms [Aα f(α)], exp[−Eα/(RTα)] and (dα/dt) can be used to identify rate-limiting steps in the overall reaction rate. The principle of this method is explained below and illustrated by Figure 7. From α = 0.02 to 0.46, (dα/dt) increases and the term [Aα f(α)] decreases, so the increase of (dα/dt) is attributed to the increase of the term exp[−Eα/(RTα)], i.e., a decrease of Eα. Thus, the increase of the rate is mainly due to a favorable energetic term. When α ≥ 0.48, (dα/dt) decreases, the term exp[−Eα/(RTα)] still increases and the term [Aα f(α)] still decreases. This indicates that the term [Aα f(α)] dominates and corresponds to a change in the rate-limiting step at this stage of the reaction. The decrease of the rate is attributed to a change in the mechanism, which may originate from a change in f(α) or from an entropic change (Aα). This entropic change may be due to a change in configuration or a decrease in the efficiency of collisions. Thus, it is identified as a change in the rate-limiting step in the overall reaction rate for α = 0.46-0.48. Table 1 gives some values of the various terms of Figure 7 used to identify a change in the rate-limiting step for 0.46 < α < 0.48.  (6)).
The comparison of the terms [A α f (α)], exp[−E α /(RT α )] and (dα/dt) can be used to identify rate-limiting steps in the overall reaction rate. The principle of this method is explained below and illustrated by Figure 7. From α = 0.02 to 0.46, (dα/dt) increases and the term [A α f (α)] decreases, so the increase of (dα/dt) is attributed to the increase of the term exp[−E α /(RT α )], i.e., a decrease of E α . Thus, the increase of the rate is mainly due to a favorable energetic term. When α ≥ 0.48, (dα/dt) decreases, the term exp[−E α /(RT α )] still increases and the term [A α f (α)] still decreases. This indicates that the term [A α f (α)] dominates and corresponds to a change in the rate-limiting step at this stage of the reaction. The decrease of the rate is attributed to a change in the mechanism, which may originate from a change in f (α) or from an entropic change (A α ). This entropic change may be due to a change in configuration or a decrease in the efficiency of collisions. Thus, it is identified as a change in the rate-limiting step in the overall reaction rate for α = 0.46-0.48. Table 1 gives some values of the various terms of Figure 7   Once the change in the rate-limiting step was identified for α = 0.46, an analysis of Figure 4 showed that ln(Aα/s -1 ) = 5.11 and Eα = 44.0 kJ·mol −1 for α = 0.46, which are close to the values used in the simulation for the catalyzed reaction, i.e., ln(A2/s -1 ) = 6.21 and E2 = 45.0 kJ·mol −1 . The change in the curvature of the Eα dependence occurred exactly in this region of α = 0.46-0.48, as seen in Figures 4, 5 and 6. This result shows that, in addition to giving information on the rate-limiting steps, the Eα and ln(Aα) dependencies can be used as estimate values of kinetic parameters to be used in a non-linear fitting procedure [9]. Figure 8 shows the variation of the rate coefficients kD, kC and kef with the extent of conversion (α) kD always decreased with α, while kC always increased. Initial values of kD were high, while they were low for kC. This is the opposite at the end where the values of kD were lower than the values of kC. This is in perfect agreement with a chemical control at the beginning of the reaction (kC << kD) followed by a diffusion control at the end (kD << kC). The variation of kef is more complex, as reflected by Equation (10), showing an increasing trend at the beginning and a decreasing tendency at the end. The previously reported value of α = 0.46-0.48 (T~183 °C), attributed to a change in the rate-limiting step, corresponded to the point at which kC started to become higher than kef and kD/kC < 10 (9.6). Generally, it is estimated that a factor of 100 is the minimum required to neglect one reaction to another. This value was obtained for α = 0.24 (T~173 °C), kD/kC ≈ 100 (102.7).  Table 1. Some values of the various terms used to identify a change in the rate-limiting step. Once the change in the rate-limiting step was identified for α = 0.46, an analysis of Figure 4 showed that ln(A α /s −1 ) = 5.11 and E α = 44.0 kJ·mol −1 for α = 0.46, which are close to the values used in the simulation for the catalyzed reaction, i.e., ln(A 2 /s −1 ) = 6.21 and E 2 = 45.0 kJ·mol −1 . The change in the curvature of the E α dependence occurred exactly in this region of α = 0.46-0.48, as seen in Figures 4-6. This result shows that, in addition to giving information on the rate-limiting steps, the E α and ln(A α ) dependencies can be used as estimate values of kinetic parameters to be used in a non-linear fitting procedure [9]. Figure 8 shows the variation of the rate coefficients k D , k C and k ef with the extent of conversion (α) k D always decreased with α, while k C always increased. Initial values of k D were high, while they were low for k C . This is the opposite at the end where the values of k D were lower than the values of k C . This is in perfect agreement with a chemical control at the beginning of the reaction (k C << k D ) followed by a diffusion control at the end (k D << k C ). The variation of k ef is more complex, as reflected by Equation (10), showing an increasing trend at the beginning and a decreasing tendency at the end. The previously reported value of α = 0.46-0.48 (T~183 • C), attributed to a change in the rate-limiting step, corresponded to the point at which k C started to become higher than k ef and k D /k C < 10 (9.6). Generally, it is estimated that a factor of 100 is the minimum required to neglect one reaction to another. This value was obtained for α = 0.24 (T~173 • C), k D /k C ≈ 100 (102.7).

Variation of the Overall Rate Coefficient k(T) and of the Effective Rate Coefficient kef(T) with Reciprocal Temperature
The variations of ln k(T) and ln kef (T) with reciprocal temperature are presented in Figure 9. The values of α = 0.46-0.48 corresponded to the temperature at which the two rate constants were very close. The closest values were reached for α = 0.54.

Fit of the Eα Dependence with the Sourour and Kamal and Diffusion Models
The isoconversional principle (Equation (7)) was applied to the autocatalytic equation of Sourour and Kamal (Equation (9)). This model is used to describe the initial stages of the reaction when it is chemically controlled [4,5,18]: with temperature, red squares: variation of k D with temperature, green lozenges: variation of k C with temperature, magenta circles: variation of k ef with temperature.

Variation of the Overall Rate Coefficient k(T) and of the Effective Rate Coefficient k ef (T) with Reciprocal Temperature
The variations of ln k(T) and ln k ef (T) with reciprocal temperature are presented in Figure 9. The values of α = 0.46-0.48 corresponded to the temperature at which the two rate constants were very close. The closest values were reached for α = 0.54.

Variation of the Overall Rate Coefficient k(T) and of the Effective Rate Coefficient kef(T) with Reciprocal Temperature
The variations of ln k(T) and ln kef (T) with reciprocal temperature are presented in Figure 9. The values of α = 0.46-0.48 corresponded to the temperature at which the two rate constants were very close. The closest values were reached for α = 0.54.

Fit of the Eα Dependence with the Sourour and Kamal and Diffusion Models
The isoconversional principle (Equation (7)) was applied to the autocatalytic equation of Sourour and Kamal (Equation (9)). This model is used to describe the initial stages of the reaction when it is chemically controlled [4,5,18]:

Fit of the E α Dependence with the Sourour and Kamal and Diffusion Models
The isoconversional principle (Equation (7)) was applied to the autocatalytic equation of Sourour and Kamal (Equation (9)). This model is used to describe the initial stages of the reaction when it is chemically controlled [4,5,18]: The isoconversional principle (Equation (7)) was also applied to the diffusion model (Equation (8)) used to describe the end of the reaction [4,5,18]: Note that, according to Equation (15), at the lowest extent of conversion (α → 0) E α tends toward the activation energy of the uncatalyzed reaction (E α → E 1 ), in agreement to what was reported in the analysis of Figure 4.
The results of the fit of Equations (15) and (17) are given in Table 2. Some differences were observed between the simulated values and the values resulting from the non-linear fit. A higher discrepancy was obtained for A 1 /A 2 (5918.03). This value should be 41.56. Nevertheless, a good agreement was found between the simulated and experimental data resulting from the non-linear fit for the other parameters. The result of this fit is presented in Figure 10.   (18) and (19).
Equations (15) and (17) can be expanded to allow computations of the pre-exponential factors, as proposed by Sbirrazzuoli in [4]. The resulting Equations (18) and (19) can be fitted to estimate A 1 , A 2 and D 0 . The results are given in Table 3.
It can be seen that the parameters are in very good agreement with the reference values and the values of A 1 , A 2 are very close the reference values in this case. For the initial part of the reaction, the restriction of the fit to the interval 0.02 < α < 0.24 resulted in an improved accuracy (lower MSSD). This confirms the previous statement that when k D /k C ≈ 100 it is possible to neglect the diffusion reaction, while for higher temperatures (or values of α) it is not completely negligible. Although the fit was better for parameters of Table 3 in comparison with those of Table 2, it is impossible to see the difference on the graph (inset of Figure 10). For the autocatalytic model of Sourour and Kamal, the accuracy of the fit was improved by addition of more flexibility when moving from Equation (15) to Equation (18). Nevertheless, this increased the possibilities of reaching local minima, resulting in several sets of parameters leading to an accurate fit of the data. The use of parameters estimated by the advanced isoconversional method, as initial values of the non-linear fit, greatly facilitated the achievement of meaningful parameters and not only fitting parameters, especially for the initial part of the reaction. However, the existence of local minima is a problem that cannot be underestimated and that is difficult to avoid when fitting complex mechanisms which involve many kinetic parameters to be determined using non-linear fits. This is less problematic for the end of the reaction, i.e., for the fit of Equation (8). The use of a genetic algorithm could be an efficient method to avoid this kind of problem [21]. This confirms the previous statement that when kD/kC ≈ 100 it is possible to neglect the diffusion reaction, while for higher temperatures (or values of α) it is not completely negligible. Although the fit was better for parameters of Table 3 in comparison with those of Table 2, it is impossible to see the difference on the graph (inset of Figure 10). For the autocatalytic model of Sourour and Kamal, the accuracy of the fit was improved by addition of more flexibility when moving from Equation (15) to Equation (18). Nevertheless, this increased the possibilities of reaching local minima, resulting in several sets of parameters leading to an accurate fit of the data. The use of parameters estimated by the advanced isoconversional method, as initial values of the non-linear fit, greatly facilitated the achievement of meaningful parameters and not only fitting parameters, especially for the initial part of the reaction. However, the existence of local minima is a problem that cannot be underestimated and that is difficult to avoid when fitting complex mechanisms which involve many kinetic parameters to be determined using non-linear fits. This is less problematic for the end of the reaction, i.e., for the fit of Equation (8). The use of a genetic algorithm could be an efficient method to avoid this kind of problem [21].  Table 2. Inset: same plot for the data of Table 3. Figure 11 shows the variation of the reaction rate and of the extent of conversion with temperature for three heating rates. In comparison with what was observed for data set 1, which include an autocatalytic step (Figure 2), the shift to higher temperatures upon increasing the heating rate was much lower in this case. This shows the difference between the reaction order and the autocatalytic mechanism for non-isothermal data. Figure 12 shows the variation of the reaction rate and of the extent of conversion with time for four temperatures (isothermal data). As expected, the maximum of the reaction rate was obtained for α = 0 in this case and the isothermal α-t curves presented the characteristic shapes of a reaction order model.  Table 2. Inset: same plot for the data of Table 3. Figure 11 shows the variation of the reaction rate and of the extent of conversion with temperature for three heating rates. In comparison with what was observed for data set 1, which include an autocatalytic step (Figure 2), the shift to higher temperatures upon increasing the heating rate was much lower in this case. This shows the difference between the reaction order and the autocatalytic mechanism for non-isothermal data. Figure 12 shows the variation of the reaction rate and of the extent of conversion with time for four temperatures (isothermal data). As expected, the maximum of the reaction rate was obtained for α = 0 in this case and the isothermal α-t curves presented the characteristic shapes of a reaction order model.  The dependence of the effective activation energy (Eα) with extent of conversion (α) is presented in Figure 13 (left axis) and the dependence with temperature is presented in Figure 14. Note that FR and NLN methods gave similar results in this case. The complexity of the mechanism is perfectly reflected by the important Eα dependence observed. The first value obtained was around 50 kJ·mol −1 , which is in perfect agreement with the value of the activation energy of the reaction order reaction E1 (50.0 kJ·mol −1 ). Figure 13 (right axis) gives the results obtained for the dependence of the logarithm of the pre-exponential factor (lnAα/s −1 ) with the extent of conversion (α). ln(Aα/s −1 ) was found to be 6.61 and Eα = 50.2 kJ·mol −1 for α = 0.02, which are close to the values used in the simulation, i.e., ln(A1/s -1 ) = 7.47 and E1 = 50.0 kJ·mol −1 . For α = 0.98, ln(Aα/s −1 ) was found to be −5.65 and Eα = 5.5 kJ·mol −1 , which are also close to the values used in the simulation, i.e., ln(D0/s -1 ) = 0.36 and ED = 4.4 kJ·mol −1 .

First-Order Reaction with Diffusion-Controlled Part (Data Set 2)
Application of the compensation parameters method [11] in the interval 0.05 < α < 0.25 (Achar-Brindley-Sharp's differential method, models F1, A3, A2, D3, r 2 = 0.999999) permit the identification of the F1 model (Mampel, first order) and the kinetic parameters of the reaction order  The dependence of the effective activation energy (Eα) with extent of conversion (α) is presented in Figure 13 (left axis) and the dependence with temperature is presented in Figure 14. Note that FR and NLN methods gave similar results in this case. The complexity of the mechanism is perfectly reflected by the important Eα dependence observed. The first value obtained was around 50 kJ·mol −1 , which is in perfect agreement with the value of the activation energy of the reaction order reaction E1 (50.0 kJ·mol −1 ). Figure 13 (right axis) gives the results obtained for the dependence of the logarithm of the pre-exponential factor (lnAα/s −1 ) with the extent of conversion (α). ln(Aα/s −1 ) was found to be 6.61 and Eα = 50.2 kJ·mol −1 for α = 0.02, which are close to the values used in the simulation, i.e., ln(A1/s -1 ) = 7.47 and E1 = 50.0 kJ·mol −1 . For α = 0.98, ln(Aα/s −1 ) was found to be −5.65 and Eα = 5.5 kJ·mol −1 , which are also close to the values used in the simulation, i.e., ln(D0/s -1 ) = 0.36 and ED = 4.4 kJ·mol −1 .
Application of the compensation parameters method [11] in the interval 0.05 < α < 0.25 (Achar-Brindley-Sharp's differential method, models F1, A3, A2, D3, r 2 = 0.999999) permit the identification of the F1 model (Mampel, first order) and the kinetic parameters of the reaction order The dependence of the effective activation energy (E α ) with extent of conversion (α) is presented in Figure 13 (left axis) and the dependence with temperature is presented in Figure 14. Note that FR and NLN methods gave similar results in this case. The complexity of the mechanism is perfectly reflected by the important E α dependence observed. The first value obtained was around 50 kJ·mol −1 , which is in perfect agreement with the value of the activation energy of the reaction order reaction E 1 (50.0 kJ·mol −1 ). Figure 13 (right axis) gives the results obtained for the dependence of the logarithm of the pre-exponential factor (lnA α /s −1 ) with the extent of conversion (α). ln(A α /s −1 ) was found to be 6.61 and E α = 50.2 kJ·mol −1 for α = 0.02, which are close to the values used in the simulation, i.e., ln(A 1 /s −1 ) = 7.47 and E 1 = 50.0 kJ·mol −1 . For α = 0.98, ln(A α /s −1 ) was found to be −5.65 and E α = 5.5 kJ·mol −1 , which are also close to the values used in the simulation, i.e., ln(D 0 /s −1 ) = 0.36 and E D = 4.4 kJ·mol −1 .
In the case of two reactions with similar activation energies, the identification of the contribution of each individual reaction using the Eα dependency could be more difficult. Nevertheless, it is highly probable that these reactions would have different pre-exponential factors. In this case, we proposed to identify the complexity of the overall process by analyzing the Aα dependency computed with the compensation method.  Figure 13. Dependence of the effective activation energy (Eα) and pre-exponential factor (lnAα) with the extent of conversion (α). Non-isothermal data. Open lozenges: pre-exponential factor computed using compensation parameters method, red lozenges: Eα computed with FR method, green circles: Eα computed with NLN method. Inset: isothermal data. Red triangles: Eα computed with FR method, green stars: Eα computed with NLN method.

Conclusions
A kinetic model has been proposed that simulates complex polymerizations with high accuracy. The first stages of the reaction were described by an autocatalytic process, followed by epoxy-amine addition and a diffusion-controlled model at the end of the reaction. Isoconversional Figure 13. Dependence of the effective activation energy (E α ) and pre-exponential factor (lnA α ) with the extent of conversion (α). Non-isothermal data. Open lozenges: pre-exponential factor computed using compensation parameters method, red lozenges: E α computed with FR method, green circles: E α computed with NLN method. Inset: isothermal data. Red triangles: E α computed with FR method, green stars: E α computed with NLN method. used in the simulation were ln(A1/s −1 ) = 7.47 and E1 = 50.0 kJ·mol −1 .
In the case of two reactions with similar activation energies, the identification of the contribution of each individual reaction using the Eα dependency could be more difficult. Nevertheless, it is highly probable that these reactions would have different pre-exponential factors. In this case, we proposed to identify the complexity of the overall process by analyzing the Aα dependency computed with the compensation method.  Figure 13. Dependence of the effective activation energy (Eα) and pre-exponential factor (lnAα) with the extent of conversion (α). Non-isothermal data. Open lozenges: pre-exponential factor computed using compensation parameters method, red lozenges: Eα computed with FR method, green circles: Eα computed with NLN method. Inset: isothermal data. Red triangles: Eα computed with FR method, green stars: Eα computed with NLN method.

Conclusions
A kinetic model has been proposed that simulates complex polymerizations with high accuracy. The first stages of the reaction were described by an autocatalytic process, followed by epoxy-amine addition and a diffusion-controlled model at the end of the reaction. Isoconversional Application of the compensation parameters method [11] in the interval 0.05 < α < 0.25 (Achar-Brindley-Sharp's differential method, models F1, A3, A2, D3, r 2 = 0.999999) permit the identification of the F1 model (Mampel, first order) and the kinetic parameters of the reaction order reaction. The values found were ln(A/s −1 ) = 7.45 and E = 49.9 kJ·mol −1 , while the kinetic parameters used in the simulation were ln(A 1 /s −1 ) = 7.47 and E 1 = 50.0 kJ·mol −1 .
In the case of two reactions with similar activation energies, the identification of the contribution of each individual reaction using the E α dependency could be more difficult. Nevertheless, it is highly probable that these reactions would have different pre-exponential factors. In this case, we proposed to identify the complexity of the overall process by analyzing the A α dependency computed with the compensation method.

Conclusions
A kinetic model has been proposed that simulates complex polymerizations with high accuracy. The first stages of the reaction were described by an autocatalytic process, followed by epoxy-amine addition and a diffusion-controlled model at the end of the reaction. Isoconversional methods-based on the assumption of the hypothesis of a single-step process only for each α value and the application of the Arrhenius equation to a very narrow temperature region related to this α value-give important insights into the change in the rate-limiting steps by analysis of the E α dependency and its variations. In addition to this, the comparisons of the terms (dα/dt), [A α f (α)] and exp[−E α /(RT α )] evaluated by Friedman's method can help to identify the change in the rate-limiting steps of the overall mechanism as measured by thermoanalytical techniques. It was also concluded that the assumption of the validity of a single-step equation, when restricted to a given α value, holds for complex reactions.