Next Article in Journal
Phytochemical and Biological Studies of Nepeta asterotricha Rech. f. (Lamiaceae): Isolation of Nepetamoside
Next Article in Special Issue
Thermal Decomposition Kinetics and Mechanism of In-Situ Prepared Bio-Based Poly(propylene 2,5-furan dicarboxylate)/Graphene Nanocomposites
Previous Article in Journal
Well-Wrapped Li-Rich Layered Cathodes by Reduced Graphene Oxide towards High-Performance Li-Ion Batteries
Previous Article in Special Issue
Kinetic Analysis of Digestate Slow Pyrolysis with the Application of the Master-Plots Method and Independent Parallel Reactions Scheme
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

by
Nicolas Sbirrazzuoli
University Côte d’Azur, Institute of Chemistry of Nice, UMR CNRS 7272, 06100 Nice, France
Molecules 2019, 24(9), 1683; https://doi.org/10.3390/molecules24091683
Submission received: 9 April 2019 / Revised: 24 April 2019 / Accepted: 25 April 2019 / Published: 30 April 2019
(This article belongs to the Special Issue Thermal Analysis Kinetics for Understanding Materials Behavior)

Abstract

:
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.

Graphical Abstract

1. 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 (Tg). 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 in a previous experimental work [4]. This complex cure is of peculiar interest as it involves several steps consisting of an autocatalytic step, a first-order reaction and a diffusion-controlled part at the end of the reaction. This mechanism was frequently observed in the crosslinking of various epoxy–amine systems and was selected because it represents a good example of reaction complexity [5]. The first stages of the reaction are often controlled by an autocatalytic process, followed by epoxy–amine addition and a diffusion-controlled part after gelation when the viscosity of the reaction medium reaches high values as a result of the high increase of the molecular weight due to crosslinking. The specific polymerization system used as an example was made of 1,3-phenylenediamine (m-PDA) and diglycidyl ether of bisphenol A (DGEBA), obtained from Sigma-Aldrich and used as received. DGEBA has a molecular weight of about 355 g mol−1, a glass transition temperature Tg of about −20 °C (midpoint DSC) and an epoxy equivalent (determined by 1 H NMR) of about 175 g equiv−1. The example proposed here to illustrate the method corresponds to a polymerization reaction which is well described by an autocatalytic model for the chemically controlled part of the reaction and by a diffusion model for the end of the reaction. Nevertheless, this procedure can be extended to other complex processes, such as parallel, competitive (consecutive) or other mechanisms.

2. 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:
α t = t i t ( d Q / d t ) d t t i t f ( d Q / d t ) d t = Q t Q t o t
where ti represents the first integration bound of the DSC signal and tf is the last integration bound selected when the reaction is finished. (dQ/dt) is the heat flow measured by DSC at time t, Qtot is the total heat released (or absorbed) by the reaction and Qt is the current heat change.
The general form of the basic rate equation is usually written as [5]:
d α d t = A exp ( E R T ) f ( α )
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):
Φ ( E α ) = i = 1 n j i n J [ E α , T i ( t α ) ] J [ E α , T j ( t α ) ]
J [ E α , T ( t α ) ] = t α Δ α t α exp [ E α R T ( t ) ] d t
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 Ti(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 E0 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]:
ln ( d α d t ) α , i = ln [ A α f ( α ) ] E α R T α , i
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):
( d α d t ) α = [ A α f ( α ) ] exp ( E α R T α )
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.
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):
[ d ln ( d α / d t ) d T 1 ] α = [ d ln k ( T ) d T 1 ] α + [ d ln f ( α ) d T 1 ] α = E α R
In this equation, the derivative of the term containing f(α) is zero because each computation is performed for a constant value of α (isoconversional methods).

3. 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]:
k D ( T , α ) = D 0 exp ( E D R T + K α )
k C ( T , α ) = A 1 exp ( E 1 R T ) + A 2 exp ( E 2 R T ) α m   [ 19 , 20 ]
k e f = k C 1 + k D 1   [ 18 ]
d α d t = k e f ( 1 α ) n
α i + 1 = α i + t i t i + 1 ( d α d t ) i d t
and using the following parameters [4]: A1 = 20739.00 s−1, E1 = 70.0 kJ·mol−1, m = 1, A2 = 499.00 s−1, E2 = 45.0 kJ·mol−1, n = 1, D0 = 1.43 s−1, ED = 4.4 kJ·mol−1, K = −7.06. Here kD, kC and kef respectively represent the specific rate constant for diffusion, the rate constant for the chemically controlled reaction and the effective rate constant. D0 represents the pre-exponential factor of the diffusion-controlled reaction, A1 the pre-exponential factor for the non-catalyzed reaction and A2 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. ED, E1 and E2 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):
k C ( T ) = A 1 exp ( E 1 R T )
and using the following parameters: A1 = 1762.24 s−1, E1 = 50.0 kJ·mol−1, n = 1, D0 = 1.43 s−1, ED = 4.4 kJ·mol−1, K = −7.06.

4. Results

4.1. Autocatalytic Reaction with Diffusion-Controlled Part (Data Set 1)

4.1.1. Reaction Rate and Extent of Conversion for Non-Isothermal and Isothermal Conditions

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.

4.1.2. 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 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 (r2 = 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(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.

4.1.3. 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.
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.
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 Figure 4, Figure 5 and Figure 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].

4.1.4. Variation of the Rate Coefficients with Extent of Conversion

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).

4.1.5. 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.

4.1.6. 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]:
E α = k 1 ( T ) E 1 + k 2 ( T ) E 2 α m k 1 ( T ) + k 2 ( T ) α m
E α = ( A 1 / A 2 ) exp ( E 1 / R T ) E 1 + exp ( E 2 / R T ) E 2 α m ( A 1 / A 2 ) exp ( E 1 / R T ) + exp ( E 2 / R T ) α m
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]:
E α = k ( T ) E D + k D ( T , α ) E 2 k ( T ) + k D ( T , α )
E α = ( A 2 / D 0 ) exp ( E 2 / R T ) E D + exp ( E D / R T + K α ) E 2 ( A 2 / D 0 ) exp ( E 2 / R T ) + exp ( E D / R T + K α )
Note that, according to Equation (15), at the lowest extent of conversion (α → 0) Eα tends toward the activation energy of the uncatalyzed reaction (EαE1), 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 A1/A2 (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.
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 A1, A2 and D0. The results are given in Table 3.
E α = A 1 exp ( E 1 / R T ) E 1 + A 2 exp ( E 2 / R T ) E 2 α m A 1 exp ( E 1 / R T ) + A 2 exp ( E 2 / R T ) α m
E α = A 2 exp ( E 2 / R T ) E D + D 0 exp ( E D / R T + K α ) E 2 A 2 exp ( E 2 / R T ) + D 0 exp ( E D / R T + K α )
It can be seen that the parameters are in very good agreement with the reference values and the values of A1, A2 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 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].

4.2. First-Order Reaction with Diffusion-Controlled Part (Data Set 2)

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.
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, r2 = 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(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.

5. 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.

Funding

This research received no external funding.

Conflicts of Interest

The author declares no conflict of interest.

References

  1. Pascault, J.P.; Sautereau, H.; Verdu, J.; Williams, R.J.J. Thermosetting polymers; Marcel Dekker Inc.: New York, NY, USA, 2002. [Google Scholar]
  2. Hale, A.; Macosko, C.W.; Bair, H.E. Glass Transition Temperature as a Function of Conversion in Thermosetting Polymers. Macromolecules 1991, 24, 2610–2621. [Google Scholar] [CrossRef]
  3. Pascault, J.P.; Williams, R.J.J. Glass Transition Temperature Versus Conversion Relationships for Thermosetting Polymers. J. Polym. Sci. Part B Polym. Phys. 1990, 28, 85–95. [Google Scholar] [CrossRef]
  4. Alzina, C.; Sbirrazzuoli, N.; Mija, A. Hybrid Nanocomposites: Advanced Nonlinear Method for Calculating Key Kinetic Parameters of Complex Cure Kinetics. J. Phys. Chem. B 2010, 114, 12480–12487. [Google Scholar] [CrossRef] [PubMed]
  5. Vyazovkin, S.; Sbirrazzuoli, N. Isoconversional kinetic analysis of thermally stimulated processes in polymers. Macromol. Rapid Comm. 2006, 27, 1515–1532. [Google Scholar] [CrossRef]
  6. Vyazovkin, S.; Burnham, A.K.; Criado, J.M.; Pérez-Maqueda, L.A.; Popescu, C.; Sbirrazzuoli, N. ICTAC kinetics committee recommendations for performing kinetic computations on thermal analysis data. Thermochim. Acta 2011, 520, 1–19. [Google Scholar] [CrossRef]
  7. Vyazovkin, S. Isoconversional Kinetics of Thermally Stimulated Processes; Springer: Berlin, Germany, 2015. [Google Scholar]
  8. Vyazovkin, S. Evaluation of activation energy of thermally stimulated solid-state reactions under arbitrary variation of temperature. J. Comput. Chem. 1997, 18, 393–402. [Google Scholar] [CrossRef]
  9. Sbirrazzuoli, N.; Vincent, L.; Vyazovkin, S. Comparison of several computational procedures for evaluating the kinetics of thermally stimulated condensed phase reactions. Chemometr. Intell. Lab 2000, 54, 53–60. [Google Scholar] [CrossRef]
  10. Vyazovkin, S. Modification of the integral isoconversional method to account for variation in the activation energy. J. Comput. Chem. 2001, 22, 178–183. [Google Scholar] [CrossRef]
  11. Sbirrazzuoli, N. Determination of pre-exponential factors and of the mathematical functions f(α) or G(α) that describe the reaction mechanism in a model-free way. Thermochim. Acta 2013, 564, 59–69. [Google Scholar] [CrossRef]
  12. Sbirrazzuoli, N. Is the Friedman method applicable to transformations with temperature dependent reaction heat? Macromol. Chem. Phys. 2007, 208, 1592–1597. [Google Scholar] [CrossRef]
  13. Sbirrazzuoli, N.; Vincent, L.; Vyazovkin, S. Electronic solution to the problem of a kinetic standard for DSC measurements. Chemometr. Intell. Lab 2000, 52, 23–32. [Google Scholar] [CrossRef]
  14. Sbirrazzuoli, N.; Brunel, D.; Elegant, L. Different kinetic equations analysis. J. Therm. Anal. 1992, 38, 1509–1524. [Google Scholar] [CrossRef]
  15. Sbirrazzuoli, N.; Girault, Y.; Elegant, L. Simulations for evaluation of kinetic methods in differential scanning calorimetry. Part 3—peak maximum evolution methods and isoconversional methods. Thermochim. Acta 1997, 293, 25–37. [Google Scholar] [CrossRef]
  16. Falco, G.; Guigo, N.; Vincent, L.; Sbirrazzuoli, N. FA polymerization disruption by protic polar solvent. Polymers 2018, 10, 529. [Google Scholar] [CrossRef] [PubMed]
  17. Friedman, H.L. Kinetics of thermal degradation of char-forming plastics from thermogravimetry. Application to a phenolic plastic. J. Polym. Sci. Part C 1964, 6, 183–195. [Google Scholar] [CrossRef]
  18. Vyazovkin, S.; Sbirrazzuoli, N. Mechanism and kinetics of epoxy-amine cure studied by differential scanning calorimetry. Macromolecules 1996, 29, 1867–1873. [Google Scholar] [CrossRef]
  19. Ryan, M.E.; Dutta, A. Kinetics of epoxy cure: A rapid technique for kinetic parameter estimation. Polymer 1979, 20, 203–206. [Google Scholar] [CrossRef]
  20. Sourour, S.; Kamal, M.R. Differential scanning calorimetry of epoxy cure: Isothermal cure kinetics. Thermochim. Acta 1976, 14, 41–59. [Google Scholar] [CrossRef]
  21. Yamini, G.; Shakeri, A.; Vafayan, M.; Zohuriaan-Mehr, M.J.; Kabiri, K.; Zolghadr, M. Cure kinetics of modified lignosulfonate/epoxy blends. Thermochim. Acta 2019, 675, 18–28. [Google Scholar] [CrossRef]
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.
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.
Molecules 24 01683 g001
Figure 2. Non-isothermal data. Variation of the reaction rate (dα/dt) and the extent of conversion (α) with temperature (T) for data set 1. The heating rate of each experiment (in K min−1) is indicated by each curve.
Figure 2. Non-isothermal data. Variation of the reaction rate (dα/dt) and the extent of conversion (α) with temperature (T) for data set 1. The heating rate of each experiment (in K min−1) is indicated by each curve.
Molecules 24 01683 g002
Figure 3. Isothermal data. Variation of the reaction rate (dα/dt) with time (t) for data set 1. Inset: Variation of the extent of conversion (α) with time (t). The temperature of each experiment (in °C) is indicated by each curve.
Figure 3. Isothermal data. Variation of the reaction rate (dα/dt) with time (t) for data set 1. Inset: Variation of the extent of conversion (α) with time (t). The temperature of each experiment (in °C) is indicated by each curve.
Molecules 24 01683 g003
Figure 4. 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, red solid lozenges: Eα computed with FR method, green circles: Eα computed with NLN method. Inset: isothermal data. Solid triangles: Eα computed with FR method, solid stars: Eα computed with NLN method.
Figure 4. 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, red solid lozenges: Eα computed with FR method, green circles: Eα computed with NLN method. Inset: isothermal data. Solid triangles: Eα computed with FR method, solid stars: Eα computed with NLN method.
Molecules 24 01683 g004
Figure 5. Non-isothermal data. Dependence of the effective activation energy (Eα) with temperature. Red solid lozenges: Eα computed with FR method, green open circles: Eα computed with NLN method.
Figure 5. Non-isothermal data. Dependence of the effective activation energy (Eα) with temperature. Red solid lozenges: Eα computed with FR method, green open circles: Eα computed with NLN method.
Molecules 24 01683 g005
Figure 6. Dependence of the effective activation energy (Eα) with the extent of conversion (α). Red lozenges: Eα computed with FR method, green circles: Eα computed with NLN method, blue triangles: (dα/dt) computed as the product of [Aα f(α)] by exp[−Eα/(RTα)] (Equation (6)).
Figure 6. Dependence of the effective activation energy (Eα) with the extent of conversion (α). Red lozenges: Eα computed with FR method, green circles: Eα computed with NLN method, blue triangles: (dα/dt) computed as the product of [Aα f(α)] by exp[−Eα/(RTα)] (Equation (6)).
Molecules 24 01683 g006
Figure 7. Red squares: variation of ln [Aα f(α)] with α, green circles: variation of −Eα/(RTα) with α, blue triangles: variation of (dα/dt) computed as the product of [Aα f(α)] by exp[−Eα/(RTα)] with α.
Figure 7. Red squares: variation of ln [Aα f(α)] with α, green circles: variation of −Eα/(RTα) with α, blue triangles: variation of (dα/dt) computed as the product of [Aα f(α)] by exp[−Eα/(RTα)] with α.
Molecules 24 01683 g007
Figure 8. Blue triangles: variation of (dα/dt) computed as the product of [Aα f(α)] by exp[−Eα/(RTα)] with temperature, red squares: variation of kD with temperature, green lozenges: variation of kC with temperature, magenta circles: variation of kef with temperature.
Figure 8. Blue triangles: variation of (dα/dt) computed as the product of [Aα f(α)] by exp[−Eα/(RTα)] with temperature, red squares: variation of kD with temperature, green lozenges: variation of kC with temperature, magenta circles: variation of kef with temperature.
Molecules 24 01683 g008
Figure 9. Green lozenges: variation of ln k(T) computed as ln k(T) = ln(Aα)–Eα/(RTα) as a function of reciprocal temperature, blue circles: variation of ln kef (T) computed using Equation (10) as a function of reciprocal temperature.
Figure 9. Green lozenges: variation of ln k(T) computed as ln k(T) = ln(Aα)–Eα/(RTα) as a function of reciprocal temperature, blue circles: variation of ln kef (T) computed using Equation (10) as a function of reciprocal temperature.
Molecules 24 01683 g009
Figure 10. Circles: Eα dependency obtained with NLN method, line: fit using the autocatalytic model for 0.02 < α < 0.46 and the diffusion model for 0.48 < α < 0.98 with the parameters of Table 2. Inset: same plot for the data of Table 3.
Figure 10. Circles: Eα dependency obtained with NLN method, line: fit using the autocatalytic model for 0.02 < α < 0.46 and the diffusion model for 0.48 < α < 0.98 with the parameters of Table 2. Inset: same plot for the data of Table 3.
Molecules 24 01683 g010
Figure 11. Non-isothermal data. Variation of the reaction rate (dα/dt) and extent of conversion (α) with temperature (T) for data set 2. The heating rate of each experiment (in K min−1) is indicated by each curve.
Figure 11. Non-isothermal data. Variation of the reaction rate (dα/dt) and extent of conversion (α) with temperature (T) for data set 2. The heating rate of each experiment (in K min−1) is indicated by each curve.
Molecules 24 01683 g011
Figure 12. Isothermal data. Variation of the reaction rate (dα/dt) with time (t) for data set 2. Inset: Variation of the extent of conversion (α) with time (t). The temperature of each experiment (in °C) is indicated by each curve.
Figure 12. Isothermal data. Variation of the reaction rate (dα/dt) with time (t) for data set 2. Inset: Variation of the extent of conversion (α) with time (t). The temperature of each experiment (in °C) is indicated by each curve.
Molecules 24 01683 g012
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.
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.
Molecules 24 01683 g013
Figure 14. Non-isothermal data. Dependence of the effective activation energy (Eα) with temperature (T). Red lozenges: Eα computed with FR method, green circles: Eα computed with NLN method.
Figure 14. Non-isothermal data. Dependence of the effective activation energy (Eα) with temperature (T). Red lozenges: Eα computed with FR method, green circles: Eα computed with NLN method.
Molecules 24 01683 g014
Table 1. Some values of the various terms used to identify a change in the rate-limiting step.
Table 1. Some values of the various terms used to identify a change in the rate-limiting step.
Tα/°CEα (FR)/kJ·mol−1αln (Aα/s−1)[Aα f(α)]/s−1exp(−Eα/(RTα)(dα/dt)/s−1
181.5645.050.425.47130.256.632 × 10−68.638 × 10−4
182.3444.390.445.30108.308.076 × 10−68.746 × 10−4
183.1143.630.465.1187.431.007 × 10−58.802 × 10−4
183.8942.760.484.8968.301.289 × 10−58.801 × 10−4
184.6641.780.504.6351.451.699 × 10−58.740 × 10−4
185.4540.680.524.3537.262.313 × 10−58.618 × 10−4
Table 2. Parameters obtained by fitting Equations (15) and (17).
Table 2. Parameters obtained by fitting Equations (15) and (17).
2 < α < 46%A1/A2E1/kJ·mol−1E2/kJ·mol−1mMSSDa
Autocatalytic5918.0379.838.20.970.4435
48 < α < 98%A/D0E2/kJ·mol−1ED/kJ·mol−1KMSSDa
Diffusion362.6748.94.9−7.840.0023
a Mean of the sum of squared deviations M S S D = ( 1 / n ) i = 1 n ( E c a l c E r e f ) 2 / E r e f .
Table 3. Parameters obtained by fitting Equations (18) and (19).
Table 3. Parameters obtained by fitting Equations (18) and (19).
A1/s−1A2/s−1E1/kJ·mol−1E2/kJ·mol−1mMSSDa
2 < α < 46%20756.17498.5267.642.61.20.6011
2 < α < 24%20258.31510.8476.146.71.30.0345
A2/s−1D0/s−1E2/kJ·mol−1ED/kJ·mol−1KMSSDa
48 < α < 98%498.91.4348.84.9−7.870.0024
a Mean of the sum of squared deviations M S S D = ( 1 / n ) i = 1 n ( E c a l c E r e f ) 2 / E r e f .

Share and Cite

MDPI and ACS Style

Sbirrazzuoli, N. Advanced Isoconversional Kinetic Analysis for the Elucidation of Complex Reaction Mechanisms: A New Method for the Identification of Rate-Limiting Steps. Molecules 2019, 24, 1683. https://doi.org/10.3390/molecules24091683

AMA Style

Sbirrazzuoli N. Advanced Isoconversional Kinetic Analysis for the Elucidation of Complex Reaction Mechanisms: A New Method for the Identification of Rate-Limiting Steps. Molecules. 2019; 24(9):1683. https://doi.org/10.3390/molecules24091683

Chicago/Turabian Style

Sbirrazzuoli, Nicolas. 2019. "Advanced Isoconversional Kinetic Analysis for the Elucidation of Complex Reaction Mechanisms: A New Method for the Identification of Rate-Limiting Steps" Molecules 24, no. 9: 1683. https://doi.org/10.3390/molecules24091683

APA Style

Sbirrazzuoli, N. (2019). Advanced Isoconversional Kinetic Analysis for the Elucidation of Complex Reaction Mechanisms: A New Method for the Identification of Rate-Limiting Steps. Molecules, 24(9), 1683. https://doi.org/10.3390/molecules24091683

Article Metrics

Back to TopTop