Oxidation of Silicon Carbide by O2 and H2O: A ReaxFF Reactive

Jul 3, 2012 - Three-site cluster analysis of oxidation of SiC for oxidation with H2O at 3000 K (a) clusters that track the progress of the primary spe...
0 downloads 8 Views 3MB Size
Download PDF
Article pubs.acs.org/JPCC

Oxidation of Silicon Carbide by O2 and H2O: A ReaxFF Reactive Molecular Dynamics Study, Part I David A. Newsome,† Debasis Sengupta,*,† Hosein Foroutan,‡ Michael F. Russo,‡ and Adri C. T. van Duin‡ †

CFD Research Corporation, 215 Wynn Dr., Huntsville, Alabama 35805, United States Department of Mechanical and Nuclear Engineering, Pennsylvania State University, 136 Research East Building, Bigler Road, University Park, Pennsylvania 16802, United States

‡

S Supporting Information *

ABSTRACT: Simulations of the initial oxidation process of a SiC surface exposed to O2 and H2O molecules was studied with ReaxFF, an atomically detailed reactive molecular dynamics method that naturally models the breaking and forming of bonds. In this work, the ReaxFF forcefield was first expanded by training it with new quantum mechanics data of the binding energy, equation of state, and heat of formation of the SiC crystal, along with data from earlier studies that describes Si − Si, Si − O, and Si − H interactions. This expanded ReaxFF forcefield is capable of simultaneously describing both Si−C−O and Si−O−H bonding interactions. Using the forcefield, oxidation simulations were performed at various temperatures (in the range of 500 to 5000 K), and the trajectories were analyzed. The analyses showed that SiC gradually transforms into the oxides of silicon with simultaneous formation of a graphite-like layer. In presence of excess O2, the graphite-like layer was further oxidized to CO and CO2. We also analyzed the trajectories with two-atom and three-atom clusters to quantitatively track the oxidation progress. This analysis clearly showed Si−O and C−C bond formation at the expense of O−O and Si−C bond consumption indicating SiC oxidation with simultaneous formation for carbon-like structure. Consumption of SiC with O2 was found to be faster than that with H2O. We have also reported the oxidation simulation of SiC with a mixture of H2O and O2. Oxidation proceeded effectively as a twopart sequence, with O2 first oxidizing the SiC, followed then by H2O. ature and high voltage electronics.7,8 Therefore, interaction of SiC with oxidizing gases is of a more general interest. The focus of this article is therefore to study the initial stage of the SiC oxidation reactions of SiC using computational chemical simulation tools. The interaction of the C−SiC−SiO2−rubber composite with an oxidizing environment, such as in the presence of oxygen and water vapor, is an extremely complex process. There are several experimental and modeling (both engineering and molecular modeling) reports on pure SiC oxidation available in the literature. Experimental research efforts on SiC can be further divided into two categories:9 (1) thermal oxidation of SiC by water vapor for electronics application9−12 and (2) oxidation of SiC for high temperature structure materials.13−16 During the oxidation process, oxide of silica (SiOX) is formed, which creates a protective layer on the SiC and hinders further oxidation. Generally, oxidation of SiC is discussed in terms of the two processes:13,14,16−18 (1) passive oxidation and (2)

1. INTRODUCTION Space vehicles endure an extremely harsh oxidative environment during the re-entry due to the friction with the atmosphere. This generates very high temperatures and pressures, requiring a thermal coating on the external surfaces of the space vehicles in order to withstand such extreme conditions. Over the years, a variety of reinforced composite materials have been developed as thermal protection system (TPS) materials that undergo a controlled ablation process. Composite materials, based on carbon, silicon carbide (SiC), silica (SiO2), and carbon fiber, typically combined with ethylene−propylene−diene monomer (EPDM) rubber, have been of particular interest as candidate thermal protection materials by aerospace agencies, such as Navy and NASA.1−6 One of the current interests of the scientific and engineering community is to understand the response and behavior of this class of materials when exposed to very high temperatures and pressures. A deeper and more fundamental understanding of the degradation mechanism of thermal protection materials is necessary in order to develop the next generation of TPS. In addition to the space vehicles application, SiC is also used extensively in the semiconductor industries for high temper© 2012 American Chemical Society

Received: June 28, 2012 Revised: July 2, 2012 Published: July 3, 2012 16111

dx.doi.org/10.1021/jp306391p | J. Phys. Chem. C 2012, 116, 16111−16121

The Journal of Physical Chemistry C

Article

metrized against a large number of high-level ab initio calculations and experimental results that include molecular structures, heats of formation, bond energies, activation barriers, transition state structures, crystal structures, etc. The ReaxFF potential has been applied in various areas involving pure gas phase reactions, gas-surface reactions, and liquid-phase reactions. In the present article, we first outline the ReaxFF potential development for atoms involved in the oxidation of SiC, i.e., the Si/C/H/O system. We then examine the high temperature oxidation of SiC with pure O2, pure H2O, and with a mixture of H2O and O2 using ReaxFF MD simulations in the temperature range of 500 to 5000 K. Each simulation trajectory was then analyzed to understand the intermediate steps during the oxidation process via computing the numbers of various types of bonds and species along the trajectory. The ReaxFF MD simulations therefore provided a fundamental insight into the high temperature oxidation of SiC at the initial stages of reactions. The TPS materials containing SiC are exposed to highly oxidative environment at very high temperatures (>3000 K) that are difficult to achieve at laboratory experiments. Therefore, our primary interest in this article is to understand the thermal response of SiC in oxidative environment especially at very high temperatures. In another article (to be communicated separately), we have utilized the ReaxFF trajectory data reported in this article to compute the rate constants of oxidation. We developed rate equations and demonstrated how the rate theory can be used in conjunction with time-consumption profile data to compute the growth of the silica layer upon SiC oxidation. We then computed activation barriers and pre-exponential factors for the oxidation reactions from Arrhenius analysis.

active oxidation. The passive oxidation is observed at high partial pressures of the oxidizing agent in which the protective oxide layer prevents further oxidation, and the oxide is formed along with the formation of CO as SiC(s) + 3/2O2 (g) → SiO2 (s) + CO(g)

Active oxidation occurs at low partial pressure of the oxidizing agent, and no protective films are expected due to the formation of SiO vapor formation: SiC(s) + O2 (g) → SiO(g) + CO(g)

There are many experimental reports on the growth of the oxide layer and are mostly performed below 2000 °C with oxygen, water vapor, carbon dioxide, or a mixture of two or more oxidizing gases,16 and the overall rate constants of the oxidation are evaluated using the standard Deal−Grove model.19 However, none of these efforts have assisted in developing a molecular level understanding of the initial stages of oxidation and the intermediate steps leading to the formation of CO. Apart from experimental efforts, several modeling papers have also appeared that considered both the macroscopic and molecular level. The macroscopic models were primarily based on a thermodynamic approach2,17 and transport model.20,21 In the thermodynamic approach, a prespecified number of chemical species were considered and Gibb’s free energy of the system was minimized to obtain the composition of the oxide and other gas species. The model reported by Hijikata et al.20,21 was based on a transport theory where Si and C emission from the oxidation surface was assumed to be a process during the oxidation. In this model, the rate parameters were then extracted by fitting the experimental growth rate vs oxide thickness plot. However, the molecular level studies of SiC oxidation are primarily based on first-principles ab initio methods.22−24 Gavrikov et al.23 reported some detailed studies of SiC oxidation and proposed a mechanism for CO formation on a clean SiC surface. These authors also performed a kinetic Monte Carlo simulation using the proposed mechanism. Soon et al.24 have performed ONIOM QM/MM calculation using a cluster model of SiC to examine the initial stages of oxidation with a clean surface. Although these calculations produced some valuable information of SiC oxidation, the dynamics of the oxidation process at different temperatures and pressures could not be obtained. In addition, static ab initio calculations performed in their report required a priori knowledge of the possible reactions. In contrast, dynamical calculations require no prior knowledge of the reactions, such as the reactive molecular dynamics simulations reported here, where reactions are taken into account inherently. The reactions occur naturally with time when initial conditions, such as temperature, pressure. and number of atoms, are set. Therefore, performing reactive dynamics simulations is perhaps the most natural way of theoretically studying a chemically reacting system. The dynamic evolution of a chemically reacting system can be studied most accurately via ab initio based MD methods, such as CPMD.25 The forces in this technique are calculated on the fly without any parametrization and approximation. However, such calculations are computationally too expensive and can be applied only to small sized atomic systems, typically containing a few tens of atoms. During the past decade, the group of van Duin developed a sophisticated bond-order based analytical potential, the ReaxFF method, which is capable of handling bond-breaking and bond-making.26−37 The parameters involved in the ReaxFF potential functions were para-

2. REAXFF POTENTIAL DEVELOPMENT FOR SI/C/H/O SYSTEM The ReaxFF method26 is an empirical force-field that allows for fully reactive atomistic scale molecular dynamics simulations of chemical reactions. The parameters used in the ReaxFF method are usually derived by fitting against a training set comprising both quantum mechanics (QM) and experimental data. The ReaxFF method is unique among reactive force field approaches in that it obtains good accuracy for both reaction energies and reaction barriers. Furthermore, the ReaxFF method is several orders of magnitude faster than QM-based simulations, which enables larger-scale (≫1000 atoms), fully dynamical simulations of chemical reactions. This combination between accurate energies and low computational expense makes the ReaxFF method ideally suited for the evaluation of complex reacting systems using molecular dynamics. To enable application of the ReaxFF method to silicon carbides, we combined previously derived ReaxFF descriptions for hydrocarbon oxidation26 and silicon/silicon oxides.37,38 To extend the set of existing ReaxFF force field parameters for C/ H/O and Si/O/H materials to the silicon carbide materials, we combined the Si/C, Si/O, Si/Si, and Si/H quantum data used previously for the silicon,28,29 silicon oxide,37 and polydimethylsiloxide ReaxFF applications.30 To these data sets, we added QM/DFT data describing the equation of state and heat of formation of the silicon carbide crystal (in the diamond configuration) as well as DFT data reported by Wang et al.22 describing the binding energy of oxygen at various locations on the silicon carbide surface. Using this combined Si/Si, Si/O, Si/ H, and Si/C QM training set, we determined ReaxFF Si−C 16112

dx.doi.org/10.1021/jp306391p | J. Phys. Chem. C 2012, 116, 16111−16121

The Journal of Physical Chemistry C

Article

Figure 1. QM (DFT) and ReaxFF equation of state of silicon carbide in the diamond- configuration.

Figure 2. QM and ReaxFF binding energies of four SiC[001] surface sites22.

article, we first evaluated its performance for the overall reaction energies related to the materials in question. These ReaxFF parameters predict a heat of formation of −215.5 kcal/ mol for the SiO2 alfa-quartz phase, which is in good agreement with the experimental value reported at the NIST website (−216.4 kcal/mol). Furthermore, the ReaxFF parameters predict a heat of formation for the SiC−diamond phase of +8.3 kcal/mol, which is in reasonable agreement with the DFT number (+13.6 kcal/mol) but somewhat removed from the experimental value (−17.1 kcal/mol39). Finally, these parameters predict heats of formation for CO and CO2 of −13.3 and −88.1 kcal/mol, respectively, which are in reasonable agreement with the values reported by Cox et al.40 of −26.4 and −94.1 kcal/mol, respectively. As such, for the overall reaction studied in this work, which is expressed below, ReaxFF predicts an energy release of 220.5 kcal/mol, in comparison with 225.7 kcal/mol energy release based on the experimental heats of formation of the reactants and products.

bond, Si−C−O angle, and C−Si−O angle parameters. Since the QM training set also contained all of our SiO2, Si, and C/ H/O data from earlier work, this approach ensures that the new S/C/H/O parameters can not only describe silicon carbide materials but also the silicon oxides, diamond, and graphite.26 It should be pointed out that the parameters are slightly modified due to the inclusion of an additional training data set. The complete forcefields are provided as Supporting Information. Figures 1 and 2 compare the ReaxFF results for the silicon carbide equation of state and for the oxygen binding to the QM data. We observed that ReaxFF provided a very accurate reproduction of the silicon carbide material around equilibrium. At high expansion (>50% expansion), ReaxFF overestimated the energy penalty; however, even at elevated temperature, we would not anticipate observing a >50% expansion in crystal volume, indicating that these highly expanded crystal states are not directly relevant for the current study. We also included the heat of formation of the SiC−diamond phase in our training. To validate the ReaxFF Si/C/H/O parameters described in this 16113

dx.doi.org/10.1021/jp306391p | J. Phys. Chem. C 2012, 116, 16111−16121

The Journal of Physical Chemistry C

Article

SiC(s) + 3/2O2 (g) → SiO2 (s) + CO(g)

The ReaxFF parameters used in this work are a combination between previously published ReaxFF descriptions for Si, SiO2, SiC, and C/O/H systems. For all these subsystems, this parameter set has demonstrated quite reasonable accuracy for observables directly related to reaction barriers. For example, Neyts and co-workers used the Si/O parameter set to study the kinetics of silicon oxidation41−44 and found good agreement with experimental growth mechanisms directly related to reaction kinetics. The ReaxFF Si-description was used by Buehler and co-workers29,45,46 to study crack-propagation in alpha-silicon and found excellent agreement between crackspeeds of ReaxFF and experiments, which indicates a good description of reaction barriers. Chenoweth et al.30 used similar ReaxFF parameter sets to study the reaction kinetics of PDMSpyrolysis, a polymer containing Si−C bonds, and found good agreement between ReaxFF and available experimental results for this polymer. The C/O/H parameters employed here were trained against a significant number of relevant reaction barriers,27 typically reproducing reaction barriers within 5 kcal/mol of their DFT values, and were successfully applied for studying hydrocarbon pyrolysis and combustion kinetics by a number of authors.32,47−53 For this reason, we believe that the parameter set applied here is likely to give a reasonable description of the kinetics associated with silicon carbide oxidation. Figure 2 demonstrates that ReaxFF gives an accurate description of the oxygen binding energy to the surface of a silicon carbide slab at four distinct positions. This oxygen binding can be seen as the initial step of the oxidation process. After this development stage, we performed a series of simulations in which a hydroxyl-passivated silicon carbide slab was exposed to a high-pressure O2 and H 2O gas at temperatures ranging from 500 to 5000 K. The force field parameters used in this study are available through the Supporting Information. These force field parameters can be used with the open-source LAMMPS code, which includes the ReaxFF method.54 2.1. ReaxFF MD Simulation Setup. To explore the structural modifications, we prepared a periodic silicon carbide slab continuous in the x- and y-orientations and interfaced with a high pressure bath of O2 and/or H2O molecules, as depicted in Figure 3 for the O2 case. Since the bare silicon carbide surface is highly reactive, we populated both the silicon carbide slab surfaces with hydroxyl groups, thus ensuring that all carbon and silicon atoms have their normal coordination of four at the beginning of the simulation. This system was first equilibrated at 100 K, primarily to remove short contacts in the O2 gas phase related to their random placement, after which the system temperature was raised quickly to target temperatures of 500, 750, 1000, 1250, 1500, 1750, 2000, 2500, 3000, 3500, 4000, 4500, 5000, and 5500 K. The temperature was controlled using a Berendsen thermostat,55 with a temperature damping constant of 100 fs. A time step of 0.25 fs was used in all simulations up to 3500 K, above which the time step was reduced to 0.1 fs. All simulations were performed in the NVT (Canonical) ensemble, using a fixed, orthogonal, periodic cell of 15.92 × 12.07 × 50 Å and containing 680 atoms. For comparison, the cutoff radius is 10 Å. The slab itself is more than twice this length, at around 25 Å. For the sake of demonstrating how ReaxFF can model surface oxidation of SiC, this structure size was deemed sufficient to model a sufficient

Figure 3. RMD simulation setup used in our simulations. The surface atoms of the initial structure were passivated by an −OH group.

number of oxidation events. Increasing the size may improve the statistics but would not change our overall observations. Though two of the periodic box lengths are less than twice the actual cutoff radius, the oxidation is actually dominated by the covalent bonds forming and breaking, so this effect would be minor, as the covalent bonds only cover a few angstroms. In addition, ReaxFF considers all image interactions, not only the minimum image, as such, apart from a small artificial symmetry effect, the relatively small system size will not cause any computational inaccuracies.

3. RESULTS AND DISCUSSION 3.1. Oxidation Simulation of SiC with O2. Figures 4 and 5 show the potential energy and the O2 removal kinetics, as observed during our O2/SiC oxidation simulations. These figures indicate that reactions involving O2(g) removal were not the only exothermic events in these simulations, as the potential energy kept decreasing even after all O2 molecules were consumed. For example, at 2000 K, all O2 were consumed at 200 000 fs (or 200 ps); however, the potential energy continues to decrease until the very end of the simulation. The potential energy, in principle, should reach a constant value once the system reaches equilibrium. However, our objective is to study the transient behavior of the oxidation processes, and therefore, no attempt has been made to run the simulation until the equilibrium is reached. This observation can be explained by looking at the structural evolution of the SiC/O2 system during the simulation. Figure 6 shows this evolution during the 3000 K simulation. We 16114

dx.doi.org/10.1021/jp306391p | J. Phys. Chem. C 2012, 116, 16111−16121

The Journal of Physical Chemistry C

Article

Figure 4. Potential energy as observed during the ReaxFF MD silicon carbide oxidation simulations.

Figure 5. O2 consumption as observed during the ReaxFF MD silicon carbide oxidation simulations.

H2O, which were formed very early during the simulation and, as such, can be described as early kinetic products. The kinetics of SiC slab oxidation can be roughly subdivided into a fast surface oxidation, coupled with a slower phase-separation and associated graphitization event. The solid-state chemistry was driven by the higher oxidation potential of silicon compared to carbon, causing the silicon to out-compete the carbon for the oxygen atoms. By studying a series of temperatures, we have specifically aimed at addressing the temperature influence on oxidation kinetics. Our observations seem to suggest the formation of a noncrystalline SiO2 even at lower oxidation temperatures, which seems reasonable given the tendency of this material to form glasses. Given the glassy nature of this material, no strong melting transition or related modified

observed a very fast oxidation of the SiC material, resulting in the formation of an amorphous or liquid-like silica suboxide. Although we are unable to conclude the phase of the oxide at this stage, considering the melting point of silicon dioxide (which is