Monte Carlo Simulation of Linear Polymer Thermal Depolymerization under Isothermal and Dynamic Modes

Abstract

Kinetics of linear polymer thermal depolymerization under isothermal and dynamic TGA modes was simulated by the Monte Carlo method. The simulation was carried out on model arrays having the same initial degree of polymerization 𝑃_𝑛=100 and different width (polydispersity index, PDI=𝑃_𝑤/𝑃_𝑛=1∼3) at three constant temperatures and five heating rates. Kinetics of the process in both modes is described by the Avrami equation, the exponent in which decreasing as the distribution width increases. Treatment of the model kinetic curves of degradation using the nonlinear regression method by the Avrami equation, under both isothermal and dynamic modes, gives correct activation energy and pre-exponential factor values independently of the initial PDI. Data obtained in the dynamic mode were also treated by two isoconversion methods, widely applied to kinetic analysis of TGA curves (Flynn-Wall-Ozawa method and Kissinger-Akahira-Sunose (KAS) method).

1. Introduction

Thermogravimetric analysis (TGA) method is widely applied to the investigation of various polymers and polymeric composites. This is a simple and reliable method that allows estimation of thermal resistance of various materials and obtaining of information on kinetics and the mechanisms of thermal degradation processes. At present, different investigators are using about 20 various methods for quantitative treatment of the TGA thermograms for determination of thermal degradation kinetic parameters, the activation energy, first of all.

The most general description of a single-stage thermal degradation of a polymer under the dynamic TGA mode looks as follows:𝑑𝛼=1𝑑𝑇𝛽⋅𝐴⋅exp−𝐸(𝛼)𝑅𝑇⋅𝑓, (1) where 𝛼 is the conversion determined from the expression 𝛼=(𝑚₀−𝑚)/(𝑚₀−𝑚_∞); 𝑚, 𝑚₀, and 𝑚_∞ are current, initial, and final masses of the sample, respectively; 𝛽 is the heating rate; 𝑇 is the absolute temperature; 𝑅 is the gas constant; 𝐴 and 𝐸 are the pre-exponential factor and the activation energy of degradation, respectively. The kinetic function 𝑓(𝛼) depends on a particular degradation mechanism.

The majority of the methods applied to kinetic analysis of TGA thermograms use different approximate solutions of (1) in the form of linear approximations. TGA data may also be treated by the method of nonlinear regression [1, 2], (1) being solved numerically. The main kinetic functions 𝑓(𝛼), used for description of thermal processes in the solid phase, are classified in [3]. Unfortunately, all these kinetic descriptions are rather formal and can hardly be bound to the real physicochemical processes proceeding in the polymers.

One of the frequently observed kinetic mechanisms of thermal degradation of polymers is depolymerization in which monomer molecules are severed one by one from the chain ends. In accordance with the literature, degradation of polymers, for example, as poly (methyl methacrylate) [4], polystyrene [5], and so forth, proceeds by this mechanism. Depolymerization of a single macromolecule having the initial degree of polymerization 𝑛 may be described by a sequence of elementary acts of a monomer severing;(𝑀)_𝑛⟶(𝑀)_𝑛(𝑀)−1+𝑀,_𝑛−1⟶(𝑀)_𝑛(𝑀)−2+𝑀,⋮,₂⟶𝑀+𝑀. (2)

In the simplest case, the monomer severing rate constant, 𝑘, is independent of the chain length. Since severing proceeds by ends only, the monomer formation rate is proportional to the chain concentration in the system, 𝑁:𝑑𝑀𝑑𝑡=𝑘⋅2𝑁. (3)

As a consequence, the conversion 𝛼=𝑀/𝑃 may be written down as follows:𝑑𝛼=𝑘𝑑𝑡𝑃⋅2𝑁, (4) where 𝑃 is the initial degree of polymerization. At the same time, the chain death happens only at the latest stage, when a dimer decays into two molecules of the monomer.

The problem is that the real polymer consists of chain of different lengths. Total kinetic scheme, which includes macromolecules of all possible sizes, is observed extremely bulky. Nevertheless, is seems obvious that integral kinetics of depolymerization shall strictly depend on PDI of the initial polymer.

The Monte Carlo simulation method is widely applied in the polymeric physicochemistry, to kinetics of polymer polymerization [6–9] and degradation [10], in particular. We have used this method for simulation of linear polymer depolymerization kinetics under isothermal and dynamic modes. For this purpose, 8 model arrays, 1000 chains each having the same number average degree of polymerization is 100 and PDI of 1 to 3 were synthesized on a computer. On these arrays, stochastic simulation of the depolymerization process at three temperatures and five heating rates was performed. Moreover, arrays with the highest probable PDI equal 2 and the degree of polymerization is 200 to 500 were synthesized, on which simulation in the isothermal mode were performed. The kinetic curves obtained were analyzed by the nonlinear regression method. Data obtained in the dynamic mode were also treated by two isoconversion methods, widely applied in the literature (Flynn-Wall-Ozawa method [11, 12] and Kissinger-Akahira-Sunose (KAS) method [13, 14]).

2. The Procedure of Simulation of the Initial Arrays of Macromolecules

The initial arrays, each containing 1000 polymeric chains, were simulated on the basis of the lognormal distribution. The array of the integers was separated into fractions from 1 to 10, 11 to 20, and so forth … from 10(𝑖 − 1)+1 to 10𝑖, where 𝑖 is the integer from 1 to 1000. For each 𝑥=10𝑖, density of the lognormal distribution was determined as follows:𝐹(𝑥)=1√2𝜋𝜎⋅𝑒^{−(log𝑥−𝜇)2/2𝜎2}, (5) where 𝜎 and 𝜇—parameters of lognormal distribution.

The number of chains in each fraction, having the degree of polymerization from 𝑥+1 to 𝑥+10, was calculated by the formula:𝑞_𝑖=(𝐹(𝑥+10)−𝐹(𝑥))⋅1000, (6) and was rounded to the nearest integer. The length of each chain within the fraction was determined as a random integer, uniformly distributed in the appropriate interval. The first interval included chains of the degree of polymerization from 2 to 10, that is, it is suggested that the initial polymer contains no monomer. The values of 𝜎 and 𝜇 were selected so that the given average degree of polymerization by the array 𝑃_𝑛 and the required PDI would be obtained. The described procedure is realized as the MS Excel macros. It was used to synthesize the arrays having 𝑃_𝑛=100 and the index PDI equal 1.14, 1.25, 1.5, 1.75, 2, 2.5, and 3. Figure 1 shows distribution histograms by the degree of polymerization for these arrays. Arrays with PDI = 2 and the degree of polymerization equal 200, 300, 400, and 500 were also obtained. Moreover, a monodispersity array (PDI = 1), in which all chains had the same length equal 100, was used in the kinetic experiments.

Figure 1 

Histograms of distribution 𝑟(𝑛) of the initial arrays of chains. Digits in the legend show values of PDI, 𝑃𝑤/𝑃𝑛.

3. The Procedure of Depolymerization Kinetics Simulation

For simulation of depolymerization, the Monte Carlo kinetic method was used [15]. The time corresponded to a single step in the isothermal experiment was calculated by the following formula:𝑡₁=−ln𝑟𝑘⋅𝑁, (7) where 0<𝑟<1 is a random number; 𝑁 is the current number of chains in the array,𝑘=𝐴₀⋅𝑒^{−𝐸/𝑅𝑇}, (8) is the rate constant of the elementary act of the monomer severing; 𝐴₀ and 𝐸 are its pre-exponential factor and the activation energy, respectively; 𝑅 is the gas constant; 𝑇 is the absolute temperature, at which the experiment was simulated. In all experiments, both isothermal and dynamic, 𝐴₀ was set equal 10⁷ min⁻¹, and 𝐸=100 kJ/mol.

In dynamic experiments, the temperature step 𝛿𝑇=𝑇_𝑖−𝑇_𝑖−1 was calculated by the formula:𝛿𝑇=−𝛽ln𝑟𝑘_𝑖⋅𝑁, (9) where 𝑘_𝑖 constant was calculated from (8) at temperature 𝑇_𝑖. If 𝛿𝑇 value exceeded 1°, the current temperature 𝑇_𝑖 was increased by one, and the procedure was repeated.

4. Data Treatment

The results of simulation in the isothermal and the dynamic modes were treated by the nonlinear regression method, using Fitter software [2]. Moreover, for TGA thermograms in the dynamic mode, two widely applied methods based on linear approximations were used.

The Flynn-Wall-Ozawa method [11, 12] represents an isoconversion method, which uses the following equation: ln𝛽=ln0.0048𝐴𝐸𝑔(𝛼)𝑅𝐸−1.0516𝑅𝑇, (10) where 𝑔(𝛼) is an arbitrary conversion function. The activation energy is determined by the temperature dependence, at which the given conversion is reached, on the heating rate. This method allows composing of a degradation profile that is dependence of the activation energy on conversion.

The KAS (Kissinger-Akahira-Sunose) method [13, 14] also relates to isoconversion methods. The relation between temperature, at which the given conversion is reached, and the heating rate is described by the following equation: 𝛽ln𝑇²=ln𝐴𝑅𝑔(𝛼)𝐸−𝐸𝑅𝑇. (11)

The slope of the line composed in the coordinates of (11) gives the activation energy value.

5. Results and Discussion

Simulation of the degradation kinetics under isothermal conditions was performed at three temperatures of 280, 300, and 320°C. Figure 2 shows kinetic curves of the relative mass loss, 𝑦=1−𝛼, at 300°C for three arrays with different initial PDI. It appears that the wider the initial distribution is, the slower degradation proceeds.

Figure 2 

Mass loss kinetic curves at 300°C for arrays of the initial PDI of 1.14 (1), 2 (2), and 3 (3). Points show values obtained by simulation, and continuous lines show treatment results by the nonlinear regression method.

For monodispersion distribution, conversion of 0.7 is reached during 137 hours, and for the widest distribution of 𝑃_𝑤/𝑃_𝑛=3 for 357 hours. Meanwhile, the initial degradation rate is the same for all arrays, because it depends on the initial number of chains only.

For each array, the curves at three temperatures were treated by the nonlinear regression method, using common models of the solid phase kinetics [3]. It has been found that the best description is given by the Avrami equation [16], which we used in the following form:𝑦=exp−(𝑘⋅𝑡)^𝑞. (12)

This form of record is somewhat different from the traditional one, in which only time 𝑡 has the exponent 𝑞, and the effective constant k is factorized. However, such presentation seems to be more correct [17], because in this case, dimensionality of the rate constant equals to the reverse time, independently of the exponent value.

Figure 2 and Tables 1 and 2 show treatment results of the depolymerization isothermal kinetics by the Avrami (continuous lines) and estimates of kinetic parameters, respectively. The activation energy of depolymerization is virtually independent of the PDI of the initial array and the initial degree of polymerization, and correlates well with the value given during simulation. The frequency factor 𝑘₀, obtained during the treatment, is bound to the pre-exponential factor of the elementary act of the monomer severing 𝐴₀ by the following expression:𝑘₀=𝐴₀𝑃_𝑛. (13)

Table 1 

Estimations of the Avrami equation parameters for arrays with the initial degree of polymerization equal 100 at isothermal depolymerization.

Table 2 

Estimations of the Avrami equation parameters for arrays with the initial PDI equal 2 at isothermal depolymerization.

This expression yields from (4). Values of 𝐴₀ calculated from this expression are also shown in Tables 1 and 2. Generally, these values are close to that of 10⁻⁷, given during simulation, although observed some trend growth 𝐴₀ when the PDI of the initial array is rising.

The exponent 𝑞 is independent of the initial degree of polymerization, but monotone decreases, as the initial sequence distribution increases.

Simulation under the dynamic mode was performed at 5 heating rates: 2, 5, 10, 20, and 40°C/min. Figure 3 shows model TGA thermograms for arrays with the initial PDI of 1, 1.5, and 3. In this case, the same situation is observed, as it is for the isothermal degradation: the process decelerates with the distribution width increase. As transited from the narrowest distribution to a wide one, temperature, at which conversion reaches 0.5, increases by 38°C.

Figure 3 

TGA thermograms at the heating rate of 10°/min for arrays of the initial PDI of 1.14 (1), 2 (2), and 3 (3). Points show values obtained by simulation, and continuous lines show treatment results by the nonlinear regression method.

Series of the curves modeled for each array at 5 heating rates were treated by the nonlinear regression, by the following model: 𝑑𝑦𝑑𝑡=−𝑔⋅𝑘^𝑔⋅𝑡^𝑔−1⋅𝑦𝑘=𝑘₀⋅𝑒^{−𝐸/𝑅𝑇}𝑇=𝑇₀+𝛽⋅𝑡 (14) representing differential form of the Avrami equation (12), Figure 3 shows the treatment results (continuous lines), and Table 3 shows estimates of parameters.

Table 3 

Estimations of the Avrami equation parameters for arrays with different initial PDI under dynamic TGA mode of depolymerization.

Thus, the Avrami equation describes well kinetics of depolymerization in both isothermal and dynamic modes. Activation energy estimates obtained from model thermograms are rather close to the value of 100 kJ/mol, given at simulation. The values of other parameters obtained in both isothermal and dynamic modes also correlate well.

Note that the Avrami equation was for many times applied by different investigators to kinetic analysis of the TGA data [18–21]. This equation was deduced for the purpose of describing the crystallization process [16], and is considered to be one of the models of the solid phase kinetics [3]. The mechanism of the processes described by the Avrami equation is often denoted as “nucleation and growth.” In the present work, we simulated the most usual kinetics, following the law of acting masses; no features of the reaction in the solid phase were taken into consideration. In this case, the specific shape of kinetic curves is stipulated by the mechanism of depolymerization itself, representing a series of consecutive acts of the monomer severing from macromolecule ends. In these processes, the chain ends represent active centers, and they die only at the last act of dimer decomposition.

Of course, the Avrami equation gives just an approximate description of depolymerization. The accuracy of this description increases with the distribution width. Obviously, exponent decrease with the increase of sequence distribution width in the Avrami equation is related to an increase of a contribution of longer macromolecules into the general process. Thermograms show this as an increase of thermal resistance, although, truly, kinetic parameters of elementary act of the monomer severing do not change. Values of q parameter may change with respect to particular shape of the sequence distribution. Nevertheless, we suggest that the tendency of q decrease with the increase of PDI of the initial polymer shall be preserved. Figure 4 shows examples of the degradation profiles calculated by data, which were obtained in the dynamic mode by two isoconversion methods. Obviously, the Flynn-Wall-Ozawa method gives a noticeably overestimated activation energy, which also increases with the conversion. Calculation by the KAS method allows obtaining of 𝐸 value, which is more reliable and independent of the conversion. Table 4 shows activation energy estimates for arrays with different initial PDI. These estimates represent the average values in the conversion range of 0.05 to 0.9.

Table 4 

Estimations of activation energy, obtained by isoconversion methods for arrays with different initial PDI.

Figure 4 

Degradation profiles for the array with 𝑃𝑤/𝑃𝑛 = 2 calculated by the Flynn-Wall-Ozawa (1) and the KAS (2) methods. Continuous horizontal line shows 𝐸 value given at the simulation.

Both isoconversion methods give the activation energy values independent of the sequence distribution width. However, estimates obtained by the KAS method are, nevertheless, closer to the true value of 100 kJ/mol, compared with the value obtained by the Flynn-Wall-Ozawa method.

6. Conclusion

The data obtained by stochastic simulation show that the kinetics of thermal polymerization is strictly affected by the sequence distribution width. The increase in PDI causes a seeming increase of thermal stability by the following indices: time of reaching the given conversion in the isothermal mode, or temperature of reaching the given conversion in the dynamic TGA mode, although the activation energy of the process remains unchanged.

The best kinetic model for depolymerization of a linear polymer in both isothermal and dynamic modes is the Avrami equation. Treatment of the kinetic curves for this model by the method of nonlinear regression, under both modes, gives the true activation energy irrespective of the initial PDI. The exponent, q, in the Avrami equation regularly decreases with increasing width of the sequence distribution.

Among two widely used isoconversion methods of TGA data treatment based on linear approximations, at thermal degradation by the depolymerization mechanism, the best results are shown by the KAS method. Treatment by the Flynn-Wall-Ozawa method gives somewhat overestimated activation energies compared with the value set out for simulation.