Molecular Dynamics Simulation of C–C Bond Scission in Polyethylene

Feb 26, 2014 - The reaction of C–C bond scission in polyethylene chains of various lengths was studied using molecular dynamics under the conditions...
0 downloads 0 Views 492KB Size
Download PDF
Article pubs.acs.org/JPCA

Molecular Dynamics Simulation of C−C Bond Scission in Polyethylene and Linear Alkanes: Effects of the Condensed Phase Konstantin V. Popov and Vadim D. Knyazev* Research Center for Chemical Kinetics, Department of Chemistry, The Catholic University of America, Washington, DC 20064, United States S Supporting Information *

ABSTRACT: The reaction of C−C bond scission in polyethylene chains of various lengths was studied using molecular dynamics under the conditions of vacuum and condensed phase (polymer melt). A method of assigning meaningful rate constant values to condensed-phase bond scission reactions based on a kinetic mechanism accounting for dissociation, reverse recombination, and diffusional separation of fragments was developed. The developed method accounts for such condensed-phase phenomena as cage effects and diffusion of the decay products away from the reaction site. The results of C−C scission simulations indicate that per-bond rate constants decrease by an order of magnitude as the density of the system increases from vacuum to the normal density of a polyethylene melt. Additional calculations were performed to study the dependence of the rate constant on the length of the polymer chain under the conditions of the condensed phase. The calculations demonstrate that the rate constant is independent of the degree of polymerization if polyethylene samples of different lengths are kept at the same pressure. However, if instead molecular systems of different polyethylene chain lengths decompose under the conditions of the same density, shorter chains result in higher pressures and lower rate constants. The observed effect is attributed to a higher degree of molecular crowding (lower fraction of free intermolecular space available for molecular motion) in the case of shorter molecules.

1. INTRODUCTION Combustion and flammability of plastics are important topics of practical interest directly related to fire safety. Even though polymers are solid or liquid, their combustion is generally a gasphase phenomenon, where gaseous fuel is supplied by the decomposing solid or liquid polymer. Thus, knowledge of the key features of polymer pyrolysis is essential for understanding of the chemistry and physics of polymer combustion and flammability. Efforts of many researchers are being directed at the elucidation of the main features of the chemical mechanisms of polymer pyrolysis; an important component of these efforts is kinetic modeling of pyrolysis based on assembling chemical models consisting of numerous elementary reactions (e.g., refs 1−5 and references therein). Such detailed kinetic modeling of complex chemical processes has achieved considerable success in recent decades, particularly in its application to such systems as combustion or chemistry of the atmosphere. However, existence of a developed and reliable database of kinetic rate coefficients is a necessary condition for the success of such modeling. Most of the experimental kinetic data available now are limited to reactions that involve relatively small molecules in the gas phase due to the comparative ease of their experimental investigation. Similarly, many theoretical studies concentrate on gas-phase reactions of small molecules, which enables isolating the reaction under study from any solvent influences. Even though theories aimed at quantitative description of solvent effects on reaction rate constants have achieved significant progress (e.g., see a review in ref 6 and references cited therein), the majority of theoretical studies aiming at precise prediction of rate constants are performed for © 2014 American Chemical Society

the gas-phase conditions. Therefore, in the kinetic modeling of such processes as combustion of real fuels (relatively large longchained, branched, and cyclic hydrocarbon molecules) and pyrolysis of polymers kinetic parameters are usually approximated by the gas-phase data obtained experimentally for chemically similar but smaller species. However, the applicability of such an approximation may appear doubtful for a number of types of chemical reactions. One example of such reactions is the backbone scission, i.e., dissociation of C−C bonds forming the polymer chain, one of the most important reaction types in polymer pyrolysis. These reactions are generally believed to be responsible for initiation of polymer pyrolysis. In evaluating the condensed-phase rate constants for polymer C−C bond dissociation, an approach based on analogy with gas-phase small molecule chemistry (decomposition of small alkanes) leads to preexponential factors of approximately 1016 to 1017 s−1 and activation energy values close to the average strength of the C−C bond, ∼350 kJ mol−1. The adequacy of this approach may appear questionable when calculated temperature dependences of the rate constant are considered in light of the known experimental data on polymer pyrolysis. For example, the above Arrhenius parameters predict that the per-bond C−C scission rate constant, kCC, will achieve the value of 1 s−1 at T ≅ 1070− 1140 K and that of 1 × 10−3 s−1 at T ≅ 910−960 K. At the same time the typical temperature range for the onset of polyethylene pyrolysis is 300−400 °C (573−673 K),1,7 where Received: November 21, 2013 Revised: February 25, 2014 Published: February 26, 2014 2187

dx.doi.org/10.1021/jp411474u | J. Phys. Chem. A 2014, 118, 2187−2195

The Journal of Physical Chemistry A

Article

the above Arrhenius parameters would give kCC ≅ 2 × 10−16 to 1 × 10−10 s−1. This disparity between the temperatures of the pyrolysis onset measured experimentally and estimated from the kinetics of decomposition of small alkanes seems to suggest that the rate constants of the reaction of C−C bond scission may have significantly higher values (also, see a discussion of rate constants in ref 1). In addition, under the condensed-phase conditions of the polymer melt, reactions can be expected to be affected by the cage effect,8 which may considerably influence their rates compared to the rates of similar reactions in the gas phase. This is another concern that dictates that extreme caution must be applied when kinetic parameters determined in the gas phase are used for modeling liquid-phase processes. These considerations and concerns make the study of the kinetics of dissociation of long-chain hydrocarbons and polymer backbones an interesting topic of research. Recently, a series of reactive molecular dynamics (RMD) studies of the thermal degradation of alkane-based polymers reported effective activation energies for C−C bond scission that are significantly lower than the corresponding bond strength (∼350 kJ mol−1). An investigation of the pyrolysis of polyethylene, polypropylene, and polyisobutylene by Nyden et al.9 reported activation energies in the range of 40−190 kJ mol−1. Later, in a RMD study of polyisobutylene (PIB) pyrolysis by Stoliarov et al.10 a correlation between the activation energy of C−C bond scission and the length of the PIB molecule was observed. The authors reported that the activation energy reduced from 239 kJ mol−1 for PIB-4 to 170 kJ mol−1 for PIB-150 (here and below the number in the polymer notation represents the number of monomers that the polymer molecule consists of). More recently, a similar effect was observed in an RMD work by Smith et al.11 for the reaction of the C−C bond decay in the polyethylene chain in the condensed phase: the activation energy values lowered from 349 kJ/mol (for PE-1) to 251 kJ/mol (for PE-100). The latter study reported a corresponding significant increase of the perbond rate constant upon the increase in the polyethylene chain length. In a work by Knyazev12 a molecular dynamics study of the reaction of C−C bond decay was performed for a series of linear alkanes and polyethylene (PE) molecules in vacuum, with the alkane chain length ranging from PE-1 to PE-1000. During the study it was discovered that the rate of the reaction increases by up to an order of magnitude with the increase in the alkane chain length. This effect was tentatively attributed to the centrifugal forces of torsional motions of molecular fragments separated by a C−C bond subject to rupture: the forces act on the bond pulling it apart and making the bond’s decay easier. No effects of the polyethylene chain length on the activation energy of C−C bond scission were observed. The study of ref 12 was limited to the gas-phase conditions (single molecule in vacuum). In the current study, we expand our molecular dynamics investigation of the C−C bond scission in alkanes and linear polyethylene to include the effects of the condensed phase. Simulations are performed for the conditions of polyethylene melt, as well as for isolated molecules in vacuum. Effects of the density of the liquid polymer melt on the observed rate constants are studied, as are those of the polyethylene chain length and the pressure.

(PE) molecules with the number of individual monomer units of 1 (ethane), 2 (n-butane), 2.5 (pentane), 25, 100, and 500. To investigate the effects of the condensed phase on the temperature dependence of the PE decomposition rate constants, the reaction of PE-500 was studied at the densities of 0.75 g cm−3 (the density of molten polyethylene13) and 0.50 g cm−3 (henceforth referred to as “reduced density”), as well as under vacuum conditions. These calculations covered the temperature range 2300−3100 K. To study the effects of the PE chain length on the rate constant values, calculations were also performed for shorter molecules, at the single temperature of 3100 K. The reaction of the C−C scission for the shorter molecules was studied under two sets of conditions: (1) the constant density of 0.75 g cm−3 and (2) variable density adjusted for each molecular system to keep a constant pressure of 1.8 × 104 bar, the pressure observed in the PE-500 system at 3100 K, and the density of 0.75 g cm−3. These conditions are referred to henceforth as “constant density” and “constant pressure,” respectively. 2.2. Force Field and Molecular Dynamics Simulations Procedures. Molecular dynamics (MD) simulations were carried out for isolated single molecules in vacuum, as well as using periodic boundary conditions for simulation of the processes in the condensed phase. Two types of MD calculations were performed in this work: the “reaction” and the “diffusion” simulations. The reasoning for using these two calculation types is provided below in section 2.3. In the “reaction” MD simulations, the all-atom OPLS-AA force field14,15 was used for the parametrization of potential energy surfaces (PES). Either harmonic potential or (to enable C−C scission) Morse potential (see below) were employed for the description of the C−C stretch terms. MD calculations were performed using GROMACS program package, version 3.3.16 Each simulation consisted of three stages. At the first stage, potential energy minimization via the steepest descent method was performed. The second stage consisted of the simulation of the thermal motion of the molecule at a given temperature for the duration of 10−50 ps with randomized initial atomic velocities set according to the Maxwell− Boltzmann distribution. To avoid any chemical reaction during this stage, the Morse terms in the force field parametrization of the potential energy of interatomic interactions were substituted with harmonic potentials. The goal of the second stage was to create randomized starting molecular configurations for the third stage. At the third stage, the simulation of the thermal motion started during the second stage was continued at the same temperature. The transition between the stages involved transferring the coordinates and the velocities of all atoms obtained at the end of the second stage and using them as the initial configuration for the third stage. The difference in the dynamics during the second and the third stages was thus in the interaction potentials between carbon atoms in C−C bonds, which were changed from harmonic to Morse to enable C−C bond scission. The “diffusion” MD simulations were performed to model separation of the two fragments of a polymer chain resulting from a broken C−C bond, with the attractive Morse potential between the fragments removed. These simulations were conducted in three stages, the first two being identical to those of the “reaction” simulations described above. The third stage of the “diffusion” simulations differed from that of the “reaction” simulations by a modification of the force field, in which the strength of one C−C bond in the PE chain was

2. MODELING AND COMPUTATIONAL METHODS 2.1. Molecular Systems and Conditions. The reaction of C−C bond scission was studied for the linear polyethylene 2188

dx.doi.org/10.1021/jp411474u | J. Phys. Chem. A 2014, 118, 2187−2195

The Journal of Physical Chemistry A

Article

reduced to 0 kJ/mol (the “broken” bond). It is worth mentioning that the term “diffusion” is applied here for the lack of a better term: the process of separation of the molecular fragments connected by a “zero-strength” C−C bond most likely does not have a completely uniform Brownian character because the bending and the torsional terms of the OPLS-AA force field do not decrease as the distance between carbon atoms in the bond increases. This means that the motions of the produced chain fragments conserve some of the constraints from the time when the fragments constituted a single molecule, which does not allow their motions to be completely independent from each other. The rate of separation of fragments can be expected to be influenced by the position of the broken bond relative to the end of the molecular chain. Therefore, in the cases of butane (PE-2) and pentane the rates of separation of fragments produced due to the rupture of the middle and the terminal C− C bonds were determined separately. For the longer alkanes (PE-25, PE-100, and PE-500), only the bond in the middle of the chain was chosen as ruptured in such “diffusion” MD simulations, under the assumption of negligible influence of the terminal C−C bonds on the kinetics of the polymer decay. Potential effects of neglecting the effects of the broken bond position in PE-25 and longer chains are discussed below, in the end of section 2.3. Temperatures used in modeling ranged from 2300 to 3100 K; such high temperatures were necessary to make the MD reaction rates observable within the practicable range of times of molecular dynamics simulations. Although these temperatures are significantly higher than those of the experimental studies of polyethylene pyrolysis (typically 600−700 K), there is no alternative to using such high temperatures in molecular dynamics simulations; thus, temperatures in the 2000−3000 K range have been used virtually invariably in all earlier reactive molecular dynamics studies of polymer pyrolysis (e.g., refs 9−12 and 17). The temperature of the systems was controlled using the Berendsen thermostat algorithm18 with the coupling time 0.5 ps. Applicability of the Berendsen temperature control method to MD rate constant calculations for large molecular systems has been investigated and confirmed earlier in ref 12. The time integration step size of 0.25 fs was used in all simulations. A further decrease of the integration time step did not influence the determined values of the rate constants. Overall translational motion of the center of mass and rotation around the center of mass were removed in all calculations. Periodic boundary conditions with square periodic boxes were used in the simulations performed for the conditions of the condensed phase. 2.3. Determination of the Rate Constants. In any molecular dynamics determination of reaction rates, a criterion of reaction occurrence must be used to determine when a reactive event takes place. In earlier bond dissociation studies, separation of the two carbon atoms forming a bond by a critical distance (denoted as RC henceforth) was used as such criterion.9,11,12 The rate constant is then calculated as the reciprocal of the average time between the beginning of the reactive simulation (the “third stage” in our computational procedure) and the moment when any of the C−C bonds is stretched beyond the critical distance. Even though the choice of such a critical distance is rather arbitrary, in the gas-phase (vacuum) simulations the resultant rate constant values are not sensitive to the value of RC.12 However, in the condensed phase the situation is quite different. As described below, obtained

effective rate constants (keff) display noticeable dependences on the values of the critical distance RC. The upper part of Figure 1

Figure 1. Dependences of the per-bond C−C scission rate constants on the values of the critical distance RC (bond dissociation criterion) obtained in molecular dynamics simulations of decomposition of PE500 at 2900 K: open triangles, gas-phase conditions; filled circles, condensed phase with the density of 0.50 g cm−3; open circles, condensed phase with the density of 0.75 g cm−3. The Morse function plotted on the same scale of interatomic distances illustrates different behavior of the keff vs RC dependences in the regions dominated by strong attractive potential (short distances) and diffusional separation of fragments (long distances, weak potential).

demonstrates such keff(RC) dependences obtained in the simulations of PE-500 decomposition at T = 2900 K and three values of density (0.75 and 0.50 g cm−3 and vacuum); the dependence of the Morse potential on the bond distance is shown in the lower part of the figure. At short bond separations, where the attractive influence of the Morse potential is strong, all curves show strong dependences of keff on RC. At longer separations, the vacuum and the condensedphase results are qualitatively different. In the case of vacuum, the effective rate constant practically loses its dependence on the critical distance at separations beyond approximately 0.43 nm, which roughly corresponds to the Morse potential energy being within 1/4RT from its long-range asymptotic value. In the two cases corresponding to the condensed-phase simulations, however, the dependence of keff on RC continues well beyond this point. To illustrate the differences, linear regression fit lines are drawn through the five points of the long-range end of each of the dependences in Figure 1. For the condensed-phase simulations these linear fit lines have considerable nonzero slopes, unlike the case of the vacuum simulations. One can distinguish two distinct regions in the qualitative behavior of the keff vs RC curves in Figure 1. The left-hand 2189

dx.doi.org/10.1021/jp411474u | J. Phys. Chem. A 2014, 118, 2187−2195

The Journal of Physical Chemistry A

Article

region is potential-controlled whereas the limiting behavior of the right-hand region can be characterized as diffusioncontrolled. The behavior observed in the right-hand region for condensed-phase simulations can also be described as influenced by the cage effect, where the presence of the surrounding parts of molecular chains limits the rate of separation of the fragments connected initially by the dissociated C−C bond and thus promotes the reverse fragment recombination. Cage effects present a challenge to determination of the rate constant; although the two limiting potentialand diffusion-controlled regions are clearly distinguishable, the placement of the border between them is ambiguous and so is the associated value of the rate constant. In the current work, we present a method of “factoring out” the contributions of the cage effects and assigning meaningful rate constant values to condensed-phase carbon−carbon bond scission reactions in reactive molecular dynamics simulations. The only existing alternative is based on using an arbitrarily chosen critical separation distance, which introduces a significant degree of uncertainty in the resultant rate constant values. The current method uses the dependences of keff on RC such as those shown in Figure 1 and the results of “diffusional” simulations described in section 2.2. The process of bond dissociation and fragment separation is represented with the following kinetic model consisting of two reactions. R − R ⇄ [R···R] [R··· R] → 2R(at distance R C)

Figure 2. Dependence of the diffusion time tD on the square of the distance between carbon atoms in the “broken” bond (see text) at T = 2900 K, density of 0.75 g cm−3. Error bars are 2σ. Upper left inset: the corresponding dependence of the diffusional separation rate constant on RC. Lower right inset: a magnified part of the main plot at (RC)2 between 0.1 and 0.25 nm2; line extrapolated from the linear dependence at longer separations is shown to emphasize deviation from linearity at low values of RC.

(1,−1) (2)

The first reaction in the mechanism represents reversible dissociation of a linear alkane chain into two primary alkyl radicals that stay in the vicinity of each other (in a cage). The second reaction represents the process of diffusional separation of these radicals to the distance RC. This process of fragment separation is characterized by a rate constant k2 (denoted as kD henceforth for the clarity of its association with the diffusion process) that depends on the value of RC whereas the rate constants of reaction 1,−1 k1 and k−1 are independent of RC. The mechanism of reactions 1 and 2 is devised in an attempt to uncouple the effects of bond dissociation (reaction 1) and diffusional separation of fragments (reaction 2). Thus, k1 here is the rate constant of alkane dissociation that is not influenced by cage effects; the analysis seeks to determine its value. The process of diffusional separation of fragments (reaction 2) was modeled separately, as described in section 2.2. For each set of conditions, the values of kD(RC) were determined for a set of RC values as the reciprocal of the average time tD between the beginning of the third stage of the molecular dynamics simulation and the time when the fragments reach the distance RC. A typical dependence of kD on RC and the corresponding dependence of tD on the square of RC is shown in Figure 2. The latter dependence displays a linear character at long separations (0.25−0.65 nm2) but curves at shorter (less than 0.25 nm2) distances. The linearity of the long-distance part of the tD vs RC2 dependence indicates a certain degree of similarity between separation of alkane chain fragments in molecular dynamics simulations and the random motion of a Brownian particle. For a 3-dimentional Brownian motion, the dependence of the average first-passage time on the distance traveled r (i.e., the time it takes for a Brownian particle to move to distance r from its original position for the first time) is described by the equation

⟨t ⟩ =

r2 6D

(I)

where D is the diffusion coefficient.19 It is interesting to note that this dependence has the same functional form as that relating the average square of the displacement ⟨r2⟩ with time, t = ⟨r2⟩/6D. The nonzero intercept between the linear fit line in Figure 2 and the X-axis is expected because the initial separation between the chain fragments is not zero and, on average, equals the C−C bond length. Considering the firstpassage times t0 and t corresponding to two different distances r0 and r, from eq I one obtains ⟨t ⟩ − ⟨t0⟩ =

r2 r2 − 0 6D 6D

(II)

where the second term in the right-hand side of the equation is a constant. The deviation from linearity at short separations is likely to be caused by the influence of the angular and the torsional terms in the potential describing the interaction between the chain fragments; even though the Morse bonding term is removed in diffusional simulations, other terms remain, like in the potential used for reactive simulations. Using the mechanism of reactions 1,−1, and 2 and the steady-state approximation for [R···R], one obtains the following equation for the effective rate constant of the overall process of dissociation and separation of fragments:

keff = 2190

k1kD k −1 + k D

(III)

dx.doi.org/10.1021/jp411474u | J. Phys. Chem. A 2014, 118, 2187−2195

The Journal of Physical Chemistry A

Article

Here, k1 and k−1 are the forward and the reverse rate constants of reaction 1,−1 and kD is the diffusion rate constant of reaction 2. Equation III can be transformed to the following form: kD k k = −1 + D keff k1 k1

0.2−0.4 nm, with shorter separations at higher temperatures. This agrees with the conclusion that at RC = 0.5−0.8 nm the reaction is on the products side of the variational transition state and the process of further separation of fragments is diffusion-controlled. The values of k1 and k−1 for all conditions (temperature and PE density) used in RMD simulations were obtained from linear fits to the dependences of kD/keff on kD analogous to that demonstrated in Figure 3. Only the left-hand sides of such dependences corresponding to the diffusion-controlled range of RC fragment separations were used in the fits, for the reasons described above. The values of the rate constants obtained via the linear regression analysis are presented in Table 23S in the Supporting Information. As mentioned above, diffusion rate constants for PE-25 and longer chains were determined in MD simulations only for midchain bond scissions, although for PE-2 and pentane both the midchain and the terminal bond scission values of kD were obtained. For these short-chain alkanes, the ratios of terminal to nonterminal kD values averaged at 1.38 and did not depend on the RC value within the uncertainties in the RC = 0.5−0.8 nm range. Terminal bonds constitute 8% of all C−C bonds in PE25 and 2% in PE-100. Including their contributions to the average diffusion rate constants may slightly increase the values of the latter; however, the factor of this increase is likely to be independent of RC, based on analogy with PE-2 and pentane. Such an increase in kD will affect the RC-dependent terms on both sides of eq IV in the same way and will have no effect on the slope of the line of the kD/keff vs kD dependence; thus, the derived values of k1 will not be affected.

(IV)

Because kD and keff depend on the critical distance RC and k1 and k−1 do not, eq IV facilitates analysis of the data on keff and kD obtained in molecular dynamics simulations using different values of RC. When these data are plotted in the kD/keff vs kD coordinates, one can expect a linear dependence with a slope equal to the reciprocal of k1 and a Y-axis intercept equal to the reciprocal of the equilibrium constant for the reaction 1,−1. An example of a dependence of kD/keff on kD is presented in Figure 3. As can be seen from the plot, the kD/keff vs kD

3. RESULTS The results of molecular dynamics simulations of the reaction of C−C bond dissociation in PE-500 are presented in Figure 4 and Table 23S (in the Supporting Information). As can be seen from the Arrhenius plots, the reaction of C−C bond dissociation in the condensed phase is noticeably slower than that in the gas phase, with the differences of up to an order of magnitude. The condensed-phase rate constants also demonstrate a strong dependence on the density, with much of the difference between the gas-phase rate constants and those obtained for the “normal” density of 0.75 g cm−3 recovered if simulations are performed with the reduced density of 0.5 g cm−3. It is worth emphasizing that the observed slowing of the scission reaction with increasing density is not a consequence of the cage effect because the latter has been accounted for in the kinetic analysis described above. The Arrhenius parameters for the C−C backbone scission rate constant for PE-500 are given in Table 1. One can see that despite the differences in the central values of the obtained activation energies, their confidence intervals (determined as two standard errors of the mean) overlap. Because the value of RT for the temperature interval of the current study is about 22 kJ mol−1, one can assert that the obtained activation energies differ by roughly 2RT from 348 kJ mol−1, the value of the threshold of the Morse potential used in the force field to describe the C−C bond, which is similar to a typical activation energy for C−C scission for small alkane molecules in the gas phase.22−24 The results of the study of C−C scission in PE molecules of varying length are presented in Table 2 and Figure 5a. In the figure, open diamonds represent the results obtained in simulations performed under the conditions of constant density

Figure 3. Dependence of the kD/keff ratio on the RC-dependent diffusion constant kD for T = 3100 K, density of 0.75 g cm−3. The fitted part of the dependence corresponds to longer separations of the C−C bond (0.5 to 0.8 nm), where the influence of the C−C interaction potential on the motions of the radicals is negligible. Error bars are 1σ.

dependence has a linear character only in the region of low values of kD, which corresponds to larger separations of the radical fragments (RC = 0.5−0.8 nm). This is the same range of RC values where the dependences of keff on RC are controlled by diffusion (as opposed to the region controlled by the bonding potential, Figure 1) and the diffusional separation of the radical fragments demonstrates a Brownian behavior (i.e., a linear dependence in the coordinates of Figure 2). We interpret this combination of data in Figures 1−3 as indicating that the twostep kinetic mechanism of reactions 1,−1, and 2 and eqs III and IV are applicable in the region of RC = 0.5−0.8 nm, whereas at shorter fragment separations the values of k1 and k−1 are not truly independent of RC. The effects of RC on k1 and k−1 at short distances between fragments are likely to be due to the influence of the bonding potential; these effects result in the nonlinearity of the kD/keff vs kD dependences (the right-hand side of Figure 3). In variational transition-state theory calculations utilizing quantum chemistry-based potential energy surfaces (e.g., refs 20 and 21), the transition state for homologous C−C bond scission is usually located at RC = 2191

dx.doi.org/10.1021/jp411474u | J. Phys. Chem. A 2014, 118, 2187−2195

The Journal of Physical Chemistry A

Article

Figure 4. Arrhenius plots of the per-bond C−C dissociation rate constants k1 obtained in molecular dynamics simulation of PE-500 decomposition. Lines are Arrhenius fits through the data. Filled squares and solid line: gas-phase conditions. Open circles and dashed line: condensed phase with the density of 0.5 g cm−3. Open triangles and dotted line: condensed phase with the density of 0.75 g cm−3. Error bars are 2σ.

Figure 5. Values of C−C scission rate constant obtained in molecular dynamics simulations under the constant density conditions with the density kept at 0.75 g cm−3 (diamonds) and under the constant pressure conditions (pressure kept at 1.8 × 104 bar, circles). Both sets of simulations were performed at 3100 K. Error bars are 2σ. (a) Rate constants obtained using the kinetic mechanism described in section 2.3. (b) Rate constants obtained using the value of the critical separation RC = 0.45 nm as the criterion of reaction occurrence. This RC value corresponds to the energy of the C−C bond stretch equal to 99% of its dissociation limit, as in the study of Smith et al.;11 see section 4.

Table 1. Arrhenius Parametersa for the Reaction of C−C Backbone Scission in PE-500 Determined for Three Values of Density

a

density/g cm−3

A/s−1

Eact/kJ mol−1

0.75 (normal) 0.5 (reduced) 0 (vacuum)

(4.21 ± 0.98) × 1012 (1.31 ± 0.68) × 1014 (2.06 ± 0.46) × 1014

275 ± 28 326 ± 60 323 ± 26

As can be seen in the figure, the rate constant obtained under the constant density conditions decrease by about an order of magnitude as the length of the polymer chain decreases from 500 monomers to one monomer. At the same time, the rate constants determined under the constant pressure conditions do not show any statistically meaningful dependence on the length of the chain.

Uncertainties are 2σ.

Table 2. Values of the C−C Bond Decay Rate Constant Obtained for the Systems of Linear Alkanes of Varying Chain Lengths under the Conditions of Constant Density and Constant Pressure at 3100 K log(k1/s

a

molecular system

constant density

500PE-1 250PE-2 200C-5 20PE-25 5PE-100 PE-500

6.91 ± 0.24 7.19 ± 0.06 8.06 ± 0.21 7.90 ± 0.12 8.00 ± 0.08

4. DISCUSSION An earlier systematic molecular dynamics study of the C−C bond scission in vacuum for a series of linear alkanes and polyethylene macromolecules by Knyazev12 demonstrated that increasing chain length results in significant acceleration of the reaction, with an increase of up to an order of magnitude in the values of the per-bond rate constant, kCC. This effect was tentatively attributed to centrifugal effects of torsional motions about the C−C bonds in the alkane chains. The differences between the gas-phase and the condensed-phase values of kCC (per-bond k1 of the mechanism of reactions 1,−1, and 2 in section 2.3) observed in the current study are similar in value to those found in ref 12 between the smallest molecule considered (C2H6) and long-chain alkanes; both are approximately a factor of 10. It appears tempting to surmise that the effect of reaction deceleration in the condensed phase is caused by the obstruction of the torsional degrees of freedom of the alkane

−1a

)

constant pressure 7.79 8.02 8.15 8.05

± ± ± ±

0.02 0.27 0.28 0.28

8.00 ± 0.08

Uncertainties are 2σ.

(0.75 g cm−3). Filled circles represent the data obtained in simulations where densities of the PE with different chain lengths were varied to maintain the same pressure, 1.8 × 104 bar, the pressure of PE-500 at the “normal” density of 0.75 g cm−3 at 3100 K. 2192

dx.doi.org/10.1021/jp411474u | J. Phys. Chem. A 2014, 118, 2187−2195

The Journal of Physical Chemistry A

Article

torsional degrees of freedom and the results of Smith et al. The results obtained in the current work for the constant density conditions are in qualitative agreement with those of Smith et al.: the per-bond rate constant increases in the sequence PE-1, PE-2, PE-25 and tends to saturation for longer chains (Figure 5a). At the same time, the rate constant values obtained under the constant pressure conditions do not show any discernible dependence on the chain length, as can be seen from the same graph. We therefore conjecture that the influence of the chain length on the C−C scission rate constant is caused by the differences in the degree of molecular crowding, or the availability of intermolecular space for molecular motion, including the motion required for extension of the C−C bonds during scission. The concept of molecular crowding is intuitively clear but difficult to define quantitatively. This term refers to how much of the physical space is occupied by molecules and thus how much (or rather how little) of the remaining space is available to molecular motion. It is different from the cage effect, i.e., limitations on the rate of separation of dissociated fragments due to obstruction of molecular motion by surrounding molecules, yet the cage effect is expected to be affected by the degree of molecular crowding. One way to quantify molecular crowding can be based on the pressure in the system: a higher degree of crowding would correspond to higher pressure. Figure 7 illustrates the dependence of the

chains. Thus, within the framework of this hypothesis, reaction acceleration observed in the gas phase with the chain length increasing from PE-1 (C2H6) to PE-25 and longer chains (up to PE-1000) should not occur in the condensed phase of polymer melt. At a first glance, this conclusion appears to contradict the findings of Smith et al.,11 who observed significant differences between the per-bond C−C scission rates in short and long PE molecules in their condensed-phase RMD study. In their work, the reaction accelerated significantly (up to 2 orders of magnitude) as the length of a PE chain increased from one to one hundred monomers. The PE molecules used in the simulations of Smith et al. had the chemical formulas C2nH4n, where n is the number of monomers in the chain. The difference with the normal alkane formula C2nH4n+2 is due to the use of the CH2 end-groups with sp3 hybridized carbon atoms. Thus, the smallest PE-1 molecule of Smith et al. had the chemical formula C2H4, with a single (not double) C−C bond and sp3 carbons. All calculations were performed with the density of PE melt maintained between 0.80 and 0.83 g cm−3. The temperature dependences of the C−C scission rate constant obtained in ref 11 for PE-1 and PE-100 are plotted in Figure 6 together with the results obtained in the current work

Figure 7. Values of pressure obtained in the simulations with the density kept equal to 0.75 g cm−3 in all the systems of linear alkanes studied in the current work.

pressure obtained in RMD simulations performed in this work for the constant density conditions on the polyethylene chain length. The pressure increases by approximately a factor of 2 as the number of monomers in the chain decreases from 500 to 1, demonstrating the corresponding increase in molecular crowding. The reason for this pressure/crowding increase is in the different dependences of the bonding and nonbonding potentials on interatomic distances. Figure 8 illustrates such dependences for carbon−carbon interactions, for the force fields used in the current work as well as by Smith et al. Nonbonding potentials in general result in larger interatomic distances, particularly at the translational energies in the thermal range. Replacing a long-chain molecule with a set of shorter ones while keeping the same number of carbon atoms results thus in higher degree of molecular crowding due to larger effective space occupied by these molecules. Another contribution to the higher degree of molecular crowding in the

Figure 6. Results obtained in the current work for PE-500 in vacuum (circles) and at the density of 0.75 g cm−3 (squares) plotted together with the results obtained by Smith et al.11 for PE-1 fragments (heavy dashed line) and PE-100 molecule (heavy solid line). Heavy dotted line: experimental k(T) dependence for ethane decomposition from Kiefer et al.25

for PE-500. The difference between the PE-1 and the PE-100 rate constants reach 2 orders of magnitude between 2500 and 2840 K, where the temperature ranges used by Smith et al. for these two chain lengths overlap. An experimental k(T) dependence reported by Kiefer et al.25 for the thermal decomposition of ethane is shown for comparison. The main motivation for our study of the dependence of the per-bond C−C scission rate constant on the polyethylene chain length was to elucidate the apparent disagreement between the hypothesis of the effects of the condensed phase on the 2193

dx.doi.org/10.1021/jp411474u | J. Phys. Chem. A 2014, 118, 2187−2195

The Journal of Physical Chemistry A

Article

separation criteria (critical distances), such as those illustrated in Figure 1. In the study of Smith et al.,11 a critical distance of 0.45 nm (corresponding to the bond stretch beyond 99% of the dissociation energy) was used as the criterion for the occurrence of the scission reaction. Although the values of the rate constants obtained using the two methods necessarily differ, the important qualitative conclusions of the current work are not method-specific. Figure 1 demonstrates that even if a fixed value of the critical C−C distance is used as the criterion of bond dissociation, the effective rate constants corresponding to the condensed-phase conditions are approximately an order of magnitude lower than those obtained in vacuum. Also, Figure 5b shows the different dependences of the per-bond C− C scission rate constants on the number of monomers in the polyethylene chain obtained under the constant density and constant pressure conditions using a fixed critical distance of 0.45 nm. As can be seen from the plots and from comparison with Figure 5a, the effect attributed to molecular crowding persists regardless of the method used to calculate rate constant values. The results of the current study do not provide for an unequivocal support or rejection of the hypothesis12 that acceleration of C−C scission in the gas phase with increasing molecular chain length is caused by centrifugal effects of torsional motions. Indeed, the differences between the gasphase and the condensed-phase values of per-bond rate constant observed in the current study for PE-500 are of the same value (an order of magnitude) as the differences between the gas-phase decomposition rates of the smallest PE molecule (PE-1, C2H6) and long-chain alkanes. This seems to support the above hypothesis because torsional degrees of freedom of long alkane chains can be expected to be obstructed in the condensed phase. However, as Figure 5 and the data in Table 2 demonstrate, in the condensed phase increase in pressure results in higher degree of molecular crowding and decreasing values of the rate constant, which is not explained by the above hypothesis. Because transition from the gas phase to the condensed phase is also accompanied by increase in pressure, the two effects may or may not be related and may be caused by the same or different phenomena. Finally, the results of the current study provide a positive answer to the question whether Arrhenius expressions for cognate gas-phase reactions of small molecules can serve as an adequate approximation to the rate constants of the C−C scission reactions in long-chain polyethylene molecules under the conditions of polymer melt. The condensed-phase rate constants are an order of magnitude lower than the gas-phase constants for long polyethylene chains. Yet the latter are an order of magnitude higher than the gas-phase constants for small molecules such as C2H6.12 Therefore, in the absence of any significant differences in activation energies, as confirmed by the current study, per-bond rate constants of C−C scission of long alkanes in the condensed phase is likely to be similar to those of small alkanes in the gas phase.

Figure 8. Bonding and nonbonding potential terms for C−C interaction used in the current study (solid lines) and by Smith et al.11 (dashed lines). Nonbonding potentials result in larger effective distances between carbon atoms. The horizontal dotted line corresponds to 1/2RT at 3100 K.

short-molecule systems used in the current work is due to the presence of additional hydrogen atoms. It is thus likely that the differences between the rates of C−C scission obtained for a long-chain (PE-100) molecule and a collection of short (PE-1) molecules observed by Smith et al.11 are caused, at least partially, by the differences in molecular crowding. As a chain consisting of one hundred monomers is replaced by one hundred of C2H4 fragments, the bonding interaction potential in every other C−C bond in the system effectively gets replaced by a nonbonding van der Waals potential, with larger effective atomic sizes. This leaves less free intermolecular space available for molecular motion. As the results of the current work are compared with those of Smith et al., it is necessary to highlight some of the important differences in the methodologies used in the two studies. The force field used by the authors of ref 11 differs from that used in the current work. One important difference is in the more realistic description of the angular and the torsional potential terms applied in ref 11, enabled by the software developed by these authors earlier.26 The force constants of these terms gradually diminish as the adjacent C−C bond is stretched, thus resulting in more “loose”27 transition states for C−C scission and in larger values of the rate constants (Figure 6). These rate constants are, on average, in better agreement with the experimental gas-phase values obtained for the thermal decomposition of ethane by Kiefer et al.,25 as illustrated in the same figure. The publicly available OPLS-AA force field used in the current study does not allow for such a decrease in the bending and the torsional force constants as the carbon atoms in the bond separate. This difference in the force fields may also have an influence on the rate with which the mutual orientation of the fragments of the broken bond relaxes during the separation of the carbon atoms forming the bond. Another important difference is in the method of determining the values of the rate constants. In the current study, a method based on the kinetic mechanism of reactions 1,−1, and 2 was developed (section 2.3, eqs III and IV, Figure 3). The method takes into account cage effects specific to the condensed phase of the polymer melt, which in this case means the influence of the finite rates of diffusional separation of radicals on the observed effective rate constants (Figure 1). The values of rate constants obtained after “factoring out” the cage effects are not affected by the influences of the choice of the

■

ASSOCIATED CONTENT

S Supporting Information *

Tables 1S-23S of detailed results of molecular dynamics calculations: keff, tD, and per-bond C-C scission rate constant values. This material is available free of charge via the Internet at http://pubs.acs.org. 2194

dx.doi.org/10.1021/jp411474u | J. Phys. Chem. A 2014, 118, 2187−2195

The Journal of Physical Chemistry A

■

Article

(16) Van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, Flexible, and Free. J. Comput. Chem. 2005, 26 (16), 1701−1718. (17) Chenoweth, K.; van Duin, A. C. T.; Goddard, W. A. ReaxFF Reactive Force Field for Molecular Dynamics Simulations of Hydrocarbon Oxidation. J. Phys. Chem. A 2008, 112, 1040−1053. (18) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81 (8), 3684−3690. (19) Stadie, F. Probleme der Brownschen Molecularbewegung. Ann. Phys. 1928, 391 (13), 751−797. (20) Truhlar, D. G.; Garrett, B. C.; Klippenstein, S. E. Current Status of Transition-State Theory. J. Phys. Chem. 1996, 100, 12771−12800. (21) Klippenstein, S. E.; Georgievskii, Yu.; Harding, L. B. Predictive theory for the combination kinetics of two alkyl radicals. Phys. Chem. Chem. Phys. 2006, 8, 1133−1147. (22) Baulch, D. L.; Cobos, C. J.; Cox, R. A.; Esser, C.; Frank, P.; Just, T.; Kerr, J. A.; Pilling, M. J.; Troe, J.; Walker, R. W.; Warnatz, J. Evaluated Kinetic Data for Combustion Modeling. J. Phys. Chem. Ref. Data 1992, 21 (3), 411−734. (23) Tsang, W.; Hampson, R. F. Chemical Kinetic Data Base for Combustion Chemistry. I. Methane and Related Compounds. J. Phys. Chem. Ref. Data 1986, 15 (3), 1087−1279. (24) Warnatz, J. Combustion Chemistry; Springer-Verlag: New York, 1984. (25) Kiefer, J. H.; Santhanam, S.; Srinivasan, N. K.; Tranter, R. S.; Klippenstein, S. J.; Oehlschlaeger, M. A. Dissociation, Relaxation and Incubation in the High-Temperature Pyrolysis RRKM of Ethane, and a Successful Modeling. Proc. Combust. Inst. 2005, 30, 1129−1135. (26) Smith, K. D.; Stoliarov, S. I.; Nyden, M. R.; Westmoreland, P. R. RMDff: a smoothly transitioning, forcefield-based representation of kinetics for reactive molecular dynamics simulations. Mol. Simul. 2007, 33, 361−368. (27) Benson, S. W. Thermochemical Kinetics; Wiley-Interscience: New York, NY, 1976.

AUTHOR INFORMATION

Corresponding Author

*V. D. Knyazev: e-mail, [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported, in part, by National Institute of Standards and Technology, Fire Research Grants Program, via Grant No. 70NANB4H1128 and, in part, by the Ministry of Science and Education of the Russian Federation under Agreement 8532 of March 7, 2012. The authors thank Drs. K. D. Smith, S. I. Stoliarov, M. R. Nyden, and P. R. Westmoreland for useful discussions and sharing of data.

■

REFERENCES

(1) Poutsma, M. L. Reexamination of The Pyrolysis of Polyethylene: Data Needs, Free-Radical Mechanistic Considerations, and Thermochemical Kinetic Simulation of Initial Product-Forming Pathways. Macromolecules 2003, 36, 8931−8957. (2) Ranzi, E.; Dente, M.; Faravelli, T.; Bozzano, G.; Fabini, S.; Nava, R.; Cozzani, V.; Tognotti, L. Kinetic Modeling of Polyethylene and Polypropylene Thermal Degradation. J. Anal. Appl. Pyrol. 1997, 40− 41, 305−319. (3) Faravelli, T.; Bozzano, G.; Scassa, C.; Perego, M.; Fabini, S.; Ranzi, E.; Dente, M. Gas Product Distribution from Polyethylene Pyrolysis. J. Anal. Appl. Pyrol. 1999, 52, 87−103. (4) Nemeth, A.; Blazso, M.; Baranyai, P.; Vidoczy, T. Thermal Degradation of Polyethylene Modeled on Tetracontane. J. Anal. Appl. Pyrol. 2008, 81, 237−242. (5) Levine, S. I.; Broadbelt, L. J. Detailed Mechanistic Modeling of High-Density Polyethylene Pyrolysis: Low Molecular Weight Product Evolution. Polym. Degrad. Stab. 2009, 94, 810−822. (6) Jalan, A.; Ashcraft, R. W.; West, R. H.; Green, W. H. Predicting solvation energies for kinetic modeling. Annu. Rep. Prog. Chem., Sect. C 2010, 106, 211−258. (7) Madorsky, S. L. Thermal Degradation of Organic Polymers; Interscience Publishers: Washington, DC, 1964. (8) Denisov, E. T.; Sarkisov, O. M.; Likhtenshtein, G. I. Chemical Kinetics: Fundamentals and New Developments; Elsevier: Amsterdam, 2003. (9) Nyden, M. R.; Stoliarov, S. I.; Westmoreland, P. R.; Guo, Z. X.; Jee, C. Applications of Reactive Molecular Dynamics to the Study of the Thermal Decomposition of Polymers and Nanoscale Structures. Mater. Sci. Eng. 2004, A365, 114−121. (10) Stoliarov, S. I.; Lyon, R. E.; Nyden, M. R. A Reactive Molecular Dynamics Model of Thermal Decomposition in Polymers. II. Polyisobutylene. Polymer 2004, 45, 8613−8621. (11) Smith, K. D.; Bruns, M.; Stolarov, S. I.; Nyden, M. R.; Ezekoye, O. A.; Westmoreland, P. R. Assessing the Effect of Molecular Weight on the Kinetics of Backbone Scission Reactions in Polyethylene Using Reactive Molecular Dynamics. Polymer 2011, 52, 3104−3111. (12) Knyazev, V. D. Effects of Chain Length on the Rates of C-C Bond Dissociation in Linear Alkanes and Polyethylene. J. Phys. Chem. A 2007, 111 (19), 3875−3883. (13) Olabisi, O.; Simha, R. A Semiempirical Equation of State for Polymer Melts. J. Appl. Polym. Sci. 1977, 21 (1), 149−163. (14) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118 (45), 11225−11236. (15) Price, M. L. P.; Ostrovsky, D.; Jorgensen, W. L. Gas-Phase and Liquid-State Properties of Esters, Nitriles, and Nitro Compounds with the OPLS-AA Force Field. J. Comput. Chem. 2001, 22 (13), 1340− 1352. 2195

dx.doi.org/10.1021/jp411474u | J. Phys. Chem. A 2014, 118, 2187−2195