Article pubs.acs.org/EF
Initial Mechanisms for an Overall Behavior of Lignin Pyrolysis through Large-Scale ReaxFF Molecular Dynamics Simulations Tingting Zhang,†,‡ Xiaoxia Li,*,† Xianjie Qiao,† Mo Zheng,† Li Guo,*,† Wenli Song,† and Weigang Lin† †
State Key Laboratory of Multiphase Complex Systems, Institute of Process Engineering, Chinese Academy of Sciences, No. 1 Zhongguancun North Second Street, Beijing 100190, People’s Republic of China ‡ University of Chinese Academy of Sciences, Beijing 100049, People’s Republic of China S Supporting Information *
ABSTRACT: Initial reaction mechanisms of lignin pyrolysis were studied by large-scale ReaxFF molecular dynamics simulations (ReaxFF MD) facilitated by the first GPU-enabled code (GMD-Reax) and the unique reaction analysis tool (VARxMD). Simulations were performed over wide temperature ranges both for heat up at 300−2100 K and for NVT at 500−2100 K with a large lignin model, which contained 15920 atoms and was constructed based on Adler’s softwood lignin model. By utilizing the relatively continuous observation for pyrolysate evolution in slow heat up simulations, three stages for lignin pyrolysis are proposed by pyrolysate fractions. The underlying mechanisms for the three stages are revealed by analyzing the species structure evolution and the reactions of linkages, aryl units, propyl chains, and methoxy substituents. Stage I is characterized with the complete decomposition of source lignin molecules at low temperatures dominated by breaking of α-O-4 and β-O-4 linkages. The temperature in stage II is relatively high where cracking of all the linkages occurs, accompanied by conversion of propyl chains and methoxy substituents. Stage III mapping to high temperature shows the formation of heavy pyrolysates by recombination reactions of five-, six-, or seven-membered aliphatic rings. The heterocyclic oxygen-containing rings are revealed as important intermediates for the aryl monomer ring opening into aliphatic rings of five-membered, seven-membered, or even larger. The pathways for small molecule formation observed in this work are broadly in agreement with the literature. This work demonstrates a new methodology for investigating the overall behaviors and the underlying complex mechanisms of lignin pyrolysis.
1. INTRODUCTION
It is now well acknowledged that molecular modeling can be useful in understanding the reaction mechanisms of lignin pyrolysis at molecular level. As a major approach for investigating chemical reaction mechanisms with high accuracy, quantum mechanics (QM) method has been employed to study the thermal decomposition behavior of lignin models, such as the bond dissociation enthalpies (BDEs) for the β-5 linkage models19 and the α/β selectivity in the pyrolysis of phenethyl phenyl ethers.20 Unfortunately, the very expensive computational cost of QM would limit the model size down to ∼100 atoms; thus it can hardly be applied for large-scale lignin models that are vital to capture the complexity of lignin pyrolysis. Classical molecular dynamics (classical MD or MD) is a useful method for simulating large molecular model (∼10000 atoms or larger) to predict its properties, for example to predict the thermodynamic properties of lignin by simulating a lignin model of about 3600 atoms.21 But MD is not feasible for exploring chemical reactions in lignin pyrolysis due to the fixed chemical bonding and atom charge in simulation. As a bridge between quantum chemistry and classical molecular dynamics, the reactive force field (ReaxFF) developed by van Duin et al.22 is a bond-order based reactive potential that can be used to manage the complexity of chemical structures and the diversity of reaction pathways in
Biomass is a promising alternative and an abundant carbonneutral renewable resource for the production of bioenergy and biomaterials.1,2 As one of the three main components of woody biomass, lignin accounts for 18−40 wt % of the dry wood3 but remains not well utilized when compared with the applications of cellulose and hemicellulose in bioethanol production, pulp, and paper refineries.4,5 Lignin is an amorphous and highly branched polymer of phenylpropane units connected by ether or carbon to carbon linkages, making it probably the only viable source to produce renewable aromatic compounds.6 The challenge of lignin utilization is how lignin can be broken into smaller molecular units and be converted to more valueadded products.7 Numerous efforts8−12 have been devoted to the determination of apparent kinetics of lignin pyrolysis, mostly based on the first-order reaction model. The kinetic parameters obtained vary widely with lignin source, its isolation method, and the reactor employed for pyrolysis. The product constitutions of lignin pyrolysis are also affected by its feedstock type,13 heating rate, and reaction temperature.14 Because of the chemical complexity of condensed-phase species, it is challenging to characterize the solid- and liquid-phase intermediates in biomass pyrolysis with available experimental techniques.1 Therefore, monomeric and dimeric lignin model compounds are often employed in understanding the conversion pathways of lignin pyrolysis experimentally.7,15−18 © 2016 American Chemical Society
Received: January 30, 2016 Revised: March 23, 2016 Published: March 25, 2016 3140
DOI: 10.1021/acs.energyfuels.6b00247 Energy Fuels 2016, 30, 3140−3150
Article
Energy & Fuels
observations of density and solubility.21 In this work, a unimolecular model of lignin is constructed based on the well accepted Adler model, which is a typical softwood lignin model consisting of 16 C9 units linked by C−C or C−O bonds. The C9 unit refers to the phenylpropane skeleton having three types:44 p-hydroxyphenyl (H), guaiacyl (G), and syringyl (S) units, as shown in Figure 1. The
large reactive molecular systems. ReaxFF employs distancedependent bond-order functions22,23 to represent the contributions of bonding relevant energy terms in its potential. Both the nonbonded interactions of van der Waals and Coulomb are calculated between every atom pair irrespective of their connectivity. ReaxFF allows a smooth transition from nonbonded to single-, double-, and triple-bonded state and ensures that the energy contributions of all connectivitydependent interactions disappear upon bond dissociation. Moreover, the atomic charges are updated at each time step using the electron equilibration method (EEM)24 approach, which is similar to the QEq scheme.25,26 More details about ReaxFF potential function can be found elsewhere.22,23,27 ReaxFF has been applied in exploring chemical reactions of complicated pyrolysis systems, such as Gloeocapsomorpha prisca microfossils,28 coal,29,30 and high density polyethylene.31 Attempts of using ReaxFF MD to lignin pyrolysis study appear only very recently, including the study of beech-wood lignin as an antioxidant for petroleum asphalt by Pan et al.32 and the simulation of lignin thermal decomposition in oxygen environment by Beste.33,34 In addition to the initial mechanisms of formaldehyde formation revealed in Beste’s work, the stabilization reactions are identified for cyclic and rigid connections that are necessary for further processing into carbon fibers. The lignin oxidation model employed contains about 3000−4000 atoms, which is relatively small probably due to the concerns for computational cost. The new methodology development for ReaxFF molecular dynamics combining GPUenabled computing and reaction analysis, particularly the attempts for large-scale simulations of coal pyrolysis30 and polymer31 in the authors’ group,35 suggests that it would be promising to explore the complex chemical reactions at the initial stage of lignin pyrolysis with ReaxFF molecular dynamics. In this work, large-scale softwood lignin pyrolysis simulations were performed using the first GPU-enabled ReaxFF MD code (GMD-Reax).36 The large lignin model contains 15920 atoms and is constructed based on Adler’s lignin model,37 which allows for a more meaningful description of the lignin pyrolysis process and the underlying chemistry in the perspective of statistics. Temperature effects on the pyrolysis behavior of lignin were investigated by performing both for the heat up (300−2100 K) and isothermal simulations (500−2100 K). Heat up simulations were used to analyze the different stages in lignin pyrolysis in a view of linkage behavior and benzene ring structure evolution. NVT simulations were conducted to investigate the reaction pathways for formation of small molecules and heavy fractions. Reaction and substructure analysis of the species were investigated by VARxMD which is the first and unique tool for visualization and reaction analysis for ReaxFF MD simulations.38 Three stages for lignin pyrolysis were proposed based on the pyrolysate evolution and the underlying mechanisms of linkages, aryl rings, propyl chains, and methoxy substituents of the monomer units.
Figure 1. Chemical structures of lignin units, from left to right: phydroxyphenyl (H) unit, guaiacyl (G) unit, and syringyl (S) unit. proportions of the three units vary in different wood sources. The guaiacyl unit accounts for approximately 90% of the softwood lignin,6 while small amounts of the other two units are also found in softwood lignin.37 Among the 16 C9 units in the Adler model, there are 14 guaiacyl units, one p-hydroxyphenyl unit, and one syringyl unit. The Adler softwood lignin model represents almost all of the common linkages in lignin.6,45 There are four types of ether bonds (α-O-4, β-O4, γ-O-α, and 4-O-5) and five types of C−C bonds (β-5, β-2, β-1, β-β, 5-5) as listed in Table 1. In addition to the oxygen in methoxyl groups and ether linkages, other oxygenated groups such as hydroxyl, carbonyl and aldehyde functions are contained in the propyl chains of the model, as listed in Table 2. The construction of the lignin model containing 15920 atoms started with the creation of a unimolecule model that was optimized by using the Materials Studio (MS)46 forcite module with the DREIDING force field.47 The optimized 3D unimolecular model is shown in Figure 2. Then a cubic periodic simulation box containing 40 optimized unimolecules was built by using the MS amorphous cell module. The initial density was set to a low density of 0.2 g/cm3 to avoid overlapping of aromatic rings and other important functional groups. The box was further compressed by rescaling the box dimension to 1.3 g/cm3. After the geometry optimization, the system was equilibrated using NPT ensemble at 1 atm and 298 K for 200 ps. This procedure led to a system expansion with a density of about 1 g/ cm3. Finally, the system underwent an annealing dynamics with the recycle from 298 to 800 K using NVT ensemble to get more reasonable structure samples of the softwood lignin model. The time step for all of the model equilibration with NPT and NVT ensemble was set to 1 fs. The lignin model containing 15920 atoms in the cubic cell with an edge length of 58 Å is shown in Figure 3. 2.2. Simulation Strategy and Details. The lignin model constructed in this work contains 15920 atoms, which is about 4.6 times the value of the lignin oxidation model of 3456 atoms employed by Beste.33 Simulation of the large lignin model would require the simulation time to be at least 4.6 times longer when the same code and computational nodes are used. GMD-Reax,36 the first implementation of GPU-enabled code of ReaxFF MD created by taking advantage of GPU’s high parallelism and memory bandwidth, can offer up to one magnitude speeding up over CPU implementation of ReaxFF MD. Therefore, simulations of the large lignin model were performed efficiently with GMD-Reax running on a CentOS 5.4 server with 2*Intel Xeon E5-2650 at 2.6 GHz, 8 GB RAM, and a K20C GPU card attached. In ReaxFF MD simulations, NVT ensembles are often employed to investigate effects of temperature and reaction time on product evolution for pyrolysis process. Compared with the isothermal simulations, heat up simulations of ReaxFF MD allow for investigating the temperature effects in a wide range with one simulation, resulting in a general picture of product evolution in lignin pyrolysis process. Heat up simulations can be more computationally economical than isothermal simulations. However, a high heating rate can limit the
2. MODELING METHODOLOGY 2.1. Softwood Lignin Model Construction. Lignin models or schematic depictions of Freudenberg and Neish,39 Adler,37 Nimz,40 and Glasser and Glasser41 proposed in the last century have been widely used.42 Freudenberg’s and Neish’s lignin schematic depiction was used by Klein and Virk43 for predicting the temporal evolution of lignin thermolysis products. Both the Adler and Nimz models were adopted by Zhang and LeBoeuf for performing conventional molecular dynamics simulations, and the property prediction results of Adler’s softwood lignin model were in good agreement with experimental 3141
DOI: 10.1021/acs.energyfuels.6b00247 Energy Fuels 2016, 30, 3140−3150
Article
Energy & Fuels Table 1. Linkages and Their Amounts in the Softwood Lignin Model with 15920 Atoms
Table 2. Oxygen-Containing Functional Groups in the Softwood Lignin Model with 15920 Atoms type
number
methoxyl hydroxyl phenolic hydroxyl carbonyl aldehyde group
640 840 200 40 80
Figure 3. Softwood lignin model with 15920 atoms, containing 40 unimolecules of C160H180O58, with one unimolecule highligted.
Table 3. Simulation Details for Lignin Pyrolysis Using a Model of 15920 Atoms simulation
temperature (K)
Figure 2. Optimized 3D unimolecule model of Alder’s softwood lignin model (C160H180O58).
heat up heat up heat up
duration time at each temperature, thus leading to the slower product evolution than what is obtained in the isothermal simulations of lignin pyrolysis. To have a better understanding of the effects of temperature and reaction time on lignin pyrolysis, both the isothermal simulations at 500−2100 K for 1 ns and heat up simulations at 300−2100 K with heating rates of 0.4, 4, and 10 K/ps were performed. The ReaxFF parameters used in this study can be found in the “reax” package in LAMMPS.48,49 A detailed description of simulation conditions is listed in Table 3. During all ReaxFF MD simulations, periodic boundary conditions
NVT
300−2100 300−2100 300−2100 500−2100 (100 K as interval)
heating rate (K/ps)
simulation time (ps)
machine hours (h)
0.4 4 10
4500 450 180 1000
714 69 28 2584
were applied in all directions to the cubic simulation box with an edge length of 58 Å. The velocity-Verlet algorithm was used to integrate Newton’s equation of motion using a time step of 0.25 fs, and the temperature was controlled using the Berendsen thermostat50 with a 0.1 ps damping constant. In addition to the large lignin model employed, three parallel simulations of each temperature condition 3142
DOI: 10.1021/acs.energyfuels.6b00247 Energy Fuels 2016, 30, 3140−3150
Article
Energy & Fuels
Figure 4. Evolution of pyrolysate fractions with temperature obtained from the three parallel heat up and NVT ReaxFF MD simulations of a large softwood lignin model with 15920 atoms: (a) C160, (b) C0−C5, (c) C6−C23, (d) C24−C159, and (e) C160+. were performed by using varied initial conditions to ensure a better description of lignin pyrolysis process. 2.3. Cheminformatics Approach for Trajectory Analysis. Long time simulations of a large lignin pyrolysis system can make the analysis of ReaxFF MD simulation trajectories quite challenging because large numbers of species, radical fragments, and chemical reactions are involved. VARxMD38 was created to tackle this challenge by using a cheminformatics approach based on 3D chemical structure analysis of the simulation trajectory. VARxMD is capable of examining the complexity of the chemical reaction network in the ReaxFF MD simulations by analyzing the 3D coordinates and bond orders of atoms in a trajectory. Therefore, VARxMD was employed to reveal detailed chemical reactions and representative structure evolutions of species involved in pyrolysis of the softwood lignin model and to further uncover the relevant complex reaction mechanisms as well.
mechanisms of lignin pyrolysis, ReaxFF MD simulations both for heat up at 300−2100 K and NVT at 500−2100 K were performed. Different heating rates were employed in heat up simulations. Up to 1 ns isothermal simulations at 17 temperatures in the range of 500−2100 K were simulated. The product evolutions with temperature, heating rate, and reaction time observed in simulations are displayed in Figure 4 with the lignin pyrolysates shown in lumped fractions as (a) C160, (b) C0−C5, (c) C6−C23, (d) C24−C159, and (e) C160+. C160 refers to the source molecule in the lignin model constructed. As shown in Figure 4a, the source lignin molecules decompose almost completely when temperature is around 1100 K in the isothermal simulations and between 1100 and 1500 K in the heat up simulations at different heating rates. The light volatiles of C0−C5 in Figure 4b increase continuously with temperature indicating that high temperature favors the
3. RESULTS AND DISCUSSION 3.1. Comparison of Heat up and Isothermal Simulations. To investigate temperature effects and initial 3143
DOI: 10.1021/acs.energyfuels.6b00247 Energy Fuels 2016, 30, 3140−3150
Article
Energy & Fuels formation of C0−C5 pyrolysates. The weight percents of both C6−C23 and C24−C159 shown in Figure 4c,d have peaks in the simulated temperature range. However, as shown in Figure 4e, little production of C160+ can be observed at both lower temperatures and higher temperatures in fast heat up simulations and short time NVT simulations. Rapid mass increase of C160+ with temperature increasing are present at a low heating rate of 0.4 K/ps for heat up simulations or after 500 ps for NVT simulations when temperature is greater than 1500 K. This observation indicates that high temperatures with low heating rate or long time duration favor the production of heavy pyrolysates in lignin pyrolysis. The dependence of volatile yields and heavy fractions on the heating rate and reaction time (duration time) agrees broadly with the well accepted views for lignin pyrolysis experimentally.51 Compared with the overall pyrolysate profiles of the isothermal simulations, the product profile evolution of the heat up simulations is slower at high heating rates of 4 and 10 K/ps due to the relatively shorter reaction time (duration time) at each temperature than in NVT simulations. The higher the heating rate, the later the profile of the pyrolysates appears. The product evolution with temperature obtained at the low heating rate of 0.4 K/ps falls in the range obtained in the NVT simulations. It can be concluded that the product evolution as a function of temperature can almost resemble the NVT simulations when at a heating rate of as low as 0.4 K/ps. One such slow heat up simulation spanning with a reaction time of 4500 ps will take about 714 machine hours by GMDReax, which is significantly faster than over 2500 machine hours for a series of NVT simulations to cover the similar temperature range as shown in Table 3. Moreover, heat up simulations at a low enough heating rate allow for continuous observations of product evolution in lignin pyrolysis process with temperature and reaction time when compared with the isothermal simulations. Therefore, heat up simulations at low heating rate can be very useful in exploring the initial mechanism of lignin pyrolysis. 3.2. Three Stages of Lignin Pyrolysis. By taking advantage of the capability of continuous observations for product evolution in slow heat up simulations, three parallel heat up simulations with a heating rate of 0.4 K/ps for lignin pyrolysis were performed. The overall picture of the product evolution is shown in Figure 5. The pyrolysates obtained in the ReaxFF MD simulations are lumped using the same fraction of carbon number as in the previous section. As shown in Figure 5, C24−C159 are pyrolysates obtained directly from the cleavage of lignin molecules (C160) by linkage breaking between the lignin monomer unit structures, which is considered as “stage I”. This stage occurs in the temperature range of 700−1100 K, where the weight percent of the source lignin molecules (C160) decreases constantly to almost zero. Meanwhile, the mass of C24−C159 increases continuously to its peak accompanied by a small generation amount of C6−C23 in the same period, indicating that the lignin molecules decompose mostly into C24−C159. Next is “stage II” where the C24−C159 pyrolysates begin to crack into C6−C23 molecules by breaking of more linkages, along with the formation of a small fraction of C0−C5 pyrolysates. This stage spans at 1100− 1500 K and produces the maximal amount of C 6−C23 pyrolysates. “Stage III” starts at around 1500 K, which is characterized by further decomposition of C6−C23 into small C0−C5 pyrolysates, resulting in the constant increase of C0−C5. Meanwhile, the observations for formation and constant
Figure 5. Overall evolution picture of three stages in lignin pyrolysis obtained from the three parallel slow heat up (0.4 K/ps) simulations of a large softwood lignin model containing 15920 atoms at 300−2100 K via ReaxFF MD.
increasing of C160+ pyrolysates as well as the slight increase of C24−C159 to its second peak at around 1600 K suggest the emerging of char formation process. These observations indicate that heat up simulations of a large lignin model with a heating rate as low as 0.4 K/ps can roughly describe the overall evolution picture of lignin pyrolysis process. 3.3. Linkage Behaviors Revealed. 3.3.1. Linkage Consumption. Because of the complexity of initial reactions in lignin pyrolysis, it is challenging to understand the initial mechanism, particularly the linkage behavior of lignin pyrolysis with ReaxFF MD simulations. Attempt was made to explore the linkage behavior underlying the three stages of lignin pyrolysis. The linkages in the simulated system are identified by substructure searching of lignin monomer unit structures with the aid of VARxMD.38 The query structures used for linkage identification and the maximum consumption of each linkage in stages I and II can be found in Supporting Information Table S1. The linkage behaviors for the slow heat up ReaxFF MD simulations at 0.4 K/ps are summarized in Figure 6, where the
Figure 6. Linkage consumption in the three stages of lignin pyrolysis obtained from the three parallel slow heat up (0.4 K/ps) simulations of a large softwood lignin model with 15920 atoms at 300−2100 K via ReaxFF MD, where β-N represents the linkages of β-5, β-1, and β-2.
linkage number percentage is defined as the ratio of the unbroken linkage number to its initial number in the source lignin molecules. For simplicity, the linkages of β-5, β-1, and β2 are grouped as β-N in the following discussions, where N refers to the position number of C atom in the lignin monomer aromatic ring. 3144
DOI: 10.1021/acs.energyfuels.6b00247 Energy Fuels 2016, 30, 3140−3150
Article
Energy & Fuels As shown in Figure 6, in stage I pyrolysis of lignin (
