Kinetic Analysis of the Curing of a Partially Biobased Epoxy Resin Using Dynamic Differential Scanning Calorimetry

This research presents a cure kinetics study of an epoxy system consisting of a partially bio-sourced resin based on diglycidyl ether of bisphenol A (DGEBA) with amine hardener and a biobased reactive diluent from plants representing 31 wt %. The kinetic study has been carried out using differential scanning calorimetry (DSC) under non-isothermal conditions at different heating rates. Integral and derivative isoconversional methods or model free kinetics (MFK) have been applied to the experimental data in order to evaluate the apparent activation energy, Ea, followed by the application of the appropriate reaction model. The bio-sourced system showed activation energy that is independent of the extent of conversion, with Ea values between 57 and 62 kJ·mol−1, corresponding to typical activation energies of conventional epoxy resins. The reaction model was studied by comparing the calculated y(α) and z(α) functions with standard master plot curves. A two-parameter autocatalytic kinetic model of Šesták–Berggren [SB(m,n)] was assessed as the most suitable reaction model to describe the curing kinetics of the epoxy resins studied since it showed an excellent agreement with the experimental data.


Introduction
Epoxy resins are low molecular weight pre-polymers with at least one epoxide group in their structure. These epoxy rings are usually located in terminal positions due to the high reactivity of this position [1]. Currently epoxy resins are the most used thermosetting resins due to the huge range of properties they can cover such as tensile strength, high adhesion strength, low shrinkage, good chemical resistance, and low volatile emission, among others [2][3][4][5]. These outstanding properties make these resins suitable for high technological sectors such as the automotive and aerospace industries. In addition, their versatility makes them useful for industrial applications such as adhesives, paints, surface coatings, advanced composite materials, electrical and electronic components, and high-performance membranes for separation and filtration, among others [3,[6][7][8][9][10][11][12][13][14]. The final properties of a structural epoxy system are greatly influenced by several factors, which include the chemical structure of both the epoxy resin and the curing agent and external factors related to the curing procedure.
Despite the huge range of epoxy resins available worldwide, undoubtedly the most used resin is that based on diglycidyl ether of bisphenol A (DGEBA), which is obtained by reacting = K(T) f (α) (1) where K(T) is the rate coefficient, which depends on temperature T, and f (α) represents the reaction model, which takes different forms depending on the mechanism of the process. The temperature dependence of K(T) could be parameterized through the Arrhenius expression shown in Equation (2): where A is the pre-exponential factor, E a is the apparent activation energy, R refers to the universal gas constant, and T is the absolute temperature. By substituting the Arrhenius temperature coefficient, Equation (1) can be written as Equation (3), which represents the general rate equation for the kinetic study. As one can realize, this is a time-dependent equation, and subsequently, it can be applied to any curing cycle.
This time-dependent expression can be converted into a temperature domain by considering Equation (4), which shows the heating rate β definition: By substituting Equation (4) into Equation (3), we obtain the general equation for kinetic analysis for dynamic heating cycles: Isoconversional Methods: Model-Free Kinetics (MFK) Isoconversional methods are derived from the basic isoconversional principle that assumes that the reaction rate at a constant conversion (α c ) depends only on temperature. If we take natural logarithms on both sides of Equation (1), we obtain: ln dα dt = ln K(T) + ln f (α) (6) Then, taking partial derivatives with respect to the inverse temperature (1/T) on Equation (6) at a constant α = α i , Equation (7) is derived: Isoconversional conditions, meaning that it has a constant α value of α i such that its corresponding f (α) value is also constant. Hence, the second term on the right-hand side of Equation (7) is zero. The first term on the right-hand side of Equation (7) can be obtained by taking the partial derivative of the logarithmic Arrhenius coefficient in Equation (2) with respect to 1/T, as shown in Equation (8): By combining Equations (7) and (8), the isoconversional methods indicate that it is possible to obtain E a (at a particular α i = constant) without assuming any reaction model, as shown in Equation (9). For this reason, sometimes these methods are called model-free kinetic methods (MFK).
Isoconversional methods can give accurate values of the apparent activation energy E a . These methods can be classified as differential, such as Friedman's method, or integral methods, such as the Flynn-Wall-Ozawa (FWO), Kissinger-Akahira-Sunose (KAS), and Starink methods. The Friedman method [44] is one of the most common differential isoconversional methods used to evaluate the apparent activation energy. By assuming that the reaction model f (α) does not change with the conversion α, the basic differential expression (Equation (3)) can be written in a natural logarithm form, as shown in Equation (10): As suggested by Equation (10), by plotting ln β dα dT α,i versus 1 , the E a can be calculated through the slope of the linear fitting. This method must be evaluated for each value of the conversion α.
As suggested by Sergey Vyazovkin et al. [43], to check its possible variation with the conversion α, it is recommended to start from α i = 0.05 up to α i = 0.95 with an increasing step size of 0.05. The Friedman method does not make any other assumption, and in a first approach, it gives a quite accurate value of E a . Isoconversional integral methods are based on integration of Equation (3) or Equation (5), as observed in Equation (11): The integral of the temperature-dependent expression from the Arrhenius constant is the so-called "temperature integral" and does not have a direct analytical solution. For this reason, different approximations have been proposed and used, leading to different methods, and obviously, to different accuracies on the E a estimation. The Flynn-Wall-Ozawa (FWO) method shows very low accuracy as it uses a rude approximation of the temperature integral. Despite this, the FWO is one of the most commonly used. It follows Equation (12) [45]. This suggests that a plot of ln(β i ) versus the inverse temperature (1/T α ) for a constant conversion α gives a straight line whose slope is −1.052/R.
The Kissinger-Akahira-Sunose [46] method gives more accurate E a values as it uses a more accurate solution of the temperature integral. The basic expression for the KAS method is shown in Equation (13), which suggests a linear correlation between ln β i T 2 α,i and the inverse temperature (1/T α ) for a particular conversion α value.
A better evaluation of the activation energy (E a ) can be done using the Starink expression, which is similar to that of the KAS method, but the particular parameters are optimized, as observed in Equation (14). In a similar way to the above-mentioned methods, a plot of ln β i T 1.92 α,i against (1/T α ) for a particular conversion α value gives a linear correlation whose slope is −1.0008/R [47].
Although it is not considered an isoconversional method, the Kissinger method (not to be confused with the Kissinger-Akahira-Sunose method) is based on the peak temperature (T p ) observed for the conversion rate at different heating rates β i [48]. This method is widely used because of its simplicity.
It assumes that d dt dα dt = 0 at the temperature for the maximum conversion rate (T p ). The assumptions of the Kissinger method are only applicable if the conversion for this peak temperature α m is almost constant for a series of heating rates, β i . The Kissinger method follows Equation (15). One important drawback of the Kissinger method is that it only provides a single E a value (at α m ) and cannot follow the evolution of E a with the conversion α. Despite this, it is widely used due to its simplicity.
By plotting ln β T 2 p versus 1 T p , the activation energy can be obtained through the slope of the linear fit. One of the restrictions of the Kissinger method is that it assumes a single-step process.
To check the existence of a single-step degradation process, different methods have been proposed. In particular, it is worthy to note the methodology described by Farjas et al. [49]. This method is not only based on the peak temperature, but also on the peak width measured as the time corresponding to the full width at half maximum (∆t FWHM ), which is more sensitive to the existence of multiple processes. If a plot of the ∆t FWHM against 1/T p gives a linear correlation, and the obtained E a by linear fitting is similar to that obtained by the Kissinger method, then it is possible to conclude that the curing process takes place in a single step (despite the complexity of the occurring reactions).
Isoconversional methods are very useful for giving an accurate estimation of the E a and do not assume any reaction model f (α). Polymerization processes are very complex, as indicated by D.R. D'hooge et al. [50]. In particular they use multi-scale modelling to unveil some of these complex processes governed by free radical polymerization (FRP) or by controlled radical polymerization (CRP). They also pay special attention to chain mobility, which is restricted by crosslinking. To obtain the kinetic triplet, it is necessary to determine the reaction model. One of the most suitable methods is the use of the Málek method [51], which allows determining the reaction model by experimentally calculating two functions defined as y(α) and z(α), as indicated in Equations (16) and (17), with x = E a /RT, i.e., the reduced temperature: The particularity of these two functions is the following. Regarding y(α), if we substitute Equation (3) into Equation (16), then Equation (18) is obtained. The relevance of this function is that it can be calculated numerically using Equation (16), and as suggested by Equation (18), it provides a plot representation proportional to f (α) by the factor A. Therefore, y(α) gives a clear idea of the reaction model. y(α) = A f (α) (18) Regarding the z(α) function, the solution to the integral expression (Equation (11)) can be defined as: where P(x) represents a particular approximation of the temperature integral. The relation between P(x) and π(x) is shown in Equation (20): By combining Equations (3), (19), and (20), Equation (17) results in the following expression: which suggests that the calculated z(α) using Equation (17), can be compared with different theoretical z(α) curves with well known f (α) and g(α) functions as shown in Table 1. Table 1. Summary of some f (α) and g(α) functions corresponding to crosslinking of polymers.

Reaction Model
Random nucleation with two nuclei of individual particle Senum and Yang have proposed different approximations of π(x) in a rational form [52] and Perez-Maqueda et al. [53] have indicated that the use of a rational approximation of π(x) with a fourth degree polynomial is enough to give good accuracy. Flynn [54] suggested a correction for the fourth rational approximation expression, as we can see at the Equation (14). These functions y(α) and z(α) can be calculated from experimental data.

Materials
The commercial epoxy resin used was a Resoltech ® 1070 ECO and the hardener was an amine-based type Resoltech ® 1074 ECO, both of them supplied by Castro Composites (Pontevedra, Spain). This partially bio-sourced epoxy resin is based on a mixture of a diglycidyl ether of bisphenol A (DGEBA) and a plant-based epoxy reactive diluent. As indicated by the supplier, the cured resin contains 31% biobased content (according to ASTM D6866-12) and maintains good transparency. The eco-epoxy resin and the hardener were mixed under the stoichiometric ratio 100:35 (%wt epoxy:%wt hardener), following the manufacturer's recommendations.
Scheme 1 shows the chemical structures of the components of the epoxy resin as well as the hardener.
The main component on the epoxy resin is the reaction product of bisphenol A and epichlorohydrin, leading to a typical diglycidyl ether of bisphenol A (DGEBA) resin. With regard to the reactive diluent, it is mainly derived from plant or vegetable oils, in particular, epoxidized vegetable oils. Finally, the hardener is an amine type, mainly composed of 3-aminomethyl-3,5,5trimethylcyclohexylamine and the reaction product of propoxylated propylidynetrimethanol and ammonia. The main component on the epoxy resin is the reaction product of bisphenol A and epichlorohydrin, leading to a typical diglycidyl ether of bisphenol A (DGEBA) resin. With regard to the reactive diluent, it is mainly derived from plant or vegetable oils, in particular, epoxidized vegetable oils. Finally, the hardener is an amine type, mainly composed of 3-aminomethyl-3,5,5trimethylcyclohexylamine and the reaction product of propoxylated propylidynetrimethanol and ammonia. A preliminary FTIR characterization was carried out. The spectra were recorded using a spectrometer Spectrum BX from Perkin-Elmer (Madrid, Spain) together with a special holder for liquid samples. Data were collected from 20 scans between 5500 and 450 cm −1 at a spectral resolution of 4 cm −1 . Figure 1 shows the FTIR spectra of both the epoxy resin ( Figure 1a) and the hardener (Figure 1b). The FTIR spectrum of the epoxy base epoxy resin shows the typical peaks and bands of a DGEBPA resin. These peak/band assignments are as follows: O-H stretching at about 3500 cm −1 ; stretching of the C-H of the oxirane ring at 3057 cm −1 ; stretching of C-H of CH 2 and CH, both aromatic and aliphatic, at 2965-2873 cm −1 ; stretching of C=C of aromatic rings at 1608 cm −1 ; stretching of the C-C of aromatic rings at 1509 cm −1 ; stretching of the C-O-C of ethers at 1036 cm −1 ; asymmetric stretching of the C-O in oxirane rings at 1241 cm −1 ; stretching of the C-O of oxirane ring at 915 cm −1 ; stretching of the C-O-C of oxirane ring at 831 cm −1 ; and rocking of CH 2 at 772 cm −1 .
With regard to the reactive diluent, an epoxidized vegetable oil contained overlapping peaks and bands of CH 2 , CH, and oxirane rings, while the ester groups typical of the triglyceride appeared with a weak signal as it was not the main component. These ester peaks were characterized by a weak signal at 645 cm −1 attributable to the O-C=O bending; a peak at 1180 cm −1 , which overlaps with ester groups and corresponds to the antisymmetric stretch of C-O-C; and a weak signal at 1765-1720 cm −1 , which is characteristic of the C=O stretch. These weak signals typical of the ester groups indicated the presence of triglycerides as reactive diluents.
Regarding the hardener, it was an amine type ( Figure 1b). The peak/band assignment was the following: NH stretching of -NH 2 at 3370-3280 cm −1 , symmetric and antisymmetric stretching of CH in -CH 2 at 2950-2890 cm −1 , symmetric and antisymmetric stretching of CH in -CH 2 -, -NH 2 deformation at 1590 cm −1 , scissors vibration of -CH 2 -at 1370 cm −1 , symmetric deformation of -CH 3 at 1370 cm −1 , stretching of C-O-C ether groups at 1250 cm −1 , stretching of -NH in C-NH 2 primary amines located at 1100 cm −1 , and wagging of -NH in R-NH 2 primary amines. A preliminary FTIR characterization was carried out. The spectra were recorded using a spectrometer Spectrum BX from Perkin-Elmer (Madrid, Spain) together with a special holder for liquid samples. Data were collected from 20 scans between 5500 and 450 cm −1 at a spectral resolution of 4 cm −1 . Figure 1 shows the FTIR spectra of both the epoxy resin ( Figure 1a) and the hardener ( Figure  1b). The FTIR spectrum of the epoxy base epoxy resin shows the typical peaks and bands of a DGEBPA resin. These peak/band assignments are as follows: O-H stretching at about 3500 cm −1 ; stretching of the C-H of the oxirane ring at 3057 cm −1 ; stretching of C-H of CH2 and CH, both aromatic and aliphatic, at 2965-2873 cm −1 ; stretching of C=C of aromatic rings at 1608 cm −1 ; stretching of the C-C of aromatic rings at 1509 cm −1 ; stretching of the C-O-C of ethers at 1036 cm −1 ; asymmetric stretching of the C-O in oxirane rings at 1241 cm −1 ; stretching of the C-O of oxirane ring at 915 cm −1 ; stretching of the C-O-C of oxirane ring at 831 cm −1 ; and rocking of CH2 at 772 cm −1 .
With regard to the reactive diluent, an epoxidized vegetable oil contained overlapping peaks and bands of CH2, CH, and oxirane rings, while the ester groups typical of the triglyceride appeared with a weak signal as it was not the main component. These ester peaks were characterized by a weak signal at 645 cm -1 attributable to the O-C=O bending; a peak at 1180 cm −1 , which overlaps with ester groups and corresponds to the antisymmetric stretch of C-O-C; and a weak signal at 1765-1720 cm −1 , which is characteristic of the C=O stretch. These weak signals typical of the ester groups indicated the presence of triglycerides as reactive diluents.
(a) Complementary to this preliminary characterization of the epoxy system, some additional properties are described. Some physical properties of the base epoxy, the hardener, and the uncured mixture are summarized in Table 2. With regard to the reactivity of this system, it is worthy to note that the exothermic peak of the epoxy system during curing at 23 °C had a value of 185 °C and occurred at 31 min. It showed a gel time of 28 min when the mixture did not exceed a volume of 70 mL, as indicated by the supplier.
In addition to these properties, the glass transition temperature, Tg, of the cured resin was obtained via dynamic mechanical thermal analysis (DMTA). The curing cycle was 1 h at 90 °C and had a post-curing stage of 30 min at 150 °C. DMTA was carried out in an oscillatory rheometer AR-G2 supplied from TA Instruments (New Castle, DE, USA), using a special clamp system made for solid samples (4 × 10 × 40 mm 3 ) working in torsion-shear conditions. Cured samples were subjected to a temperature sweep from 30 °C to 140 °C at a heating rate of 2 K min −1 , with a maximum deformation () of 0.1%, and a frequency of 1 Hz. Figure 2 shows the plot evolution of the storage modulus (G') and the damping factor (tan ). By taking the Tg as the peak maximum of tan , a value of 95.3 °C was obtained for this cured epoxy system. Complementary to this preliminary characterization of the epoxy system, some additional properties are described. Some physical properties of the base epoxy, the hardener, and the uncured mixture are summarized in Table 2. With regard to the reactivity of this system, it is worthy to note that the exothermic peak of the epoxy system during curing at 23 • C had a value of 185 • C and occurred at 31 min. It showed a gel time of 28 min when the mixture did not exceed a volume of 70 mL, as indicated by the supplier.
In addition to these properties, the glass transition temperature, T g , of the cured resin was obtained via dynamic mechanical thermal analysis (DMTA). The curing cycle was 1 h at 90 • C and had a post-curing stage of 30 min at 150 • C. DMTA was carried out in an oscillatory rheometer AR-G2 supplied from TA Instruments (New Castle, DE, USA), using a special clamp system made for solid samples (4 × 10 × 40 mm 3 ) working in torsion-shear conditions. Cured samples were subjected to a temperature sweep from 30 • C to 140 • C at a heating rate of 2 K min −1 , with a maximum deformation (γ) of 0.1%, and a frequency of 1 Hz. Figure 2 shows the plot evolution of the storage modulus (G ) and the damping factor (tan δ). By taking the T g as the peak maximum of tan δ, a value of 95.3 • C was obtained for this cured epoxy system. Finally, the Shore D hardness of the cured epoxy material was obtained in a 676-D durometer from J. Bot Instruments (Barcelona, Spain) at room temperature following ISO 868 standard, resulting in a Shore D value of 83.4 ± 1.7.

Differential Scanning Calorimetry (DSC) Kinetic Measurements
The curing kinetics of the eco-epoxy system was analyzed using differential scanning calorimetry (DSC) in a 821 DSC calorimeter from Mettler-Toledo, Inc. (Schwerzenbach, Switzerland). The sample weight was between 10-13 mg, following the recommendations of Vyazovkin et al. [55], suggesting that the sampleʹs mass should be an inverse proportion to the heating rate. The mixtures were prepared at room temperature using the stoichiometric ratio, and then the liquid mixture was stirred until homogenization. Subsequently, the corresponding mass was placed into a standard aluminium pan (40 L) and sealed with a press. Finally, the samples were subjected to the following thermal program: first, an isothermal stage at 20 °C for 1 min was programmed to stabilize the sample temperature, then a heating process up to (250 °C-280 °C) at different heating rates,  (2.5, 5, 10, and 20 K min −1 . This was followed by a cooling step down to 0 °C at a cooling rate of −20 K min −1 , and finally, a second heating step up to 180 °C at 10 K min -1 was scheduled. All tests were carried out under a dry atmosphere with a constant nitrogen flow of 30 mL min −1 . DSC was used to obtain the heat flow as a function of temperature (from the first heating), and the glass temperature (Tg) (from the second heating).
Processing the thermal data from DSC is necessary to obtain the extent of conversion . It could be expressed as indicated in Equation (23), where ∆H t is the released enthalpy or heat at a particular time t, and ∆H T is the total enthalpy released during the complete reaction. Equation (24) shows how the conversion rate can be obtained from the heat flow of the DSC thermograms. Finally, the Shore D hardness of the cured epoxy material was obtained in a 676-D durometer from J. Bot Instruments (Barcelona, Spain) at room temperature following ISO 868 standard, resulting in a Shore D value of 83.4 ± 1.7.

Differential Scanning Calorimetry (DSC) Kinetic Measurements
The curing kinetics of the eco-epoxy system was analyzed using differential scanning calorimetry (DSC) in a 821 DSC calorimeter from Mettler-Toledo, Inc. (Schwerzenbach, Switzerland). The sample weight was between 10-13 mg, following the recommendations of Vyazovkin et al. [55], suggesting that the sample's mass should be an inverse proportion to the heating rate. The mixtures were prepared at room temperature using the stoichiometric ratio, and then the liquid mixture was stirred until homogenization. Subsequently, the corresponding mass was placed into a standard aluminium pan (40 µL) and sealed with a press. Finally, the samples were subjected to the following thermal program: first, an isothermal stage at 20 • C for 1 min was programmed to stabilize the sample temperature, then a heating process up to (250 • C-280 • C) at different heating rates, β (2.5, 5, 10, and 20 K min −1 . This was followed by a cooling step down to 0 • C at a cooling rate of −20 K min −1 , and finally, a second heating step up to 180 • C at 10 K min −1 was scheduled. All tests were carried out under a dry atmosphere with a constant nitrogen flow of 30 mL min −1 . DSC was used to obtain the heat flow as a function of temperature (from the first heating), and the glass temperature (T g ) (from the second heating).
Processing the thermal data from DSC is necessary to obtain the extent of conversion α. It could be expressed as indicated in Equation (23), where ∆H t is the released enthalpy or heat at a particular time t, and ∆H T is the total enthalpy released during the complete reaction. Equation (24) shows how the conversion rate can be obtained from the heat flow of the DSC thermograms.

Estimation of the Apparent Activation Energy, E a
The curing reaction of the partially bio-based epoxy resin at different heating rates was studied using DSC. The first heating cycle was useful for determining the total enthalpy (∆H T ) and the temperature for the maximum curing rate (T p ), while the second heating cycle was used to study the glass transition temperature (T g ) of the cured resin. These values are gathered in Table 3 for different heating rates between 2.5 and 20 K min −1 . It is important to remark that the T g decreased with increasing heating rate from 95.8 • C (2.5 K min −1 ) down to 91.1 • C (for 20 K min −1 ) which was also in accordance with the curing enthalpy (∆H T ). Low curing rates allow diffusion and crosslinking reactions to occur more readily. At high curing rates, the overall speed is high and does not allow free diffusion as a crosslinked structure is immediately formed [56].  Figure 3 shows the conversion α as a function of the absolute temperature T. As the heating rate β increased, the characteristic sigmoidal curve moved toward higher values.

Estimation of the Apparent Activation Energy, Ea.
The curing reaction of the partially bio-based epoxy resin at different heating rates was studied using DSC. The first heating cycle was useful for determining the total enthalpy (∆H T ) and the temperature for the maximum curing rate (Tp), while the second heating cycle was used to study the glass transition temperature (Tg) of the cured resin. These values are gathered in Table 3 for different heating rates between 2.5 and 20 K min −1 . It is important to remark that the Tg decreased with increasing heating rate from 95.8 °C (2.5 K min -1 ) down to 91.1 °C (for 20 K min -1 ) which was also in accordance with the curing enthalpy (∆H T ). Low curing rates allow diffusion and crosslinking reactions to occur more readily. At high curing rates, the overall speed is high and does not allow free diffusion as a crosslinked structure is immediately formed [56].  Figure 3 shows the conversion α as a function of the absolute temperature T. As the heating rate β increased, the characteristic sigmoidal curve moved toward higher values.   Figure 4 shows a plot of the curing rate ( dα dT ) as a function of the conversion α obtained using Equation (24). The geometry of these plots is almost the same, thus suggesting there is no change in the reaction model with the heating rate. Equation (24). The geometry of these plots is almost the same, thus suggesting there is no change in the reaction model with the heating rate. One of the most important parameters of a kinetic study is the apparent activation energy, Ea. As indicated previously, the Kissinger method allows for estimating a single value of Ea for the curing process by using Equation (15). One important condition needed to apply the Kissinger method is that the conversion at the maximum reaction rate, denoted as , must be very similar for all the heating rates. As can be observed in Table 3, p showed an average value of 0.41 ± 0.01, and all the values corresponding to different heating rates are quite similar. Therefore, in a first instance, the Kissinger method could be applied. Figure 5 shows the typical plot representation of the Kissinger method, showing good linear fitting. In particular, the Ea obtained from the slope was 57.3 ± 3.1 kJ mol −1 . As it has been indicated previously, the peak width of the curing rate is also representative for a single step process. For this reason, a plot of ln tFWHM versus 1/Tp has been used to check this. The experimental data and the corresponding linear fit can be seen in Figure 6. The Ea value obtained from the peak width at half maximum was 57.4 ± 3.1 kJ mol −1 . This value is in good agreement with that obtained using the Kissinger method, thus suggesting the curing process took place in a single step process. The reactions that take place during a polymerization process could be quite complex due to the polymer chemical structure, diffusion phenomena, crosslinking density, monomer functionality, and chain mobility, among others [50]. Despite this, the overall process for this particular system, analysed using differential scanning calorimetry, appears to be a single step process (quite homogeneous and symmetric exothermic peak). Therefore, although there could be different activation energies for different processes occurring in the crosslinking, by using DSC, it is possible to obtain a unique "apparent" activation energy, Ea, representative for all the processes that are overlapped. One of the most important parameters of a kinetic study is the apparent activation energy, E a . As indicated previously, the Kissinger method allows for estimating a single value of E a for the curing process by using Equation (15). One important condition needed to apply the Kissinger method is that the conversion at the maximum reaction rate, denoted as α p , must be very similar for all the heating rates. As can be observed in Table 3, α p showed an average value of 0.41 ± 0.01, and all the values corresponding to different heating rates are quite similar. Therefore, in a first instance, the Kissinger method could be applied. Figure 5 shows the typical plot representation of the Kissinger method, showing good linear fitting. In particular, the E a obtained from the slope was 57.3 ± 3.1 kJ mol −1 . As it has been indicated previously, the peak width of the curing rate is also representative for a single step process. For this reason, a plot of ln ∆t FWHM versus 1/T p has been used to check this. The experimental data and the corresponding linear fit can be seen in Figure 6. The E a value obtained from the peak width at half maximum was 57.4 ± 3.1 kJ mol −1 . This value is in good agreement with that obtained using the Kissinger method, thus suggesting the curing process took place in a single step process. The reactions that take place during a polymerization process could be quite complex due to the polymer chemical structure, diffusion phenomena, crosslinking density, monomer functionality, and chain mobility, among others [50]. Despite this, the overall process for this particular system, analysed using differential scanning calorimetry, appears to be a single step process (quite homogeneous and symmetric exothermic peak). Therefore, although there could be different activation energies for different processes occurring in the crosslinking, by using DSC, it is possible to obtain a unique "apparent" activation energy, E a , representative for all the processes that are overlapped.     Although the Kissinger method is easy to use and gives an idea of the apparent activation energy of the curing process, it gives a single E a value. Therefore, the possible change of E a with the conversion cannot be evaluated using the Kissinger method. To this end, several isoconversional methods have been used. Both differential and integral methods have been used and the corresponding plots are gathered in Figure 7. The Friedman method (Figure 7a) shows the calculated data at different conversions and several heating rates. As indicated by Equation (10), the E a can be obtained through the slope of the linear fits. The Pearson's correlation coefficient r for all the linear fits at different conversions was about −0.995. The average E a obtained with the Friedman method was 61.2 ± 2.6 kJ mol −1 . This value is within the same range of the E a obtained using the Kissinger method. In a first approach, the Friedman method did not use any approximation, so it should be more accurate than integral methods since they use different approximations of the temperature integral. However, due to a noisy signal, errors related to the baseline of DSC curves, and so on, there is no great difference between the Friedman method and other integral methods [43]. The Flynn-Wall-Ozawa (FWO) method is an integral method that uses a quite crude approximation of the temperature integral. Nevertheless, it is widely used and gives quite accurate values of E a as seen in Figure 7b. After applying Equation (12), the average E a obtained using the FWO method was 62.8 ± 3.1 kJ mol −1 with excellent correlation coefficients of r ≈ −0.995. The low standard deviation was representative of a slight change in E a with the conversion. A more accurate approximation of the temperature integral is used by the Kissinger-Akahira-Sunose (KAS) method (see Figure 7c). The correlation factor was also close to −0.995 for all linear fits. The average E a obtained using the KAS method was 59.8 ± 3.3 kJ mol −1 , which is similar to all the previously reported values. The Starink method is an integral method that uses a very accurate approximation of the temperature integral and gave r values of about −0.995 (Figure 5d). By using the linear fits, the average E a was 60.0 ± 3.3 kJ mol −1 . Although the Kissinger method is easy to use and gives an idea of the apparent activation energy of the curing process, it gives a single Ea value. Therefore, the possible change of Ea with the conversion cannot be evaluated using the Kissinger method. To this end, several isoconversional methods have been used. Both differential and integral methods have been used and the corresponding plots are gathered in Figure 7. The Friedman method (Figure 7a) shows the calculated data at different conversions and several heating rates. As indicated by Equation (10), the Ea can be obtained through the slope of the linear fits. The Pearson's correlation coefficient r for all the linear fits at different conversions was about −0.995. The average Ea obtained with the Friedman method was 61.2 ± 2.6 kJ mol −1 . This value is within the same range of the Ea obtained using the Kissinger method. In a first approach, the Friedman method did not use any approximation, so it should be more accurate than integral methods since they use different approximations of the temperature integral. However, due to a noisy signal, errors related to the baseline of DSC curves, and so on, there is no great difference between the Friedman method and other integral methods [43]. The Flynn-Wall-Ozawa (FWO) method is an integral method that uses a quite crude approximation of the temperature integral. Nevertheless, it is widely used and gives quite accurate values of Ea as seen in Figure 7b. After applying Equation (12), the average Ea obtained using the FWO method was 62.8 ± 3.1 kJ mol −1 with excellent correlation coefficients of r ≈ −0.995. The low standard deviation was representative of a slight change in Ea with the conversion. A more accurate approximation of the temperature integral is used by the Kissinger-Akahira-Sunose (KAS) method (see Figure 7c). The correlation factor was also close to −0.995 for all linear fits. The average Ea obtained using the KAS method was 59.8 ± 3.3 kJ mol −1 , which is similar to all the previously reported values. The Starink method is an integral method that uses a very accurate approximation of the temperature integral and gave r values of about −0.995 (Figure 5d). By using the linear fits, the average Ea was 60.0 ± 3.3 kJ mol −1 . All the above-mentioned methods give an apparent activation energy between 57 and 62 kJ·mol −1 , which is in total accordance with typical E a values for DGEBA-based epoxy resins, thus suggesting the plant-derived reactive diluent did not affect in a remarkable way the curing process α of the resin [11]. As suggested by the standard deviation of E a from different methods, there is very low variation of E a with α, as can be seen in Figure 8. All the above-mentioned methods give an apparent activation energy between 57 and 62 kJ•mol −1 , which is in total accordance with typical Ea values for DGEBA-based epoxy resins, thus suggesting the plant-derived reactive diluent did not affect in a remarkable way the curing process of the resin [11]. As suggested by the standard deviation of Ea from different methods, there is very low variation of Ea with , as can be seen in Figure 8. Due to this low variability of Ea, it can be considered that the curing process of this partially biobased epoxy resin follows a single step process with a global "apparent" activation energy defined by the average for all the conversion range [57]. Once the Ea has been accurately obtained, it is important to obtain the remaining kinetic parameters to complete the kinetic triplet. This means obtaining the reaction model f( ), and the pre-exponential factor, A. As indicated previously, Málek and Criado-Maqueda described a methodology to compare the calculated y(α) and z(α) , as indicated by Equation (16) and Equation (17), respectively, with different master plots typical of thermally activated processes, as summarized in Table 1. Figure 9 shows the plot of y(α) as a function of the conversion α for different heating rates. As it can be seen, the maximum of the y(α), denoted as , was not zero ( = 0 could be representative for an n th order reaction model), thus suggesting another reaction model. In fact, this maximum was located close to 0.1 and the typical geometry suggested an autocatalytic reaction model as obtained by comparison with different master plots [51,55,58]. The maximum of the calculated y(α) functions for different heating rates are summarized in Table 4. Due to this low variability of E a , it can be considered that the curing process of this partially bio-based epoxy resin follows a single step process with a global "apparent" activation energy defined by the average for all the conversion range [57]. Once the E a has been accurately obtained, it is important to obtain the remaining kinetic parameters to complete the kinetic triplet. This means obtaining the reaction model f (α), and the pre-exponential factor, A. As indicated previously, Málek and Criado-Maqueda described a methodology to compare the calculated y(α) and z(α), as indicated by Equation (16) and Equation (17), respectively, with different master plots typical of thermally activated processes, as summarized in Table 1. Figure 9 shows the plot of y(α) as a function of the conversion α for different heating rates. As it can be seen, the maximum of the y(α), denoted as α M , was not zero (α M = 0 could be representative for an nth order reaction model), thus suggesting another reaction model. In fact, this maximum was located close to 0.1 and the typical geometry suggested an autocatalytic reaction model as obtained by comparison with different master plots [51,55,58]. The maximum of the calculated y(α) functions for different heating rates are summarized in Table 4. Figure 10 gathers the z(α) plots for different heating rates. As can be seen, the z(α) functions show the same geometry, and some of them are overlapped. In can be clearly seen that the maximum of the z(α) function, denoted as α ∞ p , was between 0.4 and 0.5. The actual values for the different heating rates were obtained from the experimental values of z(α) and are shown in Table 4. By comparison of the experimental z(α) plots with some generalized master plots [58], this particular geometry is consistent with an autocatalytic model with a reasonable weight of the autocatalytic effect.  Figure 10 gathers the z(α) plots for different heating rates. As can be seen, the z(α) functions show the same geometry, and some of them are overlapped. In can be clearly seen that the maximum of the z(α) function, denoted as α p ∞ , was between 0.4 and 0.5. The actual values for the different heating rates were obtained from the experimental values of z(α) and are shown in Table 4. By comparison of the experimental z(α) plots with some generalized master plots [58], this particular geometry is consistent with an autocatalytic model with a reasonable weight of the autocatalytic effect.   Figure 10 gathers the z(α) plots for different heating rates. As can be seen, the z(α) functions show the same geometry, and some of them are overlapped. In can be clearly seen that the maximum of the z(α) function, denoted as α p ∞ , was between 0.4 and 0.5. The actual values for the different heating rates were obtained from the experimental values of z(α) and are shown in Table 4. By comparison of the experimental z(α) plots with some generalized master plots [58], this particular geometry is consistent with an autocatalytic model with a reasonable weight of the autocatalytic effect.  The peak for y(α) functions, α M , was 0.090 ± 0.011, and this means there was a certain autocatalytic effect on the curing reaction [59,60]. In fact, one of the conditions for an autocatalytic process requires α M > 0. The autocatalytic effect will be more intense with higher α M values. With regard to the maximum of the z(α) plots, the average α ∞ p was 0.436 ± 0.011. Another condition for an autocatalytic process is that α M < α ∞ p , and obviously, this is true. Therefore, it is possible to use an autocatalytic reaction model to obtain the remaining parameters for the kinetic triplet. For the curing of epoxies, the two parameters Šesták-Berggren reaction model [SB(m,n)] [51] can be used as suggested by both y(α) and z(α) functions. The Šesták-Berggren function is shown in Equation (25) and considers two parameters, n and m. The n parameter represents the typical nth order reaction models as it indicates the reaction rate is proportional to the unreacted material (1 − α), while the m exponent is represents the autocatalytic effect as it indicates, the conversion rate is proportional to the reacted material, α: Once the reaction model has been checked, Equation (25) can be substituted into Equation (3), thus leading to Equation (26): As the apparent activation energy is known (isoconversional methods), it is possible to write Equation (26) as: It is possible to obtain the n and m exponents using an iterative linear fit, but in this work, a non-linear curve fitting has been used. By applying natural logarithms to both sides of Equation (27), we obtain: ln β dα dT + E a RT = ln(A) + m ln α + n ln(1 − α) The term on the left-hand side of Equation (28) is ln(y(α)), such that a non-linear curve fitting, as indicated in Equation (29), can be used to obtain the optimum m and n exponents, as well as the pre-exponential factor, A: The parameters obtained after the non-linear curve fitting are summarized in Table 5. The average value of ln(A) was 18.47 ± 0.05, and regarding the reaction model exponents, m = 0.15 ± 0.01 and n = 1.76 ± 0.08. The low standard deviation indicates these values did not change in a remarkable way with the heating rate.
For an autocatalytic reaction model, the maximum of y(α) is related to the n and m exponents via Equation (30) [51,61]: Substitution of the average n and m values obtained with the non-linear curve fitting led to a theoretical α M of about 0.08, which is in total accordance with the calculated value from experimental data. Once the kinetic triple had been determined, the model wa integrated. Figure 11 shows a comparison of the experimental data (symbols) and the corresponding theoretical models (lines). As it can be seen, Equation (3), with the obtained kinetic triplet (E a , A, and f (α)), gave good accuracy relative to the experimental data for all the heating rates, thus giving consistency to the obtained parameters.
fitting, corresponding to the crosslinking of a partially bio-based epoxy resin for different heating rates. For an autocatalytic reaction model, the maximum of y(α) is related to the n and m exponents via Equation (30) [51,61]: Substitution of the average n and m values obtained with the non-linear curve fitting led to a theoretical α M of about 0.08, which is in total accordance with the calculated value from experimental data.
Once the kinetic triple had been determined, the model wa integrated. Figure 11 shows a comparison of the experimental data (symbols) and the corresponding theoretical models (lines). As it can be seen, Equation (3), with the obtained kinetic triplet (Ea, A, and f( )), gave good accuracy relative to the experimental data for all the heating rates, thus giving consistency to the obtained parameters.

Conclusions
The cure kinetics of a partially bio-based epoxy resin derived from DGEBA and a plant-based reactive diluent (representing 31% of the cured resin) was analyzed using dynamic DSC at different heating rates. The apparent activation energy, Ea, was determined using different isoconversional

Conclusions
The cure kinetics of a partially bio-based epoxy resin derived from DGEBA and a plant-based reactive diluent (representing 31% of the cured resin) was analyzed using dynamic DSC at different heating rates. The apparent activation energy, E a , was determined using different isoconversional methods, and it was found to be between 57 and 62 kJ mol −1 . The low dispersion of E a values for different conversion, α, values suggested a relatively independence of E a from α.
The reaction model was obtained via comparison of the calculated y(α) and z(α) with some well-established reaction models of thermally activated processes. It was concluded that the two-parameter [SB(m,n)] autocatalytic kinetic model of Šesták-Berggren was the most suitable for the description of the curing process of this partially bio-based epoxy resin. Therefore, this work provides a methodology to obtain the kinetic triplet corresponding to the curing/crosslinking of an epoxy resin.