Problems with Applying the Ozawa–Avrami Crystallization Model to Non-Isothermal Crosslinking Polymerization

Ozawa has modified the Avrami model to treat non-isothermal crystallization kinetics. The resulting Ozawa–Avrami model yields the Avrami index (n) and heating/cooling function (χ(T)). There has been a number of recent applications of the Ozawa–Avrami model to non-isothermal crosslinking polymerization (curing) kinetics that have determined n and have used χ(T) in place of the rate constant (k(T)) in the Arrhenius equation to evaluate the activation energy (E) and the preexponential factor (A). We analyze this approach mathematically as well as by using simulated and experimental data, highlighting the following problems. First, the approach is limited to the processes that obey the Avrami model. In cases of autocatalytic or decelerating kinetics, commonly encountered in crosslinking polymerizations, n reveals a systematic dependence on temperature. Second, χ(T) has a more complex temperature dependence than k(T) and thus cannot produce exact values of E and A via the Arrhenius equation. The respective deviations can reach tens or even hundreds of percent but are diminished dramatically using the heating/cooling function in the form [χ(T)]1/n. Third, without this transformation, the Arrhenius plots may demonstrate breakpoints that leads to questionable interpretations. Overall, the application of the Ozawa–Avrami model to crosslinking polymerizations appears too problematic to be justified, especially considering the existence of well-known alternative kinetic techniques that are flexible, accurate, and computationally simple.


Introduction
The Avrami model (also known as the Johnson-Mehl-Avrami-Kolmogorov model [1,2]) has been developed to describe the kinetics of the crystalline phase formation via the nucleation mechanism. The most common application area of the model is the kinetics of crystallization [1,2]. The classical Avrami model applies to isothermal conditions and is typically used in the following form: where t is the time, α is the fractional volume transformed into the crystalline phase, n is the Avrami exponent (or index), and k A (T) is the Avrami rate constant. An essential feature of the model is that it imitates the so-called sigmoid kinetics, i.e., the situation when the process rate initially increases, reaches a maximum, and finally decreases. Respectively, the cumulative growth of the new phase demonstrates a sigmoid profile with respect to time. It is noteworthy that apart from crystallization, there are a multitude of processes that demonstrate the sigmoid kinetics. This fact has inspired a variety of applications of the Avrami model that have stretched far beyond the process of crystallization. Some examples include: gelation [3][4][5], the adsorption of solutes [6][7][8] and

Avrami and Ozawa-Avrami Models
One of the most common techniques used to monitor crystallization is differential scanning calorimetry (DSC) and/or differential thermal analysis (DTA). The application is based on the fact that the volume of the crystalline phase that is formed is proportional to the heat released during crystallization.
In kinetic terms, it is assumed that the rate of crystallization is directly proportional to the heat flow dQ/dt: (2) where Q 0 is the total heat released, k(T) is the Arrhenius rate constant, E is the activation energy, A is the preexponential factor, R is the gas constant, and f (α) is the reaction model. Note that because the heat flow is not species specific, the respective rate generally has an overall nature, i.e., it can correspond to more than a single reaction step. For this reason, the parameters of Equation (2) also generally have an overall nature. The aforementioned proportionality between the rate and heat flow was introduced in an early work by Borchardt and Daniels [33]. Their analysis of the heat flow components as measured by DTA or heat flux DSC suggests that one can neglect the so-called thermal inertia term. This, however, may not always be justified. In particular, when one uses larger sample masses and faster heating rates, neglecting the thermal inertia term can result in significant systematic errors in the activation energy [34,35]. In this regard, it is important to stress that accounting for thermal inertia can be accomplished via a relatively simple procedure [35].
Per Equation (1), the extent of conversion is determined experimentally as the fractional area of the DSC peak. Then, the parameters of the Avrami model can be estimated from Equation (3) ln by plotting its left-hand side against the lnt. Assuming that k A (T) has an Arrhenian temperature dependence permits evaluating the Arrhenius parameters, i.e., E and lnA [36]. However, the direct substitution of k A (T) into the Arrhenius equation as follows: creates some problems with correctly estimating its parameters. The issue has been addressed in several publications [37][38][39][40][41]. Briefly, one needs to recognize that the rate in Equation (2) invariably has the units of (time) −1 . Since f (α) is dimensionless, k(T) also has to have the units of (time) −1 . However, estimating k A (T) from the Avrami equation (Equations (1) or (3)) yields the value that has the units of (time) -n . For that reason, it should not be substituted directly into the Arrhenius Equation (4). The problem is resolved by a simple transformation that converts the units of k A (T) into (time) −1 : Without this transformation, the direct substitution of k A (T) in Equation (4) results in estimating apparent Arrhenius parameters that are n times larger than the correct values. That is: [37,39,41] The problem is avoided altogether if the Avrami equation is used in the modified form: [37][38][39][40][41][42] 1 The resulting Avrami rate constant, k mA always has the units of (time) −1 , and thus, its direct substitution into Equation (4) yields the correct values of E and lnA. Overall, depending on the form of the Avrami equation used, the correct values of E and lnA are determined by defining the Arrhenius rate constant as follows: Ozawa has developed a method of adjusting the Avrami model for non-isothermal conditions [20]. The basic equation of the resulting Ozawa-Avrami model is: where β is the cooling or heating rate and χ(T) is, respectively, the cooling or heating function. The parameter n is sometimes confusingly called the Ozawa exponent, but in fact is the Avrami exponent [43]. The method makes use of the isothermal sections of the non-isothermal data. That is, one selects temperature and then determines what values of α correspond to this temperature at different values of β. Then, n and χ(T) are readily evaluated using Equation (10): by plotting the left-hand side against lnβ. Although Ozawa has never mentioned that χ(T) can be applied to estimating the Arrhenius parameters, such an application has been used in numerous publications devoted to the kinetic analysis of crosslinking polymerization [21][22][23][24][25][26][27]. In essence, they have proposed substituting χ(T) for k(T) in the Arrhenius Equation (4) to determine E and lnA. In our opinion, this is a questionable proposition. The first clue of this being in question is that χ(T) has units of (heating or cooling rate) −n , whereas k(T) of (time) −1 , i.e., these are two different physical quantities. The actual relationship between χ(T) and k(T) is established to be: [44][45][46][47] d The derivations of Equation (11) are not presented here as they are readily available in the original publications [44][45][46][47]. Note that the negative sign in the right-hand side of Equation (11) applies only to the conditions of cooling, i.e., when β < 0. This is because on cooling χ(T) is a decreasing function of temperature [20], so its derivative is negative, whereas k(T) must always be positive. For the conditions of heating, the derivative is positive, so the negative sign is dropped. Clearly, χ(T) is nothing else but the rate constant integrated over temperature. Perhaps some empirical support for using χ(T) in the Arrhenius equation comes from the fact that lnχ(T) tends to show a reasonably linear dependence on T −1 and yields the E values rather similar to the activation energies. These issues are analyzed further by using some simulated data.

Simulations
To look closer into the use of χ(T) for estimating the Arrhenius parameters as well as the Ozawa-Avrami model in general, we simulate processes of two types: one that obeys the Avrami model with n = 3 and another that follows the second-order kinetics model. These models are commonly denoted as A3 and F2, respectively. Note that the reaction-order kinetics are routinely encountered in crosslinking polymerizations as seen in examples publicized elsewhere [28][29][30]. Both processes have the same Arrhenius parameters: E = 100 kJ mol −1 and A = 10 12 min −1 (i.e., lnA = 27.63) and are simulated in the form of α vs. T curves at five heating rates: 1, 1.5, 2, 3, and 5 K min −1 . This is accomplished by solving the following equation for α: 3 for the A3 model or (1 − α) −1 − 1 for the F2 model [31,39,48]. The temperature integral is solved with the aid of a highly accurate approximation [49].

Analysis of Simulated A3 Data
The simulated kinetic curves are displayed in Figure 1. Figure 2 shows a set of the Ozawa plots (Equation (10)) obtained from the simulated A3 data. As expected, the plots are perfectly linear (the coefficient of linear correlation in all cases is 1) and yield the exact values of n (Table 1). These plots also permit estimating the lnχ(T) values. Indeed, the Arrhenius plots of lnχ(T) vs. T −1 demonstrate excellent linearity ( Figure 3). The activation energy estimated directly from such a plot is 319.5 kJ mol −1 , which is more than three times larger than the value used in simulations. The source of this large difference is probably the same as in the case of using lnk A (T) instead of lnk A (T)/n (see Equation (6)). Expectedly, plotting lnχ(T)/n vs. T −1 yields a three times smaller value, E = 106.5 kJ mol −1 . It may seem to be reasonably close to the simulated value. The respective systematic error is 6.5%, which could be considered acceptable in the case of experimental data but not for the case of simulated data. For lnA, the error is even larger: 16%. At any rate, for the simulated data, we have to be obtaining the exact values. For comparison, the estimated n values are exact, at least to the third decimal place (Table 1). That is, the error is less than 0.1%. Thus, the errors in E and lnA should be expected to be on the level of a fraction of a percent.         Of course, the question that arises is "Why are we getting linear Arrhenius plots for the χ(T) although it represents a quantity that differs entirely from k(T)?" The basic reason is quite simple. The integration of the exponential function in k(T) (see Equation (11)) would necessarily include a similar exponential function in χ(T), which simply is a property of the exponential function. To illustrate the point, consider the simplest case of n = 1. Then, for the conditions of heating, Equation (11) can be transformed into The temperature integral in Equation (13) does not have an analytical solution. Yet, it has a number of accurate approximations [50], the simplest of which is as follows: [51] As we can see, the integrated form, χ(T), does include the same exponential function as the Arrhenius rate constant, k(T). Taking the logarithm of both sides of Equation (14), we obtain: Let us now substitute the simulated (true) values of E = 100 kJ mol −1 and A = 10 12 min −1 into Equation (15) and plot the lnχ(T) vs. T −1 dependence for the temperature range 394-414 K (i.e., the same temperature range that we used for the Ozawa analysis (Table 1)). In Figure 3, we can see that this dependence is practically the same as the dependence of lnχ(T)/n vs. T −1 plotted by using the data from Table 1. This obviously indicates that Equation (15) accurately imitates the behavior of the lnχ(T)/n vs. T −1 dependence obtained from the Ozawa analysis (Table 1). Indeed, fitting the lnχ(T) vs. T −1 dependence obtained by Equation (15) to a straight line yields practically the same values of the intercept and E as the corresponding values obtained by fitting the lnχ(T)/n vs. T −1 data from Table 1 to the Arrhenius equation (Equation (4); see Figure 3). The resulting absolute value of the coefficient of linear correlation for the Arrhenius plot set by Equation (15) is 1, just as for the two other plots in Figure 3.
In general, Equation (15) suggests that a strong linearity for lnχ(T)/n vs. T −1 data should be expected as long as the 2lnT term does not vary much with the temperature. In the temperature range 394-414 K, this term increases by only 0.8%, which is practically negligible. By taking the derivative of Equation (15) with respect to T −1 , one can easily find that the E derived from χ(T) is larger than the true value by 2RT. That is, the relative systematic deviation of the activation energy estimated with the plot lnχ(T)/n vs. T −1 from the true value should increase for processes occurring at a higher temperature and with a lower activation energy. Furthermore, it should be stressed that one should not stretch the predictions of Equation (15) too far because the approximation used (Equation (14)) is reasonably accurate within the range 28 < E/RT < 50 [50].
To further illustrate the inaccuracy of using the χ(T) value for estimating the Arrhenius parameters, we treat the simulated A3 data with an advanced isoconversional method [52]. In it, the activation energy E α is found at multiple conversions with the step ∆α (usually 0.01 or 0.02). As a result, one obtains a dependence of E α on α. For each α, E α is found as a minimum of the function: where the integral: is evaluated by using the trapezoid rule. The preexponential factor is then estimated by using the compensation effect equation [31,53]. The resulting dependencies of E α and lnA α on α have demonstrated the constancy of both values, as expected for a single-step process. The mean values are found to be E = 100.05 kJ mol −1 and lnA = 27.59. They deviate, respectively, from the exact values by 0.05 and 0.1%, as opposed to the 6.5 and 16% deviations found for E and lnA determined from using χ(T).
The above examples clearly indicate that the expected accuracy of estimating E and lnA from using the model data should, respectively, be on the scale of a fraction of a percent, i.e., definitely not 6.5 and 16%.

Analysis of Simulated F2 Data
The simulated kinetic curves are displayed in Figure 4. Figure 5 presents a set of the Ozawa plots (Equation (10)) obtained from the simulated F2 data. The plots are nearly perfectly linear. The lowest coefficient of linear correlation is 0.997. An obvious difference between these plots and the ones for A3 is that the slope appears to change with temperature. Indeed, the values of n decrease systematically with temperature ( Table 2). This is a clear indication that the Avrami model is inapplicable to the data that do not represent the sigmoid kinetics. It is noteworthy that a similar decrease in n with increasing temperature is found in experimental data on non-isothermal crosslinking [21,23,27], which thus may be a sign of the inapplicability of the Avrami model. Naturally, if the Avrami model does not apply to the data, one faces the problem of the applicability of the Ozawa-Avrami analysis altogether. the sigmoid kinetics. It is noteworthy that a similar decrease in n with incre ature is found in experimental data on non-isothermal crosslinking [21,23,27 may be a sign of the inapplicability of the Avrami model. Naturally, if the A does not apply to the data, one faces the problem of the applicability of th rami analysis altogether.     Figure 6 shows the Arrhenius plots built by employing the heating func ble 2). We should take notice that lnχ(T) vs. T −1 is distinctly nonlinear an described as a combination of two linear segments or, to put it differently, as plot with a breakpoint. Such a situation is commonly interpreted as eviden   Figure 6 shows the Arrhenius plots built by employing the heating function data ( Table 2). We should take notice that lnχ(T) vs. T −1 is distinctly nonlinear and thus can be described as a combination of two linear segments or, to put it differently, as an Arrhenius plot with a breakpoint. Such a situation is commonly interpreted as evidence of the reaction mechanism change. Of course, in the case of our simulated data, we know perfectly well Polymers 2022, 14, 693 9 of 13 that the process obeys a single mechanism with a single set of the Arrhenius parameters. Needless to say, the activation energy values, 89 and 50 kJ mol −1 , estimated from the two linear segments (see Figure 6) deviate significantly from the simulated value 100 kJ mol −1 . The importance of this example is that it clearly indicates that a breakpoint observed in the lnχ(T) vs. T −1 plot can be a trivial computational artifact. The latter can simply be due to the inapplicability of the Avrami model and can thus have no physical significance. Again, we should notice that the breakpoints have been customarily seen [21,23,[25][26][27] in experimental lnχ(T) vs. T −1 plots for crosslinking.
ble 2). We should take notice that lnχ(T) vs. T −1 is distinctly nonlinear and thus described as a combination of two linear segments or, to put it differently, as an Ar plot with a breakpoint. Such a situation is commonly interpreted as evidence of t tion mechanism change. Of course, in the case of our simulated data, we know p well that the process obeys a single mechanism with a single set of the Arrhenius eters. Needless to say, the activation energy values, 89 and 50 kJ mol −1 , estimated f two linear segments (see Figure 6) deviate significantly from the simulated valu mol −1 . The importance of this example is that it clearly indicates that a breakpoint o in the lnχ(T) vs. T −1 plot can be a trivial computational artifact. The latter can sim due to the inapplicability of the Avrami model and can thus have no physical signi Again, we should notice that the breakpoints have been customarily seen [21,23,2 experimental lnχ(T) vs. T −1 plots for crosslinking. As found in the case of the simulated A3 data, the Arrhenius parameters are better estimated by using the lnχ(T)/n vs. T −1 plots. The corresponding plot is depicted in Figure 6. It is remarkable that the aforementioned breakpoint does not appear in this plot. This is yet another clue that the breakpoint seen in the lnχ(T) vs. T −1 plot can be no more than a computational artifact. The plot is linear (the coefficient of linear correlation 0.9978), as expected. It yields the following Arrhenius parameters: E = 96.9 kJ mol −1 and lnA = 28.8. The respective deviations from the correct values are 3 and 10%. On the other hand, an advanced isoconversional method yields E = 100.2 kJ mol −1 and lnA = 27.88, which, respectively, deviate from the correct values by only 0.2 and 0.9%.

Analysis of Experimental Data
As an experimental example, we employ our previously published data on the crosslinking polymerization kinetics of cyanate ester based on a dimer of 4-tert-butylphenol [54]. Figure 7 shows the DSC curves representing this process. The details of the experiments and kinetic treatment are provided in a previous publication [54]. A comprehensive kinetic analysis has demonstrated that the process follows single-step kinetics. The Arrhenius parameters have been determined to be: E = 85 kJ mol −1 and ln(A/min −1 ) = 18.72. The kinetics of the process have been established to be autocatalytic, i.e., of the sigmoid type. The actual model is f(α) = α 1.20 (1 − α) 1.33 . As the kinetics are sigmoid, the Avrami model can be used as a reasonable approximation.
The Ozawa plots are displayed in Figure 8. The points appear to fall on the straight lines reasonably well. The n and χ(T) values derived from the fit lines are collected in Table 3. It is seen that the n values decrease systematically with an increasing temperature. Although the trend is not nearly as strong as in the case of the F2 data (Table 2), this still may be indicative that the Avrami model is not a proper representation of the process under study. [54]. Figure 7 shows the DSC curves representing this process. The details of the experiments and kinetic treatment are provided in a previous publication [54]. A comprehensive kinetic analysis has demonstrated that the process follows single-step kinetics. The Arrhenius parameters have been determined to be: E = 85 kJ mol −1 and ln(A/min -1 ) = 18.72. The kinetics of the process have been established to be autocatalytic, i.e., of the sigmoid type. The actual model is f(α) = α 1.20 (1-α) 1.33 . As the kinetics are sigmoid, the Avrami model can be used as a reasonable approximation.  Table 3. It is seen that the n values decrease systematically with an increasing temperature. Although the trend is not nearly as strong as in the case of the F2 data (Table 2), this still may be indicative that the Avrami model is not a proper representation of the process under study.   Figure 9 presents the Arrhenius plots based on the heating function data ( Table 3). The lnχ(T) vs. T −1 plot yields the activation energy 141 kJ mol −1 , which is 65% larger than the value found in our previous work by using an advanced isoconversional method (Equations (16) and (17)). As established in the analysis of the simulated A3 and F2 data, more reasonable Arrhenius parameters are obtained while employing the lnχ(T)/n vs. T −1 plots. Indeed, the resulting E is 90 ± 3 kJ mol −1 and lnA is 20.9 ± 0.6. In agreement with the simulated results, the estimated E value is larger than the correct one by ~6%. In turn, the estimated lnA value is larger by ~12%.   Figure 9 presents the Arrhenius plots based on the heating function data ( Table 3). The lnχ(T) vs. T −1 plot yields the activation energy 141 kJ mol −1 , which is 65% larger than the value found in our previous work by using an advanced isoconversional method (Equations (16) and (17)). As established in the analysis of the simulated A3 and F2 data, more reasonable Arrhenius parameters are obtained while employing the lnχ(T)/n vs. T −1 plots. Indeed, the resulting E is 90 ± 3 kJ mol −1 and lnA is 20.9 ± 0.6. In agreement with the simulated results, the estimated E value is larger than the correct one by~6%. In turn, the estimated lnA value is larger by~12%. the value found in our previous work by using an advanced isoconversional method (Equations (16) and (17)). As established in the analysis of the simulated A3 and F2 data, more reasonable Arrhenius parameters are obtained while employing the lnχ(T)/n vs. T −1 plots. Indeed, the resulting E is 90 ± 3 kJ mol −1 and lnA is 20.9 ± 0.6. In agreement with the simulated results, the estimated E value is larger than the correct one by ~6%. In turn, the estimated lnA value is larger by ~12%.

Conclusions
Indubitably, the Ozawa-Avrami model can be applied to precisely determine the Avrami exponent as long as the process is established to obey the Avrami model. In the case of crosslinking polymerization, one cannot simply assume that the Avrami model applies. This is because the process can follow non-sigmoid (decelerating) reaction-order kinetics or more diverse autocatalytic kinetics. Only in the latter case the Avrami model may apply as a reasonable approximation. The most straightforward way to establish the Avrami model applicability is to perform at least one isothermal run and see whether its data can be fit by this model. An indirect way of testing this is to determine the value of conversion related to the rate (DSC) maximum, α p , under non-isothermal conditions. For the Avrami models, the respective value should be within the range 0.61-0.63 [55]. This criterion would help to differentiate the Avrami models from the second-order reaction model, for which α p is in the range 0.45-0.50. However, this criterion cannot distinguish the Avrami models from the first-order reaction model, for which α p is in the range 0.59-0.63. In addition, one should watch for a systematic dependence of the Avrami exponent on temperature as it can be a warning sign of the inapplicability of the Avrami model.
Concerning using the χ(T) function for determining the Arrhenius parameters, we have demonstrated that this function has a more complex temperature dependence than the regular rate constant, k(T). Yet, within certain limits, χ(T) can at best provide a mediocre approximation for k(T). However, if one is to use lnχ(T) in the Arrhenius plot, the values must be necessarily scaled by the values of the respective Avrami indexes (i.e., used as lnχ(T)/n). Without such a transformation, the errors in the activation energy can easily reach tens or even hundreds of percent. Additionally, not performing this transformation may produce breakpoints in the Arrhenius plots that confuse kinetic interpretations. Even if the transformation is performed, the Arrhenius plots of lnχ(T)/n still fail to produce the exact values of E and lnA.
Overall, the application of the Ozawa-Avrami model to non-isothermal polymerization is rather difficult to justify. It is limited to a single model, prone to computational artifacts, not very accurate, and not particularly easy to use. More importantly, there is a good number of time-proven kinetic techniques that are flexible, very accurate, and quite simple computationally [31]. Moreover, they have a solid track record of successful applications to crosslinking polymerization as well as many other processes.