Molecular and Kinetic Models for High-Rate Thermal Degradation of

Subscriber access provided by MT ROYAL COLLEGE

Article

Molecular and Kinetic Models for HighRate Thermal Degradation of Polyethylene J. Matthew D. Lane, and Nathan W Moore J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.7b11180 • Publication Date (Web): 01 Feb 2018 Downloaded from http://pubs.acs.org on February 1, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry A is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Molecular and Kinetic Models for High-Rate Thermal Degradation of Polyethylene J. Matthew D. Lane∗ and Nathan W. Moore Sandia National Laboratories, Albuquerque, NM 87185 E-mail: [email protected]

∗

To whom correspondence should be addressed

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract Thermal degradation of polyethylene is studied under the extremely high-rate temperature ramps expected in laser-driven and x-ray ablation experiments – from 1010 to 1014 K/s in isochoric, condensed phases. The molecular evolution and macroscopic state variables are extracted as a function of density from reactive molecular dynamics simulations using the ReaxFF potential. The enthalpy, dissociation onset temperature, bond evolution and observed cross-linking are shown to be rate dependent. These results are used to parameterize a kinetic rate model for the decomposition and coalescence of hydrocarbons as a function of temperature, temperature ramp rate, and density. The results are contrasted to first-order random-scission macrokinetic models often assumed for pyrolysis of linear polyethylene under ambient conditions.

Introduction Recent extreme environments experiments at both Sandia and Livermore National Laboratories challenge our understanding of the high-rate degradation of polymer materials. These state-of-the-art experiments involve high-rate energy deposition, which drives both a thermal ramp and a mechanical shock wave, resulting in an evolving mixture of nonequilibrium solid and vapor phases. At LLNL’s National Ignition Facility (NIF), inertial confinement fusion experiments are using precision-manufactured ablative polymer coatings to optimize implosion sphericity during fuel capsule compression. 1 During compression, heating rates exceed 1014 K/s in the hydrocarbon coated fuel capsules. 2 At Sandia, recent experiments on the Z Machine have heated polymeric materials, using ∼3 keV x-rays, to 100+ µm depths over tens of ns rise times, thereby driving temperatures to thousands of Kelvin at heating rates of 1012 K/s. By contrast, widely available polymer degradation data are generally measured from thermogravimetric decomposition experiments, in which thermal ramps rarely exceed 1 K/s. 3,4 Here, we apply molecular dynamics simulation to investigate the role of thermal ramp 2


Page 2 of 28

Page 3 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


rate and density variation in the degradation processes of polyethylene. We then use these molecular simulations to construct a kinetic model which incorporates these processes. More generally, understanding the response of polymeric materials in ablation conditions, which span solid to vapor, can also aid in extending our equations-of-state (EOS) models to account for response due to physical and chemical processes which are coupled to dynamics. Molecular dynamics (MD) simulations may present an avenue for developing and/or evaluating EOS over a large phase space. Polyethylene (PE) is a useful test case because of its structural simplicity. Beyler and Hirschler, 5 in their review of thermal decomposition of polymers, identify four primary mechanisms for thermal degradation of polymers: (1) random-chain scission, (2) end-chain scission, (3) chain stripping, and (4) cross-linking. In linear polyethylene, chain stripping is not a process, since there are no side-chains or branching structures to cleave. Most decomposition models, which are based on experiments involving rates reminiscent of combustion or conventional heating processes, focus on low-rates and emphasize bond scission processes and associated volatilization of low-mass end-products. The pyrolysis of linear alkanes has been long-studied, 5–7 with recent molecular-level investigations highlighting, for example, the roles of chain structure and length, 8–10 chain length distributions, 11 evaporation rates, 12 density, 10,13,14 and hydrogen abstraction 15 on the overall conversion rate of solid polyethylene to gaseous products using a wide range of experimental and computational tools. We show here that, at extremely high thermal-ramp rates, cross-linking effects should not be ignored in kinetic models and energy balance.

Results and Discussion Molecular Dynamics Models Molecular dynamics has been used to study chemical decomposition of polymers. Nyden et al. 16 used single-chain reactive molecular dynamics simulation to study scission products for 3



Figure 1: Pressure dependence on temperature for 1013 K/s ramps in polyethylene samples of density ranging from 0.3 to 0.96 g/cm3 (full density). several polymers. Stoliarov et al. 17 used similar methodology to study PMMA. Popov and Knyazev also used single-molecule MD to study density effects on the effective bond-energy for first scission events. 14 Knyazev studied the effect of chain length, as well, showing that chain-length effects saturate quickly for PE chains of over 5 CH2 monomer units. 8 We extend these past studies to explicitly study the effect of energy-deposition rate on the degradation, and observe cross-linking and molecular recombination in addition to scission processes. These non-scission events are found to contribute significantly to the energy balance of the degradation. Further, we propose a kinetic model based on our molecular models which generalizes the effect of rate. Past kinetic models, 18 which have focussed on bond scission, appear to be insufficient to capture degradation at the extremely high rates encountered in laser-driven and pulsed x-ray heating experiments. Large-scale classical molecular dynamics, using Sandia’s LAMMPS code, 19,20 were conducted on systems of 168 linear polyethylene chains of length 44 monomers. The simple linear polyethylene structure readily forms semi-crystalline solids. Initial full-density systems were constructed in a hexagonal lattice with dimensions of 5.2 nm × 6.0 nm × 6.0 nm, and 22,176 atoms. From this 0.93 g/cm3 initial density a range of densities were equilibrated from 0.3 g/cm3 to 0.96 g/cm3 . These relatively modest system sizes were chosen to allow 4


Page 4 of 28

Page 5 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


extremely long simulations, which span 4 decades of thermal ramp rate from 1010 to 1014 K/s. A small number of systems with twice the number of PE chains were run to verify that results were not system size dependent. An all-atom representation was used, with the reactive ReaxFF interatomic potential. 21,22 ReaxFF is a bond-order force field which allows chemical reactions (through variable charge equilibration), and dynamic bond formation and breaking. ReaxFF has been used successfully to model polymer chemistry, and specifically has been shown to capture response in extreme environments. For example, one of us has shown, previously, that ReaxFF quantitatively agrees with DFT and experiment in extreme environments of polyethylene up to 50 GPa and greater than 10,000 K in shock compression simulations. 23–25 Thermal ramps were applied with a linear temperature profile from 300 to 5000 K at rates of 1010 to 1014 K/s. Constant pressure and temperature (NPT) simulations were used to equilibrate the initial systems. Then constant volume simulations with constrained temperatures (NVT) were used during the thermal ramps. Constant volume simulations are used to approximate the inertially confined state of experiments at these extreme rates. A Nose-Hoover thermostat was used with a 0.1 fs timestep. Short timesteps were necessary due to the extreme temperature conditions and the high-frequency oscillations of the carbonhydrogen bond. Total simulation times ranged from 7 ps to 12.75 ns, or 70 000 to 127 500 000 timesteps, respectively. Following the definition in the ReaxFF potential, bonding is defined by instantaneous atom separation distances. One of the strengths of classical molecular dynamics is the ability to efficiently model lowdensity states present in ablating polyethylene vapor. In these states, chemical dissociation plays a significant role in the energy balance, and therefore in the macroscopic relationships between pressure, temperature, and density (i.e. the EOS). Figure 1 shows the pressure dependence on temperature for a single ramp rate of 1013 K/s and for several densities ranging from full density (0.96 g/cm3 ) down to 0.3 g/cm3 . At all densities, we see a characteristic change in P -T slope between 2600 and 2800 K, which indicates the onset of chemical decom-

5



Figure 2: Pressure-density isotherms predicted by ReaxFF for polyethylene are smooth except in the low-temperature and high-density regime where solids undergo a crystal/amorphous transition. position. At the higher densities, we also see, near 700 K, the signature of an order-disorder transition of the initially aligned molecular crystal. This melt transition is observable for several densities, but the effect is less dramatic for lower densities and disappears entirely below 0.7 g/cm3 . The same (ρ, P, T ) data is plotted, in Figure 2, as isotherms in pressure versus density space. At 1000 K and above, the isotherms show smooth trends, while below 1000 K, we see increased fluctuation in the data, especially at densities above 0.7 g/cm3 . This observation is consistent with the order/disorder (i.e. melt) transition seen in the P -T plot, as well as in structure visualizations of these systems. These data for 1013 K/s indicate that the macroscopic state is capturing the effect of the polyethylene chemistry, but point to a dissociation onset temperature much higher than typically reported from experiments. Low-rate experiments typically measure dissociation onset below 1000 K. 3 In order to explore the dependence on thermal ramp rate, we conducted simulations with rates ranging from 1014 down to 1011 K/s. Figure 3 shows a plot of the system enthalpy as a function of temperature for our two extremum densities, 0.96 and 0.3 g/cm3 . The dissociation process in these plots is seen by a shift in the otherwise linear slope. The higher slope region, which indicates chemical activity, shifts to lower tempera-

6


Page 6 of 28

Page 7 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 3: Enthalpy as a function of temperature for ramp rates from 1011 to 1014 K/s for 0.3 g/cm3 and 0.96 g/cm3 initial density. Enthalpy shows two increases in 0.96 g/cm3 samples, indicating both the solid disordering transition near 600 K, and the dissociation transition above 2000 K. In the 0.3 g/cm3 sample, only the dissociation transition is seen. Degradation threshold temperature is highly dependent on temperature ramp rate. tures for slower rates, and also seems to slightly steepen, indicating a narrower dissociation temperature range. This trend would be consistent with the observation that slower rate experiments have lower dissociation onset temperatures. Figure 4 plots the dissociation temperature range as a function of thermal ramp rate, along with experimental data taken with rates nearly 8 decades lower. When plotted on log-log scales, this data is connected with a surprisingly simple power law relationship across these widely different rates. The fit in the plot extrapolates from the low rate experiments, and can be seen to pass through the high-rate simulations. The power law exponent is ∼ 0.04. The fit agreement appears to be weakest at the 1014 K/s rate, and it should be noted that a downward extrapolation from the simulations would not intersect the low-rate data, perhaps indicating new physics is at play at the very highest-rate simulations. It should be noted that while we have shown that the onset temperature of the degradation process is rate dependent, this does not require or imply that the activation energy is rate dependent. In fact, our molecular model and our kinetic analysis (see below) explicitly assume that activation energies are not affected by energy deposition rate. The shift in the 7



Figure 4: The initial dissociation threshold range vs. thermal ramp rate, plotted for MD simulations (blue) and previously published experiments 3 (red) of full density polyethylene. A simple power law fit connects the data across approximately 8 decades of ramp rate. onset temperature is likely the result of the relative time scales of chemical processes to the applied ramp rate, not to any change in the energy required to break bonds. While confirmation of macroscopic state is important, the molecular simulations allow us to probe more deeply into the processes which determine this state. Specifically, we have the ability to answer questions about the details of chemical dissociation, such as the role of random vs. end-chain scission. At lower rates, the nature of the scission processes helps determine the mechanical properties of the remaining material, and determines the volatile species which are outgassed. Kinetic models have been proposed which predict the chain length distributions given by different processes. 8–11 At high-rates the distinction between random and end-chain scission is not the dominant characteristic of chemistry. And, bond scission-based macrokinetic model predictions of chain length distributions do not compare well with our MD results. A more detailed investigation of the polyethylene bonding shows that a random bond scission model cannot account for the degradation processes at high ramp rates. The solid lines plotted in Figure 5 look very much like the trend lines predicted from scission processes, in which the initial system’s carbon-carbon (CC) bonds (solid) drop to

8


Page 8 of 28

Page 9 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5: Plot indicating the number of systemwide carbon-carbon bonds as a function of temperature during 1011 K/s thermal ramp in polyethylene. Initially scission dominates, reducing the number of bonds, but recombination ultimately increases the total number of bonds in the system. The initial bonds drop to zero indicating that all the originally chain bonds have been broken, while others have formed. zero over a relatively narrow temperature range. However, a closer look shows that, while the initial bonds (i.e. the initial chain backbone bonds) in the system are dissociated, new bonds are formed, such that the total number of bonds in the system drops by less than 50%, before recovering. This plot points to the importance of recombination and crosslinking. At high temperatures, the degradation of polyethylene can be expected to lead to char formation, especially at low thermal ramp rates. In that case, gases are able to diffuse to a surface and then out of the material, leaving a condensation of dense aromatic structures. 27 Our bulk simulations do not have a free surface, and therefore do not allow for removal of volatiles (e.g. H2 ). A variety of metrics for characterizing the progression of organic chars have been developed, typically related to the degree of aromaticity (total proportion of aromatic C) and/or the degree of aromatic condensation (proportion of aromatic C in condensed phases). 28 Aromatic formation is difficult to observe directly in our isochoric simulations, due to the short timescales. However, since many indirect metrics for aromaticity are measured through spectroscopic sensitivities, we take as a simpler metric for the extent of char formation, the ratio, χHC = (nH − nH2 )/nC , which is analogous to the H/C ratio 9



Figure 6: Plot indicating the degree of dehydrogenation of the backbone chain, as measured by the ratio χHC = (nH − nH2 )/nC as a function of temperature during thermal ramps at two rates in polyethylene. used in experiments. nH is the total number of hydrogen, nH2 is the molecular hydrogen and nC is the total carbon. In pyrolysis experiments, where outgassing from the reaction vessel keeps nH2 ∼ 0, a decrease in χHC indicates the formation of unsaturated carbons such as aromatic rings. 29 Here, the expectation is χHC ∼ 2 for aliphatic hydrogen, χHC ∼ 1 for aromatic hydrogen and χHC ∼ 0 for condensed phases. This degree of dehydrogenation of the carbon backbones might indicate the longer-time formation of ring structures. Figure 6 shows the evolution of the ratio of bound hydrogen to carbon χHC through ramps at 1011 and 1012 K/s. At low temperatures, the value is near 2, indicating that there are 2 hydrogen bound to each carbon in the system. As the temperature rises, the ratio drops as the chains dehydrogenate, and radicals are formed. For 0.3 g/cm3 systems the final value is approximately 0.75, while from the 0.96 g/cm3 systems the value is approximately 1.25. This indicates a mixed state of char formation in both cases, and more dehydrogenation for lower initial densities. Further work, including the introduction of a free surface, would likely be necessary to make more conclusive statements about char formation at these high thermal ramp rates. We propose here the population concentration of coordination number as a method for

10


Page 10 of 28

Page 11 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 7: Populations of carbon atoms with z=1, 2 or 3 coordination number as a function of temperature for three rates of thermal ramp in 0.96 g/cm3 (full-density) polyethylene. tracking the average structure in high-rate degradation. Figure 7 shows an example of the thermal population distributions of carbon atoms which have only one, two, or three carbon bonds for thermal ramp in full-density 0.96 g/cm3 polyethylene. We note that z = 1 indicates a terminal carbon, z = 2 indicates a backbone carbon, and z = 3 indicates a branch point or 3D structure of condensed carbon. Higher coordination numbers are not seen in high percentages, but can be followed. The atomistic visualizations in Figure 8 illustrate the populations at various temperatures for the case of our slowest ramp rate of 1011 K/s. Here, the carbon atoms are colored to represent C-C coordination, C1 (orange), C2 (blue), and C3 (green). The trends in Figure 7 are common across the different rates shown, but shifted in onset temperature. It should be noted that since the temperature is increasing in time, the population plots can be viewed as a temporal plot, by replacing the temperature on the x-axis with time. In the plots, we see C2 dominates the systems at lower temperature, as one would expect for the initial long linear chains. As temperatures rise, the C2 population drops and C1 increases indicating some combination of random and/or end scission. But, as temperatures increase further, we note a drop in C1 population, which likely accounts for an increase in C3 and C2 local coordination. Comparing initial (backbone) bond population 11



Figure 8: Molecular snapshots showing the degradation state of 0.3 g/cm3 polyethylene at 1900 K (top left) , 2100 K (top right), 2300 K (bottom left), and 2500 K (bottom right) during a 1011 K/s thermal ramp. Carbon atoms are color coded by coordination number C1 is orange, C2 is blue, C3 and higher are green. Hydrogen atoms are white, unless they form molecular hydrogen, H2 yellow. Images visualized using Ovito. 26 to the total CC bond population shows scission processes are quickly overwhelmed by CC recombination, c.f. Figure 5. This had not been previously observed at lower rates and temperatures. The role of dehydrogenation will be discussed in regards to the macrokinetic model, but we see that dehydrogenation is a step in the process of recombination. In Figure 8 we see an increase in H2 molecular hydrogen (colored yellow) with temperature. These molecular simulation results indicate that at high thermal ramp rates, additional processes are necessary to model the nonequilibrium evolution of hydrocarbon structure. We believe that the difference is due to the fact that at high rates, energies rise to the threshold of recombination before the volatile byproducts of scission can escape the system through

12


Page 12 of 28

Page 13 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


outgassing and evaporation. At low rates these processes are separable and scission leads to mass loss through vaporization before the onset of the other processes. In the next section we propose a macrokinetic model which captures the complex rate dependence of both scission, dehydrogenation, and recombination.

Macrokinetic Model We seek an accurate, but fast general continuum-scale model for polyethylene thermal degradation consistent with our MD simulations to illustrate the population dynamics in this unique temperature-time regime. In contrast to the typical assumption that gaseous products are removed from the reacting material through bulk diffusion processes (e.g., Refs. 11,12 ), our simulations show that recombination of molecular fragments must be accounted for in modeling the thermal degradation of volumetrically confined polyethylene. The consequence of volumetric confinement of the molecular ensemble appears similar to the effect of molecularlevel confinement, or “cage effects”, believed to promote recombination of polyethylene fragments in the solid-phase when there is insufficient mobility for fragments to diffuse far from scission sites. 14 However, here the effect arises strictly from rapid heating precluding bulk diffusion, which is implicit in the isochoric simulations. We chose to use a structure-independent model based on carbon coordination to simplify the rich chemistry. This circumvents the complexity of enumerating the vast range of possible chemical species, 7 while taking advantage of the ease of calculating carbon coordination from atomic positions (e.g., Fig. 7). We assumed an abbreviated number of elementary reactions: 1) carbon-carbon scission, either randomly along the chain backbone or at chain ends, 2) dehydrogenation and 3) recombination of radicals, leading to branched structures and carbon clusters. 5,15 The carbon coordination number (n) was the independent variable, assumed to change stepwise by ±1. Further, recombination was only allowed between two carbon radicals, which could be produced either from carbon-carbon scission or dehydrogenation. Reaction rates were assumed 13



to be independent of the carbon coordination (z) and charge number of carbons neighboring the reacting species. This framework is structure-independent; i.e., highly branched or even cyclical molecules are allowed but are not differentiated in the model. As well, the observation that random scission and cleavage of terminal carbons may proceed at different rates, 16 is accounted for by allowing terminal carbons to dissociate at unique rates. As a simplification, however, carbons bound to terminal carbons, are considered internal to the carbon backbone and not further differentiated. These approximations have the advantage of reducing the number of reactive species to only eight carbon entities. The result is a system of ordinary differential equations that can be solved numerically for the concentrations of carbon coordination at a given temperature and ramp rate. Arrhenius rate expressions were used to represent binary recombination and unary decomposition, incorporating temperature-dependent diffusion. Activation energies were determined from MD simulations of single-molecule decomposition and re-association, leaving the rate constants as the only fit parameters. The range of possible carbon entities in the structure-independent macrokinetic model are notated C1 , C2 , C3 , C4 , C˙ 1 , C˙ 2 , and C˙ 3 , with the subscript indicating the carbon-carbon coordination number, and a dot indicating a radical species. Note that neighboring carbons need not be notated in the chemical formulae under these simplifications. The assumed range of possible reactions are enumerated in Table 1. Whereas the recombination reactions are binary and depend on the concentrations of the two combining species, the decomposition and dehydrogenation reactions are assumed to be unary and to depend only on the concentration of the decomposing species. This simplification greatly reduces the number of free parameters and is considered sufficient for illustrating the essential mechanisms of the thermal degradation under volumetric confinement. As a further simplification, dehydrogenation and scission of bound radicals are assumed irreversible such that reaction rates correspond to net forward conversion rates. Although one could expand the set to include reactions of the type Cn C˙ m → C˙ n−1 C˙ m−1 , such species are

14


Page 14 of 28

Page 15 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


likely unstable and short-lived, and are thus ignored in the model. As a final simplification, reaction rate constants are assumed to depend only on the coordination number. Following mass conservation across the reactions in Table 1, the following set of first-order differential equations can be written to describe the population evolution: C˙ 0 dt C˙ 1 dt C˙ 2 dt C˙ 3 dt C1 dt C2 dt C3 dt C4 dt

+ ˙ 2 + ˙ + ˙ + ˙ = k1− ([C1 ] + [C˙ 1 ]) − k00 [C0 ] + k10 [C1 ][C˙ 0 ] + k20 [C2 ][C˙ 0 ] + k30 [C3 ][C˙ 0 ] + ˙ + ˙ 2 + ˙ + ˙ = k2− ([C2 ] + [C˙ 2 ]) − k1− [C˙ 1 ] + k1H [C1 ] − k10 [C1 ][C˙ 0 ] + k11 [C1 ] + k21 [C2 ][C˙ 1 ] + k31 [C3 ][C˙ 1 ] + ˙ + ˙ + ˙ 2 + ˙ = k3− ([C3 ] + [C˙ 3 ]) − k2− [C˙ 2 ] + k2H [C2 ] − k20 [C2 ][C˙ 0 ] + k21 [C2 ][C˙ 1 ] + k22 [C2 ] + k32 [C3 ][C˙ 2 ] =

k4− [C4 ]

−

k3− [C˙ 3 ]

+

k3H [C3 ]

−

+ ˙ k30 [C3 ][C˙ 0 ]

+

+ ˙ k31 [C3 ][C˙ 1 ]

+

+ ˙ k32 [C3 ][C˙ 2 ]

+

+ ˙ 2 k33 [C3 ]

+ ˙ 2 + ˙ + ˙ + ˙ k00 [C0 ] + k10 [C1 ][C˙ 0 ] + k20 [C2 ][C˙ 0 ] + k30 [C3 ][C˙ 0 ] − (k1− + k1H )[C1 ] + ˙ + ˙ 2 + ˙ + ˙ = k10 [C1 ][C˙ 0 ] + k11 [C1 ] + k21 [C2 ][C˙ 1 ] + k31 [C3 ][C˙ 1 ] − (k2− + k2H )[C2 ] + ˙ + ˙ + ˙ 2 + ˙ = k20 [C2 ][C˙ 0 ] + k21 [C2 ][C˙ 1 ] + k22 [C2 ] + k32 [C3 ][C˙ 2 ] − (k3− + k3H )[C3 ] + ˙ + ˙ + ˙ + ˙ 2 [C3 ][C˙ 0 ] + k31 [C3 ][C˙ 1 ] + k32 [C3 ][C˙ 2 ] + k33 [C3 ] − k4− [C4 ] = k30 =

+ where brackets denote concentration, knm is the reaction rate constant for recombination

involving species n and m, kn− is the net forward reaction rate constant for scission, and knH is the net forward reaction rate constant for dehydrogenation. Note that in this structure+ + , and kn− and knH do not vary with radicalization. independent model, knm = kmn

The reaction rate constants are assumed to have the Arrhenius form:

ki = ki0 fD (T ) eEi /RT

(1)

where ki0 is the number of reaction attempts per unit concentration at 298K, R is the ideal gas constant, T is temperature, and Ei is the activation energy. We have included the factor fD (T ) to account for the temperature-dependence of the reaction attempt frequency due to the temperature dependence of molecular diffusion rates, which are proportional to D1/2 for

15



Page 16 of 28

Table 1: Reactions assumed in the macro kinetic model.

Scission & Recombination

Dehydrogenation

C1 C1 * ) 2C˙ 0

C1 → C˙ 1

C2 C1 * ) C˙ 1 + C˙ 0

C2 → C˙ 2

C2 C2 * ) 2C˙ 1

C3 → C˙ 3

C3 C1 * ) C˙ 2 + C˙ 0 C3 C2 * ) C˙ 2 + C˙ 1

Scission of Bound Radicals

C3 C3 * ) 2C˙ 2

C˙ 1 → C˙ 0

C4 C1 * ) C˙ 3 + C˙ 0

C˙ 2 → C˙ 1

C4 C2 * ) C˙ 3 + C˙ 1

C˙ 3 → C˙ 2

C4 C3 * ) C˙ 3 + C˙ 2 C4 C4 * ) 2C˙ 3 a Fickian diffusion constant D. In the Chapman-Enskog theory of gases, D ∝ T 3/2 , whereas in the Stokes-Einstein theory, D ∝ T for dilute liquids with constant viscosity. 30 For our simulations spanning gas to liquid densities and a variety of molecular species, we assume the more general form: fD (T ) =

T T0

θ (2)

for scaling the reaction attempt frequency with respect to an initial temperature, T0 = 298K, and assume an intermediary value of θ = 5/8 for the binary recombination reactions and of θ = 0 for the unary decomposition reactions. To reduce the number of free parameters, activation energies were taken from a combination of literature sources and single-molecule simulations. Activation energies were taken from bond energies noting this is approximate in assuming equivalence in reaction mechanisms between isolated molecular fragments and long molecular chains at higher density. 14,31 A detailed comparison of these discrepancies may be important in some physical situations but is beyond the scope of this study.

16


Page 17 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 2: Activation energies determined from single molecule MD simulations.

Reaction

Ea (kJ/mol)

–C0

Scission

Recombination

C1 C1 * ) 2C˙ 0

–

350

18

C2 C1 * ) C˙ 1 + C˙ 0

–CH2

390

32

–CH3

339

21

C2 C2 * ) 2C˙ 1

–

339

34

C3 C1 * ) C˙ 2 + C˙ 0

–CH2

379

25

–CH3

350

22

C3 C2 * ) C˙ 2 + C˙ 1

–

350

24

C3 C3 * ) 2C˙ 2

–

333

43

C4 C1 * ) C˙ 3 + C˙ 0

–CH2

379

38

–CH3

429

80

C4 C2 * ) C˙ 3 + C˙ 1

–

507

163

C4 C3 * ) C˙ 3 + C˙ 2

–

417

113

C4 C4 * ) 2C˙ 3

-

603

289

Experimental values for bond energies for dehydrogenation of carbon structures range from 403 to 462 kJ/mol, depending on the coordination state, being among the largest for methyl groups. 32 A characteristic value of 439 kJ/mol was taken to account for the variety of structures present during polyethylene decomposition. Activation energies for the decomposition and recombination reactions were determined from MD simulations of small molecular clusters. Representative molecular fragments were first made reactive by removing hydrogen at C-C bonding sites, and then repositioned at a constant rate of approach along a linear trajectory until complete bonding of a single site occurred. Decomposition was simulated by separating the fragments along the reverse trajectory. Activation energies were determined from the difference in total simulation cell energy just before and after complete bonding. For reactions involving C0 , runs were repeated 17



Page 18 of 28

Table 3: Macrokinetic model parameters (s−1 )

Scission

Recombination

k1− = 4.0 × 1017

+ k00 =

k2− = 3.6 × 1017

+ k10 = 2.5 × 1015

k3− = 4.5 × 1017

+ k20 =

0

k4− =

+ k30 =

0

0

0

+ k11 = 1.0 × 1012

Dehydrogenation

+ k21 = 1.0 × 1016

k1H = 1.5 × 1017

+ k31 =

k2H = 5.0 × 1016

+ k22 = 5.0 × 1016

k3H = 1.0 × 1017

+ k32 =

0

+ k33 =

0

0

using both –CH2 and –CH3 termini to explore differences in bond energy. Activation energies determined in this manner are listed in Table 2. For recombination reactions, the values for individual reactions can be used unambiguously. However, because our macrokinetic model aggregates scission reactions by coordination number, activation energies characteristic of the range of reactions are required. Because the activation energy for dehydrogenation of a methyl group (438.5 kJ/mol 32 ) is significantly higher than the values for carbon scission in nearly all cases, scission of terminal –CH3 is expected to be more preponderant than scission of terminal –CH2 for initially saturated polyethylene chains. For reactions involving C1 , an exception is scission with C4 , which has an Ea value for –CH3 loss comparable to that of dehydrogenation, such that scission of CH2 is more probable. The average value, Ea = 355 kJ/mol, of the four most probable reactions involving C1 scission was used in the model. It follows that for reactions involving C2 , scission of terminal –CH2 can be ignored, as can scission of C4 due to the unusually large activation energy. The average value, Ea = 343 kJ/mol, of the lowest-energy scission pathways for C2 was used in the model. Similarly, for reactions involving C3 , scission with C4 and terminal –CH2 are 18


Page 19 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ignored, and the average value, Ea = 344 kJ/mol, was taken for C3 scission. The model is expected to be relatively insensitive to parameters for C4 , owing to their comparatively small number fraction (0.015 − 0.053) for the conditions studied. To test this sensitivity, the lower value of 379 kJ/mol for terminal –CH2 loss from C4 C1 scission was included in calculating the average value, Ea = 477 kJ/mol, for C4 scission. To compare to the MD simulations with 44-carbon chains, equations 2 through 10 were solved for the initial conditions [C1 ] = 0.045, [C2 ] = 0.955, and [Cn ] = [C˙ n ] = 0 for n = 0, 3, and 4, with concentrations normalized to the initial number of carbons. The temperature is T (t) = T0 + rt for ramp rate r. The solutions were fit to the coordination number histories from the MD simulations (Figure 7) for r = 1011 K/s and 0.3 g/cm3 density by varying the Arrhenius prefactors, ki0 . First, a grid-based least-squares algorithm was used over 1800 < T < 2750 with a wide range of initial values to allow for the possibility that not all of the reactions enumerated in Table 1 may occur, i.e., in an attempt to find the global minimum on the error surface. The parameter set with the fewest number of free parameters that gave reasonable fits to the data were then refined to best represent the population dynamics over 1800 < T < 4200 simultaneously for 1011 < r < 1013 K/s. The final parameters are listed in Table 3. Parameters with value ‘0’ indicate the model was insensitive to these reaction pathways, owing to numerical resolution of the molecular simulations. It should be noted that the rate constants, while within typical reported ranges, do not represent a universal depiction of polyethylene degradation; indeed, differing assumptions about the dominant reaction mechanisms for a given process, and the associated activation energies, is a well-known source of variability. 14 Rather, the results are used to show that an abbreviated set of scission, recombination and dehydrogenation reactions can produce the ensemble dynamics observed in the molecular simulations. The populations for C1 , C2 and C3 determined with the macrokinetic model are shown in Fig. 9. The results capture the essential population dynamics seen in the MD simulations (c.f. Fig. 7), viz., initial conversion of C2 to C1 on heating followed by increased C2 and C3

19



Page 20 of 28

concentration due to recombination of radical fragments. Although not shown for clarity, a smaller population of C4 also develops as in the simulation. The fraction of carbon bonds, θB , relative to the total carbon concentration, is: 1 θB = 2

4 X

, i[Ci ]

i=1

4 X

! [Ci ]

(3)

i=0

As shown in Figure 10, this fraction initially decreases due to scission, then increases above the original concentration due to recombination and branching reactions, reaching a peak before converging as the solutions reach equilibrium at higher temperatures where smaller species are thermodynamically favored. These dynamics markedly differ from random scission models (e.g., 8,11,12 ), wherein the number of carbon bonds decreases monotonically with temperature under thermal ramp, or with time at constant temperature, due to the absence of recombination reactions. It may also be noted that this model correctly describes the macroscopic behavior of polyethylene in terms of carbon coordination, independent of structure and without enumeration of detailed reaction pathways beyond scission, dehydrogenation and recombination. Although this framework will be too simplistic for every conceivable application, we hope that this model, which nonetheless efficiently captures the dynamics of polyethylene degradation, will find use for other polymer systems, including possibly fluorocarbons and branched

Figure 9: Macrokinetic model for population of carbon atoms with 1, 2 or 3 coordination number as a function of temperature for three rates of thermal ramp (0.3 g/cm3 density).

20


Page 21 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 10: Fraction of carbon bonds relative to the initial concentration predicted by the macrokinetic model for initial density of 0.3 g/cm3 polyethylene. alkanes under high-rate heating. The macrokinetic model is currently being implemented in a continuum solver at Sandia. In combination with a realistic EOS model, a continuum-scale multi-physics approach, incorporating nonequilibrium chemistry, will be extremely useful in comparing to high-rate experiments conducted at the national labs. The results of that work will be reported elsewhere.

Conclusions We have examined the high heating rate degradation of polyethylene molecules using simulated temperature ramps spanning 1010 − 1014 K/s. In contrast to conventional experiments, here the thermal degradation is expected to proceed isochorically, leading to substantially different degradation chemistry. Specifically, recombination of molecular fragments dominates the reaction landscape once enough radicals have formed through carbon-carbon scission or dehydrogenation, leading to molecular branching. A continuum macrokinetic model capturing these essential processes approximated the bond population dynamics seen in the MD simulations at several ramp rates. In contrast to random scission models, the concentration of carbon-carbon bonds reached a peak value higher than the initial concentration

21



due to branching, then decreased with increased temperature as smaller fragments became thermodynamically favored, highlighting essential differences in degradation chemistry of polyethylene at heating rates relevant to high energy density experiments. Finally, we note that it is our hope that the macrokinetic model proposed here will serve as a framework for broad application to other polymeric systems.

Acknowledgement Sandia National Laboratories is a multi-mission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA0003525.

References (1) Craxton, R.; Anderson, K.; Boehly, T.; Goncharov, V.; Harding, D.; Knauer, J.; McCrory, R.; McKenty, P.; Meyerhofer, D.; Myatt, J.; et al., Direct-drive Intertial Confinement Fusion: A Review. Physics of Plasmas 2015, 22, 110501–1–110501–153. (2) Moore, A. S.; Prisbrey, S.; Baker, K. L.; Celliers, P. M.; Fry, J.; Dittrich, T. R.; Wu, K.J. J.; Kervin, M. L.; Schoff, M. E.; Farrell, M.; et al., Off-Hugoniot Characterization of Alternative Inertial Confinement Fusion Ablator Materials. J. Phys.: Conf. Series 2016, 717, 012038. (3) Conesa, J. A.; Marcilla, A.; Font, R.; Caballero, J. Thermogravimetric Studies on the Thermal Decomposition of Polyethylene. J. Anal. Appl. Pyrolysis 1996, 36, 1–15. (4) Vyazovkin, S. Characterization of Materials, 2 ed.; Wiley & Sons, Inc., 2012; pp 1–12.

22


Page 22 of 28

Page 23 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(5) Beyler, C. L.; Hirschler, M. M. Thermal Decomposition of Polymers. SFPE handbook of fire protection engineering 2002, 2, 110–131. (6) Compton, R.; Bamford, C.; Tipper, C. Degradation of Polymers; Comprehensive Chem. Kinetics; Elsevier: Amsterdam, 1975; Vol. 14. (7) Levine, S.; Broadbelt, L. Detailed Mechanistic Modeling of High-Density Polyethylene Pyrolysis: Low Molecular Weight Product Evolution. Polymer Degradation and Stability 2009, 98, 810–822. (8) Knyazev, V. Effects of Chain Length on the Rates of C-C Bond Dissociation in Linear Alkanes and Polyethylene. J. Phys. Chem. A 2007, 111, 3875–3883. (9) Smith, K.; Bruns, M.; Stoliarov, S.; Nyden, M.; Ezekoye, O.; Westmoreland, P. Assessing the Effect of Molecular Weight on the Kinetics of Backbone Scission Reactions in Polyethylene Using Reactive Molecular Dynamics. Polymer 2011, 52, 3104–3111. (10) Popov, K.; Knyazev, V. Initial Stage of Pyrolysis of Polyethylene. J. Phys. Chem. A 2015, 119, 11737–11760. (11) Koptelov, A.; Milekhin, Y.; Baranets, Y. Simulation of Thermal Decomposition of a Polymer at Random Scissions of C-C Bonds. Russian J. Phys. Chem. B 2012, 6, 626– 633. (12) Koptelov, A.; Koptelov, I. Statistical Model of Thermal Degradation of Linear Polymers. Polymer Sci. B 2009, 51, 313–319. (13) Popov, K.; Knyazev, V. Molecular Dynamics Study of C-C Bond Dissociation in Linear Alkanes and Polyethylene: Effects of Condensed Phase. Eastern State Fall Technical Meeting, Chemical & Physical Processes in Combustion. 2007. (14) Popov, K.; Knyazev, V. Molecular Dynamics Simulation of C-C Bond Scission in

23



Polyethylene and Linear Alkanes: Effects of the Condensed Phase. J. Phys. Chem. A 2014, 118, 2187–2195. (15) Liu, X.; Li, X.; Liu, J.; Kong, B.; Gong, X.; Yang, X.; Lin, W.; Guo, L. Study of High Density Polyethylene (HDPE) Pyrolysis with Reactive Molecular Dynamics. Polymer Degradation and Stability 2014, 104, 62–70. (16) Nyden, M.; Stoliarov, S.; Westmoreland, P.; Guo, Z.; Jee, C. Applications of Reactive Molecular Dynamics to the Study of the Thermal Decomposition of Polymers and Nanoscale Structures. Mater. Sci. Engin. A 2004, 365, 114–121. (17) Stoliarov, S. I.; Westmoreland, P. R.; Nyden, M. R.; Forney, G. P. A Reactive Molecular Dynamics Model of Thermal Decomposition in Polymers: I. Poly (Methyl Methacrylate). Polymer 2003, 44, 883–894. (18) Bruns, M.; Koo, J.; Ezekoye, O. Population-Based Models of Thermoplastic Degradation: Using Optimization to Determine Model Parameters. Polymer degradation and stability 2009, 94, 1013–1022. (19) LAMMPS code available at http://lammps.sandia.gov. (20) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comp. Phys. 1995, 117, 1. (21) van Duin, A. C. T.; Dasgupta, S.; Lorant, F.; Goddard III, W. A. J. Phys. Chem. A 2001, 105, 9396. (22) Chenoweth, K.; van Duin, A. C. T.; Goddard III, W. A. ReaxFF Reactive Force Field for Molecular Dynamics Simulation of Hydrocarbon Oxidation. J. Phys. Chem. A 2008, 112, 1040. (23) Mattsson, T. R.; Lane, J. M. D.; Cochrane, K. R.; Desjarlais, M. P.; Thompson, A. P.;

24


Page 24 of 28

Page 25 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Pierce, F.; Grest, G. S. First-Principles and Classical Molecular Dynamics Simulation of Shocked Polymers. Phys. Rev. B 2010, 81, 054301. (24) Root, S.; Haill, T. A.; Lane, J. M. D.; Thompson, A. P.; Grest, G. S.; Schroen, D. G.; Mattsson, T. R. Shock Compression of Hydrocarbon Foam to 200 GPa: Experiments, Atomistic Simulations, and Mesoscale Hydrodynamic Modeling. J. Appl. Phys. 2013, 114, 103502. (25) Lane, J. M. D.; Grest, G. S.; Mattsson, T. R. Hot Spot and Temperature Analysis of Shocked Hydrocarbon Polymer Foams Using Molecular Dynamics Simulation. Comput. Mater. Sci. 2013, 79, 873–876. (26) Stukowski, A. Visualization and analysis of atomistic simulation data with OVITO–the Open Visualization Tool. Modelling and Simulation in Materials Science and Engineering 2009, 18, 015012. (27) Sogancioglu, M.; Yel, E.; Ahmetli, G. Pyrolysis of Waste High Density Polyethylene (HDPE) and Low Density Polyethylene (LDPE) Plastics and Production of Epoxy Composites with Their Pyrolysis Chars. J. Cleaner Production 2017, 165, 369–381. (28) Wiedemeier, D. B.; Abiven, S.; Hockaday, W. C.; Keiluweit, M.; Kleber, M.; Masiello, C. A.; McBeath, A. V.; Nico, P. S.; Pyle, L. A.; Schneider, M. P. Aromaticity and Degree of Aromatic Condensation of Char. Organic Geochem. 2015, 78, 135–143. (29) Baldock, J. A.; Smernik, R. J. Chemical Composition and Bioavailability of Thermally Altered Pinus Resinosa (Red pine) Wood. Organic Geochem. 2002, 33, 1093–1109. (30) Bird, R.; Stewart, W.; Lightfoot, E. Transport Phenomena (2nd Ed.); John Wiley & Sons, 2007. (31) Gilliom, R. Activation Energies from Bond Energies. A Modification. J. Amer. Chem. Soc. 1977, 99, 8399–8402. 25



(32) Blanksby, S.; Ellison, G. Bond Dissociation Energies of Organic Molecules. Accounts Chem. Res. 2003, 36, 255–263.

26


Page 26 of 28

Page 27 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Graphical TOC Entry

27



TOC entry 82x44mm (300 x 300 DPI)


Page 28 of 28

Molecular and Kinetic Models for High-Rate Thermal Degradation of

Recommend Documents