Kinetic Analysis of the Thermal Degradation of Recycled Acrylonitrile-Butadiene-Styrene by non-Isothermal Thermogravimetry

This work presents an in-depth kinetic study of the thermal degradation of recycled acrylonitrile-butadiene-styrene (ABS) polymer. Non-isothermal thermogravimetric analysis (TGA) data in nitrogen atmosphere at different heating rates comprised between 2 and 30 K min−1 were used to obtain the apparent activation energy (Ea) of the thermal degradation process of ABS by isoconversional (differential and integral) model-free methods. Among others, the differential Friedman method was used. Regarding integral methods, several methods with different approximations of the temperature integral were used, which gave different accuracies in Ea. In particular, the Flynn-Wall-Ozawa (FWO), the Kissinger-Akahira-Sunose (KAS), and the Starink methods were used. The results obtained by these methods were compared to the Kissinger method based on peak temperature (Tm) measurements at the maximum degradation rate. Combined Kinetic Analysis (CKA) was also carried out by using a modified expression derived from the general Sestak-Berggren equation with excellent results compared with the previous methods. Isoconversional methods revealed negligible variation of Ea with the conversion. Furthermore, the reaction model was assessed by calculating the characteristic y(α) and z(α) functions and comparing them with some master plots, resulting in a nth order reaction model with n = 1.4950, which allowed calculating the pre-exponential factor (A) of the Arrhenius constant. The results showed that Ea of the thermal degradation of ABS was 163.3 kJ mol−1, while ln A was 27.5410 (A in min−1). The predicted values obtained by integration of the general kinetic expression with the calculated kinetic triplet were in full agreement with the experimental data, thus giving evidence of the accuracy of the obtained kinetic parameters.


Introduction
Thermal degradation of polymers is a topic of great relevance. On one hand, due to their high sensitivity to temperature, degradation or aging at moderate temperatures is responsible for a dramatic loss in overall performance [1,2]. On the other hand, thermal degradation at high temperatures, studied by thermogravimetric analysis (TGA), gives a sound idea of the thermal stability of polymers and composites [3,4]. In addition to conventional analytical techniques, theoretical studies are

Theoretical Background
Kinetic studies are based on the following general kinetic expression: dα dt = k f (α) (1) This expression indicates that the reaction or conversion rate (dα/dt) is directly related to the rate or kinetic constant (k), and a function of the conversion f (α), which takes several forms, depending on the main mechanism, and subsequently, is representative for the reaction model. The rate constant (k) is temperature-dependent and it is assumed to follow the Arrhenius equation in thermally activated reactions in the solid state (Equation (2)): where A (in min −1 ) corresponds to the pre-exponential or frequency factor, E a (in J mol −1 ) stands for the apparent activation energy, R is the gas constant (8.314 J mol −1 K −1 ), and T is the absolute temperature (in K). By combining Equations (1) and (2), the following expression is obtained: As Equation (3) suggests kinetic analysis assists on determining the kinetic triplet: A, E a , f (α).

Isoconversional Methods
Isoconversional methods have been extensively used to estimate E a accurately by using differential or integral methods (in either isothermal or non-isothermal conditions). Differential methods are based on Equation (3), which can be used for whatever thermal program. This time-domain expression can be converted into a temperature-domain by assuming the heating rate (β).
Then, by substituting Equation (4) in Equation (3) yields Equation (5). This temperaturedependent equation is applicable to dynamic thermal programs. β dα dT = A exp(−E a /RT) f (α) (5) Equation (5) can also be represented in its integral form to give Equation (6), which is the base for the integral methods.
Isoconversional methods derive from the basic isoconversional principle that assumes that the reaction rate at a constant conversion (α c ) is only dependent on temperature. By taking natural logarithms on both sides of Equation (1), it yields: ln dα dt = ln k(T) + ln f (α) (7) Polymers 2019, 11, 281 4 of 23 Then, by taking partial derivatives with respect to (1/T) on Equation (7) at a constant α = α i , the following expression is obtained: In isoconversional conditions, at a constant α value α i , its corresponding f (α) value is also constant, and therefore, the second term on the right-hand side of Equation (8) equals zero. The first term on the right-hand side of Equation (8) can be obtained by taking the partial derivative of the natural logarithmic form of Equation (2) with respect to 1/T. Therefore, the following expression can be obtained: Then, by combining Equation (9) and Equation (8), the isoconversional methods are based on the fact that it is possible to obtain an estimation of E a (at a particular α i = constant) without assuming any reaction model, as shown in Equation (10). Thus, these methods are habitually called Model-Free Kinetic methods (MFK). These methods are widely used to study the apparent activation energy of thermal activated processes (degradation and cross-linking) [23].
Although isoconversional methods can give quite accurate values of E a by using differential methods, such as Friedman, or integral methods, such as the Flynn-Wall-Ozawa (FWO), Kissinger-Akahira-Sunose (KAS), and Starink, among others, no additional information is obtained about the reaction model and further methods are then required.

Differential Isoconversional Methods
The Friedman method is a differential method based on Equation (5). By taking natural logarithms and re-arranging terms, Equation (11) is obtained, which suggests that a plot of ln β dα dT against 1/T, at a constant α, for different heating rates β, gives a linear plot whose slope is E a R . As f (α) is constant at a constant α value, and A is also a constant, this approach can give an accurate value of E a without assuming any particular reaction model.

Integral Isoconversional Methods
Integral methods are based on the integration of Equation (6) to give the integral form, g(α) of the reaction model, f (α).
The temperature integral (right-hand side) does not have an analytical solution, so that if we consider the reduced temperature as x = E a /RT, a polynomial P(x) must be obtained by numerical methods. g(α) = AE a βR P(x) (13) Polymers 2019, 11, 281 5 of 23 As indicated previously, by taking x = E a /RT as the reduced temperature, a series expansion can give different approximations of P(x). Asymptotic expansions after a single integration-by-parts are interesting approximations of P(x) [24].
Another asymptotic expansion series is that proposed by Schlomilch, as shown in Equation (15): Or the Lyon's expansion shown in Equation (16) [25]: There are several integral methods which use different solutions of the temperature integral, such as the Broido method [26]. The Flynn-Wall-Ozawa (FWO) method allows estimating the apparent activation energies (E a ) over the entire α range from TGA data at different heating rates (β). This method is widely used because it provides an accurate value of E a for a particular α (isoconversional) value. The FWO method is an integral method that uses the Doyle's approximation of the temperature integral P(x). The empirical equation proposed by Doyle for the temperature integral is shown in Equation (17) [27][28][29][30]: where B(x) is comprised between −1.195 and −1.034 over the x domain. Usually, an average value of −1.052 is used for integral methods for x values comprised between 10 and 60 with a relatively low error. By substituting Equation (17) in Equation (13), it yields: By taking natural logarithms on both sides of the equality and rearranging terms, Equation (19) is obtained. This equation is the base of the FWO integral method: As it can be seen from Equation (19), a plot of log(β) versus 1/T (or ln(β) versus 1/T) allows a linear fit, whose slope is directly related to the E a at that particular α value [31].
Different approximations give different expressions, which can be useful to estimate the values of E a . The FWO method is not the most accurate, as it makes use of a quite essential approximation of P(x). Among others, the Kissinger-Akahira-Sunose method (KAS) employs the Murray and White approximation of P(x) as indicated in Equation (20), which gives more accurate results: Polymers 2019, 11, 281 6 of 23 By using the KAS method (Equation (21)), a plot representation of ln β T 2 versus 1/T for a particular conversion α, gives a straight line with a slope consisting on E a /R [32].
Another approximation is that given by Starink [24] in Equation (22), which allows an even more accurate estimation of E a .
As one can realize, Equation (19), Equation (21), and Equation (22) show similar mathematical structures as Starink [24] already demonstrated, thus leading to a generalized expression: where B and C are experimental values related to a particular approximation of the temperature integral, so that Doyle's approximation considers B = 0 and C = 1.052. The Murray and White approximation leads to B and C values of 2 and 1, respectively.

Combined Kinetic Analysis (CKA)
The above-mentioned methods allow estimating the apparent activation energies without any assumption of the reaction model. As it has been described previously, kinetic analysis offers several tools to determine the so-called kinetic triplet, that is, E a , A, and f(α). One interesting approach to get more information and simultaneously obtain E a and A is the Combined Kinetic Analysis (CKA). Equation (3) can be written in the following form by taking natural logarithms on both sides and rearranging terms: By choosing the appropriate kinetic model, f (α), a plot of the left-hand side of Equation (24) versus 1/T, allows determining E a from the slope of the linear fit while the pre-exponential factor, A, can be obtained from the intercept once E a has been obtained. Obviously, the quality of the linear fit will depend on the selected reaction model, f (α). Table 1 summarizes several kinetic models with their typical f (α) functions and their corresponding integral forms, g(α). Controlled reaction on the surface (migration volume; contracting sphere) Random nucleation with one nucleus of individual particle. Mample first order.
(1−α) −ln(1−α) F 2 Random nucleation with two nuclei of individual particle (1−α) 2 (1−α) −1 F 3 Random nucleation with three nuclei of individual particle In addition to the correct selection of the kinetic model, experimental data can deviate from the theoretical kinetic models due to sample size, geometry, shape, and so on. These drawbacks can be overcome or minimized by the procedure described by Perez-Maqueda et al. [33], which suggests the use of a generic equation derived from the generalized Sestak-Berggren equation shown in Equation (25) [34]. Depending on the combination of the superscripts, m, n and p, Equation (25) can represent almost any reaction model.
The generalized expression proposed by Perez-Maqueda et al. [33] is shown in Equation (26), and as it has been demonstrated, it can fit most of the typical reaction models by simply adjusting the c, m, and n parameters.
If Equation (26) is included into Equation (24) and the terms are rearranged, it gives Equation (27). An optimization of the Pearson's linear correlation coefficient can be considered as the objective to obtain the optimum m and n values that allow a linear fit between the left-hand side part of Equation (27) against 1/T. Once m and n are obtained, the E a and ln cA values can be obtained from the slope and intercept of the optimum linear fit, respectively, as reported by Sánchez-Jiménez et al. [35].

Criado and Maqueda Method
Criado-Maqueda methodology is based on the use of master plots corresponding to different reaction models (Table 1). By comparison with experimental data, it is possible to elucidate the reaction model if it is unique, or in some cases, if it is based on two or more consecutive reaction models, depending on the conversion, α. As previously stated, isoconversional models allow estimating the E a in an accurate way but they do not give any information about the reaction model. Once E a has been obtained, a y(α) function can be defined as shown in Equation (28), with x = E a /RT representing the reduced temperature: Substituting Equation (5) in Equation (28) and rearranging terms, it yields Equation (29). Since A is a constant value, the plot of the second member in Equation (29) against α should be proportional to f (α). Therefore, y(α) gives the geometry of f (α) and it can be very useful to elucidate the reaction model. It is possible to normalize y(α) to enable comparison with the normalized theoretical functions in Table 1. (29) To give more consistency to this method, Criado-Maqueda also proposed the use of a z(α)-type function to represent the master curves corresponding to different reaction models that take into account both f (α) and its integral form g(α). The z(α) function is defined in Equation (30), in which π(x) is related to the Arrhenius temperature integral: As observed in Equation (14), Equation (15), and Equation (16), the approximations of the temperature integral are defined by a series of asymptotic functions. Nevertheless, all these equations can be expressed for a particular expansion, as indicated in Equation (31): Senum and Yang have proposed different approximations of π(x) in a rational form, as summarized in Table 2 [36]. Perez-Maqueda et al. [37] validated the accuracy of these approximations up to eight degrees and they concluded that a π(x) rational approximation with a fourth degree gives an error lower than 10 −5 , therefore the fourth degree is enough for kinetic studies although they have demonstrated increasing accuracy with increasing degree. Table 2 shows the corrected value for the first power x term of the numerator of the fourth degree of 86 instead of 88 as reported by Flynn [38].
Then, by substituting Equation (3) and Equation (31) in Equation (30), the following expression is derived: As the E a of thermally activated processes can be determined by several non-isothermal tests independently of the reaction mechanism f (α), the theoretical master plots obtained by Equation (32) for different well-established reaction mechanisms, as summarized in Table 2, can be easily obtained. These plots can be very useful in determining the main reaction mechanism or even the existence of overlapped reaction mechanisms, as the z(α) corresponding to experimental results can be calculated by using Equation (33), which is derived from including π(x) from Equation (31) on Equation (20) and Or by taking into account the heating rate, β, Equation (33) can be written as: As in the case of y(α), a comparison of the calculated z(α) from the experimental data using Equation (33) or Equation (34) and the different normalized z(α) master plot curves corresponding to different reaction mechanisms for solid state reactions, as described in Table 2, and calculated by Equation (32), can be very useful to assess the reaction model or a combination of them.

Materials
A recycled acrylonitrile-butadiene-styrene (ABS) obtained from waste electric-electronics equipment (WEEE) was used as the base material for kinetic analysis. This ABS grade is a binary blend of SAN and BR. Its characteristic glass transition temperature (T g ) is located at 104.4 • C. The melt flow index (MFI) is 7.7 g/10 min, measured at 210 • C and a load of 5 kg, and 20.5 g/10 min, measured at 230 • C and a load of 5 kg in a ATF Faar S.p.A. melt flow index (MFI) (Novegro-Tregarezzo, Italy). Scheme 1 shows a plot of the structure of this polymer blend.

Thermal Degradation Analysis
Kinetic analysis of the thermal degradation of recycled ABS was carried out by using dynamic thermogravimetry (TGA). To this end, a TGA/SDTA 851 thermobalance from Mettler-Toledo (Schwerzenbach, Switzerland) was used. An average sample mass of 8-10 mg was subjected to a dynamic program from 30 °C up to 700 °C at different heating rates (2,5,10,15,20,25, and 30 K min −1 ) in nitrogen atmosphere with a constant nitrogen flow of 80 mL min −1 . Nitrogen atmosphere was selected to simulate typical "pyrolysis" conditions, which could be helpful for practical uses, such as energy valorization of plastic wastes. Standard graphite crucibles were used. The conversion α at each temperature was calculated using Equation (35): Where w0 denotes the initial mass of the sample, wf stands for the residual mass, and wT represents the mass at a particular temperature (T).

Determination of the Apparent Activation Energy of the Thermal Degradation of ABS
The typical conversion curves for the degradation of recycled ABS are shown in Figure 1a with varying α. Figure 1b shows the corresponding (dα/dT) plots versus T, including the characteristic temperatures for the maximum degradation rates (Tm). A simple method to obtain the Ea from a set of experiments at different heating rates is the Kissinger method [32]. This is widely used because of its simplicity, and it is based on the assumption that d dt dα dt = 0 at Tm. It is assumed that the conversion at the maximum degradation rate, m, is almost constant for a series of heating rates, β. Derivation of Equation (3) gives the following expression: Taking into account the definition of the heating rate β, as indicated in Equation (4) Scheme 1. Schematic plot of the acrylonitrile-butadiene-styrene (ABS) structure, mainly composed of poly(styrene-co-acrylonitrile) (SAN) and poly(butadiene) rubber (BR).

Thermal Degradation Analysis
Kinetic analysis of the thermal degradation of recycled ABS was carried out by using dynamic thermogravimetry (TGA). To this end, a TGA/SDTA 851 thermobalance from Mettler-Toledo (Schwerzenbach, Switzerland) was used. An average sample mass of 8-10 mg was subjected to a dynamic program from 30 • C up to 700 • C at different heating rates (2,5,10,15,20,25, and 30 K min −1 ) in nitrogen atmosphere with a constant nitrogen flow of 80 mL min −1 . Nitrogen atmosphere was selected to simulate typical "pyrolysis" conditions, which could be helpful for practical uses, such as energy valorization of plastic wastes. Standard graphite crucibles were used. The conversion α at each temperature was calculated using Equation (35): where w 0 denotes the initial mass of the sample, w f stands for the residual mass, and w T represents the mass at a particular temperature (T).

Determination of the Apparent Activation Energy of the Thermal Degradation of ABS
The typical conversion curves for the degradation of recycled ABS are shown in Figure 1a with varying α. Figure 1b shows the corresponding (dα/dT) plots versus T, including the characteristic temperatures for the maximum degradation rates (T m ). A simple method to obtain the E a from a set of experiments at different heating rates is the Kissinger method [32]. This is widely used because of its simplicity, and it is based on the assumption that d dt dα dt = 0 at T m . It is assumed that the conversion at the maximum degradation rate, α m , is almost constant for a series of heating rates, β. Derivation of Equation (3) gives the following expression: As indicated previously, the Kissinger methods is based on the assumption d dt dα dt = 0, therefore, by considering Equation (37) and rearranging terms, the following expressions are obtained: Then, by taking natural logarithms on both sides of Equation (39), it yields: For many reaction models, the second term in the right-hand side of Equation (40) is constant [39], thus suggesting that a simple plot of ln β T m 2 versus 1/Tm, at different heating rates β, yields a linear plot whose slope is −Ea/R. If the term - dα is dependent on the heating rate, then the linear relationship cannot be inferred, which represents a limitation of this method [40]. Table 3 summarizes the main results obtained in relation to Tm and the conversion for each Tm value ( m). As it can be seen, the conversion at the peak temperature is almost identical with an average value for the different heating rates of 0.531 ± 0.004.
As one can observe in Figure 2, the correlation between ln β T m 2 versus 1/Tm is very accurate and it allows for estimating a single value of Ea. In particular, the Ea value obtained was 163.1±1.0 kJ mol −1 with a correlation coefficient r = −0.9998. One of the assumptions of this method is that f(α) does not change with α. In fact, it is giving Ea at a particular conversion αm = 0.531. The Kissinger method gives a reliable value of Ea if αm does not change in a significant way with the heating rate, β. On the other hand, the systematic error decreases with increase of the reduced temperature x = Ea/RT. In this particular case, x was comprised between 27.40 and 29.95 for the corresponding Tm values. As it has been reported, x values higher than 10 lead to an error of less than 5% in Ea [41,42]. Taking into account the definition of the heating rate β, as indicated in Equation (4), arranging terms in Equation (36) and considering the definition of dα dt in Equation (3), Equation (36) can be written as: As indicated previously, the Kissinger methods is based on the assumption d dt dα dt = 0, therefore, by considering Equation (37) and rearranging terms, the following expressions are obtained: Then, by taking natural logarithms on both sides of Equation (39), it yields: For many reaction models, the second term in the right-hand side of Equation (40) is dependent on the heating rate, then the linear relationship cannot be inferred, which represents a limitation of this method [40]. Table 3 summarizes the main results obtained in relation to T m and the conversion for each T m value (α m ). As it can be seen, the conversion at the peak temperature is almost identical with an average value for the different heating rates of 0.531 ± 0.004. Table 3. Summary of the thermal parameters corresponding to the temperature for the maximum degradation rate or peak maximum (T m ), and their corresponding conversion (α m ) and conversion rate (dα/dT) m values, for the thermal degradation of recycled ABS at different heating rates β. As one can observe in Figure 2, the correlation between ln β T 2 m versus 1/T m is very accurate and it allows for estimating a single value of E a . In particular, the E a value obtained was 163.1±1.0 kJ mol −1 with a correlation coefficient r = −0.9998. One of the assumptions of this method is that f (α) does not change with α. In fact, it is giving E a at a particular conversion α m = 0.531. The Kissinger method gives a reliable value of E a if α m does not change in a significant way with the heating rate, β. On the other hand, the systematic error decreases with increase of the reduced temperature x = E a /RT. In this particular case, x was comprised between 27.40 and 29.95 for the corresponding T m values. As it has been reported, x values higher than 10 lead to an error of less than 5% in E a [41,42].   On the other hand, due to its simplicity, the Kissinger method can only represent a single-step kinetic process as it is derived from Equation (3). Additional tests are needed to validate the applicability of the Kissinger method. In particular, as indicated by Farjas et al. [43], the peak width of the degradation rate dα dT , is a very sensitive parameter to the existence of multiple transformations.
As indicated, all the thermo-analytical curves differ by a time scale factor τ, which can be defined as the inverse of the Arrhenius constant that is shown in Equation (41): There is a direct relationship between the peak width, measured as the time corresponding to the full width at half maximum (∆tFWHM) of the dα dT curve that is shown in Figure 3a, and the time scale factor corresponding to the peak temperature (Tm), τm (Equation (42)). On the other hand, due to its simplicity, the Kissinger method can only represent a single-step kinetic process as it is derived from Equation (3). Additional tests are needed to validate the applicability of the Kissinger method. In particular, as indicated by Farjas et al. [43], the peak width of the degradation rate dα dT , is a very sensitive parameter to the existence of multiple transformations.
As indicated, all the thermo-analytical curves differ by a time scale factor τ, which can be defined as the inverse of the Arrhenius constant that is shown in Equation (41): There is a direct relationship between the peak width, measured as the time corresponding to the full width at half maximum (∆t FWHM ) of the dα dT curve that is shown in Figure 3a, and the time scale factor corresponding to the peak temperature (T m ), τ m (Equation (42)).
where ∆t FW HM represents the full width at half maximum of a curve obtained when time is normalized to τ m = 1. Therefore, ∆t FW HM depends only on the reaction model and is a constant value, as reported by Farjas et al. [43] Then, substituting τ m in Equation (42), applying natural logarithms and rearranging terms, Equation (42) results in Equation (43), which suggests a linear relationship between ln ∆t FW HM and 1/T m . If ABS degradation follows a single-step kinetic process, then the E a value obtained by Equation (43)   The Ea value obtained after the linear fit between ∆t FWHM vs. 1/Tm at different heating rates was 163.7 ± 1.1 kJ mol −1 , which is in total agreement with the Ea value obtained by applying the Kissinger method. Then, it can be concluded that ABS degradation occurs in a single-step process and the assumptions of the Kissinger method are appropriate, thus giving a reliable Ea value.
Although the premises of the Kissinger method have been checked and the Ea value was accurate, it represents the Ea for a particular conversion αm = 0.531. To evaluate the possible changes in the reaction model during the conversion, isoconversional methods were applied at different conversion values, since they can give much more information about the variability (or not) of Ea with the conversion, α . Figure 4 shows the typical linear fit plots corresponding to different isoconversional methods, both differential and integral. For all the four methods, isoconversional values ranged between α = 0.05 and α = 0.95 at an increasing step of ∆α = 0.05. Figure 4a shows the linear fits in the above mentioned α range according to the Friedman method, which is a differential method. In contrast to the integral methods, which rely on some approximation of the temperature integral, the Friedman method seems to be more accurate, since it does not need any approximation. The E a value obtained after the linear fit between ∆t FW HM vs. 1/T m at different heating rates was 163.7 ± 1.1 kJ mol −1 , which is in total agreement with the E a value obtained by applying the Kissinger method. Then, it can be concluded that ABS degradation occurs in a single-step process and the assumptions of the Kissinger method are appropriate, thus giving a reliable E a value.
Although the premises of the Kissinger method have been checked and the E a value was accurate, it represents the E a for a particular conversion α m = 0.531. To evaluate the possible changes in the reaction model during the conversion, isoconversional methods were applied at different conversion values, since they can give much more information about the variability (or not) of E a with the conversion, α. Figure 4 shows the typical linear fit plots corresponding to different isoconversional methods, both differential and integral. For all the four methods, isoconversional values ranged between α = 0.05 and α = 0.95 at an increasing step of ∆α = 0.05. Figure 4a shows the linear fits in the above mentioned α range according to the Friedman method, which is a differential method. In contrast to the integral methods, which rely on some approximation of the temperature integral, the Friedman method seems to be more accurate, since it does not need any approximation. By using the Friedman method (see Equation (11)), the slope of the linear fit corresponds to −E a /R. The average E a obtained by the Friedman method was 162.9 ± 1.3 kJ mol −1 , which indicates a very low variability of E a with α. All the Pearson's correlation coefficients for the linear fits were at least −0.9995. This narrow value of E a , however, does not necessarily means that ABS degradation is a simple degradation process. In fact, it is a complex process consisting of numerous overlapped consecutive and competitive reactions, though the overall process can be considered as a single-step process with unique leading "apparent" activation energy.  As previously indicated, differential methods do not need any approximation, so that in the first instance they could seem more accurate than integral methods. Nevertheless, the application of these methods to integral data, such as those obtained by TGA, involves calculating the first derivative by numerical methods, with the corresponding inaccuracy being due to its noisy signal. With regard to integral methods, although they use more or less accurate approximations of the temperature integral, they can help in obtaining reliable values for Ea [44]. As indicated by Nikolaidis et al. [45], application of differential and integral methods usually gives differences in Ea but the integral methods become more accurate as the accuracy of the temperature integral increases. The FWO is one of the widely used integral methods, though it uses a quite rude approximation of the temperature integral. As indicated in Equation (19), by using the FWO method, the linear fit of ln (β) As previously indicated, differential methods do not need any approximation, so that in the first instance they could seem more accurate than integral methods. Nevertheless, the application of these methods to integral data, such as those obtained by TGA, involves calculating the first derivative by numerical methods, with the corresponding inaccuracy being due to its noisy signal. With regard to integral methods, although they use more or less accurate approximations of the temperature integral, they can help in obtaining reliable values for E a [44]. As indicated by Nikolaidis et al. [45], application of differential and integral methods usually gives differences in E a but the integral methods become more accurate as the accuracy of the temperature integral increases. The FWO is one of the widely used integral methods, though it uses a quite rude approximation of the temperature integral. As indicated in Equation (19), by using the FWO method, the linear fit of ln(β) vs. 1/T yields a straight line whose slope is −1.052 E a /R. Using this equation, the average E a obtained was 165.7 ± 0.8 kJ mol −1 (with good Pearson's correlation factor r of at least −0.9998), which is in total accordance with the previous results. As the FWO method uses a quite rude approximation of the temperature integral, different corrections have been proposed [24,[46][47][48]. As indicated by Venkatesh et al. [49], the Doyle's approximation can lead to an error of more than 10% in the E a values. If E a is independent of the conversion α, the FWO method gives similar values to those obtained by the Friedman method. As pointed out by Flynn [48], the error related to the Doyle's approximation is less than 1% for reduced temperature values, x = E a /RT, comprised between 21 and 81. The reduced temperature values x for the average temperatures corresponding to an isoconversional α value at different heating rates were between 27 (α iso = 0.95) and 31 (α iso = 0.05), so that the expected error in the FWO method was less than 1%. By applying the corrections described by Flynn [48], the corrected activation energy E a was 163.4 ± 0.8 kJ mol −1 , which is close to the value given by the Friedman method. The KAS method ( Figure 4c) uses a more accurate solution for the temperature integral, thus, in a first instance it gives more accurate values. In a similar way to the above-mentioned methods, the KAS gave an almost independent E a value with the conversion α, with an average value of E a of 163.0 ± 0.8 kJ mol −1 , with r values close to −0.9999 for all isoconversional α values considered. The value obtained for E a is in total agreement with the previous values reported by the Friedman method and it is almost the same value than that obtained by using the Starink method, which uses an even better approximation for the temperature integral. In fact, the linear fits using the Starink method (Figure 4d) resulted in an average value of E a = 163.3 ± 0.8 kJ mol −1 (r values equal or better than −0.9998), thus giving more consistency to the previous differential and integral methods. In addition, the non-dependency observed for E a with the conversion suggests that although the degradation reactions could be quite complex, the overall degradation takes place in a single-step process, and then all the assumptions that have been previously made with the different methods lead to coherent values.
It is not easy to elucidate the reaction model by using isoconversional methods. Huang [50] used the Coats and Redfern integral method with different g(α) functions and this allowed an approach to the reaction model. Nevertheless, in the last years, with the development of advanced computational methods, new approaches have been proposed. One of them is the so-called Combined Kinetic Analysis (CKA) proposed by Criado-Maqueda [33,35,[51][52][53][54][55]. This method makes it possible to obtain the E a value by assuming a general expression for f (α) and maximizing the linear fit between ln as a function of 1/T as indicated in Equation (27). As the slope of these linear fits is negative, then the goal has been minimizing the Pearson's correlation factor (r) by varying the n and m exponents. The results obtained by using this optimization process are summarized in Table 4. It is worthy to note that the linear fit has been applied in a wide conversion range, comprised between α = 0.01 and α = 0.99, and the optimization procedure was possible by using the "minimize" function in MathCad 15 from PTC (Boston, MA, USA).  Figure 5 shows the typical plots, which suggest an excellent linear fit (r = −1.0000) between the combination of experimental data β( dα dT ) and a generic mathematical expression that covers almost all potential reaction models in thermally activated processes (cα m (1 − α) n ). It is worthy to note that the average E a obtained by Combined Kinetic Analysis (CKA) is very similar to those obtained with the previous methods (E a = 163.7 ± 0.2 kJ mol −1 ). Although the use of this modified Sestak-Berggren equation (Equation (26)) is only considered for mathematical purposes, the obtained results for n = 1.4926 ± 0.0005 and m = 0.0032 ± 0.0009 suggest that the (1−α) n term of this mathematical expression contributes in a major extent than the α m term. Therefore, although this method is based only on a mathematical expression that covers all the possible (or a wide variety of) reaction models, it can also suggest the reaction model for the thermal degradation of recycled ABS. One can considerer the n th reaction model as the main reaction model. The extremely low SD obtained for E a , is also representative for its independency from the heating rate, β, which indicates increased consistency to the assumptions of the model, and subsequently, the obtained results are representative. In fact, the E a value given by the CKA is within the same range of all isoconversional methods and the Kissinger method as well.

Determination of the Reaction Model of the Thermal Degradation of ABS
As previously indicated, all the above-mentioned methods give coherent and almost identical Ea values from different approaches. Nevertheless, one step beyond includes determining the remaining kinetic parameters to give a complete kinetic triplet (Ea, A, f(α)) that will be helpful to predict the thermal degradation of recycled ABS and also to check the accuracy between the experimental data

Determination of the Reaction Model of the Thermal Degradation of ABS
As previously indicated, all the above-mentioned methods give coherent and almost identical E a values from different approaches. Nevertheless, one step beyond includes determining the remaining kinetic parameters to give a complete kinetic triplet (E a , A, f (α)) that will be helpful to predict the thermal degradation of recycled ABS and also to check the accuracy between the experimental data and the values predicted by the given kinetic model. To this end, the Criado-Maqueda method [17,[56][57][58] was employed. This considers a comparison of the calculated y(α) and z(α) plots with several well-known reaction models. For this purpose, y(α) and z(α) were calculated as indicated in Equations (28) and (34), respectively. These functions have been calculated for each heating rate and are shown in Figures 6 and 7. As it can be seen in Figure 6, corresponding to y(α) curve plots, the maximum α M for this plot is close to 0 and its geometry suggests a n th order reaction model. Moreover, the different y(α) curves for the different heating rates are almost coincident; y(α) suggests a n th reaction model and the plot of z(α) corroborates this, as can be seen in Figure 7. For a n th order reaction model, the maximum for the z(α) plot, α ∞ p is located at 1 − n ( 1 1−n ) [59], which resulted in an average theoretical n value of 1.4810 by using the values described in Table 5.     By assuming a n th order reaction model, with f(α) = (1−α) n , both the reaction order and ln A values can be simultaneously obtained by linear fit of ln y(α) as a function of ln(1−α). In fact, the slope of this linear fit corresponds to the reaction order n, while the slope represents ln A. In addition to the n th reaction order, the conditions for an autocatalytic process are covered by the calculated y(α) and z(α) plots with α M > 0 and α p ∞ > α M , using the characteristic peak values shown in Table 5. For this reason, the n and m values of the autocatalytic reaction model were calculated by the non-linear fitting of ln y(α) as a function of the conversion α (curve fitting to k + m ln α + n ln(1−α)). Table 6 summarizes the results obtained for both reaction models and their corresponding fitting processes. Application of both reaction models suggests that the n th order reaction model is the best. In fact, when the autocatalytic reaction model is applied, the autocatalytic exponent, m was very low. This can be clearly observed in Figure 8, which shows a comparison of the calculated y(α) functions for different heating rates with different normalized master plots. Figure 8a shows a comparison with By assuming a n th order reaction model, with f (α) = (1−α) n , both the reaction order and ln A values can be simultaneously obtained by linear fit of ln y(α) as a function of ln(1−α). In fact, the slope of this linear fit corresponds to the reaction order n, while the slope represents ln A. In addition to the n th reaction order, the conditions for an autocatalytic process are covered by the calculated y(α) and z(α) plots with α M > 0 and α ∞ p > α M , using the characteristic peak values shown in Table 5. For this reason, the n and m values of the autocatalytic reaction model were calculated by the non-linear fitting of ln y(α) as a function of the conversion α (curve fitting to k + m ln α + n ln(1−α)). Table 6 summarizes the results obtained for both reaction models and their corresponding fitting processes. Application of both reaction models suggests that the n th order reaction model is the best. In fact, when the autocatalytic reaction model is applied, the autocatalytic exponent, m was very low. This can be clearly observed in Figure 8, which shows a comparison of the calculated y(α) functions for different heating rates with different normalized master plots. Figure 8a shows a comparison with typical master plots of the n th order reaction model with different n values from 1.00 (linear) to 2.00 (curved line). The master plot for n = 1.5 gave high coincidence with the calculated y(α) plots. With regard to the autocatalytic master plots, Figure 8b shows different exponent combinations (with n + m = 1.5). As one can observe, if n = m, the maximum is located at a conversion α of 0.5. If m > n, the maximum is moved to higher α values thus indicating prevalence of the autocatalytic term. On the contrary, if m < n, the maximum is displaced to lower α values, thus indicating a minor effect of the autocatalytic contribution. As m becomes lower (e.g., n = 1.4995, m = 0.005), the characteristic plot looks like the typical master plot of the n th order reaction model. Table 6. Main kinetic parameters obtained for the thermal degradation of recycled ABS by using the n th order (linear fit) and the autocatalytic (non-linear curve fit) reaction models obtained for different heating rates. (curved line). The master plot for n = 1.5 gave high coincidence with the calculated y(α) plots. With regard to the autocatalytic master plots, Figure 8b shows different exponent combinations (with n + m = 1.5). As one can observe, if n = m, the maximum is located at a conversion α of 0.5. If m > n, the maximum is moved to higher α values thus indicating prevalence of the autocatalytic term. On the contrary, if m < n, the maximum is displaced to lower α values, thus indicating a minor effect of the autocatalytic contribution. As m becomes lower (e.g. n = 1.4995, m = 0.005), the characteristic plot looks like the typical master plot of the n th order reaction model.  Once the kinetic triplet was obtained, integration of the conversion rate with the corresponding parameters was carried out and its results are shown in Figure 9. The integration was carried out in Mathcad using the fourth-order Runge-Kutta fixed-step method for both n th order and autocatalytic model. Figure 9a shows a comparison plot of the experimental data and theoretical (model) curves obtained by integration of the kinetic model using the n th reaction model with average values obtained from Table 6: n = 1.4950 ± 0.0002, ln A = 27.5410 ± 0.0011. Regarding the autocatalytic model (actually a n th order model too), the average parameters used for the integration were: n = 1.4939 ± 0.0002, m = 0.0050 ± 0.0009, and ln A = 27.5654 ± 0.0009. As it can be seen, there is very good accuracy of the experimental data with both models (actually, the n th order reaction model since the autocatalytic model has a very low contribution of the autocatalytic component). (curved line). The master plot for n = 1.5 gave high coincidence with the calculated y(α) plots. With regard to the autocatalytic master plots, Figure 8b shows different exponent combinations (with n + m = 1.5). As one can observe, if n = m, the maximum is located at a conversion α of 0.5. If m > n, the maximum is moved to higher α values thus indicating prevalence of the autocatalytic term. On the contrary, if m < n, the maximum is displaced to lower α values, thus indicating a minor effect of the autocatalytic contribution. As m becomes lower (e.g. n = 1.4995, m = 0.005), the characteristic plot looks like the typical master plot of the n th order reaction model.  Once the kinetic triplet was obtained, integration of the conversion rate with the corresponding parameters was carried out and its results are shown in Figure 9. The integration was carried out in Mathcad using the fourth-order Runge-Kutta fixed-step method for both n th order and autocatalytic model. Figure 9a shows a comparison plot of the experimental data and theoretical (model) curves obtained by integration of the kinetic model using the n th reaction model with average values obtained from Table 6: n = 1.4950 ± 0.0002, ln A = 27.5410 ± 0.0011. Regarding the autocatalytic model (actually a n th order model too), the average parameters used for the integration were: n = 1.4939 ± 0.0002, m = 0.0050 ± 0.0009, and ln A = 27.5654 ± 0.0009. As it can be seen, there is very good accuracy of the experimental data with both models (actually, the n th order reaction model since the autocatalytic model has a very low contribution of the autocatalytic component).

K min −1 , and
(curved line). The master plot for n = 1.5 gave high coincidence with the calculated y(α) plots. With regard to the autocatalytic master plots, Figure 8b shows different exponent combinations (with n + m = 1.5). As one can observe, if n = m, the maximum is located at a conversion α of 0.5. If m > n, the maximum is moved to higher α values thus indicating prevalence of the autocatalytic term. On the contrary, if m < n, the maximum is displaced to lower α values, thus indicating a minor effect of the autocatalytic contribution. As m becomes lower (e.g. n = 1.4995, m = 0.005), the characteristic plot looks like the typical master plot of the n th order reaction model. Once the kinetic triplet was obtained, integration of the conversion rate with the corresponding parameters was carried out and its results are shown in Figure 9. The integration was carried out in Mathcad using the fourth-order Runge-Kutta fixed-step method for both n th order and autocatalytic model. Figure 9a shows a comparison plot of the experimental data and theoretical (model) curves obtained by integration of the kinetic model using the n th reaction model with average values obtained from Table 6: n = 1.4950 ± 0.0002, ln A = 27.5410 ± 0.0011. Regarding the autocatalytic model (actually a n th order model too), the average parameters used for the integration were: n = 1.4939 ± 0.0002, m = 0.0050 ± 0.0009, and ln A = 27.5654 ± 0.0009. As it can be seen, there is very good accuracy of the experimental data with both models (actually, the n th order reaction model since the autocatalytic model has a very low contribution of the autocatalytic component).
Once the kinetic triplet was obtained, integration of the conversion rate with the corresponding parameters was carried out and its results are shown in Figure 9. The integration was carried out in Mathcad using the fourth-order Runge-Kutta fixed-step method for both n th order and autocatalytic model. Figure 9a shows a comparison plot of the experimental data and theoretical (model) curves obtained by integration of the kinetic model using the n th reaction model with average values obtained from Table 6: n = 1.4950 ± 0.0002, ln A = 27.5410 ± 0.0011. Regarding the autocatalytic model (actually a n th order model too), the average parameters used for the integration were: n = 1.4939 ± 0.0002, m = 0.0050 ± 0.0009, and ln A = 27.5654 ± 0.0009. As it can be seen, there is very good accuracy of the experimental data with both models (actually, the n th order reaction model since the autocatalytic model has a very low contribution of the autocatalytic component).

Conclusions
This work shows a systematic kinetic study of the thermal degradation of recycled ABS in nitrogen atmosphere. The obtained data by dynamic TGA at different heating rates, comprised between 2 and 30 K min −1 , were used to apply different model-free kinetic methods (MFK) based on the isoconversional principle. With regard to the differential methods, the Friedman approach gave an Ea value of 162.9 ± 1.3 kJ mol −1 , which is quite similar to those obtained by applying several integral methods, such as FWO (Ea = 163.4 ± 0.8kJ mol −1 ), KAS (Ea = 163.0 ± 0.8 kJ mol −1 ), and the Starink method (Ea = 163.3 ± 0.8 kJ mol −1 ). In addition, the variation of Ea with the conversion α was almost negligible, thus suggesting the degradation process of ABS takes place in a single-step process. A different approach obtained by CKA gave an Ea value of 163.7 ± 0.2 kJ mol −1 , which is in total agreement with other isoconversional methods, thus giving clear evidence of the consistency of the CKA method. Once the consistency of the Ea values was confirmed, the other kinetic parameters were obtained by calculating the y(α) and z(α) characteristic plots and comparing them with the master curve plots. Both y(α) and z(α) plots suggested a n th order reaction model. Then, both n and ln A were obtained by linear fitting resulting in n = 1.4950 ± 0.0002 and ln A = 27.5410 ± 0.0011 (A in min −1 ). Therefore, an excellent agreement between the experimental data and the values predicted by the kinetic model was found. In addition, the autocatalytic reaction model was tested, and the corresponding n and m exponents were: n = 1.4939 ± 0.0002, m = 0.0050 ± 0.0009, and ln A = 27.5654 ± 0.0009. These results indicated that the autocatalytic component has a very low contribution in the reaction model (extremely low m values) and it also corroborated the good agreement with a simple n th order reaction model.    Once the kinetic triplet was obtained, integration of the parameters was carried out and its results are shown in Figu Mathcad using the fourth-order Runge-Kutta fixed-step met model. Figure 9a shows a comparison plot of the experimen obtained by integration of the kinetic model using the n obtained from Table 6: n = 1.4950 ± 0.0002, ln A = 27.5410 ± 0.0 (actually a n th order model too), the average parameters use 0.0002, m = 0.0050 ± 0.0009, and ln A = 27.5654 ± 0.0009. As it of the experimental data with both models (actually, th autocatalytic model has a very low contribution of the autoca 5 K min −1 , and netic triplet was obtained, integration of the conversion rate with the corresponding carried out and its results are shown in Figure 9. The integration was carried out in he fourth-order Runge-Kutta fixed-step method for both n th order and autocatalytic a shows a comparison plot of the experimental data and theoretical (model) curves egration of the kinetic model using the n th reaction model with average values able 6: n = 1.4950 ± 0.0002, ln A = 27.5410 ± 0.0011. Regarding the autocatalytic model der model too), the average parameters used for the integration were: n = 1.4939 ± 0 ± 0.0009, and ln A = 27.5654 ± 0.0009. As it can be seen, there is very good accuracy ental data with both models (actually, the n th order reaction model since the del has a very low contribution of the autocatalytic component).

Conclusions
This work shows a systematic kinetic study of the thermal degradation of recycled ABS in nitrogen atmosphere. The obtained data by dynamic TGA at different heating rates, comprised between 2 and 30 K min −1 , were used to apply different model-free kinetic methods (MFK) based on the isoconversional principle. With regard to the differential methods, the Friedman approach gave an E a value of 162.9 ± 1.3 kJ mol −1 , which is quite similar to those obtained by applying several integral methods, such as FWO (E a = 163.4 ± 0.8kJ mol −1 ), KAS (E a = 163.0 ± 0.8 kJ mol −1 ), and the Starink method (E a = 163.3 ± 0.8 kJ mol −1 ). In addition, the variation of E a with the conversion α was almost negligible, thus suggesting the degradation process of ABS takes place in a single-step process. A different approach obtained by CKA gave an E a value of 163.7 ± 0.2 kJ mol −1 , which is in total agreement with other isoconversional methods, thus giving clear evidence of the consistency of the CKA method. Once the consistency of the E a values was confirmed, the other kinetic parameters were obtained by calculating the y(α) and z(α) characteristic plots and comparing them with the master curve plots. Both y(α) and z(α) plots suggested a n th order reaction model. Then, both n and ln A were obtained by linear fitting resulting in n = 1.4950 ± 0.0002 and ln A = 27.5410 ± 0.0011 (A in min −1 ). Therefore, an excellent agreement between the experimental data and the values predicted by the kinetic model was found. In addition, the autocatalytic reaction model was tested, and the corresponding n and m exponents were: n = 1.4939 ± 0.0002, m = 0.0050 ± 0.0009, and ln A = 27.5654 ± 0.0009. These results indicated that the autocatalytic component has a very low contribution in the reaction model (extremely low m values) and it also corroborated the good agreement with a simple n th order reaction model.