Toward a Process-Based Molecular Model of SiC Membranes. 2

Jan 28, 2013 - Because free hydrogen radicals are in high energy state, they have a strong .... the average coordination number of Si and/or C atoms, ...
0 downloads 0 Views 5MB Size
Download PDF
Article pubs.acs.org/JPCC

Toward a Process-Based Molecular Model of SiC Membranes. 2. Reactive Dynamics Simulation of the Pyrolysis of Polymer Precursor To Form Amorphous SiC Saber Naserifar,† William A. Goddard, III,*,§ Lianchi Liu,‡ Theodore T. Tsotsis,† and Muhammad Sahimi*,† †

Mork Family Department of Chemical Engineering & Materials Science, University of Southern California, Los Angeles, California 90089-1211, United States § Materials and Process Simulation Center, California Institute of Technology, Pasadena, California 91125, United States ‡ Computational Chemistry Group, School of Chemistry and Chemical Engineering, Shanghai Jiao Tong University, Min Hang, Shanghai 200240, China ABSTRACT: In part 1 of this series we developed the reactive force field ReaxFF, choosing the parameters adjusted to fit quantum mechanics description of prototypical reactions. In the present paper we use ReaxFF for reactive dynamics (RD) simulation of thermal decomposition of a siliconcontaining polymer, hydridopolycarbosilane (HPCS) over a wide range of temperature. Many properties of the material during its pyrolysis were computed and were found to be in good agreement with the experimental data. We then simulated pyrolysis of the HPCS under conditions that closely mimic the conditions of the fabrication of nanoporous SiC membranes to produce amorphous SiC material. Again, the computed properties of the SiC ceramic were found to be in good agreement with the experimental data. In particular, the ReaxFF RD simulations generate ceramic SiC that contains mostly Si−C bonds with some Si−Si and C−C bonds, consistent with experimental data. This indicates again the accuracy of the parameters of the ReaxFF determined in Part 1.

1. INTRODUCTION In Part 1 of this series1 we described the development of the ReaxFF reactive force field, in order to study the pyrolysis of hydridopolycarbosilane (HPCS), [SiH2CH2]n. As described in Part 1, the pyrolysis process is an essential step in fabrication of silicon-carbide (SiC) nanoporous membranes. In the present paper we utilize ReaxFF with reactive dynamics (RD) simulation to study the temperature-dependence of the decomposition of the HPCS. We carry out RD simulation of the pyrolysis process under the conditions that mimic closely the experimental conditions in which the SiC nanoporous membranes are fabricated. This provides a strong test of the accuracy of ReaxFF, by further comparing the properties of the predicted material - amorphous SiC - with the relevant experimental data. The rest of this paper is organized as follows. Section 2 briefly describes the generation of an atomistic model of the HPCS. How the polymer model is utilized in the MD simulations of its pyrolysis is described in Section 3. Molecular dynamics simulation of thermal decomposition of the HPCS is described in Section 4, where the results are presented, dicussed, and compared with the experimental data. In section 5 the formation of amorphous SiC by MD simulations, under the conditions that closely mimic the experimental conditions in which nanoporous SiC membranes are fabricated, is © 2013 American Chemical Society

described. Section 6 compares the properties of the resulting SiC ceramic with the experimental data, while the paper is summarized in Section 7.

2. ATOMISTIC MODEL OF THE POLYMER As the first step we generate an atomistic model of the HPCS, in order to carry out MD simulations of its pyrolysis using ReaxFF, and study the effect of pyrolysis on the final atomistic model and composition of the SiC material. The molecular model of the HPCS was generated using the polymerconsistent force field (PCFF). The procedure for the generation of the polymer’s molecular model as well as the parameters of the PCFF were given in Part 11 and, thus, are not repeated here. 3. ENERGY MINIMIZATION, PRE-EQUILIBRATION, AND HEATING After generating the atomistic structure of the HPCS, we carried out reactive molecular dynamics (RMD) calculations using the ReaxFF force field to describe the thermal degradation process over a wide range of temperature. The Received: August 6, 2012 Revised: January 23, 2013 Published: January 28, 2013 3320

dx.doi.org/10.1021/jp307799p | J. Phys. Chem. C 2013, 117, 3320−3329

The Journal of Physical Chemistry C

Article

K, 3200 K, and 3400 K for which the simulations lasted 60 ps. The time step, the minimum absolute temperature deviation, and the rescaling coefficient were the same as before. We then computed the temperature-dependence of the populations of the types of the bonds formed in the final configuration that were obtained from the (NVT)-MD simulation after 20 ps. The results for two temperature intervals are shown in Figures 1 and 2. Although the pyrolysis in the fabrication of SiC nanoporous membranes is at about 1000 K,2,3 our RMD simulations indicated that, below 1800 K, there were very small changes in the number of bonds over time scales of a few nanoseconds, practical for the RMD. Thus, higher temperatures are required for the pyrolysis process to commence on the practical computational time scale. At and above 1800 K we find that Si−H and C−H bonds break often to form hydrogen radicals. Then, these radicals form H−H bonds to produce hydrogen gas. As Figure 2 indicates, above 2400 K the RMD leads to breaking of the Si−C bonds, which leads to the formation of CH2 and SiH2 radicals. They are responsible for the production of methane, silane, and methyl silane gases by forming bonds with the available H radicals. At and above 2400 K only a very small number of Si−Si bonds remains in the system. With increasing temperature the numbers of H−H and C−C bonds also increase, resulting in an increase in the number of C2 hydrocarbons and hydrogen molecules, which agrees with the experimental data. Analysis of the experimental data for the gaseous products of pyrolysis of the HPCS indicated4 that H2 is the major product of the pyrolysis, while methane, silane, and methyl silanes gases are also produced in small amounts. The final configuration of the system at 3600 K, obtained by (NVT)-RMD simulation, is shown in Figure 3. The complete list of the products at the same temperature was provided in Table 11 of Part 1.1 Figure 4 compares the populations of the gaseous products in the final configuration after RMD simulations for 20 ps, in the range 2400 K−4000 K. This indicates that H2 is the major product of the pyrolysis. Similar results were obtained for the populations of the CHn and SiHn (n = 1−4) radicals and molecules in the same system. To understand better how the populations of hydrogen radicals and hydrogen molecules H2 form over time, we present in Figure 5 time- and temperaturedependence of their populations. Thus, one may view each pair of hydrogen radicals as being equivalent to one H2 molecule, because over longer periods of time (not currently practical in RMD simulations) the hydrogen atoms in high energy states bond together and produce H2. This insight, which greatly facilitates the MD simulations, will be used later for removing the hydrogen atoms from the system during pyrolysis simulation of the HPCS at a specified temperature. Similarly, time- and temperature-dependence of the populations of other species were also computed during the (NVT)-MD simulations. For example, Figure 6 presents the evolution of production of CHn (n = 1−4) radicals and the final methane molecules at 3200 K. At 3200 K the population of CH3 increases with the time and produces methane by reacting with hydrogen, consistent with the increase in the number of CH4 molecules during the same period. The trends are similar at 3600 and 4000 K (not shown). The effect of high temperatures on cleavage of C−H bonds becomes clear at 4000 K, at which they break more easily, and, therefore, after about 5 ps the numbers of CH3 radicals and CH4 molecules begin decreasing. Similar trends were obtained for SiH3 radicals and silane molecules.

RMD simulations of the pyrolysis process were carried out under the conditions that closely mimic the experimental ones under which the SiC nanoporous membranes are fabricated. The procedure for the cook-off simulations of the polymer was described in Part 1.1 Before heating up the polymer, energy minimization and preequilibration simulations were carried out using ReaxFF, in order to further relax the molecular structure of the HPCS. All the simulations in this section were implemented using GRASP - General Reactive Atomistic Simulation Program - a MD simulator developed by Sandia National Laboratories. The GRASP simulator is written in FORTRAN language and can be executed in parallel computations. The total energy E of the polymer was minimized using the steepest-descent method with relaxation in all the directions. The maximum displacement allowed in any direction during any time step was set to 0.01 Å, while the time step was 0.1 fs. To further check whether the true equilibrium structure of the polymer was obtained, the total energies of the atoms were also examined along the dynamics. The results indicated the constancy of all such energies. Following energy minimization, we equilibrated the structure using isothermal-isochoric (NVT)-MD simulations for 5 ps at 50 K, during which temperature was first rescaled to 50 K and then held fixed. The time step was 0.1 fs. The minimum absolute temperature deviation for which rescaling was performed was 10 K, and after each 10 time steps the temperature was checked and compared with the target value. The rescaling coefficient (which lies between 0 and 1) was 0.7, so selected carefully to rescale the velocities to obtain the target temperature.

4. TEMPERATURE-DEPENDENCE OF DECOMPOSITION OF HPCS To analyze the kinetics and details of the pyrolysis process, three series of (NVT)-RMD simulations were carried out. The

Figure 1. Population of the bonds in the final configuration obtained from the (NVT)-MD simulations. Also given is the population for the HPCS in the bulk before raising the temperature.

first series consisted of six simulations in the range 300−1800 K, divided into 300 K intervals, while the second series included three simulations in the range 2000−2800 K with 400 K intervals. Eighteen distinct simulations were carried out in the last series, which were in the range 3000 K−4000 K with 200 K intervals. They were repeated twice for each temperature. Duration of each simulation was 20 ps, except for those at 3000 3321

dx.doi.org/10.1021/jp307799p | J. Phys. Chem. C 2013, 117, 3320−3329

The Journal of Physical Chemistry C

Article

Figure 2. Population of the bonds in the final configuration obtained from the (NVT)-MD simulations, for the temperature range 2400 K−4000 K. Also given is the population for the HPCS in the bulk before raising the temperature.

Figure 4. Populations of gaseous molecules in the final configurations and their dependence on temperature. Figure 3. The final configuration of the material at 3600 K, obtained by (NVT)-MD simulation. The density is 0.98 g/cm3.

4.1. Kinetics of HPCS Decomposition. The hydrogen radicals represent an intermediate product in the decomposition of the polymer and proceed in a fast secondary reaction to form H2 as the primary product of the pyrolysis. Therefore, we utilized production of H2 as a basis for studying the kinetics of the HPCS decomposition. Temperaturedependence of the rate k of production of H2 in the interval 3000 K−4000 K was computed. The results are presented in Figure 7. Three trials runs were used to obtain a measure of the

As for the CHnSi (n = 1−6) radicals and molecules, our simulations indicated that increasing temperature increases cleavage of Si−H and C−H bonds, and, thus, more CHSi and CH2Si radicals are produced. There is a competition between the CHnSi (n > 2) radicals and breaking of the bonds in Si−H and C−H, resulting in large fluctuations in the population of the radicals around a mean value. 3322

dx.doi.org/10.1021/jp307799p | J. Phys. Chem. C 2013, 117, 3320−3329

The Journal of Physical Chemistry C

Article

Figure 5. Time- and temperature-dependence of production of hydrogen radical and hydrogen molecule during (NVT)-MD simulations.

Figure 6. Time-dependence of production of CHn, during (NVT)-MD simulation at 3200 K.

errors in the estimated values. The results were then fitted to the classical Arrhenius functional form

k(T ) = k 0 exp( −Ea /RT ) 3323

(1)

dx.doi.org/10.1021/jp307799p | J. Phys. Chem. C 2013, 117, 3320−3329

The Journal of Physical Chemistry C

Article

Figure 9. Time-dependence of the pressure and temperature during MD simulations for heating up four HPCS strands at a density of 0.98 g/cm3.

Figure 7. Temperature-dependence of rate of production of H2. The results for three separate (NVT)-MD simulation runs with distinct initial configurations are shown, in order to provide a measure of the variations. Straight line represents fit of the data from trial 1.

where Ea is the activation energy, and R is the gas constant. We obtained, k0 ≃ 1.28 × 1012 mol·L−1·s−1 and Ea ≃ 39.7 kcal/mol.

5. FORMATION OF AMORPHOUS SIC The final step of modeling the process of formation of amorphous SiC is to carry out MD simulations under the

Figure 10. Time-dependence of the potential energy during (NVT)MD simulations. The target temperature is 2000 K.

Figure 8. The schematic of shell script for removing the hydrogen radicals from the system.

conditions that mimic the experimental setup for fabrication of SiC nanoporous membranes.2,3 Selecting the temperature at which the MD simulations are to be carried out is a crucial issue, because the time scale of the simulation, being at most on the order of nanosecond, is much shorter than the actual experiments where the pyrolysis of the polymer take several hours to complete. In order to describe the dramatic changes in geometries and charges, and large changes in density and atom connectivity that require valence angle and torsion potential during the simulation of formation of amorphous/crystalline SiC, it is necessary to use short time steps, 0.1 fs. Therefore, the temperature for the RMD simulations must be selected high enough to accelerate the computations to allow the chemical reactions to occur on the computational time scale, though the quality of nanoporous SiC membrane, produced as a result of pyrolysis of the polymer

Figure 11. The increase in the density of the system during the (NPT)-MD simulation to compress the material.

on the pores’ surface of the membrane’s support, might be slightly affected with varying the pyrolysis temperature. On the other hand, experimental data5 indicate that for temperatures below 1175 K the SiC ceramic that is produced by the pyrolysis process is in the amorphous state, in which Si− H and C−H bonds break readily, and cross-linking between hydrogen radicals results in the production of H2 gas, the main 3324

dx.doi.org/10.1021/jp307799p | J. Phys. Chem. C 2013, 117, 3320−3329

The Journal of Physical Chemistry C

Article

Figure 12. The increase in the temperature and pressure. The target temperature is 2000 K.

Figure 13. Time-dependence of the temperature and pressure during (NVT)-MD simulation to lower the temperature from 2000 to 300 K.

product of the pyrolysis. Below 1175 K experiments show6 that the breaking of Si−C bonds is very limited. As the populations of the bonds and the analysis of the gasesous products over a broad range of temperatures revealed (see above and Part 1), the RMD simulations show that decomposition of the polymer begins at 1000 K by breaking the Si−H and C−H bonds and is enhanced by increasing the temperature. Below 2400 K the RMD simulations find that Si−C bonds remain unbroken, though few of them are very mobile and continuously form and break bonds at temperatures close to 2400 K, but above which the Si−C bonds begin breaking considerably. Therefore, the simulation temperatures of 2400 K and lower are, roughly speaking, equivalent to the experimental temperatures of 1175 K and lower at which the ceramic product is in the amorphous state and the decomposition process occurs mainly by breaking the Si−H and C−H bonds. Since in the fabrication of SiC nanoporous membranes2,3 the pyrolysis of the HPCS is done at 1025 K, we decided to carry out the (NVT)-MD simulations at 2000 K. Moreover, during the membrane fabrication, the support’s pores, coated by the HPCS, are heated at 1025 K for several hours in a tubular furnace, into which argon flows to sweep all the gaseous products of the pyrolysis out of the furnace. Thus, similar to the experiments, all the gaseous products of the pyrolysis must be deleted from the MD simulation cell, but doing so requires some additional considerations, in order to reduce the computation time to a reasonable level. As pointed out earlier, gas chromatography and mass spectroscopy have indicated4 that H2 is the major product of the pyrolysis, with traces of methane, silane, and methyl silanes gases also present, which are in agreement with our MD simulation results presented earlier. Because free hydrogen radicals are in high energy state, they have a strong tendency to

form H−H bonds, implying that after long enough simulation all the hydrogen radicals react and form H2. Thus, every two hydrogen atoms may be viewed as being equivalent to one H2 molecule, and, therefore, deletion of the hydrogen radicals during the pyrolysis simulations is equivalent with deleting H2 gas that would form at much later computational time scales. This is an important consideration, as it is computationally very expensive to carry out RMD simulations for long enough times to allow all the hydrogen radicals to form H2 molecules. Thus, we used shell programming to integrate the computer program with the GRASP package to automatically delete the H radicals at 2−5 ps intervals. The schematic of the shell script for removing the H atoms is shown in Figure 8. The external program evaluates the distance of every H atom in the system from all the neighboring atoms - C, H, and Si - and if it exceeds the equilibrium bond length between H and the corresponding neighboring atom by a factor of 1.2, it considers the H atom to be unbounded and deletes it from the simulation cell. Then, the program creates a new input file for the next cycle of the (NVT)-RMD simulation. The total time used to remove all the H radicals from the simulation cell was 1.16 ns. Therefore, extensive MD simulations were carried out in the (NVT) ensemble at 2000 K. Following energy minimization and pre-equilibration (see above and Part 1), the HPCS was heated up to 2000 K. The density of the polymer at the beginning of the simulation was 0.98 g/cm3, consistent with experimental data. The simulations were carried out using a heating rate of 0.1 K/fs. The pressure change during the simulation was 4.24 GPa; see Figure 9. The time step, minimum absolute temperature deviation, and the rescaling coefficient were selected similar to those for the (NVT)-MD simulations at 50 K (see above). 3325

dx.doi.org/10.1021/jp307799p | J. Phys. Chem. C 2013, 117, 3320−3329

The Journal of Physical Chemistry C

Article

using a time step of 0.1 fs with the Nosé-Hoover thermostat (with a 10 fs damping constant) and the Rahman-Parrinello barostat (with a 100 fs pressure damping constant). The structure was compressed to 1.1ρ, where ρ is the density of the ceramic, 2.2 g/cm3, and then the simulation cell was allowed to relax and expand at ambient temperature and pressure, reaching the right density. The procedure allows movement and arrangement of the atoms in the structure. Figure 11 presents the results from the (NPT)-MD simulation at 2000 K for the density, while Figure 12 presents the results for the variations of the pressure and temperature. The pressure change during the simulation was set to 8 GPa to reach the required density after 84 ps of simulation (selected carefully based on extensive preliminary simulations), in order to obtain the required density. Next, the system was cooled down to the ambient temperature with a cooling rate of 0.1 K/fs (using the same time step of 0.1 fs and thermostats). The results are shown in Figure 13. As shown, the system required further equilibration at ambient temperature and pressure, to release the extra pressure and achieve a structure in equilibrium. Therefore, additional (NPT)-MD simulations with the same parameters as before were utilized to relax the MD cell. To avoid sudden expansion of the cell, the pressure was lowered gradually from 8 GPa to ambient pressure, which also provided enough time for the atoms to rearrange themselves in the system and achieve a better equilibrium state. Figure 14 indicates that after 60 ps of simulation the pressure has reached its ambient value, with the volume of the system expanding slowly with time. The system was relaxed for an extra 60 ps of simulation using (NPT)-MD simulation to make sure that it had reached its true equilibrium state. Figure 15 presents the structural evolution of the system after removing the H radicals, carrying out 818 ps of (NVT)MD simulations, and finally compressing, cooling down, and relaxing the system at ambient condition. Subsequently, simulated annealing8 in the (NVT) ensemble was utilized to provide the system with the opportunity for surmounting the energetic barriers, in a search for conformations with energies lower than the possibly localminimum energy identified by the energy minimization step. The system was annealed between 300 and 3000 K. It was first heated to 3000 K using a heating rate of 0.5 K/fs, stayed at 3000 K for 5 ps, after which it was cooled off to 300 K with the rate of 0.1 K/fs; see Figure 16. The state of this system at this point was taken to be the true equilibrium state of the material.

Figure 14. The change in the density and pressure during (NPT)-MD simulation. The pressure is released from its high value, and temperature is kept at 300 K.

After removing all the 1288 H atoms, the system consisted of 320 Si and 320 C atoms, which were relaxed by carrying out (NVT)-MD simulations using the GRASP at 2000 K for a sufficient time to allow for the rearrangement of the atoms and cross-linking between Si and C atoms. All the parameters of the simulations were similar to those at 50 K. The potential energy shown in Figure 10 remained unchanged after 820 ps, implying that the system had reached its true equilibrium state. During this time period the Si and C atoms begin to bond together, with Si−C rings appearing in the structure. Experimentally, the final product of pyrolysis of the HPCS at 1025 K is an amorphous ceramic composed of nearstoichiometric SiC, which also contains a small percentage of hydrogen. The density of the ceramic resulting from the pyrolysis of the polymer has been reported7 to be 2.2 g/cm3. The volume of the simulation cell obtained from the last step of the (NVT)-MD simulation at 2000 K is the same as the one corresponding to the density of the HPCS at the beginning of the simulation, namely, 0.98 g/cm3. Therefore, the simulation cell should be compressed enough to achieve the density of the amorphous SiC ceramic. Thus, we carried out isothermalisobaric (NPT)-MD simulations to compress the cell at 2000 K,

6. FURTHER COMPARISON WITH THE EXPERIMENTAL DATA In addition to the comparisons with the experimental data and knowledge described above, we analyzed the final SiC ceramic and computed its radial distribution function, g(r), in order to compare it with the experimental data. Figure 17 presents the results at 300 K. The location of the largest peak corresponds to the Si−C covalent bond and is in very good agreement with the experimental data.9 In addition, experimental studies indicated the existence of some homonuclear bonds, Si−Si and C−C, in amorphous SiC ceramic. The presence of such bonds is indicated by the computed g(r). The peaks corresponding to C−C and Si−Si covalent bonds are consistent with experimental data.10 These results confirm the accuracy of the MD simulation and, hence, ReaxFF, in producing the true arrangement of Si and C atoms in the material. 3326

dx.doi.org/10.1021/jp307799p | J. Phys. Chem. C 2013, 117, 3320−3329

The Journal of Physical Chemistry C

Article

Figure 15. The structure of the SiC ceramic after (a) removing all the hydrogen radicals from the system during (NVT)-MD simulation at 2000 K; (b) (NVT)-MD simulation and equilibration of the system at 2000 K, and (c) (NPT)-MD simulation for compressing the system at 300 K. For better clarity, the size of the simulation cell in (c), which is smaller than those in (a) and (b), was set equal to them.

Figure 17. The computed radial distribution function g(r) of the SiC ceramic, as well as those for the individual bonds in the material, obtained by MD simulation. The experimental data were adopted after Musumeci et al.14

The radial distribution function may also be used to gain information about the coordination number Z of each atom in the material.11 The relation between Z and g(r) is given by

Figure 16. Variations of the properties of the system during further energy minimization, using simulated annealing. 3327

dx.doi.org/10.1021/jp307799p | J. Phys. Chem. C 2013, 117, 3320−3329

The Journal of Physical Chemistry C

Article

temperature. The time step was 0.1 fs, and the Andersen thermostat was used. In the next step the Si atoms were saturated and energy minimization with 10 ps long (NVT)-MD simulations were carried out again. This resulted in better arrangement of the newly added hydrogen and avoided sudden changes in the energy and pressure of the system. The final structure contained 29.3 wt % C, 68.4 wt % Si, and 2.3 wt % H, a composition in agreement with the experimental data4,12 for the amorphous structure of SiC, obtained by pyrolysis at 1025 K. This was followed by MD simulations in the (NVT) ensemble at 300 K for 50 ps, in order to further relax the saturated structure. Then, MD simulations were carried out in the (NPT) ensemble for several ns at 1 atm and 300 K. At the end of the simulations, the structure of the SiC ceramic was validated by computing the X-ray diffraction (XRD) pattern of the material. The results are shown in Figure 18 and compared with the experimental data.14 The computed locations of the peaks, which correspond to the arrangements of the Si and C atoms in the amorphous phase of the material, are close to the corresponding experimental angles. Thus, the computed results are in good agreement with the data, hence further confirming the accuracy of the results and, hence, the ReaxFF force field. The next step of this work is to use the predicted SiC ceramic for simulation of transport of gaseous mixtures through it, in order to test the accuracy of the simulated material in terms of its separation properties.

7. SUMMARY The reactive force field ReaxFF, which was developed in Part 11 for simulating thermal decomposition of polymers that contain Si, was utilized to carry out extensive and detailed MD simulation of thermal decomposition of hydridopolycarbosilane (HPCS) over a wide range of temperature. Various properties of the system during the pyrolysis were computed. Whenever a comparison with experimental data was possible, we found good agreement between the results provided by the MD simulations and the data. We then simulated the pyrolysis of the HPCS under the conditions that closely matched the experimental one used in the fabrication of nanoporous SiC membranes, in order to generate the amorphous SiC material used in the membranes. Once again, several properties of the ceramic were computed and were found to be in good agreement with the experimental data. In particular, consistent with the experimental data, the MD simulations generate ceramic SiC that also contains some Si−Si and C−C bonds. Thus, the parameters of the ReaxFF that we have determined provide an accurate force field for simulating thermal decomposition of an important class of Si-based materials, which have many applications in science and technology.

Figure 18. The computed XRD pattern of the SiC (top) and (bottom) experimental XRD patterns of SiC at various temperatures.14

Z = 4π

∫0

rm

ρg (r )r 2dr

(2)

where rm is the first minimum point after the largest peak in g(r), related to the Si−C bond. This equation provides an estimate of the average coordination number of Si and/or C atoms, where ρ = N/V, with N being the number of the atoms, V is the volume (in Å3), and r is the radial distance (in Å). We obtained an average coordination number of 3.6 for Si/C atoms, indicating that most of the atoms must have a coordination number of 4, with the rest having Z < 4. More precisely, bond distance analysis of the final SiC ceramic produced by the simulation indicated that about 48.1% of the Si atoms and 70.0% of the C atoms have a coordination number of 4. The rest of the atoms have a coordination number that is mainly 3. Experimental data4,12−14 have revealed that the amorphous SiC, prepared by various methods, contains some hydrogen with a concentration of up to 40−50 atomic percent. Therefore, we may conclude that the remaining Si and C atoms with coordination numbers less than 4 are saturated with hydrogen atoms. Hence, we added hydrogen atoms to the system such that each of them had the longest possible distance from the other atoms connected to the same central atom. This was done by computing first the directional vectors that connect other bonded atoms to the central atom and then setting the H atoms in the reverse direction of the sum of all the directional vectors with the same magnitude of the X-H bond length with X = Si and C. Adding hydrogen atoms was done step by step. First, the carbon atoms were saturated with hydrogen, followed by minimization and 10 ps long (NVT)-MD simulation at room

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The support of the Office of Energy Sciences, Chemical Sciences, Geosciences & Biosciences Division of the Department of Energy is gratefully acknowledged. T.T.T. also acknowledges the support of the National Science Foundation. 3328

dx.doi.org/10.1021/jp307799p | J. Phys. Chem. C 2013, 117, 3320−3329

The Journal of Physical Chemistry C

Article

W.A.G. received support from the Department of Transportation.

■

REFERENCES

(1) Naserifar, S.; Liu, L.; Goddard, W. A.; Tsotsis, T. T.; Sahimi, M. J. Phys. Chem. A 2013, DOI: doi.org/10.1021/jp3078002. (2) Elyassi, B.; Sahimi, M.; Tsotsis, T. T. J. Membr. Sci. 2007, 288, 290. (3) Elyassi, B.; Sahimi, M.; Tsotsis, T. T. J. Membr. Sci. 2008, 316, 73. (4) Interrante, L. V.; Whitmarsh, C. W.; Sherwood, W. Mater. Res. Soc. Symp. Proc. 1994, 346, 593. (5) Ly, H. Q.; Taylor, R.; Day, R. J.; Heatley, F. J. Mater. Sci. 2001, 4045. (6) Tang, X.; Zhang, L.; Tu, H.; Gu, H.; Chen, L. J. Mater. Sci. 2010, 45, 5749. (7) Starfire Systems, Inc. http://www.starfiressystems.com (accessed Jan 31, 2013). (8) Kirkpatrick, S.; Gelatt, C. D.; Vecchi, M. P. Science 1983, 220, 671. (9) CRC Handbook of Chemistry and Physics, 92nd ed.; 2012. http:// www.hbcpnetbase.com/ (accessed Jan 31, 2013). (10) Ishimaru, M.; Bae, I.-T.; Hirotsu, Y.; Matsumura, S.; Sickfus, K. E. Phys. Rev. Lett. 2002, 055502. (11) McQuarrie, D. A. Statistical Mechanics; University Science Books: Sausalito, CA. (12) Hasegawa, Y.; Okamura, K. J. Mater. Sci. 1983, 18, 3633. (13) Soraru, G. D.; Babonneau, J. D.; Mackenzie, J. D. J. Mater. Sci. 1990, 25, 3886. (14) Musumeci, P.; Calcagno, L.; Makhtari, A. Mater. Sci. Eng., A 1998, 253, 296.

3329

dx.doi.org/10.1021/jp307799p | J. Phys. Chem. C 2013, 117, 3320−3329