Reactive Molecular Dynamics Simulations of Polystyrene Pyrolysis

1. Introduction

Activated carbon is an essential industrial product with a wide range of energy, environment, biology, and medicine applications. The most commonly used raw materials for preparing activated carbon are coconut shells [1], fruit shells, anthracite, bitumen, petroleum coke, etc [2-8]. However, even with rigorous modification methods, the preparation of activated carbon from these raw materials still faces the problem of high impurity content, poor performance retention, and uncontrollable pore size. Recently, there has been an increasing number of experimental reports on the fine control of the pore structure of activated carbon using polymers and copolymers, such as polypropylene (PP), polyethylene (PE), polystyrene (PS), polyacrylonitrile (PAN), and styrene-divinylbenzene (S-DVB), as pyrolysis precursors [9, 10]. The resulting activated carbon has high purity, regular geometric morphology, and adjustable compositional content [11-15]. Nevertheless, preparing activated carbon from pyrolytic polymer precursors is plagued by unclear mechanisms and difficulties regulating pyrolysis tissues and products.

Polymers are rich in carbon, so converting polymers into high-value-added functional carbon materials is low-cost. PS is a typical polymeric material that is widely available, inexpensive, and structurally controllable, and the pyrolysis of PS precursor into functional carbon materials is considered an attractive option for the treatment of PS solids [1, 16-18]. In attempts to prepare activated carbon by pyrolysis of PS, it has been found that the heating rate and the target temperature are critical elements in determining the intermediate and final products [13, 1, 19-21]. Since microstructures of carbon materials play a vital role in the properties and characterization, the controlled pyrolysis of polymers has been studied by Tang et al. [14]. There are only a few reports on preparing carbonaceous materials using PS as a pyrolysis precursor [1, 13]. Experimentally, although PS thermal decomposition yields relatively simple products compared to other commodity polymers such as PE [22] and PP [23, 24], it is still unsatisfactory in controlling pyrolysis product type, morphology, and pore size [11].

Towards the controllable preparation of PS pyrolysis, this work adopts the reactive molecular dynamics (ReaxFF-MD) method to study the effects of heating rate and target temperature on the pyrolysis mechanism. The evolution of the microstructure and the generation of pyrolysis products is investigated in detail to provide a clear cognition of the pyrolysis mechanism from the atomic level. It is found that the heating rate would mainly affect the diversity of pyrolysis intermediates, while the target temperature would affect the tissue morphology of the final equilibrium products. Carbon ring analysis bridges theoretical simulations and experimental observations, making ReaxFF-MD simulations practically relevant.

2. Computational Details

2.1. Modeling of Pyrolysis Precursor

The pyrolysis precursor is modeled by a cubic amorphous polystyrene cell, which contains five 5 PS molecular chains, and the polymerization degree of a single chain is 20, as Figure 1(a) shows. According to the actual density of polystyrene (1.05g/cm3), the length of the cubic cell is set to 25.45Å. To make the initial structure reasonable, we constructed 10 cell configurations and optimized each configuration with the COMPASS II force field. Then, the configuration with the lowest energy was selected as the precursor model, as depicted in Figure 1(b). Subsequently, an energy minimization of the model was performed before the ReaxFF-MD simulations. In this process, the Forctice force field is utilized for the geometric optimization and annealing simulations. Specifically, 10 cycles of annealing simulations were carried out with the NVT ensemble in the temperature range 200K~500K. Thus, the precursor model was further optimized, as Figure 1(c) shows.

2.2. ReaxFF-MD Details for Polystyrene Pyrolysis

This work employs NVT integration and high-temperature techniques in ReaxFF-MD pyrolysis simulations on PS with different heating rates and target temperatures. Generally speaking, the time scale in experiments is in minutes, while a significant number of reactions in molecular simulations can occur within nanoseconds timescale. Therefore, to accelerate the reaction progress in ReaxFF-MD simulations, higher temperatures are often used to facilitate atomic motion and molecular collisions [25, 26], facilitating the observation of reactions [27]. Despite the differences in time and temperature between MD simulations and experiments, numerous studies have shown that increasing the temperature in ReaxFF-MD yields reasonable and consistent results compared to experimental outcomes, and the method has been applied to investigate combustion, pyrolysis, and explosion in larger systems at high temperatures [27-31].

ReaxFF MD simulations were carried out with the Large-scale Atomic Molecular Massively Parallel Simulator (LAMMPS) package. In all NVT-ReaxFF-MD simulations, periodic boundary conditions are applied, and the time step is set to 0.2 fs. Temperature control is achieved by the Nose-Hoover thermostat method. Structural energy minimization is performed before the ReaxFF-MD simulations, eliminating any potential strong van der Waals interactions that could lead to local structure distortions and unstable simulations. The minimization (or optimization) procedure is an iterative process in which atomic coordinates and possible cell parameters are adjusted to minimize the energy of the structure. Snapshots of the simulation results are generated by Visual Molecular Dynamics (VMD) [32] and OVITO [33] software at different time points during the process. The pore size distribution is analyzed by the Zeo++ [34] tool. The pyrolysis simulation employs the ReaxFF force field (CHON-2019). This force field was developed to characterize the dissociation and formation of carbon bonds, showing a good agreement with experimental results [35] and thus can provide a sound foundation for PS pyrolysis.

3. Results and Discussions

Preparing pyrolyzed carbon materials with well-defined microstructures requires a detailed understanding of the relationship of the pyrolysis mechanism with the heating rate and the target temperature. As widely recognized, the pyrolysis process consists of two main stages: the cleavage of polymers into small organic molecules and the carbonization of degradation products [36, 37], and we discuss both aspects below.

3.1. Effects of Heating Rate on the Pyrolysis Products

Firstly, the temperature of the PS precursor was raised to 3000K with different heating rates (5K/ps, 10K/ps, 50K/ps, 100K/ps, and 200K/ps), and we examined the effect of different heating rates on pyrolysis intermediates. We found that the PS pyrolysis occurs after 2000K, so Figure 2(a) shows the trends in the number and types of pyrolysis species during the heating process above 2000K. According to the changes in pyrolysis intermediates, the heating process can be divided into low-rate (5K/ps, 10K/ps) and high-rate heating (50K/ps, 100K/ps, 200K/ps). Overall, the number and type of pyrolysis intermediates at the low-rate heating are more significant than those at the high-rate heating. A peak occurs in the number and type during the low-rate heating, which involves the further cross-linking transformation of intermediate products obtained from the early pyrolysis of PS into larger carbon-hydrogen macromolecules. However, the number and type of pyrolysis intermediates show a negative correlation with the heating rate during the high-rate heating, which indicates that a higher heating rate suppresses the thermal decomposition of the precursor.

Then, 6ns isothermal simulations at 3000K were performed for both low- and high-rate heating systems. As Figure 2(b) shows, in this stage, the total number of pyrolysis intermediates increases while types decrease for every heating rate, which suggests that diverse molecules during the heating process undergo secondary cross-linking reactions that eventually converge to produce the same final pyrolysis equilibrium products. Therefore, the heating rate will not affect the equilibrium carbon products at 3000K after the isothermal process.

A more detailed pyrolysis products and tissue analysis further prove the above discussion. As shown in Figure 3, hydrocarbon intermediates increase for high-rate heating, while the number and type of hydrocarbon intermediates peak and gradually decrease for low-rate heating. As a result, the heating rate affects the type and quantity of pyrolysis intermediates and their appearance time. In experiments, Onwudili et al. [21] reported the pyrolysis of PS in an intermittent autoclave, where PS was degraded at around 350°C, mainly into viscous, dark-colored oil. Carbon formation increased slightly until 425°C. However, carbon formation increased significantly at 450~500°C. In the absence of a catalyst, PS undergoes pyrolysis at temperatures above 500 degrees Celsius, resulting in the production of styrene and styrene-like compounds, as well as the Formation of coke [19, 21]. Westerhout et al. concluded that even small temperature changes could lead to large changes in reaction rates during PS pyrolysis [38]. Thus, the experiments are highly consistent with our simulations.

3.2. Effects of Target Temperature on the Pyrolysis Products

To test the effects of target temperatures on the equilibrium pyrolysis products, we first heated the system to 2500K, 2750K, 3000K, 3250K, 3500K, and 3750K. As described in the previous section, the heating rate hardly affects the equilibrium pyrolysis products, so the 50K/ps heating rate was implied for every simulated system. Subsequently, a 6ns isothermal simulation was followed at the respective target temperature and all simulated systems reached equilibrium, as proved in Figure S2 in SI.

The evolution of pyrolysis tissues and products at the isothermal simulation process is shown in Figure 4. According to the changes in microscopic morphology and pyrolysis products, these simulated systems can be classified into low-temperature system (LTs, 2500K), proper-temperature system (PTs, 2750K, 3000K, 3250K), high-temperature system (HTs, 3500K), and ultra-high-temperature system (UHTs, 3750K).

Comparing the four sets of simulations, we found that the time required for pyrolysis to reach equilibrium is proportional to the target temperature, i.e., the higher the target temperature, the shorter the time for pyrolysis to reach equilibrium. However, the process of pyrolytic evolution and the final equilibrium products are quite different. In the LTs, the equilibrium pyrolysis products are predominantly organized in cross-linked networks with discrete layers of graphite flakes of small molecules. In the PTs, the pyrolysis products show large graphite lamellar tissues and further evolve to randomly folded carbon networks when the pyrolysis equilibrium is reached. In the HTs, large graphite lamellar tissues appear more significantly during the isothermal simulation and evolve into a regular and hierarchical graphitized carbon network at the end of the simulation. However, if the temperature is extremely high, no graphitized tissues are formed, and cracked small molecules are haphazardly dispersed throughout the system at the end of the simulation, as Figure 4 shows. Thus, excessively low and high target temperatures are not conducive to generating large graphitized tissues.

3.3. Tissue Characterization and Carbon Ring Analysis of Pyrolysis Products

In the previous two sections, we discussed the effects of heating rate and target temperature on the pyrolysis intermediates and final products. However, we need a quantitative characterization that can connect with the experiments to discover the microstructure's evolution during the pyrolysis process, which is crucial for fine-tuning the pyrolysis products in experiments. The previous study showed that the type and quantity of carbon rings in pyrolysis tissues can be fitted by the component peaks in X-ray photoelectron spectroscopy (XPS) [39]. Therefore, in the following, carbon ring analysis is employed to elaborate on the pyrolysis mechanism.

3.3.1. Carbon Ring Evolution with Different Heating Rates

As described in section 3.1, the temperature was raised to 3000K at 5K/ps, 10K/ps, 50K/ps, 100K/ps, and 200K/ps heating rates, respectively. The changes in the various types of carbon rings during the heating process at different heating rates are shown in Figure 5. There are three main types of carbocycles generated during pyrolysis, i.e., 5-, 6-, and 7-membered rings, which we denote by 5r, 6r, and 7r, respectively. Since there are only 6r in the PS pyrolysis precursor, the change in its number represents the pyrolysis process, while the presence of 5r and 7r represents the generation of new pyrolysis products.

Figure 5 shows that higher heating rates can inhibit the early cleavage of the 6r and the production of 5r and 7r. Consistent with the previous analysis, this suggests that the low-rate heating (5K/ps, 10K/ps) contributes to the sufficient cleavage of the precursor, resulting in the production of type-rich pyrolysis intermediates, as 5r and 7r are significantly observed. In sharp contrast, 6r changes extremely slowly in the high-rate heating plots, with a very small number of pyrolysis intermediates, and no new product formation is even observed at a heating rate of 200 K/ps. The analysis of carbon rings further proves that low-rate heating helps yield more diverse pyrolysis intermediates.

The changes in carbon rings during the next thermostatic simulation at 3000K are shown in Figure 6. As can be seen, there is a clear inflection point in 6r plots. At the early stage of the thermostatic simulation, the cracking of the precursor continues and 6r keeps decreasing. During this period, the magnitude of the decrease in the 6r plots is proportional to the rate of heating, i.e., the higher the heating rate, the greater the magnitude of the decrease in the 6r plots, indicating insufficient cleavage at the heating stage is compensated in the thermostatic stage.

Subsequently, accompanying the completion of the precursor pyrolysis and the secondary cross-linking reaction, graphite-type tissues are generated in large quantities, as evidenced by the fact that the amount of 6r starts to rise rapidly and becomes the main structure in the final product, outnumbering both 5r and 7r. Since the heating rate hardly affects the equilibrium product, the same equilibrium tissue is obtained at 3000 K for the five heating rates.

3.3.2. Carbon Ring Evolution at Different Target Temperatures

The evolutions of carbon rings in six thermostats with different target temperatures, as defined in section 3.2, were plotted in Figure 7. We can see that the 6r curves for the other five systems, except for the UHTs (3750K), show inflection points in the early stage of the isothermal simulations, indicating that further cleavage of the tissues from the heating process is still required at the beginning of the isothermal. It was found that the higher the target temperature, the faster the further cleavage of the precursor, which is due to the fact that the 6r decline curve becomes steeper with increasing temperature during the initial period of isothermal, as shown in the individual subplots. However, due to the excessively high temperature in UHTs, the precursor is completely cleaved and no inflection point is observed in the system.

After the inflection points, the secondary cross-linking reaction happens, accompanied by a rapid increase in the number of 6r in LTs (2500K), PTs (2750K, 3000K, 3250K), and HTs (3500K) and a stable evolution of 5r and 7r. The number of 6r in PTs and HTs is obviously larger than that in LTs. Since graphite has a typical 6-membered ring structure, this change of 6r in the system after the inflection point means that the tissues undergo a graphitizing process, which is sufficient in the PTs and HTs. Regarding the main tissue morphology and the number of carbon rings, the final equilibrium products in the PTs and HTs are similar, consisting of large graphitized lamellae. However, as will be discussed in the next chapter, the pore size distributions of these tissues are quite different. In the UHTs, the number of 6r is almost the same as 5r and 7r, indicating that there is a large amount of discrete fragmented organization in the system and no large-scale graphite lamellar structure is generated. Therefore, the evolution of the carbon ring under different conditions can be used as a descriptor and measured value by experiments in the controlled pyrolysis process of polystyrene.

3.3.3. Pore Size Analysis

We analyze the pore sizes and their distribution in the respective pyrolysis products to further differentiate the difference in equilibrium tissues at different target temperatures. The pore size distribution in the equilibrium product at each target temperature is shown in Table 1. The detailed evolution of pore distribution can be found in Figure S7 in SI.

According to the previous discussion, we know that the LTs and UHTs contain a large number of discretized fragmented molecules, so the conductivity will not be good. We focus on the better-graphitized PTs and HTS systems. Inside the tissues of these systems, 5~11Å pore size dominates, about 75%. In particular, the 11 Å apertures increase while 5~9 Å apertures decrease with temperature in the PTs. This phenomenon helps to discern tissue nuances in the PTs.

4. Conclusions

Understanding the heat mechanism is essential to realize controlled pyrolysis of polystyrene. In this paper, we systematically investigated the effect of heating rate and target temperature on the pyrolysis tissues of polystyrene using ReaxFF-MD simulations. The study found that the heating rate would mainly affect the diversity of pyrolysis intermediates, while the target temperature would affect the tissue morphology of the final equilibrium products. From a computational point of view, a temperature increase to 3000~3500 K at a rate of 5~10 K/ps followed by constant temperature pyrolysis yields abundant pyrolysis intermediates as well as an extensive graphitized lamellar tissue. However, constrained by the computational methodology and model size, this by no means implies that the experiments should strictly follow these numbers. Importantly, we have revealed through theoretical simulations the processes that inevitably occur in experiments and are challenging to track in real time. In this way, by comparing the experimental observations with the theoretical simulations, we can know the pyrolysis stage and decide the direction of the next step, thus achieving controlled pyrolysis. In the future, we hope to apply this method to prepare amorphous pyrolyzed carbon materials with metallic and nonmetallic loadings, thus promoting the development of new energy materials.