Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
A process-based molecular model of nano-porous silicon carbide membranes
(USC Thesis Other)
A process-based molecular model of nano-porous silicon carbide membranes
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
A PROCESS-BASED MOLECULAR MODEL OF NANO-POROUS SILICON CARBIDE MEMBRANES by Saber Naserifar A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulllment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (CHEMICAL ENGINEERING) August 2013 Copyright 2013 Saber Naserifar Acknowledgments First and foremost, I wish to express my gratitude to my advisors, Professors Muhammad Sahimi and Theodore T. Tsotsis who were abundantly helpful and oered invaluable assistance, support and guidance. I am also grateful to Professor William A. Goddard, III, of Caltech and his group for their fruitful collaboration, generous help, and very useful advice at various stages of this work. My deepest gratitude are also due to the member of my Ph.D. supervisory committee, Professor Stephan Haas and also members of qualifying committee Professors Andrea M. Armani, Stephan Haas and Katherine Shing without whose knowledge and assistance this study would not have been successful. Special thanks also to my graduate student friends both at USC and Caltech for sharing experience and assistance, as well as the sta members and my friends at Mork family department for making a friendly and enjoyable environment. I would also like to thank the National Science Foundation and the Department of Energy for providing the nancial support. Last, but not the least, I wish to express my love and gratitude to my beloved family; for their support, understanding and endless love, through the duration of my studies. ii Contents Acknowledgments ii List of Figures vi List of Tables xv Abstract xvii 1 Introduction 1 1.1 Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Fabrication Methods of Silicon Carbide (SiC) Membranes . . . . . . 2 1.3 Pyrolysis Mechanism and Film Growth . . . . . . . . . . . . . . . . 3 1.4 Modeling and Simulation . . . . . . . . . . . . . . . . . . . . . . . . 7 1.5 The Plan of the Dissertation . . . . . . . . . . . . . . . . . . . . . . 10 2 Background and Computational Methods 11 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.2 Quantum Chemistry . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.2.1 The Hartree-Fock approximation . . . . . . . . . . . . . . . 15 2.2.2 The Roothaan equations . . . . . . . . . . . . . . . . . . . . 16 2.2.3 Electron correlation . . . . . . . . . . . . . . . . . . . . . . . 21 2.2.4 General energy models . . . . . . . . . . . . . . . . . . . . . 22 2.2.5 Density functional theory . . . . . . . . . . . . . . . . . . . 23 2.3 Molecular Dynamics and Monte Carlo Simulations . . . . . . . . . . 26 2.4 ReaxFF: A Reactive Force Field . . . . . . . . . . . . . . . . . . . . 26 2.4.1 Bond order calculation . . . . . . . . . . . . . . . . . . . . . 27 2.4.2 Charge calculation . . . . . . . . . . . . . . . . . . . . . . . 29 2.4.3 Potential energy . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.4.4 One-body interactions . . . . . . . . . . . . . . . . . . . . . 32 2.4.5 Two-body interactions . . . . . . . . . . . . . . . . . . . . . 34 2.4.6 Three-body interactions . . . . . . . . . . . . . . . . . . . . 37 2.4.7 Four-body interactions . . . . . . . . . . . . . . . . . . . . . 40 iii 2.5 Dual Control Volume-Grand Canonical Molecular Dynamics . . . . 42 2.5.1 The DCV-GCMD method . . . . . . . . . . . . . . . . . . . 43 2.5.2 Potential models of the gases and the membrane . . . . . . . 47 2.5.3 Transport properties . . . . . . . . . . . . . . . . . . . . . . 48 3 Development of a Reactive Force Field 50 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.2 Quantum Chemical Calculations . . . . . . . . . . . . . . . . . . . . 51 3.3 Optimization of ReaxFF . . . . . . . . . . . . . . . . . . . . . . . . 53 3.4 Atomistic Model of the Polymer . . . . . . . . . . . . . . . . . . . . 60 3.5 Cook-o Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 4 Amorphous SiC Formation via Reactive Dynamics Simulation of the HPCS 79 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 4.2 Atomistic Model of the Polymer . . . . . . . . . . . . . . . . . . . . 80 4.3 Energy Minimization, Pre-equilibration, and Heating . . . . . . . . 80 4.4 Temperature-Dependence of Decomposition of HPCS . . . . . . . . 82 4.5 Kinetics of HPCS Decomposition . . . . . . . . . . . . . . . . . . . 86 4.6 Formation of Amorphous SiC . . . . . . . . . . . . . . . . . . . . . 88 4.7 Further Comparison with the Experimental Data . . . . . . . . . . 101 5 Diusion of Gases in the Amorphous SiC Model 109 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 5.2 Molecular Modeling and Simulation . . . . . . . . . . . . . . . . . . 111 5.2.1 Gas insertion in the SiC model . . . . . . . . . . . . . . . . 112 5.2.2 The force eld . . . . . . . . . . . . . . . . . . . . . . . . . . 113 5.3 Accessible Free Volume and Solubility . . . . . . . . . . . . . . . . . 114 5.4 Estimation of the Self-diusivities . . . . . . . . . . . . . . . . . . . 115 5.5 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . 116 5.5.1 Morphological properties . . . . . . . . . . . . . . . . . . . . 116 5.5.2 Self-diusion coecients . . . . . . . . . . . . . . . . . . . . 117 6 Application to Transport and Separation of Gaseous Mixtures in SiC Membranes 129 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 6.2 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . 130 6.2.1 Details of the simulations . . . . . . . . . . . . . . . . . . . 130 6.2.2 Steady-state condition . . . . . . . . . . . . . . . . . . . . . 132 6.2.3 Eect of the temperature . . . . . . . . . . . . . . . . . . . . 137 6.2.4 Eect of the external pressure drop . . . . . . . . . . . . . . 139 6.2.5 Comparison with the experimental data . . . . . . . . . . . 141 iv 7 Conclusions 156 References 159 A Bond and Angle Distribution for HPCS 171 B Temperature-Dependence Production of H 2 175 B.1 Smoothing the Population Data . . . . . . . . . . . . . . . . . . . . 175 B.2 Reaction Rate at Various Temperature . . . . . . . . . . . . . . . . 176 C Separation Factor Fitting and Extrapolation 195 v List of Figures 1.1 Schematic illustration of SiC fabrication process with a SEM picture of thin SiC lm on the surface of the support. . . . . . . . . . . . . 4 1.2 Schematic illustration of the eect of temperature change on the pyrolysis products of the polycarbosilane precursor. . . . . . . . . . 6 2.1 Schematic illustration of the atomic congurations in energy terms: (a) one-body, (b) two-body, (c) three-body, (d) four-body, (e) hydrogen- bonding, and (f) non-covalent iterations. A red sphere indicates the position of the primary atom in each energy term. Solid lines indi- cate covalent bonds, while non-covalent bonds are represented by dotted lines. Picture courtesy of [1]. . . . . . . . . . . . . . . . . . . 32 2.2 Schematic illustration of the SiC amorphous membrane used in the DCV-GCMD simulations. The h and ` areas are, respectively, the high- and low-pressure control volumes. . . . . . . . . . . . . . . . . 43 3.1 Comparison of the energies computed by ReaxFF and QC for (a) Si{C and (b) double-bond dissociations. . . . . . . . . . . . . . . . 54 3.2 Comparison of distortion energies computed by ReaxFF and QC for the angles in (a) Si{C{Si; (b) Si{C{H; (c) C{Si{H; (d) C{C{Si; (e) C{Si{C, and (f) C{Si{Si triplets in the HPCS. . . . . . . . . . . . . 55 3.3 Comparison of the Mulliken charge distributions computed by the DFT (red) and ReaxFF (blue). . . . . . . . . . . . . . . . . . . . . 56 3.4 The equations of state (compression and expansion) for the three crystal structures of SiC. (a) -SiC; (b) 4H-SiC, and (c) -SiC, as computed by the DFT and ReaxFF. . . . . . . . . . . . . . . . . . 58 3.5 Time-dependence of the potential energy of the HPCS during (NPT )- MD simulation at 300 K. The force eld used was the PCFF. . . . 64 3.6 The computed radial distribution function of the HPCS. . . . . . . 65 vi 3.7 The amorphous structure of the HPCS, as generated by the PCFF. 65 3.8 Distribution of (a) Si{C distances; (b) Si{C{Si angles, and (c) H{ C{Si{C dihedral angles in the HPCS polymer, as generated by the PCFF. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 3.9 Time-dependence of the energy of the HPCS and the system's pres- sure during minimization using the ReaxFF. . . . . . . . . . . . . . 68 3.10 Energies of all the atoms of the polymer at the end of minimization, using ReaxFF. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 3.11 Time-dependence of the energy and temperature of the HPCS, and the system's pressure during pre-equilibration at 50 K, using ReaxFF. 70 3.12 Energies of all atoms of each type at the end of pre-equilibration (at 50 K), using ReaxFF. . . . . . . . . . . . . . . . . . . . . . . . . . . 71 3.13 Time-dependence of the temperature and pressure during the cook- o simulation of the HPCS, using ReaxFF. . . . . . . . . . . . . . . 73 3.14 The evolution of the number of bonds during the cook-o simulation of the HPCS, using ReaxFF. . . . . . . . . . . . . . . . . . . . . . . 74 3.15 The evolution of the number of the hydrogen radicals during the cook-o simulation of the HPCS, using ReaxFF. . . . . . . . . . . . 75 3.16 The evolution of the number of the CH n (top), SiH n (middle) and CSiH n (bottom) radicals during the cook-o simulation of the HPCS, using ReaxFF. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 3.17 The evolution of the C 2 hydrocarbons during the cook-o simulation of the HPCS, using ReaxFF. . . . . . . . . . . . . . . . . . . . . . . 77 3.18 The evolution of various compounds during cook-o simulation of the HPCS, using ReaxFF. . . . . . . . . . . . . . . . . . . . . . . . 78 4.1 Population of the bonds in the nal conguration obtained from the (NVT )-MD simulations. Also given is the population for the HPCS in the bulk before raising the temperature. . . . . . . . . . . . . . . 83 4.2 Population of the bonds in the nal conguration obtained from the (NVT )-MD simulations, for the temperature range 2400{4000 K. Also given is the population for the HPCS in the bulk before raising the temperature. . . . . . . . . . . . . . . . . . . . . . . . . 84 vii 4.3 The nal conguration of the material at 3600 K, obtained by (NVT )-MD simulation. The density is 0.98 g/cm 3 . . . . . . . . . . 85 4.4 Populations of gaseous molecules in the nal congurations and their dependence on temperature. . . . . . . . . . . . . . . . . . . . . . . 87 4.5 Populations of CH n radicals and gaseous molecules in the nal con- gurations and their dependence on temperature. . . . . . . . . . . 88 4.6 Populations of SiH n radicals and gaseous molecules in the nal con- gurations and their dependence on temperature. . . . . . . . . . . 89 4.7 Populations of hydrogen radicals in the nal congurations and their dependence on temperature. . . . . . . . . . . . . . . . . . . . . . . 90 4.8 Time- and temperature-dependence of production of hydrogen rad- ical and hydrogen molecule during (NVT )-MD simulations. . . . . . 91 4.9 Time-dependence of production of CH n , during (NVT )-MD simu- lation at 3200 K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 4.10 Time-dependence of production of CH n , during (NVT )-MD simu- lation at 3600 K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 4.11 Time-dependence of production of CH n , during (NVT )-MD simu- lation at 4000 K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 4.12 Time-dependence of production of SiH n , during (NVT )-MD simu- lation at 3200 K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 4.13 Time-dependence of production of SiH n , during (NVT )-MD simu- lation at 3600 K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 4.14 Time-dependence of production of SiH n , during (NVT )-MD simu- lation at 4000 K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 4.15 Time-dependence of production of CH n Si, during (NVT )-MD sim- ulation at 4000 K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 4.16 Temperature-dependence of rate of production of H 2 . The results for three separate (NVT )-MD simulation runs with distinct initial congurations are shown, in order to provide a measure of the vari- ations. Straight line represents t of the data from trial 1. . . . . . 99 4.17 The schematic of shell script for removing the hydrogen radicals from the system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 viii 4.18 Time-dependence of the pressure and temperature during MD sim- ulations for heating up four HPCS strands at a density of 0.98 g/cm 3 .100 4.19 Time-dependence of the potential energy and pressure during (NVT )- MD simulations. The target temperature is 2000 K. . . . . . . . . . 101 4.20 The increase in the density of the system during the (NPT )-MD simulation to compress the material. . . . . . . . . . . . . . . . . . 102 4.21 The change in the temperature and pressure of the system during the (NPT )-MD simulation to compress the material. The target temperature is 2000 K. . . . . . . . . . . . . . . . . . . . . . . . . . 102 4.22 Time-dependence of the density, T and P during (NVT )-MD sim- ulation to lower the temperature from 2000 K to 300 K. . . . . . . . 103 4.23 The change in the density and pressure during (NPT )-MD simula- tion. The pressure is released from its high value, and temperature is kept at 300 K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 4.24 The change in the density and pressure during (NPT )-MD simula- tion for an extra 60 ps of MD simulation. The temperature is kept at 300 K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 4.25 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. . . . . . 106 4.26 Variations of the properties of the system during further energy minimization, using simulated annealing. . . . . . . . . . . . . . . . 106 4.27 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. . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 4.28 The experimental (top) and the computed (bottom) XRD pattern of the SiC, computed at the end of the nal MD simulation. The experimental data were adopted after Musumeci et al. [2]. . . . . . 108 5.1 The accessible free volume fraction versus time with a probe diam- eter of 2.5 A. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 ix 5.2 Free volume percentage as a function of the probe size for the SiC model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 5.3 The distribution of the free volume versus time. . . . . . . . . . . . 120 5.4 The ensemble of the trajectories of H 2 , CH 4 and CO 2 . The simula- tion cells have been rotated to provide a better view of trajectories. 122 5.5 Displacements of three gases from their last positions. . . . . . . . . 123 5.6 Logarithmic plot of the MSD versus time of the three gases, com- puted in the (NPT ) ensemble. . . . . . . . . . . . . . . . . . . . . . 126 5.7 The MSDs of the gases versus time computed in the (NPT ) ensem- ble at 300 K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 5.8 The MSDs of H 2 versus time, computed in the (NPT ) ensemble at 300 K and 400 K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 6.1 Time-dependance of the separation factor for the H 2 /CO 2 mixture for pressures drops of 20, 40, and 60 atm at 475 K. . . . . . . . . . 134 6.2 Time-dependance of the separation factor for the H 2 /CH 4 mixture for pressures drops of 20, 40, and 60 atm at 475 K. . . . . . . . . . 135 6.3 The dimensionless density proles of the H 2 /CO 2 mixture in the axial (X ) direction, in the membrane (middle) and the two con- trol volumes. The temperature is 375 K and the applied external pressure drop is 20 atm. . . . . . . . . . . . . . . . . . . . . . . . . 137 6.4 The dimensionless density proles of the H 2 /CO 2 mixture in the axial (X ) direction, in the membrane (middle) and the two con- trol volumes. The temperature is 475 K and the applied external pressure drop is 20 atm. . . . . . . . . . . . . . . . . . . . . . . . . 138 6.5 The dimensionless density proles of the H 2 /CO 2 mixture in the axial (X ) direction, in the membrane (middle) and the two con- trol volumes. The temperature is 475 K and the applied external pressure drop is 60 atm. . . . . . . . . . . . . . . . . . . . . . . . . 139 6.6 The dimensionless density proles of the H 2 /CO 2 mixture in the axial (X ) direction, in the membrane (middle) and the two con- trol volumes. The temperature is 475 K and the applied external pressure drop is 100 atm. . . . . . . . . . . . . . . . . . . . . . . . . 140 x 6.7 The dimensionless density proles of the H 2 /CH 4 mixture in the axial (X ) direction, in the membrane (middle) and the two con- trol volumes. The temperature is 375 K and the applied external pressure drop is 20 atm. . . . . . . . . . . . . . . . . . . . . . . . . 141 6.8 The dimensionless density proles of the H 2 /CH 4 mixture in the axial (X ) direction, in the membrane (middle) and the two con- trol volumes. The temperature is 475 K and the applied external pressure drop is 20 atm. . . . . . . . . . . . . . . . . . . . . . . . . 142 6.9 The dimensionless density proles of the H 2 /CH 4 mixture in the axial (X ) direction, in the membrane (middle) and the two con- trol volumes. The temperature is 475 K and the applied external pressure drop is 60 atm. . . . . . . . . . . . . . . . . . . . . . . . . 143 6.10 The dimensionless density proles of the H 2 /CH 4 mixture in the axial (X ) direction, in the membrane (middle) and the two con- trol volumes. The temperature is 475 K and the applied external pressure drop is 100 atm. . . . . . . . . . . . . . . . . . . . . . . . . 144 6.11 The average dimensionless temperature T in the membrane (in the middle) and the two control volumes for the H 2 /CO 2 mixture. Numbers in the parentheses indicate the set temperatures. The applied external pressure drop is 20 atm. . . . . . . . . . . . . . . . 145 6.12 The average dimensionless temperature T in the membrane (in the middle) and the two control volumes for the H 2 /CH 4 mixture. Numbers in the parentheses indicate the set temperatures. The applied external pressure drop is 20 atm. . . . . . . . . . . . . . . . 146 6.13 Eect of the temperature on the separation factor for the two gaseous mixtures. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147 6.14 Eect of the temperature on the uxes for the H 2 /CO 2 gaseous mixtures. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147 6.15 Eect of the temperature on the uxes for the H 2 /CH 4 gaseous mixtures. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 6.16 Eect of the applied external pressure drop P on the separation factors for a membrane of thickness 64.2 A. . . . . . . . . . . . . . . 148 6.17 Eect of the membrane's thickness on its separation factor for the H 2 /CO 2 (top) and H 2 /CH 4 (bottom) mixtures at T = 475 K and P = 20 atm. The data are tted to the S d (X) function. . . . . . . 150 xi 6.18 Eect of the membrane's thickness on its separation factor for the H 2 /CO 2 (top) and H 2 /CH 4 (bottom) mixtures at T = 475 K and P = 30 atm. The data are tted to the S d (X) function. . . . . . . 151 6.19 Eect of the membrane's thickness on its separation factor for the H 2 /CO 2 (top) and H 2 /CH 4 (bottom) mixtures at T = 475 K and P = 40 atm. The data are tted to the S d (X) function. . . . . . . 152 6.20 The tting functions,S d (X) andS l (X), approaching the maximum separation factor for the H 2 /CO 2 (top) and H 2 /CH 4 (bottom) mix- tures at T = 475 K and P = 20 atm. . . . . . . . . . . . . . . . . 153 6.21 The tting functions,S d (X) andS l (X), approaching the maximum separation factor for the H 2 /CO 2 (top) and H 2 /CH 4 (bottom) mix- tures at T = 475 K and P = 30 atm. . . . . . . . . . . . . . . . . 154 6.22 The tting functions,S d (X) andS l (X), approaching the maximum separation factor for the H 2 /CO 2 (top) and H 2 /CH 4 (bottom) mix- tures at T = 475 K and P = 40 atm. . . . . . . . . . . . . . . . . 155 A.1 Distribution of (a) Si{C ; (b) Si{H, and (c) C{H distances in the HPCS polymer, as generated by the PCFF. . . . . . . . . . . . . . 172 A.2 Distribution of (a) Si{C{Si ; (b) Si{C{H, (c) C{Si{H and (d) C{Si{C triplet angles in the HPCS polymer, as generated by the PCFF. . . 173 A.3 Distribution of (a) H{C{Si{C ; (b) H{Si{C{Si, (c) H{C{Si{H and (d) Si{C{Si{C dihedral angles in the HPCS polymer, as generated by the PCFF. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 174 B.1 Temperature-dependence of the rate of production of H 2 in both original and smoothed data (top) at 3000 K for trial 1, tted f(t) against original (middle) and smoothed data (bottom). . . . . . . . 177 B.2 Temperature-dependence of the rate of production of H 2 in both original and smoothed data (top) at 3000 K for trial 2, tted f(t) against original (middle) and smoothed data (bottom). . . . . . . . 178 B.3 Temperature-dependence of the rate of production of H 2 in both original and smoothed data (top) at 3000 K for trial 3, tted f(t) against original (middle) and smoothed data (bottom). . . . . . . . 179 B.4 Temperature-dependence of the rate of production of H 2 in both original and smoothed data (top) at 3200 K for trial 1, tted f(t) against original (middle) and smoothed data (bottom). . . . . . . . 180 xii B.5 Temperature-dependence of the rate of production of H 2 in both original and smoothed data (top) at 3200 K for trial 2, tted f(t) against original (middle) and smoothed data (bottom). . . . . . . . 181 B.6 Temperature-dependence of the rate of production of H 2 in both original and smoothed data (top) at 3200 K for trial 3, tted f(t) against original (middle) and smoothed data (bottom). . . . . . . . 182 B.7 Temperature-dependence of the rate of production of H 2 in both original and smoothed data (top) at 3400 K for trial 1, tted f(t) against original (middle) and smoothed data (bottom). . . . . . . . 183 B.8 Temperature-dependence of the rate of production of H 2 in both original and smoothed data (top) at 3400 K for trial 2, tted f(t) against original (middle) and smoothed data (bottom). . . . . . . . 184 B.9 Temperature-dependence of the rate of production of H 2 in both original and smoothed data (top) at 3400 K for trial 3, tted f(t) against original (middle) and smoothed data (bottom). . . . . . . . 185 B.10 Temperature-dependence of the rate of production of H 2 in both original and smoothed data (top) at 3600 K for trial 1, tted f(t) against original (middle) and smoothed data (bottom). . . . . . . . 186 B.11 Temperature-dependence of the rate of production of H 2 in both original and smoothed data (top) at 3600 K for trial 2, tted f(t) against original (middle) and smoothed data (bottom). . . . . . . . 187 B.12 Temperature-dependence of the rate of production of H 2 in both original and smoothed data (top) at 3600 K for trial 3, tted f(t) against original (middle) and smoothed data (bottom). . . . . . . . 188 B.13 Temperature-dependence of the rate of production of H 2 in both original and smoothed data (top) at 3800 K for trial 1, tted f(t) against original (middle) and smoothed data (bottom). . . . . . . . 189 B.14 Temperature-dependence of the rate of production of H 2 in both original and smoothed data (top) at 3800 K for trial 2, tted f(t) against original (middle) and smoothed data (bottom). . . . . . . . 190 B.15 Temperature-dependence of the rate of production of H 2 in both original and smoothed data (top) at 3800 K for trial 3, tted f(t) against original (middle) and smoothed data (bottom). . . . . . . . 191 xiii B.16 Temperature-dependence of the rate of production of H 2 in both original and smoothed data (top) at 4000 K for trial 1, tted f(t) against original (middle) and smoothed data (bottom). . . . . . . . 192 B.17 Temperature-dependence of the rate of production of H 2 in both original and smoothed data (top) at 4000 K for trial 2, tted f(t) against original (middle) and smoothed data (bottom). . . . . . . . 193 B.18 Temperature-dependence of the rate of production of H 2 in both original and smoothed data (top) at 4000 K for trial 3, tted f(t) against original (middle) and smoothed data (bottom). . . . . . . . 194 C.1 Eect of the membrane's thickness on its separation factor for the H 2 /CO 2 (top) and H 2 /CH 4 (bottom) mixtures at T = 475 K and P = 50 atm. The data are tted to the S d (X) function. . . . . . . 196 C.2 Eect of the membrane's thickness on its separation factor for the H 2 /CO 2 (top) and H 2 /CH 4 (bottom) mixtures at T = 475 K and P = 60 atm. The data are tted to the S d (X) function. . . . . . . 197 C.3 The tting functions,S d (X) andS l (X), approaching the maximum separation factor for the H 2 /CO 2 (top) and H 2 /CH 4 (bottom) mix- tures at T = 475 K and P = 50 atm. . . . . . . . . . . . . . . . . 198 C.4 The tting functions,S d (X) andS l (X), approaching the maximum separation factor for the H 2 /CO 2 (top) and H 2 /CH 4 (bottom) mix- tures at T = 475 K and P = 60 atm. . . . . . . . . . . . . . . . . 199 xiv List of Tables 2.1 General models: HF, Hartree-Fock; MP2, MP3, MP4, Meller- Plesset models; QCI, quadratic conguration interaction [QCISD(T)]; FCI, full conguration interaction. Table courtesy of [3]. . . . . . . 23 2.2 The parameters of the bond order functions. . . . . . . . . . . . . . 28 2.3 Classication of the interactions in ReaxFF. . . . . . . . . . . . . . 31 2.4 The parameters of one-body energy functions. . . . . . . . . . . . . 34 2.5 Values of Tap in the taper function equation. . . . . . . . . . . . . 36 2.6 The parameters of two-body energy functions. . . . . . . . . . . . . 37 2.7 The parameters of three-body energy functions. . . . . . . . . . . . 40 2.8 The parameters of four-body energy functions. . . . . . . . . . . . . 42 2.9 The Lennard-Jones parameters for the various gases and membrane atoms, used in the DCV-GCMD simulations. . . . . . . . . . . . . . 44 2.10 The conversion between the actual and reduced units. Subscript 1 refers to the value of the parameters for component 1, which is CH 4 in this dissertation. Table courtesy of [4]. . . . . . . . . . . . . . . . 44 3.1 Comparison of the reaction energies E (in kcal/mol), computed by the DFT and ReaxFF . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.2 The atom parameters of ReaxFF, computed for the HPCS. The rst three from the right are the EEM parameters, while the rest are the bond-order correction coecients. . . . . . . . . . . . . . . . . . . . 59 3.3 Bond energiesD e (in kcal/mol) and bond-order parameters of ReaxFF, computed for the HPCS. The parameter r is in A. . . . . . . . . . 59 3.4 van der Waal parameters of ReaxFF, computed for the HPCS. . . . 59 xv 3.5 Parameters for the triplets in ReaxFF, computed for the HPCS. The angles are in degrees. . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.6 Torsion parameters of ReaxFF for quadruplets, computed for the HPCS. V i is in kcal/mol. . . . . . . . . . . . . . . . . . . . . . . . . 60 3.7 Values of the force constants K i and the equilibrium lengths l 0 for various pairs, used with the PCFF. . . . . . . . . . . . . . . . . . . 62 3.8 Values of the force constants H i and the equilibrium angles 0 for various triplets, used with the PCFF. . . . . . . . . . . . . . . . . . 62 3.9 Values of the force constants D i for the quadruplets, used with the PCFF. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.10 Values of the size and energy parameters of the non-bond part of the PCFF. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.11 The pyrolysis products after heating up the polymer to 3600 K. N is the number of the produced product, and MW its molecular weight. 72 5.1 The computed self-diusivities of H 2 , CH 4 and CO 2 . All numbers are in m 2 /s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 6.1 The predicted values of the membrane thickness for the H 2 /CO 2 and H 2 /CH 4 mixtures for the pressure drops P = 20, 30, and 40 atm. The temperature is 475 K. . . . . . . . . . . . . . . . . . . . 149 xvi Abstract A broad class of important materials, such as silicon carbide (SiC), are fabricated by temperature-controlled pyrolysis of pre-ceramic polymers. In particular, fabri- cation of SiC membranes by pyrolysis of a polymer precursor that contains Si is quite attractive for separation of hydrogen from other gases. It has been quite dif- cult to extract atomistic-scale information about such SiC membranes, since they are amorphous. The research presented in this dissertation extends the ReaxFF reactive force eld to describe the processes involved in the thermal decomposition of hydridopolycarbosilane (HPCS) to form SiC nanoporous membranes. First, we carry out quantum mechanical calculations on models meant to cap- ture the important reaction steps and structures. Then, we develop a model of the HPCS polymer and utilize ReaxFF to describe the thermal degradation and decomposition of the polymer as the system is heated by molecular dynamics (MD) simulations. In the next step, we use ReaxFF for reactive dynamics of HPCS over a wide range of temperature. We then simulate pyrolysis of the HPCS under conditions that closely mimic the conditions of the fabrication of nanoporous SiC membranes to produce amorphous SiC material. The pyrolysis results and the computed properties of the SiC ceramic are shown to be in good agreement with the experimental data. xvii To further test the validity of the model and to provide insight into improving the performance of the membrane, extensive MD simulations were carried out to compute the self-diusivities of H 2 , CO 2 and CH 4 in the molecular model of SiC. The results indicate higher values of the diusivity for H 2 . The morphology of the amorphous SiC layer is characterized by computing its accessible free volumes and the cavity distributions. Finally, we use the molecular model of SiC and non-equilibrium MD simula- tions in order to study transport and separation in the membrane of two binary gaseous mixtures, namely, H 2 /CO 2 and H 2 /CH 4 at various temperatures and pres- sure drops, applied externally to the membrane. When compared with our own experimental data, the model is demonstrated to provide accurate predictions for various properties of interest and, in particular, for the separation factors of the mixtures. The model can be used to determine the optimal membrane's thickness. xviii Chapter 1 Introduction 1.1 Outlook Many important materials are fabricated by temperature-controlled pyrolysis of pre-ceramic polymers. Examples include carbon molecular sieves, silicon carbide (SiC), silicon nitride (Si 3 N 4 ), and aluminumnitride (AlN) [5{9]. Such materials have important applications in the design and synthesis of fuels and chemicals as separation and reactive media, such as adsorbents, membranes and catalysts, as well as various types of sensors. The fabrication of such materials involves a series of tailored reactions via which the original polymer precursor is converted into the nal ceramic. Consider, for example, fabrication of nanoporous membranes. Growing inter- est in the hydrogen economy has motivated research on inorganic, hydrogen- permselective membranes for use in the processes related to H 2 production that take place at high temperatures and pressures. Due to its many unique properties, such as high thermal conductivity [10], thermal shock resistance [11], biocompat- ibility [12], resistance to acidic and alkali environments [13], chemical inertness, and high mechanical strength [14, 15], a promising material for fabricating such inorganic membranes is SiC. This chapter is organized as follows. In Sections 1.2 and 1.3 we review, respec- tively, the experimental methods for the fabrication of SiC membranes and the 1 mechanism of the degradation of the polymer precursor as a function of temper- ature. Section 1.4 describes the importance, as well as the approaches for the modelings and simulations that are used throughout this dissertation. The last section describes the plan of this dissertation. 1.2 Fabrication Methods of Silicon Carbide (SiC) Membranes Two distinct approaches have been developed for fabrication of SiC nanoporous membranes. One involves the use of chemical-vapor deposition (CVD)/chemical- vapor inltration (CVI) techniques, while the second approach is based on a dip- coating technique combined with the pyrolysis of a polymer precursor. Both approaches have been pursued by our group at USC [5{8, 16, 17], but the rst one suers from certain diculties [16,17] that make it less attractive. The second method utilizes the pyrolysis of a polymer precursor, hydri- dopolycarbosilane (HPCS), or a partially allyl-substituted hydridopolycarbosilane (AHPCS). The role of the allyl group is to slightly increase the carbon content of the nal ceramic, and enable the polymer to cross link at lower temperatures. The AHPCS is curable in the presence of argon, rather than oxygen that introduces Si{O{C bonds in the polycarbosilane (PCS)-derived ceramics that are thermally and hydrothermally unstable. The primary pyrolysis product of both the HPCS and AHPCS is a SiC ceramic with Si{C bonds, as well as some Si{Si and C{C ones (see below). To fabricate SiC nanoporous membrane using the dip-coating technique and the pyrolysis of the AHPCS [5{8], porous SiC support tubes are prepared using uniaxial cold-pressing of -SiC powder, together with the appropriate sintering aids. The 2 pores of the support are, however, too large to be eective for separating gaseous (and even liquid) mixtures. To endow the support with high selectivity for the separation of uid mixtures, one deposits on the support (via dip-coating) a thin layer of the pre-ceramic polymer precursor, which leaves behind a thin nanoporous SiC lm after pyrolysis. Figure 1.1 represents a schematic of SiC fabrication process with a scanning electron microscope (SEM) picture of the SiC lm on the surface of the support. During pyrolysis bonds in the polymer break, releasing various atoms and radicals that, in turn, react with each other and produce such gases as H 2 , CH 4 and SiH 4 . The dip-coating step is repeated a number of times, in order to \repair" some of the cracks and imperfections that develop during the formation of the ceramic layer. Although pyrolysis of the AHPCS is supposed to result in a SiC ceramic with a near stoichiometric Si:C ratio, one cannot exclude the possibility that there may remain a small amount of carbon after the pyrolysis, which may in uence the membrane's performance. 1.3 Pyrolysis Mechanism and Film Growth Understanding equilibrium and non-equilibrium phenomena, such as adsorption, reaction, and transport of uid mixtures in SiC membranes entails having an accu- rate model of the materials' pore space, which can be utilized for designing the optimal pore structure for the separation process. One of the most important ingre- dients for developing such a model is a fundamental understanding of the pyrolysis of the polymer precursors, such as the HPCS and AHPCS, and of the evolution of the structure of the thin lm deposited on the support. It is, therefore, instructive to brie y review the experimental results for pyrolysis of the PCS to SiC ceramics, 3 Figure 1.1: Schematic illustration of SiC fabrication process with a SEM picture of thin SiC lm on the surface of the support. as they will be compared with the results of the molecular simulations presented in following chapters. Hasegawa and Okamura [18] studied pyrolysis of several types of the PCS using various techniques, and suggested six stages for their conversion to SiC ceramic. In the rst stage, up to about 625 K, gas generation is very small and degra- dation of low-molecular weight molecules occurs. The average molecular weight of the PCS begins to change during the second stage from 625 K to 825 K, but the gas production is still small. Dehydrogenation of the Si{H bonds occurs to some extent and, due to cross-linking, the Si{Si bonds also appear and some H 2 is produced. Depending on the type of the PCS, some alkanes may form that are the result of dehydrocarbonation condensation. In the third stage from 825 K to about 1075 K the Si{H and C{H bonds break readily, resulting in the formation of many free hydrogen radicals that bond together and generate H 2 . For the PCSs 4 that have methyl group as side chains, CH 4 may also be produced as a result of the decomposition of the methyl group and the subsequent reaction with free hydro- gen radicals. With increasing temperature the connectivity of the Si{C{Si bonds increases as a result of dehydrogenation of the Si{H and C{H bonds and demetha- nation of SiCH 3 , producing eventually an inorganic amorphous residue. At about 975 K the residue includes C=C bonds that are produced by dehydrogenation of the {CH{CH bonds [19,20], as well as hydrogen atoms. The fourth and fth stages happen between 1075 K and 1475 K, with the former known as a transition stage in which the structure is still amorphous. The remaining hydrogen atoms decompose with increasing temperature to 1275 K, above which nucleation of Si{C begins and the amorphous SiC is transformed into a crystalline structure. The carbon content of the system in this stage is more than its Si content and results in a nonstoichiometric crystalline structure with crystal sizes of 2{3 nm. Finally, above 1475 K crystal growth occurs in which a stoichiometric ratio between Si and C is obtained by burning the extra carbon. A schematic representation that summarizes the eect of the temperature change on the pyrolysis products of the PCS polymer precursor is given in Figure 1.2. Tang et al. [21] pyrolyzed the PCS in H 2 atmosphere, and reported the initiation of the process by cleavage of the Si{H bonds below 825 K, and of the Si{CH 3 bonds that break signicantly above 1175 K. The properties of SiC membranes and bers at high temperatures are aected negatively by the presence of extra Si. Therefore, Tang et al. used a hydrogen atmosphere to achieve a Si/C ratio close to one in the nal ceramic ber via reaction of the free radicals with H 2 and generation of CH 4 . Interrante and Moraes [22] studied the pyrolysis chemistry of the AHPCS and a hyperbranched polycarbosilane. The results indicated that H 2 was the main product of the pyrolysis whose loss caused the polymers to cross-link above 525 5 Figure 1.2: Schematic illustration of the eect of temperature change on the pyrol- ysis products of the polycarbosilane precursor. K. Amorphous SiC was produced after the polymer was heated up to 1275 K, which was then transformed into a nanocrystalline SiC at about 1875 K. A quite similar process was reported by Interrante et al. [23] for the pyrolysis of the HPCS. Its related and compositionally equivalent linear polymer analog, polysilaethylene (PSE), [SiH 2 CH 2 ] n , was studied by Liu et al. [24]. Both polymers produced H 2 above 575 K by a (1,1)-elimination process to form a cross-linked network. The process, also called unimolecular elimination, is a model for explaining a particular type of chemical elimination reaction that consists of ionization and deprotonation. More H 2 and an amorphous SiC were produced between 725 K and 1275 K. Mass spectrometric analysis of the gas production during the pyrolysis of the PCS, such as the HPCS, indicated that although most of the hydrogen production occurs below 1075 K, a signicant amount is also produced above this temperature [24]. The loss of hydrogen results in further condensation of the amorphous material. 6 A reaction pathway was postulated for cross-linking and pyrolysis of the HPCS in which both (1,1)-H 2 elimination and intramolecular H-transfer reactions lead to highly reactive silylene intermediates, producing eventually Si{Si bonds that rapidly rearrange to Si{C at such temperatures and form Si{C interchain cross- links. 1.4 Modeling and Simulation To obtain a fundamental understanding of such pyrolysis processes we carried out detailed atomistic-scale simulations. In principle, ab-initio quantum mechanics (QM) can provide information about the structure of the amorphous systems. However, to determine the structure of the SiC membrane layer one should capture in the simulations the various reactive processes involved in forming the layer. This requires QM simulations on systems with about 3000 atoms per cell at temperature of 1200 K for microseconds, which are far beyond the current QM capabilities. Instead, this dissertation extends the ReaxFF reactive force eld, validated for high temperature reactions of other materials, to describe the processes involved in the thermal decomposition of HPCS to form SiC nanoporous membranes. First, we carry out QM calculations on models meant to capture important reaction steps and structures. Then, we use ReaxFF to carry out reactive molecular dynamics (MD) simulation of the pyrolysis processes. There are currently several MD methods that can simulate reactive systems. These are reviewed by Farah et al. [25]. Three of the better known of such methods are as follows. ReaxFF reactive force eld was developed by van Duin, Goddard and coworker [26], and applied to a large number of reactions in metallic, ceramic, semiconductor, and polymer systems that we study here. Another method is the 7 reactive empirical bond-order potentials developed by Brenner and coworkers [27] and applied to carbon materials. Nyden and Noid [28] used a conventional non- reactive FF, but dynamically varied it to describe chemical reactions. Chemical reactions can also be simulated by the tight binding and density functional meth- ods, but the focus of this dissertation is on MD simulation of chemical reactions. ReaxFF bridges the gap between the computational methods based on quan- tum chemical (QC) and empirical force elds. Although the QC methods are, in principle, applicable to all materials, their computational expense makes them inapplicable for large systems (with more than a few hundred atoms). ReaxFF, on the other hand, enables one to carry out MD simulations of large-scale reactive chemical systems with accuracy near QC at a fraction of the cost. It also allows for accurate description of bond breaking and bond formation and, therefore, it may be used in the MD simulation of pyrolysis of a polymer. The work in this dissertation develops the ReaxFF aimed at simulating pyrol- ysis of the HPCS. ReaxFF has been previously developed for hydrocarbons and hydrocarbon oxidation [26, 29], metaloxides [30] and metal hydrides [31], Si/Si oxides [32], aluminum/aluminum oxides [33], nitramines [34], and a variety of complex chemical systems [35{38]. There has also been a number of studies for simulating pyrolysis and thermal degradation of polymers using ReaxFF. Chenoweth et al. [38] used ReaxFF to investigate the thermal decomposition of poly(dimethylsiloxane). Desai et al. [39] studied the pyrolysis of phenolic polymer in the presence of carbon nanotubes, while Saha and Schatz studied [40] thermal decomposition of polyacrylonitrile in the carbonization process. As described earlier, the pyrolysis process is an essential step in the fabrication of nanoporous SiC membranes. Therefore, we utilize ReaxFF with reactive dynam- ics (RD) simulation to study the temperature-dependence of the decomposition of 8 the HPCS. We carry out RD simulation of the pyrolysis process under the condi- tions 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. To further test the validity of the model, we carry out extensive MD simulations to compute the self diusivities of H 2 , CO 2 and CH 4 in amorphous SiC generated by inserting the gases in the SiC. Atomistic models of the composites are gen- erated using energy minimizations, MD simulations, and the polymer-consistent force eld. The morphology of the SiC model is characterized by its density, the radial distribution function, and the accessible free volume, which we compute with probe molecules of dierent sizes. The distributions of the cavity volumes are computed using the Voronoi tessellation method. Additionally, we utilize the generated amorphous SiC model in non-equilibrium MD simulations in order to study transport and separation in the membrane of two binary gaseous mixtures, namely, H 2 /CO 2 and H 2 /CH 4 , and test the accuracy of the results by comparing them with the experimental data. The eect of the temperature and the pressure drop applied externally to the membrane is described. Computations for the work described in this dissertation were supported by the University of Southern California Center for High-Performance Computing and Communications (hpcc.usc.edu), San Diego Supercomputer Center at the Uni- versity of San Diego (sdsc.edu), and Material and Process Simulation Center at California Institute of Technology (wag.caltech.edu). 9 1.5 The Plan of the Dissertation The rest of this dissertation is organized as follows. In the next chapter we pro- vide the reader with detailed information and references about the techniques and tools that are used for carrying out the simulations throughout this work. They include quantum-chemistry (QC), MD simulation, details of ReaxFF, and dual control volume Grand-Canonical MD (DCV-GCMD) simulation technique. Chap- ter 3 explains the QC computations and force eld optimization for estimating the parameters of ReaxFF for the pyrolysis of HPCS. We also describe in that chapter the generation of a molecular model of the HPCS, as well as the results of MD simulation of thermal decomposition of the polymer. In Chapter 4, we describe how the polymer model is utilized in the MD simulations of its pyrolysis. Molec- ular dynamics simulation of thermal decomposition of the HPCS is described in that chapter, where the results are presented, discussed, and compared with the experimental data. 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 also described in that chapter. The chapter ends with comparison of the properties of the resulting SiC ceramic with the experimental data. In Chapter 5, the self-diusivity of various gases and physical properties of the SiC material are described. Chapter 6 discusses the non-equilibrium MD simu- lations of gas mixtures transport within the generated SiC layer. The dissertation is summarized in the last chapter. 10 Chapter 2 Background and Computational Methods 2.1 Introduction In this chapter, we provide detailed information and references about the tech- niques and tools that are used for carrying out the simulations throughout this work. Section 2.2 explains the essentials of the quantum chemistry with the focus on the topics that are related to this work. The materials for that section were taken carefully from dierent text books and journal papers to clearly explain the concepts [3,41{43]. Section 2.3 explains very brie y about molecular dynamics (MD) and Monte-Carlo (MC) simulation and provides references for further details on these topics. In Section 2.4, details of ReaxFF reactive force eld are discussed. Section 2.5 explains the dual control volume Grand-Canonical MD (DCV-GCMD) simulation technique. 2.2 Quantum Chemistry The concepts of quantum mechanics (QM) and quantum chemistry (QC) have been discussed, from elementary to advanced levels, by many quantum chemists [3,41{ 43]. In this section we provide a short overview of the QC with focus on the methods and techniques that will be used in the following chapters. 11 The goal of QC is to provide an accurate explanation of chemical processes at its most fundamental level. Erwin Schr odinger proposed the basic nonrelativistic wave equation governing the motion of nuclei and electrons in molecules, which provides the theoretical foundation for computational chemistry, ^ H =E ; (2.1) where is the wavefunction of a particular state, which is a function of the spin and cartesian coordinates of the included particles, ^ H is the Hamiltonian operator associated with the observable energy, and E is a scalar that represents the total energy of the system [41{43]. The Hamiltonian of the system ( ^ H) is an operator that contains all the required terms that contribute to the energy of the system, ^ H = ^ T + ^ V ; (2.2) where ^ T and ^ V are the kinetic and potential energy operators, respectively. The contributions of electron and nuclei are included in each of these terms, ^ T = ^ T e + ^ T n ; (2.3) ^ V = ^ V nn + ^ V ne + ^ V ee ; (2.4) ^ T e = X i h 2 2m e r 2 r i ; (2.5) ^ T n = X i h 2 2M i r 2 R i ; (2.6) ^ V nn = X i X j>i Z i Z j e 2 4 0 jR i R j j ; (2.7) 12 ^ V ne = X i X j Z i e 2 4 0 jR i r j j ; (2.8) ^ V ee = X i X j>i e 2 4 0 jr i r j j ; (2.9) where M i is the mass of the nucleus i, Z i is the atomic number of nucleus i, m e is the mass of the electron,R i is the nuclear coordinate,r i is the electron coordinate, and h is the reduced Planck constant. As an example, consider the hydrogen atom. Schr odinger's equation can be easily solved by considering the nucleus as the center of the coordinates and, therefore, the Hamiltonian is given by, ^ H = 1 2 r 2 r 1 r : (2.10) Note that the wavefunction is simply a function ofr. The results for hydrogen atom are found to be in perfect agreement with experimental spectroscopic data [41{43]. In addition to the total energy of the system, many chemical properties can be obtained from derivatives of the energy with respect to some external param- eters. External parameters include geometric parameters (bond lengths, angles, dihedrals), external electric eld, external magnetic eld (NMR experiments), and so on. Examples of computed properties contain bond energies, reaction energy, structures of ground-, excited- and transition-states, atomic charges and electro- static potentials, vibrational frequencies (IR and Raman), transition energies and intensities for UV and IR spectra, NMR chemical shifts, dipole moments, polar- isabilities, hyperpolarisabilities and reaction pathways and mechanisms. Most of such properties are needed to be computed in this dissertation in order to develop the force eld ReaxFF. However, the exact solution of the Schr odinger's equation exists only for a small number of systems, such as the rigid rotor, the harmonic oscillator, a particle in a box and the hydrogenic ions (H, He + , Li 2+ ,...), but for 13 many other systems, even for a simple molecule such as H + 2 that consists of three particles, the Schr odinger's equation cannot be solved analytically, hence posing a challenge to the quantum chemists for many years. In order to solve the problem, approximate mathematical models must be devel- oped. Although the theoretical models are not exact, they should be well-dened and accurate enough to assist the qualitative interpretation of chemical phenom- ena and provide predictive capability. In the development of such models one must take into account the target accuracy, precise formulation of approximate mathe- matical models, reasonable time and cost of implementation, verication against known chemical facts and, nally, predictions of the model by applying to unknown chemical problems [3]. The Born-Oppenheimer (BO) approximation is often adopted that takes care of the great dierence in masses of nuclei and electrons. Nuclei travel much more slowly than electrons because of being much heavier (the mass of a proton 2000 times that of an electron) and, hence to a good approximation, one can consider the electrons in a molecule to be moving in the eld of xed nuclei. In other words, the electrons can respond almost instantaneously to the displacement of the nuclei. Therefore, the wavefunction of a molecule is broken into its electronic and nuclear (vibrational, rotational) components: total (R;r) = nuclei (R) electron (r;R); (2.11) where the \; " notation indicates a parametric dependence. The potential energy surface is a direct consequence of the BO approximation. The focus is the solution of the electronic Schr odinger's equation, H (R;r) =E(R) (r;R); (2.12) 14 where R is a set of xed locations of the nuclei and E(R) is the electronic energy. The Hamiltonian is given by, H = h 2 2m e X i r 2 i X i X j Z i e 2 4 0 jR i r j j + X i X j>i e 2 4 0 jr i r j j : (2.13) It is conventional not to write the Hamiltonian with the nucleus-nucleus repulsion term, but at the end of calculations it is added as a classical term [41{43]. There are two approaches to solving the Schr odinger's equation. One approach is based on ab-initio calculation (the term comes from the Latin word for \from the beginning"), where models utilize only the fundamental constants of physics. The other approach is semi-empirical that uses a simplied version of the Hamilto- nian, and experimental data are used to obtain the adjustable parameters. Accu- rate ab-initio calculations for large molecules are computationally expensive and, therefore, semi-empirical methods have been developed to deal with a larger range of chemical species [41,42]. 2.2.1 The Hartree-Fock approximation Development of approximate solutions to the electronic Schr odinger's equation has been a major preoccupation of quantum chemists since the start of quantum mechanics. Quantum chemists are faced with many-electron problems, except for the very simple cases like H + 2 . The critical problem is the electron-electron potential energy, which depends on electron-electron separationr ij as given by the last term of the Equation 2.13. The Hartree-Fock (HF) approximation has been a solution to such problems [41{43]. 15 One should take into account the contribution of spin orbitals. For a given molecule, if one assigns the 2n electrons to a set of n molecular orbitals i (i = 1; ;n), the related many electron wave function are given as = (n!) (1=2) det a (1) a (2) a (3) z (n) : (2.14) where and are spin functions and i are considered to be orthonormal. Equa- tion 2.14 is called a Slater determinant. The energy E is evaluated as the expec- tation value of the of the full Hamiltonian, E =h jHj i : (2.15) If the energy is minimized by varying i , then it is completely dened and, accord- ing to the variational principle, is an upper bound for the exact Schr odinger energy (Equation 2.1). Using this methodology, Fock derived a set of related dierential equations for i . Due to some early applications of the method on atoms by Hartree, the method is known as the Hartree-Fock theory [3]. 2.2.2 The Roothaan equations The HF procedure for atoms and their spherical symmetry is relatively straight- forward, which means that HF equations can be solved numerically for the spinor- bitals, but such solution for molecular systems is suciently computationally complex that the technique is required to be improved. A major improvement was proposed by Roothaan [44]. Roothaan introduced a set of N basis function 16 ( = 1; 2; ;N; N > n) and expressed each spatial wavefunction i as a linear combination of these basis functions, i = N X =1 c i ; (2.16) where c i are as yet unknown coecients. The N linearly-independent spatial wavefunctions can be obtained from a set of N basis functions. Therefore, the problem of solving the wavefunctions has been transferred to one of computing the coecients. This leads to a set of equations that can be written in matrix form, FC = SCE ; (2.17) where F =H + X P [(j) (j)=2]; (2.18) H = Z H d ; (2.19) S = Z d ; (2.20) E ij =" i ij ; (2.21) P = 2 n X i c i c i ; (2.22) (j) = Z Z (1) (1) (1=r 12 ) (2) (2)d 1 d 2 ; (2.23) where F , S, and P are, respectively, the Fock, overlap and density matrices. H is the core Hamiltonian, describing the movement of an electron in the empty eld of the nuclei. The eigenvalues " i are the Fock energies for one-electron, the lowest n belonging to the occupied molecular orbitals 1; 2; ;n [3]. 17 If the given functions are exclusively determined according to the nuclear positions, then a full mathematical model is provided by the above nonlinear equa- tions. They are often called self-consistent eld (SCF) equations. In the primitive versions of the molecular orbital theory, the were selected to be the atomic orbitals of the included atoms, with the method known as LCAOSCF for \linear combination of atomic orbitals." The setf g is usually known as the basis set. The basis functions are normally chosen to center at the nuclei and to depend solely on the atomic numbers of that nucleus [3]. In principle, to represent spinorbitals exactly, a complete set of basis functions must be used. However, this is not feasible because of being computationally expensive. Therefore, a nite basis set is always used and an error is computed due to approximation of the basis set which is called basis-set truncation error. The basis sets are usually real functions, and there are probably as many basis sets dened as there are quantum chemists. The choice of a basis set is not nearly random as it may rst appear. The best practice to make a basis set is to have a low number of basis functions (to minimize the number of evaluation of two- electron integrals) and to select them wisely (to minimize the computational eort for the evaluation of each integral), while achieving a small basis-set truncation error [41]. One choice of the basis functions are the Slater-type orbitals (STOs) [45], which are a set of functions that decay exponentially with distance from the nuclei. The introduction of Gaussian-type orbitals (GTOs) by Boys [46] made the ab-initio calculations computationally feasible. Because it is easier to calculate the overlap 18 and other integrals with Gaussian basis functions, this led to huge computational savings. The Cartesian Gaussians are functions of the form, ijk (r 1 r c ) = (x 1 x c ) i (y 1 y c ) j (z 1 z c ) k exp jr 1 rcj 2 ; (2.24) where (x c ;y c ;z c ) are the Cartesian coordinates of the center of the Gaussian func- tion at r c , (x 1 ;y 1 ;z 1 ) are the Cartesian coordinates of an electron at r 1 , is a positive exponent and i, j, and k are non-negative integers. When i =j =k = 0, i +j +k = 1, and i +j +k = 2, the Cartesian Gaussians are s-type Gaussian, p-type Gaussian, and d-type Gaussian, respectively [41]. Despite the fact that the computations are feasible with Gaussian-type orbitals (GTO), it gives a poorer representation of the orbitals at the atomic nuclei and, therefore, a larger basis set must be used to reach accuracies close to the STOs. The problem is alleviated by grouping several GTOs to form what are known as contracted Gaussian functions, j = X i d ji g i ; (2.25) where is contracted Gaussian function, d is the contraction coecient, and g is the primitive Gaussian function. The spatial orbitals are then shown as a linear combination of the contracted Gaussians, i = X j c ji j : (2.26) Today, there are hundreds of basis sets composed of the GTOs. The smallest of these are called minimal basis set in which one function is used to represent each of the orbital elementary valance theory [41]. 19 There are many ways to construct contracted Gaussian basis sets. One approach is using a least-squares t of N primitive Gaussian to a set of STOs. Using this approach, the earlier full Slater results have been reproduced. A com- mon choice is N = 3, which proved to give satisfactory results and is referred as STO-3G basis and the general theoretical model HF/STO-3G [47]. However, the minimal HF/STO-3G model exhibits very high stability of single bonds, when compared with multiple bonds. In the case of anisotropic atoms, this can be con- sidered as the failure of the minimal basis. This problem was solved by using two basis functions per valence atomic orbital, instead of one, which is called 6-31G basis. For the inner shell, it has a single contracted six-Gaussian basis function, while for the valence shell of each atom a set of outer uncontracted and a set of inner three-contracted Gaussians are used. This is a sample of a split-valence basis. There are also other similar basis set such as double-zeta (DZ) and triple-zeta (TZ) basis sets, in which there are, respectively, two and three basis functions, instead of each basis function in the minimal basis set [3,41]. Split-valence bases have several considerable failures that require careful atten- tion. First, such bases prefer structures with high symmetry. Second, it emphasizes too much on the polarity using electric dipole moments, which can be due to lim- itation of lone pair to pure p-orbital types [3]. The HF models were considerably improved by the introduction of 6-31G , also indicated as 6-31G(d). Here, the polarization functions in the form of six d-type functions for each atom (except hydrogen) are added to the 6-31G basis [3,48,49]. Addition of the polarization func- tions in the form of three p-type functions for each hydrogen atom would result in another star, an additional polarization function, 6-31G or 6-31G(d,p) [3]. The basis set can be extended by adding higher polarization functions. The 6- 31G(2df;p) consists ofd functions (two sets),f functions (one set) on heavy atoms, 20 and p functions (one set) on hydrogen [3]. In addition, diuse functions are used that are especially important for anions and electronic states, which are denoted by a \+" as in 6-311+G(d) [3]. The Hartree-Fock models HF/6-31G and 6-31G have been shown to be very eective in the description of molecular conforma- tions [3]. These models have been used frequently in this dissertation for the quantum chemistry computations of dierent properties of the polymeric system. Their overall performance and other aspects have been documented elsewhere [50]. 2.2.3 Electron correlation The HF ground-state wavefunction is still not the \exact" wavefunction, regardless of how accurate it can be. The method works with the averages and the major problem with it is neglecting the electron correlation between the motion of elec- tron in antiparallel spin ( correlation). This can result in considerable errors in, for example, bond dissociation energies. These wavefunctions can be improved by taking into account the electron correlation. The procedure is to start with the HF determinant and form linear combinations with other determinants. Usually, it is done by forming extra determinants from the virtual molecular orbitals (unoc- cupied), which are the higher eigenfunctions of the Fock operator. In general, a basis withN Cartesian functions and 2n electron hasNn virtual orbitals, which may be occupied by or electrons. A general wavefunction with combination of several determinants can be written as =a 0 0 + X ia a a i a i + X ijab a ab ij ab ij +:::; (2.27) 21 wherei;j;k;::: are occupied spinorbitals; a;b;::: are virtual, 0 , a i and ab ij are HF, singly substituted and doubly substituted determinants, respectively. This method is called the conguration interaction (CI) [3,41]. There are also other methods of incorporating electron correlation [3], such as perturbation theory, where a perturbed Hamiltonian is dened as H() =F 0 +fHF 0 g ; (2.28) where F 0 is the Fock Hamiltonian. The computed energy is expressed in powers of , E() =E 0 +E 1 + 2 E 2 + 3 E 3 +::: ; (2.29) where at some level the series is cut and then set = 1. This method is known as Meller-Plesset many-body perturbation theory [51], and it is often indicated by MPn if nished at order n. Another method for correlation theory is the coupled- cluster method whose details can be found elsewhere [52{54]. 2.2.4 General energy models The eort over the past decades has been focused on developing models that can produced results in terms of dierent chemical properties (for example bond and angle energy) close to experimental results. To achieve this goal appropriate basis set and correlation must be used. Table 2.1 gives a summary of some of the options. On the top row dierent correlation methods are given, for which the sophistication increases from left to right. The most left column shows the basis sets that becomes more exible from top to bottom. Full conguration interaction (FCI) considers complete solution within the nite space dened by the basis. The results of applying the complete basis set (only theoretically) is given at the bottom 22 Table 2.1: General models: HF, Hartree-Fock; MP2, MP3, MP4, Meller-Plesset models; QCI, quadratic conguration interaction [QCISD(T)]; FCI, full congura- tion interaction. Table courtesy of [3]. Basis HF MP2 MP3 MP4 QCI FCI STO-3G 6-31G 6-31G* 6-31+G* 6-311+G* 6-311+G(2df) 1 S-eqn of the table. Application of a complete basis set with full conguration interaction to complete solution of the Schr odinger equation is given at the bottom right. The empty boxes in the table describe a well-dened size-consistent theoretical model. The procedure is to start from top-left of the table and then move on to the bottom-right until acceptable agreement between theory and experiment is obtained. Some comparisons between dierent basis and correlation are reviewed by others [3]. 2.2.5 Density functional theory So far and in the ab-initio methods the computations start with the HF approx- imation, where the HF equations are rst solved to construct conguration state functions by nding spinorbitals. These methods suer from computational dif- culty when performing calculations with large basis sets on molecules that con- tain many electrons and atoms. An alternative for this, which is widely used by quantum chemists, is the density functional theory (DFT), which begins with the concept of electron probability density, in contrast to the previous methods that use CSFs. The advantage of the DFT method is that it is less computationally expensive, while it takes into account the electron correlation. It can be used to 23 do calculations on molecules of 100 or more atoms in much less time than the HF methods [41]. The main idea is to write the energy of the system in terms of the electron probability density (r), which represents the total electron density of a system withn electron at a special point r in space. The electronic energy is also dened as a functional of (denoted byE []), which means that for a given function(r) there is a single corresponding energy [41]. For example, consider the Kohn-Sham (KS) equation [55], E [] = h 2 2m e n X i=1 Z i (r 1 )r 2 1 i (r 1 )dr 1 j 0 N X I=1 Z I r I1 (r 1 )dr 1 + 1 2 j 0 Z (r 1 )(r 2 ) r 12 dr 1 dr 2 +E XC [] ; (2.30) where the one electron spatial orbitals i (i = 1; 2;:::;n) are the KS orbitals. The exact ground-state electron density is given by (r) = n X i j i (r)j 2 : (2.31) Once the KS orbitals are computed (see below), is also known. In Equation 2.30, the rst term is the kinetic energy, the second is electron-nucleus attraction (N is the number of nuclei with indexI and atomic numberZ I ), the third is the Coulomb interaction, while the last term is the exchange-correlation energy of the system which is also a function of the density and considers all the non-classical electron- electron interactions. The analytical form of E XC is not known and approximate 24 forms are used. The KS orbitals are determined by solving the KS equations. For one electron orbitals the KS equation is ( h 2 2m e r 2 i j 0 N X I=1 Z I r I1 +j 0 Z (r 2 ) r 12 dr 2 +V XC (r 1 ) ) i (r 1 ) =" i i (r 1 ); (2.32) where " i are the orbital energies of KS, and the V XC is the exchange-correlation potential which is related to exchange-correlation energy by V XC [] = E XC [] : (2.33) By making an approximation for E XC , the V XC and nally the KS equations can be solved [41]. The simplest E XC approximation is known as the local density approximation (LDA) [55], which is mostly used for solids. Another approxima- tion is generalized gradient approximations (GGAs) [56], useful for chemical cal- culations. Hybrids [57] is another approximation that replaces a fraction of GGA exchange [58] with the HF exchange, which leads to the most popular and widely used B3LYP [59] in chemistry. This is used in this dissertation. There is also the PBE GGA approximation [60], which is probably the most popular approximation in materials. These three approximation methods (LDA, PBE, and B3LYP) are known as the standard approximations, because of being used dominantly, compared with other approximation methods [61]. 25 2.3 Molecular Dynamics and Monte Carlo Sim- ulations Molecular dynamics (MD) can be described as brute force integration of Newton's equations of motion to obtain the time evolution of a set of interacting particles. Given the positions, masses and velocities of all the particles in the system and the forces on the particles, the motion of all (individual) particles can be followed in time by calculating the (deterministic) single particle trajectories. On the other hand, the Monte Carlo (MC) method can be described as repeated random sam- pling to obtain numerical results. In this method, the simulations are run many times in order to heuristically compute the similar probabilities. The details of the MD simulations under dierent ensembles, as well as the MC simulations are extensively discussed elsewhere [62{67] and, therefore, due to limited size of this dissertation are not repeated here. Instead, the author would like to provide more details about the background and technical steps of other concepts, which represent the central core of this work. They are given in previous (Quantum-Chemistry) and also subsequent (ReaxFF and GCMD) sections of this chapter. 2.4 ReaxFF: A Reactive Force Field The word \classical" in MD simulations implies that each particle in the system experiences a force and its motion is determined according to Newton's second law. However, one can adopt the eective forces to take into account some quantum mechanical properties of the system and then use them with MD simulations. Clearly, the forces constructed in this way are system dependent and should be tuned (or even reconstructed) when considering new systems. ReaxFF is one of 26 such force eld and is being widely used today (see section 1.4). In ReaxFF, the energy of a system is calculated based on atomic coordinates, bond order (BO) between atomic pairs, and atomic charges. Most of the materials in this section were adopted from the References [1,68]. 2.4.1 Bond order calculation In ReaxFF, the nature of covalent bonds between two atoms i and j is described by a bond order BO ij . There are two steps in calculating the bond order. In the rst step, the BO 0 ij , referred as the raw bond orders, are calculated as a function of interatomic distances, BO 0 ij =BO 0 ij +BO 0 ij +BO 0 ij = exp p bo1 r ij r ij p bo2 + exp p bo3 r ij r ij p bo4 + exp p bo5 r ij r ij p bo6 ; (2.34) where r i is the position of atom i and r ij =jr i r j j is the interatomic distance between atoms i and j. In the second step, the BO 0 ij are improved to provide better atomic valences as BO ij =BO ij +BO ij +BO ij ; (2.35) BO ij =BO 0 ij f 1 ( 0 i ; 0 j )f 4 ( 0 i ;BO 0 ij )f 5 ( 0 j ;BO 0 ij ); (2.36) BO ij =BO 0 ij f 1 ( 0 i ; 0 j )f 1 ( 0 i ; 0 j )f 4 ( 0 i ;BO 0 ij )f 5 ( 0 j ;BO 0 ij ); (2.37) BO ij =BO 0 ij f 1 ( 0 i ; 0 j )f 1 ( 0 i ; 0 j )f 4 ( 0 i ;BO 0 ij )f 5 ( 0 j ;BO 0 ij ); (2.38) 0 i = X j2n(i) BO 0 ij Val i ; (2.39) 27 Table 2.2: The parameters of the bond order functions. p bo1 p bo6 bond parameters r 0 ,r 0 ,r 0 bond radius parameters Val i ,Val boc i valencies of atom i p boc1 ,p boc2 overcoordination parameters p boc3 p boc5 13 bond order corrections f 1 ( 0 i ; 0 j ) = 1 2 Val i +f 2 ( 0 i ; 0 j ) Val i +f 2 ( 0 i ; 0 j ) +f 3 ( 0 i ; 0 j ) + Val j +f 2 ( 0 i ; 0 j ) Val j +f 2 ( 0 i ; 0 j ) +f 3 ( 0 i ; 0 j ) ; (2.40) f 2 ( 0 i ; 0 j ) = exp(p boc1 0 i ) + exp(p boc1 0 j ); (2.41) f 3 ( 0 i ; 0 j ) = 1 p boc2 ln 1 2 exp(p boc2 0 i ) + exp(p boc2 0 j ) ; (2.42) f 4 ( 0 i ;BO 0 ij ) = 1 1 + exp(p boc3 (p boc4 BO 0 ij BO 0 ij 0boc i ) +p boc5 ) ; (2.43) f 4 ( 0 j ;BO 0 ij ) = 1 1 + exp(p boc3 (p boc4 BO 0 ij BO 0 ij 0boc j ) +p boc5 ) ; (2.44) 0boc i =Val boc i + X j2n(i) BO 0 ij ; (2.45) wheren(i) is a group of atoms close to atomi, which are placed within the atom's single covalent bond cuto distance. The list of all the parameters in this section is provided in Table 2.2. The corrected overcoordination of atomi, after correcting bond orders, is computed as i =Val i + X j2n(i) BO ij ; (2.46) 28 2.4.2 Charge calculation Atomic charges are computed in ReaxFF using an electronegativity equalization method (EEM) [69, 70]. Sanderson's electronegativity equalization principle [71{ 73] considers all included atoms in a molecule to have equal electronegativity. In addition, one can use DFT [69,74,75] to calculate the electronegativity of atom i in a molecule as i =A i +B i q i + N X j=1(j6=i) q i R ij ; (2.47) A i = 0 i + i ; (2.48) B i = 2( 0 i + i ); (2.49) where N is the number of atoms in the molecule, q i is charge, i is electronega- tivity, i is hardness of atom i, R ij is the distance between atoms i and j. , i and i are parameters to be determined empirically and describe the molecular environment. Sanderson's electronegativity equalization principle sets the elec- tronegativity of all atoms in a molecule equal to molecule's electronegativity, 1 = 2 = = N = ; (2.50) where is molecular electronegativity. Conservation of charge gives N X i=1 q i =Q; (2.51) 29 where Q is the total charge of the molecule. Finally, atomic charges can be calcu- lated in terms offA i g,fB i g, and Q as 2 6 6 6 6 6 6 6 6 6 6 4 B 1 R 12 R 1N 1 R 21 B 2 R 2N 1 . . . . . . . . . . . . . . . R N1 R N2 B N 1 1 1 1 0 3 7 7 7 7 7 7 7 7 7 7 5 2 6 6 6 6 6 6 6 6 6 6 4 q 1 q 1 . . . q 1 3 7 7 7 7 7 7 7 7 7 7 5 = 2 6 6 6 6 6 6 6 6 6 6 4 A 1 A 1 . . . A N Q 3 7 7 7 7 7 7 7 7 7 7 5 : (2.52) In this dissertation, the electrostatic energy is minimized with a charge neu- trality constraint, equivalent to the EEM method. Consider the energy to be an arbitrary function of the atomic charges, E = E(fq i g) with the charge neu- trality constraint (i.e. P i q i = 0). The Lagrange multiplier method is utilized to minimize the energy subject to the charge neutrality constraint. Therefore, E = E(fq i g) P i q i can be minimized without a constraint, where is the Lagrange multiplier. One can see that all atoms have same electronegativity when the system is minimized, g i =: @E @q i = @E @q i = 0; (2.53) where @E=@q i is the chemical potential of atom i, or negative of electronegativity of atom i. The conjugate-gradient method is used to minimize energy. 30 Table 2.3: Classication of the interactions in ReaxFF. Number of explicitly involved atoms Energy terms 1 E lp +E over +E under 2 E bond +E vdWaals +E Coulomb 3 E val +E pen +E coa +E hbond 4 E tors +E conj 2.4.3 Potential energy There are twelve energy terms in ReaxFF that make the overall potential energy, E ReaxFF =E bond +E lp +E over +E under +E val +E pen +E coa +E tors +E conj +E hbond +E vdWaals +E Coulomb ; (2.54) whereE bond is bond energy,E lp is lone pair energy,E over is overcoordination energy, E under is undercoordination energy, E val is valence-angle energy, E pen is penalty energy, E coa is three-body conjugation energy, E tors is torsion-angle energy, E conj is four-body conjugation energy, E hbond is hydrogen bonding energy, E vdWaals is van der Waals energy, and E Coulomb is Coulomb energy. Only the last two terms are the energy of non-covalent interactions, while the other terms are the energy of valence interactions. These interactions can be classied as one-, two-, three-, and four-body inter- actions, based on the number of explicitly-included atoms; see Table 2.3 and Fig- ure 2.1. The presence of bond orders in ReaxFF ensures smooth vanishing of valance interaction energies as the bonds break. 31 Figure 2.1: Schematic illustration of the atomic congurations in energy terms: (a) one-body, (b) two-body, (c) three-body, (d) four-body, (e) hydrogen-bonding, and (f) non-covalent iterations. A red sphere indicates the position of the primary atom in each energy term. Solid lines indicate covalent bonds, while non-covalent bonds are represented by dotted lines. Picture courtesy of [1]. 2.4.4 One-body interactions As it was mentioned above, the lone, overcoordination, and undercoordination energies are classied for single atoms. The number of lone pairs around atom i is computed as n lp;i = e i 2 + exp $ p lp1 2 + e i 2 e i 2 2 % ; (2.55) where e i calculates the dierence between the sum of bond orders around atom i and the total number of outer shell electrons (e.g. 1 for hydrogen, 6 for oxygen or 4 for silicon), e i =Val e i + X j2n(i) BO ij : (2.56) 32 If for an atom i the assigned bond order became more than its normal value, then Equation 2.55 would let a lone pair to break down smoothly. This causes a deviation lp i from the optimum number of loan pairs n lpopt;i , computed as lp i =n lpopt;i n lp;i : (2.57) Having computed lp i , the penalty energy function for atoms with lone electron pairs exceeding the optimal number is dened by E lp = X i p lp2 lp i 1 + exp 75 lp i : (2.58) Even after correcting the bond orders as explained in Section 2.4.1, those of some atoms may still be more than what the theory of valence electrons antici- pates. In that case energy penalties are assigned to overcoordinations. However, as mentioned earlier, the lone pair electrons break up through Equation 2.55 and the resulting overcoordinations should not be punished. Therefore, rst the corrected overcoordination lpcorr i is computed as lpcorr i = i lp i 1 +p ovun3 exp h p ovun4 P j2n(i) j lp j BO ij +BO ij i ; (2.59) after obtaining corrected overcoordinations, each over coordinated atom i is assigned an energy penalty given by E over = X i P j2n(i) p ovun1 D e BO ij lpcorr i +Val i lpcorr i 2 4 1 1 + exp p ovun2 lpcorr i 3 5 : (2.60) 33 Table 2.4: The parameters of one-body energy functions. Val e i valency of atom i n lpopt;i normal number of lone pair electrons for atom i parameters p lp1 lone pair parameter p lp2 lone pair energy parameter p ovun1 p ovun4 overcoordination parameters p ovun5 undercoordination energy parameter p ovun6 p ovun8 undercoordination parameters D e -bond energy parameter If two undercoordinated atoms i and j are bonded together and the bond between them has partially -bond character, the energy contribution resulting from the resonance of-electrons between centers of two connected atoms is taken into account as E under = X i (p ovun5 ) 1 exp p ovun6 lpcorr i 1 + exp p ovun2 lpcorr i 1 1 +p ovun7 exp h p ovun8 P j2n(i) j lp j BO ij +BO ij i : (2.61) The list of all the parameters for the one-body energy functions is given in Table 2.4. 2.4.5 Two-body interactions For all pairs of atoms, the energies for the bond, van der Waals, and Coulomb are calculated, when the distance between two atoms is less than the cuto length of 34 the related interaction. The energy contributions from -, -, and -bonds are added together in calculating the bond energy as E bond = X i X j D e BO ij exp p be1 1 BO ij p be2 D e BO ij D e BO ij : (2.62) As mentioned earlier, ReaxFF takes into account two types of reactions. One is valance interactions, while the other one is non-covalent (non-bonded) interactions. The rst one results from the overlapping between the electron wave functions, whereas the second one is based on van der Waals and Coulomb interactions, where Pauli exclusion principle for repulsion at short distances and attraction due to dispersion at long distances are implemented. In ReaxFF, to screen the van der Waals and Coulomb forces, a taper function, which is a distance-dependent 7 th order polynomial, is used, Tap (r ij ) = 7 X =0 Tap r ij r cnb : ; (2.63) where r ij is the distance between atoms i and j, and r cnb is the cuto length, which is used when non-covalent energies are computed for all pairs of atoms with distances less than r cnb . Tap are coecients and their values are given in Table 2.5. The above taper function is multiplied by the energy and the force when non- covalent energy and force are computed between a pair of atoms. The Tap are selected to result in continuous derivatives (1 st , 2 nd , and 3 rd order) with respect to r ij , and allow them to vanish smoothly at the cuto distance. 35 Table 2.5: Values of Tap in the taper function equation. Tap Values Tap 0 1 Tap 1 0 Tap 2 0 Tap 3 0 Tap 4 -35 Tap 5 84 Tap 6 -70 Tap 7 20 In ReaxFF, a shielded Morse potential is used to account for the van der Waals interaction, as given by E vdWaals = X i X j D ij Tap (r ij ) exp ij 1 f 13 (r ij ) r vdW 2 exp 1 2 ij 1 f 13 (r ij ) r vdW : (2.64) f 13 (r ij ) = r p vdW1 ij + 1 w p vdW1 1=p vdW1 : (2.65) Large repulsions between connected atoms and also between atoms sharing a valence angle are avoided by using a shielded potential. The Coulomb interac- tions are also calculated by a shielded potential as E Coulomb = X i i q i + 1 2 X i X j q i H (r ij )q i ; (2.66) H (r ij ) =J i ij + Tap (r ij ) r 3 ij + (1= ij ) 3 1=3 : (2.67) The list of all the parameters for the two-body energy functions is given in Table 2.6. 36 Table 2.6: The parameters of two-body energy functions. D e ;D e ;D e -, - and -bond energy parameters p be1 ;p be2 bonding energy parameters D ij van der Waals energy parameter ij van der Waals parameter r vdW van der Waals parameter p vdW1 van der Waals shielding w van der Waals parameter ij Coulomb parameter 2.4.6 Three-body interactions The valence angle, hydrogen bonding, three-body conjugation, and penalty are taken into account for triplets of the atoms. The valence angle considers an angle formed across two bonds with three atoms in which one of them is common between the two bonds. For instance, the valence angle in a water molecule is about 104:5 . Over/under coordination, hybridization, and the number of lone pair electrons of central atom in the triplet are the important factors that aect the value of equilibrium valence angle 0 . The hybridization of the central atom itself is a function of the sum of -bond orders (SBO) of the central atom. When lone pair electrons and over/under coordination are not present, 0 is equal to 109:5 , 120 , and 180 for sp 3 (0 -bond), sp 2 (1 -bond), and sp (2 -bonds) hybridization, respectively. In ReaxFF, 0 is calculated as 0 (SBO) = 0;0 f1 exp [p val10 (2SBO2)]g ; (2.68) SBO2 = 8 > > > > > > > < > > > > > > > : 0 (SBO 0) SBO p val9 (0<SBO 1) 2 (2SBO) p val9 (0<SBO 2) 2 (SBO> 2) ; (2.69) 37 SBO j = X m2n(j) BO jm +BO jm + 2 4 1 Y m2n(j) exp BO 8 jm 3 5 angle j p val8 n lp;j ; (2.70) angle j =Val angle j + X m2n(i) BO jm : (2.71) The list of all the parameters in the equations of this section is given in Table 2.7. Energy contribution resulting from the deviation of valence angle ijk from its equilibrium value 0 is calculated by E val = X i X j X k f 7 (BO ij )f 7 (BO jk )f 8 ( j ) p val1 p val1 exp p val2 ( 0 (SBO j ) ijk ) 2 ; (2.72) f 7 (BO ij ) = 1 exp p val3 BO p val4 ij ; (2.73) f 8 ( j ) =p val5 (p val5 1) 2 + exp p val6 angle j 1 + exp p val6 angle j + exp p val7 angle j : (2.74) When one of the bonds forming the valence angles breaks up, the energy contri- bution from a valence angle approaches zero in the function E val . When two double bonds are included in a valence angle, as in, for example, propadiene, an additional energy term is used to penalize deviations from the equilibrium angle E pen = X i X j X k p pen1 f 9 ( j ) exp p pen2 (BO ij 2) 2 exp p pen2 (BO jk 2) 2 ; (2.75) f 9 ( j ) = 2 + exp (p pen3 j ) 1 + exp (p pen3 j ) + exp (p pen4 j ) : (2.76) 38 In chemistry, conjugation means the delocalization of-bonds in a group. This would result in lowering the overall energy and increasing the stability of the group. In ReaxFF, conjugation is used for two terms, the four- and three-body conjugation terms. The four-body conjugation takes into account the stability of conjugated hydrocarbon systems, for example benzene. On the other hand, three- body conjugation term considers the stability of the groups, for example the nitro group ({NO 2 ). The three-body conjugation energy E coa is given by E coa = X i X j X k p coa1 1 1 + exp p coa2 val j exp 2 4 p coa3 0 @ BO ij + X m2n(i) BO im 1 A 2 3 5 exp 2 4 p coa3 0 @ BO jk + X m2n(k) BO km 1 A 2 3 5 exp p coa4 (BO ij 1:5) 2 exp p coa4 (BO jk 1:5) 2 : (2.77) A hydrogen bond is usually dened as electromagnetic attractive interaction between a hydrogen atom with other atoms, such as a nitrogen, oxygen or uorine atom. However, in ReaxFF, for any triplet in which hydrogen is the central atom, a hydrogen bonding energy is dened, and is given by E hbond = X i X j X k p hb1 [1 exp (p hb2 BO ij )] exp p hb3 r 0 hb r jk + r jk r 0 hb 2 sin 8 ijk 2 : (2.78) 39 Table 2.7: The parameters of three-body energy functions. p val1 valence angle energy parameter p val2 p val7 valence angle parameters p val8 valency/lone pair parameter p val9 ;p val10 valence angle parameters 0;0 equilibrium valence angle parameter Val angle i atom valency p pen1 penalty energy parameter p pen2 double bond/angle parameter p pen3 ;p pen4 double bond/angle parameters for overcoordination p coa1 three-body conjugation energy parameter p coa2 p coa4 valency angle conjugation parameters p hb1 hydrogen bond energy parameter p hb2 ;p hb3 hydrogen bond parameters r 0 hb hydrogen bond radius 2.4.7 Four-body interactions Four-body conjugation and torsion angle energies are dened for a connected group of four atoms. For each quartet, with atoms being connected consecutively, the torsion or dihedral angles is dened as the angle between the planes formed by the rst three and the last three atoms. In other words, a torsion or dihedral angle ! ijkl is dened between two triplet of atoms (i;j;k) and (j;k;l) as ! ijkl = cos 1 (^ n ijk ^ n jkl ) ; (2.79) where ^ n ijk and ^ n jkl are unit vectors perpendicular to the planes dened by triplets of atoms (i;j;k) and (j;k;l), respectively. In ReaxFF, to ensure that dihedral 40 angles are always near their equilibrium values, angle energy contribution E tors is dened as E tors = 1 2 X i X j X k X l f 10 (BO ij ;BO jk ;BO kl ) sin ijk sin jkl V 1 (1 + cos! ijkl ) +V 2 exp n p tor1 2BO jk f 11 ( j ; k ) 2 o (1 cos 2! ijkl ) +V 3 (1 + cos 3! ijkl ) : (2.80) f 10 (BO ij ;BO jk ;BO kl ) = n [1 exp (p tor2 BO ij )] [1 exp (p tor2 BO jk )] [1 exp (p tor2 BO kl )] o ; (2.81) f 11 ( j ; k ) = ( 2 + exp h p tor3 angle j + angle k i ) ( 1 + exp h p tor3 angle j + angle k i + exp h p tor4 angle j + angle k i ) 1 : (2.82) If one of the valence angles ( ijk or jkl ) approaches , the contribution of the torsion angles disappears. The list of all the parameters in the equations of this section is given in Table 2.8. As mentioned above, the four-body conjugation energy takes into account the stability of hydrocarbon conjugated-systems. The energy E conj is calculated by E conj = X i X j X k X l p cot1 f 12 (BO ij ;BO jk ;BO kl ) 1 + cos 2 ! ijkl 1 sin ijk sin jkl ; (2.83) 41 Table 2.8: The parameters of four-body energy functions. p val1 torsion energy parameters p val2 p val7 torsion parameter p val8 torsion/BO parameter p val9 ;p val10 torsion overcoordination parameters 0;0 four-body conjugation energy parameter Val angle i four-body conjugation parameter f 12 (BO ij ;BO jk ;BO kl ) = exp p cot2 (BO ij 1:5) 2 exp p cot2 (BO jk 1:5) 2 exp p cot2 (BO kl 1:5) 2 : (2.84) The largest contribution of the four-body conjugation is when consecutive bonded atoms have bond orders equal to 1.5 (e.g benzene or other aromatics). 2.5 Dual Control Volume-Grand Canonical Molecular Dynamics The dual control volume-grand canonical molecular dynamics (DCV-GCMD) is considered as one of the nonequilibrium MD simulation methods that combines MD and MC simulations [76,77]. This method is ideal for nonequilibrium systems under an external driving force, such as pressure gradient or chemical potential [4,78,79]. Using this method, it is possible to study the transport and separation of gas mixtures in membranes, such as H 2 /CO 2 and H 2 /CH 4 in the present dissertation (see Chapter 6). In this section, the essentials of this method with focus on the simulation techniques that we used in our home-made code are discussed. 42 Figure 2.2: Schematic illustration of the SiC amorphous membrane used in the DCV-GCMD simulations. The h and ` areas are, respectively, the high- and low- pressure control volumes. 2.5.1 The DCV-GCMD method In the DCV-GCMD method, the system under consideration is divided into three regions, similar to the schematic representation shown in Figure 2.2. Here, the origin of the coordinates is located at the center of the system. The external driving force is a pressure gradient (see Chapter 6) applied in the direction of the uid ow (the -x direction). The control volumes (CV) chambers on the left (h) and right side (`) represent, respectively, the high- and low-pressure regions, in equilibrium with the bulk condition. The middle region is allocated for the atomistic structure of the model under study { the amorphous SiC membrane in this dissertation [4,78]. The membrane's length is nL with n being an integer. In our calculations (Chapter 6) we used a variety of membranes lengths from n = 1 to 10. Periodic boundary conditions were employed in the y- and z-directions and, therefore, it is only x-direction that is used to study the eect of the membranes' length on transport properties such separation factor and permeance [4,78]. 43 Table 2.9: The Lennard-Jones parameters for the various gases and membrane atoms, used in the DCV-GCMD simulations. Gases and atoms ( A) =k B (K) Mass C 3.400 28.0 12.00 Si 3.742 23.6 28.09 H 2 2.827 59.7 2.04 CO 2 3.794 225.3 44.01 CH 4 3.810 148.1 16.04 Table 2.10: The conversion between the actual and reduced units. Subscript 1 refers to the value of the parameters for component 1, which is CH 4 in this disser- tation. Table courtesy of [4]. Variable Reduced form Length L L =L= 1 Energy U U =U=" 1 Mass M M =M=M 1 Density = 3 1 Temperature T T =k B T=" 1 Pressure P P =P 3 1 =" 1 Time t t =t(" 1 =M 1 2 1 ) 1=2 Flux J J =J 3 1 (M 1 =" 1 ) 1=2 Permeability K K =K(M 1 " 1 ) 1=2 = 1 The uids and atoms of the membrane are represented as Lennard-Jones (LJ) spheres, characterized by eective LJ size and energy parameters, and ". The values of and" for the gaseous molecules (H 2 , CH 4 , and CO 2 ) and the membrane atoms (Si and C), which are used in this dissertation, are given in Table 2.9 [78,80, 81]. For unlike pairs of atoms, the Lorentz-Berthelot [82, 83] rules, ij = ( i j ) 1=2 , and ij = 1 2 ( i + j ), were used to calculate the size and energy of the pairs. All the quantities of interest were made dimensionless with the help of and " of one of the uids (CH 4 in this dissertation); the dimensionless groups are listed in Table 2.10. 44 As mentioned above, in the DCV-GCMD method the MD moves are combined with the GCMC insertions and deletions in the h and ` control volumes (see Fig- ure 2.2). The Verlet velocity algorithm was implemented in the MD simulations to solve the equations of motion [84]. The main purpose of using the DCV-GCMD technique in this dissertation is to study the eect of an external pressure drop on the transport properties of gaseous mixtures and, therefore, the contribution of temperature gradient to the transport must be eliminated by holding the temperature constant in the CVs, as well as within the membrane. This was achieved by rescaling the velocity components independently in all the three directions during the entire MD simulation [4,78,79]. The results are shown in Chapter 6. The density of each gaseous component in the two CVs must be held at some xed values. These values, which are in equilibrium with two bulk phases, are directly related to the required chemical potentials and pressures in the correspond- ing CVs. To obtain constant values of the densities, sucient number of the GCMC insertions and deletions of the molecules were carried out in the CVs [4,78,79]. The probabilities of inserting and deleting a particle of componenti are, respec- tively, given by p + i = min Z i V c N i + 1 exp(E=k B T ); 1 ; (2.85) and p i = min N i Z i V c exp(E=k B T ); 1 ; (2.86) where E is the potential energy change resulting from creating or deleting a particle,V c the volume of the CV,N i the number of atoms of componenti in each 45 CV, k B the Boltzmann's constant and Z i the absolute activity at temperature T . The values of Z i are calculated by Z i = exp( i =k B T ) 3 i ; (2.87) where i and i are, respectively, the de Broglie wavelength and chemical potential of component i. The pressures were converted to equivalent chemical potentials using Soave-Redlich-Kwong (SRK) equation of state [85]. The Maxwell-Boltzmann distribution was used to assign a thermal velocity to every newly inserted particle in a CV at the given T . In the DCV-GCMD simulations, it is crucial to select an appropriate value for the ratioR, which is the number of GCMC insertions and deletions in each CV relative to the number of the MD steps between successive GCMC steps. This ratio must be chosen such that the correct density and pressure in the CVs are created and maintained during simulations [4,78]. In our simulations (Chapter 6) the best value forR was found to be 10:1. During the MD simulations the molecules crossing the outer boundaries of the CVs inx-direction were removed, but the periodic boundary condition was used in the other two directions, and the image of particles were found in the corresponding CVs. The number of such particles was, however, very small compare to the total number of molecules that were removed during the GCMC simulations. In addition, for each component and in the membrane region a nonzero streaming velocity (the ratio of the ux to the concentration of each component) was allowed, consistent with the presence of bulk pressure/chemical potential gradients along the ow direction. The absence of nonzero streaming velocity results in severely underestimated uxes in the transport region, used in many of the previous works. 46 The actual streaming velocities of the molecules in the membrane region were calculated by the MD simulations [4,78]. 2.5.2 Potential models of the gases and the membrane The standard cut-and-shifted LJ potential was used to model the molecule{ molecule interactions, given by U(r) = 8 < : U LJ (r)U LJ (r c ); rr c ; 0; r>r c : (2.88) Here, U LJ (r) is the standard 6-12 LJ potential, U LJ (r) = 4" r 12 r 6 ; (2.89) where r is the distance between the interacting pair, and r c is the cut-o distance which was chosen to be 2:5 CH 4 . As discussed in Chapter 6, the DCV-GCMD simulations for the present study are computationally expensive. Therefore, the simulation time for computing the interactions between gas molecules and the Si and C atoms in the membrane was reduced by discretizing the simulation box into n x n y n z grid points, with the distance between two consecutive grid points being equal to 0.43 A. For example, for a length of 64 A, 149 grid points were used in corresponding direction. The interactions were calculated at each of the grid points, and 3D piecewise cubic Hermite interpolation [86] was used for interpolating the results between the grid points. The results were then saved and used in the simulations. 47 2.5.3 Transport properties Several quantities of interest, such as the density prole, ux, and the separation factor, were computed for the gaseous mixtures. The density prole z i (x) of the gas was computed by dividing the simulation box in the x-direction into grids of size CH 4 , and averaging the number of particles of component i over the distance CH 4 for each MD step [4,78]. As discussed in Chapter 6, the density proles are important to understanding the transport of the gases within the two CVs and the membrane region. The uxJ i of each componenti was calculated by counting the net number of its particles crossing a given yz plane, J i = N LR i N RL i NtA yz ; (2.90) whereN LR i andN RL i are the number of the particles of type i displacing from the left to the right and vice versa, respectively, A yz is the cross-sectional area, t is the MD time step (t = 1 10 3 , i.e., t' 0:00137 ps), and N is the number of the MD steps over which the average was taken (e.g every 1 million steps). The ux J i was calculated at three locations along the membrane (i.e. entrance, middle, and the exit of the membrane), and was then averaged. The system was considered to have reached steady state when the uxes computed at various yz planes were within 5% from the averaged value [4,78]. The equations of motion were integrated over large number of steps (up to 100 million steps) within which the steady state was achieved. The permeabilityK i of species i was computed using K i = J i P i =nL = nLJ i P i ; (2.91) 48 where P i =x i P is the partial pressure drop for species i along the membrane, withx i being the mole fraction of component i (x i = 0:5 in this dissertation), and P the total pressure drop applied externally to the membrane [4, 78]. A most important property that is computed for gaseous mixtures and compared with experimental data is the dynamic separation factor S 21 , dened as S 21 = K 2 K 1 : (2.92) 49 Chapter 3 Development of a Reactive Force Field 3.1 Introduction As discussed in Chapter 2, the original ReaxFF was developed [26, 32] to make practical molecular dynamics (MD) simulation of large-scale reactive compounds with thousands of atoms (see Section 1.4). The conventional empirical force elds used in the standard MD simulations assume rigid connectivity of the atoms and cannot describe the complex chemistry of real materials, but ReaxFF may be used to examine such chemical processes as polymer decomposition. It uses a general relationship between bond distance and bond order, on the one hand, and between bond order and bond energy, on the other hand, which lead to proper dissociation of the bonds and atoms. Its development is based on the fundamental assumption that the bond order between a pair of atoms can be obtained directly from the interatomic distance r ij between atoms i and j. Other valence terms present in ReaxFF are also dened in terms of the same bond order, so that all the terms vanish smoothly as the bonds break. In addition, ReaxFF contains Coulomb and van der Waals potentials for describing nonbond interactions between all the atoms. They are shielded at short range, so that they become constant as the distance between two atoms vanishes. 50 In this chapter we describe the development of ReaxFF [87] in order to simulate the thermal decomposition of the polymeric precursor that is used in the fabrication of SiC nanoporous membranes. In Sections 3.2 and 3.3, we explain the quantum chemical (QC) computations for estimating the parameters of ReaxFF for the pyrolysis process of the polymeric precursor, hydridopolycarbosilane (HPCS). In Section 3.4, we describe the generation of a molecular model of the polymer and the details of minimization of its energy and its pre-equilibration, before it is heated up. Section 3.5 describes the results of MD simulation of thermal decomposition of the polymer [87]. 3.2 Quantum Chemical Calculations Similar to nonreactive empirical force elds, ReaxFF computes the total energy E of a system as the sum of partial energies, each representing a particular contribu- tion to E, E =E bond +E lp +E over +E under +E val +E pen +E coa +E tors +E conj +E hbond +E vdWaals +E Coulomb : (3.1) Each partial energy term is parameterized by a number of parameters. The description of each term and their parameterized functional forms are given in Section 2.4. The overall number of parameters for the HPCS is 438, of which 39 are general, 50 are computed specically in this dissertation for the HPCS, while the remaining 349 parameters are those computed previously [38] for ReaxFF for poly(dimethylsiloxane). Since we have three types of atoms in the HPCS, namely, Si, C, and H, then, there are 26 combinations of torsion angles. The torsion energy, 51 however, contributes weakly to ReaxFF and, thus, only 9 terms out of 26 are used in the ReaxFF that is developed here. The parameters of ReaxFF were estimated such that they correctly reproduce the rst principles quantum mechanical interactions in the HPCS, and provide a way of extending the rst-principles' accuracy to large length- and time-scale sim- ulations. This was achieved by adding the SiC crystal chemistry and the QC data for systems relevant to the HPCS to the ReaxFF training set for Si materials [32]. The QC data for nonperiodic systems were obtained from the density-functional theory (DFT) calculations, carried out using the Jaguar (version 5.5) software [88], B3LYP functional [59, 89], and Pople's 6-311G** basis set [90, 91]. The Mulliken charges were also computed using the 6-31G** basis set developed by Pople and co-workers [48,49]. To obtain the equation of state for the -SiC crystal structure, the DFT cal- culations were carried out with the software Quantum-Espresso [92]. The general- ized gradient approximation (GGA), proposed by Perdew et al. [60] (see Chapter 2) was used for the exchange-correlation energy, and ultrasoft pseudopotentials were utilized to replace the core electrons. We used the Perdew's implementa- tion of the GGA [93] with a kinetic energy cuto of 320 eV. The Monkhorst-Pack scheme [94] was used to generate the k-space grid with a 0.1 A 1 spacing. In gen- eral, the accuracy of the results are similar to, or better than that provided by the parameterized model number 3 (PM3) [95], a semi-empirical method for quantum mechanical computation of the electronic structure of molecules. ReaxFF is about 100 times faster than PM3, which in turn is about 100 times faster than the QC calculations. The DFT calculations were then utilized to compute the geometries and ener- gies of the species and complexes in the HPCS, as well as the products of the 52 pyrolysis. The ReaxFF parameters were determined by combining the Si training set with the results from the DFT calculations, followed by optimization of the parameters, in order to minimize the dierences between the energies calculated by the DFT and ReaxFF. 3.3 Optimization of ReaxFF To optimize the Si{C bond energy, we used the dissociation pathway for the bond in H 3 Si{CH 3 , shown in Figure 3.1(a), and the Si=C bond in H 2 Si=CH 2 , depicted in Figure 3.1(b). To optimize the ReaxFF parameters, the DFT results for the singlet case were used around the equilibrium bond distance. The valence angle parameters were optimized to represent the angle-bending energies (between three atoms), obtained from the DFT calculations on various clusters of atoms. There are six valence angles in the HPCS that need to be accounted for, namely, Si{C{Si, Si{C{H, C{Si{H, C{C{Si, C{Si{C, and C{Si{Si. For the rst three the triplets H 3 Si{CH 2 {SiH 3 and for the rest H 3 C{CH 2 {SiH 3 , H 3 C{SiH 2 {CH 3 , and H 3 Si{SiH 2 { CH 3 were used, respectively, as the representative molecules for valence angles. For each angle, the geometry of the representative molecule was optimized with the angles of interest xed, in order to compute the angle-bending energies relative to the optimal geometry. The results obtained by the DFT computations are compared in Figure 3.2 with those that ReaxFF yielded. The charge distributions in ReaxFF were calculated using the electronegativity equalization method (EEM) [95] (see Chapter 2). The Mulliken charge distribu- tions obtained from the DFT calculations were used to optimize the EEM charge parameters. The partial charges for SiH 3 CH 2 SiH 3 , CH 3 SiH 2 CH 3 , and C{(SiH 2 { CH 2 ) 3 (six-membered ring), computed by the DFT and ReaxFF, are compared 53 Figure 3.1: Comparison of the energies computed by ReaxFF and QC for (a) Si{C and (b) double-bond dissociations. 54 Figure 3.2: Comparison of distortion energies computed by ReaxFF and QC for the angles in (a) Si{C{Si; (b) Si{C{H; (c) C{Si{H; (d) C{C{Si; (e) C{Si{C, and (f) C{Si{Si triplets in the HPCS. in Figure 3.3. They all indicate good agreement between the two sets of results, indicating the accuracy of ReaxFF. 55 Figure 3.3: Comparison of the Mulliken charge distributions computed by the DFT (red) and ReaxFF (blue). To sample the possible reactions in the Si{C materials relevant to the decompo- sition of the HPCS, the DFT- and ReaxFF-computed reaction energies were com- puted and compared. The reactions include not only the energetically-favorable clusters, but also the strained and energetically less favorable clusters that test the ability of ReaxFF to represent both cases. Cleavage of the Si{H and C{H bonds are important in the degradation process of the HPCS, when exposed to heat. The energies computed by ReaxFF and DFT for these systems are listed in Table 3.1. As pointed out in Chapter 1, the goal in this dissertation is to develop ReaxFF in order to carry out MD simulation of the structure that results from pyrolysis of the HPCS. The structure is a condensed phase of SiC that may be either in crystalline or amorphous form, particularly at high pressures and temperatures. 56 Table 3.1: Comparison of the reaction energies E (in kcal/mol), computed by the DFT and ReaxFF Reaction E (DFT) E (ReaxFF) SiH 2 (CH 3 ) 2 + CH 2 (SiH 3 ) 2 ! c-(SiH 2 {CH 2 ) 3 + 2H 2 17.15 12.75 SiH 2 (CH 3 ) 2 + CH 2 (SiH 3 ) 2 ! (CH 3 ) 2 SiH{CH(SiH 3 ) 2 + H 2 7.86 5.15 CH 3 (SiH 2 CH 2 ) 2 SiH 3 ! c-(SiH 2 {CH 2 {SiH 2 {CH 2 {CH){SiH 3 + H 2 27.84 22.08 CH 3 (SiH 2 CH 2 ) 2 SiH 3 ! c-(SiH 2 {CH 2 {SiH 2 {CH){SiH 2 CH 3 + H 2 24.43 25.81 CH 3 (SiH 2 CH 2 ) 2 SiH 3 ! c-(SiH 2 {CH 2 {SiH 2 {CH 2 {SiH){CH 3 + H 2 13.00 17.88 There are, in fact, more than 250 crystalline structures of SiC, and three major SiC polytypes known as () 3C-SiC, () 6H-SiC, and 4H-SiC [96]. Polytypes refer to a family of similar crystalline structures, and represent variations of the same chemical compound that are identical in two dimensions, but dier in the third one. Thus, the ability of ReaxFF for representing the equation of state of the SiC crystal was tested against the aforementioned crystalline structures of SiC. The Crystal Builder module [97] in Cerius 2 and the unit cell parameters reported by others [98,99] were used to construct the SiC crystals. For each crystalline form of SiC, the quantum energies were computed for a broad range of volumes (densities), describing both expansion and compression. Figure 3.4 compares the ReaxFF results against those obtained by the DFT calculations, indicating that ReaxFF correctly describes the equation of state for the SiC structures. There is, however, a rather large dierence between the energies computed by ReaxFF and DFT for 4H-SiC at the expansion part of -SiC. This does not, however, cause a problem for the particular conditions that we study in this dissertation, because the most relevant and physically interesting parts are the equilibrium volumes from which the relative phase stabilities are deduced. 57 Figure 3.4: The equations of state (compression and expansion) for the three crystal structures of SiC. (a) -SiC; (b) 4H-SiC, and (c) -SiC, as computed by the DFT and ReaxFF. Optimization of -SiC crystal with ReaxFF yielded a density of 3.16 g/cm 3 , only 1.6% smaller than the experimental density of 3.21 g/cm 3 , measured at 300 K [99]. Values of the ReaxFF parameters that we computed for the atoms, bond ener- gies and bond orders, van der Waals, angle, and torsion contributions to the HPCS and SiC crystals are given, respectively, in Tables 3.2 to 3.6. They may be used directly in any MD simulator, and in particular those that are freely available, such as the LAMMPS package [100]. A major goal of the development of ReaxFF, in addition to its use in the MD simulation of thermal decomposition of HPCS, is to estimate the parameters that are transferable to other Si-based materials and polymers. Since the complete training set has been developed for the Si-containing polymers, the new parameters of ReaxFF that we have computed, in addition to 58 Table 3.2: The atom parameters of ReaxFF, computed for the HPCS. The rst three from the right are the EEM parameters, while the rest are the bond-order correction coecients. Atom (eV) (eV) ( A) p ov=uv p val3 p boc3 p boc4 p boc5 C 5.725 6.923 0.871 -2.898 4.782 33.563 5.713 11.996 H 3.819 9.883 0.762 -15.768 3.350 3.246 3.865 1.000 Si 2.912 7.427 0.629 -5.849 2.285 8.0583 4.511 0.838 Table 3.3: Bond energiesD e (in kcal/mol) and bond-order parameters of ReaxFF, computed for the HPCS. The parameter r is in A. Pair D e; D e; p be1 p bo1 p bo2 p bo3 p bo4 r r C{H 159.8520 0.000 -0.4646 -0.0098 8.5954 1.0386 Si{H 238.3492 0.000 -0.7467 -0.0438 6.8851 1.3266 C{Si 100.3210 25.714 0.1000 -0.1100 6.0000 -0.290 8.80 1.640 1.570 be relevant to the HPCS, make the FF general and capable of describing other polymers that contain other combinations of Si with C and H. We point out that the important point in using ReaxFF is to enable chemical reactions to proceed with appropriate energy barriers. Thus, ReaxFF uses a bond energy-to-bond order-to-bond distance relation that can capture double and single bonds, as well as lower-order bonds during reactions. As a result, the compress- ibility curves computed by ReaxFF are generally not as accurate as one would like them to be. For covalent systems calculations of accurate second-order properties are best carried out by using ReaxFF to obtain the structure, and then using a nonreactive FF to analyze the distortions near an equilibrium structure. Table 3.4: van der Waal parameters of ReaxFF, computed for the HPCS. Pair ( A) (kcal/mol) vdW C{H 1.7207 0.0431 10.3632 3.2757 Si{H 1.4723 0.0640 13.8391 2.6229 C{Si 1.9000 0.0600 11.0000 1.7008 59 Table 3.5: Parameters for the triplets in ReaxFF, computed for the HPCS. The angles are in degrees. Triplet 0;0 p val1 p val2 p val4 p val7 p coa1 p pen1 C{C{Si 72.52 22.3583 2.0393 1.040 1.0031 0.000 0.000 C{Si{C 69.33 18.9270 2.1333 1.040 1.0031 0.000 0.000 Si{Si{C 69.33 19.6964 2.0703 1.040 1.0031 0.000 0.000 Si{C{Si 69.33 18.9270 2.0703 1.040 1.0031 0.000 0.000 Si{C{H 72.59 13.8347 2.4952 1.040 1.0000 0.000 0.000 C{Si{H 72.59 14.8347 2.4952 1.040 1.0000 0.000 0.000 Table 3.6: Torsion parameters of ReaxFF for quadruplets, computed for the HPCS. V i is in kcal/mol. Quadruplet V 1 V 2 V 3 p tor1 p cot1 X{Si{C{X 0.0000 0.6675 0.0000 -8.2352 0.0000 C{C{C{C 0.0000 38.9174 0.3649 -8.2931 -2.0127 C{C{C{H 0.0000 49.1001 0.2713 -8.5284 -1.5309 H{C{C{H 0.0000 34.0265 0.3804 -6.3917 -0.9965 H{Si{Si{H 0.0000 0.0000 0.0640 -2.4426 0.0000 H{Si{Si{Si 0.0000 0.0000 0.1587 -2.4426 0.0000 3.4 Atomistic Model of the Polymer To test the accuracy of ReaxFF, we study thermal decomposition of the HPCS. The goal is to see whether the reaction products that MD simulation using ReaxFF determines agree with what is known experimentally, described in Chapter 1. To do so, we rst generate an atomistic model of the HPCS. A modied version of the self-avoiding walk (SAW) method of Theodorou and Suter [101{103] was used to generate the initial structure of the HPCS. A cubic simulation cell was used in which three connected atoms of the HPCS' backbone were placed in random orientations. The polymer was then constructed by adding one atom at a time in a stepwise fashion, using the modied SAW algorithm. The allowed rotational states of the successive bonds between adjacent atoms were determined from the 60 probability distribution functions that were governed by energy considerations (see below). To suppress the surface eects, periodic boundary conditions were used. The polymer-consistent force eld (PCFF) was utilized in order to generate the atomistic model of the HPCS in which the total potential energyE of the polymer is given by E =E s +E +E +E nonbond ; (3.2) with E s = 4 X i=1 K i (ll 0 ) i ; E = 4 X i=2 H i ( i ) i ; E = 3 X i=1 D i (1 cosi); (3.3) and E nonbond = X i;j q i q j r ij + X ij ij " 2 ij r ij 9 3 ij r i j 6 # : (3.4) Here, E s is the energy associated with bond stretching, with l 0 being a bond's equilibrium length and l its actual length at any time during the simulations, and K i a stretching-force constant. Values of K i and the equilibrium lengths l 0 are listed in Table 3.7 for all the atoms. E represents the energy associated with the changes in the bonds' angles, where 0 is the equilibrium angle of a bond, its angle at any given time during the simulations, and H i the corresponding force constant. Numerical values of H i and 0 are listed in Table 3.8. The contribution of the torsional forces is represented by E , with D i being a force constant, and i the dihedral angle. Values of the parameters are given in Table 3.9 for all the quadruplets. The parameters i and i are, respectively, the Lennard-Jones (LJ) size and energy parameters for atom i, and q i is the partial charge of i. The non- bond interactions were cut o at an interatomic distance of 12.5 A, while a spline 61 Table 3.7: Values of the force constants K i and the equilibrium lengths l 0 for various pairs, used with the PCFF. Pair l 0 ( A) K 2 K 3 K 4 C{C 1.417 470.836 -627.618 1327.635 C{H 1.101 345.000 -691.890 844.60 C{Si 1.899 189.653 -279.421 307.513 Si{H 1.478 202.779 -305.360 280.268 Table 3.8: Values of the force constantsH i and the equilibrium angles 0 for various triplets, used with the PCFF. Triplet 0 (deg) H 2 H 3 H 4 C{Si{C 113.185 36.206 -20.394 20.017 C{Si{H 112.098 36.483 -12.809 0.000 H{C{H 107.660 39.641 -12.921 -2.432 H{Si{H 108.605 32.541 -8.316 0.000 S{C{H 112.035 28.772 -13.952 0.000 switching function was used between 9.5 and 12.5 A. The cutos were selected care- fully, based on extensive preliminary simulations, as well as our previous extensive experience with modeling of Coulombic interactions using the multipole expan- sion [104] and the Ewald summation [105, 106] techniques. Other methods of simulating Coulombic interactions are also available [107]. It would perhaps be somewhat more accurate to also include corrections due to the introduction of the cutos, but for the main purpose of the present work, namely, creating the initial atomistic model of the HPCS to be used with ReaxFF for modeling its pyrolysis, the cutos without any corrections seemed to suce (see below). The parameters associated with E nonbond are listed in Table 3.10. We emphasize that the PCFF was used only to generate the initial structure of the HPCS. To obtain the equilibrated atomistic model of the polymer, energy minimization and MD simulations were utilized. The density of the polymer at the beginning 62 Table 3.9: Values of the force constants D i for the quadruplets, used with the PCFF. Quadruplet D 1 D 2 D 3 C{Si{C{H 0.00 0.00 -0.0657 H{C{Si{H 0.00 0.00 -0.0657 Table 3.10: Values of the size and energy parameters of the non-bond part of the PCFF. Atom i ( A) i (kcal/mol) C 4.010 0.054 H 2.995 0.020 Si 4.450 0.190 of the energy minimization was set to be 0.15 g/cm 3 , so as to avoid ring spearings and catenations. The total energyE of the polymer was then minimized using the conjugate-gradient method [108]. The resulting structure was then compressed using MD simulation in the (NPT ) ensemble at pressures between 0.5 and 0.6 GPa, in order to increase the polymer density and bring it to as close to the experimental value of 0.98 g/cm 3 as possible. Energy minimization and short- duration MD simulations in the (NVT ) ensemble at 1000 K were then followed, in order to further relax the material. Then, MD simulations were carried out in the (NPT ) ensemble for several ns at 1 atm and 300 K. The time step was 1 fs. The pressure was controlled by the Andersen method [109], while the temperature was held xed using the Nos e-Hoover thermostat [110,111]. We used W = 20 amu for masslike parameter of the Andersen barostat. For the Nos e-Hoover thermostat, the parameter that was used to scale the ctitious mass was 1. Figure 3.5 presents the change with the time of the potential energy of the HPCS at 300 K during MD simulations in the (NPT ) ensemble. We note that 63 Figure 3.5: Time-dependence of the potential energy of the HPCS during (NPT )- MD simulation at 300 K. The force eld used was the PCFF. while the constant value of the total energy of the polymer may be a good indi- cation of whether it has attained its true equilibrium conguration, it might not be sucient under certain conditions and, therefore, we also computed its radial distribution functiong(r). Figure 3.6 presents the results at 300 K. The computed g(r), and in particular the locations of the various peaks that correspond to the three types of the covalent bonds in the polymer, are in very good agreement with the experimental data [23], hence conrming the accuracy of the equilibrium con- guration of the HPCS. Figure 3.7 presents the structure of the HPCS at the end of MD simulation, which is composed of 1928 atoms within a cell of dimensions of 28:9 28:9 28:9 ( A) 3 . The cell contained 4 separate chains, each composed of 482 atoms. The computed density of the polymer is 0.973 gr/cm 3 , which should be com- pared with the experimental value of 0.980 gr/cm 3 . To further check whether the 64 Figure 3.6: The computed radial distribution function of the HPCS. Figure 3.7: The amorphous structure of the HPCS, as generated by the PCFF. generated polymer structure is accurate, we also computed the bond distances, as well as the triple and dihedral angles that exist in the polymer. The results for 65 Si{C bonds, Si{C{Si angles, and H{C{Si{C dihedral are presented in Figure 3.8. The most probable Si{C distance was computed to be 1.9 A, very close to the actual distance of 1.89 A. The most probable angle in the Si{C{Si triplets is 109 , which is in the expected range for such triplets. As for the dihedral angles, the most probable angles for the cis and trans congurations were computed to be 58 and 180 , respectively, which are again very close to the theoretical values of 60 and 180 . To our knowledge, there are no experimental data for the angles in our polymer. Similar results were obtained for other bond distances, as well as triplet and dihedral angles and are given in Appendix A. The atomistic model of the HPCS was then utilized with ReaxFF to carry out MD simulation of its thermal decomposition. Before heating the polymer up, energy minimization and pre-equilibration simulations using ReaxFF were utilized to further relax the atomistic structure of the polymer. All the simulations in this section were implemented using the MD simulator GRASP { General Reactive Atomistic Simulation Program { which can be utilized in parallel computations. The total energyE of the polymer was minimized with relaxation in all the direc- tions. The maximum displacement allowed in any direction was 0.01 A, while the time step was 0.1 fs. Figure 3.9 presents the change with the time of the total energy and pressure of the system during minimization. It is clear that the energy of the polymer has reached its minimum after 0.025 ps. To further check whether the true equilib- rium was reached, the atoms were numbered and their energies were monitored in the last time frame of the minimization simulation. The results are shown in Figure 3.10, indicating that all the energies are constant. To further equilibrate the polymer structure after energy minimization, MD simulation was carried out in the (NVT ) ensemble for 5 ps at 50 K. The minimum absolute temperature 66 0 1 2 3 4 5 1.6 1.7 1.8 1.9 2 2.1 2.2 2.3 2.4 Probability (1/Å) Si−C Bond Length (Å) (a) 0 0.01 0.02 0.03 0.04 60 80 100 120 140 160 Probability (1/degrees) Si-C-Si Angle (degrees) (b) 0 0.001 0.002 0.003 0.004 0.005 -200 -150 -100 -50 0 50 100 150 200 Probability (1/degrees) H-C-Si-C Angle (degrees) (c) Figure 3.8: Distribution of (a) Si{C distances; (b) Si{C{Si angles, and (c) H{C{ Si{C dihedral angles in the HPCS polymer, as generated by the PCFF. 67 Figure 3.9: Time-dependence of the energy of the HPCS and the system's pressure during minimization using the ReaxFF. deviation for which the rescaling of the velocities was performed was 10 K, and after each 10 time steps the temperature was compared with its target value. The rescaling coecient was 0.7, so selected carefully to rescale the velocities to attain the target temperature. The dependence on time of the total energy, temperature, and pressure of the system is shown in Figure 3.11, which indicate that temper- ature was kept close to its target value, while the energy reached its equilibrium value after a few ps. In addition, Figure 3.12 indicates that the energies of each type of the atoms did not deviate abruptly from their average values. 3.5 Cook-o Simulation After determining the parameters of ReaxFF and generating a molecular model of the HPCS, cook-o simulation of the polymer was carried out to obtain an 68 -300 -250 -200 -150 -100 -50 0 50 0 482 964 1446 1928 Energy (kcal/mol) Atom Number Si C H Figure 3.10: Energies of all the atoms of the polymer at the end of minimization, using ReaxFF. overview of its thermal degradation as a result of heating it up. Molecular dynamics simulations were carried out in the (NVT ) ensemble using a heating rate of 0.1 K/fs, during which the temperature was increased from 50 K to 5000 K. The variations of the pressure during the simulation are shown in Figure 3.13. The parameters of the simulations were the same as before. For temperatures less than but close to 1000 K, the Si{H and C{H bonds were very mobile, and continuously formed and broke before the initiation of the polymer decomposition. The simulations indicated that the degradation of the HPCS is initiated at about 1000 K by releasing a few H atoms from the Si{H and C{H bonds; see Figure 3.14. The Si{H bond is about 20 kcal/mol weaker than the C{H bond. Therefore, its cleavage initiated the decomposition process and happened more frequently, consistent with the experimental results [23]. At higher temperatures, extensive hemolytic Si{H and C{H bond breaking occurred, 69 -178 -177.9 -177.8 -177.7 -177.6 -177.5 0 1 2 3 4 5 Total Energy (kcal/mol/1000) Time (ps) 20 30 40 50 60 70 0 1 2 3 4 5 Temperature (K) Time (ps) -1 -0.5 0 0.5 1 0 1 2 3 4 5 Pressure (GPa) Time (ps) Figure 3.11: Time-dependence of the energy and temperature of the HPCS, and the system's pressure during pre-equilibration at 50 K, using ReaxFF. which led to the formation of a large number of H radicals, shown in Figure 3.15. A secondary reaction involving bond formation between the H radicals led to the formation of hydrogen gas, the main pyrolysis product. The loss of H from SiH n (n = 1; 2) led to the formation of reactive silylene (=Si:) intermediates, resulting in the formation of the Si{Si bonds due to the initial chain cross-linking. But, above 1000 K it quickly broke because there is tendency for Si atoms to form more stable bonds with C, hence explaining the large uctuations for the number of the Si{Si bonds shown in Figure 3.14. Cleavage of the Si{C bonds initiated at about 2000 K and accelerated at 2400 K and higher, which led to the formation of CH 2 , SiH 2 and CH 4 Si radicals, shown in Figure 3.16. They reacted with the H radicals to form CH 3 , SiH 3 and CH 5 Si radicals, and eventually methane, silane, and methylsilane gases. In addition, 70 -300 -250 -200 -150 -100 -50 0 50 0 482 964 1446 1928 Energy (kcal/mol) Atom Number Si C H Figure 3.12: Energies of all atoms of each type at the end of pre-equilibration (at 50 K), using ReaxFF. various C 2 hydrocarbons, including C 2 H 6 and C 2 H 4 , and such radicals as C 2 H 5 and C 2 H 3 , as well as C 2 H 2 were rst formed, as shown in Figure 3.17. The time evolution of the gaseous compounds is shown in Figure 3.18. Table 3.11 lists all the products that are obtained after heating up the HPCS to 3600 K. Analysis of the gaseous products by mass spectroscopy and gas chromatog- raphy has indicated that the simulation results are consistent, both quantitatively and qualitatively, with the experimental data [23]. Reference [23] reported that most of the produced gas was hydrogen, followed in much smaller quantities by methane, methyl silanes and C 2 hydrocarbons. These are also what we obtained, as well as other gases that may not be possible to measure experimentally with any precision (due to their trace amounts). In fact, this is an advantage that ReaxFF oers. 71 Table 3.11: The pyrolysis products after heating up the polymer to 3600 K. N is the number of the produced product, and MW its molecular weight. N Product MW N Product MW 122 H 1.00 1 C 2 H 5 Si 3 113.32 26 H 2 2.00 1 C 2 H 7 Si 3 115.33 3 CH 3 15.03 3 C 2 H 10 Si 3 118.36 11 CH 4 16.04 1 C 2 H 11 Si 3 119.36 2 C 2 H 4 28.05 1 C 5 H 9 Si 2 125.29 3 HSi 29.09 3 C 3 H 8 Si 3 128.35 7 H 2 Si 30.10 1 C 3 H 9 Si 3 129.36 5 H 3 Si 31.11 2 C 3 H 10 Si 3 130.37 10 H 4 Si 32.12 1 C 4 H 12 Si 3 144.39 1 CSi 40.09 1 C 4 H 13 Si 3 145.40 4 CH 2 Si 42.11 1 C 2 H 12 Si 4 148.46 14 CH 3 Si 43.12 1 C 3 H 10 Si 4 158.45 16 CH 4 Si 44.13 1 C 3 H 11 Si 4 159.46 10 CH 5 Si 45.13 1 C 4 H 12 Si 4 172.48 5 CH 6 Si 46.14 1 C 4 H 13 Si 4 173.49 1 CH 7 Si 47.15 1 C 5 H 14 Si 4 186.50 1 C 2 H 2 Si 54.12 1 C 3 H 15 Si 5 191.58 1 C 2 H 3 Si 55.13 1 C 4 H 13 Si 5 201.57 3 C 2 H 5 Si 57.14 1 C 4 H 14 Si 5 202.58 3 C 2 H 6 Si 58.15 1 C 4 H 15 Si 5 203.59 1 C 2 H 7 Si 59.16 1 C 5 H 15 Si 5 215.60 1 CH 4 Si 2 72.21 1 C 5 H 19 Si 5 219.63 2 CH 5 Si 2 73.22 1 C 6 H 16 Si 5 228.62 4 CH 6 Si 2 74.23 2 C 7 H 16 Si 5 240.63 1 CH 7 Si 2 75.24 1 C 5 H 15 Si 6 243.69 3 C 2 H 5 Si 2 85.23 1 C 6 H 19 Si 6 259.73 3 C 2 H 6 Si 2 86.24 1 C 6 H 19 Si 7 287.81 4 C 2 H 7 Si 2 87.25 1 C 8 H 25 Si 7 317.88 1 C 2 H 8 Si 2 88.25 1 C 9 H 20 Si 8 352.94 2 C 2 H 9 Si 2 89.26 1 C 14 H 37 Si 12 542.47 1 C 2 H 10 Si 2 90.27 1 C 14 H 39 Si 13 572.57 1 C 3 H 8 Si 2 100.27 1 C 17 H 42 Si 13 611.63 2 C 3 H 9 Si 2 101.27 72 Figure 3.13: Time-dependence of the temperature and pressure during the cook-o simulation of the HPCS, using ReaxFF. 73 0 20 40 60 80 100 120 0 10 20 30 40 50 Number of Bonds Time (ps) (a) C-C 700 800 900 1000 1100 1200 1300 1400 0 10 20 30 40 50 Number of Bonds Time (ps) (b) C-H 0 20 40 60 80 100 120 140 0 10 20 30 40 50 Number of Bonds Time (ps) (c) H-H 700 800 900 1000 1100 1200 1300 1400 0 10 20 30 40 50 Number of Bonds Time (ps) (d) Si-H 700 800 900 1000 1100 1200 1300 1400 0 10 20 30 40 50 Number of Bonds Time (ps) (e) Si-C 0 10 20 30 40 50 0 10 20 30 40 50 Number of Bonds Time (ps) (f) Si-Si Figure 3.14: The evolution of the number of bonds during the cook-o simulation of the HPCS, using ReaxFF. 74 Figure 3.15: The evolution of the number of the hydrogen radicals during the cook-o simulation of the HPCS, using ReaxFF. 75 0 5 10 15 20 25 0 10 20 30 40 50 0 1000 2000 3000 4000 5000 Number of Radicals Temperature (K) Time (ps) (a) CH 2 CH 3 Temperature (K) 0 5 10 15 20 25 0 10 20 30 40 50 0 1000 2000 3000 4000 5000 Number of Radicals Temperature (K) Time (ps) (b) SiH 2 SiH 3 Temperature (K) 0 5 10 15 20 25 30 35 0 10 20 30 40 50 0 1000 2000 3000 4000 5000 Number of Radicals Temperature (K) Time (ps) (c) CH 4 Si CH 5 Si Temperature (K) Figure 3.16: The evolution of the number of the CH n (top), SiH n (middle) and CSiH n (bottom) radicals during the cook-o simulation of the HPCS, using ReaxFF. 76 Figure 3.17: The evolution of the C 2 hydrocarbons during the cook-o simulation of the HPCS, using ReaxFF. 77 0 10 20 30 40 50 60 25 30 35 40 45 50 2500 3000 3500 4000 4500 5000 Number of Molecules Temperature (K) Time (ps) H 2 CH 4 SiH 4 CH 6 Si C 2 Hydrocarbons Temperature (K) Figure 3.18: The evolution of various compounds during cook-o simulation of the HPCS, using ReaxFF. 78 Chapter 4 Amorphous SiC Formation via Reactive Dynamics Simulation of the HPCS 4.1 Introduction In Chapter 3 we described the development of the ReaxFF reactive force eld, in order to study the pyrolysis of hydridopolycarbosilane (HPCS), [SiH 2 CH 2 ] n . As described in that chapter, the pyrolysis process is an essential step in fabrica- tion of silicon-carbide (SiC) nanoporous membranes. In the present chapter we utilize ReaxFF with reactive dynamics (RD) simulation to study the temperature- dependence of the decomposition of the HPCS [112]. 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 [112]. This pro- vides 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 organization of this chapter is as follows. Section 4.2 brie y describes the generation of an atomistic model of the HPCS. How the polymer model is utilized in the molecular dynamics (MD) simulations of its pyrolysis is described in Section 4.3. Molecular dynamics simulation of thermal decomposition of the HPCS is described in Section 4.4, where the results are presented, discussed, and 79 compared with the experimental data. Section 4.5 describes how the generation of H 2 gas is used to nd the kinetics of the HPCS decomposition. In Section 4.6, 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 described. Section 4.7 compares the properties of the resulting SiC ceramic with the experimental data. 4.2 Atomistic Model of the Polymer As the rst step we generate an atomistic model of the HPCS, in order to carry out MD simulations of its pyrolysis using ReaxFF, and study the eect of pyrolysis on the nal atomistic model and composition of the SiC material. The molecular model of the HPCS was generated using the polymer-consistent force eld (PCFF). The procedure for the generation of the polymer's molecular model, as well as the parameters of the PCFF were given in Section 3.4. 4.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 eld to describe the thermal degradation process over a wide range of temperature. The 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-o simulations of the polymer was described in Section 3.5. 80 Before heating up the polymer, energy minimization and pre-equilibration sim- ulations 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 devel- oped by Sandia National Laboratories. The GRASP simulator is written in FOR- TRAN language and can be executed in parallel computations. The total energyE 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 A, 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, similar to Figure 3.10. Following energy minimization, we equilibrated the structure using isothermal- isochoric (NVT )-MD simulations for 5 ps at 50 K, during which temperature was rst rescaled to 50 K and then held xed. 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 temperature was checked and compared with the target value. The rescaling coecient (which lies between 0 and 1) was 0.7, so selected carefully to rescale the velocities to obtain the target temperature. Similar to minimization step, the constancy of all atom energies were observed (see Figure 3.12). 81 4.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 rst 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{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 K, 3200 K, and 3400 K for which the simulations lasted 60 ps. The time step, the minimum absolute temperature deviation, and the rescaling coecient were the same as before. We then computed the temperature-dependence of the populations of the types of the bonds formed in the nal conguration that were obtained from the (NVT )- MD simulation after 20 ps. The results for two temperature intervals are shown in Figures 4.1 and 4.2. Although the pyrolysis in the fabrication of SiC nanoporous membranes is at about 1000 K [5, 6], our RMD simulations indicated that, below 1800 K, there were very small changes in the numbers of the bonds in nanosecond time scale, pratical 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 nd 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 4.2 indicates, above 2400 K the RMD leads to breaking of the Si{C bonds, which leads to the formation of CH 2 and SiH 2 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 82 1200 1220 1240 1260 1280 1300 Bulk 300 600 900 1200 1500 1800 Number of Bonds Temperature (K) SiC SiH CH Figure 4.1: Population of the bonds in the nal conguration obtained from the (NVT )-MD simulations. Also given is the population for the HPCS in the bulk before raising the temperature. 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 C 2 hydrocarbons and hydrogen molecules, which agrees with the experimental data. Analysis of the experimental data for the gaseous products of pyrolysis of the HPCS indicated [23] that H 2 is the major product of the pyrolysis, while methane, silane, and methyl silanes gases are also produced in small amounts. The nal conguration of the system at 3600 K, obtained by (NVT )-RMD simulation, is shown in Figure 4.3. The complete list of the products at the same temperature was provided in Table 3.11. Figure 4.4 compares the populations of the gaseous products in the nal cong- uration after RMD simulations for 20 ps, in the range 2400{4000 K. This indicates 83 900 950 1000 1050 1100 1150 1200 1250 1300 Bulk 2400 2800 3200 3600 4000 Number of Bonds Temperature (K) SiH 900 950 1000 1050 1100 1150 1200 1250 1300 Bulk 2400 2800 3200 3600 4000 Number of Bonds Temperature (K) SiC 0 10 20 30 40 50 2400 2800 3200 3600 4000 Number of Bonds Temperature (K) SiSi 900 950 1000 1050 1100 1150 1200 1250 1300 Bulk 2400 2800 3200 3600 4000 Number of Bonds Temperature (K) CH 0 20 40 60 80 100 2400 2800 3200 3600 4000 Number of Bonds Temperature (K) CC 0 20 40 60 80 100 2400 2800 3200 3600 4000 Number of Bonds Temperature (K) HH Figure 4.2: Population of the bonds in the nal conguration obtained from the (NVT )-MD simulations, for the temperature range 2400{4000 K. Also given is the population for the HPCS in the bulk before raising the temperature. that H 2 is the major product of the pyrolysis. Figures 4.5 and 4.6 present the pop- ulations of the CH n and SiH n (n=1{4) radicals and molecules in the same system. Figure 4.7, which presents temperature-dependence of the number of H radicals, indicates that their population is a strong function of the temperature and, due to their large number at and above 3200 K, H 2 molecules begin to form. 84 Figure 4.3: The nal conguration of the material at 3600 K, obtained by (NVT )- MD simulation. The density is 0.98 g/cm 3 . To understand better how the populations of hydrogen radicals and hydrogen molecules H 2 form over time, we present in Figure 4.8 time- and temperature- dependence of their populations. Thus, one may view each pair of hydrogen rad- icals as being equivalent to one H 2 molecule, because over longer periods of time (not currently practical in RMD simulations) the hydrogen atoms in high energy states bond together and produce H 2 . 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 specied temperature. Similarly, time- and temperature-dependence of the populations of other species were also computed during the (NVT )-MD simulations. For example, Figures 4.9 85 to 4.11 present the evolution of production of CH n (n=1{4) radicals and the nal methane molecules at 3200 K, 3600 K and 4000 K, respectively, while the corre- sponding results for SiH n (n=1{4) are shown in Figures 4.12 to 4.14. At 3200 K the population of CH 3 increases with the time, and produces methane by reacting with hydrogen, consistent with the increase in the number of CH 4 molecules dur- ing the same period, shown in Figure 4.9. The trends are similar at 3600 K (see Figure 4.10). The eect 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 CH 3 radicals and CH 4 molecules begin decreasing, as indicated by Figure 4.11. Figures 4.12 to 4.14 indicate similar trends for SiH 3 radicals and silane molecules. As for the CH n Si (n = 1{6) radicals and molecules, Figure 4.15 presents the results at 4000 K. Increasing temperature increases cleavage of Si{H and C{H bonds and, thus, more CHSi and CH 2 Si radicals are produced. There is a com- petition between the CH n Si (n > 2) radicals and breaking of the bonds in Si{H and C{H, resulting in large uctuations in the population of the radicals around a mean value. This is clearly indicated by the results shown in Figure 4.15. 4.5 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 H 2 as the primary product of the pyrolysis. Therefore, we utilized production of H 2 as a basis for studying the kinetics of the HPCS decomposition. Temperature-dependence of the ratek of production of H 2 in the interval 3000{4000 K was computed. The results are presented in Figure 4.16. Three trial runs were used to obtain a measure of 86 0 5 10 15 20 25 30 35 40 2400 2800 3200 3600 4000 Number of Molecules Temperature (K) SiH 4 CH 4 H 2 CH 6 Si Figure 4.4: Populations of gaseous molecules in the nal congurations and their dependence on temperature. the errors in the estimated values. The results were then tted to the classical Arrhenius functional form, k(T ) =k 0 exp E a RT ; (4.1) where E a is the activation energy, and R is the gas constant. We obtained, k 0 ' 1:28 10 12 (molL 1 s 1 ), and E a ' 39:7 (kcal/mol). See Appendix B for the change of H 2 population with time for all temperatures and trials and also the details of the procedure for tting and smoothing data. 87 0 2 4 6 8 10 12 14 16 18 2800 3200 3600 4000 Number of Radicals/Molecules Temperature (K) CH 1 CH 2 CH 3 CH 4 Figure 4.5: Populations of CH n radicals and gaseous molecules in the nal cong- urations and their dependence on temperature. 4.6 Formation of Amorphous SiC The nal step of modeling the process of formation of amorphous SiC is to carry out MD simulations under the conditions that mimic the experimental setup for fabrication of SiC nanoporous membranes [5,6]. 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 valance angle and torsion potential during the simulation of formation of amorphous/crystalline SiC, it is necessary to use short 88 0 2 4 6 8 10 12 14 2800 3200 3600 4000 Number of Radicals/Molecules Temperature (K) SiH 1 SiH 2 SiH 3 SiH 4 Figure 4.6: Populations of SiH n radicals and gaseous molecules in the nal cong- urations and their dependence on temperature. 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 on the pores' surface of the membrane's support, might be slightly aected with varying the pyrolysis temperature. On the other hand, experimental data [23] indicate that for temperatures below 1175 K the SiC ceramic that is produced by the pyrolysis process is in the amor- phous state, in which Si{H and C{H bonds break readily, and cross-linking between hydrogen radicals result in the production of H 2 gas, the main product of the pyrol- ysis. Below 1175 K experiments show [21] that the breaking of Si{C bonds is very limited. As the populations of the bonds and the analysis of the gaseous products 89 0 20 40 60 80 100 120 140 160 180 2000 2400 2800 3200 3600 4000 Number of H Radicals Temperature (K) Figure 4.7: Populations of hydrogen radicals in the nal congurations and their dependence on temperature. over a broad range of temperatures revealed (see above and Section 3.5), 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 nd 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 speak- ing, 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 membranes [5, 6] 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 90 0 20 40 60 80 100 120 140 0 10 20 30 40 50 60 Number of Radicals/Molecules Time (ps) H 3200K 0 5 10 15 20 25 30 35 0 10 20 30 40 50 60 Number of Radicals/Molecules Time (ps) H 2 3200K 0 20 40 60 80 100 120 140 160 180 0 5 10 15 20 Number of Radicals/Molecules Time (ps) H 3600K 0 5 10 15 20 25 30 35 40 0 5 10 15 20 Number of Radicals/Molecules Time (ps) H 2 3600K 0 50 100 150 200 250 0 5 10 15 20 Number of Radicals/Molecules Time (ps) H 4000K 0 10 20 30 40 50 60 70 0 5 10 15 20 Number of Radicals/Molecules Time (ps) H 2 4000K Figure 4.8: Time- and temperature-dependence of production of hydrogen radical and hydrogen molecule during (NVT )-MD simulations. 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 ows 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. 91 0 5 10 15 20 25 0 10 20 30 40 50 60 Number of Radicals/Molecules Time (ps) 3200K CH 0 5 10 15 20 25 0 10 20 30 40 50 60 Number of Radicals/Molecules Time (ps) 3200K CH 2 0 5 10 15 20 25 0 10 20 30 40 50 60 Number of Radicals/Molecules Time (ps) 3200K CH 3 0 5 10 15 20 25 0 10 20 30 40 50 60 Number of Radicals/Molecules Time (ps) 3200K CH 4 Figure 4.9: Time-dependence of production of CH n , during (NVT )-MD simulation at 3200 K. As pointed out earlier in this dissertation, gas chromatography and mass spec- troscopy have indicated [23] that H 2 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 hydro- gen 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 H 2 . Thus, every two hydrogen atoms may be viewed as being equivalent to one H 2 molecule and, therefore, deletion of the hydrogen radicals during the pyrolysis simulations is equivalent with deleting H 2 gas that would form at much later computational time scales. This is an important consideration, as it is com- putationally very expensive to carry out RMD simulations for long enough times to allow all the hydrogen radicals to form H 2 molecules. 92 0 5 10 15 20 25 0 5 10 15 20 Number of Radicals/Molecules Time (ps) 3600K CH 0 5 10 15 20 25 0 5 10 15 20 Number of Radicals/Molecules Time (ps) 3600K CH 2 0 5 10 15 20 25 0 5 10 15 20 Number of Radicals/Molecules Time (ps) 3600K CH 3 0 5 10 15 20 25 0 5 10 15 20 Number of Radicals/Molecules Time (ps) 3600K CH 4 Figure 4.10: Time-dependence of production of CH n , during (NVT )-MD simula- tion at 3600 K. 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 4.17. 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 le 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. 93 0 5 10 15 20 25 0 5 10 15 20 Number of Radicals/Molecules Time (ps) 4000K CH 0 5 10 15 20 25 0 5 10 15 20 Number of Radicals/Molecules Time (ps) 4000K CH 2 0 5 10 15 20 25 0 5 10 15 20 Number of Radicals/Molecules Time (ps) 4000K CH 3 0 5 10 15 20 25 30 0 5 10 15 20 Number of Radicals/Molecules Time (ps) 4000K CH 4 Figure 4.11: Time-dependence of production of CH n , during (NVT )-MD simula- tion at 4000 K. Therefore, extensive MD simulations were carried out in the (NVT ) ensemble at 2000 K. Following energy minimization and pre-equilibration (see above), the HPCS was heated up to 2000 K. The density of the polymer at the beginning of the simulation was 0.98 g/cm 3 , 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 4.18). The time step, minimum absolute temperature deviation, and the rescaling coecient were selected similar to those for the (NVT )-MD simulations at 50 K (see above). 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 sucient time to allow for the rearrangement of the atoms and cross-linking between Si and C atoms. All the parameters of the simulations 94 0 5 10 15 20 25 0 10 20 30 40 50 60 Number of Radicals/Molecules Time (ps) 3200K SiH 0 5 10 15 20 25 0 10 20 30 40 50 60 Number of Radicals/Molecules Time (ps) 3200K SiH 2 0 5 10 15 20 25 0 10 20 30 40 50 60 Number of Radicals/Molecules Time (ps) 3200K SiH 3 0 5 10 15 20 25 0 10 20 30 40 50 60 Number of Radicals/Molecules Time (ps) 3200K SiH 4 Figure 4.12: Time-dependence of production of SiH n , during (NVT )-MD simula- tion at 3200 K. were similar to those at 50 K. The potential energy shown in Figure 4.19 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 nal product of pyrolysis of the HPCS at 1025 K is an amorphous ceramic composed of near-stoichiometric SiC, which also contains a small percentage of hydrogen. The density of the ceramic resulted from the pyrol- ysis of the polymer has been reported [113] to be 2.2 g/cm 3 . 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/cm 3 . Therefore, the simulation cell should be compressed enough to achieve the density of the amorphous SiC ceramic. Thus, 95 0 5 10 15 20 25 0 5 10 15 20 Number of Radicals/Molecules Time (ps) 3600K SiH 0 5 10 15 20 25 0 5 10 15 20 Number of Radicals/Molecules Time (ps) 3600K SiH 2 0 5 10 15 20 25 0 5 10 15 20 Number of Radicals/Molecules Time (ps) 3600K SiH 3 0 5 10 15 20 25 0 5 10 15 20 Number of Radicals/Molecules Time (ps) 3600K SiH 4 Figure 4.13: Time-dependence of production of SiH n , during (NVT )-MD simula- tion at 3600 K. we carried out isothermal-isobaric (NPT )-MD simulations to compress the cell at 2000 K, using a time step of 0.1 fs with the Nos e-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/cm 3 , and then the simulation cell was allowed to relax and expand at ambi- ent temperature and pressure, reaching the right density. The procedure allows movement and arrangement of the atoms in the structure. Figure 4.20 presents the results from the (NPT )-MD simulation at 2000 K for the density, while Fig- ure 4.21 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 96 0 5 10 15 20 25 0 5 10 15 20 Number of Radicals/Molecules Time (ps) 4000K SiH 0 5 10 15 20 25 0 5 10 15 20 Number of Radicals/Molecules Time (ps) 4000K SiH 2 0 5 10 15 20 25 0 5 10 15 20 Number of Radicals/Molecules Time (ps) 4000K SiH 3 0 5 10 15 20 25 0 5 10 15 20 Number of Radicals/Molecules Time (ps) 4000K SiH 4 Figure 4.14: Time-dependence of production of SiH n , during (NVT )-MD simula- tion at 4000 K. 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 4.22. 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 4.23 indicates that after 60 ps of simulation the pressure has reached its ambient value, with the volume of 97 0 5 10 15 20 0 5 10 15 20 Number of Radicals/Molecules Time (ps) 4000K CHSi 0 5 10 15 20 0 5 10 15 20 Number of Radicals/Molecules Time (ps) 4000K CH 2 Si 0 5 10 15 20 25 30 35 40 45 0 5 10 15 20 Number of Radicals/Molecules Time (ps) 4000K CH 3 Si 0 5 10 15 20 25 30 35 40 45 0 5 10 15 20 Number of Radicals/Molecules Time (ps) 4000K CH 4 Si 0 5 10 15 20 0 5 10 15 20 Number of Radicals/Molecules Time (ps) 4000K CH 5 Si 0 5 10 15 20 0 5 10 15 20 Number of Radicals/Molecules Time (ps) 4000K CH 6 Si Figure 4.15: Time-dependence of production of CH n Si, during (NVT )-MD simu- lation at 4000 K. 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. The results are shown in Figure 4.24, indicating constant density, pressure and temperature. Figure 4.25 presents the structural evolution of the system after removing the H radicals, carrying out 818 ps of (NVT )-MD 98 Figure 4.16: Temperature-dependence of rate of production of H 2 . The results for three separate (NVT )-MD simulation runs with distinct initial congurations are shown, in order to provide a measure of the variations. Straight line represents t of the data from trial 1. simulations, and nally compressing, cooling down and relaxing the system at ambient condition. Subsequently, simulated annealing [114] in the (NVT ) ensemble was utilized to provide the system with the opportunity for surmounting the energetic bar- riers, in a search for conformations with energies lower than the possibly local- minimum energy identied by the energy minimization step. The system was annealed between 300 K and 3000 K. It was rst 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 o to 300 K with the rate of 0.1 K/fs (see Figure 4.26). The state of this system at this point was taken to be the true equilibrium state of the material. 99 Figure 4.17: The schematic of shell script for removing the hydrogen radicals from the system. Figure 4.18: Time-dependence of the pressure and temperature during MD simu- lations for heating up four HPCS strands at a density of 0.98 g/cm 3 . 100 -86 -84 -82 -80 -78 -76 -74 -72 0 100 200 300 400 500 600 700 800 Potential Energy (kcal/mol/1000) Time (ps) -6 -4 -2 0 2 4 6 0 100 200 300 400 500 600 700 800 Pressure (GPa) Time (ps) 1900 1950 2000 2050 2100 0 100 200 300 400 500 600 700 800 Temperature (K) Time (ps) Figure 4.19: Time-dependence of the potential energy and pressure during (NVT )- MD simulations. The target temperature is 2000 K. 4.7 Further Comparison with the Experimental Data In addition to the comparisons with the experimental data and knowledge described above, we analyzed the nal SiC ceramic and computed its radial dis- tribution function, g(r), in order to compare it with the experimental data. Fig- ure 4.27 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 experimen- tal data [115]. 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 101 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 0 10 20 30 40 50 60 70 80 Density (g/cm 3 ) Time (ps) Figure 4.20: The increase in the density of the system during the (NPT )-MD simulation to compress the material. -4 -2 0 2 4 6 8 10 12 14 0 10 20 30 40 50 60 70 80 Pressure (GPa) Time (ps) 1600 1700 1800 1900 2000 2100 2200 2300 2400 0 10 20 30 40 50 60 70 80 Temperature (K) Time (ps) Figure 4.21: The change in the temperature and pressure of the system during the (NPT )-MD simulation to compress the material. The target temperature is 2000 K. and Si{Si covalent bonds are consistent with experimental data [116]. These results conrm the accuracy of the MD simulation and, hence, ReaxFF, in producing the true arrangement of Si and C atoms in the material. 102 2.2 2.25 2.3 2.35 2.4 0 2 4 6 8 10 12 14 16 Density (g/cm 3 ) Time (ps) 0 2 4 6 8 10 12 14 0 2 4 6 8 10 12 14 16 Pressure (GPa) Time (ps) 200 400 600 800 1000 1200 1400 1600 1800 2000 2200 0 2 4 6 8 10 12 14 16 Temperature (K) Time (ps) Figure 4.22: Time-dependence of the density, T and P during (NVT )-MD simu- lation to lower the temperature from 2000 K to 300 K. The radial distribution function may also be used to gain information about the coordination number Z of each atom in the material [117]. The relation between Z and g(r) is given by, Z = 4 Z rm 0 g(r)r 2 dr; (4.2) where r m is the rst 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 A 3 ), and r is the radial distance (in A). 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 nal SiC ceramic produced by the simulation indicated that about 48.1% of the Si atoms and 70.0% of the C atoms 103 2.1 2.14 2.18 2.22 2.26 2.3 2.34 2.38 2.42 2.46 2.5 0 10 20 30 40 50 60 Density (g/cm 3 ) Time (ps) -4 -2 0 2 4 6 8 10 12 0 10 20 30 40 50 60 Pressure (GPa) Time (ps) 270 280 290 300 310 320 330 340 0 10 20 30 40 50 60 Temperature (K) Time (ps) Figure 4.23: 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. have a coordination number of 4. The rest of atoms have a coordination number that is mainly 3. Experimental data [2,18,20,23], 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 rst 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, 104 2.1 2.14 2.18 2.22 2.26 2.3 0 10 20 30 40 50 60 Density (g/cm 3 ) Time (ps) -4 -3 -2 -1 0 1 2 3 4 5 0 10 20 30 40 50 60 Pressure (GPa) Time (ps) 270 280 290 300 310 320 330 340 0 10 20 30 40 50 60 Temperature (K) Time (ps) Figure 4.24: The change in the density and pressure during (NPT )-MD simulation for an extra 60 ps of MD simulation. The temperature is kept at 300 K. the carbon atoms were saturated with hydrogen, followed by minimization and 10 ps long (NVT )-MD simulation at room 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 nal structure contained 29.3 wt% C, 68.4 wt% Si, and 2.3 wt% H, a composition in agreement with the experimental data [18, 23] 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. 105 Figure 4.25: 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. 0 500 1000 1500 2000 2500 3000 3500 0 5 10 15 20 25 30 35 40 Temperature (K) Time (ps) -6 -4 -2 0 2 4 6 8 10 0 5 10 15 20 25 30 35 40 Pressure (GPa) Time (ps) Figure 4.26: Variations of the properties of the system during further energy min- imization, using simulated annealing. 106 Figure 4.27: The computed radial distribution functiong(r) of the SiC ceramic, as well as those for the individual bonds in the material, obtained by MD simulation. At the end of the simulations, the structure of the SiC ceramic was validated by computing the x-ray diraction (XRD) pattern of the material. The results are shown in Figure 4.28 and compared with the experimental data [2]. 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 conrming the accuracy of the results and, hence, the ReaxFF force eld. The next chapters use the predicted SiC ceramic for simulation of transport and diusion of gaseous mixtures through it, in order to test the accuracy of the simulated material in terms of its separation properties. 107 Figure 4.28: The experimental (top) and the computed (bottom) XRD pattern of the SiC, computed at the end of the nal MD simulation. The experimental data were adopted after Musumeci et al. [2]. 108 Chapter 5 Diusion of Gases in the Amorphous SiC Model 5.1 Introduction In Chapters 3 and 4, we described the development of a new model for amor- phous SiC membranes, with a pore space that consists of interconnected pores of irregular shapes, sizes, and connectivity. To test the validity of the model, we studied the evolution of the gases during the pyrolysis of the polymer precursor, hydridopolycarbosilane (HPCS), using the reactive force eld ReaxFF. We then utilized the molecular model to compute several of its properties, such as the radial distribution function and the x-ray diraction pattern. Whenever possible, we also compared the computed properties with the experimental data, which were found to be in good agreements, both quantitatively and qualitatively, with the data. Studying ow and transport of uids in conned media are currently of fun- damental and practical interest [4, 118, 119]. Its importance becomes even more signicant when one deals with nanoporous materials. The reason is threefold. First, the energy dissipation caused by friction in conned uid can change the chemistry of the process and, in addition to forming new components, phase tran- sitions as well as drastic changes in the uid's static and dynamical properties also occur, which are phenomena of fundamental importance that can be very dierent from those under the bulk conditions. Second, it is clearly of practical interest to 109 understand the mechanism of transport in the pore space of nanoporous materi- als [118{121]. Third, the morphology of the pore space can be characterized and better understood by studying transport and adsorption in the nanoporous mate- rials [118{121], which greatly in uences the transport properties of the matrix and the materials [4]. The nanoporous material of interest to us in this dissertation contains pores with an average size of about a few A. Because of such exceedingly small pore sizes, which is comparable with those of the gases that we consider here, the main resis- tance to the uid ow in nanoporous material is due to the nano- and mesopores and, therefore, the molecular interactions between the gases' molecules, as well as between them and the pores' walls, are important and cannot be ignored. Hence, one must utilize atomistic modeling and simulation to study the transport processes that take place in the pore space of such materials. A better understanding of the eect of connement on the behavior of the uids has been achieved by a large num- ber of studies on the ow and transport of uid mixtures in nanopores [4,122,123]. In recent years, our group has been developing atomistic models of various types of materials, such as carbon molecular-sieve membranes, nanotube-polymer composites, and silicon-carbide (SiC) nano-porous membranes. The models have been used in the study of sorption and transport of low molecular-weight gasses in the membranes, using molecular dynamics (MD) simulations [4, 78, 79, 124{ 129]. Our models and studies have provided not only estimates of the macroscopic properties of interest, such as the eective diusivity and solubility of various gases, but also detailed information on the mechanisms of transport processes in such materials. In this chapter we use the molecular model of amorphous SiC material that was developed in Chapter 4 to carry out extensive MD simulations and compute 110 the self-diusivities of H 2 , CO 2 and CH 4 gases in the material. In addition, the morphology of the material is characterized by its densities and the accessible free volume, which we compute with probe molecules of various sizes. Such a study is important not only on its own, as it provides further validation of the molecular model and the simulation methodologies that were developed in Chapter 4, but also provide insight into improving the performance of SiC nanoporous membranes for which the polymer precursor is the HPCS, as well as any other polymer precursor that contains Si, C, and H. The rest of this chapter is organized as follows. In Section 5.2 we present the details of the molecular modeling and computations, including modeling of the molecular interactions between the gases and the SiC model, MD simulation of gas diusion in the material, and the force eld. Sections 5.3 and 5.4 describe the computations of the free-volume distribution and the self-diusivity, respectively. Section 5.5 contains the results and discussion of the study and their implications. 5.2 Molecular Modeling and Simulation The MD simulations consist of two steps. First, we generate the molecular model of the system, consisting of the SiC model with each of H 2 , CO 2 and CH 4 gases by the methods described below, and compute the free-volume distribution of the model. Second, we carry out equilibrium MD simulations of diusion of H 2 , CO 2 and CH 4 in the material. In what follows, we explain the details of the computational techniques used in this chapter. 111 5.2.1 Gas insertion in the SiC model To compute the self-diusivity of H 2 , CH 4 and CO 2 in the SiC model, 10 molecules each gases were inserted in the SiC model separately. For convenience, we use MH2, MCH4 and MCO2 to denote the SiC material that contains H 2 , CH 4 and CO 2 , respectively. In each case, the system was constructed by adding one molecule at a time, in a stepwise fashion, into the SiC model. The molecules were ran- domly inserted in the free space of the SiC model and any physically-unacceptable contact and overlap between the gaseous molecules and atoms of the SiC model was avoided. During the construction and after the insertion of every molecule, the total energy of the system was minimized using the conjugate-gradient algo- rithms [108, 130]. This process was followed until the desired number of gas molecules were reached for each case in the SiC model. To suppress the surface eects and the small size of the system, the periodic boundary conditions were used. After generation of the systems, the energy of the structures were rst mini- mized using the conjugate-gradient algorithm until it achieved its true minimum. To further equilibrate the structures after energy minimization, MD simulations were carried out in the (NVT ) ensemble for 20 ps at 300 K. Then, another series of MD simulations was carried out in the (NPT ) ensemble for several ns at 1 atm and 300 K. The time step was 1 fs. The pressure was controlled by the Andersen method [109], while the temperature was held xed using the Nos e-Hoover thermo- stat [110,111]. We used W = 20 amu for the masslike parameter of the Andersen barostat. For the Nos e-Hoover thermostat, the parameter that was used to scale the ctitious mass was 1. 112 5.2.2 The force eld All the interaction potentials for the SiC model and its gaseous mixtures were computed based on the universal force eld (UFF) [131], which is an all atom potential in which every atom has its own parameters. The parameters of the UFF are estimated according to general rules based only on the element, its connectivity and hybridization. The potential energy is expressed as the sum of the bonded interactions and non-bonded interactions. The total energy E is given by E =E s +E +E +E ! +E vdWaals +E Coulomb ; (5.1) with E s = 1 2 S ij (ll ij ) 2 ; (5.2) E =T ijk m X n=1 C n cosn; (5.3) E =P ijkl m X n=1 C n cosn ijkl ; (5.4) E ! =G ijkl (C 0 +C 1 cos! ijkl +C 2 cos 2! ijkl ) ; (5.5) E vdWaals =V ij 2 r ij r 6 + r ij r 12 ; (5.6) E Coulomb = 332:0637 q i q j r ij : (5.7) Here, E s is the energy associated with bond stretching, with l ij being the equilib- rium or natural bond length,l its actual length at any time during the simulation, and S ij a stretching-force constant. E represents the energy associated with the changes in the angles of the bonds, where the coecients C n are chosen to satisfy the appropriate, physically justied boundary conditions where the function has a 113 minimum at an equilibrium bond angle 0 , is the angle at any given time during the simulation, and T ijk the corresponding force constant. The contribution of the torsional forces is represented by E , with P ijkl being a force constant. The values of P ijkl and the coecients C n are determined by the rotational barrier V , the periodicity of the potential, and the equilibrium angle. For a given central j{k bond, all torsions about this bond are considered, with each torsional barrier being divided by the number of torsions present about the bond. E ! represents the energy associated with the inversion, which is a one- or two-term cosine Fourier expansion, with G ijkl being a force constant, and ! ijkl is the angle between the il axis and the ijk plane. E vdWaals is a Lennard-Jones (LJ) 6-12 type energy expression for the non-bond interactions, withV ij being the well depth, andr ij the van der Waals bond length [131]. The general r ij and V ij are obtained from the homonuclear parameters through the use of combination rules r ij = (r i r j ) 1 2 ; (5.8) V ij = (V i V i ) 1 2 : (5.9) E Coulomb represents the energy associated with the electrostatic interactions, where the q i is the atomic charge [131]. In our work the potential energy for the non- bonded interactions was cut o at 6 A. We used the Ewald summation tech- nique [130] and the multiple expansion method [132] in order to account for the Coulombic interactions. 5.3 Accessible Free Volume and Solubility To compute the free-volume distribution, MD simulations in the NPT ensemble were carried out. Twenty congurations of the SiC model structures were taken 114 from the simulation trajectories after every 50 fs. The accessible free volumes were then computed using the following method based on a hard-sphere probe. First, the cubic simulation cell was rst partitioned into a three-dimensional mesh of 100 100 100 subcells. Then, a hard sphere of a given diameter was inserted as the probe at the center of each subcell, and the distance to the nearest atom of the material was calculated. If the sum of the material's atom and the van der Waals radii of the probe were smaller than the computed distance, the subcell was considered as contributing to the accessible free volume [124]. The cavity volume was then computed using a Voronoi tessellation of the space according to the algorithm of Tanemura et al. [133]. The centers of the atoms in the system were connected to each other by vectors. The vectors were then bisected perpendicularly and a large number of intersecting planes were formed, which result in construction of polyhedrons around the atoms using the algorithm. The Voronoi polyhedron around an atom represents its own available space. This can be related to the void volume of the SiC model, which was then used to study the evolution of the free volume [124]. 5.4 Estimation of the Self-diusivities In a system at equilibrium, the particles will move in accordance with the equations of motion that dene the system and, in general, will tend to diuse away from their original location. To compute the self-diusivity of H 2 , CO 2 and CH 4 in the SiC model, the systems described in Section 5.2.1 were used. The self-diusivities 115 were computed from the mean-square displacements (MSD) of the gas molecules based on the Einstein relation, D = 1 6 lim t!1 d dt N X n=1 jR(t)R(0)j 2 ; (5.10) whereD is the diusion coecient, andR(t) denotes the Cartesian position vector at time t. Thehi represents averages over all the gas molecules and over all possible time origins for the movement of the gas molecules. The limiting slope of the mean-square displacement as a function of time is then used to evaluate the self-diusion coecient of a molecule undergoing random Brownian motion in three dimensions [124]. 5.5 Results and Discussion The cell length and density of the SiC model were about 21.4 A and 2.2 g/cm 3 , respectively. The cell lengths and densities, after energy minimization and relax- ation, were about 23.7 A and 2.22 g/cm 3 for MH2, 24.2 A and 2.25 g/cm 3 for MCH4, and 23.6 A and 2.29 g/cm 3 for MCO2. Given the void space that exists in the pure SiC model (see below), their densities appear reasonable, considering the size of the cells and the mass of the gases. 5.5.1 Morphological properties Figure 5.1 presents the results of the free-volume fractions f v versus time for a probe of size = 2:5 A, averaged over 20 congurations taken from the MD simulation in the (NPT ) ensemble. As shown in the gure, the uctuations of f v were very small, which are due to long characteristic times in the SiC model. 116 Figure 5.2 presents the free volume fractionf v as a function of the probe diam- eter, represented by its LJ size parameter , for the SiC model. As the gure indicates, the values of f v decrease rapidly as the diameter increases. The value of the f v is less than 0.5% when the probe diameter reaches 4 A. This indicates that the pores of the SiC model are large enough for the gaseous molecules (such as H 2 , CH 4 , and CO 2 ) to diuse through them. As described above, the Voronoi tessellation method was used to calculate the distribution of the cavity volumes for the unit cells of the SiC model, taking all the atoms into account. Figure 5.3 presents the resulting cavity volume distributions at three dierent times. They indicate that the distributions do not vary much with time, hence indicating the stability of the entire system. It is worth noticing that there is a small but signicant fraction of large cavity volumes, as shown in the graph. Considering the size and the volume of the gaseous molecules, it is the tail of the cavity volume distributions that are responsible for gas permeation, diusion and separation, when used as a membrane or in other applications. 5.5.2 Self-diusion coecients We now present the results for the self-diusivities of H 2 , CH 4 and CO 2 . In the current SiC model, which represents an amorphous material, the accessible free volume for the diusion of small gas molecules is distributed inhomogeneously throughout the system. In such a material, the diusion coecients may vary signicantly locally. That is, they will be dierent if computed for dierent parts of the SiC model. Doolittles law explains the dependence of the diusivities on the local free-volume fraction by D =A exp B max ; (5.11) 117 1 1.5 2 2.5 3 0 200 400 600 800 1000 1200 Accessible Free Volume Fraction (%) Time (ps) SiC Figure 5.1: The accessible free volume fraction versus time with a probe diameter of 2.5 A. or a modication of it [134], whereA andB are constants, and max are, respec- tively, the packing fraction of the material and its maximum value. One may also consider the Cohen-Turnbull statistical mechanical theory dened by D =M exp N v m v f ; (5.12) where M and N are two constants, v f is a fraction of the available free volume over the volume of the involved component, and v m is the minimum free volume needed for a diusing molecule to getaway from its nearby atoms [124,135]. In order to decrease the degree of inhomogeneity for calculation of the self- diusivities, 10 molecules of each gas were inserted at random positions in the SiC model. Therefore, the eective self-diusivities that were computed from the slope 118 0 2 4 6 8 10 0 1 2 3 4 5 6 Accessible Free Volume Fraction (%) Diameter of Probe (Å) SiC Figure 5.2: Free volume percentage as a function of the probe size for the SiC model. of the plots of the MSDs versus time should be an accurate indicator of the eect of the SiC model's morphology on the diusion of small gas molecules through it. The process of random insertion of the gases in the SiC model and the calculation of the self-diusivities were repeated several times for each case, and similar results were obtained each time. The ensembles of the H 2 , CH 4 and CO 2 trajectories during the simulation are displayed in Figure 5.4, respectively. As shown in the gures, all the gases move around in the structure, showing that the pores of the SiC model are large enough for the gases to diuse through them. However, there is an important result with practical implication, which is the longer range of the H 2 motion compared to the other two gases. This means that the generated SiC model is more favorable for the diusion and transport of H 2 and, consequently, yielding favorable separation of 119 0 0.05 0.1 0.15 0.2 0.25 0 20 40 60 80 100 Distribution Cavity Volume (ų) 0.2 ns 0 0.05 0.1 0.15 0.2 0.25 0 20 40 60 80 100 Distribution Cavity Volume (ų) 0.5 ns 0 0.05 0.1 0.15 0.2 0.25 0 20 40 60 80 100 Distribution Cavity Volume (ų) 1 ns Figure 5.3: The distribution of the free volume versus time. 120 H 2 when dealing with the gaseous mixtures. The smallest range of motion belongs to CO 2 and this can be explained not only based on being a (relatively) heavy molecule, but also due to the presence of the C atoms on the pores' surface, as CO 2 does have some anity for adsorption on the pores' surface and, hence, more resistance for its diusion. To better understand the gas movement, time-dependence of the displacement of the three gases is shown in Figure 5.5. The displacement of H 2 uctuates from 1 A to 35 A, while that of CH 4 uctuates between 1 A and 10 A, and CO 2 between 1 A and 8 A. The result is consistent with what we obtained from the trajectory ensemble, and clearly indicates better diusion of H 2 in the material and a relatively restricted morphology of the model for diusion of the CO 2 . The above results are all consistent with the hopping mechanism that is often advocated to describe diusion of small gas molecules in the pores of amorphous materials. Based on the hopping mechanism, the gas molecules uctuates in a cavity until they nd a way to hop to the adjacent cavities that are not occupied by other gas molecules [124]. The amount of the uctuations and the frequency of the hops depend on material's morphology and density and dier from one cavity to another. The SiC model of this study represents an amorphous material with irregular pore shapes. Therefore, heterogenous structures give rise to increased uctuations, hence leading to more frequent hops or jumps of the gas molecules between the pores. Such uctuations then explain the larger oscillations and the jump frequencies of H 2 due to its size that is smaller than other two gases. 121 Figure 5.4: The ensemble of the trajectories of H 2 , CH 4 and CO 2 . The simulation cells have been rotated to provide a better view of trajectories. Figure 5.6 presents the plot of log(MSD) versus log(time) for the three gases at 300 K, computed in the (NPT ) ensemble. In the Fickian regime of diusion, the mean-square displacements varies proportionally with time according to lim t!1 jR(t)R(0)j 2 = 6Dt; (5.13) 122 0 5 10 15 20 25 30 0 200 400 600 800 1000 Displacement R(t) (Å) Time (ps) H 2 0 2 4 6 8 10 12 14 0 200 400 600 800 1000 Displacement R(t) (Å) Time (ps) CH 4 0 2 4 6 8 10 12 14 0 200 400 600 800 1000 Displacement R(t) (Å) Time (ps) CO 2 Figure 5.5: Displacements of three gases from their last positions. 123 Table 5.1: The computed self-diusivities of H 2 , CH 4 and CO 2 . All numbers are in m 2 /s. gas type computed D 10 8 H 2 1.81 CH 4 0.39 CO 2 0.09 and therefore, log lim t!1 jR(t)R(0)j 2 = log (6D) + log (t) ; (5.14) whereD is the diusion coecient, andR(t) denotes the Cartesian position vector at time t. As shown in the Equation 5.14, a plot of log(MSD) versus log(t) must be linear with a slope of 1. Figure 5.6 also represent such plots, and in all cases the slopes are' 1. Moreover, in all the cases diusion of the gases reaches the Fickian regime almost quickly. Figure 5.7 depicts the time-dependence of the MSD of the three gases, com- puted in the (NPT ) ensemble. The computed self-diusivities of the H 2 , CH 4 and CO 2 in the corresponding systems are listed in Table 5.1. Figure 5.8 com- pares the time-dependence of the MSD of H 2 at temperatures of 300 K and 400 K. As shown, at higher temperatures H 2 molecules diuse faster, resulting in larger diusion coecients. Although in the next chapter we compute and compare the transport properties of the SiC model for gaseous mixtures, under an external pressure drop, with the corresponding experimental properties obtained by our group, we are aware of no experimental data for the diusion of the three gases in materials similar to amorphous SiC, for which we developed a molecular model in this study. However, the computed self-diusivities provide a very important result for the practical 124 purpose of fabricating nano-porous SiC membranes for separation of H 2 from the other gases (mainly CH 4 and CO 2 ) [5,6]. As Table 5.1 indicates, the self-diusivity of H 2 is larger than both CH 4 by a factor of 5 times and CO 2 by a factor of 20, which means that the porous SiC materials is more amenable to diusion of H 2 , due to its small size, which then results in its separation from the other gases, hence further conrming the appropriate structure of the amorphous SiC and its performance as a membrane layer. However, these results can only provide a qualitative explanation of gas separation, because under experimental condition, mixtures of gases are passed through the membrane under an external pressure drop that can change the ow pattern and separation of gases from each other. In the next chapter we will describe the performance of the molecular model of amorphous SiC in terms of the transport properties of the gases under conditions similar to experiment. 125 1.5 2 2.5 3 3.5 4 4.5 0 0.5 1 1.5 2 2.5 3 log(MSD) log(time) H 2 1 1.5 2 2.5 3 3.5 0 0.5 1 1.5 2 2.5 3 log(MSD) log(time) CH 4 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3 log(MSD) log(time) CO 2 Figure 5.6: Logarithmic plot of the MSD versus time of the three gases, computed in the (NPT ) ensemble. 126 0 2 4 6 8 10 12 14 0 200 400 600 800 1000 Mean Square Displacement/1000 (Ų) Time (ps) H 2 0 0.5 1 1.5 2 2.5 0 200 400 600 800 1000 Mean Square Displacement/1000 (Ų) Time (ps) CH 4 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 200 400 600 800 1000 Mean Square Displacement/1000 (Ų) Time (ps) CO 2 Figure 5.7: The MSDs of the gases versus time computed in the (NPT ) ensemble at 300 K. 127 0 5 10 15 20 0 200 400 600 800 1000 Mean Square Displacement/1000 (Ų) Time (ps) T=300 K T=400 K Figure 5.8: The MSDs of H 2 versus time, computed in the (NPT ) ensemble at 300 K and 400 K. 128 Chapter 6 Application to Transport and Separation of Gaseous Mixtures in SiC Membranes 6.1 Introduction In previous chapter, we described the importance of using atomistic modeling and simulation to study the transport processes that take place in the pore space of the nanoporous materials such as the SiC membrane model that was generated and discussed in Chapter 4. We computed the self diusivities of H 2 , CO 2 and CH 4 gases in the SiC model and studied the morphology of it. In this chapter we use atomistic simulation to study the transport and sep- aration of two binary mixtures through the SiC membrane model. Due to the external pressure that is applied to the membrane in order to facilitate transport and separation of the gases in the SiC membrane, what we study in this chap- ter represent non-equilibrium phenomena and, therefore, require a non-equilibrium molecular dynamics (MD) method (see below). As described in Chapter 2, the simulation method that we utilize is the dual control volume-grand-canonical MD (DCV-GCMD) method, the essentials of which were already described in Chap- ter 2. Thus, in this chapter we use the method to simulate the transport and 129 separation of the binary gas mixtures in the molecular model of the developed SiC membrane. This chapter is organized as follows. In Section 6.2, the results of non- equilibrium MD simulations are presented, while the criterion for reaching the steady-state condition is described in Section 6.2.2. In Sections 6.2.3 and 6.2.4 we study the eect of the various important parameters, such as the temperature and the magnitude of the applied pressure drop on the performance of the model mem- brane. The eect of membrane's thickness on the separation factors of the model membrane, as well as the methodology for comparing the results of the computed separation factors with the membrane's thickness are described in Section 6.2.5, where the results are also compared with the experimental data. 6.2 Results and Discussion We carried out extensive DCV-GCMD simulations of two equimolar binary gas mixtures in the model membranes, namely, the H 2 /CO 2 and H 2 /CH 4 mixtures, under a variety of conditions. We carefully studied and compared the simulation results from more than 300 series of simulations that were utilized under a variety of conditions, such as various magnitudes of the applied pressure drop, temperatures, and the membrane's thicknesses in order to identify the most accurate simulation strategy that can yield the most information about and insights into the transport properties in the SiC membrane model. 6.2.1 Details of the simulations The most important contributing factors to the separation properties of a mem- brane are its thickness, as well as the magnitude of the applied pressure gradient 130 and the temperature at which the separation operation is carried out. Thus, we vary all three factors in order to study their eect. They include simulations under external pressure drops ranging from 2 atm to 10 atm with increments of 2 atm, those ranging from 10 atm to 100 atm with increments of 10 atm, simulations at 300 K and at several other temperatures ranging from 375 K to 675 K with incre- ments of 100 K, the membrane's thickness ranging from about 21.4 A to 128.4 A with increments of 21.4 A and from 128.4 A to 1284.2 A with increments of 128.4 A. The membrane's linear size in the yz plane perpendicular to the direction of the applied pressure gradient was also increased from 21.4 A by increments of 21.4 A to 64.2 A. The number of atoms in the simulations varied from about 640 to 1:710 5 atoms. Such a large number of simulation runs not only makes it possible to identify a computational strategy that yields results that can be compared with the experimental data, but also makes it feasible to study the behavior and trans- port properties of the membrane under conditions that are not easily achievable experimentally. As we describe shortly, it was required in all the simulations to increase the size of the original molecular model of the SiC membrane, which was developed in Chapter 4, in order to improve the statistics of the simulations by allowing to ow a reasonable number of gas molecules in the system and to dynamically reach the steady-state condition, where the transport properties of the gases need to be studied. The membrane layer with increasing sizes was constructed by replicating the original model along the -x, -y and -z directions until the required size in each direction and each simulation was achieved. After extensive simulations, we found that the most accurate cross section of the membrane in the yz plane perpendicu- lar to the uid ow in -x direction is 64.2 A by 64.2 A. Dierent membrane lengths along the -x direction (see above) was then constructed using the identied cross 131 section. We emphasize that one of our goals is to study how the transport prop- erties of the gases in the membrane vary as the membrane length along the uid ow direction (the -x direction) changes. In all simulations the pressure drop across the membrane was xed by setting the upstream pressure in upstream control volume (CV) to the target pressure drop and the downstream pressure to zero (see Figure 2.2). In practice this provides a better comparison with the experiment [6], where the gases exit to a large vacuum at zero pressure. 6.2.2 Steady-state condition As the rst step, it is essential to ensure that simulations are carried out correctly in the sense that, the system has reached the steady-state condition, and the pres- sure, temperature, density, and uxes in the CVs, and consequently, the separation factors are xed. Reaching the steady-state condition using the GCMD method, especially for the present application, proved to be very dicult. Experimentally, the transport properties are measured for membranes with cross sections and thick- nesses of several cm 2 andm, respectively, representing a molecularly-large system with billions of atoms. Moreover, extensive experiments must be carried out for at least several minutes before collecting any data, in order to ensure that the system has reached the steady-state condition. Therefore, it is quite clear that it would be computationally prohibitive to simulate a system with the size and time scale that are comparable with the experimental conditions. They are, in fact, well beyond the capabilities of today's supercomputers. In addition, one more point requires care, which relates to the nature of GCMD simulations. During GCMD simulations, the Monte Carlo (MC) moves occur after every several MD steps to construct the required pressures and densities in the two CVs. As explained in 132 Chapter 2, the method uses a random insertion and deletion of the gas particles in the two CVs. This, in principle causes the uid ow to experience disturbances every several time steps, as a result of which the simulations slow down while trying to reach the steady state. The solution to this problem is to have large enough CVs and to apply higher external pressures than in the experiment. The method provides a large enough number of gas molecules in the CVs, which in turn provides better statistics after averaging over specied period of times. More importantly, the disturbance caused by the MC moves is minimized. On the other hand, it should be noted that there is a limit for enlarging the system to avoid expensive computations. In addition, to reach the steady-state condition, the simulations in most of the cases were run for more than 100 ns, with the averages taken every 1 million steps. Figures 6.1 and 6.2 present the dependence of the separation factors on the time for the binary mixtures of H 2 /CO 2 and H 2 /CH 4 for several pressure drops. As shown, in all the cases, after about 30 million time steps, the separation factors attain constant values, which is an indication of reaching the steady-state condition. The density of each gaseous component in the two CVs must be held at some xed values. These values, which are in equilibrium with two bulk phases, are directly related to the required chemical potentials and pressures in the corre- sponding CVs. To obtain constant values of the densities, a sucient number of the GCMC insertions and deletions of the molecules were carried out in the CVs. The GCMC and method of computing the chemical potential and pressure in each CV were described in Chapter 2. The time-averaged dimensionless density proles i (x) for the binary gas mix- ture of H 2 /CO 2 as a function of X = x= 1 along the SiC membrane model and CVs are presented in Figures 6.3 to 6.6 at several temperatures and pressures. 133 1 2 3 4 5 6 7 8 0 10 20 30 40 50 Separation Factor Number of Steps/10 6 H2/CO2 T=475 K 20 atm 40 atm 60 atm Figure 6.1: Time-dependance of the separation factor for the H 2 /CO 2 mixture for pressures drops of 20, 40, and 60 atm at 475 K. In this and subsequent gures, the dashed lines represent the boundaries of the membrane model. Same proles are shown for the binary gas mixture of H 2 /CH 4 in Figures 6.7 to 6.10. The pressures at the upstream CVs are 20 atm, 60 atm and 100 atm. The density proles follow the same trends at the specied pressures and temperatures. The only dierences are the values of the densities that are higher at the higher P with the same T and are lower at the higher T with the same P . These dierences are expected because when T is constant andP is increased, more particles are required to be inserted in the CVs to construct the new pressure and, consequently, a higher density is obtained. On the other hand, for the same P , fewer particles are required when T is increased, because the particles move faster at higher T and, therefore, more collisions between them that contribute to the pressure increase in the CVs. 134 1 1.5 2 2.5 3 3.5 4 4.5 5 0 10 20 30 40 50 Separation Factor Number of Steps/10 6 H2/CH4 T=475 K 20 atm 40 atm 60 atm Figure 6.2: Time-dependance of the separation factor for the H 2 /CH 4 mixture for pressures drops of 20, 40, and 60 atm at 475 K. The density proles are essentially at in the upstream (left) CV, i.e. in the region25<X <8:33, and are very negligible in the downstream (right) CV, i.e. in the region 8:33<X < 25 as the pressure there is set to zero. Although due to the numerical noise there are small uctuations in the proles, the averages of the numerical values are similar to those computed with the GCMC method under the same conditions. This indicates that the chemical potentials and pressures are properly computed and maintained during the DCV-GCMD simulations. There is a sharp increase in the density of CO 2 and CH 4 at the entrance of the membrane,X =8:33, which is due to the adsorption anity of the gaseous molecules for the SiC surface. The density's peak at the entrance of the membrane is sharper for CO 2 than CH 4 , which is due to the larger energy parameter of CO 2 - membrane interactions. In other words, the SiC membrane is more attractive to 135 CO 2 than to CH 4 and, therefore, more molecules of CO 2 accumulate at the pores' surface of the membrane than do the CH 4 molecules. On the other hand, H 2 has very little anity for adsorption on the surface of SiC and, compared to CO 2 and CH 4 , the energy parameter for H 2 interactions with the membrane is low. Therefore, as shown in the gures, the accumulation of H 2 at the entrance of the membrane is smaller than the other two gases. Within the membrane (the transport region),8:33<X < 8:33, the density proles for all the components decrease along the membrane's length, which is expected. The density proles are neither linear nor smooth. The total ux for every component is the sum of the diusive and convective contributions and, due to the external pressure drop applied to the membrane, the overall streaming velocity is nonzero. Therefore, the density proles are not linear. In addition, the membrane is an amorphous material that is highly heterogeneous with pores of dierent sizes and connectivities (see Chapter 5). Thus, as a result, unlike a single pore the density proles are not smooth. They change broadly in the membrane. The temperature gradient eect on the transport properties must be elimi- nated by holding the temperature constant, when the eect of the applied external pressure gradient on the transport of a mixture is studied. Therefore, special care was taken to hold the temperature constant by rescaling independently the velocity components in all the three directions during the entire MD simulation. Figures 6.11 and 6.12 present the dimensionless temperatures T for the binary mixtures of H 2 /CO 2 and H 2 /CH 4 at two temperatures and for the applied external pressure drop of 20 atm, respectively. Each T at an axial position X presents an average over a cross section of the membrane at X . As Figures 6.11 and 6.12 136 0 5 10 15 20 25 -25 -15 -5 5 15 25 Density x 10 3 X * T=375 K ΔP=20 atm CO2 H2 Figure 6.3: The dimensionless density proles of the H 2 /CO 2 mixture in the axial (X ) direction, in the membrane (middle) and the two control volumes. The temperature is 375 K and the applied external pressure drop is 20 atm. indicate, the temperatures in the two CVs are completely constant. The uctua- tions in the middle part (membrane layer) are small and due to the pore space's heterogeneity. 6.2.3 Eect of the temperature Increasing the temperature of the molecules causes them to diuse much faster and, consequently, result in increased uxes. One should not also ignore the possibility of a small amount of adsorption of gases, if any, on the surface of SiC pores at higher temperatures. Higher uxes are obtained as a result of the combination of the two factors. Since the hydrogen is the lightest gas among other gases in the present study, higher temperatures would result in higher separation factors 137 0 5 10 15 -25 -15 -5 5 15 25 Density x 10 3 X * T=475 K ΔP=20 atm CO2 H2 Figure 6.4: The dimensionless density proles of the H 2 /CO 2 mixture in the axial (X ) direction, in the membrane (middle) and the two control volumes. The temperature is 475 K and the applied external pressure drop is 20 atm. for H 2 over both CO 2 and CH 4 . The results of the MD simulations are presented in Figure 6.13, which conrms the trends. As shown in the gure, the separation factor is always an increasing function of temperature. Figures 6.14 and 6.15 present the change of the uxes with increasing temper- ature for the mixtures of H 2 /CO 2 and H 2 /CH 4 under an external pressure drop of P = 20 atm, respectively. The decrease of the uxes with increasing the temperature is due to the fact that a fewer number of molecules are required at higher temperatures to build the same pressure, which should not be confused with the fact that molecules diuse faster at higher temperature. Instead, one should consider the ratio of the uxes (i.e. separation factors), as shown in Figures 6.14 138 0 5 10 15 20 25 30 35 40 45 -25 -15 -5 5 15 25 Density x 10 3 X * T=475 K ΔP=60 atm CO2 H2 Figure 6.5: The dimensionless density proles of the H 2 /CO 2 mixture in the axial (X ) direction, in the membrane (middle) and the two control volumes. The temperature is 475 K and the applied external pressure drop is 60 atm. and 6.15. It is clear in Figure 6.13 that the separation factor increases when the temperature rises. 6.2.4 Eect of the external pressure drop Figure 6.16 presents the dependence on the external pressure drop P on the separation factor for the SiC membrane for the two mixtures that we studied. The separation factor slightly decreases at higher pressures for the H 2 /CO 2 mixture, but it remains almost unchanged for the H 2 /CH 4 mixture. In overall the separation factors are slightly higher at the lower pressure dierences. It is due to the tight and small pores of the membrane, where increasing the external pressure drop generates a mixture that is increasingly liquid-like. Hence, the contribution of the 139 0 10 20 30 40 50 60 70 -25 -15 -5 5 15 25 Density x 10 3 X * T=475 K ΔP=100 atm CO2 H2 Figure 6.6: The dimensionless density proles of the H 2 /CO 2 mixture in the axial (X ) direction, in the membrane (middle) and the two control volumes. The temperature is 475 K and the applied external pressure drop is 100 atm. diusion of the lighter and faster H 2 to its overall ux is decreased, which results in lower H 2 uxes and, therefore, smaller separation factors. The H 2 /CO 2 mixture exhibits more or less similar trends as for H 2 /CH 4 mix- ture. The only dierence is that CO 2 has some anity for adsorption on the pores' surface due to existence of the C atoms, whereas H 2 and CH 4 do not have any sig- nicant anity for adsorption on the SiC surfaces. Hence, surface ow exists in the membrane for the H 2 /CO 2 mixture, which can also be seen through the energy interaction parameters of the gases with the pores' surface. The results of this section indicate that, in general, applying higher exter- nal pressure drop would result in lower separation factor. However, as shown in Figure 6.16 the separation factor changes slightly by varying the pressure, and 140 0 5 10 15 20 -25 -15 -5 5 15 25 Density x 10 3 X * T=375 K ΔP=20 atm CH4 H2 Figure 6.7: The dimensionless density proles of the H 2 /CH 4 mixture in the axial (X ) direction, in the membrane (middle) and the two control volumes. The temperature is 375 K and the applied external pressure drop is 20 atm. remains almost unchanged for the external pressure drops studied in this disser- tation. Therefore, one may conclude that the separation factor for the membrane in this study is not a strong function of the external pressure drop. The same behaviour of the membrane was also observed experimentally. 6.2.5 Comparison with the experimental data As the next step, we carried out a series of simulations in which the membrane's thickness was increased and the corresponding separation factorsS were computed. The experimental thickness of the membrane varies from 2 to 6 m [5, 6]. As discussed earlier, simulating a system in the range of the experimental size is computationally impossible. Therefore, one way to compare the separation factor 141 0 5 10 15 -25 -15 -5 5 15 25 Density x 10 3 X * T=475 K ΔP=20 atm CH4 H2 Figure 6.8: The dimensionless density proles of the H 2 /CH 4 mixture in the axial (X ) direction, in the membrane (middle) and the two control volumes. The temperature is 475 K and the applied external pressure drop is 20 atm. resulting from the simulations with the experimental data would be to extrapolate the available data for smaller thicknesses to the range of experimental thickness of a membrane. For this purpose one needs to know how, qualitatively, the separation factor changes with the thickness of the membrane, which means that there exits an optimal thickness such that, the separation factor will not improve, if not decline, for any thickness larger than the optimal value. In other words, the separation factor increases nonlinearly with the membrane thickness but, at the same time, during its increase it also decays to an optimal point where the optimal membrane thickness has also been achieved. Therefore, one may conclude that it would not 142 0 5 10 15 20 25 30 35 40 45 -25 -15 -5 5 15 25 Density x 10 3 X * T=475 K ΔP=60 atm CH4 H2 Figure 6.9: The dimensionless density proles of the H 2 /CH 4 mixture in the axial (X ) direction, in the membrane (middle) and the two control volumes. The temperature is 475 K and the applied external pressure drop is 60 atm. be unrealistic to t the separation factor data resulted from dierent lengths to an exponential function, S d (X) =S 1 +A exp (X=B) ; (6.1) where X is the membrane thickness, A and B are constants to be determined, and S 1 is the asymptotic separation factor that is set equal to the corresponding experimental separation factors, 93 for the H 2 /CO 2 and 104 for the H 2 /CH 4 mix- tures [5, 6]. The experimental values were measured for applied external pressure drops of up to P = 10 atm, at T = 475 K. It was also shown in the previous section, and has been also conrmed experimentally, that separation factor is not 143 0 10 20 30 40 50 60 -25 -15 -5 5 15 25 Density x 10 3 X * T=475 K ΔP=100 atm CH4 H2 Figure 6.10: The dimensionless density proles of the H 2 /CH 4 mixture in the axial (X ) direction, in the membrane (middle) and the two control volumes. The temperature is 475 K and the applied external pressure drop is 100 atm. a very strong function of P . Hence, one would expect similar or close separation factors at higher P . Figure 6.17 presents the change in the separation factor with the thickness, as well as the t of the results for the H 2 /CO 2 and H 2 /CH 4 mixtures with P = 20 atm. The thickness of membrane was increased from 64.2 A to 385.2 A, the largest thickness for which the steady state can be obtained based on a reasonable com- putational time. In all the cases, the separation factor increases with increasing the thickness of the membrane. Having obtained all the required parameters after tting the data to Equation 6.1, the optimal thickness of the membrane for each gaseous mixture is computed, which is an important result with practical implica- tion. If the maximum experimental separation factors were to remain more or less 144 0 1 2 3 4 5 6 -25 -15 -5 5 15 25 T* X * T=375 K (T*=2.51) T=475 K (T*=3.19) H2/CO2 Figure 6.11: The average dimensionless temperature T in the membrane (in the middle) and the two control volumes for the H 2 /CO 2 mixture. Numbers in the parentheses indicate the set temperatures. The applied external pressure drop is 20 atm. unchanged at higher pressure drops (as it was shown in previous section), similar plots may be obtained for other pressure drops, such as P = 30 atm and 40 atm as presented, respectively, in Figures 6.18 and 6.19 and, consequently, the optimal membrane thickness can be computed. Appendix C provides results for P = 50 atm and 60 atm. Table 6.1 summarizes the computed optimal thicknesses at dierent pressure drops, which indicate that for higher pressure drops the optimum thickness of the membrane would be also larger, which is expected and consistent with the results 145 0 1 2 3 4 5 6 -25 -15 -5 5 15 25 T* X * T=375 K (T*=2.51) T=475 K (T*=3.19) H2/CH4 Figure 6.12: The average dimensionless temperature T in the membrane (in the middle) and the two control volumes for the H 2 /CH 4 mixture. Numbers in the parentheses indicate the set temperatures. The applied external pressure drop is 20 atm. presented above. The anity of CO 2 for adsorption on the pores' surface due to existence of the C atoms causes the optimum thickness for the H 2 /CO 2 mixture to be smaller than the one for H 2 /CH 4 mixture at the same pressure drops. Another signicance of the above results is shown in Figures 6.20 to 6.22, where the results of the t to Equation 6.1 along with a t to a linear function of the form, S l (X) =M +NX are presented for each case. Clearly, the linear function, when it crosses the S 1 line, results in the minimum optimal thickness, compared with any other tting function, yet the estimated values from the linear line would 146 0 1 2 3 4 5 300 400 500 600 700 Separation Factor Temperature (K) H2/CH4 H2/CO2 ΔP=20 atm Figure 6.13: Eect of the temperature on the separation factor for the two gaseous mixtures. 0 10 20 30 300 400 500 600 700 Flux x 10 5 Temperature (K) H2 CO2 ΔP= 20 atm Figure 6.14: Eect of the temperature on the uxes for the H 2 /CO 2 gaseous mix- tures. 147 0 10 20 30 300 400 500 600 700 Flux x 10 5 Temperature (K) H2 CH4 ΔP= 20 atm Figure 6.15: Eect of the temperature on the uxes for the H 2 /CH 4 gaseous mix- tures. 0 2 4 6 8 10 0 20 40 60 80 100 Separation Factor ΔP (atm) H2/CO2 H2/CH4 X=64.2 Å T=475 K Figure 6.16: Eect of the applied external pressure drop P on the separation factors for a membrane of thickness 64.2 A. 148 Table 6.1: The predicted values of the membrane thickness for the H 2 /CO 2 and H 2 /CH 4 mixtures for the pressure drops P = 20, 30, and 40 atm. The tempera- ture is 475 K. Mixture P (atm) Optimal thickness (m) H 2 /CO 2 20 1.79 H 2 /CO 2 30 1.95 H 2 /CO 2 40 2.43 H 2 /CH 4 20 3.83 H 2 /CH 4 30 4.87 H 2 /CH 4 40 9.68 result in a reasonable value for the membrane thickness. This is another indication of the validity of the developed SiC model. 149 0 2 4 6 8 10 12 14 16 0 0.01 0.02 0.03 0.04 0.05 Separation Factor X (μm) S d (X)= S ∞ - Aexp[-X/B] S ∞ = 93.0000 A = 91.7219 B = 0.3901 H2/CO2 ΔP=20 atm S data S d 0 2 4 6 8 10 12 14 16 0 0.01 0.02 0.03 0.04 0.05 Separation Factor X (μm) S d (X)= S ∞ - Aexp[-X/B] S ∞ = 104.0000 A = 101.9310 B = 0.8360 H2/CH4 ΔP=20 atm S data S d Figure 6.17: Eect of the membrane's thickness on its separation factor for the H 2 /CO 2 (top) and H 2 /CH 4 (bottom) mixtures at T = 475 K and P = 20 atm. The data are tted to the S d (X) function. 150 0 2 4 6 8 10 12 14 16 0 0.01 0.02 0.03 0.04 0.05 Separation Factor X (μm) S d (X)= S ∞ - Aexp[-X/B] S ∞ = 93.0000 A = 91.9552 B = 0.4236 H2/CO2 ΔP=30 atm S data S d 0 2 4 6 8 10 12 14 16 0 0.01 0.02 0.03 0.04 0.05 Separation Factor X (μm) S d (X)= S ∞ - Aexp[-X/B] S ∞ = 104.0000 A = 101.7440 B = 1.0627 H2/CH4 ΔP=30 atm S data S d Figure 6.18: Eect of the membrane's thickness on its separation factor for the H 2 /CO 2 (top) and H 2 /CH 4 (bottom) mixtures at T = 475 K and P = 30 atm. The data are tted to the S d (X) function. 151 0 2 4 6 8 10 12 14 16 0 0.01 0.02 0.03 0.04 0.05 Separation Factor X (μm) S d (X)= S ∞ - Aexp[-X/B] S ∞ = 93.0000 A = 91.8586 B = 0.5290 H2/CO2 ΔP=40 atm S data S d 0 2 4 6 8 10 12 14 16 0 0.01 0.02 0.03 0.04 0.05 Separation Factor X (μm) S d (X)= S ∞ - Aexp[-X/B] S ∞ = 104.0000 A = 101.4411 B = 2.1130 H2/CH4 ΔP=40 atm S data S d Figure 6.19: Eect of the membrane's thickness on its separation factor for the H 2 /CO 2 (top) and H 2 /CH 4 (bottom) mixtures at T = 475 K and P = 40 atm. The data are tted to the S d (X) function. 152 0 20 40 60 80 100 120 140 0 0.5 1 1.5 2 2.5 Separation Factor X (μm) S d (X)= S ∞ - Aexp[-X/B] S ∞ = 93.0000 A = 91.7219 B = 0.3901 S l (X)= M + NX M = 1.3994 N = 221.0276 H2/CO2 ΔP=20 atm S data S d S l S ∞ 0 20 40 60 80 100 120 140 0 1 2 3 4 5 Separation Factor X (μm) S d (X)= S ∞ - Aexp[-X/B] S ∞ = 104.0000 A = 101.9310 B = 0.8360 S l (X)= M + NX M = 2.0987 N = 118.4696 H2/CH4 ΔP=20 atm S data S d S l S ∞ Figure 6.20: The tting functions, S d (X) and S l (X), approaching the maximum separation factor for the H 2 /CO 2 (top) and H 2 /CH 4 (bottom) mixtures atT = 475 K and P = 20 atm. 153 0 20 40 60 80 100 120 140 0 0.5 1 1.5 2 2.5 Separation Factor X (μm) S d (X)= S ∞ - Aexp[-X/B] S ∞ = 93.0000 A = 91.9552 B = 0.4236 S l (X)= M + NX M = 1.1519 N = 204.9233 H2/CO2 ΔP=30 atm S data S d S l S ∞ 0 20 40 60 80 100 120 140 0 1 2 3 4 5 Separation Factor X (μm) S d (X)= S ∞ - Aexp[-X/B] S ∞ = 104.0000 A = 101.7440 B = 1.0627 S l (X)= M + NX M = 2.2770 N = 93.4779 H2/CH4 ΔP=30 atm S data S d S l S ∞ Figure 6.21: The tting functions, S d (X) and S l (X), approaching the maximum separation factor for the H 2 /CO 2 (top) and H 2 /CH 4 (bottom) mixtures atT = 475 K and P = 30 atm. 154 0 20 40 60 80 100 120 140 0 0.5 1 1.5 2 2.5 3 3.5 4 Separation Factor X (μm) S d (X)= S ∞ - Aexp[-X/B] S ∞ = 93.0000 A = 91.8586 B = 0.5290 S l (X)= M + NX M = 1.2072 N = 165.9645 H2/CO2 ΔP=40 atm S data S d S l S ∞ 0 20 40 60 80 100 120 140 0 2 4 6 8 10 Separation Factor X (μm) S d (X)= S ∞ - Aexp[-X/B] S ∞ = 104.0000 A = 101.4411 B = 2.1130 S l (X)= M + NX M = 2.5650 N = 47.3986 H2/CH4 ΔP=40 atm S data S d S l S ∞ Figure 6.22: The tting functions, S d (X) and S l (X), approaching the maximum separation factor for the H 2 /CO 2 (top) and H 2 /CH 4 (bottom) mixtures atT = 475 K and P = 40 atm. 155 Chapter 7 Conclusions The research in this dissertation dealt with various areas of science, such as quan- tum chemistry (QC), molecular dynamics (MD) simulation, and transport phe- nomena. In order to make it more convenient for the reader to follow the work, an overview of the important concepts and methods, which were used in this disser- tation, was provided in Chapter 2. In Chapter 3, the reactive force eld ReaxFF was developed in order to carry out MD simulation of pyrolysis of hydridopolycarbosilane (HPCS), which is a cru- cial step in the fabrication of nanoporous SiC membranes. The parameters of ReaxFF that are specically associated with the HPCS were estimated using QC computations based on density-functional theory. An atomistic model of the HPCS was then generated and utilized with ReaxFF in the MD simulation of its thermal decomposition. According to the MD simulations the degradation of the HPCS was initiated at about 1000 K. Cleavage of the Si{H bonds initiated the decom- position. At higher temperatures a large number of H radicals were formed, along with a secondary reaction that involved bond formation between the radicals that produced hydrogen gas, the main pyrolysis product. Reactive silylene intermedi- ates were also produced, which led to the formation of the Si{Si bonds that quickly broke above 1000 K, as Si tends to form more stable bonds with C. Cleavage of the Si{C bonds was initiated at about 2000 K and accelerated at 2400 K and higher, which led to the formation of some intermediate gaseous reaction products that eventually form methane, silane, and methylsilane gases, as well as some C 2 156 hydrocarbons. The results were consistent, both quantitatively and qualitatively, with the experimental data [23]. In Chapter 4 we 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 generated ceramic SiC that also contained some Si{Si and C{C bonds. Thus, the parameters of the ReaxFF that we have determined provide an accurate force eld for simulating thermal decomposition of an important class of Si-based materials, which have many applications in science and technology. In Chapter 5 the morphology of the molecular model of amorphous SiC that was developed in Chapter 4 was studied by computing its free-volume and cavity-size distributions. It was shown that the SiC model is a porous material with a range of pore sizes up to several A, suitable for diusion of small gaseous molecules. In addition, atomistic simulations were used to generate models of the SiC material containing various gases, by randomly inserting the gases in the SiC material. The self-diusivities of H 2 , CH 4 and CO 2 at 300 K were estimated using MD simulations, and the eect of temperature on the diusion of H 2 was studied. The simulations indicate that, among the three gases that we considered H 2 has the highest diusivity, which is important from the practical point of view of using the SiC material as a membrane layer for separation of H 2 from other gases. The results conrm the capability of the molecular model of SiC that we have developed in this dissertation. It was also shown that with increasing the temperature the diusivity of H 2 increases, which is expected. 157 The generated model was then utilized in Chapter 6 using non-equilibrium MD simulations to study transport and separation of two binary gaseous mixtures, namely H 2 /CO 2 and H 2 /CH 4 , under a variety of conditions. The model provided predictions for the separation properties of the SiC membrane that were in rea- sonable agreement with the experimental data [6]. In all the cases, the separation factor increased with increasing the thickness of the membrane. It was shown that for the two gas mixtures the separation factor reduces with increasing the external pressure drop. It is due to the formation of a liquid-like mixture at high external pressure drops that lower uxes of H 2 are obtained. The only dierence between the mixtures was that CO 2 has some anity for adsorption on the pores' surface, due to existence of the C atoms, hence the surface ow exists in the membrane for the H 2 /CO 2 mixture. Increasing the temperature caused the molecules to diuse much faster and, consequently, resulted in higher uxes and separation factors. Although simulating a system with the size used in practice is computationally prohibitive at this point, the eect of the membrane's thickness on the separation factor was studied by extrapolating the results for smaller thicknesses and tting them to a mathematical function. The optimal thickness of the membrane for each gaseous mixture was computed from the tted functions, which is an important result with practical implication. In summary, having obtained an accurate model for SiC nanoporous mem- branes, we believe that the proposed methodology using a reactive force led for the pyrolysis simulation of the polymer precursor is quite accurate, and applicable to the modeling of not only SiC membranes, but also a wide variety of inorganic membranes. The models can be used to study the transport and separation of mixtures in such membranes, and provide insight into improving the performance of the membrane as well their fabrication. 158 References [1] K. Nomura, R. K. Kalia, A. Nakano, and P. Vashishta, \A scalable paral- lel algorithm for large-scale reactive force-eld molecular dynamics simula- tions," Computer Physics Communications, vol. 178, no. 2, pp. 73{87, 2008. [2] P. Musumeci, L. Calcagno, and A. Makhtari, \Relaxation phenomena in keV- ion implanted hydrogenated amorphous silicon carbide," Materials Science and Engineering: A, vol. 253, no. 1, pp. 296{300, 1998. [3] J. A. Pople, \Nobel lecture: Quantum chemical models," Reviews of Modern Physics, vol. 71, no. 5, pp. 1267{1274, 1999. [4] M. Firouzi, T. T. Tsotsis, and M. Sahimi, \Nonequilibrium molecular dynam- ics simulations of transport and separation of supercritical uid mixtures in nanoporous membranes. I. Results for a single carbon nanopore," The Jour- nal of Chemical Physics, vol. 119, no. 13, pp. 6810{6822, 2003. [5] B. Elyassi, M. Sahimi, and T. T. Tsotsis, \Silicon carbide membranes for gas separation applications," Journal of Membrane Science, vol. 288, no. 1, pp. 290{297, 2007. [6] B. Elyassi, M. Sahimi, and T. T. Tsotsis, \A novel sacricial interlayer- based method for the preparation of silicon carbide membranes," Journal of Membrane Science, vol. 316, no. 1, pp. 73{79, 2008. [7] B. Elyassi, T. W. Kim, M. Sahimi, and T. T. Tsotsis, \Eect of polystyrene on the morphology and physical properties of silicon carbide nanobers," Journal of Materials Chemistry and Physics, vol. 118, no. 1, pp. 259{263, 2009. [8] B. Elyassi, M. Sahimi, and T. T. Tsotsis, \Inorganic membranes," in Ency- clopedia of Chemical Processing (S. Lee, ed.), pp. 1{16, Taylor Francis Group, LLC, 2009. 159 [9] M. Kotani, K. Nishiyabu, S. Matsuzaki, and S. Tanaka, \Processing of polymer-derived porous SiC body using allylhydridopolycarbosilane (AHPCS) and PMMA microbeads," Journal of the Ceramic Society of Japan, vol. 119, no. 1391, pp. 563{569, 2011. [10] Y. Takeda, K. Nakamura, K. Maeda, and Y. Matsushita, \Eects of elemen- tal additives on electrical resistivity of silicon carbide ceramics," Journal of the American Ceramic Society, vol. 70, no. 10, pp. C266{C267, 1987. [11] K. Schulz and M. Durst, \Advantages of an integrated system for hot gas ltration using rigid ceramic elements," Filtration & Separation, vol. 31, no. 1, pp. 25{28, 1994. [12] A. J. Rosenbloom, D. M. Sipe, Y. Shishkin, Y. Ke, R. P. Devaty, and W. J. Choyke, \Nanoporous SiC: A candidate semi-permeable material for biomed- ical applications," Biomedical Microdevices, vol. 6, no. 4, pp. 261{267, 2004. [13] Z. Li, K. Kusakabe, and S. Morooka, \Preparation of thermostable amor- phous Si{C{O membrane and its application to gas separation at elevated temperature," Journal of Membrane Science, vol. 118, no. 2, pp. 159{168, 1996. [14] C. A. Zorman, A. J. Fleischman, A. S. Dewa, M. Mehregany, C. Jacob, S. Nishino, and P. Pirouz, \Epitaxial growth of 3C{SiC lms on 4 in. diam (100) silicon wafers by atmospheric pressure chemical vapor deposition," Journal of Applied Physics, vol. 78, no. 8, pp. 5136{5138, 1995. [15] S. H. Kenawy and W. M. N. Nour, \Microstructural evaluation of thermally fatigued SiC{reinforced Al 2 O 3 /ZrO 2 matrix composites," Journal of Mate- rials Science, vol. 40, no. 14, pp. 3789{3793, 2005. [16] F. Chen, R. Mourhatch, T. T. Tsotsis, and M. Sahimi, \Experimental stud- ies and computer simulation of the preparation of nanoporous silicon-carbide membranes by chemical-vapor inltration/chemical-vapor deposition tech- niques," Chemical Engineering Science, vol. 63, no. 6, pp. 1460{1470, 2008. [17] R. Mourhatch, T. T. Tsotsis, and M. Sahimi, \Network model for the evo- lution of the pore structure of silicon-carbide membranes during their fabri- cation," Journal of Membrane Science, vol. 356, no. 1, pp. 138{146, 2010. [18] Y. Hasegawa and K. Okamura, \Synthesis of continuous silicon carbide bre," Journal of Materials Science, vol. 18, no. 12, pp. 3633{3648, 1983. [19] G. D. Soraru, F. Babonneau, and J. D. Mackenzie, \Structural concepts on new amorphous covalent solids," Journal of Non-Crystalline Solids, vol. 106, no. 1, pp. 256{261, 1988. 160 [20] G. D. Soraru, F. Babonneau, and J. D. Mackenzie, \Structural evolutions from polycarbosilane to SiC ceramic," Journal of Materials Science, vol. 25, no. 9, pp. 3886{3893, 1990. [21] X. Tang, L. Zhang, H. Tu, H. Gu, and L. Chen, \Decarbonization mech- anisms of polycarbosilane during pyrolysis in hydrogen for preparation of silicon carbide bers," Journal of Materials Science, vol. 45, no. 21, pp. 5749{ 5755, 2010. [22] L. V. Interrante, K. Moraes, L. MacDonald, and W. Sherwood, \Mechanical, thermochemical, and microstructural characterization of AHPCS-derived SiC," Advanced SiC/SiC Ceramic Composites: Developments and Applica- tions in Energy Systems, vol. 144, pp. 123{140, 2002. [23] L. V. Interrante, C. W. Whitmarsh, W. Sherwood, H. J. Wu, R. Lewis, and G. Maciel, \High yield polycarbosilane precursors to stoichiometric SiC. Syn- thesis, pyrolysis and application," in Materials Research Society Proceedings, vol. 346, pp. 593{603, Cambridge Univ Press, 1994. [24] Q. Liu, H. J. Wu, R. Lewis, G. E. Maciel, and L. V. Interrante, \Investigation of the pyrolytic conversion of poly(silylenemethylene) to silicon carbide," Chemistry of Materials, vol. 11, no. 8, pp. 2038{2048, 1999. [25] K. Farah, F. M uller-Plathe, and M. C. B ohm, \Classical reactive molecular dynamics implementations: State of the art," ChemPhysChem, vol. 13, no. 5, pp. 1127{1151, 2012. [26] A. C. T. van Duin, S. Dasgupta, F. Lorant, and W. A. Goddard, III, \ReaxFF: A reactive force eld for hydrocarbons," The Journal of Physi- cal Chemistry A, vol. 105, no. 41, pp. 9396{9409, 2001. [27] D. W. Brenner, \Empirical potential for hydrocarbons for use in simulating the chemical vapor deposition of diamond lms," Physical Review B, vol. 42, no. 15, pp. 9458{9471, 1990. [28] M. R. Nyden and D. W. Noid, \Molecular dynamics of initial events in the thermal degradation of polymers," The Journal of Physical Chemistry, vol. 95, no. 2, pp. 940{945, 1991. [29] K. Chenoweth, A. C. T. van Duin, and W. A. Goddard, III, \ReaxFF reactive force eld for molecular dynamics simulations of hydrocarbon oxidation," The Journal of Physical Chemistry A, vol. 112, no. 5, pp. 1040{1053, 2008. [30] K. Chenoweth, A. C. T. van Duin, P. Persson, M.-J. Cheng, J. Oxgaard, and W. A. Goddard, III, \Development and application of a ReaxFF reactive 161 force eld for oxidative dehydrogenation on vanadium oxide catalysts," The Journal of Physical Chemistry C, vol. 112, no. 37, pp. 14645{14654, 2008. [31] S. Cheung, W. Q. Deng, A. C. T. van Duin, and W. A. Goddard, III, \ReaxFF MgH reactive force eld for magnesium hydride systems," The Jour- nal of Physical Chemistry A, vol. 109, no. 5, pp. 851{859, 2005. [32] A. C. T. van Duin, A. Strachan, S. Stewman, Q. Zhang, X. Xu, and W. A. Goddard, III, \ReaxFF SiO reactive force eld for silicon and silicon oxide systems," The Journal of Physical Chemistry A, vol. 107, no. 19, pp. 3803{ 3811, 2003. [33] Q. Zhang, T. C a gn, A. C. T. van Duin, W. A. Goddard, III, Y. Qi, and L. G. Hector Jr, \Adhesion and nonwetting-wetting transition in the Al/{Al 2 O 3 interface," Physical Review B, vol. 69, no. 4, p. 045423, 2004. [34] A. Strachan, A. C. T. van Duin, D. Chakraborty, S. Dasgupta, and W. A. Goddard, III, \Shock waves in high-energy materials: The initial chemical events in nitramine RDX," Physical Review Letters, vol. 91, no. 9, p. 098301, 2003. [35] A. C. T. van Duin, Y. Zeiri, F. Dubnikova, R. Koslo, and W. A. Goddard, III, \Atomistic{scale simulations of the initial chemical events in the ther- mal initiation of triacetonetriperoxide," Journal of the American Chemical Society, vol. 127, no. 31, pp. 11053{11062, 2005. [36] K. D. Nielson, A. C. T. van Duin, J. Oxgaard, W. Q. Deng, and W. A. Goddard, III, \Development of the ReaxFF reactive force eld for describing transition metal catalyzed reactions, with application to the initial stages of the catalytic formation of carbon nanotubes," The Journal of Physical Chemistry A, vol. 109, no. 3, pp. 493{499, 2005. [37] L. Liu, A. Jaramillo-Botero, W. A. Goddard, III, and H. Sun, \Development of a ReaxFF reactive force eld for ettringite and study of its mechanical failure modes from reactive dynamics simulations," The Journal of Physical Chemistry A, vol. 116, no. 15, pp. 3918{3925, 2012. [38] K. Chenoweth, S. Cheung, A. C. T. Van Duin, W. A. Goddard, III, and E. M. Kober, \Simulations on the thermal decomposition of a poly(dimethylsiloxane) polymer using the ReaxFF reactive force eld," Jour- nal of the American Chemical Society, vol. 127, no. 19, pp. 7192{7202, 2005. [39] T. G. Desai, J. W. Lawson, and P. Keblinski, \Modeling initial stage of phenolic pyrolysis: Graphitic precursor formation and interfacial eects," Polymer, vol. 52, no. 2, pp. 577{585, 2011. 162 [40] B. Saha and G. C. Schatz, \Carbonization in polyacrylonitrile (PAN) based carbon bers studied by ReaxFF molecular dynamics simulations," The Journal of Physical Chemistry B, vol. 116, no. 15, pp. 4684{4692, 2012. [41] P. W. Atkins and R. S. Friedman, Molecular quantum mechanics, vol. 3. Oxford University Press, 1997. [42] A. Szabo and N. S. Ostlund, Modern quantum chemistry: Introduction to advanced electronic structure theory. Courier Dover Publications, 1989. [43] R. Shankar, Principles of quantum mechanics, vol. 2. Plenum Press, 1994. [44] C. C. J. Roothaan, \New developments in molecular orbital theory," Reviews of Modern Physics, vol. 23, no. 2, pp. 69{89, 1951. [45] J. C. Slater, \Atomic shielding constants," Physical Review, vol. 36, no. 1, pp. 57{64, 1930. [46] S. F. Boys, \Electronic wave functions. I. A general method of calculation for the stationary states of any molecular system," Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, vol. 200, no. 1063, pp. 542{554, 1950. [47] W. J. Hehre, R. F. Stewart, and J. A. Pople, \Self-consistent molecular- orbital methods. I. Use of gaussian expansions of slater-type atomic orbitals," The Journal of Chemical Physics, vol. 51, no. 6, pp. 2657{2664, 1969. [48] P. C. Hariharan and J. A. Pople, \The in uence of polarization functions on molecular orbital hydrogenation energies," Theoretica Chimica Acta, vol. 28, no. 3, pp. 213{222, 1973. [49] M. M. Francl, W. J. Pietro, W. J. Hehre, J. S. Binkley, M. S. Gordon, D. J. DeFrees, and J. A. Pople, \Self-consistent molecular orbital methods. XXIII. A polarization-type basis set for second-row elements," The Journal of Chemical Physics, vol. 77, no. 7, pp. 3654{3665, 1982. [50] W. J. Hehre, L. Radom, P. V. R. Schleyer, and J. A. Pople, Ab initio molec- ular orbital theory, vol. 33. Wiley-Interscience, 1986. [51] C. Mller and M. S. Plesset, \Note on an approximation treatment for many- electron systems," Physical Review, vol. 46, no. 7, pp. 618{622, 1934. [52] J. C zek, \On the correlation problem in atomic and molecular systems. Calculation of wavefunction components in ursell-type expansion using quantum-eld theoretical methods," The Journal of Chemical Physics, vol. 45, no. 11, pp. 4256{4266, 1966. 163 [53] P. R. Taylor, G. B. Bacskay, N. S. Hush, and A. C. Hurley, \The coupled- pair approximation in a basis of independent-pair natural orbitals," Chemical Physics Letters, vol. 41, no. 3, pp. 444{449, 1976. [54] R. J. Bartlett and G. D. Purvis, \Many-body perturbation theory, coupled- pair many-electron theory, and the importance of quadruple excitations for the correlation problem," International Journal of Quantum Chemistry, vol. 14, no. 5, pp. 561{581, 1978. [55] W. Kohn and L. J. Sham, \Self-consistent equations including exchange and correlation eects," Physical Review, vol. 140, no. 4A, pp. A1133{A1138, 1965. [56] J. P. Perdew, \Density-functional approximation for the correlation energy of the inhomogeneous electron gas," Physical Review B, vol. 33, no. 12, pp. 8822{8824, 1986. [57] A. D. Becke, \Density-functional thermochemistry. III. The role of exact exchange," The Journal of Chemical Physics, vol. 98, no. 7, pp. 5648{5652, 1993. [58] A. D. Becke, \Density-functional exchange-energy approximation with cor- rect asymptotic behavior," Physical Review A, vol. 38, no. 6, pp. 3098{3100, 1988. [59] C. Lee, W. Yang, and R. G. Parr, \Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density," Physical Review B, vol. 37, no. 2, pp. 785{789, 1988. [60] J. P. Perdew, K. Burke, and M. Ernzerhof, \Generalized gradient approx- imation made simple [Phys. Rev. Lett. 77, 3865 (1996)]," Physical Review Letter, vol. 78, no. 7, pp. 1396{1396, 1997. [61] K. Burke, \Perspective on density functional theory," The Journal of Chem- ical Physics, vol. 136, no. 15, p. 150901, 2012. [62] M. P. Allen and D. J. Tildesley, Computer simulation of liquids. Oxford University Press, 1989. [63] D. Frenkel and B. Smit, Understanding molecular simulation: From algo- rithms to applications. Academic Press, 2001. [64] J. M. Haile, Molecular dynamics simulation: Elementary methods. John Wiley & Sons, Inc., 1992. 164 [65] W. R. Gilks, S. Richardson, and D. J. Spiegelhalter, Markov chain Monte Carlo in practice, vol. 2. Chapman & Hall/CRC, 1996. [66] K. Binder and A. Baumg artner, The Monte Carlo method in condensed mat- ter physics. Springer, 1995. [67] R. Y. Rubinstein, Simulation and the Monte Carlo method, vol. 190. Wiley- Interscience, 2009. [68] M. H. Vedadi, Shock-induced nanobubble collapse and its applications. PhD thesis, University of Southern California, 2013. [69] W. J. Mortier, S. K. Ghosh, and S. Shankar, \Electronegativity{equalization method for the calculation of atomic charges in molecules," Journal of the American Chemical Society, vol. 108, no. 15, pp. 4315{4320, 1986. [70] G. O. A. Janssens, B. G. Baekelandt, H. Toufar, W. J. Mortier, and R. A. Schoonheydt, \Comparison of cluster and innite crystal calculations on zeo- lites with the electronegativity equalization method (EEM)," The Journal of Physical Chemistry, vol. 99, no. 10, pp. 3251{3258, 1995. [71] R. T. Sanderson, Chemical bonds and bond energy. Academic Press, 1971. [72] R. G. Parr and R. G. Pearson, \Absolute hardness: Companion parameter to absolute electronegativity," Journal of the American Chemical Society, vol. 105, no. 26, pp. 7512{7516, 1983. [73] R. T. Sanderson, Polar covalence. Academic Press, 1983. [74] W. J. Mortier, K. A. van Genechten, and J. Gasteiger, \Electronegativity equalization: Application and parametrization," Journal of the American Chemical Society, vol. 107, no. 4, pp. 829{835, 1985. [75] K. A. van Genechten, W. J. Mortier, and P. Geerlings, \Intrinsic framework electronegativity: A novel concept in solid state chemistry," The Journal of Chemical Physics, vol. 86, no. 9, pp. 5063{5071, 1987. [76] G. S. Heelnger and F. van Swol, \Diusion in Lennard-Jones uids using dual control volume grand canonical molecular dynamics simulation (DCV- GCMD)," The Journal of Chemical Physics, vol. 100, no. 10, pp. 7548{7552, 1994. [77] J. M. D. MacElroy, \Nonequilibrium molecular dynamics simulation of dif- fusion and ow in thin microporous membranes," The Journal of Chemical Physics, vol. 101, no. 6, pp. 5274{5280, 1994. 165 [78] L. Xu, M. Sahimi, and T. T. Tsotsis, \Nonequilibrium molecular dynam- ics simulations of transport and separation of gas mixtures in nanoporous materials," Physical Review E, vol. 62, no. 5, pp. 6942{6948, 2000. [79] N. Rajabbeigi, T. T. Tsotsis, and M. Sahimi, \Molecular pore-network model for nanoporous materials. II: Application to transport and separation of gaseous mixtures in silicon-carbide membranes," Journal of Membrane Sci- ence, vol. 345, no. 1, pp. 323{330, 2009. [80] R. C. Reid, J. M. Prausnitz, and B. E. Poling, The properties of gases and liquids. McGraw Hill, 1987. [81] Y. Kobayashi, S. Takami, M. Kubo, and A. Miyamoto, \Non-equilibrium molecular simulation studies on gas separation by microporous membranes using dual ensemble molecular simulation techniques," Fluid Phase Equilib- ria, vol. 194, pp. 319{326, 2002. [82] H. A. Lorentz, \ Uber die anwendung des satzes vom virial in der kinetischen theorie der gase," Annalen der Physik, vol. 248, no. 1, pp. 127{136, 1881. [83] A. des Sciences (Paris), Comptes rendus hebdomadaires des s eances de l'Acad emie des sciences, vol. 33. Gauthier-Villars, 1851. [84] L. Verlet, \Computer \experiments" on classical uids. I. Thermodynamical properties of Lennard-Jones molecules," Physical Review, vol. 159, no. 1, pp. 98{103, 1967. [85] G. Soave, \Equilibrium constants from a modied Redlich-Kwong equation of state," Chemical Engineering Science, vol. 27, no. 6, pp. 1197{1203, 1972. [86] A. Spitzbart, \A generalization of Hermite's interpolation formula," The American Mathematical Monthly, vol. 67, no. 1, pp. 42{46, 1960. [87] S. Naserifar, L. Liu, W. A. Goddard, III, T. T. Tsotsis, and M. Sahimi, \Toward a process-based molecular model of SiC membranes. I. Development of a reactive force eld," The Journal of Physical Chemistry C, vol. 117, no. 7, pp. 3308{3319, 2013. [88] S. M. Bachrach, \Jaguar 5.5," Journal of the American Chemical Society, vol. 126, no. 15, pp. 5018{5019, 2004. [89] A. D. Becke, \A new mixing of Hartree-Fock and local density-functional theories," The Journal of Chemical Physics, vol. 98, no. 2, pp. 1372{1377, 1993. 166 [90] R. Krishnan, J. S. Binkley, R. Seeger, and J. A. Pople, \Self-consistent molec- ular orbital methods. XX. A basis set for correlated wave functions," The Journal of Chemical Physics, vol. 72, no. 1, pp. 650{654, 1980. [91] T. Clark, J. Chandrasekhar, G. W. Spitznagel, and P. v. R. Schleyer, \E- cient diuse function-augmented basis sets for anion calculations. III. The 3-21+G basis set for rst-row elements, Li{F," Journal of Computational Chemistry, vol. 4, no. 3, pp. 294{301, 1983. [92] P. Giannozzi et al., \QUANTUM ESPRESSO: A modular and open-source software project for quantum simulations of materials," Journal of Physics: Condensed Matter, vol. 21, no. 39, p. 395502, 2009. [93] J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Peder- son, D. J. Singh, and C. Fiolhais, \Atoms, molecules, solids, and surfaces: Applications of the generalized gradient approximation for exchange and correlation," Physical Review B, vol. 46, no. 11, pp. 6671{6687, 1992. [94] H. J. Monkhorst and J. D. Pack, \Special points for Brillouin-zone integra- tions," Physical Review B, vol. 13, no. 12, pp. 5188{5192, 1976. [95] J. J. P. Stewart, \Optimization of parameters for semiempirical methods. I Method," Journal of Computational Chemistry, vol. 10, no. 2, pp. 209{220, 1989. [96] A. Taylor and R. M. Jones, \Silicon Carbide-a high temperature semicon- ductor," in Proceedings of the Conference on Silicon Carbide (Boston, MA 1959), pp. 147{154, 1960. [97] Cerius 2 , version 4.0; Accelyrs: San Diego, California, 1999. [98] M. E. Levinshtein, S. L. Rumyantsev, and M. S. Shur, Properties of advanced semiconductor materials: GaN, AIN, InN, BN, SiC, SiGe. Wiley- Interscience, 2001. [99] P. Patnaik, Handbook of inorganic chemicals, vol. 28. McGraw-Hill, 2003. [100] S. Plimpton et al., \Fast parallel algorithms for short-range molecular dynamics," Journal of Computational Physics, vol. 117, no. 1, pp. 1{19, 1995. [101] D. N. Theodorou and U. W. Suter, \Detailed molecular structure of a vinyl polymer glass," Macromolecules, vol. 18, no. 7, pp. 1467{1478, 1985. 167 [102] E. Tocci, D. Hofmann, D. Paul, N. Russo, and E. Drioli, \A molecular simula- tion study on gas diusion in a dense poly(ether{ether{ketone) membrane," Polymer, vol. 42, no. 2, pp. 521{533, 2001. [103] M. M. Ostwal, M. Sahimi, and T. T. Tsotsis, \Water harvesting using a conducting polymer: A study by molecular dynamics simulation," Physical Review E, vol. 79, no. 6, p. 061801, 2009. [104] A. R. Mehrabi and M. Sahimi, \Diusion of ionic particles in charged disor- dered media," Physical Review Letters, vol. 82, no. 4, pp. 735{738, 1999. [105] N. Kim, A. Harale, T. T. Tsotsis, and M. Sahimi, \Atomistic simulation of nanoporous layered double hydroxide materials and their properties. II. Adsorption and diusion," The Journal of Chemical Physics, vol. 127, no. 22, p. 224701, 2007. [106] A. R. Mehrabi and M. Sahimi, \Cluster conformations and multipole dis- tributions in ionic uids. I. Two-dimensional systems of mobile ions," The Journal of Chemical Physics, vol. 128, no. 23, p. 234503, 2008. [107] D. Wolf, P. Keblinski, S. R. Phillpot, and J. Eggebrecht, \Exact method for the simulation of Coulombic systems by spherically truncated, pairwise r 1 summation," The Journal of Chemical Physics, vol. 110, no. 17, pp. 8254{ 8282, 1999. [108] M. Sahimi, Heterogeneous materials II: Nonlinear and breakdown properties and atomistic modeling, vol. 23. Springer, 2003. [109] H. C. Andersen, \Molecular dynamics simulations at constant pressure and/or temperature," The Journal of Chemical Physics, vol. 72, no. 4, pp. 2384{2391, 1980. [110] S. Nos e, \A unied formulation of the constant temperature molecular dynamics methods," The Journal of Chemical Physics, vol. 81, no. 1, pp. 511{519, 1984. [111] W. G. Hoover, \Canonical dynamics: Equilibrium phase-space distribu- tions," Physical Review A, vol. 31, no. 3, pp. 1695{1697, 1985. [112] S. Naserifar, W. A. Goddard, III, L. Liu, T. T. Tsotsis, and M. Sahimi, \Toward a process-based molecular model of SiC membranes. II. Reactive dynamics simulation of the pyrolysis of polymer precursor to form amorphous SiC," The Journal of Physical Chemistry C, vol. 117, no. 7, pp. 3320{3329, 2013. 168 [113] Starre Systems, Inc., http://www.starressystems.com. [114] S. Kirkpatrick, D. G. Jr., and M. P. Vecchi, \Optimization by simulated annealing," Science, vol. 220, no. 4598, pp. 671{680, 1983. [115] D. R. Lide and T. J. Bruno, CRC handbook of chemistry and physics. CRC Press, 2012. [116] M. Ishimaru, I.-T. Bae, Y. Hirotsu, S. Matsumura, and K. E. Sickafus, \Structural relaxation of amorphous silicon carbide," Physical Review Let- ters, vol. 89, no. 5, p. 055502, 2002. [117] D. A. McQuarrie, Statistical mechanics. University Science Books, 2000. [118] M. Sahimi, Flow and transport in porous media and fractured rock. Wiley- VCH, 2012. [119] S. Torquato, Random heterogeneous materials: Microstructure and macro- scopic properties, vol. 16. Springer, 2001. [120] T. J. Pinnavaia and M. F. Thorpe, Access in nanoporous materials. Plenum Press, 1995. [121] H. Kral, J. Rouquerol, K. S. W. Sing, and K. K. Unger, Characterization of porous solids, vol. 39. Elsevier Science, 1988. [122] L. D. Gelb, K. E. Gubbins, R. Radhakrishnan, and M. Sliwinska-Bartkowiak, \Phase separation in conned systems," Reports on Progress in Physics, vol. 62, no. 12, pp. 1573{1659, 1999. [123] M. Sahimi and T. T. Tsotsis, \Computational methods for statistical meth- ods of nanoporous materials and their properties," in Handbook of Theoreti- cal and Computational Nanotechnology (W. Schommers and M. Reith, eds.), ch. 74, American Scientic Publications, 2005. [124] S. Y. Lim, M. Sahimi, T. T. Tsotsis, and N. Kim, \Molecular dynamics simulation of diusion of gases in a carbon-nanotube-polymer composite," Physical Review E, vol. 76, no. 1, p. 011810, 2007. [125] J. Ghassemzadeh, L. Xu, T. T. Tsotsis, and M. Sahimi, \Statistical mechan- ics and molecular simulation of adsorption in microporous materials: Pil- lared clays and carbon molecular sieve membranes," The Journal of Physical Chemistry B, vol. 104, no. 16, pp. 3892{3905, 2000. [126] J. Ghassemzadeh and M. Sahimi, \Molecular modelling of adsorption of gas mixtures in montmorillonites intercalated with Al13-complex pillars," Molec- ular Physics, vol. 102, no. 13, pp. 1447{1467, 2004. 169 [127] L. Xu, T. T. Tsotsis, and M. Sahimi, \Statistical mechanics and molecular simulation of adsorption of ternary gas mixtures in nanoporous materials," The Journal of Chemical Physics, vol. 114, no. 16, pp. 7196{7210, 2001. [128] S. Y. Lim, T. T. Tsotsis, and M. Sahimi, \Molecular simulation of diusion and sorption of gases in an amorphous polymer," The Journal of Chemical Physics, vol. 119, no. 1, pp. 496{504, 2003. [129] N. Rajabbeigi, B. Elyassi, T. T. Tsotsis, and M. Sahimi, \Molecular pore- network model for nanoporous materials. I: Application to adsorption in silicon-carbide membranes," Journal of Membrane Science, vol. 335, no. 1, pp. 5{12, 2009. [130] N. Kim, Y. Kim, T. T. Tsotsis, and M. Sahimi, \Atomistic simulation of nanoporous layered double hydroxide materials and their properties. I. structural modeling," The Journal of Chemical Physics, vol. 122, no. 21, p. 214713, 2005. [131] A. K. Rappe, C. J. Casewit, K. S. Colwell, W. A. Goddard, III, and W. M. Ski, \UFF, a full periodic table force eld for molecular mechanics and molecular dynamics simulations," Journal of the American Chemical Society, vol. 114, no. 25, pp. 10024{10035, 1992. [132] H.-Q. Ding, N. Karasawa, and W. A. Goddard, III, \Atomic level simula- tions on a million particles: The cell multipole method for Coulomb and London nonbond interactions," The Journal of Chemical Physics, vol. 97, no. 6, pp. 4309{4315, 1992. [133] M. Tanemura, T. Ogawa, and N. Ogita, \A new algorithm for three- dimensional Voronoi tessellation," Journal of Computational Physics, vol. 51, no. 2, pp. 191{207, 1983. [134] J. Budzien, J. D. McCoy, and D. B. Adolf, \Solute mobility and packing fraction: A new look at the Doolittle equation for the polymer glass tran- sition," The Journal of Chemical Physics, vol. 119, no. 17, pp. 9269{9273, 2003. [135] M. H. Cohen and D. Turnbull, \Molecular transport in liquids and glasses," The Journal of Chemical Physics, vol. 31, no. 5, pp. 1164{1169, 1959. 170 Appendix A Bond and Angle Distribution for HPCS Figures A.1 to A.3 present the distributions of possible bonds, triplet angles and dihedral angles for the HPCS atomistic model, respectively (see section 3.4). 171 0 1 2 3 4 5 1.6 1.7 1.8 1.9 2 2.1 2.2 2.3 2.4 Probability (1/Å) Si−C Bond Length (Å) (a) 0 1 2 3 4 5 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 Probability (1/Å) Si−H Bond Length (Å) (b) 0 1 2 3 4 5 6 7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 Probability (1/Å) C−H Bond Length (Å) (c) Figure A.1: Distribution of (a) Si{C ; (b) Si{H, and (c) C{H distances in the HPCS polymer, as generated by the PCFF. 172 0 0.01 0.02 0.03 0.04 60 80 100 120 140 160 Probability (1/degrees) Si-C-Si Angle (degrees) (a) 0 0.01 0.02 0.03 0.04 60 80 100 120 140 160 Probability (1/degrees) Si-C-H Angle (degrees) (b) 0 0.01 0.02 0.03 0.04 60 80 100 120 140 160 Probability (1/degrees) C-Si-H Angle (degrees) (c) 0 0.01 0.02 0.03 0.04 60 80 100 120 140 160 Probability (1/degrees) C-Si-C Angle (degrees) (d) Figure A.2: Distribution of (a) Si{C{Si ; (b) Si{C{H, (c) C{Si{H and (d) C{Si{C triplet angles in the HPCS polymer, as generated by the PCFF. 173 0 0.001 0.002 0.003 0.004 0.005 -200 -150 -100 -50 0 50 100 150 200 Probability (1/degrees) H-C-Si-C Angle (degrees) (a) 0 0.001 0.002 0.003 0.004 0.005 -200 -150 -100 -50 0 50 100 150 200 Probability (1/degrees) H-Si-C-Si Angle (degrees) (b) 0 0.001 0.002 0.003 0.004 0.005 -200 -150 -100 -50 0 50 100 150 200 Probability (1/degrees) H-C-Si-H Angle (degrees) (c) 0 0.001 0.002 0.003 0.004 0.005 0.006 -200 -150 -100 -50 0 50 100 150 200 Probability (1/degrees) Si-C-Si-C Angle (degrees) (d) Figure A.3: Distribution of (a) H{C{Si{C ; (b) H{Si{C{Si, (c) H{C{Si{H and (d) Si{C{Si{C dihedral angles in the HPCS polymer, as generated by the PCFF. 174 Appendix B Temperature-Dependence Production of H 2 B.1 Smoothing the Population Data Because of high temperature during (NVT )-MD simulations, as well as the high energy state of the H radicals, the bond formed between two H radicals can break and form several times, which causes uctuations in the population of H 2 molecules. A lter known as exponential moving average (EMA), or an exponentially weighted moving average (EWMA), was used to smooth the H 2 population data. The l- ter has the type of innite impulse response with weighting factors that decrease exponentially. The EMA can be calculated for the H 2 population data recursively, and the weighting decreases exponentially for each previous data point, but never vanishes. It has the form of, t> 1;P t =D t + (1)P t1 ; (B.1) where, is the degree of weighting decrease (between 0 and 1), which was selected to be 0.9, D t is the population of H 2 molecules at time t, and P t is the value of EMA at timet. The value ofP 1 was initialized by the rst value in H 2 population data set (i.e. at t = 0). 175 B.2 Reaction Rate at Various Temperature The rate of H 2 production increase during (NVT )-MD simulations strongly depends on the temperature and can be tted to an exponential function, f(t) =N f (N f N 0 ) exp(kt); (B.2) where, N f and N 0 are initial and asymptotic number of H 2 , respectively, and k is reaction rate (i.e. an inverse of average characteristic time). The exponential function,f(t), was tted against original and smoothed data for H 2 production at various temperatures (3000-4000 K) and for three trials (dierent initial congura- tions). The results are presented in Figures B.1 to B.18. The asymptotic standard errors for all the ttings were found to be less than 1 percent. The computed reaction rates from smoothed data were tted using the classical Arrhenius form for the temperature-dependence of a rst-order reaction rate constant, k(T ) =k 0 exp E a RT ; (B.3) wherek is reaction rate (s 1 ),E a is the activation energy (kcal/mol), andA is the frequency prefactor (molL 1 s 1 ). 176 0 2 4 6 8 10 12 14 16 18 0 5 10 15 20 25 30 35 40 45 50 55 60 Number of Hydrogen Molecules Time (ps) 3000K, Trial 4 3000K 3000K smoothed 0 2 4 6 8 10 12 14 16 0 5 10 15 20 25 30 35 40 45 50 55 60 Number of Hydrogen Molecules Time (ps) 3000K, Trial 4 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.0262444 ps -1 3000K f(t) 0 3 6 9 12 15 0 5 10 15 20 25 30 35 40 45 50 55 60 Number of Hydrogen Molecules Time (ps) 3000K, Trial 4 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.0173742 ps -1 3000K smoothed f(t) Figure B.1: Temperature-dependence of the rate of production of H 2 in both orig- inal and smoothed data (top) at 3000 K for trial 1, tted f(t) against original (middle) and smoothed data (bottom). 177 0 3 6 9 12 15 0 5 10 15 20 25 30 35 40 45 50 55 60 Number of Hydrogen Molecules Time (ps) 3000K, Trial 5 3000K 3000K smoothed 0 3 6 9 12 15 0 5 10 15 20 25 30 35 40 45 50 55 60 Number of Hydrogen Molecules Time (ps) 3000K, Trial 5 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.0174332 ps -1 3000K f(t) 0 3 6 9 12 15 0 5 10 15 20 25 30 35 40 45 50 55 60 Number of Hydrogen Molecules Time (ps) 3000K, Trial 5 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.0173742 ps -1 3000K smoothed f(t) Figure B.2: Temperature-dependence of the rate of production of H 2 in both orig- inal and smoothed data (top) at 3000 K for trial 2, tted f(t) against original (middle) and smoothed data (bottom). 178 0 2 4 6 8 10 12 14 16 18 20 0 5 10 15 20 25 30 35 40 45 50 55 60 Number of Hydrogen Molecules Time (ps) 3000K, Trial 6 3000K 3000K smoothed 0 2 4 6 8 10 12 14 16 18 20 0 5 10 15 20 25 30 35 40 45 50 55 60 Number of Hydrogen Molecules Time (ps) 3000K, Trial 6 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.0274611 ps -1 3000K f(t) 0 2 4 6 8 10 12 14 16 18 20 0 5 10 15 20 25 30 35 40 45 50 55 60 Number of Hydrogen Molecules Time (ps) 3000K, Trial 6 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.0273507 ps -1 3000K smoothed f(t) Figure B.3: Temperature-dependence of the rate of production of H 2 in both orig- inal and smoothed data (top) at 3000 K for trial 3, tted f(t) against original (middle) and smoothed data (bottom). 179 0 4 8 12 16 20 24 28 32 0 5 10 15 20 25 30 35 40 45 50 55 60 Number of Hydrogen Molecules Time (ps) 3200K, Trial 4 3200K 3200K smoothed 0 4 8 12 16 20 24 28 32 0 5 10 15 20 25 30 35 40 45 50 55 60 Number of Hydrogen Molecules Time (ps) 3200K, Trial 4 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.0321686 ps -1 3200K f(t) 0 4 8 12 16 20 24 28 32 0 5 10 15 20 25 30 35 40 45 50 55 60 Number of Hydrogen Molecules Time (ps) 3200K, Trial 4 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.032039 ps -1 3200K smoothed f(t) Figure B.4: Temperature-dependence of the rate of production of H 2 in both orig- inal and smoothed data (top) at 3200 K for trial 1, tted f(t) against original (middle) and smoothed data (bottom). 180 0 4 8 12 16 20 0 5 10 15 20 25 30 35 40 45 50 55 60 Number of Hydrogen Molecules Time (ps) 3200K, Trial 5 3200K 3200K smoothed 0 4 8 12 16 20 0 5 10 15 20 25 30 35 40 45 50 55 60 Number of Hydrogen Molecules Time (ps) 3200K, Trial 5 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.0303268 ps -1 3200K f(t) 0 4 8 12 16 20 0 5 10 15 20 25 30 35 40 45 50 55 60 Number of Hydrogen Molecules Time (ps) 3200K, Trial 5 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.0302331 ps -1 3200K smoothed f(t) Figure B.5: Temperature-dependence of the rate of production of H 2 in both orig- inal and smoothed data (top) at 3200 K for trial 2, tted f(t) against original (middle) and smoothed data (bottom). 181 0 4 8 12 16 20 24 0 5 10 15 20 25 30 35 40 45 50 55 60 Number of Hydrogen Molecules Time (ps) 3200K, Trial 6 3200K 3200K smoothed 0 4 8 12 16 20 24 0 5 10 15 20 25 30 35 40 45 50 55 60 Number of Hydrogen Molecules Time (ps) 3200K, Trial 6 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.0323923 ps -1 3200K f(t) 0 4 8 12 16 20 24 0 5 10 15 20 25 30 35 40 45 50 55 60 Number of Hydrogen Molecules Time (ps) 3200K, Trial 6 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.0322618 ps -1 3200K smoothed f(t) Figure B.6: Temperature-dependence of the rate of production of H 2 in both orig- inal and smoothed data (top) at 3200 K for trial 3, tted f(t) against original (middle) and smoothed data (bottom). 182 0 4 8 12 16 20 24 28 32 0 5 10 15 20 25 30 35 40 45 50 55 60 Number of Hydrogen Molecules Time (ps) 3400K, Trial 4 3400K 3400K smoothed 0 4 8 12 16 20 24 28 32 0 5 10 15 20 25 30 35 40 45 50 55 60 Number of Hydrogen Molecules Time (ps) 3400K, Trial 4 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.041055 ps -1 3400K f(t) 0 4 8 12 16 20 24 28 32 0 5 10 15 20 25 30 35 40 45 50 55 60 Number of Hydrogen Molecules Time (ps) 3400K, Trial 4 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.0408728 ps -1 3400K smoothed f(t) Figure B.7: Temperature-dependence of the rate of production of H 2 in both orig- inal and smoothed data (top) at 3400 K for trial 1, tted f(t) against original (middle) and smoothed data (bottom). 183 0 5 10 15 20 25 30 35 40 45 0 5 10 15 20 25 30 35 40 45 50 55 60 Number of Hydrogen Molecules Time (ps) 3400K, Trial 5 3400K 3400K smoothed 0 5 10 15 20 25 30 35 40 45 0 5 10 15 20 25 30 35 40 45 50 55 60 Number of Hydrogen Molecules Time (ps) 3400K, Trial 5 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.041055 ps -1 3400K f(t) 0 5 10 15 20 25 30 35 40 45 0 5 10 15 20 25 30 35 40 45 50 55 60 Number of Hydrogen Molecules Time (ps) 3400K, Trial 5 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.0408728 ps -1 3400K smoothed f(t) Figure B.8: Temperature-dependence of the rate of production of H 2 in both orig- inal and smoothed data (top) at 3400 K for trial 2, tted f(t) against original (middle) and smoothed data (bottom). 184 0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 30 35 40 45 50 55 60 Number of Hydrogen Molecules Time (ps) 3400K, Trial 6 3400K 3400K smoothed 0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 30 35 40 45 50 55 60 Number of Hydrogen Molecules Time (ps) 3400K, Trial 6 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.0384296 ps -1 3400K f(t) 0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 30 35 40 45 50 55 60 Number of Hydrogen Molecules Time (ps) 3400K, Trial 6 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.0382548 ps -1 3400K smoothed f(t) Figure B.9: Temperature-dependence of the rate of production of H 2 in both orig- inal and smoothed data (top) at 3400 K for trial 3, tted f(t) against original (middle) and smoothed data (bottom). 185 0 5 10 15 20 25 30 35 0 2 4 6 8 10 12 14 16 18 20 Number of Hydrogen Molecules Time (ps) 3600K, Trial 4 3600K 3600K smoothed 0 5 10 15 20 25 30 35 0 2 4 6 8 10 12 14 16 18 20 Number of Hydrogen Molecules Time (ps) 3600K, Trial 4 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.0809067 ps -1 3600K f(t) 0 5 10 15 20 25 30 35 0 2 4 6 8 10 12 14 16 18 20 Number of Hydrogen Molecules Time (ps) 3600K, Trial 4 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.0799084 ps -1 3600K smoothed f(t) Figure B.10: Temperature-dependence of the rate of production of H 2 in both original and smoothed data (top) at 3600 K for trial 1, tted f(t) against original (middle) and smoothed data (bottom). 186 0 5 10 15 20 25 30 0 2 4 6 8 10 12 14 16 18 20 Number of Hydrogen Molecules Time (ps) 3600K, Trial 5 3600K 3600K smoothed 0 5 10 15 20 25 30 0 2 4 6 8 10 12 14 16 18 20 Number of Hydrogen Molecules Time (ps) 3600K, Trial 5 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.0835207 ps -1 3600K f(t) 0 5 10 15 20 25 30 0 2 4 6 8 10 12 14 16 18 20 Number of Hydrogen Molecules Time (ps) 3600K, Trial 5 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.0825157 ps -1 3600K smoothed f(t) Figure B.11: Temperature-dependence of the rate of production of H 2 in both original and smoothed data (top) at 3600 K for trial 2, tted f(t) against original (middle) and smoothed data (bottom). 187 0 5 10 15 20 25 30 35 0 2 4 6 8 10 12 14 16 18 20 Number of Hydrogen Molecules Time (ps) 3600K, Trial 6 3600K 3600K smoothed 0 5 10 15 20 25 30 35 0 2 4 6 8 10 12 14 16 18 20 Number of Hydrogen Molecules Time (ps) 3600K, Trial 6 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.0882821 ps -1 3600K f(t) 0 5 10 15 20 25 30 35 0 2 4 6 8 10 12 14 16 18 20 Number of Hydrogen Molecules Time (ps) 3600K, Trial 6 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.0872054 ps -1 3600K smoothed f(t) Figure B.12: Temperature-dependence of the rate of production of H 2 in both original and smoothed data (top) at 3600 K for trial 3, tted f(t) against original (middle) and smoothed data (bottom). 188 0 5 10 15 20 25 30 35 40 45 50 0 2 4 6 8 10 12 14 16 18 20 Number of Hydrogen Molecules Time (ps) 3800K, Trial 4 3800K 3800K smoothed 0 5 10 15 20 25 30 35 40 45 50 0 2 4 6 8 10 12 14 16 18 20 Number of Hydrogen Molecules Time (ps) 3800K, Trial 4 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.0970347 ps -1 3800K f(t) 0 5 10 15 20 25 30 35 40 45 50 0 2 4 6 8 10 12 14 16 18 20 Number of Hydrogen Molecules Time (ps) 3800K, Trial 4 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.118124 ps -1 3800K smoothed f(t) Figure B.13: Temperature-dependence of the rate of production of H 2 in both original and smoothed data (top) at 3800 K for trial 1, tted f(t) against original (middle) and smoothed data (bottom). 189 0 5 10 15 20 25 30 35 0 2 4 6 8 10 12 14 16 18 20 Number of Hydrogen Molecules Time (ps) 3800K, Trial 5 3800K 3800K smoothed 0 5 10 15 20 25 30 35 0 2 4 6 8 10 12 14 16 18 20 Number of Hydrogen Molecules Time (ps) 3800K, Trial 5 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.119635 ps -1 3800K f(t) 0 5 10 15 20 25 30 35 0 2 4 6 8 10 12 14 16 18 20 Number of Hydrogen Molecules Time (ps) 3800K, Trial 5 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.118124 ps -1 3800K smoothed f(t) Figure B.14: Temperature-dependence of the rate of production of H 2 in both original and smoothed data (top) at 3800 K for trial 2, tted f(t) against original (middle) and smoothed data (bottom). 190 0 5 10 15 20 25 30 35 40 45 50 0 2 4 6 8 10 12 14 16 18 20 Number of Hydrogen Molecules Time (ps) 3800K, Trial 6 3800K 3800K smoothed 0 5 10 15 20 25 30 35 40 45 50 0 2 4 6 8 10 12 14 16 18 20 Number of Hydrogen Molecules Time (ps) 3800K, Trial 6 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.0970347 ps -1 3800K f(t) 0 5 10 15 20 25 30 35 40 45 50 0 2 4 6 8 10 12 14 16 18 20 Number of Hydrogen Molecules Time (ps) 3800K, Trial 6 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.0959031 ps -1 3800K smoothed f(t) Figure B.15: Temperature-dependence of the rate of production of H 2 in both original and smoothed data (top) at 3800 K for trial 3, tted f(t) against original (middle) and smoothed data (bottom). 191 0 10 20 30 40 50 60 0 5 10 15 20 Number of Hydrogen Molecules Time (ps) 4000K, Trial 4 4000K 4000K smoothed 0 5 10 15 20 25 30 35 40 45 50 55 60 0 2 4 6 8 10 12 14 16 18 20 Number of Hydrogen Molecules Time (ps) 4000K, Trial 4 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.136001 ps -1 4000K f(t) 0 5 10 15 20 25 30 35 40 45 50 55 60 0 2 4 6 8 10 12 14 16 18 20 Number of Hydrogen Molecules Time (ps) 4000K, Trial 4 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.134419 ps -1 4000K smoothed f(t) Figure B.16: Temperature-dependence of the rate of production of H 2 in both original and smoothed data (top) at 4000 K for trial 1, tted f(t) against original (middle) and smoothed data (bottom). 192 0 10 20 30 40 50 60 0 5 10 15 20 Number of Hydrogen Molecules Time (ps) 4000K, Trial 5 4000K 4000K smoothed 0 5 10 15 20 25 30 35 40 45 50 55 60 0 2 4 6 8 10 12 14 16 18 20 Number of Hydrogen Molecules Time (ps) 4000K, Trial 5 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.126645 ps -1 4000K f(t) 0 5 10 15 20 25 30 35 40 45 50 55 60 0 2 4 6 8 10 12 14 16 18 20 Number of Hydrogen Molecules Time (ps) 4000K, Trial 5 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.125042 ps -1 4000K smoothed f(t) Figure B.17: Temperature-dependence of the rate of production of H 2 in both original and smoothed data (top) at 4000 K for trial 2, tted f(t) against original (middle) and smoothed data (bottom). 193 0 10 20 30 40 50 0 5 10 15 20 Number of Hydrogen Molecules Time (ps) 4000K, Trial 6 4000K 4000K smoothed 0 5 10 15 20 25 30 35 40 45 50 0 2 4 6 8 10 12 14 16 18 20 Number of Hydrogen Molecules Time (ps) 4000K, Trial 6 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.128283 ps -1 4000K f(t) 0 5 10 15 20 25 30 35 40 45 50 0 2 4 6 8 10 12 14 16 18 20 Number of Hydrogen Molecules Time (ps) 4000K, Trial 6 f(t) = N f - (N f - N 0 ) exp(-kt) k = 0.126858 ps -1 4000K smoothed f(t) Figure B.18: Temperature-dependence of the rate of production of H 2 in both original and smoothed data (top) at 4000 K for trial 3, tted f(t) against original (middle) and smoothed data (bottom). 194 Appendix C Separation Factor Fitting and Extrapolation Figures C.1 to C.4, similar to gures in Section 6.2.5, present the change in the separation factor with the thickness, as well as the t of the results for the H 2 /CO 2 and H 2 /CH 4 mixtures with P = 50 and 60 atm. 195 0 2 4 6 8 10 12 14 16 0 0.01 0.02 0.03 0.04 0.05 Separation Factor X (μm) S d (X)= S ∞ - Aexp[-X/B] S ∞ = 93.0000 A = 91.4562 B = 0.8790 H2/CO2 ΔP=50 atm S data S d 0 2 4 6 8 10 12 14 16 0 0.01 0.02 0.03 0.04 0.05 Separation Factor X (μm) S d (X)= S ∞ - Aexp[-X/B] S ∞ = 104.0000 A = 101.4193 B = 1.9095 H2/CH4 ΔP=50 atm S data S d Figure C.1: Eect of the membrane's thickness on its separation factor for the H 2 /CO 2 (top) and H 2 /CH 4 (bottom) mixtures at T = 475 K and P = 50 atm. The data are tted to the S d (X) function. 196 0 2 4 6 8 10 12 14 16 0 0.01 0.02 0.03 0.04 0.05 Separation Factor X (μm) S d (X)= S ∞ - Aexp[-X/B] S ∞ = 93.0000 A = 91.4562 B = 0.8790 H2/CO2 ΔP=60 atm S data S d 0 2 4 6 8 10 12 14 16 0 0.01 0.02 0.03 0.04 0.05 Separation Factor X (μm) S d (X)= S ∞ - Aexp[-X/B] S ∞ = 104.0000 A = 101.7198 B = 1.4486 H2/CH4 ΔP=60 atm S data S d Figure C.2: Eect of the membrane's thickness on its separation factor for the H 2 /CO 2 (top) and H 2 /CH 4 (bottom) mixtures at T = 475 K and P = 60 atm. The data are tted to the S d (X) function. 197 0 20 40 60 80 100 120 140 0 1 2 3 4 5 Separation Factor X (μm) S d (X)= S ∞ - Aexp[-X/B] S ∞ = 93.0000 A = 91.4562 B = 0.8790 S l (X)= M + NX M = 1.5717 N = 101.0587 H2/CO2 ΔP=50 atm S data S d S l S ∞ 0 20 40 60 80 100 120 140 0 2 4 6 8 10 Separation Factor X (μm) S d (X)= S ∞ - Aexp[-X/B] S ∞ = 104.0000 A = 101.4193 B = 1.9095 S l (X)= M + NX M = 2.5873 N = 52.4068 H2/CH4 ΔP=50 atm S data S d S l S ∞ Figure C.3: The tting functions, S d (X) and S l (X), approaching the maximum separation factor for the H 2 /CO 2 (top) and H 2 /CH 4 (bottom) mixtures atT = 475 K and P = 50 atm. 198 0 20 40 60 80 100 120 140 0 1 2 3 4 5 Separation Factor X (μm) S d (X)= S ∞ - Aexp[-X/B] S ∞ = 93.0000 A = 91.4562 B = 0.8790 S l (X)= M + NX M = 1.5717 N = 101.0587 H2/CO2 ΔP=60 atm S data S d S l S ∞ 0 20 40 60 80 100 120 140 0 1 2 3 4 5 6 7 8 Separation Factor X (μm) S d (X)= S ∞ - Aexp[-X/B] S ∞ = 104.0000 A = 101.7198 B = 1.4486 S l (X)= M + NX M = 2.2890 N = 69.1188 H2/CH4 ΔP=60 atm S data S d S l S ∞ Figure C.4: The tting functions, S d (X) and S l (X), approaching the maximum separation factor for the H 2 /CO 2 (top) and H 2 /CH 4 (bottom) mixtures atT = 475 K and P = 60 atm. 199
Abstract (if available)
Abstract
A broad class of important materials, such as silicon carbide (SiC), are fabricated by temperature-controlled pyrolysis of pre-ceramic polymers. In particular, fabrication of SiC membranes by pyrolysis of a polymer precursor that contains Si is quite attractive for separation of hydrogen from other gases. It has been quite difficult to extract atomistic-scale information about such SiC membranes, since they are amorphous. The research presented in this dissertation extends the ReaxFF reactive force field to describe the processes involved in the thermal decomposition of hydridopolycarbosilane (HPCS) to form SiC nanoporous membranes. ❧ First, we carry out quantum mechanical calculations on models meant to capture the important reaction steps and structures. Then, we develop a model of the HPCS polymer and utilize ReaxFF to describe the thermal degradation and decomposition of the polymer as the system is heated by molecular dynamics (MD) simulations. In the next step, we use ReaxFF for reactive dynamics of HPCS over a wide range of temperature. We then simulate pyrolysis of the HPCS under conditions that closely mimic the conditions of the fabrication of nanoporous SiC membranes to produce amorphous SiC material. The pyrolysis results and the computed properties of the SiC ceramic are shown to be in good agreement with the experimental data. ❧ To further test the validity of the model and to provide insight into improving the performance of the membrane, extensive MD simulations were carried out to compute the self-diffusivities of H₂, CO₂ and CH₄ in the molecular model of SiC. The results indicate higher values of the diffusivity for H₂. The morphology of the amorphous SiC layer is characterized by computing its accessible free volumes and the cavity distributions. ❧ Finally, we use the molecular model of SiC and non-equilibrium MD simulations in order to study transport and separation in the membrane of two binary gaseous mixtures, namely, H₂/CO₂ and H₂/CH₄ at various temperatures and pressure drops, applied externally to the membrane. When compared with our own experimental data, the model is demonstrated to provide accurate predictions for various properties of interest and, in particular, for the separation factors of the mixtures. The model can be used to determine the optimal membrane's thickness.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Silicon carbide ceramic membranes
PDF
Fabrication of silicon carbide sintered supports and silicon carbide membranes
PDF
Molecular modeling of silicon carbide nanoporous membranes and transport and adsorption of gaseous mixtures therein
PDF
Continuum and pore netwok modeling of preparation of silicon-carbide membranes by chemical-vapor deposition and chemical-vapor infiltration
PDF
Fabrication of nanoporous silicon carbide membranes for gas separation applications
PDF
Experimental studies and computer simulation of the preparation of nanoporous silicon-carbide membranes by chemical-vapor infiltration/chemical-vapor deposition techniques
PDF
Exploring properties of silicon-carbide nanotubes and their composites with polymers
PDF
Fabrication of silicon-based membranes via vapor-phase deposition and pyrolysis of organosilicon polymers
PDF
Fabrication of nanoporous silicon oxycarbide materials via a sacrificial template technique
PDF
Development of carbon molecular-sieve membranes with tunable properties: modification of the pore size and surface affinity
PDF
Methanol synthesis in a membrane reactor
PDF
A study of the application of membrane-based reactive separation to the carbon dioxide methanation
PDF
Process intensification in hydrogen production via membrane-based reactive separations
PDF
The use of carbon molecule sieve and Pd membranes for conventional and reactive applications
PDF
Molecular-scale studies of mechanical phenomena at the interface between two solid surfaces: from high performance friction to superlubricity and flash heating
PDF
On the use of membrane reactors in biomass utilization
PDF
Quantum molecular dynamics and machine learning for non-equilibrium processes in quantum materials
PDF
Modeling and simulation of complex recovery processes
PDF
Dynamics of water in nanotubes: liquid below freezing point and ice-like near boiling point
PDF
Reactive and quantum molecular dynamics study of materials: from oxidation and amorphization to interfacial charge transfer
Asset Metadata
Creator
Naserifar, Saber
(author)
Core Title
A process-based molecular model of nano-porous silicon carbide membranes
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Chemical Engineering
Publication Date
07/01/2013
Defense Date
06/17/2013
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
amorphous SiC model,diffusion of gases,OAI-PMH Harvest,pyrolysis simulation of HPCS,ReaxFF force filed development,separation of gaseous mixtures,silicon-carbide membranes
Format
application/pdf
(imt)
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Sahimi, Muhammad (
committee chair
), Haas, Stephan W. (
committee member
), Tsotsis, Theodore T. (
committee member
)
Creator Email
naseirfar_saber@yahoo.com,naserifa@usc.edu
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-c3-283510
Unique identifier
UC11294000
Identifier
etd-NaserifarS-1725.pdf (filename),usctheses-c3-283510 (legacy record id)
Legacy Identifier
etd-NaserifarS-1725.pdf
Dmrecord
283510
Document Type
Dissertation
Format
application/pdf (imt)
Rights
Naserifar, Saber
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Tags
amorphous SiC model
diffusion of gases
pyrolysis simulation of HPCS
ReaxFF force filed development
separation of gaseous mixtures
silicon-carbide membranes