University of Southern California Dissertations and Theses

Thermal properties of silicon carbide and combustion mechanisms of aluminum nanoparticle

(USC Thesis Other)

Thermal properties of silicon carbide and combustion mechanisms of aluminum nanoparticle

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content THERMAL PROPERTIES OF SILICON CARBIDE AND COMBUSTION MECHANISMS OF ALUMINUM NANOPARTICLE by Weiqiang Wang A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (MATERIAL SCIENCE) December 2008 Copyright 2008 Weiqiang Wang ii ACKNOWLEDGEMENTS My first and foremost gratitude goes to my Ph.D. advisor, Dr. Priya Vashishta, for his expertise, kindness, and most of all, for his mentoring in my research and life throughout all these years. For all his patient support, guidance, and insight view of my projects, it is difficult to overstate my gratitude to him. I would also like to sincerely thank the other two professors, Dr. Rajiv K. Kalia, and Dr. Aiichiro Nakano, in the group of Collaboratory for Advanced Computing and Simulations (CACS), for continuously encouraging, advising and inspiring me in my research. Without them, this thesis would not have been completed or written. I would also like to thank many of my colleagues in CACS for helping me and providing me a stimulating and fun environment to work. In particular, I would like to thank Dr. Nomura, who has been very kind to help me and many other CACS members and provided us great research and dinning information. Richard Clark, as my close collaborator on the research project, has provided me numerous helps in editing my thesis. My grateful thanks also go to Dr. Paulo Branicio, Dr. Satyavani Vemparala, Dr. Ashish Sharma, Dr. Cindy L. Rountree, Dr. Cheng Zhang, Dr. Zhen Lu, Yi-Chun Chen, Hsiu-Pin Chen, Richard Seymour, Hikmet Dursun, Liu Peng, and Zaoshi Yuan for wonderful discussions, collaborations, and for making the group more cheerful. Special thanks to Sabrina Feeley and Patricia Wong for their coordinating efforts in CACS. Lastly, and most importantly, I would like to express my love and gratitude to my family and friends. In particular, I sincerely thank my parents for their boundless love, unconditional support, and unflinching faith in the way I took to fulfill my graduate study. It is to them that I dedicate this work. iii TABLE OF CONTENTS ACKNOWLEDGEMENTS ............................................................................................ ii LIST OF TABLES .......................................................................................................... v LIST OF FIGURES ....................................................................................................... vi Abstract ......................................................................................................................... xv Chapter 1 ......................................................................................................................... 1 1.1 Silicon Carbide Ceramic ................................................................................. 1 1.2 Thermal Conductivity in Theory, Experiments and Simulations ................... 4 1.3 Combustion of Energetic Materials ................................................................ 8 1.4 Thermite ........................................................................................................ 11 1.5 Metastable Intermolecular Composite .......................................................... 14 Chapter 2 ....................................................................................................................... 23 2.1 Classical Mechanics in Molecular Dynamics ............................................... 23 2.1.1 Formulations of Equations of Motions ................................................. 25 2.1.2 Numerical Techniques .......................................................................... 27 2.2 Statistical Mechanics in Molecular Dynamics .............................................. 40 2.2.1 Statistical Ensembles ............................................................................ 41 2.2.2 Physical Quantities................................................................................ 46 2.3 Potential Forms Used in Molecular Dynamics ............................................. 51 2.3.1 An Effective Interatomic Interaction Potential for Ceramics ............... 52 2.3.2 EAM Potential for Metals ..................................................................... 57 2.3.3 Coupling of Two Different Potential Forms ......................................... 62 Chapter 3 ....................................................................................................................... 73 3.1 Thermal Conductivity Simulation Methodologies ........................................ 73 3.1.1 Equilibrium MD Method ...................................................................... 74 3.1.2 Non-Equilibrium MD Method .............................................................. 75 3.2 Thermal Properties of 3C-SiC ...................................................................... 79 3.2.1 Specific Heat of 3C-SiC ........................................................................ 80 3.2.2 Thermal Conductivities of 3C-SiC ....................................................... 81 3.3 Thermal Properties of a-SiC ......................................................................... 88 3.3.1 Preparation of amorphous SiC .............................................................. 89 3.3.2 Specific Heat of amorphous SiC ........................................................... 90 3.3.3 Thermal Conductivity of amorphous SiC ............................................. 91 Chapter 4 ....................................................................................................................... 97 4.1 Combustion of Aluminum Nanoparticles with Crystalline Shell ................. 98 4.1.1 Preparation of Aluminum Nanoparticles .............................................. 98 iv 4.1.2 Results and Discussion ....................................................................... 100 4.2 Combustion of Aluminum Nanoparticles with Amorphous Shell .............. 124 4.2.1 Preparation of Aluminum Nanoparticle with Amorphous Shell ......... 124 4.2.2 Results and Discussion ....................................................................... 127 Chapter 5 ..................................................................................................................... 154 References: .................................................................................................................. 159 v LIST OF TABLES Table 2.1 Comparison of physical properties of SiC by different potential models .......54 Table 2.2 Parameters in two- and three-body potential terms used in the MD simulation of SiC. ...........................................................................................54 Table 2.3 Comparison of physical properties for a-Al 2 O 3 between MD model and experiments ....................................................................................................55 Table 2.4 Parameters in two- and three-body potential terms used in MD simulation of Al 2 O 3 ..........................................................................................................55 Table 2.5 Comparison of physical properties for Al between EAM model and experiments ....................................................................................................61 Table 2.6 Comparison of bond length ( Ǻ) and bond angles (deg) in Al X O Y molecules ........................................................................................................69 Table 2.7 Comparison of bond length ( Ǻ) in Al 2 O molecules in MD simulations with different value of α R ..............................................................................72 Table 3.1 The zero pressure volumes of the 3C-SiC systems at various temperatures. ..................................................................................................86 Table 3.2 Comparison of bulk thermal conductivities at different temperatures from extrapolation results by using different strategies ..........................................87 vi LIST OF FIGURES Figure 1-1 The stacking sequence of double layers of different SiC polytypes9. ............3 Figure 1-2 High-resolution TEM images of α–SiC before (a) and after (b) irradiation with ions at 300K10. .....................................................................4 Figure 1-3. Schematic of shock and equilibrium Hugoniot curve and Rayleigh line on the P-V plane32. P-V plane is divided into four areas by the dashed lines, where the upper left .............................................................................11 Figure 1-4 Photos of agglomerates (a) and original aluminum (b) (exerted from the research result by Frolov, etc35). ..................................................................13 Figure 1-5 (a) Drainage of aluminum from the oxide “skins” on a platinum wire above the metal melting point; (b) Exudation of aluminum and collapse of the original oxide skin during the subsequent cooling in laboratory heating tests34. ..............................................................................................13 Figure 1-6 Evidence of a detached flame envelope and accumulation of molten oxide on the aluminum droplet surface observed for quenched droplet34. ..14 Figure 1-7 (a) HRTEM of Fe2O3/Al xerogel nanocomposite; (b) SAED pattern of the labeled Al particle in (a) 51. ....................................................................17 Figure 1-8 High-resolution transmission electron micrograph of nano-aluminum53. ...18 Figure 1-9 Schematic of the comparison of between the diffusion mechanism for the oxidation of micro-sized aluminum particle and the melt-dispersion mechanism for the oxidation of nano-sized aluminum particle55. ...............20 Figure 1-10 (Left) Schematic of the reaction of ALEX with oxidizer and the condition for the two nanoparticle’s merge under the condition of flash heating by laser; (Right) experimental results of specimens exposed to weaker, strong and even stronger pulse (from top to bottom) resulting in complete ignition54. .....................................................................................21 Figure 2-1 Cell decomposition of a system for linked list algorithm. Atom i resides in cell 8. .........................................................................................................31 Figure 2-2 Data structure for linked list algorithm, where “E” represents empty. .........32 Figure 2-3 Schematic of periodic boundary condition algorithm with an example of a single SiO2 molecule 2-dimensional system. The silicon atom in the system is represented by a solid blue circle, and the oxygen atoms in the vii system are represented by solid red circles. Dashed squares and circles represent the replicated systems. ...................................................................35 Figure 2-4 Inter-processor communication for data exchange in 2D example. (a) At step 1 and 2, CPU 4 get boundary atom information from CPU 1 and CPU 7, respectively. (b) At step 3 and 4, CPU 3 and CPU 5 send the information of boundary atoms to CPU 4, together with the information of atoms at the corners label as Blue. ...........................................................38 Figure 2-5 Schematic of primary and secondary boundary shown as red and black solid squares. Red and Black arrow represents the primary and secondary cutoffs. .........................................................................................39 Figure 2-6. Variation of atomistic-level plane stress with respect to the increasing of strain in the ligament region between two voids. .........................................48 Figure 2-7. Energy per atom (top) as a function of temperature. The vertical dashed line is the experimental reported melting temperature. ................................56 Figure 2-8. Specific heat for alpha and amorphous alumina at two densities. ...............57 Figure 2-9 Equations of state of aluminum metal by EAM under the conditions of both NPT ensemble (a) and NVE ensemble (b). The red vertical dash line in (a) is used to indicate where the experimental melting temperature (933 K) is located. .....................................................................62 Figure 2-10. Plot of the weighting factor as a function of the oxidation degree o i n ......64 Figure 2-11. Plot of ) ( ik r Θ as the function of distance between a pair of Al and O atoms .............................................................................................................65 Figure 2-12. (a) Energy fluctuations of the aluminum nano-particle system with various time steps tested. (b) Plot of the linear proportionality between the standard deviation of energy with respect to the square of the time step size. ........................................................................................................67 Figure 2-13. Stability test for the interface between Al core and Al2O3 shell when changing the parameter α R from 3.0 Ǻ (a) to 3.5 Ǻ (b). ..............................68 Figure 3-1 Schematic description of the velocity-scaling method. System is divided into slices of equal size. Velocities of atoms in the two slices positioned at ¼ and ¾ of the system length are scaled to remove and add heat energy into the system respectively. .............................................................76 viii Figure 3-2 (a) The calculated vibrational phonon density of states of 3C-SiC by molecular dynamics using our potential. (b) Temperature dependence of specific heat (constant volume) of 3C-SiC as a result of integration using the d.o.s from (a). The classical limit of the specific heat (1.2423 J/gK) is plotted as a straight line in (b) ...................................................................81 Figure 3-3 Temperature profile along y-axis for a 1,048,576 atoms 3C-SiC system at 700K ..........................................................................................................82 Figure 3-4 System temperature profiles and their linear fits for four consecutive 0.5 ns periods. The differences among the temperature profiles can be seen from the slopes of the linear fits. The temperature gradient of the first 0.5 ns period differ by about 47% from the first 0.5 ns period. ....................83 Figure 3-5 extrapolation of thermal conductivities from four different size systems at 700K, the inverse of the intercept gives the bulk thermal conductivity as 59.88(Wm-1K-1). .....................................................................................84 Figure 3-6 The temperature dependence of thermal conductivity, extrapolated from four different sizes of systems. Also shown are experiment and simulation results from another group116 (data is exerted from ref. by J.Li et. al93) ..................................................................................................85 Figure 3-7 The estimated phonon mean free path for 3C-SiC is inversely proportional to temperature. ..........................................................................88 Figure 3-8. The flow chart of the preparation procedure for the amorphous SiC system. ..........................................................................................................89 Figure 3-9 (a) The calculated vibrational phonon density of states of a-SiC by molecular dynamics using our potential, and (b) The temperature dependence of specific heat (at constant volume) of a-SiC as a result of integration of the d.o.s from (a). The specific heat of 3C-SiC is also shown again for comparison with that of a-SiC ............................................91 Figure 3-10 Temperature profile for an a-SiC system at 700 K with steady heat flow established along the y-axis. .........................................................................92 Figure 3-11 Temperature profile along the y-axis of an a-SiC system with four heat sinks and sources placed at each 1/8 position along the y-axis. ...................92 Figure 3-12 Extrapolation of the thermal conductivities for a-SiC at 900 K and 50 K from three differently sized systems. The linear-fit result shows that the extrapolated bulk thermal conductivity for 900 K a-SiC is not very different from the result in the largest sample. However, for the system ix at 50 K, the extrapolated result is considerably different from the result in the largest sample ......................................................................................93 Figure 3-13 Comparison of the temperature gradient between two consecutive 0.5 ns simulations of a 500K a-SiC system (with four heat sources and sinks) using the first region between heat source and sink. The results show the two linear fits, y1=484.20556+3.16712x (for 2nd 0.5ns) and y2=483.17733+3.15441x (for 1st 0.5ns), have only a 0.4% difference in the temperature gradient (slope) ...................................................................94 Figure 3-14 Temperature dependence of thermal conductivities of our system with an effective length of 27.9 nm. .....................................................................95 Figure 4-1 Oxidation of aluminum nanoparticle color coded by the local temperature ...................................................................................................97 Figure 4-2. (a) Initial configuration of the system with oxygen atoms shown by red color and aluminum atoms shown as white. (b) Mix-up of the core and shell at the interface after 10ps of thermalization. ........................................99 Figure 4-3 (a) Increase of kinetic energy per aluminum atom in the process of explosion for the C3, C6 and C9 systems. (b) Survival fraction of unoxidized aluminum atoms along with time. ............................................101 Figure 4-4. Snapshots of explosion status of the systems at different time steps with different initial core temperature are shown by color coding the core aluminum with white, the shell aluminum with yellow, the shell oxygen with red and the environmental oxygen with blue. (a) expansion snapshots for C3 system at different time steps; (b) expansion snapshots for C6 system at different time steps; (c) expansion snapshots for C9 system at different time steps. .....................................................................105 Figure 4-5. (a) Plot of the outside radius of the nanoparticle in C3, C6 and C9 systems along with time. (b) Change of shell’s minimum thickness of the three systems along with time. ..............................................................106 Figure 4-6 Schematic of events happened in the process of aluminum nanoparticle’s detonation. Red solid represents the oxygen atoms and grey solid represent the aluminum atoms. Left white region is the core. Center blue is the Shell region. Right light purple is the outside region of the nanoparticle. Atoms moving are shown as the ones with shadow. Six events of exchange of atoms between shell with either core or environment are Al(s → e) , Al(s → c) , O(s → e) , O(s → c) , Al(c → s) , and O(e → s) . There are also two events Al(c → e) and O(e → c) , that represent the direct transport of atoms between core x region and environment region when there is a channel present (shown as the light green region in the shell). .........................................................108 Figure 4-7. Deterioration of the shells in C3, C6 and C9 nanoparticle system is quantitatively measured by (a) the number of shell aluminum atoms that jet out of the nanoparticle and (b) that fly inside the nanoparticle. ............109 Figure 4-8. Recovery of the shells in C3, C6 and C9 nanoparticle system is measured quantitatively by (a) counting the number of environmental oxygen atoms that oxidize the shell and (b) the number of core aluminum atoms that react with oxygen atoms and deposit on the shell. ...111 Figure 4-9. Increasing of the total number of atoms on the shell in C3, C6 and C9 systems. .......................................................................................................112 Figure 4-10. Change of the shell’s composition along with time in C3, C6 and C9 systems by calculating the ratio of number of aluminum atoms on the shell to number of oxygen atoms on the shell. ............................................113 Figure 4-11. Semi-log plot of the mixing of core aluminum atoms and environment oxygen atoms in C3, C6 and C9 system is measured quantitatively by counting (a) the number of core aluminum atoms that jet into the out- particle environment and (b) the number of environmental oxygen atoms that fly into the core. ...................................................................................114 Figure 4-12. Snapshots of shell morphology of C6 system (a) and C9 system (c) at 100ps with colors representing the number density of the shell as shown by the color bar. Snapshots of the shell morphology together with yellow color core aluminum atoms are shown in (b) for C6 system and (d) for C9 system. Environmental oxygen atoms that are not part of the shell are excluded from the snapshots. Core aluminum atoms are seen jet out of the nanoparticle from the weak areas on the shell during the expansion. ...116 Figure 4-13. (a) Statistics of size distribution of fragments in the C3 system at 200ps shows that the majority of the fragments are small size fragment that has number of atoms in the fragment less than 6. (b) Statistics of fragment composition (ratio of number of aluminum atoms to number of oxygen atoms in a fragment) for fragments with number of atoms less than 6 in C3 system at 200ps. Majority of the small fragments are found to be aluminum rich fragments. ..................................................................119 Figure 4-14. (a) Statistics of size distribution of fragments in the C6 system at 200ps shows that the majority of the fragments are small size fragment with total number of atoms in the fragment less than 6 and the number of small fragments has as much 5 times of that in C3 system (b) Statistics xi of fragment composition (ratio of number of aluminum atoms to number of oxygen atoms in a fragment) for fragments with number of atoms less than 6 in C6 system at 200ps. As in C3 system, most of the small fragments are AlO or aluminum rich fragments, but there are also many oxygen rich fragments.................................................................................120 Figure 4-15. (a) Statistics of size distribution of fragments in the C9 system at 200ps shows that most of the fragments are small sized fragments and the number is almost tripled as compared to the number in C6 system at 200ps. (b) Statistics of fragment composition (ratio of number of aluminum atoms to number of oxygen atoms in a fragment) for fragments with number of atoms less than 6 in C9 system at 200ps. The majority of the fragments are not aluminum rich fragments as contrary to the result in C3 and C6 system, and significant amount of the small fragments are oxygen rich fragments (AlO2 and AlO3). ...........................121 Figure 4-16. (a). Snapshots of the cross section of the center for C3 (a), C6 (b) and C9 (c) system with colors representing the temperature profile as shown by the color bar. The C3 system shows a medium temperature sphere localized to the nanoparticle only. The C6 system has high temperature range expanded to a neighborhood of the nanoparticle, but the high temperature regime drops to a medium temperature within twice of the nanoparticle’s radius. The C9 system has high temperature regime well expanded to a much larger sphere with radius as large as 5 times the radius of the original nanoparticle. .............................................................122 Figure 4-17. Temperature profile of the C3, C6 and C9 systems at 200ps. Temperature is averaged over a 2A thickness shell from the center of the nanoparticle. This plot quantitatively shows that with a higher initial core temperature, the system will result a larger and hotter hot region inside the system. The vertical bar is where the outside surface of the nanoparticle located, and they are used merely to guide the viewer. .........123 Figure 4-18. A center slice view for the nanoparticle with crystalline shell (a) and amorphous shell (b). Shell thickness of crystalline shell is about 4nm, and shell thickness of amorphous shell is about 2~3 nm. ...........................125 Figure 4-19. Schematic of the schedule of amorphorizing the crystalline shell of Al nanoparticles. Shell thickness is also reduced in this process to 2~3nm. ...126 Figure 4-20. Examination of the shell structure by checking the structural quantities of the shell, including the partial pair distribution (a) and the bond angle distribution (b). ...........................................................................................127 xii Figure 4-21. (a) Increase of kinetic energy per aluminum atom in the process of explosion for the A3, A6 and A9 systems. (b) Survival fraction of unoxidized aluminum atoms along with time for A3, A6 and A9 systems. (c) Comparison of energy release among ANP systems with crystalline shell and amorphous shell. (d) Comparison of the Al atom consumption among ANPAS and ANPCS systems. Difference among C3, A9 and other systems are also obvious. ...............................................130 Figure 4-22. Snapshots of explosion status of the ANPAS systems at two time steps with different initial core temperature are shown by color coding the core aluminum with white, the shell aluminum with yellow, the shell oxygen with red and the environmental oxygen with blue. (a) expansion snapshots for A3 system at 70 and 200 ps; (b) expansion snapshots for A6 system at 70 and 200 ps; (c) expansion snapshots for A9 system at 70 and 200 ps. At 200ps, the A9 system shows dramatic difference in its configuration from the other two, indicating a different oxidation mechanism. .................................................................................................133 Figure 4-23. (a) Radial plot of the number density distribution for A3, A6 and A9 systems. (b) Radial plot of the composition distribution for A3, A6 and A9 systems. .................................................................................................134 Figure 4-24. (a) Plot of the outside radius of the nanoparticle in A3, A6 and A9 systems along with time. (b) Change of shell’s minimum thickness of the three systems along with time. (c) Comparison of nanoparticle radius expansion among all ANPAS and ANPCS system. (d) Comparison of minimum shell thickness among all ANPAS and ANPCS systems. ..........137 Figure 4-25. Comparison of the total number of covalently bonded fragments among all the ANPAS and ANPCS systems. .............................................139 Figure 4-26. (a) Statistics of fragment composition (ratio of number of aluminum atoms to number of oxygen atoms in a fragment) for fragments with number of atoms less than 6 in A3 system at 200ps, where majority of the small fragments are AlO and AlO2. (b) Statistics of fragment composition for fragments with number of atoms less than 6 in C3 system at 200ps, where majority of the small fragments are found to be aluminum rich fragments. ...........................................................................140 Figure 4-27. (a) Statistics of size distribution of fragments in the 6K system at 200ps shows that the majority of the fragments are small size fragment with total number of atoms in the fragment less than 6 and the number of small fragments has as much 5 times of that in 3K system (b) Statistics of fragment composition for fragments with number of atoms less than 6 in 6K system at 200ps. As in 3K system, most of the small fragments are xiii AlO or aluminum rich fragments, but there are also many oxygen rich fragments.....................................................................................................141 Figure 4-28. (a) Statistics of size distribution of fragments in the 9K system at 200ps shows that most of the fragments are small sized fragments and the number is almost tripled as compared to the number in 6K system at 200ps. (b) Statistics of fragment composition for fragments with number of atoms less than 6 in 9K system at 200ps. The majority of the fragments are not aluminum rich fragments as contrary to the result in 3K and 6K system, and significant amount of the small fragments are oxygen rich fragments (AlO2 and AlO3). ..................................................141 Figure 4-29. Chang of the total number of atoms on the largest fragments of the three ANPAS system. .................................................................................144 Figure 4-30. Deterioration of the shells in A3, A6 and A9 nanoparticle system is quantitatively measured by (a) the number of shell aluminum atoms that jet out of the nanoparticle ) ( e s n Al → and (b) that are left inside the nanoparticle ) ( c s n Al → . ............................................................................146 Figure 4-31. Oxidation of the shells in A3, A6 and A9 systems are measured quantitatively by (a) counting the number of environmental oxygen atoms that oxidize the shell ) ( s e n O → and (b) the number of core aluminum atoms that react with oxygen atoms and deposit on the shell n Al (c → s) ...................................................................................................146 Figure 4-32. (a) Overall view of the shell structure of the A6 system at 200 ps. Color coding is as following: red represents the oxygen atoms from shell; blue is the oxygen from environment; yellow is the Al atoms from shell; white is the Al atoms from core. (b) Atomistic view of a part of the shell from the left picture (a) shows the complexity of the shell structure is consist of aluminum oxide fragments mixed with pure aluminum metal clusters. Bond length cutoff is 2.5 Å and the color of red represent all oxygen atoms and color of white represent all aluminum atoms. .........................................................................................147 Figure 4-33. Semi-log plot of the mixing of core aluminum atoms and environment oxygen atoms in A3, A6 and A9 system is measured quantitatively by counting (a) the number of core aluminum atoms that jet into the out- particle environment ) ( e c n Al → and (b) the number of environmental oxygen atoms that fly into the core ) ( c e n O → . (c) Plot of ) ( e c n Al → for both ANPCS and ANPAS systems. (d) Plot of ) ( c e n O → for both ANPCS and ANPAS systems. ....................................................................149 xiv Figure 4-34. Snapshots of shell morphology of C9 system (a) and A9 system (c) at 100ps with colors representing the number density of the shell as shown by the color bar. Snapshots of the shell morphology together with yellow color core aluminum atoms are shown in (b) for C9 system and (d) for A9 system....................................................................................................151 Figure 4-35. Snapshots of the cross section of the center for A3 (a), A6 (b) and A9 (c) system with colors representing the temperature profile as shown by the color bar. With the increase of the initial core temperature, a larger and hotter spot is found. ..............................................................................152 Figure 4-36. (a) Temperature profile of the A3, A6 and A9 systems at 200ps. (b) Temperature profile comparison for ANPCS and ANPAS system. Temperature is averaged over a 2A thickness shell from the center of the nanoparticle. The vertical bar is where the outside surface of the nanoparticle located, and they are used merely to guide the viewer. .........153 xv Abstract Silicon Carbide (SiC) and Alumina (Al 2 O 3 ) are two important engineering materials for their outstanding electronic, thermal and physical properties. Applications of SiC to high-temperature and high-power devices require knowledge of its thermal properties. We have studied the specific heat and thermal conductivities of both cubic SiC and amorphous SiC as functions of temperature. Results show different heat transport behaviors for crystalline and amorphous SiC. A layer of about 4 nm Al 2 O 3 is naturally formed on the surface of aluminum nanoparticles (ANPs). Presence of this layer affects the combustion behavior of ANPs. Using multi-million-atom molecular dynamics (MD) simulations, we have investigated the combustion behavior of a single ANP. Transition of the combustion mechanism from diffusion to ballistic transportation of atoms is found when the temperature of the preheated aluminum core of the ANP is changed from 3000K to 9000K. Higher initial temperature of the aluminum core results the mechanical breakdown of the alumina shell in the expansion phase of the ANP. Breaking of the shell provides direct oxidation paths for core aluminum atoms, which results faster oxidation reaction and then faster energy release. Thus, the mechanism of the combustion is transferred from thermodynamics to mechanochemistry regime. Furthermore, when the shell structure of the ANP is changed from thicker crystalline to thinner amorphous, this transition of combustion mechanism is found at lower temperatures. Specially, for initial core temperature of 9000K, the ANP with amorphous shell explodes and results complete shattering of the shell. No large fragment of oxide is found after the explosion except for clusters of tens of oxidized Al atoms, that form a superficial particle and reactivity of which is greatly enhanced. 1 Chapter 1 Introduction Ceramics are important materials because they have been not only utilized widely in electronic devices, fusion reactors, and armors, but also involved in energetic materials as either desensitizers or combustion products. For example, silicon carbide (SiC) is considered as one of the most prominent candidates for semiconductor electronic devices under circumstances of high temperature, high power, and high radiation 17 . In section 1.1 and 1.2, backgrounds about SiC and thermal conductivity will be introduced. An example of ceramic usage in energetic materials is the alumina skin of aluminum nanoparticles. This naturally grown alumina skin of aluminum nanoparticle desensitize the reaction of aluminum nanoparticles with oxidizers, and can affect the performance of the aluminum nanoparticles by varying its thickness 106 . Theories and experiments related to combustion, thermite and thermite on nanoscales will be introduced in section 1.3 to 1.5. 8.1 Silicon Carbide Ceramic Semiconductor devices are being widely used in people’s daily life. Among all the semiconductor materials, silicon carbide (SiC) is one of the most prominent compounds, which has significantly important potential in optoelectronic and environmental applications due to its merit of combining excellent semiconducting and mechanical properties, i.e., good chemical stability and electronic properties, high stiffness and hardness, light weight, high strength at high temperature, low thermal expansion, high thermal conductivity, wide band-gap, low level of induced radioactivity, resistance to 2 oxidation, corrosion and creep at high temperature 17,36,109,110 . As candidates for electronic device applications, SiC is prominent due to the breakdown electric field strength and the high saturated drift velocity of the electrons. In these applications, thermal conditions always have significant effect on the devices’ performance. Thus, thermal management of SiC on various scales and composite forms are of great importance to scientists now. SiC has over two hundred known polytypes, which differ from each another only in the stacking sequence of the Si-C double layer. But only a few are commonly grown and acceptable for electronic semiconductor uses, including 3C-SiC (also known as β–SiC), 4H-SiC and 6H-SiC (also known as α–SiC). Schematics of SiC polytypes are shown in Fig. 1, where 3C-SiC is in structure of cubic with ABC stacking sequence of the double layer, 2H-, 4H- and 6H- SiC in structure of hexagonal with AB, ABAC and ABCACB stacking sequences, respectively. Over years, 3C-SiC has attracted more interest for microelectronmechnical system applications than other polytypes, and experimentally, there are great interests in growing large size bulk and thick film 3C-SiC 16,73 . The floating zone method is used in bulk growth of SiC from liquid to enhance the growth rate, and chemical vapor decomposition (CVD) method is the most prominent method among epitaxial growth methods, due to higher purity of its products. Other emerging methods for growing high quality bulk 3C-SiC include physical vapor transport 16 and vapor-liquid-solid 87 methods. 3 Figure 8-1 The stacking sequence of double layers of different SiC polytypes 48 . Other than the extensive studies of crystalline SiC, much of the research interests have also been given to the amorphous counterpart of it. Technologically, change of physical properties of SiC due to its amorphization during ion-beam irradiation attracted the attention of researchers 107 . An example is given in Fig. 2, which shows the high resolution transmission electron microscopy (TEM) pictures of SiC structure before and after the irradiation. High-level deposition requirement of SiC layers in the applications of engines, turbines and reactor makes the choice of amorphous SiC (a-SiC) attractive due to the high temperature stability of semiconducting properties 83,109 . Fundamental research interest also arises for studying amorphous semiconductors in order to understand the chemical bond change, the thermal property change, the specification of the network topology, and the causes of different ordering strength in different semiconductor compounds 26,45,94 . 4 Figure 8-2 High-resolution TEM images of α–SiC before (a) and after (b) irradiation with ions at 300K 107 . Numerous studies have been carried out for the application of SiC in various electronic devices, but few of the theoretical advantages of SiC have been realized in device applications 17 . Even few efforts have been put into the study of thermal properties of SiC/a-SiC, although thermal stability of semiconductor devices using SiC requires more attention. Specially, in the high temperature applications, where SiC is superior to Si, the ability to dissipate heat is often the limiting factor that determines device performance. Furthermore, the temperature dependence of thermal conductivities in crystalline and amorphous semiconductors shows different characteristics. To address the above issues, the following section gives backgrounds of the theoretical calculations and experimental measurements of the thermal conductivities as a function of temperature. 8.2 Thermal Conductivity in Theory, Experiments and Simulations Thermal conductivity is one of the fundamental thermophysical properties of materials. For a given temperature gradient between a heat source and a heat sink, thermal conductivity determines the rate of heat transfer from the source to the sink. Different mechanisms of heat transport are found in different types of materials. In gases, 5 heat is transferred through the collisions between molecules. In metal, heat is transferred through the drift of electrons, which are also the carriers of electricity. Unlike gases or metals, heat transportation in crystalline insulators is carried out by “phonon”, a quantized motion of lattice vibration. Lattice thermal conductivity derived from the elementary kinetic theory is detailed in the book by Ashcroft and Mermin 6 , which states that the thermal conductivity κ is proportional to the product of specific heat v c , phonon velocity c and phonon mean free path l . l c c v ⋅ ⋅ = 3 1 κ (1.1) By considering both the normal and umklapp collision processes of phonons, the elementary kinetic theory explains quite nicely the temperature dependence of thermal conductivities in crystalline insulators, spanning the temperature range from very low temperature to as high as above Debye’s temperature. At extremely low temperature (about several K), the umklapp processes are almost “frozen out”, and temperature dependence of thermal conductivity only comes from the specific, thus proportional to T 3 . The increasing of κreaches maximum when the exponentially increasing umklapp process catches up and the increase of specific heat slows down, when κ starts to decrease exponentially. At temperature close and above the Debye’s temperature, κ decreases by a slow power law due to the increasing number of phonons. However, when crystalline insulators are amorphized, heat conduction shows distinctive features from that in its crystalline counterpart. First, κof amorphized insulators is usually order of magnitude lower than that of crystalline. Second, κ of glasses of a wide range of composition is similar in magnitude and shows universal 6 temperature dependent features. i.e. κ of glasses initially increases as T 2 at very low temperatures (T ≤ 1K), after which κ exhibits a “plateau” regime at temperature in between 2K to about 30K, and then it rises again smoothly to a saturated value. These anomalous features of κ of glasses in the low temperature regimes are explained by models including scattering from structural disorder 66 , tunneling interactions 47 , and soft- potential model 30 , etc. κ of glasses above the plateau is found to have a minimum value, and it has been modeled by Cahill et al. 14 with modified Einstein oscillators. With the decreasing trend of the electronic device size, thermal management at scale of nanometers creates a demand for understanding the thermophysical properties of materials at correspondingly small scales. Challenge rises when the scale of the studied material becomes comparable to the wavelength of the phonons 13 . Issues that need to be addressed include the definition of temperature in nonequilibrium nanoscale systems. Thermal conduction at solid-solid interfaces has also become important. So far, there is no systematic model that can be used to explain all these questions. Computational approaches including numerically solving Fourier’s law, calculations of Boltzmann transportation equation and atomistic level molecular dynamics (MD) appear to be useful tools in helping to answer these questions. The traditional lattice dynamics method 57 has advantages to quickly study large systems. This method considers the heat transfer and interaction of phonons explicitly, thus certain input parameters such as phonon density of states and phonon relaxation time, which might be obtained from either experiments or other simulations are required. On the other hand, atomistic simulations 43,80 are limited to the study of small system. But the application of it does not need a sophisticated understanding of the fundamental heat 7 transport process of a particular material. Given the interatomic potential of a special system, interesting physical properties can then be generated naturally through simply evolving the system by Newton’s equations. Thus, this method is more suitable for studying of thermal properties of new materials, interfaces and materials in nanoscale regime, which either has less heat transport process understood or the novel properties can not be extrapolated from the known properties of bulk systems. Two most popular methods are the equilibrium MD and the non-equilibrium MD. Details of these methods and their applications in studying the thermal conductivity of both 3C-SiC and a-SiC will be discussed in Chapter 3. Parallel to the development of modeling of the thermal conductivity, experimental methods have also been advanced to probe the thermal conductivity with more accuracy and at finer scales. Conventionally, the method of measuring thermal conductivity uses separate heat source and thermometer that requires big sample size and needs long time to establish the equilibrium heat flow state. Methods of using single element as heat source and thermometer in the time domain 34 were later developed, as well as methods of using the frequency domain technique while having separate heat source and thermometer 88 . An approach of combining these two techniques proposed by Cahill et al. is the 3 ω method 15 , which provides reliable thermal conductivity data over a wide temperature range due to the elimination of the blackbody radiation error in the measurement. Initially, the 3 ω method was developed for studying the thermal conductivities of homogenous materials, such as glasses. Later, the 3 ω method has been widely implemented for measuring the cross-plane thermal conductivity of thin films 46 , 8 and then extended for measurement of the in-plane and cross-plane thermal conductivities of anisotropic films 10 . Apart from the 3 ω method, there are many other methods that measure the thermal conductivity either directly or indirectly. For example, the axial flow method for low temperature samples and the guarded hot plate method for bigger samples measure the thermal conductivity directly 39 . An example of indirect measurement method is the flash diffusivity method 38 , which measures the thermal diffusivity and the thermal conductivity can be determined from the following relation: ρ κ aCp = (1.2) where a is the thermal diffusivity, Cp is the specific heat, and ρ is the density. In this thesis, we will present the work on simulations of thermal conductivity of both 3C-SiC and a-SiC with NEMD method, results and discussions of which are detailed in Chapter 3. 8.3 Combustion of Energetic Materials Energetic materials are the substances that store large amount of chemical energies and release them rapidly upon certain stimuli. Typically, energetic materials are classified into propellants, pyrotechnics, and explosives. Depending on the preparation processes, energetic materials may consist of either a chemically pure compound or a composite of a fuel and an oxidizer. The chemically pure, or the monomolecular, energetic materials are prepared by incorporating the oxidizing and fuel moieties into one molecule and they are generally explosives, which can be ignited and explode to produce detonation front waves of exothermic reactions. The high energy release power of these monomolecular 9 energetic materials is primarily due to the intramolecular reactions controlled by chemical kinetics. But the total achievable energy density is limited due to the requirement for a chemically stable material and the feasibility of making the density of the materials higher. Examples of the monomolecular energetic materials include trinitrotoluene, nitramine, TNT and RDX, etc. On the other hand, energetic materials of fuel and oxidizer composition are prepared through physical mixing of solid oxidizers and fuels and can be tailored through different proportions to be either explosives or pyrotechnics and propellants. Pyrotechnic compositions and propellants, like the explosives, involve exothermic reactions and produce heat, light and/or thrust, but non detonative. Energy densities in the composite energetic materials can be significantly higher than that stored in monomolecular energetic materials. Maximum energy density can be reached by balancing the ratio between the oxidizer and fuel. However, the granular nature of composite energetic materials determines the mass transport rates between reactants as the reaction kinetics, and diffusion is the mechanism of mass transport for composite with fuel size on the scale of micrometers. Hence, the energy release rate of the composites is usually lower that of monomolecular energetic materials, even though they may have higher energy densities. Burning process of the energetic materials can be explained by the well developed combustion theory, which involves not only the elementary thermodynamics, but also the conservation equations of fluid dynamics and chemical kinetics due to the flow of reacting and diffusing gases. Analysis of the combustion process would include the overall continuity equation, momentum conservation, energy conservation, continuity of 10 species, radiant heat transfer, etc. An important result of the combustion theory is the Rankine-Hugoniot relations, within which plane deflagration and detonation waves may be investigated. Deflagration and detonation are two general branches of combustion with characteristic reaction speeds of subsonic and supersonic, respectively. Deflagration combustion features laminar propagation flames with a subsonic speed and usually propagates by igniting the cold material through conducted heat from the ignited region. Detonation combustion features propagation of shock wave with a supersonic speed through the material, which is compressed and ignited and then releases the energy to further support the propagation of the shock wave. Difference of the two classes of combustion can be seen from the plot of Hugoniot curves, where two nonoverlapping regions are representations of the two combustion branches respectively (see Fig. 3). Detailed derivation of the Rankine-Hugoniot relations can be found in the book by F. Williams 108 . 11 Figure 8-3. Schematic of shock and equilibrium Hugoniot curve and Rayleigh line on the P-V plane 28 . P-V plane is divided into four areas by the dashed lines, where the upper left corner represents the detonation branch and right lower corner represents the deflagration branch. 8.4 Thermite Thermite, the composite energetic material consisting of aluminum powder and a metal oxide, can create short bursts of extremely high temperature in the aluminothermic reaction, also known as the thermite reaction. The discovery of thermite reactions is attributed to the German chemist Hans Goldschmidt 31 . In the process of thermite reactions, the oxides are reduced and the aluminum is oxidized to form aluminum oxide. This solid-state redox reaction can be formally described as the following: H A MO AO M Δ + + → + (1.3) 12 where M refers to the aluminum metal, A can be either metal or non-metal and H Δ is the heat generated. In principle, M can be any reactive metal, but aluminum is the favored one due to its merits of large availability, low melting but high boiling temperatures, and the insensitivity to ambient environment because of the passivation layer on it. Metals as fuel ingredients in explosives and propellants have been studied extensively. The most commonly used and also studied metal is again aluminum. Researchers have found that the combustion behavior of metal powders (with the size of the metal particles usually in the range of 10-40 micrometers) is characterized by the accumulation of metal on the burning surface, which results in aggregated, relatively large agglomerates that burn after leaving the burning surface. A picture of typical aluminum particles as propellant ingredient and the observed agglomerate is shown in Fig. 4. The reason of this fusion of individual metal particles is attributed to the physical properties of the metals and their oxides. The considered properties of metals and their oxides include: thermal expansion coefficients, melting points, solubility to each other, and the boiling points. Particles are quenched in the process of burning to study microscopically what happened to the nanoparticle in the burning process. Drainage or exudation of aluminum (see Fig. 5), accompanied by the collapse of the oxide skins, was observed in experiments when aluminum particles are continuously heated above the aluminum melting point in inert atmosphere 50 . 13 Figure 8-4 Photos of agglomerates (a) and original aluminum (b) (exerted from the research result by Frolov, etc 27 ). Figure 8-5 (a) Drainage of aluminum from the oxide “skins” on a platinum wire above the metal melting point; (b) Exudation of aluminum and collapse of the original oxide skin during the subsequent cooling in laboratory heating tests 50 . Ignition of single aluminum particle is manifested by a sudden and drastic increase in luminosity. The burning droplet is found often with a residue of molten “retracted” oxide on the surface and a detached flame envelope as shown in Fig. 6. The oxide found on the surface of the droplet or in the flame may be either a relic of the original oxide skin or that newly formed during the continuous burning of the droplet. 14 Figure 8-6 Evidence of a detached flame envelope and accumulation of molten oxide on the aluminum droplet surface observed for quenched droplet 50 . Applications of thermite reactions have been involved in a wide range of civilian and military uses 104 . However, the traditional energetic materials still feature a slow reaction rate and incomplete burning. The emerging of energy demanding systems has provoked researchers’ efforts into the improvement of the combustion efficiency of thermites 8,78 . Increase of porosity has long been recognized as one of the methods to trigger the burning of energetic materials from slow to violent. Transition of granular explosives from deflagration to detonation has been studied by many researchers 44 . But presence of porosity inevitably derogates the energy density of the energetic material. 8.5 Metastable Intermolecular Composite Classical thermites release enormous heat at burning, but the burning velocity is relatively slow. Recent developments in nanoscale energetic materials have attracted many researchers’ attention due to the merit of their high heats of combustion and fast energy release rates. A promising class of the energetic composite, the metastable intermolecular composite (MIC), is a subclass of thermite and consists of mixtures of fuel 15 and oxidizer particles on the scale of nanometers. Faster combustion velocities have been achieved by introducing metal powder with micro and nano sized particles into traditional energetic materials 32,49 . Specifically, researchers have found that reducing the size of the metallic particles can significantly increase the combustion velocities 42-44 . Reaction propagation can be as high as 1km/s, and the dominant mechanism of this fast propagation is proposed to be convection 5 . To explain the reduced ignition time and higher burning efficiency of MIC, the reaction mechanism for classical thermites are employed. In traditional thermites ，the metallic particles have size of micron meters, and their burning is controlled by the mass diffusion of oxidizer through the oxide shell to the metal fuel and vise versa. Similarly, when it comes to the MIC, shorter diffusion of the fuel to oxidizer is the first intuitive thought for the reason behind this enhanced performance. Modeling by Shimizu and Saitou 84 and experiments by Aumann 8 have both shown that this enhancement of burning could be attributed to the increased contact intimacy of fuel and oxidizer particles. But when the metallic particle size is continuously reduced to the nano-scale and with the improvement of nano-engineering techniques in synthesizing MIC, diffusion mechanism alone had become insufficient for the explanation of such high burning rates of MIC. It was found that the reactivity can be increased by orders of magnitude 44,45,47,48 , which can not be simply addressed by the explanation of intimate contact between fuel and oxidizer, together with the traditional diffusion and oxidization mechanism. Technologies of producing nanoparticles have been developed quickly over the years. Techniques of synthesizing ultra-fine nanopowders can be either chemically or physically. Synthesis of iron nanoparticles via chemical reduction method has been 16 reported, where sodium borohydride, polyacrylic acid, and palladium ions are used as the reducing agent, the dispersing agent, and seeds for iron nanoparticle nucleation, respectively 40 . Iron nanoparticles produced possess an average diameter of 6 nm and a geometric standard deviation of 1.3. Photochemical and radiation-chemical reduction methods are alternative ways of producing high purity nanoparticles of various sizes. Photochemical method has been reported to synthesize polygonal gold nanoparticles capable of absorbing near infrared light (NIR) radiation 41 . In the above chemical reduction methods, stabilizer is always used to prevent nanoparticles from aggregating. In the absence of stabilizers, metallic nanoparticles of high activity tend to aggregate to larger particles even without activation energy. Cryochemical synthesis method realizes the stabilization by using inert gas matrix isolation at superlow temperatures. Other than the above chemical methods, a number of physical methods are used for preparing metal nanoparticles. These methods include those employing low-temperature plasma, gas evaporation, various dispersions, etc. Development of the nano-engineering technology has made it possible to produce energetic composites with precisely controlled fuel-oxidizer compositions and uniform densities through the sol-gel synthetic approach 95 . Recently, development in this sol-gel approach has enabled the incorporation of organic functionality into the energetic materials, a potential application as pyrotechnics and propellants 18 . Figure 7 (a) shows the intimate contact of iron-oxide with ultra-fine grain (~25nm diameter) aluminum particles in the thermite nanocomposite xerogel prepared by this sol-gel method and visualized through the high resolution transmission electron microscopy (HRTEM). Figure 7 (b) shows the confirmation of the aluminum particle labeled in (a) by the selected area electron diffraction (SAED) pattern. 17 Figure 8-7 (a) HRTEM of Fe 2 O 3 /Al exerogel nanocomposite; (b) SAED pattern of the labeled Al particle in (a) 95 . High resolution TEM has also been used by Ramaswamy et al. to reveal the oxide surface and aluminum lattice down to the atomic level 75 . Figure 8, the HRTEM picture of a corner of aluminum nanoparticle, shows that the thickness of the oxide skin of the nanoparticle is about 2.5 nm and the aluminum core has a crystal lattice structure. 18 Figure 8-8 High-resolution transmission electron micrograph of nano-aluminum 75 . With the ability of looking at the structure of the MIC down to atomic levels, researchers are very interested in seeing the combustion details of the MIC. A great deal of experiments and theoretical calculations on energetic materials involving aluminum nanoparticles has been carried out in many different research groups 41,51,54-61 . In experiments of ALEX (a commercially available Al nanoparticle aggregates) mixed with various oxidizers, Ivanov and Tepper 42 noticed that exotherm could occur 19 without oxidation, thus proposed that the source of this burning behavior could be from internal energy. However, experiments by Mench et al. 60 and theoretical calculations by Desena and Kuo 21 have both proved that the amount of energy stored in process of making ALEX is realistically negligible. Researchers in Texas Technology University have studied extensively on combustion behavior of various MICs (Al+MoO 3 and Al+Fe 2 O 3 nanocomposites) and proposed a new melt dispersion mechanism for the fast reaction of Al nanoparticles during fast heating 48,63 . The melt-dispersion mechanism explains the fast reactivity of the MIC by considering the spallation of the oxide shell due to the high pressure of the melted aluminum core when the nanoparticle is heated. Spallation of the oxide shell gives space for the liquid aluminum clusters, which jet out of the core and react with neighboring oxidizers. This reaction pathway is much faster than the diffusion mechanism as seen for micron sized thermites. Figure 9 is the schematic of this mechanism by Pantoya et al. 68 They have calculated the conditions for the spallation of oxide shell, but based on the bulk properties of aluminum and alumina, and there has not been any direct proof of it. 20 Figure 8-9 Schematic of the comparison of between the diffusion mechanism for the oxidation of micro-sized aluminum particle and the melt-dispersion mechanism for the oxidation of nano-sized aluminum particle 68 . In experiments conducted by Dlott’s group in University of Illinois at Urbana- Champaign (UIUC), researchers use laser to flash heat ALEX embedded in either Teflon or NC oxidizers 2,54,64-67 . They have intensively investigated the relationship between laser fluencies and the shock induced chemical reaction range with variables such as ALEX concentration, nanoparticle size, shell thickness etc. They have found the linear relationship between the reaction propagation distance and laser energy, and explained it with hydrodynamic model instead of thermal explosion model. From the observed two- stage consumption of ONO 2 groups in the NC oxidizer, the researchers found that the 21 energy release process consists of two steps. The first fast stage (~300ps) involves the reaction of ALEX nanoparticles with their surrounding shell of NC, and the second slower stage (~2ns) is the reaction of NC in between the nanoparticles triggered by the ignition of the nanoparticles and their shells. A schematic and supporting experimental image results for their explanation of the two-stage reaction between fuel nanoparticles and surrounding oxidizers are shown in Fig. 10. Figure 8-10 (Left) Schematic of the reaction of ALEX with oxidizer and the condition for the two nanoparticle’s merge under the condition of flash heating by laser; (Right) experimental results of specimens exposed to weaker, strong and even stronger pulse (from top to bottom) resulting in complete ignition 22 . To understand the ignitability of a thermite, it is essential to firstly understand the initiation process of individual thermite nanoparticles. Knowledge of a single 22 nanoparticle’s initiation process in general will greatly enhance the understanding of burning process of MIC composites. Analysis of the experiments have suggested that mechanical breakdown of the nanoparticle’s shell plays important role in the initiation process 52 , but details of this mechanism are still unknown. Simulation of the burning process of a single nanoparticle by MD is an appropriate way of starting for the study of nanoparticle’s process at atomic level. In this thesis, we will present the result of multi-million-atom simulations of the detonation process of a single aluminum nanoparticle, where details of the detonation at atomic level have been studied. Effects of initial core temperature and shell structure are studied by a set of detonation simulations for aluminum nanoparticles of crystalline or amorphous oxide shell and the aluminum core at three different temperatures. So, the outline of this thesis is as following. We first discuss the methodologies of MD simulations in Chapter 2. Simulation results of thermal properties of SiC ceramic in both crystalline and amorphous phase are presented in Chapter 3. Detailed results and discussions about the combustion mechanism of aluminum nanoparticles and the effect of shell structures are presented in Chapter 4. In Chapter 5, we will conclude the work. 23 Chapter 2 Molecular Dynamics Methodology As a “virtual experiment” that bridges between laboratory experiments and theoretical postulations, Molecular Dynamics (MD) 1,64,71 has become a well developed interdisciplinary method. Based firmly on classical mechanics and statistical mechanics, it probes the properties of complex systems through computer simulations. Successful applications of the MD method have been reported in many research fields, mainly concentrated in materials science and biological sciences. In these areas, the details of a system’s evolution in time at the atomic scale are often of great interest to researchers. In this chapter, we first introduce the formulation of the equations of motions in the MD method and the algorithms for efficiently solving them numerically on computers. Then we describe the applications of MD to different types of systems and the physical properties it predicts. Finally, we will introduce the different forms of the interatomic potential, the essential input needed by MD, and present examples of its adaptation to our research focus. 10.1 Classical Mechanics in Molecular Dynamics Theoretically, the laws of quantum mechanics more accurately describe the true evolution of a system than classical mechanics. However, the complexity of solving the full time-dependent Schrodinger equation is prohibitive for a system of merely thousands of atoms. On the other hand, the MD method approximates the trajectory of the system by treating the motion of each atom within a classical framework. In this case, each atom 24 is represented as a point particle, the state of which is precisely defined by its position and velocity vectors. This approximation toward classical mechanics becomes reasonable when the de Broglie wavelength is much smaller than the other dimensions of the system. Popular simulation methods that apply classical mechanics include the Molecular Dynamics method, the Monte Carlo (MC) method, the Langevin dynamics method, etc. The way that the MD method applies classical mechanics is to numerically solve Newton’s equations for a statistical ensemble of N atoms, obtaining the system trajectories in the 6N phase space. The Newtonian equation for any atom i in the system is V F r m i i i i ∇ − = = v v & & v (2.1) where i m is the mass of the atom, i r & & v is the second derivative of the atom's coordinates in the system, i F v is the force exerted on the atom and V is the potential energy of the system. Given the force exerted on the atoms, the equations of motion can be integrated to yield the trajectory of the atoms, which in turn determines how the system evolves. Because the solution to the equations of motion is deterministic, given the initial position and velocity of each atom r r i 0 () , r Ý r i 0 ( ) { } , the state of the system r r i t ( ) , r Ý r i t () { } at any later or earlier time t can be predicted. The physical properties of the system can then be obtained by performing statistical averages over the trajectory of the system. Initially, the MD method was implemented in the microcanonical ensemble, which has a fixed number of atoms N, a fixed system volume V and a fixed total energy of the system E (thus the name of “NVE ensemble” used). Application of the MD method was 25 extended to the canonical ensemble by Nose and Anderson to incorporate the isothermal constraint and isobaric constraint, respectively. 10.1.1 Formulations of Equations of Motions Newton’s laws of motions that form the foundation of classical mechanics are also the underpinnings of MD method. Formulations of the equations of motions can follow three ways: the Lagrangian formulation method, the Hamiltonian formulation method, and the operator-based formulation method. 10.1.1.1 The Lagrangian Formulation Method The classical Lagrangian for a system of N particles is ( ) ∑ − = N i N i i i r r r m L ) ,... ( 2 1 2 r r & r φ (2.2) where i m is the mass of atom i , i r & r is the velocity of atom i and ) ,... ( N i r r r r φ is the potential energy of the system. The equations of motion can then be derived directly from the Euler-Lagrange equation: 0 = ∂ ∂ − ∂ ∂ i r L r L dt d r & r (2.3) By substituting equation (2.2) into (2.3), we then obtain the 3N Newton’s second law equations for the system. 26 2.1.1.2 The Hamiltonian Formulation Method Alternatively, Newton’s laws of motions can also be determined through the Hamiltonian formulation. Hamiltonian is a natural function of coordinates and their conjugate momenta defined in terms of the Lagrangian: )) ( , ( ) ( ) , ( 1 p r r L p r p r p H i N i i i r & r r r & r r r r − ⋅ = ∑ = (2.4) Now plugging in the formula for the Lagrangian (equation (2.2)) yields the form: ) ,... ( 2 1 2 N i N i i i r r m p H r r φ + = ∑ = (2.5) The equations of motion of the system are given by Hamilton’s equations: i i i i r H p p H r r & r r & r ∂ ∂ ∂ ∂ − = = , (2.6) which gives precisely the same Newton’s laws of motion as the Lagrangian formulation. 2.1.1.3 The Operator-based Formulation Because the true time trajectories of a Hamiltonian system are symplectic and microscopically reversible, it’s desirable to formulate the equations of motion by a symplectic and reversible integrator. Otherwise one obtains less stable solutions to the long-time integration of Hamiltonian systems. The evolution of a quantity ) , ( p r A r r is given as: {} iLA H A p A r H r A p H dt p r dA p r = = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ ⋅ ∂ ∂ − ∂ ∂ ⋅ ∂ ∂ = ∑ , ) , ( , r r r r r r r r (2.7) where, L is the Liouville operator, defined as, 27 iΓ = {...,H} (2.8) and {...,...} is the Poisson bracket. By plugging in either the position or its conjugate momentum, the above equation gives Hamilton’s equations of motion. Then for a state in the phase space at time t, defined as {} p r t r r , ) ( = Γ , Hamilton’s equations of motion can be written as the evolution of ) (t Γ : {} ) ( ), ( ) ( t iL H t dt t d Γ = Γ = Γ (2.9) which can then be integrated to yield ) 0 ( ) ( Γ = Γ iLt e t (2.10) where {} ) 0 ( ), 0 ( ) ( p r t r r = Γ are the initial conditions. Evolution of the system is then merely the action of the classical propagator iLt e on the state function ) (t Γ . 10.1.2 Numerical Techniques After the formulation of the equations of motion, the time evolution of the system can then be determined numerically by integrating these equations up to the desired time. In addition to the integration algorithms for calculating the trajectories, numerical techniques involved in the MD procedure also include the linked list method for reducing the computation complexity of O(N 2 ) to O(N), multiple time step (MTS) method to take advantage of the different interaction stiffness in parts of a complex potential, the parallelization algorithm to map a large system to a cluster of computing nodes and the periodic boundary condition for treating simulations where surface properties are not wanted. 28 2.1.2.1 Velocity Verlet Algorithm Historically, there have been many integration algorithms, including the Euler algorithm, leap-frog algorithm and velocity Verlet, etc. Among these algorithms, we have selected the velocity Verlet algorithm due to its symplectic and time-reversible properties. The velocity Verlet algorithm can be simply derived by approximating atomic coordinates in a Taylor series expansion: ) ( 2 ) ( ) ( ) ( ) ( ) ( 2 1 ) ( ) ( ) ( 2 2 t F m t t t v t r t t r t t r t r t t r r r r & & r & r r r ⋅ Δ + Δ ⋅ + = Δ ⋅ + Δ ⋅ + = Δ + (2.11) where higher order terms, 3 ) ( t O Δ have been neglected. By approximating the positions in the negative time direction, it becomes: ) ( 2 ) ( ) ( ) ( ) ( ) ( 2 1 ) ( ) ( ) ( 2 2 t t F m t t t t v t t r t t t r t t t r t t r t r Δ + ⋅ Δ + Δ ⋅ Δ + − Δ + = Δ ⋅ Δ + + Δ ⋅ Δ + − Δ + = r r r & & r & r r r (2.12) Solving for the velocity at time t t Δ + for the velocity, we obtain: [ ] ) ( ) ( 2 ) ( ) ( t t F t F m t t v t t v Δ + + ⋅ Δ + = Δ + r r r r (2.13) Equations (2.11) and (2.13) give the formulation of velocity Verlet algorithm. The algorithm is implemented as following: i) Update the velocity of each atom to one half time step ) ( 2 ) ( ) 2 ( t F m t t v t t v r r r ⋅ Δ + = Δ + (2.14) ii) Advance the position of each atom to a full time step 29 t t t v t r t t r Δ ⋅ Δ + + = Δ + ) 2 ( ) ( ) ( r r r (2.15) iii) Calculate the new force on each atom ) ,... ( ) ( N i r r t t F r r v r φ ∇ − = Δ + (2.16) iv) Update the velocity of each atom to the full time step ) ( 2 ) 2 ( ) ( t t F m t t t v t t v Δ + ⋅ Δ + Δ + = Δ + r r r (2.17) In the above equations, t Δ is the time step for the system being updated. The time step size is generally taken as an order of magnitude less than the Einstein period. However, in simulations of fast detonation or hyper velocity impact, it should be properly checked against the energy conservation laws prior to the start of a full simulation. Velocity Verlet algorithm can also be formally derived from the Liouville formalism by approximating the time propagator by the product of unitary exponential operators. First, the time propagator for a discrete time step is defined in: e iLt = e iLΔt [ ] n (2.18) where n t t / = Δ and t p p r r t iL e e Δ ⋅ ∂ ∂ + ∂ ∂ Δ = ) ( r & r r & r To approximate the discrete time operator while retaining the symplectic and reversibility properties, the following approximant of Trotter formula is used: 2 2 ) ( At Bt At t B A e e e e ≅ + (2.19) Similarly, the discrete time propagator can be approximated as: 30 ) ( 3 2 2 ) ( t O e e e e e t p p t r r t p p t r r p p t iL Δ + ≅ = Δ ⋅ ∂ ∂ Δ ⋅ ∂ ∂ Δ ⋅ ∂ ∂ Δ ⋅ ∂ ∂ + ∂ ∂ Δ r & r r & r r & r r & r r & r (2.20) Applying the above operator onto the state vector { } p r r r , simply yields the Taylor expansion of the state at the next time step, with the velocity Verlet form as follows: 2 ) ( ) 2 ( ) ( ) 2 ( ) 0 ( ) ( 2 ) 0 ( ) 0 ( ) 2 ( t t F t p t p t m t p r t r t F p t p Δ ⋅ Δ + Δ = Δ Δ ⋅ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ Δ + = Δ Δ ⋅ + = Δ r r r r r r r r r (2.21) Switching r r r & r ∂ ∂ and p p r & r ∂ ∂ in the above equation gives an alternative algorithm, the position Verlet algorithm, which is also symplectic and time-reversible. 2.1.2.2 Linked List Cell Method The potential energy that describes the interaction between atoms inevitably results in the two-body, three-body and even many-body forces among atoms. Even in the simplest case of a pair potential, the naive implementation of the velocity Verlet algorithm to a system of N atoms will have the computational complexity of O(N 2 ) due to the calculation of forces between each pair of atoms. To improve the scalability of the velocity Verlet algorithm, the potential is first truncated to include only a limited range of interactions for each atom. This range of interaction is predefined in the program as the cutoff length. Potential form is then shifted to eliminate the discontinuity at the cutoff. After this, the whole system is divided into cells which are slightly larger than half the cutoff length. Each cell is assigned an integer number c in between 0 and N C -1 (N C is the total number of cells of the system) as its ID. 31 Atoms in each cell are then “linked” using the linked-list algorithm, and interactions of atoms in a given cell are limited to atoms within the current cell and its neighboring 26 cells (see figure 2.1 for a 2D example). Figure 10-1 Cell decomposition of a system for linked list algorithm. Atom i resides in cell 8. The algorithm of the linked-list cell is implemented using a linked-list data structure (see Figure 2.2) 1 . In this data structure, two arrays are created. The first array is the head array with size equal to N C . Each element of the head array is mapped to a cell by matching its index in the array to the cell’s ID and stores only the head atom of linked list from the cell. The second array is the linked-list array with size equal to the number of atoms in the system. Indices of this array are mapped to the IDs of atoms. Each element of this array stores the ID of the atom that is linked to the atom whose ID is this element’s index. In case of the atom is the tail of the list, it points to the empty “E”. To look up the atoms in a given cell, one needs to firstly find the head atom in the head array and then traverse the “linked” atoms in the cell by searching the linked-list. 32 Implementation of this algorithm reduces the computation complexity to O(N), which makes MD simulations of large systems computationally feasible. Figure 10-2 Data structure for linked list algorithm, where “E” represents empty. To facilitate the searching of neighboring cells, each cell is also assigned a vector ID ) , , ( z y x c c c c = r , which is related to its sequential ID as z cz y cz cy x c L c L L c c + + = , where α c L is the cell size along the α direction. To find a cell’s neighbor, one just need to update its vector ID according to the neighbor’s relative position, and then the scalar ID can be calculated as above. In programming, computational cost can be further reduced in observing that, Newton’s third law allows the calculation of interaction forces being restrained to only half of the neighbor cells. Practically, to further reduce the computational cost, a neighbor list for each atom i is constructed in order to reduce the number of force calculations between atom i and its neighbors. The neighbor list is constructed by searching over the neighboring cells of atom i’s resident cell. Only atoms within the distance of the potential’s cutoff plus a certain skin δ are added to the neighbor list. This neighbor list does not need to be updated unless the atoms in the list move a distance of δ/2 or more. 33 2.1.2.3 Multiple Time Steps Another technique used to enhance the computational efficiency is the multiple- time-step (MTS) method proposed by Streett et al. 89 The basic idea behind the MTS method is that in a MD system, interactions of an atom i to all other atoms in the neighbor list can be further divided into groups of close (also called “primary”) neighbors and far (also called “secondary”) neighbors. Forces between atom i and its primary neighbors generate mostly high frequency motions of atom i, and thus are dominant in determining its motion. Forces between atom i and its secondary neighbors normally are softer and drive slower degrees of freedom. Consequently, the state of the system can be updated for several time steps by using the stiffer forces before it is updated by the softer forces. Formulation of the MTS algorithm can be derived elegantly from the classical integrator such as velocity Verlet. To simplify the derivation, we assume the system potential can also be divided into two terms, V 0 and V 1 , corresponding to the stiffer potential and softer potential, respectively. 1 0 V V V + = (2.22) Then, the Liouvillean can be split as p A r V iL p A r V r A r iL r r r r r & r ∂ ∂ ⋅ ∂ ∂ − = ∂ ∂ ⋅ ∂ ∂ − ∂ ∂ ⋅ = 1 1 0 0 (2.23) The discrete time propagator can then be split and approximated by the Trotter formula as 2 2 ) ( 1 1 1 0 1 1 1 0 1 t iL n n t iL t iL t L L i t iL e e e e e Δ Δ Δ Δ + Δ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = = (2.24) 34 n t t 1 0 Δ = Δ is the time step for the stiffer potential. Factorization of the Liouvillean 0 iL as in the velocity Verlet algorithm gives the MTS propagator (in this case, the multiple time step is a two-fold time step) 2 2 2 2 1 1 0 0 0 0 0 1 1 1 t p r V n t p r V t r r t p r V t p r V t iL e e e e e e Δ ⋅ ∂ ∂ ∂ ∂ − Δ ⋅ ∂ ∂ ∂ ∂ − Δ ⋅ ∂ ∂ Δ ⋅ ∂ ∂ ∂ ∂ − Δ ⋅ ∂ ∂ ∂ ∂ − Δ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ≅ r r r r r & r r r r r (2.25) and the implementation of MTS algorithm is as following ) , ; 0 ( ) 2 ( ) 0 ( 2 ) 0 ( ) 0 ( ) 2 ( 1 1 1 1 + + < = Δ = ′ Δ ⋅ + = Δ i n i i for t p p t F p t p r r r r r 2 ) 2 ) 1 (( ) ) 2 1 (( ) ) 1 (( ) ) 2 1 (( ) ( ) ) 1 (( 2 ) 2 ( ) 2 ( ) ) 2 1 (( 0 0 0 0 0 0 0 0 0 0 0 0 0 0 t t i F t i p t i p t m t i p t i r t i r t t i F t i p t i p Δ ⋅ Δ + + Δ + ′ = Δ + ′ Δ ⋅ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ Δ + ′ + Δ = Δ + Δ ⋅ Δ + Δ ′ = Δ + ′ r r r r r r r r r (2.26) 2 ) ( ) ( ) ( 1 1 1 1 1 t t F t p t p end Δ ⋅ Δ + Δ ′ = Δ r r r 2.1.2.4 Periodic Boundary Condition The periodic boundary condition (PBC) was implemented by Born and von Karman in 1912 to overcome the surface effect problem 11 . In a real three dimensional simulation, the PBC is sometimes desired only in certain directions while surface effects or other boundary conditions (such as the momentum reflection wall used in a recent RDX shock 35 simulation performed by Ken-ichi, et al.) are desired in other directions. By applying the PBC along a given axis, the MD box of the system is replicated infinitely in that direction (see Figure 2.3). In such a system, if an atom moves out of the system in the PBC direction, an image atom from the replica MD box will move into the system. Thus, the conservation of total number of atoms is naturally held by PBC. Figure 10-3 Schematic of periodic boundary condition algorithm with an example of a single SiO 2 molecule 2-dimensional system. The silicon atom in the system is represented by a solid blue circle, and the oxygen atoms in the system are represented by solid red circles. Dashed squares and circles represent the replicated systems. Implementation of the PBC algorithm is simply the minimum-image convention method, and given below using the example of a one-dimensional system. The coordinate of atom i is updated by: ⎪ ⎩ ⎪ ⎨ ⎧ ≤ − > = − − − = 0 2 0 2 ) , 2 ( ) , 2 ( ) , 2 ( i x i x i x x i x i x i i x L x L x L SignR L x L SignR x L SignR x x (2.27) 36 where x L is the size of the MD box along x-axis. Implementation of the PBC can sometimes reduce the required system size when calculating the equilibrium properties of a system that uses short-range interactions in the simulation. By taking advantage of this, the cost of computation can be significantly reduced. Meanwhile, there exist problems in using PBC algorithm instead of simulating a bigger system when studying the dynamic properties of a system. First, events in dynamics could be missed by using a small sample size. Second, the PBC could inhibit the occurrence of long wavelength fluctuations. For example, in our simulations of the thermal conductivity of 3C-SiC the phonon wavelength is much larger than the system size. This makes the simulation results size-dependent. 2.1.2.5 Parallelization Algorithm Even though the computational time of MD can scale as O(N) by using the linked list cell method as mentioned in section 2.1.2.2, in practice the number of atoms in a MD simulation is still small due to the limited memory on a single processor and the necessary storage for atomic trajectories. To achieve large scale MD simulations involving millions or billions atoms in our research, a parallel scheme of spatial decomposition has been used 63,64 . This spatial decomposition method combined with the linked-list cell method can theoretically be used for MD simulations of an arbitrarily large system size, provided enough computing resources are available. In the spatial decomposition method, a simulated physical system is partitioned into small cubes of the same size by using a 3-dimensional mesh. Each cube is assigned a unique scalar ID p, an integer index between 0 and N-1, where N is the total number of 37 cubes in the system. These N subsystem cubes are then assigned to N computing processors to carry out the MD calculations. In our simulations, since the interactions between atoms are truncated at the cutoff length, the MD simulations for each subsystem can be computed locally by each processor provided the positions and velocities of all neighboring atoms within the interaction range are first obtained. In practice, to communicate, a processor needs to know the IDs of its neighbors. As when searching neighbor cells in the linked list cell method, a vector ID r p = (p x , p y , p z ) is given to each processor. This ID is related to the sequential ID by p = p x N py N pz + p y N pz + p z , where α p N is the number of processors along the α direction. For the purpose of inter-processor communication, a boundary layer of width equal to the cutoff distance is drawn around each processor, and the phase trajectories of atoms in that layer are copied to its neighbor at each MD step. In the 3D system, there are a total of 26 neighbors for each processor. The communication with all the neighboring processors is completed in six steps by carrying on the new information from each previous step. An example in a 2D system is given in figure 2.4. 38 Figure 10-4 Inter-processor communication for data exchange in 2D example. (a) At step 1 and 2, CPU 4 get boundary atom information from CPU 1 and CPU 7, respectively. (b) At step 3 and 4, CPU 3 and CPU 5 send the information of boundary atoms to CPU 4, together with the information of atoms at the corners label as Blue. When the MTS algorithm is implemented, the boundary layers can also be further divided into several layers (the number of layers depends on the number of folds in the MTS) according to the interaction range for the different potentials in equation (2.22). In the case of the two-fold MTS, the boundary layer is divided into a primary layer and a secondary layer as shown in figure 2.5. For each MD step in the inner loop of equation (2.26), the primary boundary atoms are copied to each processor. Atom information in the secondary layer is copied only when the outside MD step is executed. 39 Figure 10-5 Schematic of primary and secondary boundary shown as red and black solid squares. Red and Black arrow represents the primary and secondary cutoffs. In the spatial decomposition method, atoms of the simulated system are assigned to each cube according to their coordinates. In MD simulations, atoms can experience large forces and might move in between different subsystems. To account for this, the list of resident atoms in each subsystem is updated at each MD step according to the newest coordinates. Migration of atoms follows a similar method as the copying scheme mentioned above. In cases where atoms do not move across the subsystems, an atom decomposition scheme can be implemented by assigning fixed groups of atoms to each processor. Alternatively, in simulations with small numbers of atoms but large force computational costs, parallelization can also be implemented by partitioning the force calculation between atoms while each processor has position and velocity information for all atoms in the system. The above schemes can be implemented on massive computer clusters, such as USC’s HPC (high performance computing) center, where, as of today, the Linux 40 computing source consists of 1,893 dual-processor nodes and 5 x 16 processor, 64GB large memory servers, interconnected with Ethernet and Myrinet. The most popular software used for communication in between processors is the Message Passing Interface (MPI). In the case of the computing source with shared memory architecture, another application programming interface OpenMP can also be used. In a hybrid model of parallel program, MPI and OpenMP can be used together to improve the computing performance. 10.2 Statistical Mechanics in Molecular Dynamics As discussed above, MD simulations can accurately calculate the time trajectory of a system on the atomistic scale using classical mechanics. However, in general the simulated system only contains thousands of atoms (sometimes billions of atoms), which is far less than in experiments. In order to validate the calculations by comparison with experimental results or to provide insights into experimental results, macroscopic observables such as pressure, energy, heat capacities, etc are often needed from the MD simulations. In MD, these macroscopic properties are connected to the microscopic simulations via statistical mechanics, which calculates the macroscopic properties rigorously by mathematical expressions using the distribution and motion of the atoms in the system. In statistical mechanics, the mathematical expression for a macroscopic property is often an ensemble average of a physical quantity taken over a large number of replicas of the system considered simultaneously: A ensemble = d r p N d r r N A r p N , r r N ( ) ρ r p N , r r N ( ) ∫∫ (2.28) 41 Where ( ) N N r p A r r , is the physical quantity, and ( ) N N r p r r , ρ is the probability density of the ensemble: () ( ) ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − = T k r p H Q r p B N N N N r r r r , exp 1 , ρ (2.29) Where T is the temperature, H is the Hamiltonian, k B is the Boltzmann constant and Q is the partition function: ( ) ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − = ∫∫ T K r p H r d p d Q B N N N N r r r r , exp (2.30) The integration in equation (2.28) is ideally over all possible states of a given system, though this is practically impossible to calculate. Instead, the MD method calculates the time average of the desired quantity and gives essentially the same result as the ensemble average. This is formally stated by the Ergodic hypothesis: A ensemble = A time = lim τ →∞ 1 τ A r p N , r r N ( ) dt ∫ (2.31) This hypothesis indicates that all the possible states of a system occur when the system is allowed to evolve for a long enough period of time. 10.2.1 Statistical Ensembles To mimic macroscopic properties observed in experiments, MD simulations of various ensembles are needed. The most commonly used is the microcanonical ensemble (NVE), which is characterized by thermodynamic state with a fixed number of atoms N, a fixed volume V, and a fixed total energy E. Other frequently used ensembles include the canonical (NVT) and Isobaric-Isothermal (NPT) ensembles. 42 The velocity Verlet algorithm discussed in the previous section was developed to preserve the symplectic property of the Hamiltonian given in equation (2.5). Thus, it can be directly applied to the evolution of a system in the NVE ensemble. The equations of motion for this ensemble are those given as equation (2.6). By modifying the Hamiltonian in equation (2.6) to include the external effect of a thermostat, Nose extended the application of the MD method from the NVE ensemble to the NVT ensemble 65 , where the temperature of the system is kept constant instead of the total energy. This formalism was proven by Nose to be able to give the canonical ensemble distribution of positions and momenta of an ergodic system. The new Hamiltonian is given as: ) ln( 2 ) ,... ( 2 2 1 2 2 s T gk Q p r r s m r H req B s N i N i i i + + + = ∑ = r r r φ (2.32) where the first two terms are the same as in the original NVE system, and the third and fourth terms are associated with the kinetic and potential energies of the thermostat. The terms s and p, are the coordinate and momentum of the thermostat, respectively. Q is a fictitious quantity representing the effective “mass” of the thermostat. Treq is the required temperature. g is the number of degrees of freedom in the system. B k is the Boltzmann constant. The equations of motion can be obtained using the newly formed Hamiltonian, and the conserved quantity of the new system is the Hamiltonian. Through the fluctuation of s, the kinetic energies can be exchanged between the system and the thermostat to effectively control the system’s temperature. The oscillation frequency of s is determined by the effective mass Q, which is often chosen according to the intrinsic thermal 43 frequency of the system. A system in which Q is very large approaches the microcanonical ensemble, and thus is not efficient in controlling the temperature. On the other hand, a system in which Q is too small will result in high frequency temperature oscillations and slow kinetic energy exchange. Because the equations of motion in Nose’s formalism use the virtual time scale dt, which is related to the real time scale t d ′ by s dt dt / '= , oscillations of s cause uneven sampling in the virtual time scale. To avoid this, Hoover 37 later extended Nose’s formalism and transformed the equations of motion into the real time scale. This Nose- Hoover formalism is the most widely used NVT ensemble MD method. The Nose-Hoover formalism was further developed by Martyna et al. to overcome its limitation in systems having few degrees of freedoms, where the assumption of ergodicity in the Nose-Hoover formalism fails 59 . This new formalism includes more than one thermostat variable, and thus is named the Nose-Hoover chain formalism. The Hamiltonian of the system with M thermostats is: i M i B req B M i i s N i N i i i Ts k s T Nk Q p r r m p H i ∑ ∑ ∑ = = = + + + + = 2 1 1 2 1 2 3 2 ) ,... ( 2 r r r φ (2.33) and the equations of motion are: i i i m p r v & r = (2.34) 1 1 ) ,... ( Q p p r r r p s i i N i i v v r r & v − ∂ ∂ − = φ (2.35) i s i Q p s i = & (2.36) 44 2 1 2 2 1 1 3 Q p p T Nk m p p s s req B N i i i s − ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ − = ∑ = & (2.37) 1 1 2 1 1 + − + − − ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ − = j s s req B j s s Q p p T k Q p p j j j j & (2.38) req B M s s T k Q p p M M − = − − 1 2 1 & (2.39) where the effective masses are: 2 3 1 s req B s T Nk Q ω = 2 s req B s T k Q i ω = (2.40) where ω s is the intrinsic thermal frequency of the MD system. Similar to the NVT ensemble implementation of MD with constant temperature constraint, Anderson proposed the NPE ensemble MD method 2 . Here the pressure is constrained to be constant and system volume is treated as a dynamic variable. Parrinello and Rahman 69 later developed Anderson’s method to incorporate the shape change of the system by describing the system with a set of scaled coordinatess i = (s ix ,s iy ,s iz ). The real coordinates of an atom i r i = (r ix ,r iy ,r iz ) are related to the scaled coordinates by the H- matrix: r r i = s i1 r h 1 + s i2 r h 2 + s i3 r h 3 = r r H • r s i , (2.41) where r r H = ( r h 1 r h 2 r h 3 ) = h 1x h 2x h 3x h 1y h 2y h 3y h 1z h 2z h 3z ⎛ ⎝ ⎜ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ ⎟ (2.42) 45 and the three lattice vectors ) h , h , h ( 3 2 1 r r r describe the shape of simulation box. The volume of the system at a given time is determined from the determinant of the H- matrix: V = det || r r H || (2.43) By combining the NPE method and NVT method, MD simulations with constrains of both constant temperature and constant pressure can be performed. The Hamiltonian of the new Isobaric-Isothermal ensemble is given as: () V P s T k N W p Q p r r m p H req req B s N i N i i i + + + + + + = ∑ = ) ln( 1 3 2 2 ) ,... ( 2 2 2 1 2 ε φ r r r (2.44) where, p e , and W are the momentum and effective mass of the barostat. P req is the desired pressure. From the Hamiltonian equations, the equations of motion are obtained as: i i i i r W p m p r r v & r ε + = (2.45) i s i i N i i p Q p p W p N d r r r p v v v r r & v 1 1 3 1 ) ,... ( − ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + − ∂ ∂ − = ε φ (2.46) Q p s s = & (2.47) ∑ = + − + = N i reqt B i i s T k N W p m p p 1 2 2 ) 1 3 ( ε & (2.48) W dVp V ε = & (2.49) () ∑ = − + − = N i s i i req in p Q p m p N d P P dV p 1 2 3 ε ε & (2.50) where P in is the internal pressure and defined as: 46 ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − ∂ ∂ ⋅ + = ∑∑ == N i N i N i i N i i i i in V r r dV r r r r m p dV P 11 2 ) ,... ( ) ( ) ,... ( 1 ∂ ∂φ φ r r v r r r r (2.51) W is the effective “mass” of barostat and normally chosen to be: W = (3N + d)k B T ext ω b 2 (2.52) whereω b is the intrinsic fluctuation frequency of the MD box. 10.2.2 Physical Quantities As stated above, MD simulations calculate the time average of a physical quantity as the ensemble average of it by using the ergodic hypothesis. Physical quantities that can be calculated by MD include thermodynamic quantities, structural quantities and dynamics quantities. 2.2.2.1 Thermodynamics Quantities The most calculated thermodynamic quantity is the temperature of the system, which can be obtained from the average of the kinetic energy <K>: B Nk K T 3 2 > < = (2.53) where k B is the Boltzmann constant, and N is the number of atoms in the system. Fluctuation of the kinetic energy is also related to another useful thermodynamic quantity, the specific heat, which in the case of the NVE ensemble is given as: ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − = V B B NVE C Nk T Nk K 2 3 1 2 3 2 2 2 δ (2.54) 47 Equation (2.53) is the well-known equipartition theory, which says that for each degree of freedom, there is an average energy of k B T/2. This equation can also be derived from the virial theorem: T k q H A B k = ∂ ∂ (2.54) From the virial theorem, the stress tensor can also be derived. The form of the stress tensor that is independent of the origin of the coordinate system is: ∑∑ ∑ < + = ij i ij ij i i i i f r v v m V β α β α αβ σ 1 (2.55) where α, β are the coordinate indices, m i and v i are the mass and the velocity of i-th atom, respectively, f ij is the force between i-th and j-th atoms. In many studies, stress variation at atomistic level is of great interest. There are mainly two ways of calculating the stress at the atomistic level 90,97 . The first is the mechanical approach, which calculates the forces and momentum across a virtual plane of interest and applies them as the stress on the plane. The second approach is the thermodynamic method where the virial expression is given by equation (2.55). The results of these two methods are equivalent when averaged over the volume in a homogeneous system. However, for an inhomogeneous system, the mechanical approach has been argued to be the more appropriate choice 97 . And example of stress calculation by the virial approach can be found in our previous paper on wing-crack simulations 55 . Below is an example of the stress calculation by direct method when studying the void coalescence mechanism. In the simulation, two voids with a given separation distance are generated in silica glass. The whole system is 48 then subjected to hydrodynamic tensile stress, which causes the coalescence of the two voids. To study the stress variation in the ligament region during the process of void coalescence, we selected a plane across the centers of both voids, and plotted the variation of the normal stress on the plane with increasing strain in Figure 2.6. Some fluctuations of the stress are to be expected and make it difficult for analysis. But, with the knowledge of the coalescence of two voids at 7.5% strain, we can correlate the stress variation in the void coalescence with the accumulation and disappearance of stress before and after the coalescence, respectively. Tensile Plane stress variation along with strain -15 -10 -5 0 5 10 15 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Strain Stress (GPa) Figure 10-6. Variation of atomistic-level plane stress with respect to the increasing of strain in the ligament region between two voids. 2.2.2.2 Structural Quantities Useful structural information can be extracted from the pair-distribution function and/or bond-angle distributions. The pair-distribution function (PDF) describes the probability of finding an atom with a given separation distance to another atom. It provides both the long-range and short-range atomic order information of the material, which can be validated by experimental PDF analysis on powder diffraction data. The mathematical expression of the PDF for a uniform system is: 49 g αβ ( r r ) = V N α N β δ r r − r r ij () j ∈ β {} i≠ j ∑ i ∈ α {} ∑ (2.56) where V is the volume of the system, N α and N β are the number of atoms of species α and β, <…> represents an ensemble average, and j i ij r r r r r r − = . For an isotropic system, equation (2.56) will be reduced to a radial pair distribution function independent of the orientation: g αβ (r) = V 4πr 2 N α N β δ r − r ij () j ∈ β {} i≠ j ∑ i ∈ α {} ∑ (2.57) This is especially useful in studying the atomic structure of amorphous and liquid phase materials. Since the probability of finding an atom at a large distance from another atom becomes a constant related to the density of the material, the PDF at large r converges to 1 by normalizing by the density. While the PDF provides useful information about how many atoms are packed within a certain distance to another atom, sometimes it’s more important to know how these atoms are packed. The Bond-angle distribution is a useful tool for providing such information through the calculation of the three-particle angular correlations. Essentially, it calculates the histogram of bond angles of a given type of three-particle triplets. The bond angle is defined as the angle between two arms of the triplet. 2.2.2.3 Dynamics Quantities In addition to thermodynamic and structural quantities, MD can also be used to study various dynamic quantities, such as thermal conductivity, the diffusion coefficient, 50 viscosity, etc. Mathematical formulations of these quantities can all be derived from the famous Green-Kubo formula, proposed by Green and Kubo in the 1950’s 24 . This Green- Kubo formula relates a physical quantity ) (t ϕ to its associated transport coefficient g through the auto-correlation function of its time derivative: ∫ ∞ = 0 ) 0 ( ) ( ϕ τ ϕ τ & & d g (2.58) where ) 0 ( ) ( ϕ τ ϕ & & is the time autocorrelation function of ) (t ϕ & . In the case of the velocity autocorrelation function, the Green-Kubo formula gives the diffusion coefficient: ∫ ∞ = 0 ) 0 ( ) ( 1 α α τ τ v v d d D r r (2.59) where α is the atom species, d is the dimension of the system, α τ ) 0 ( ) ( v v r r is the partial velocity autocorrelation function for species α . Fourier transform of the partial velocity autocorrelation function gives the partial vibrational density of states ∫ ∞ = 0 ) cos( ) 0 ( ) ( 2 ) ( ωτ τ τ π ω α α α α v v d T k m N G B r r (2.60) For the stress autocorrelation, the Green-Kubo formula gives the shear viscosity ∫ ∞ = 0 ) 0 ( ) ( 1 αβ αβ σ τ σ τ η d T Vk B (2.61) The definition of the stress tensors are given in equation (2.55). When the autocorrelation function is of the heat current, the Green-Kubo formula gives the thermal conductivity expression as: ∫ ∞ = 0 2 ) 0 ( ) ( 1 j j d T Vk B r r τ τ λ (2.62) where ) (t j r is the heat current at time t. 51 Since all the above equations are derived from the Green-Kubo formula, their implementation is limited to equilibrium MD systems. Among the above mentioned dynamic quantities, the diffusion constant is the easiest to calculate using this method. The others require long time MD simulations. This is due to the long tail decay in correlation functions of stresses and heat current. In practice, there are also non- equilibrium MD methods developed to calculate dynamic properties by mimicking the experimental setup in MD simulations. In the following chapter, we have used the non- equilibrium MD method to calculate the thermal conductivity of silicon carbide in both the crystalline and amorphous phases. 10.3 Potential Forms Used in Molecular Dynamics As seen above, the only essential input for a successful MD simulation is a satisfactory interatomic interaction potential, adjusted properly to describe the macroscopic properties of the material. Over the years, many different potential forms have been developed from the simple pair potential forms to the more complicated forms involving 3-body interactions, and the later complex bond-order forms. A pair potential form still in use today is the Lennard-Jones potential proposed in 1920’s by John Lennard-Jones 51 . It is the first realistic potential successfully implemented for the study of a number of properties of liquid Ar in the 1960’s by Rahman 74 . In contrast, potential form can also be as complex as the ReaxFF 98 , in which the chemical reactions of continuous bond breaking and forming are described with bond-order and covalent interactions. The calibration of parameters for a given potential form involves iterative modifications to fit the experimental data including lattice constants, bulk modulus, 52 melting temperature and cohesive energy, etc. In certain cases, calibrating to quantum calculation results is also desirable. While simple forms of the potential may not describe accurately all chemical reactions, a too complex form will certainly encumber the computational efficiency of the simulation. In the simulations we have carried out, we have used a potential form consisting of two-body and three-body terms, the parameters of which are successfully calibrated to yield the properties of the ceramics SiC and Al 2 O 3 100,101 . This validation is described in detail in section 2.2.1. To study the dynamics of the detonation process of aluminum nanoparticles, the embedded atom method (EAM) potential 88-90 is used for the aluminum metal (see section 2.2.2) and a coupling methodology between the two potentials (shown in section 2.2.3) is used to incorporate the aluminum oxidation process. 10.3.1 An Effective Interatomic Interaction Potential for Ceramics The semi-empirical potential form used in our simulations for SiC and Al 2 O 3 consists of two terms, namely the two-body term and the three-body term: ) , ( ) ( ) 3 ( ) 2 ( ik N j i N k j i ij jik ij ij r r V r V V r r ∑∑ << < + = (2.63) The two-body interatomic interactions consider the effects of steric repulsion, Coulomb interaction, charge-dipole interaction and Van der Waals interaction as given below: 6 4 ) 2 ( 4 1 ) ( r w e r D e r Z Z r H r V ij r r ij r r j i ij ij s s ij − − + = − − η (2.64) The three-body interatomic interaction reflects the covalent effects through the bond- bending and bond-stretching interactions: 53 () 0 2 0 2 0 0 0 ) 3 ( , ) cos (cos 1 ) cos (cos exp ) , ( r r r C r r r r B r r V ik ij ijk ijk ijk ik ij ijk ik ij jik ≤ − + − ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − + − = θ θ θ θ ξ ξ r r (2.65) In Eq. (2.64), ij H is the steric repulsion strength, i Z is the effective charge in units of the electronic charge, i α represents the charge-dipole strength and ij w is the strength of van der Waals attractions. ij η is the exponent of the steric repulsion, and ls r and s r λ are the screening lengths for the Coulomb and charge-dipole interactions, respectively. ij r is the distance between any two atoms in the system. In equation (2.65), the three-body interaction, ijk B is the strength of the three-body interaction, ijk θ is the angle formed by ij r r and ik r r , and o r is the cutoff for the three-body interaction. 2.3.1.1 Potential Validation for SiC For MD simulations of SiC, at least four empirical interatomic potentials have been proposed: the Pearson potential (PP) 70 , in a form similar to the Stillinger-Weber model; a modified embedded atom method (MEM) by Baskes 9 ; the Tersoff potential (TP) 53 , which is a many-body interatomic potential model including explicit bond-orders; and the potential used here 100 , which consists of two-body and three body terms. Parameters in our potential have been chosen to fit the experimental data on bulk modulus, melting temperature, etc. In addition, the calculated melting temperature and structural transition pressure agree well with experiment values. This potential has also been used to calculate the indentation hardness values of a-SiC 91 , and it predicts a zincblende-to-rocksalt structural transition mechanism 85 , which is later confirmed by first-principles electronic structure calculations. Table 2.1 compares the equilibrium properties of 3C-SiC calculated using the four potentials, and table 2.2 lists all the parameters for our SiC 54 potential. The most up-to-date parameters for our SiC potential can be found in the reference 100 . Table 10.1 Comparison of physical properties of SiC by different potential models Properties MD TP PP MEAM Expt. Result Lattice Constant (Å) 4.36 4.32 4.19 4.2 4.36 Elastic Constants (GPa) C 11 C 12 C 44 390 144 179 436 120 255 1095 937 606 205 390 142 150-256 Bulk Modulus B(GPa) B=(C 11 +2C 12 )/3 225 225 990 211 225 Cohesive Energy (eV) -6.34 -6.18 -7.71 -6.4 -6.34 Table 10.2 Parameters in two- and three-body potential terms used in the MD simulation of SiC. Two-Body Si C Z i (e) 1.26 -1.26 Si-Si Si-C C-C η ij 7 9 7 H ij (eV•Å η ) 474.3120 410.8439 155.4348 D ij (e 2 •Å 3 ) 1.1907 0.59535 0 W ij (eV•Å 6 ) 0 43.4129 0 ls r =5.0Å s r λ =3.75Å r c =6.0 Å Three-Body B jik (eV) θ jik (º) C jik γ (Å) r 0 (Å) Si-C-Si 16.5106 109.47 20 1.0 2.90 C-Si-C 16.5106 109.47 20 1.0 2.90 55 2.3.1.2 Potential Validation for Al 2 O 3 Similar to the simulations of SiC, a number of simulations for Al 2 O 3 have been conducted with different potential forms. Most of the theoretical studies have focused on the structural and dynamic properties of Al 2 O 3 in crystalline phase 33,35,58 . There are also a few potential forms that describe the structural properties of the liquid phase 35,79 . The potential we have used in this study has been validated for properties of Al 2 O 3 in both crystalline and liquid phases. The validation against experimental results is shown in Table 2.3 and the parameters are listed in Table 2.4. Table 10.3 Comparison of physical properties for a-Al 2 O 3 between MD model and experiments MD Experiments E/N (eV/atom) -6.35 -6.35 B (GPa) 253 255 C 11 (GPa) 523 498 C 12 (GPa) 147 163 C 13 (GPa) 129 117 C 14 (GPa) 7.5 -23 C 33 (GPa) 427 502 C 44 (GPa) 135 147 C 66 (GPa) 174 167 * * Calculated from the relation, C 66 = (C 11 -C 12 )/2. The experimental data for the elastic constants are from Ref. {Gieske, 1968 #1} Table 10.4 Parameters in two- and three-body potential terms used in MD simulation of Al 2 O 3 Two-Body Al O Z i (e) 1.5237 -1.0158 Al-Al Al-O O-O η ij 7 9 7 H ij (eV•Å η ) 12.4236 244.3038 553.3917 D ij (e 2 •Å 3 ) 0 3.4374 1.52775 W ij (eV•Å 6 ) 0 0 63.4894 56 Table 2.4: Continued λ=5.0 Å ξ=3.75 Å r c =6.0 Å Three-Body B jik (eV) θ jik (º) C jik γ (Å) r 0 (Å) Al-O-Al 8.1149 109.47 10 1.0 2.90 O-Al-O 12.4844 90.0 10 1.0 2.90 In a recent study of the structural and dynamical properties of both crystalline and amorphous phases, this potential was further validated by comparing the simulated neutron static structure factor of both alpha and liquid Al2O3 with the experimental results. Dynamic properties such as the vibrational density of states of alpha alumina are also consistent with experiments (see the detailed paper 101 ). In our current research studying the detonation process of aluminum nanoparticles, the transition of alumina to the liquid phase and the specific heat of alumina at high temperatures are both important. To illustrate the validation of our potential on these issues, plots of the melting of alpha alumina and the specific heats of alpha and amorphous alumina obtained from reference 101 are shown again in Fig. 2.7 and 2.8, respectively. Figure 10-7. Energy per atom (top) as a function of temperature. The vertical dashed line is the experimental reported melting temperature. 57 Figure 10-8. Specific heat for alpha and amorphous alumina at two densities. 10.3.2 EAM Potential for Metals In MD simulations of metals, there have been many potential forms developed, the most popular one being the embedded-atom method (EAM). The EAM potential was first proposed by Daw and Baskes [Phys. Rev. B29, 6443 (1984)] to study the properties of realistic metal systems based on density-functional theory. EAM was developed further to study the shear behavior of covalent material silicon by Baskes in 1987 and metallic alloys by Foiles et al. in 1986. Since then, the EAM has been continuously developed to study various properties over a wide selection of metallic materials, including both f.c.c, b.c.c and h.p.c structured pure metals and metallic alloys. In 1987, Voter and Chen successfully applied EAM to accurately describe the Ni3Al alloy system as well as systems of diatomic Ni 2 , diatomic Al 2 , f.c.c Ni and f.c.c Al. In our research into the detonation mechanism of aluminum nano-particles, we have adopted the Voter-Chen’s EAM formalism and coupled it with Vashishta’s potential. This will describe the not fully oxidized Al-Al interaction in the presence of oxygen. The next two subsections will briefly introduce the EAM potential and its validation for aluminum metal. The interpolation between the two potentials is discussed in following section in detail. 58 2.3.2.1 Formulation of Voter-Chen’s EAM The underpinnings of EAM are similar to the basis of density functional theory (DFT), which states that the energy of a system of atoms can be expressed exactly by a function of its electronic density. In the original EAM proposed by Daw and Baskes 20 , the energy of a metallic system is calculated by the following expression: ∑ = n i i tot E E (2.66) where n is the number of atoms in the system and Ei is the energy of atom i: ) ( ) ( 2 1 i i j ij i F r E ρ φ + = ∑ (2.67) The first term in the above equation describes the electrostatic interaction between atom i and its neighboring atoms. The second term describes the attractive interaction which models placing (embedding) a positively charged atom in the electron density due to the free valence sea of electrons created by the host system of atoms. The embedding function ) ( i i F ρ is a function of the superimposed densities i ρ , due to the charge density distribution of neighboring atoms: ∑ = j ij i r ) ( ρ ρ (2.68) where ) ( ij r ρ is the pair-wise electronic density as a function of the distance ij r between atom i and atom j, but without angular dependency. Originally, the electrostatic interaction was composed of a pure repulsive columbic term and the electron densities taken from the Hartree-Fock calculations. In different applications, however, the way the electron densities are determined and the forms of the pair potentials vary. In many cases the functional ) ( i i F ρ is nonlinear and positively 59 curved ( 0 2 2 > ρ d F d ). The physical significance of the nonlinearity is that it provides a many-body contribution to the energy. The positive curvature ensures the weakening of successive bonds, which is consistent with the traditional chemical bonding concepts (making more bonds or shorter bonds increases the bond order of an atom). In the case of Voter-Chen’s EAM expression, the pair-wise electronic interaction is taken to be a Morse potential: { } M M M M D R r D r − − − = 2 )] ( exp[ 1 ) ( α φ (2.69) where D M and R M define the depth and position of the minimum, respectively. α M is a measure of the curvature at the minimum. The density function is a hydrogenic 4s orbital with a relative normalization factor added to ensure the monotonicity: ) 2 ( ) ( 2 9 6 r r e e r r β β ρ − − + = (2.70) where β is an adjustable parameter. In order to implement the EAM in MD to scale as O(N), both ) (r φ and ) (r ρ are truncated at cutoff r cut using: cut r r m cut cut cut dr dh r r m r r h r h r h / = ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − + − = 1 ) ( ) ( ) ( (2.71) where ) (r h can be ) (r φ or ) (r ρ , and m=20. With these definitions, the only unknown in the EAM form is the function ) ( i i F ρ , which is numerically calculated by matching the universal energy form of Rose et al 77 . ∗ − ∗ ∗ + ⋅ − = a coh e a E a E ) 1 ( ) ( (2.72) 60 whre E coh is the cohesive energy per atom of the crystal, and ∗ a is the reduced lattice constant defined by the equilibrium lattice constant 0 a , the bulk modulus B, the equilibrium atomic volume Ω and E coh . Ω − = ∗ B E a a a coh 9 ) 1 ( 0 (2.73) Similar models developed for metallic systems include the Finnis-Sinclair’s model 25 that is based on the second moment approximation to Tight Binding theory. This model explicitly states that, the band width and thus the energy of a system is proportional to the square root of the electron density, instead of numerically constructing the embedding function as in EAM. This method has the advantage of being easily extended for a more accurate potential description by incorporating higher moments. It becomes increasingly difficult, however, when constructing binary and ternary alloys. The glue model developed by Ercolessi et al. uses same forms as EAM but only uses a first-nearest- neighbor cutoff distance 23 . Although this model can successfully describe various Au surface reconstructions, its application is limited because high order neighbor interactions are needed in cases involving more complex crystal structures. 2.3.2.2 Validation of Voter-Chen’s EAM for Aluminum metal As has been shown by Voter and Chen in their papers 102,103 , the Voter-Chen formalism for aluminum has been validated by comparing various properties of simulated pure aluminum metal with experimental values. The compared properties are given in Voter-Chen’s paper and listed in table 2.5 for convenience. The advantage of the EAM’s fitting procedure is that it exactly matches the lattice constant, cohesive energy and bulk modulus to experimental values. Matching of the other properties, including elastic 61 constants, vacancy formation energy ΔE f 1V , bond length R e and bond energy D e of the diatomic molecule, is based on minimizing the root-mean-square deviation χ rms between the experimental and calculated values. Table 10.5 Comparison of physical properties for Al between EAM model and experiments EAM Experiments Ecoh(eV) 3.36 a 0 ( Ǻ) 4.05 B (10 12 erg/cm 3 ) 0.79 C 11 (10 12 erg/cm 3 ) 1.07 1.14 C 12 (10 12 erg/cm 3 ) 0.652 0.619 C 44 (10 12 erg/cm 3 ) 0.322 0.316 ΔE f 1V (ev) 0.73 0.75 D e (ev) 1.54 1.60 R e ( Ǻ) 2.45 2.47 χ (rms%) 3.85 Validations of the Voter-Chen EAM also include the successful study of interfacial structures, dislocation structures, grain boundary and free surface properties in metals 103,104 . A modified version of EAM incorporating angular dependences into the electron densities has also been applied successfully in the study of binary systems, such as nickel/silicon and aluminum/sapphire systems 3 . In more recent times, Boyukata and Guvenc 12 have successfully studied the energetics, melting and isomerization of aluminum microclusters Al n ( 13 2 ≤ ≤ n ) by using the Voter-Chen EAM potential in MD simulations implemented using the predictor-corrector algorithm. In our research of a typical aluminum nano-particle, where there are about two million Al atoms contained within the metallic core, we focus on the melting and the oxidation of the aluminum as the driving force behind the detonation of the aluminum nano-particle. As such, we have further tested the equation of state of aluminum by this model for both NVE and NPT 62 ensembles. The results are shown below in figure 2.9. From figure 2.9 (a), we can see the melting temperature of aluminum at zero-pressure is very close to the experimental value of 933 K. In confined conditions (figure 2.9 (b)), the melting temperature is elevated to the higher value of 1713 K. Figure 10-9 Equations of state of aluminum metal by EAM under the conditions of both NPT ensemble (a) and NVE ensemble (b). The red vertical dash line in (a) is used to indicate where the experimental melting temperature (933 K) is located. 10.3.3 Coupling of Two Different Potential Forms The oxidized Aluminum Nano-Particle (OANP) has been studied by large scale MD simulations, which requires building a stable interface between the pure aluminum core and the alumina shell in the initial state of the system. However, a simple average of Al- Al interactions from Vashishta’s potential 101 and Voter’s potential 103 can’t stabilize the interface due to the large difference in the potential stiffness of the two potentials. Consequently, a many-body interpolation scheme is proposed and validated by examining the stability of the interface and the energy conservation. In the OANP system, the interaction between O-Al takes the form of Vashishta’s potential, in consideration of the strong oxidizing nature of aluminum. The interaction 63 between O-O also takes the form of Vashishta’s potential, which is a pure repulsive force. This is again adopted in considering that oxygen molecules dissociate easily to oxidize the not-fully-oxidized aluminum atoms. This will certainly cause the amount of energy released different during the process of the aluminum nano-particle’s detonation. However, since the main purpose of our simulation is to check the mechanism of the detonation, this effect should be minimal. As Voter’s aluminum potential is a many-body potential, the interpolation for the Al-Al interaction is done for each Al atom, such that the total energy of this Al atom due to neighboring Al atoms is the weighted average of the Vashishta’s potential and Voter’s potential. The weight factor is based on the oxidation degree of the Al atom, which is modeled by a function that is smooth in the second order of the derivatives. 2.3.3.1 Interpolation Scheme Mathematically, the total potential energy of the system can be written as the sum of the potential energy of all aluminum atoms and that of all oxygen atoms: ∑ ∑ ∑ + = = ) ( ) ( O i i Al i i i i tot E E E E (2.74) The potential energy of each Al atom i is a weighted average between Vashishta’s potential and the EAM: ) 1 ( ) ( ) ( 2 1 2 1 ) 1 ( , , , , 2 , , f r F r f V f E f E E i j j i i j j i i j j i b p i v i p i − ⋅ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + + ⋅ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = − ⋅ + ⋅ = ∑ ∑ ∑ ≠ ≠ ≠ − ρ φ (2.75) where f is the weighting factor as a function of the oxidation degree o i n of atom i: 64 () () ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ >= < < + ⋅ − <= = 1 1 1 0 1 ) 5 . 0 ( sin 2 1 0 0 ) ( o i o i o i o i o i n n n n n f π , (2.76) A plot of the weighting factor is shown below Figure 10-10. Plot of the weighting factor as a function of the oxidation degree o i n The weight factor is constructed such that it is continuous on the order of second derivative. The derivative of the weight factor is given as: () ⎪ ⎩ ⎪ ⎨ ⎧ < < ⋅ ⋅ ⋅ − = else n r d dn n r d n df o i i o i o i i o i 0 1 0 2 ) 5 . 0 ( cos ) ( r r π π (2.77) The oxidation degree o i n of atom i is: ) ( ik k o i r n ∑ Θ = (2.78) where ) ( ik r Θ is a continuous function that counts the degree of oxidation according to the distance of the considered aluminum atom to a neighboring oxygen atom: 65 ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ + ≥ + < < − + − + + − − − ≤ = Θ α α α α α α α α α α α α α α π π D R r D R r D R D D R r D D R r D R r r , 0 , 2 ) ) ( sin( 2 1 , 1 ) ( (2.79) The plot of ) ( ik r Θ as the function of the distance between Al and O atoms is: Figure 10-11. Plot of ) ( ik r Θ as the function of distance between a pair of Al and O atoms The derivative of the oxidation degree factor is given as: () ∑ ∑ ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ + ≥ + < < − + − + − − ≤ = Θ = k k i k i i o i D R r D R r D R D D R r D D R r r d r d r d n d α α α α α α α α α α α α π , 0 , ) ) ( cos( 1 2 1 0 ) ( ) ( , v v (2.80) The forces on atoms involved in the energy equation (2.75) are calculated as: 1) Force on the aluminum atom i: () i j i b p i j j i i j i j i i Al j i j i i j i b p i i p i v i i v i i p i i i r d df V F r f r d r d d dF f r d r d f r d dV r d df E E f r d dE f r d dE E r d d F r r r r r r r r ⋅ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − + + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − ⋅ ⋅ + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − ⋅ + ⋅ − = ⋅ − + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − ⋅ + ⋅ − = − = − ≠ ≠ ≠ − ∑ ∑ ∑ , , 2 , , ) ( , , , 2 , , , , 2 1 ) ( ) ( 2 1 ) 1 ( ) ( ) 1 ( ) ( 2 1 ) 1 ( ρ φ ρ ρ φ (2.81) 66 2) Force on each neighboring Al-atom j: ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − ⋅ ⋅ + − ⋅ + ⋅ − = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − ⋅ + ⋅ − = − = − ) 1 ( ) ( ) 1 ( ) ( 2 1 2 1 ) 1 ( , , , , 2 , , f r d r d d dF f r d r d f r d dV f r d dE f r d dE E r d d F j j i j j i j j i b p j i v j i p i j j r r r r r r ρ ρ φ (2.82) 3) Force on each neighboring O-atom k: () k j i b p i j j i k i p i v i k k r d df V F r r d df E E E r d d F r r r ⋅ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − + = ⋅ − = − = − ≠ ∑ , , 2 , , , 2 1 ) ( ) ( 2 1 ρ φ (2.83) 2.3.3.2 Validation for the coupled force field To validate the interpolation scheme, we first checked the energy conservation of the nano-particle system with this new potential. A typical energy fluctuation plot for various time step sizes is shown below in figure 2.12. From this figure (a), we can see that the total energy of the system conserves very well even for time step size of 2fs. Fluctuation of the total energy reduces accordingly with respect to the reducing of time step size. From figure 2.12 (b), the linear proportionality of energy fluctuation to the square of the time step is shown. In the test simulation, the time step of 1fs was used. 67 Figure 10-12. (a) Energy fluctuations of the aluminum nano-particle system with various time steps tested. (b) Plot of the linear proportionality between the standard deviation of energy with respect to the square of the time step size. Another check of this interpolated potential is if this scheme can make a stable interface between the aluminum core and the alumina shell once they are put together in the initial configuration of the aluminum nano-particle. This is because the strength of 68 oxidation can be affected by the choice of parameters in the interpolation scheme. The degree of oxidation of the Al atom can be varied using different values of the two parameters α R and α D in function ) ( ik r Θ . Specifically, for those aluminum atoms at the interface of aluminum core and alumina shell in the initial configuration, most of them are partially oxidized. The choice of α R and α D will affect how much the Al atoms from the core will be oxidized, as well as the structure of the interface. In the following graphs, we compare the interface structure resulting from changing the value of α R from 3.0 Ǻ to 3.5 Ǻ, while fixing the value of α D at 0.5 Ǻ. Comparison between these two choices shows that when the value of α R is increased, the aluminum atoms become more easily oxidized and the repulsive Al-Al force from Vashishta’s potential results in a separation at the interface. Figure 10-13. Stability test for the interface between Al core and Al 2 O 3 shell when changing the parameter α R from 3.0 Ǻ (a) to 3.5 Ǻ (b). 69 We further validated this interpolated potential by running MD simulations for several Al X O Y molecules. These molecules are quickly quenched by the steepest descent method to the energy minimum configuration from their initial state, which is similar to, but not exactly the same as, the quantum mechanical result for the stable configuration, listed in the last column of the table below. These molecular structures at the energy minimum found via MD are listed in the middle column to compare to the quantum calculations. Table 10.6 Comparison of bond length ( Ǻ) and bond angles (deg) in Al X O Y molecules MD CBS-QB3* AlO AlO 2 Al 2 O Al 2 O Al 2 O Same as above Al 2 O 2 Al 2 O 2 O 1.630 Al O Al O 1.704 O Al O 1.654 O 1.682 Al O Al Al 2.446 3.238 O Al Al 2.819 1.615 O Al Al 2.607 1.742 97 Al O Al 1.729 Al O Al 1.716 Al O Al 1.722 O 1.756 1.699 Al O Al 1.724 O 1.676 1.604 Al Al O 3.239 O 2.448 Al Al O 1.602 O 2.550 70 Table 2.6: Continued Al 2 O 2 (planar) Al 2 O 3 Al 2 O 3 (planar) Al 2 O 3 (planar) Same as above Al 2 O 4 (planar) Al 2 O 4 O Al Al 2.57 1.76 93 43 O O Al Al 2.439 1.757 88 46 O O Al O 1.699 Al 1.748 O O Al O 1.600 Al 1.678 O O Al Al 1.739 1.735 87 92 O O 1.761 133 O Al Al 1.72 1.81 94 88 O O 1.75 138 O Al 1.667 1.700 58 O O 1.726 151 Al 87 O Al Al 1.736 O O 1.739 133 O 93 94 O Al Al 1.72 O O 1.79 137 O 86 92 87 O Al Al 1.734 O O 1.760 43 O 94 1.846 O Al Al 1.74 O O 1.78 O 2.53 1.86 71 Table 2.6: Continued Al 3 O 2 (planar) Al 3 O 2 (planar) Al 3 O 3 (planar) Al 3 O 3 (planar) * CBS-QB3 involves density functional (B3LYP) geometry optimization and ab initiao “complete basis set” extrapolation 72 . From the above table, it is clear that this interpolated potential gives the correct structures for a number of molecules with varying Al and O composition. Among the above seventeen tested configurations, there are only three of which show different structures from the QM calculations. The other fourteen configurations are very similar to the quantum results except for a small variation in the bond length and bond angles. It is well known that the final stable structure of alumina is Al 2 O 3 . So, these checked structures are the intermediate molecular structures that could be present during the 124 O O 1.72 Al Al 1.75 170 Al 122 O O 1.707 Al Al 1.707 180 Al Al O O 1.77 Al Al 1.74 136 98.6 1.89 Al O O 1.762 Al Al 1.742 130 87 1.927 Al 88 O Al Al 1.679 O O 1.755 134 92 1.746 92 1.712 Al 93 O Al Al 1.76 O O 1.75 137 88 1.79 86 1.71 O Al O 1.801 Al O Al 1.694 1.870 1.854 93 47 42 O Al O Al O Al 1.75 105.3 134.5 72 process of detonation. The confirmed minimum energy structures for these intermediate molecules thus gives part of the validation of this interpolation scheme applied to the procedure of detonation of aluminum nano-particles. In the above MD configurations for Al X O Y molecules, the parameters used for α R and α D are equal to 3.0 Ǻ and 0.5 Ǻ, respectively. As also noted above, the change of these two parameters’ values, especially the value of α R , can affect the structures. An example is given for the third molecule Al 2 O listed in the above table 2.6. By changing of the value of α R , we show the large change of bond lengths for this molecule at the stable configuration (see table 2.7). These results show that the precision of describing the oxidation dynamics of aluminum by this interpolation scheme can still be improved by adjusting these variables. Table 10.7 Comparison of bond length ( Ǻ) in Al 2 O molecules in MD simulations with different value of α R Initial setup Final stable structure 0 . 3 = α R Ǻ 0 . 3 = α R Ǻ 0 . 2 = α R Ǻ 5 . 1 = α R Ǻ CBS-QB3 72 O Al Al 3.0 3.0 O Al Al 2.446 3.238 O Al Al 3.0 1.8 O Al Al 3.0 2.1 O Al Al 2.424 2.152 O Al Al No bond 1.682 O Al Al 2.819 1.615 O Al Al 3.0 1.8 O Al Al 2.375 1.72 73 Chapter 3 Thermal Conductivity of Ceramic SiC When heat accumulats at certain region in ceramics, it will dissipates gradually to neighboring cold regions along the direction of temperature gradient. To study the slow heat transportation mechanism in ceramics, we have performed extensive simulations of thermal conductivities for both crystalline and amorphous SiC. In the following sections, the simulation method used for calculating the thermal conductivities and the corresponding results are presented. 12.1 Thermal Conductivity Simulation Methodologies Thermal conductivity is defined in Fourier’s law of heat conduction as the proportionality constant which relates the heat current to the temperature gradient. ∑ ∂ ∂ − = ν ν μν μ κ x T J , μ , v = 1, 2, 3 (1) where, J is the heat current vector and ν x T ∂ ∂ is the temperature gradient. Usually, μν κ is a 3 by 3 tensor. In the case of a cubic structure crystal, it degenerates to a scalar. There are two methods generally used in molecular dynamics for studying thermal conductivity. One is the Equilibrium Molecular Dynamics (EMD) 81,96 method, which often uses the Green-Kubo formalism. The other one is the Non-Equilibrium Molecular Dynamics (NEMD) method 19,43,62,93 , which mimics the experimental direct-measurement case (thus named the “direct method”). 74 12.1.1 Equilibrium MD Method The Green-Kubo method is a result of the fluctuation-dissipation theorem, which defines transport coefficients in terms of linear response of the system to an external perturbation. Green-Kubo formalism relates the thermal conductivity to the equilibrium current-current autocorrelation function as: ∫ ∞ > < = 0 2 ) 0 ( ) ( 1 dt J t J T Vk k B ν μ μν (2) where, B k is the Boltzmann constant and V is the volume of the system. The angular bracket denotes the ensemble average and J μ (t) is the microscopic heat current at time t, defined as: J μ (t) δε μ = r iμ (ε i − ε i ) i ∑ (3) where μ i r is a coordinate component of atom i at time t, and i ε is the site energy. In the case involving a three-body term in the energy form, the resulting heat current from this derivative can have different forms with respect to different definition of the size energy. But it was tested that the definition of site energy has very little effects on the determination of thermal conductivity. The Green-Kubo method is applied for system that is in equilibrium state. Since there is no external driving force imposed in the system, the linear response assumption naturally holds. An obvious advantage of this method is that it can easily capture the anisotropic effects of thermal conductivity in a single simulation. Also, there is the advantage of having less severe finite-size effects compared to NEMD. For example, in 75 the case of calculating thermal conductivity of Si by P.K.Schelling et al 81,96 , the convergence of k is obtained with only 1728 atoms. Nevertheless, they have also found that this method needs very long time simulation in order to reduce the noise-signal ratio. Because the convergence of current-current autocorrelation function is very slow, there were works previously attempted to reduce this demanding requirement for long time simulations. Methods include either exponentially fitting the current-current autocorrelation function by using the data obtained from the beginning period of the autocorrelation function, or alternatively extrapolating the Fourier transform of the ensemble average of heat current to the zero frequency limit. However, the basic assumption of exponential decay in the current-current autocorrelation function that is used to reduce the simulation time turns out to be unjustified. A double-exponential fit has been recently proposed for solid argon 96 . This, however, will not significantly reduce the simulation time, since to fit the second part of the current-current autocorrelation function, a sufficiently long time run is still needed. 12.1.2 Non-Equilibrium MD Method Historically, the NEMD direct method was used in calculating thermal conductivities of liquids with 1-D and 2-D solid models. Variations of this method include making a rigid wall connected to a heat reservoir 93 , using thermostats in the simulated systems 19 , setting up steady heat flow by exchanging cold and hot atoms’ velocities 62 , and properly scaling velocities of atoms in heat source and sink regions of the system as proposed by P. Jund and R. Jullien 43 . In our simulations, we selected the last algorithm for its merits discussed in the following. 76 Schematic representation of the NEMD direct method is shown as Fig.1, in which the two shaded areas are heat source and sink respectively. In this algorithm, same fixed amount of heat are input and output from the system by scaling the velocities of atoms in heat source and heat sink respectively. The volume of the system and the total number of atoms in the system are not changed. As such, this method is consistent with a microcanonical ensemble system. Also, this algorithm has the advantage of being compatible with the periodic boundary condition, which is apparent in the symmetry of the system. Because it is necessary for a local temperature to be defined using velocities in a local system at equilibrium, or at least in a local system where the center of mass has zero velocity. This rescaling algorithm is also favorable as it meets these requirements (seen in the algorithm below). Figure 12-1 Schematic description of the velocity-scaling method. System is divided into slices of equal size. Velocities of atoms in the two slices positioned at ¼ and ¾ of the system length are scaled to remove and add heat energy into the system respectively. The algorithm of velocity rescaling is straightforward as seen in the work by P. Jund and R. Jullien 43 , and are shown again for reference purpose as follows: New velocities of atoms in the two slices are given by: ) ( ' G i G i V V V V − + = α (4) ¾ L Inputting Heat Taking out Heat ¼ L 77 where, i V is the velocity of each atom in the two slices before being scaled, G V is the velocity of the center of mass of the two slices, and α is the scaling factor, defined as: R C E ε α Δ ± = 1 (5) where ε Δ is the inputted or outputted heat, and R C E is the kinetic energy of each slice relative to its center of mass: ∑ ∑ − = i G i i i i R C v m v m E 2 2 2 1 2 1 (6) The temperature gradient for a fixed-length system is determined mainly by the value of ε Δ and the rate of scaling velocity. A typical simulation using this method for a system initially at equilibrium includes a beginning run of making steady heat flow state, where the total system’s average temperature is fixed, and a subsequent run of collecting temperature data with the temperature constraints removed. A linear fit is applied to the collected data of the temperature profile along the heat flow direction (y-axis in Cartesian coordinates) to estimate the temperature gradient. The heat flux in the system is calculated as: t A J y Δ ⋅ Δ − = ε (7) where A is the area of the cross section of the system, and t Δ is the time interval between each energy input to the system. The thermal conductivity is then determined from the Fourier’s law of heat conduction. This NEMD direct method, when compared with the EMD method, is believed to have advantages including shorter required simulation time, less statistical error, and 78 providing direct information on transport phenomena. Especially when dealing with heat conduction across the interface between different materials, this method is preferable. However, this method does have the disadvantage that it can only provide one component of the thermal conductivity for each simulation process, which is inconvenient when calculating anisotropic thermal conductivity. Also, this method has severe finite-size effects on the calculated thermal conductivity. This finite-size effect can be corrected by extrapolating from the inverse thermal conductivities of finite systems of different sizes, as has been shown by many other researchers’ simulations 43,81 . This extrapolation is based on the estimation of the thermal conductivity from elementary kinetic theory, and the estimation that the inverse of the phonon’s effective mean free path (m.f.p.) is proportional to the inverse of the system’s size. Thermal conductivity deduced from elementary theory is given as follows: l v Cv * * 3 1 = κ (8) where Cv is the specific heat, v is the average phonon group velocity, and l is the m.f.p. of the phonon. By considering only the contribution from acoustic modes, the estimated relationship for the size of the systems and the thermal conductivity of silicon was obtained by P.K.Schelling 81,82 , etc. as follows: ) 1 ( 2 1 Z B L b l nv k + = ∞ κ (9) where Z L is the size of systems, n is the number density of the system, v is the phonon group velocity, ∞ l is the phonon m.f.p. of infinite volume, and b is a parameter that was taken as 4 in Schelling’s work based on the assumption of ballistic movement of phonons 79 in their system. This equation is a crude estimate based on elementary kinetic theory for the thermal conductivity calculation and the classical Dulong and Peti model for the calculation of specific heat. Nevertheless, it has been shown that this equation was successful in estimating the bulk thermal conductivity of Si. Also, the slope and intercept can be used to estimate the phonon group velocity and m.f.p respectively. This extrapolation solves the problem of finite-size effect, but makes the total simulation time longer and induces more error in extrapolating the bulk thermal conductivity. This is nontrivial, especially for systems having large thermal conductivities. From the above equation, Schelling 81 also estimated that to obtain a value of κ within 10% of the correct bulk value, the size of the system has to be at least 40 times larger than the bulk m.f.p. Nevertheless, considering the isotropic characteristics of this SiC polytype and our capability of running large-scale simulations on very large system sizes with parallel multi-processors, we have chosen to use this method in calculating the thermal conductivity of 3C-SiC. 12.2 Thermal Properties of 3C-SiC Elementary kinetic theory has shown that the thermal conductivity is directly connected to the specific heat and the phonon group velocity. The phonon group velocity can be calculated directly from the lattice constants, the validation of which has been shown in Chapter 2, Table I. We will further validate the application of our potential to this study by calculating the specific heat from the density of states (d.o.s.), obtained from our simulation. 80 12.2.1 Specific Heat of 3C-SiC The density of states of 3C-SiC can be easily calculated by applying the Fourier’s transform to the velocity autocorrelation function obtained during the relaxation of an equilibrium system. In our simulation, an equilibrium 3C-SiC system of 13702 atoms was used to calculate the velocity autocorrelation function and the density of states of the system. The results are shown below in Figure 2 (a). From Figure 2 (a), we can see that the range of optical modes and acoustic modes fall in the range of less than 20THz and 22~29THZ, which is in good agreement with experiment. In Figure 2 (b), we also plot the calculated temperature dependence of the heat capacity, which can be readily calculated from the following equation once the density of states is obtained. ∫ − = ω ω ω ω ω d g T k T k T k k Cv B B B B ) ( ] 1 ) [exp( ) exp( ) ( 2 2 h h h (10) The calculated specific heat shows excellent agreement with the experimental values for 6H-SiC 61 . 6H-SiC has a different structure than 3C-SiC, but their specific heats have been shown to be close. Also, the calculated specific heat at the high temperature limit is shown to be consistent with the classical limit 3NK B , which is plotted as a straight line in Figure 2(b). The confirmation of the correct range of acoustic and optical modes, the consistence of our simulated constant volume specific heat with the experimental result, and the trend to the classical high temperature limit further validates our potential model in calculating the thermal properties of 3C-SiC. 81 Figure 12-2 (a) The calculated vibrational phonon density of states of 3C-SiC by molecular dynamics using our potential. (b) Temperature dependence of specific heat (constant volume) of 3C-SiC as a result of integration using the d.o.s from (a). The classical limit of the specific heat (1.2423 J/gK) is plotted as a straight line in (b) 12.2.2 Thermal Conductivities of 3C-SiC Thermal conductivities of 3C-SiC at various temperatures are calculated using the NEMD method with the P. Jund and R. Jullien algorithm 43 , as explained in the previous section. Since there is an artificial heat input and output in two slices of the system, the scattering of phonons in these two regions will be affected. This will result in non-linear 82 regions near the two slices. A typical temperature profile of the system at temperature 700K is given in Figure 3. In the previous work of other authors, two strategies have been used for deciding the temperature gradient either including the nonlinear region or not 43,81 . We will first show the results which include the non-linear region, and will then show that there is actually little difference in the results obtained from these two strategies. Figure 12-3 Temperature profile along y-axis for a 1,048,576 atoms 3C-SiC system at 700K There are two main factors, which affect the final results of the simulation. One is the simulation time, consisting of the time used for establishing the required steady heat flow in the system and the time used for averaging the temperature of each slice of the system. The second factor is the well-known system size effect discussed earlier. In considering the time factor effect on the simulation result, we have performed very long time simulations for all the systems at different temperatures. A typical simulation includes 50 ps relaxing the equilibrium system at each given temperature, 250 83 ps making steady heat flux status, and finally 2.0 ns of collecting temperature data. In the last 2.0 ns, the temperature of each slice of the system is calculated, and averaged over a 0.5 ns interval. This long run is necessary as the temperature gradient profile can not be stabilized in a shorter time 56 . As shown in Figure 4. , the temperature gradient calculated from the average of the four 0.5 ns periods can be different by as much as 50%. In all our simulations, we have used the temperature averaged over the last 0.5 ns to ensure the steady heat flow condition is satisfied. Figure 12-4 System temperature profiles and their linear fits for four consecutive 0.5 ns periods. The differences among the temperature profiles can be seen from the slopes of the linear fits. The temperature gradient of the first 0.5 ns period differ by about 47% from the first 0.5 ns period. 84 With regard to the size effect, researchers have shown that the bulk thermal conductivity of dielectric materials can be correctly extrapolated from small size systems 53 . In our case, because of the large thermal conductivity and relatively long phonon mean free path of 3C-SiC, it is desirable to perform this NEMD method on relatively large size systems, and extrapolate the results. Four sizes have been used for the extrapolation. All the four systems used have the same cross section area: 49.3Å by 49.3Å. The length of the four systems varies from 2511.1Å to 5580.3Å. The systems used are so large (the maximum number of atoms in a system is about 1.3 million) that parallel computing becomes a necessity. The largest system has used 160 CPUs for our simulations. For a given temperature, the thermal conductivities with different size systems are calculated by using equation (1). Finally the bulk thermal conductivity is extrapolated from the calculated finite-size values. A representative plot of the extrapolation of bulk thermal conductivity of 3C-SiC at 700K is given in Figure 5. Figure 12-5 extrapolation of thermal conductivities from four different size systems at 700K, the inverse of the intercept gives the bulk thermal conductivity as 59.88(Wm-1K- 1). 85 The thermal conductivity dependence on temperature is shown in Figure 6. These results include the nonlinear region in calculating the temperature gradients. Figure 12-6 The temperature dependence of thermal conductivity, extrapolated from four different sizes of systems. Also shown are experiment and simulation results from another group 92 (data is exerted from ref. by J.Li et. al 53 ) The results from Ju Li’ EMD simulation 53 and experiments are also shown in Figure 6. Compared to their results, we obtain smaller values along the temperature line. However, the decreasing trend of thermal conductivity with increasing temperature from our simulation is the same as that obtained by the other two means. The reason for the resulting smaller thermal conductivities might be attributed to the larger anharmonicity from our potential as compared to the values from experiments. The larger anharmonicity also produces larger linear thermal expansion coefficients for our system, compared to the experimental results. The thermal expansion coefficient has been found to be around 4.5 (10 -6 K -1 ) experimentally 54 . To calculate the thermal expansion of our system, we first determine the zero pressure volumes of our system at various temperatures. The results are shown below in Table 1. Using this data, the calculated thermal expansion coefficient of our system is 1.097 (10 -5 K -1 ), which is about 86 two or three times larger than the experimental values. It has been shown that when a perfect crystal is at temperatures greater than the Debye temperature, the thermal conductivity is inversely proportional to the square of Grüneisen parameter. This parameter is proportional to the thermal expansion coefficient, as is shown in the following equation. Cv B and α γ γ κ 3 , ) 5 . 0 ( 1 2 = + ∝ (11) where γ is the Grüneisen parameter, α is the thermal expansion coefficient, and B is the bulk modulus. Thus, thermal conductivity decreases with increasing thermal expansion coefficient. By using the thermal expansion coefficient calculated here and the specific heat and bulk modulus calculated previously, the Grüneisen parameter is found to be 2.3. Replacing the thermal expansion value with the experimental one gives the Grüneisen parameter of 0.95, which reduces the thermal conductivity by a factor of 3. Table 12.1 The zero pressure volumes of the 3C-SiC systems at various temperatures. Temperature (K) Volume (Å 3 ) 1 294283.4 500 299003.2 1000 303793.6 1500 309227 2000 314178.1 As stated before, we have also calculated the thermal conductivity by using a strategy that excludes the regions in the vicinity of the heat source and heat sink. This has been carried out using the middle 80% of the temperature profile region and the middle 60% of the temperature profile. The comparison between these results is shown in Table 2. It is clear that the difference between the two strategies is trivial. The reason for this is 87 not clear yet, but we may attribute this similarity to the fact that the middle linear region is taking up most of the temperature profile region in our fairly large systems, as shown in Figure 3 and Figure 4. By also considering the symmetry of the increased temperature in the heat source region and the decreased temperature in the heat sink region, this lack of effect from the nonlinearity may be understood. Table 12.2 Comparison of bulk thermal conductivities at different temperatures from extrapolation results by using different strategies Temperature (K) Thermal Conductivity κ (W/mK) for different data regions used all middle region middle 80% region middle 60% region 300 144.9275 149.2537 153.8462 400 107.5269 108.6957 109.8901 500 84.74576 84.74576 82.64463 600 72.9927 74.07407 74.07407 700 59.88024 60.24096 60.24096 800 55.55556 55.86592 54.64481 900 48.54369 48.78049 51.28205 In the study of the thermal conductivity of 3C-SiC, the phonon group velocity and phonon’s m.f.p. also play important roles, as we can see from the equation (8). To study these two properties of 3C-SiC, we use the “crude estimation” equation (9) by Schelling 81 . Since we have much large simulation system, the value of parameter b will differ, as the assumption of ballistic movement of phonons across the system may not be valid. Consequently, we first calculate the phonon velocity from the elastic constant 88 obtained before and then approximate the m.f.p from the intercept of the extrapolation plot. The calculated group velocity from the elastic constant is 8.64 (10 3 m/s). With this, the m.f.p is approximated and plotted below in Figure 7. estimated mean free path from extrapolation y = 7046.5x -0.9899 R 2 = 0.9974 0 5 10 15 20 25 30 0 200 400 600 800 1000 T (K) mfp (nm) Figure 12-7 The estimated phonon mean free path for 3C-SiC is inversely proportional to temperature. 12.3 Thermal Properties of a-SiC The thermal conductivity of amorphous ceramics has distinctive differences in comparison to crystalline ones. Firstly, the thermal conductivity of amorphous materials increases with increasing temperature, instead of decreasing as is seen in crystalline materials. Secondly, the increase of thermal conductivity of amorphous materials undergoes a plateau region before increasing again to a saturated value. To compare the differences in thermal transport behavior between crystalline and amorphous SiC, we have calculated the specific heat and thermal conductivity of amorphous SiC at various temperatures. 89 12.3.1 Preparation of amorphous SiC In our simulation, a-SiC systems are prepared as follows: we first slowly heat up the 3C-SiC system (49.3235 Å by 279.0144 Å by 49.3235 Å) by scaling the velocity of each atom by a factor of 1.005 every 50 time steps and relaxing the system at every 500 K temperature interval. The system is heated up to 4700K, at which temperature it is completely melted. The system is then relaxed for 50,000 time steps and the size of the system is then expanded by 10% (resulting in a reduced density of 2.893 g/cc). After this, the system is reheated again to counteract to the cooling effect of expanding the system size. Next, the system is well thermalized and relaxed in the melted phase for 50 ps and then cooled down slowly to zero temperature using the velocity scaling factor of 0.997 for every 100 time steps. After applying the conjugate gradient method (CGM) to the cooled system, the system slowly heated again to the desired temperature for calculating the thermal conductivity. A flow chart of the preparation procedure for the amorphous SiC systems is shown in Figure 8. Figure 12-8. The flow chart of the preparation procedure for the amorphous SiC system. 90 12.3.2 Specific Heat of amorphous SiC As an additional check in the study of thermal properties of a-SiC, we have calculated the density of state (d.o.s.) and the specific heat of the a-SiC system by the same procedure used for 3C-SiC. The results are consistent with previous calculations 100 and shown in Figure 9. As we can see from Figure 9 (a), the peaks of the phonon d.o.s are broadened and move down to the lower frequencies (i.e. the acoustic mode), compared to the d.o.s of 3C-SiC. Additionaly, there is no clear gap between the optical mode and the acoustic mode. This is very different than the d.o.s of 3C-SiC. This may imply that the scattering mechanism for phonon propagation in a-SiC could be very different from 3C- SiC. The resulting specific heat does not change too much from 3C-SiC to a-SiC (except at low temperature where the specific heat of a-SiC is a little higher). Also, in the high temperature region, we found the specific heat of a-SiC approaches the theoretically predicted value as shown in Figure 9 (b). At low temperature, the specific heat formula can be approximated as: ∫ = ω ω ω ω d g T k T k k Cv B B B ) ( ) exp( ) ( 2 h h (12) From this formula, it is apparent the d.o.s of a-SiC has most phonons in low frequency range, which results a higher specific heat at low temperatures. 91 Figure 12-9 (a) The calculated vibrational phonon density of states of a-SiC by molecular dynamics using our potential, and (b) The temperature dependence of specific heat (at constant volume) of a-SiC as a result of integration of the d.o.s from (a). The specific heat of 3C-SiC is also shown again for comparison with that of a-SiC 12.3.3 Thermal Conductivity of amorphous SiC The thermal conductivities of a-SiC at various temperatures are also calculated using the NEMD molecular dynamics with P. Jund and R. Jullien’s algorithm 43 . Compared to the simulation for 3C-SiC, the results of this method show much less non- linear effect in the regions near the heat source and heat sink. As is shown in Figure 10, a typical temperature profile for a 700K a-SiC system with steady heat flow shows very good linearity. This system consists of 65,536 atoms with the dimensions 51.0866 Å by 288.989 Å by 51.0866 Å. However, in some smaller size systems at low temperatures, there are singularities in heat source and heat sink. As such, in all the calculations of the temperature gradients for the a-SiC systems, we use the middle 75% of the region between the heat source and the heat sink to ensure all the calculations are consistent. 92 Figure 12-10 Temperature profile for an a-SiC system at 700 K with steady heat flow established along the y-axis. To check the size effect on the bulk thermal conductivity of a-SiC, thermal conductivities of smaller sizes are calculated. In considering the long time requirement for making amorphous SiC systems, a simpler method has been adopted to calculate the thermal conductivities in smaller a-SiC systems. We simply put two and four heat sources and heat sinks into the system’s symmetric positions (taking advantage of the symmetric nature of the original algorithm). A representative temperature profile for four heat sources of the a-SiC system at 500 K is shown in Figure 11. The temperature gradient is calculated from the region between the first heat source and heat sink. Figure 12-11 Temperature profile along the y-axis of an a-SiC system with four heat sinks and sources placed at each 1/8 position along the y-axis. 93 We found that the size effect on calculating the thermal conductivity of a-SiC is not evident at temperatures above 300K. This is apparent in the extrapolation plot for the bulk thermal conductivity of a-SiC at 900K in Figure 12. The thermal conductivity of the largest system from direct calculation is 1.9964 WK -1 m -1 , which is only 3% different than the extrapolated result of bulk thermal conductivity 2.0574 WK -1 m -1 . The limiting of the size effect may be attributed to the shorter phonon m.f.p in a-SiC systems and the large size of the system used in our simulations (~ 289 Å length). Figure 12-12 Extrapolation of the thermal conductivities for a-SiC at 900 K and 50 K from three differently sized systems. The linear-fit result shows that the extrapolated bulk thermal conductivity for 900 K a-SiC is not very different from the result in the largest sample. However, for the system at 50 K, the extrapolated result is considerably different from the result in the largest sample In contrast, at low temperature the size effect is a little more obvious. As is also shown in Figure 12, the extrapolated thermal conductivity of a-SiC at 50 K is 2.0419 WK -1 m -1 , as opposed to the directly calculated result of 1.8679 WK -1 m -1 (a 9.3% difference). While not very different, it is a larger difference than in the case of 900 K. Also, since the steady heat flow is much easier to establish in the a-SiC systems than in 3C-SiC systems, there is no need to run as long. A typical simulation for a-SiC 94 includes 50 ps relaxing the equilibrium system at each given temperature, 250 ps making steady heat flux status and finally 0.5 ns collecting the temperature data. To check the time effect, we simulate 0.5 ns on the 500K a-SiC system (with four pairs of heat sources and heat sinks) consecutively to the previous run. The resulting temperature profile shows that both almost overlap in all regions. A calculation of the temperature gradient for these two time periods is shown in Figure 13. We can see that there is only a 0.4% difference between the two measurements. Figure 12-13 Comparison of the temperature gradient between two consecutive 0.5 ns simulations of a 500K a-SiC system (with four heat sources and sinks) using the first region between heat source and sink. The results show the two linear fits, y1=484.20556+3.16712x (for 2 nd 0.5ns) and y2=483.17733+3.15441x (for 1 st 0.5ns), have only a 0.4% difference in the temperature gradient (slope) In Figure 14, we present the relationship between temperature and thermal conductivity of a-SiC, which is not extrapolated from different sizes. From Figure 14, we can see that at temperature above 300K, there is not much change of the thermal conductivity. The range of the thermal conductivity above 300K is 1.98 ±5% WK -1 m -1 . The temperature dependence of a-SiC thermal conductivity is very different from that of 95 3C-SiC. This phenomenon has been previously addressed before by both experiments and simulations 13,43 . Figure 12-14 Temperature dependence of thermal conductivities of our system with an effective length of 27.9 nm. In the above plot of thermal conductivity of a-SiC at various temperatures, we observed the saturation of thermal conductivity at temperatures above 300K. The saturation phenomenon of the thermal conductivities of amorphous solids has been addressed previously by Gahill from experimental observation 15 , while the limit was determined by using the minimum thermal conductivity model by Slack 86 . As a crude estimation, we would like to use the modified version of Slack’s minimum thermal conductivity model proposed by D. Gahill to calculate the lower bound for thermal conductivity of a-SiC at high temperature region. This minimum thermal conductivity model is also based on equation (8), with the assumptions that the phonon m.f.p. in the amorphous solid is equal to the averaged interatomic distance and that the specific heat is estimated from Dulong and Petit’s model. The phonon velocity is calculated from the elastic constants. 96 In our study, we first calculate the elastic constants of a-SiC from simulation. The averaged phonon velocity is determined to be 6.45 × 10 3 m/s. The specific heat is taken from the above ideal constant volume specific heat 1.24 Jg -1 K -1 . The average interatomic distance is estimated from the density of the system. Finally, the minimum thermal conductivity from this model is determined to be 1.76 WK -1 m -1 , which is in good agreement with the simulated result of 1.98 ±5% WK -1 m -1 . In conclusion, we have calculated the thermal properties of both crystalline and amorphous SiC. The specific heat calculated from the density of states shows that there is little difference in the amount of energy required to increase the temperature of either 3C- SiC or a-SiC by a certain interval. However, the behavior of heat transportation in 3C- SiC and a-SiC is very different, as seen in the differing results of the thermal conductivities as a function of temperature. Nevertheless, our model still successfully calculated the changing trend of thermal conductivity with respect to increasing temperatures for both 3C-SiC and a-SiC. In particular, the results of thermal conductivity of amorphous SiC conform to Slack’s minimum thermal conductivity model at high temperatures. 97 Chapter 4 Combustion Mechanisms of Al Nanoparticle In experiments, thermite reaction involves the oxidization of aluminum to form aluminum oxide, which releases large amount of heat due to the exothermic nature of the reaction. Knowledge of a single nanoparticle’s initiation process in general will greatly enhance the understanding of burning process of MIC composites. Analysis of the experiments have suggested that mechanical breakdown of the nanoparticle’s shell plays important role in the initiation process 52 , but details of this mechanism are still unknown. Previous simulations in our group have shown the oxidation process of aluminum nanoparticle (ANP) under ambient conditions. Formation of alumina shell and the structure details of it are studied. The oxidation process colored coded by temperature is illustrated in Fig. 1. However, combustion of ANP can have different mechanisms other than diffusion, evidenced by the much faster energy release. In this chapter, we will present simulations of the burning process of a single ANP. Combustion mechanisms of ANP under the effects of initial core temperatures and shell structures are discussed. Figure 14-1 Oxidation of aluminum nanoparticle color coded by the local temperature 98 14.1 Combustion of Aluminum Nanoparticles with Crystalline Shell To understand the effect of initial core temperature on the combustion mechanism of ANP, we have conducted million atom Molecular Dynamics simulations of the combustion behavior of a flash heated ANP coated by aluminum oxide shell and surrounded by oxygen atoms. Here, we consider the initial core temperature as the single variable that affects the combustion behavior of the ANP. Atomistic level details of the change of combustion mechanism due to the change of core temperature are discussed in this section. 14.1.1 Preparation of Aluminum Nanoparticles In our simulations, an ANP is prepared as following. First, we cut out an aluminum sphere of 40nm diameter simulated by EAM potential and an alumina shell of 4nm thick simulated by our potential with the inner radius of this shell about 20.5nm. Then, the ANP is prepared by putting the aluminum core and alumina shell together with a distance of about 0.5nm between the shell inner surface and the core aluminum sphere surface. A snapshot of the initial setup of the system is shown in Fig. 2 (a), where white color atoms are aluminum and red ones are oxygen. In this system, there are 1,928,931 core aluminum atoms, 1,152,528 shell aluminum atoms, 1,727,292 oxygen atoms in the shell, and 6,210,716 environmental oxygen atoms surrounding the nanoparticle. In Fig. 2 (a), most of the oxygen atoms surrounding the nanoparticle are removed for the purpose of view. A quarter of the system is cut out to show the inside of the nanoparticle. The whole system is then thermalized for 10ps to let the aluminum core and the alumina shell to form a stable interface in between them. A close up view of the interface is shown in Fig. 99 2 (b), where disorder of the inner surface of the shell is seen, and the mix-up of shell and core atoms forms stable interface. Figure 14-2. (a) Initial configuration of the system with oxygen atoms shown by red color and aluminum atoms shown as white. (b) Mix-up of the core and shell at the interface after 10ps of thermalization. This thermalized system is then embedded inside the oxygen gas environment. After that, the aluminum core is heated to a desired high temperature and thermalized for 2ps while holding the alumina shell atoms fixed. To simulate a nanoparticle with uniformly heated core, we run Langevin Dynamics at the desired temperature for the core atoms for one more picosecond. At this point, the counter for the simulation time is reset to 0ps, and the system is ready to be released for explosion. This procedure corresponds to an adiabatic flash heating of the aluminum core only in experiments, because the wavelength of laser beam can be selected in the near IR regime (~1.06µm), thus absorbed only by the metal core but leave the oxidizer shell unheated. In experiments of laser drilling of alumina by Atanasov et al. 7 , they found that using a Nd:YAG laser 100 source(λ=1.06µm), the temperature on the surface of alumina only increases about 300K in 1 ns, which is far below the experimental melting temperature of alumina 2330 K. So, we have fixed the shell temperature to room temperature at the beginning of the simulations. To study the core temperature effect, we have heated up the core of this ANP to 3000K, 6000K and 9000K respectively. These three systems will be referred to as C3 system, C6 system and C9 system from now on. Constrains on shell atoms of these nanoparticles are then released and the systems undergo subsequent thermal explosion for 200ps in NVE ensemble. 14.1.2 Results and Discussion 4.1.2.1 Energy Release Rate and the Consumption of Aluminum Atoms To compare the different energy release rate in the three systems, change of total kinetic energy per atom is plotted with respect to time in Fig. 3(a). First, we can see that among the three systems, a higher initial core temperature gives more energy release in a total 200ps explosion process. And, at a given time, the energy release rate for the C3 system is much slower than those of the other two systems. This positive relation of higher initial input heat to the aluminum core results faster burning of the nanoparticle is consistent with the experimental results from Dlott’s group 58,59,66,67 . In their experiments, the released energy is not measured directly. But since consumption of ONO 2 indicates their reaction with aluminum fuel, which will generate heat due to the exothermic characteristic of the reaction, the energy release is shown indirectly by the consumption of ONO 2 groups by ALEX fuel. It was observed that, ONO 2 groups are consumed faster with higher laser fluence. Also observed in the experiments is that the consumption of 101 ONO2 groups can be divided into two stages. This two-stage consumption is explained as a first reaction of aluminum fuel with surrounding NC oxidizer and then the slower consumption of NC oxidizer in between aluminum nanoparticles. -0.2 0 0.2 0.4 0.6 0.8 1 1.2 050 100 150 200 time (ps) 9000K 6000K 3000K (a) 0.5 0.6 0.7 0.8 0.9 1 0 50 100 150 200 time (ps) 9000K 6000K 3000K (b) Figure 14-3 (a) Increase of kinetic energy per aluminum atom in the process of explosion for the C3, C6 and C9 systems. (b) Survival fraction of unoxidized aluminum atoms along with time. Except for the positive relation between the burning rate and the initial core temperature, this two-stage burning process is also seen in our simulation results shown in Fig. 3(a). In Fig. 3(a), the total kinetic energy reaches its minimum after crossing the first peak in less than 5ps. In this 5ps, the reaction at the interface of aluminum core and alumina shell increases the kinetic energy in the first one or two picoseconds, while the total kinetic energy drops due to the increase of potential energy by the expansion. After that, the total kinetic energy keeps increasing, and it can be divided into a first fast regime and a second slower regime. In experiments, the fast stage is attributed to the reaction of aluminum core with surrounding NC shell and the slower stage is attributed to the consumption of NC in between nanoparticles. In our simulation, the surrounding oxygen atoms are not energetic chemical reactants, so it’s not directly comparable to the 102 experimental results. However, the consumption of ONO 2 in experiment result indicates its reaction with aluminum fuel, as also seen in the decay of Al emission in the same manner found in their other research result 111 . The consumption of aluminum can be measured directly in our simulation, and the result is plotted in the Fig. 3(b), which clearly shows the characteristic of two-stage consumption of aluminum atoms. These agreements show macroscopically the consistency between the experimental results and our simulation model. Based on that, we would like to study further the fundamental mechanism that causes the fast burning process of ANP. With the help of MD simulation, we can gain atomistic level details of each burning process for the three systems. Specially, even all of the three systems have shown a similar two-stage characteristic in burning process, the difference among them is still obvious. Compared to the C3 system, C6 system and C9 system have much more similarity of energy release behavior. But they still show different trend of energy release near the end of the 200ps run. To be specific, the C6 system shows a flattened trend, but the C9 system shows a continuous increasing trend. To address all these issues, we need a closer look at the status of the systems during their burning processes. And what we found are presented in the following sections. 4.1.2.2 Shell Expansion Dynamics To see the burning process of the three systems, a center slice of each system is cut out and the chronicle snapshots of them are shown in Fig. 4. In this figure, the environmental oxygen atoms are color coded as blue and shell aluminum atoms are color coded as yellow to distinguish them from the red color shell oxygen atoms and white color core aluminum atoms. 103 Figure 4(a) is a series of expansion snapshots of the C3 system at three time steps. At 5ps, the nanoparticle expands and makes rooms for the aluminum core to have some porosity inside it. At 70ps, the expanded nanoparticle shrinks back to a nanoparticle with solid aluminum core. At 200ps, the nanoparticle keeps its original shape, but with some oxygen atoms from the shell diffuse into the surface part of the aluminum core. Here diffusion speed is limited, and the diffusion depth at 200ps is only about 3.8nm. And no aluminum atom is seen outside the nanoparticle. Figure 4(b) is a series of expansion snapshots of the C6 system. At 5ps, the nanoparticle expands to a slightly larger volume as compared to that in C3 system, and leaves more pores observable inside the aluminum core; At 70ps, the nanoparticle expands to an even larger volume. Most of the core aluminum atoms are attached to the shell with only few aluminum atoms left in the center, forming a hollow cavity inside the nanoparticle; At 200ps, the nanoparticle collapsed to form the solid core again, but with lots of oxygen into the center and clusters of aluminum atoms from both shell and core jet out of the nanoparticle. Also, some environmental oxygen atoms diffuse into the nanoparticle in the process of oxidation reaction. As result of this expansion process, the original shell structure is destroyed. This mechanical expansion of the shell has greatly enhanced the diffusion of shell oxygen atoms both outwards and inwards. At 200ps, many shell oxygen atoms are already seen at the core aluminum center, which means that the inward diffusion speed of the oxygen atom could be as high as 100m/s. Meanwhile, the quick diffusion of shell oxygen atoms into the core has left the outside layer of the nanoparticle oxygen deficient and thus ready for more oxidation by environmental oxygen atoms. 104 Figure 4(c) is a series of expansion snapshots of the C9 system. At 5ps, this nanoparticle expands to an even larger volume as compared to the C3 and C6 systems at the same time step, and more pores are also seen inside the core; At 70 ps, the system expands so large that parts of the shell become much thinner than the initial 4nm shell thickness and even large size hole on the shell are observable. Shell oxygen atoms are seen both inside and outside the nanoparticle and core aluminum atoms are seen already outside the nanoparticle; At 200ps, the expanded hollow nanoparticle shrinks due to the contraction of the shell. But the inside of the nanoparticle keeps hollow and even larger holes are seen formed on the shell. Thickness of the remaining shell parts becomes thicker due to the reaction of core aluminum atoms with the shell and the oxidation by the environmental oxygen atoms. Opening of holes on the shell makes paths for core aluminum atoms to react with environmental oxygen atoms directly. Thus, compared to C6 system, mechanical expansion of the nanoparticle has not only destroyed the original shell structure to make the nanoparticle easy for oxidation, but also breaks the shell and makes free path for direct reaction between core aluminum atoms and environmental oxygen atoms. This mechanical breakdown of the shell is important because it has excluded the opportunity for the nanoparticle to be oxidized again to form an outside layer of alumina shell, which will prevent further oxidation to the core aluminum atoms as seen in the previous oxidation simulation studies. Opening of free paths of oxidation generated by the holes ensures more complete oxidation of the nanoparticle at later time, meaning more aluminum fuel would be oxidized. In other words, there will be more energy released in the C9 system compared to that in C6 system. This is consistent with our energy release 105 plot as shown in Fig. 3. In Figure 3, the amount of energy increased per atom is almost the same for both C6 and C9 systems in the first 100ps. But after that, the difference of the amount of energy increased in C6 and C9 system becomes prominent. Figure 14-4. Snapshots of explosion status of the systems at different time steps with different initial core temperature are shown by color coding the core aluminum with white, the shell aluminum with yellow, the shell oxygen with red and the environmental oxygen with blue. (a) expansion snapshots for C3 system at different time steps; (b) expansion snapshots for C6 system at different time steps; (c) expansion snapshots for C9 system at different time steps. Snapshots of the three systems have shown us their explosion behavior qualitatively. To have some quantitative knowledge of the system’s expansion behavior, we plot the expansion of the nanoparticle’s radius with respect to time and the change of the particle shell’s minimum thickness in Fig. 5(a) and 5(b) respectively. The radius of the nanoparticle is an average of the radii of the nanoparticle along the three Cartesian 106 axes. The shell of the nanoparticle is defined as following: We first search through all the atoms in the system and find all the clusters that are covalently bonded. Next, the shell is identified as the largest one among all the clusters. In this calculation, an oxygen atom and an aluminum atom are claimed to be bonded if the distance between them is less than 2.5Å or if the distance is in between 2.5Å and 3.0Å and the two atoms are moving toward each other. After identifying the shell, a virtual shell with thickness zero and radius of the nanoparticle’s radius is constructed. This virtual shell is then divided into areas of same size by a mesh. The mesh firstly discretizes the zenith angle evenly and then divides the azimuth angles accordingly to ensure that each surface area has the same solid angle. Atoms of the shell are assigned to the areas of the mesh according to their coordinates, and the shell thickness of each area can then be determined from the distances of the atoms to the center of the nanoparticle. Finally, the minimum shell thickness of the nanoparticle is the minimum thickness of all the divided areas of the shell. 200 300 400 500 600 0 50 100 150 200 Radius ( ) time (ps) 9000K 6000K 3000K (a) 0 10 20 30 40 50 0 50 100 150 200 Minimum shell thickness (Å) time (ps) 9000K 6000K 3000K (b) Figure 14-5. (a) Plot of the outside radius of the nanoparticle in C3, C6 and C9 systems along with time. (b) Change of shell’s minimum thickness of the three systems along with time. From Fig. 5(a), we can see the difference of expansion behavior in the three systems. The C3 system only experiences a small bulge at the beginning 3ps, after which 107 it starts to increase from around 30ps but in a very slow manner; The C6 system expands during the first 60ps, but it stops to expand after that; The C9 system shows dramatic expansion to more than twice of its original size, and reaches the maximum at around 80ps. Figure 5(b) gives us more insight about the shell’s expansion behavior by comparing the minimum shell thickness between the three systems. The thickness of the C9 system decreases to zero at 22ps, which indicates the appearance of holes, and the holes are never recovered since then. The C6 system also experiences a quick drop of the shell thickness, but it does not reach to zero (appearance of holes) till 60ps. More importantly, the hole disappears quickly after that, and the minimum shell thickness becomes thicker and thicker. The C3 system, on the contrary, never had a time with holes generated in the process of burning. After an initial drop of the shell thickness, it gradually becomes thicker and thicker from around 20ps. The expansion and change of minimum thickness of the shell are seen accompanied by the change of shell structure as shown in Fig. 4. The quantitative analysis of the shell structure change is presented in the following section. 4.1.2.3 Change of shell structure As can be seen from Fig. 4, during the explosion of the nanoparticle, not only the shell expands, but it might also lose atoms by jetting them out and gain atoms through oxidation. Figure 5 has quantitatively shown the expansion of the shell and its fracture. To quantitatively grasp the change of the shell structure, we have checked both the deterioration of the shell due to the loss of shell aluminum atoms and the healing of the shell due to the oxidation. A schematic is shown below for the events involved in the 108 shell structure change. These events include: Al(s → e) representing the shell aluminum atoms that move to the outside environment; Al(s → c) representing the shell aluminum atoms that move to the nanoparticle core; O(s → e) representing the shell oxygen atoms that move to the outside environment; O(s → c) represents the shell oxygen atoms that move to the nanoparticle core; Al(c → s) representing the core aluminum atoms that become part of the shell; O(e → s) representing the environmental oxygen atoms that become part of the shell. In this picture, grey solid circle represents the aluminum atom and red solid circle represents the oxygen atoms. Also shown in the schematic are the two events that atoms transport ballistically between the core and the environment. And they are: Al(c → e) representing the core aluminum atoms that jet out to the environment; O(e → c) representing the environment oxygen atoms that fly to the inside core region. Figure 14-6 Schematic of events happened in the process of aluminum nanoparticle’s detonation. Red solid represents the oxygen atoms and grey solid represent the aluminum atoms. Left white region is the core. Center blue is the Shell region. Right light purple is the outside region of the nanoparticle. Atoms moving are shown as the ones with shadow. Six events of exchange of atoms between shell with either core or environment are Al(s → e), Al(s → c), O(s → e), O(s → c), Al(c → s), and O(e → s). There are also two events Al(c → e) and O(e → c), that represent the direct transport of atoms between core region and environment region when there is a channel present (shown as the light green region in the shell). 109 First, we have checked the deterioration of the shell represented by the loss of shell aluminum atoms that either jet outside the nanoparticle (see Fig. 7 (a)) or are left inside the core regime (see Fig 7 (b)). The blue lines shown in Fig. 7 (a) and (b) show that the C3 system has no aluminum atoms that jet out of the nanoparticle, and very few shell aluminum atoms diffuse into the core gradually. The green lines in Fig. 7 (a) and (b) show that, in C6 systems there are more shell aluminum atoms that diffuse inwards than that jet outside the nanoparticle. The red lines in Fig. 7 (a) and (b) show that the shell aluminum atoms of the C9 system keeps jetting outside the nanoparticle, but the number of atoms moving inward decreases after the first 100ps of increase. 0 2 4 6 8 10 12 050 100 150 200 time(ps) 9000K 6000K 3000K (a) n Al (s J e) 0 2 4 6 8 050 100 150 200 time (ps) 9000K 6000K 3000K (b) n Al (s J c) Figure 14-7. Deterioration of the shells in C3, C6 and C9 nanoparticle system is quantitatively measured by (a) the number of shell aluminum atoms that jet out of the nanoparticle and (b) that fly inside the nanoparticle. The above plots show that hot aluminum core’s impact to the shell causes the shell’s deterioration. At the same time, due to the oxidation reaction on the shell, there is also a healing process happening during the burning process. This process is plotted in Fig. 8. Figure 8(a) shows the oxidation of the shell by environmental oxygen atoms. Here 110 again, C3 system has a very small magnitude due to the fact that its shell was not deteriorated dramatically from its original stable Al 2 O 3 structure except for few diffusion of shell oxygen atoms into the core. In return, it becomes hard for the environmental oxygen atoms to diffuse into the shell. On the contrary, C6 system and C9 system both experience more violent change of the shell structure, and become easy for environmental oxygen atoms to oxidize the shell aluminum atoms. Figure 8(b) shows the reaction of core aluminum atoms with shell. For the C3 system, this number almost kept a constant in the 100ps, and these atoms are actually the layer of core aluminum atoms that are covalently bonded to the oxygen atoms near the inside surface of the shell. There is still an increase, though, due to the slowly diffusion of shell oxygen atoms into the core. For C6 and C9 systems, there are huge amount of core aluminum atoms become part of the shell. This is partially due to their reaction with shell oxygen atoms and partially due to their reaction with environmental oxygen atoms that oxidize the shell aluminum atoms and become part of the shell first. Also, because the number of environmental oxygen atoms on the shell is much larger than the number of core aluminum atoms deposited on the shell, it suggest a possible change of the shell composition in the process of explosion, which will be shown later shortly. 111 0 5 10 15 20 25 30 050 100 150 200 time (ps) 9000K 6000K 3000K (a) n O (e J s) 0 50 100 150 200 050 100 150 200 time (ps) 9000K 6000K 3000K (b) n Al (c J s) Figure 14-8. Recovery of the shells in C3, C6 and C9 nanoparticle system is measured quantitatively by (a) counting the number of environmental oxygen atoms that oxidize the shell and (b) the number of core aluminum atoms that react with oxygen atoms and deposit on the shell. To have the complete information about the atom exchange events happened to the shell as described in Fig. 6. We have also measured the change of total number of atoms on the shell as shown in Fig. 9. Figure 9 shows that, for all three systems, the total number of atoms on the shell has been increased in a either slow or fast manner. This means that oxidation reaction happened on the shell continues to incorporate more number of atoms than that of lost due to their diffusion or jet away from the shell (see Fig. 7 (a) and (b)). Comparison between the three systems shows that, for a higher initial core temperature, the total number of atoms on the shell is higher. Specially, even we have seen that the shell becomes broken for C9 system in Fig. 4 (c), its total number of atoms on the shell is still the largest one among the three systems (shown as the red line in Fig. 9), which indicates it has more reactions happened on the shell due to the oxidation as compared to C3 and C6 systems. By looking at the magnitude of the plots in Fig. 7 and Fig. 8, we know that the major contribution of the increase of total number of 112 atoms on the shell is due to the deposition of core Al atoms to the shell in the 200ps of simulation time. 2.8 3 3.2 3.4 3.6 0 50 100 150 200 time (ps) 9000K 6000K 3000K Figure 14-9. Increasing of the total number of atoms on the shell in C3, C6 and C9 systems. Since the number of atoms on the shell are changed, and there are both aluminum and oxygen atoms deposited on the shell but with a quite different number (see Fig. 8(a) and (b)), the ratio of number of Al to the number of O on the shell might as well be changed. To check the structural change of the shell in terms of its composition, we have plotted in Fig. 10 the composition change of the shell by measuring the number of oxygen atoms on the shell divided by the number of aluminum atoms on the shell. We can see for C9 system, the ratio of number of Al to the number of O quickly increases to 1. This increase of composition is mainly due to the attack and attachment of core aluminum atom to the shell. But for C6 system, this process becomes slower, and it starts to increase at the end of 200ps, which is caused by the close of the pores on the shell and the oxidation by the environmental oxygen gas. For C3 system, the ratio slowly increases due to the diffusion of oxygen atoms to the core, which connects more and more core 113 aluminum atoms to the shell. The continuous increasing for C3 system shows that it is still at the stage of initial diffusion between shell oxygen and core aluminum atoms towards each other, and the oxidization by environmental oxygen atoms to the nanoparticle can come only after the depletion of shell oxygen reaches the outmost layer of the shell, after which more environmental oxygen can then oxidize the nanoparticle. But for the C6 and C9 system, they are already at the stage of being oxidized by large number of the environmental oxygen atoms. 0.6 0.7 0.8 0.9 1 050 100 150 200 time (ps) 9000K 6000K 3000K Figure 14-10. Change of the shell’s composition along with time in C3, C6 and C9 systems by calculating the ratio of number of aluminum atoms on the shell to number of oxygen atoms on the shell. 4.1.2.4 Direct oxidation by crossing the shell Since the energy release is a result of oxidation of aluminum atoms in our system, we have also checked the possibility of direct oxidation of core aluminum atoms by environmental atoms both outside the nanoparticle and inside the core region. This oxidation can happen by either the core aluminum atom jets outside the nanoparticle or the environmental oxygen flies inside the nanoparticle core as seen in Fig. 4 (b) and (c). These two events are also shown in Fig. 6 schematically, where n Al (c → e) represent the 114 number of core aluminum atoms that jet outside the nanoparticle and n O (e → c) represent the number of environmental oxygen atoms that are inside the nanoparticle. And these results are plotted in Fig. 11 (a) and (b), respectively. 10 0 10 1 10 2 10 3 10 4 10 5 10 6 050 100 150 200 time (ps) 9000K 6000K 3000K (a) n Al (c J e) 10 0 10 1 10 2 10 3 10 4 10 5 10 6 0 50 100 150 200 time (ps) 9000K 6000K 3000K (b) n O (e J c) Figure 14-11. Semi-log plot of the mixing of core aluminum atoms and environment oxygen atoms in C3, C6 and C9 system is measured quantitatively by counting (a) the number of core aluminum atoms that jet into the out-particle environment and (b) the number of environmental oxygen atoms that fly into the core. The above measurements shown in Fig. 11 give us the information about the possibility of direct oxidation of core aluminum atoms by environmental oxygen atoms, even though it does not exclude the possibility of these atoms react with shell atoms before they pass through the shell. From Fig. 11, we can see in C3 system, there is neither core aluminum atom nor environmental oxygen atom can pass through the shell. Combined with the Figure 6 and Figure 7, we know that the C3 system has experienced some heat enhanced diffusion and oxidation on the inside and outside surface of the shell, but the shell structure is not changed, which certainly limits the direct oxidation of core aluminum atoms. For C6 system, the number of atoms that pass the shell is in the order of hundreds, showing the limited possibility of direct oxidation reaction between core 115 aluminum and environmental oxygen. In the C9 system, this number is order of magnitude larger than that in C6 system, and we know it’s due to the large opening of holes on the shell that makes a easy path for the atoms’ crossing. To show the enhancement of the number of crossing atoms due to the opening holes on shell, we visualized the core aluminum atoms and the shell by excluding the environmental oxygen gas in Fig. 12. Figure 12 (a) shows the shell of the C6 system color coded by number density of the shell, and Figure 12 (b) shows the shell of the C6 system again with the same color code but together with the yellow colored core aluminum atoms at 100ps. In Fig. 12(c) and 12(d), we visualize the same features of the C9 system at 72ps. Here we can see that most of the core aluminum atoms first jet out from the weak points on the shell, and the C9 system at 72ps has already much more core aluminum jets out of the shell than that in C6 system at 100ps, which is also shown in the Fig. 11 (a). In the case of C9, we have seen quantitatively that the holes formed on the broken shell at early stage (~20ps) due to the expanded large volume and the loss of shell atoms at the weak areas is not able being compensated by the healing process. Because oxidation of the shell happens at localized parts of the shell, the holes are never healed. 116 Figure 14-12. Snapshots of shell morphology of C6 system (a) and C9 system (c) at 100ps with colors representing the number density of the shell as shown by the color bar. Snapshots of the shell morphology together with yellow color core aluminum atoms are shown in (b) for C6 system and (d) for C9 system. Environmental oxygen atoms that are not part of the shell are excluded from the snapshots. Core aluminum atoms are seen jet out of the nanoparticle from the weak areas on the shell during the expansion. From the above comparison, we now know the reason of the energy release rate difference for the three systems shown in Fig. 3(a). Because in our simulated systems, the extra heat can only be generated from the oxidation of aluminum atoms, the higher energy release rate for C9 and C6 system must have a larger amount of reaction for aluminum atoms compared to C3 system. This can be explained after seeing that the shell structure is very little affected in the 200ps for C3 system as shown in Fig. 4(a). Since the shell is stable itself, the oxidization of aluminum atoms happens mostly at the inside interface between the hot aluminum core and the shell due to the elevated inside temperature, and very few surrounding oxygen atoms can oxidize the shell outside surface at near room temperature (see Fig. 7 and Fig. 8). For the C6 and C9 system, there is an early expansion phase, during which structure of the shell is changed and becomes Low High 117 easy for oxidation. In the expansion phase for C6 system, stretching of the shell atoms makes some small pores for core aluminum atoms to jet out and environmental oxygen atoms to flies into the core. These pores are small compared to the holes observed in C9 system, but can still be seen from the Fig. 12(a). Also, in the expansion phase, core aluminum atoms attack the shell and react with shell oxygen atoms, and the reactions are much easier due to the expanded structure of the shell, which is the fast stage of energy release. After the expanded nanoparticle starts to shrink, small pores on the shell are closed and shell is oxidized to form the stable structure again, which limits the number of core aluminum atoms jetting outward as seen in Fig. 11 (b). When the shell becomes condensed again, it slows down the reaction of core aluminum atoms with shell oxygen atoms. Also, speed of the core aluminum atoms jetting out and environmental oxygen gas flies into the core is slowed down. Then the energy release rate becomes slow and it falls to the regime of diffusion again, even though the diffusion and oxidation are at the fast end of the speed limit due to the enhancement by the high temperature of the nanoparticle. In the C9 system, the shell structure is deteriorated by the impact of the hot aluminum core, which makes shell oxygen atoms diffuse inside the core and shell aluminum atoms jet out of the shell. And the openings on the shell give an easy path for the oxidation of core aluminum atoms by the environmental oxygen atoms. This destructive mechanical expansion of the nanoparticle in C9 system gives an even faster energy release in the C9 system compared to the C6 system. The openings are then never closed, which guarantees the continuous oxidization of the core aluminum atoms and explains the continuous increasing trend of energy release at the near end of 200ps for C9 system as compared to the flattened trend in C6 system. 118 From the comparison between the three systems, we can see that the shock induced mechanical breakage of the nanoparticle shell is the mechanism of its enhanced burning efficiency. When the shell becomes broken, it not only helps the oxidation to core aluminum atoms, but also causes the shell’s composition to change (see Fig. 10), which induces oxidation to the shell. 4.1.2.5 Fragment analysis In the previous sections, we have checked carefully the change of the shell structure and its effect on the burning process of the ANP. But the burning of the nanoparticle in experiments is found on the order of several nanoseconds 105 , and then a hot spot is formed due to the huge amount of heat released at each nanoparticle location. The nanosecond time scale is mainly due to the nascent of the most stable product Al 2 O 3 . In our simulations, we have run the explosion of the three systems only up to 200ps due to the limitation of MD timescale. At this time, many intermediate explosion products are observed in our simulation as expected. Simulations of oxidation of aluminum cluster and experiments of laser ablation of Al films have all shown the products of AlO before they form the final product Al 2 O 3 . To check the explosion products in our systems, we have run fragment analysis for the three systems at 200ps. Fragments are decided again by the atoms that are covalently bonded, using the same definition mentioned before. And Bonds between aluminum and aluminum or between oxygen and oxygen are not counted, since we are only interested in the oxide fragments. Detailed results of the fragment analysis are shown in the following sections. But the non-bonded aluminum or oxygen atoms are not shown for the purpose of conciseness. In other words, single atoms are not represented in the statistical results. 119 In Fig. 13, we show the fragment analysis result of explosion products of C3 system at 200ps. Figure 13 (a) is the size distribution of fragments, which shows that most of the fragments in the system are small sized fragments. These small fragments have number of atoms less than 6 in each fragment. There is only one fragment that has size larger than 100 atoms, and it is the shell consisting of 2,978,188 atoms, composition of which can be found in Fig. 10. Figure 13 (b) shows the detail composition of the small fragments. It’s clear that, these small fragments are all aluminum rich fragments, and the majority of the small fragments are Al 2 O. Together with the previous knowledge of the shell behavior we know that these aluminum rich small fragments are the ones formed inside the core due to the diffusion of oxygen atoms into the core and their reactions with core aluminum atoms. In this system at 200ps, the number of unoxidized aluminum atoms is 1,716,567, and they are all inside the nanoparticle. 0 0.5 1 1.5 2 020 40 60 80 100 Fragment size (a) 0 0.2 0.4 0.6 0.8 1 01 2 3 4 Fragment ratio n Al /n O AlO Al 2 O (b) Al 3 O Al 4 O Figure 14-13. (a) Statistics of size distribution of fragments in the C3 system at 200ps shows that the majority of the fragments are small size fragment that has number of atoms in the fragment less than 6. (b) Statistics of fragment composition (ratio of number of aluminum atoms to number of oxygen atoms in a fragment) for fragments with number of atoms less than 6 in C3 system at 200ps. Majority of the small fragments are found to be aluminum rich fragments. 120 Figure 14 shows the fragment analysis result for C6 system. In this case, the largest fragment is again the shell itself with 3,156,289 atoms, but its composition is reduced to 1.12 compared to the value 1.32 in C3 system. Decreasing of the shell composition is due to the deposition of aluminum atoms on the shell as shown in Figure 8. Size distribution analysis shows that majority of the fragments are still small size fragments, and the number of fragments having size smaller than 6 has more than quadrupled compared to that in C3 system. The detailed composition distribution shown in Fig. 14 (b) reveals that this system is a much richer system than the C3 system. Compared to the C3 system, the majority of the small fragments in this system is AlO, even there are still lots of aluminum rich fragments, such as Al 2 O and Al 3 O. Also, oxygen rich fragments are seen in Fig. 14 (b) as contrast to that in C3 system, which indicates the oxidation of the aluminum atoms outside the nanoparticle. As compared to the C3 system at 200ps, the number of unoxidized aluminum atoms is reduced from 1,716,567 to 1,261,128. 0 2 4 6 8 020 40 60 80 100 Fragment size (a) 0 1 2 3 0 123 4 Fragment ratio n Al /n O (b) AlO Al 2 O AlO 2 Al 3 O Figure 14-14. (a) Statistics of size distribution of fragments in the C6 system at 200ps shows that the majority of the fragments are small size fragment with total number of atoms in the fragment less than 6 and the number of small fragments has as much 5 times of that in C3 system (b) Statistics of fragment composition (ratio of number of aluminum atoms to number of oxygen atoms in a fragment) for fragments with number of atoms less than 6 in C6 system at 200ps. As in C3 system, most of the small fragments are AlO or aluminum rich fragments, but there are also many oxygen rich fragments. 121 Figure 15 (a) shows the size distribution of fragments in C9 system at 200ps. Except for the largest shell fragment (total number of atoms on shell is equal to 3,544,013), the number of small fragments is about three times of that in C6 system, indicating a more destructive impact from the higher temperature aluminum core. Figure 15 (b) shows that among all the small fragments, AlO 2 is the dominant fragments, instead of AlO as in the C6 system, indicating an advance to the complete oxidation of the small fragments. Also, at 200ps, the number of unoxidized aluminum atoms has reduced from 1,261,128 in the C6 system to 1,087,524, and all of them are exposed for later oxidation from the opened holes on the shell as seen in Fig. 4 (c). 0 5 10 15 20 25 020 40 60 80 100 Fragment size (a) 0 5 10 15 20 01 2 3 4 Fragment ratio n Al /n O AlO Al 2 O AlO 2 (b) Figure 14-15. (a) Statistics of size distribution of fragments in the C9 system at 200ps shows that most of the fragments are small sized fragments and the number is almost tripled as compared to the number in C6 system at 200ps. (b) Statistics of fragment composition (ratio of number of aluminum atoms to number of oxygen atoms in a fragment) for fragments with number of atoms less than 6 in C9 system at 200ps. The majority of the fragments are not aluminum rich fragments as contrary to the result in C3 and C6 system, and significant amount of the small fragments are oxygen rich fragments (AlO 2 and AlO 3 ). 4.1.2.6 Formation of hot spot From Fig. 5, we know the expansion of the systems have ceased due to the release of the internal pressure of the aluminum core at 200ps. At this time, the explosion process 122 is being transformed from the initiation stage of chemical reaction onset to the ignition stage of the burning process, during which the aluminum fuel will be completely oxidized and energy is to be completely released. A hot spot is then being formed at where the nanoparticle is located. In Fig. 16, we have visualized the temperature profiles of the three systems at 200ps. Figure 16 (a) shows that in C3 system the nanoparticle is at an elevated temperature around 2000K, but the propagation of the heat is strictly limited by the random collision of surrounding oxygen atoms to shell. In the C6 system shown in Fig. 16 (b), the temperature of the nanoparticle is much higher, and the heat is propagated to the nearby oxygen atoms, but the hot spot volume is much smaller compared to that in C9 system. In the C9 system shown in Fig. 16 (c), the high temperature region spreads to almost entire system. Figure 14-16. (a). Snapshots of the cross section of the center for C3 (a), C6 (b) and C9 (c) system with colors representing the temperature profile as shown by the color bar. The C3 system shows a medium temperature sphere localized to the nanoparticle only. The C6 system has high temperature range expanded to a neighborhood of the nanoparticle, but the high temperature regime drops to a medium temperature within twice of the nanoparticle’s radius. The C9 system has high temperature regime well expanded to a much larger sphere with radius as large as 5 times the radius of the original nanoparticle. 0K 6,000K 123 To quantitatively see the temperature profile, we have measured the temperature distribution of the systems at 200ps shown in Fig. 17. The temperature on y-axis of Fig. 17 is an averaged temperature of the atoms in a shell of 2 Å thickness and by a distance of r Å away from the center of the nanoparticle. It can be seen that the temperature profile of C3 system has a sharp drop within 10 Å of where the nanoparticle surface located. Hot regions in the C6 and C9 system are extended more to the outside of the nanoparticle, and the decreasing of temperature from center to nearby region in C6 and C9 system is more continuous. However, the hottest region in the C6 system is still limited to the close region of the nanoparticle, but it spreads out of the nanoparticle far away in the C9 system. The different size and magnitude of the hot spots confirm that the mechanically enhanced detonation process as seen in C9 system forms a much bigger and hotter hot spot as compared to the other two systems where the burning is dominated by the classical diffusion mechanism. Hot spot formed in the C9 system is also extended beyond the region of the nanoparticle size as compared to the localized ones in C3 and C6 systems. 0 2 4 6 8 10 12 0200 400 600 800 r (Angstrom) Figure 14-17. Temperature profile of the C3, C6 and C9 systems at 200ps. Temperature is averaged over a 2A thickness shell from the center of the nanoparticle. This plot quantitatively shows that with a higher initial core temperature, the system will result a larger and hotter hot region inside the system. The vertical bar is where the outside surface of the nanoparticle located, and they are used merely to guide the viewer. 124 14.2 Combustion of Aluminum Nanoparticles with Amorphous Shell In the previous section, we have used crystalline shell as the outside passivation layer of the ANPs, which generates a stronger confinement to the aluminum metal core than in nature. In reality, the shell structure is a more complex format. As has been shown by Ramaswamy et al., 75,76 the natural initial structure of the shell is a layer of about 2.5nm thick and mostly amorphous alumina, which is mostly porous, with the outside surface terminated by hydrogen and the inner layer crystallized. In the experiments by Gertsman and Kwok 29 , it was also observed that the amorphous shell can be partially crystalline. In order to study the shell structure effect by comparing with the previous results of aluminum nanoparticles with crystalline shell (ANPCS), we have run a set of detonation simulations for aluminum nanoparticles with amorphous shells (ANPAS). We have also studied the temperature effect on these ANPAS systems, and compared to the results of ANPCS systems. 14.2.1 Preparation of Aluminum Nanoparticle with Amorphous Shell Setup of the ANPAS systems is made similar to that in section 4.1, except that the shell is amorphous instead of crystalline (see Fig. 18). Figure 18 (b) shows the amorphous structure of the shell as compared to the crystalline shell shown in Fig. 18 (a), which has only disorder at the interface of core and shell. It’s well known that the amorphous alumina is less dense with smaller elastic constants as compared to those values of crystalline alumina. In order to emphasize the effect of weakening of the shell, we also purposely reduced the thickness of the shell from the 4nm in the crystalline case to 2~3 nm in the amorphous case. 125 Figure 14-18. A center slice view for the nanoparticle with crystalline shell (a) and amorphous shell (b). Shell thickness of crystalline shell is about 4nm, and shell thickness of amorphous shell is about 2~3 nm. To make the amorphous shell we amorphorized the crystalline shell by using the following schedule shown in Fig. 19. In this process, the core aluminum atoms are all hold frozen to keep their crystalline structure. Melting temperature of alumina with Vashishta’s potential has been tested to be around 2500K for system with free surfaces 99 . Since for our nanoparticles, the shell is in an open environment, we could heat up the alumina shell up to a temperature of 3500K and relax for a long time to obtain the amorphous structure, as long as making sure no atoms from the shell are overheated too much to break the bonds and fly away. 126 Figure 14-19. Schematic of the schedule of amorphorizing the crystalline shell of Al nanoparticles. Shell thickness is also reduced in this process to 2~3nm. After the shell is amorphorized, structure quantities of the shell are checked carefully to make sure it obtains amorphous structure. The checked structural quantities include the pair distribution and bond angel distribution as shown in Fig. 20. In Fig. 20 (a), the first peak positions for bond length of Al-O, Al-Al, and O-O are 1.83 Å, 2.67 Å, and 3.07 Å, respectively. Figure 20 (b) shows that the two bond angle peaks for O-Al-O are: 92º and 170º, and the two bond angle peaks for Al-O-Al are 95º and 112º. The broadening of the peaks in Fig. 20 (a) and (b) show that the structure of the shell is amorphous. Also, the average coordination number of the Al atoms on the shell is calculated to be 4.8, which is consistent with the experimental results. 127 0 5 10 15 02 468 10 r (Angstrom) O-O Al-O Al-Al (a) 060 120 180 θ (degree) Al-O-Al O-Al-O (b) Figure 14-20. Examination of the shell structure by checking the structural quantities of the shell, including the partial pair distribution (a) and the bond angle distribution (b). The prepared aluminum nanoparticles are then ready to be flash heated to various temperatures for the purpose of studying detonation process as we have done in section 4.1. Here again, while holding the shell atoms frozen, the core of the nanoparticles are heated to 3000K, 6000K and 9000K respectively and these system are then referred to as A3 system, A6 system and A9 system respectively. These systems are then immersed into the well thermalized sea of oxygen at the 10 times the atmosphere pressure and room temperature. After the simulation time is reset to 0 ps, constrains on shell atoms of these nanoparticles are released and the systems are let go subsequent thermal explosion for 200ps in NVE ensemble. 14.2.2 Results and Discussion 4.2.2.1 Energy Release Rate and the Aluminum Atoms Consumption In section 4.1, the simulation of ANPCS systems at three different temperatures, we have seen that with increasing of the initial core temperature, the energy release rate increases accordingly. This proportionality is consistent with the experimental results 2,67 . 128 And in these simulations of ANPAS systems, we have found the same trend as shown in Fig. 21 (a). But compared to those results in ANPCS systems, the energy release in the ANPAS systems is much faster. For a given initial core temperature, the gained kinetic energy per atom in ANPAS systems is much larger than that in ANPCS systems at a given time. As in ANPCS systems, at the beginning of the detonation, there is a quick decrease of the gained kinetic energy for each atom due to the fast expansion of the shell. This is attributed again to the fast expansion of the shell that causes the sudden increase of the potential energy of the shell. Also, among these three systems, the increasing difference of the amount of energy released per atom along with time is similar to that in the ANPCS systems. But in the current simulation result, the magnitude of the energy difference is much bigger than that among the ANPCS systems. Similarly, we have found the corresponding pure Al atom consumption in the current simulations as shown in Fig. 21 (b). It is worth noting that in the A9 system, the consumption of Al atoms is much faster than those in the A6 and A3 systems, and it is almost completely consumed within only 200 ps. To compare the effect of different shells on the energy release rate, we have plotted the energy release in both ANPCS and ANPAS systems together in Fig. 21 (c). Figure 21 (c) shows that, in general the energy release in ANPAS systems is larger than that in ANPCS systems. It is interesting to note that, the upper energy release limit in the ANPCS system (the C9 system represented as the dashed red line in Fig. 21 (c) ) is close to and finally exceeds the lower limit in the ANPAS system (the A3 system represented as the solid blue line in Fig. 21 (c) ). A possible explanation will be addressed in the following section. A comparison of the consumption of pure Al atoms among ANPAS 129 and ANPCS systems is plotted in Fig. 21 (d), which shows the distinction of the Al atom consumption behavior among the C3, A9 and other systems. Among all the simulations, the C3 system has the slowest Al atom consumption in the 200 ps, and the A9 system has the fastest one. Al atom consumptions for all other systems fall between them and can be seen clearly different from those of the C3 and A9 systems. Even though the hierarchy of the magnitude of the released energy in Fig. 21 (c) does not resemble that in Fig. 21 (d) due to the complex oxidation process of Al atoms, Figure 21 (c) and (d) both have shown the distinctive behaviors among the C3, A9 and other systems. Since the previous study in section 4.1 have shown the transition of detonation mechanism from diffusion mechanism (for the C3 system) to the ballistic atom transport mechanism (for the C9 system), here the special behavior of A9 system might also indicate a new mechanism for the fast burning of the ANPs. -0.5 0 0.5 1 1.5 2 2.5 3 3.5 0 50 100 150 200 time (ps) 9000K 6000K 3000K (a) 0 0.2 0.4 0.6 0.8 1 050 100 150 200 time (ps) 9000K 6000K 3000K (b) Figure 14-21. (a) Increase of kinetic energy per aluminum atom in the process of explosion for the A3, A6 and A9 systems. (b) Survival fraction of unoxidized aluminum atoms along with time for A3, A6 and A9 systems. 130 -0.5 0 0.5 1 1.5 2 2.5 3 3.5 A3 A6 A9 C3 C6 C9 0 50 100 150 200 time (ps) (c) 0 0.2 0.4 0.6 0.8 1 050 100 150 200 A3 A6 A9 C3 C6 C9 time (ps) (d) Figure 4-21: Continued (c) Comparison of energy release among ANP systems with crystalline shell and amorphous shell. (d) Comparison of the Al atom consumption among ANPAS and ANPCS systems. Difference among C3, A9 and other systems are also obvious. 4.2.2.2 Shell expansion To see the burning process of the three systems in the current simulation, a center slice of each system is cut out and two snapshots of them at 70 ps and 200 ps are shown in Fig. 22. In this figure, the environmental oxygen atoms are color coded as blue and shell aluminum atoms are color coded as yellow to distinguish them from the red color shell oxygen atoms and white color core aluminum atoms. Figure 22 (a) shows the expansion of the A3 system at 70 and 200 ps. From the initial ANP solid spherical shape, the nanoparticle in the A3 system expands large enough at 70 ps to generate a hollow center inside the nanoparticle. At 200ps, the expansion phase of the nanoparticle is already ended and changed to the shrinking phase. By that time, the shell of the nanoparticle has become unclear from the core, and oxygen can be seen over all regions of the nanoparticle. Figure 22 (a) shows that this expansion process of the A3 system is similar to that in the previous simulated C6 system (see Fig. 4 (b)). 131 Figure 22 (b) shows the expansion of the A6 system at 70 and 200 ps. Expansion of the nanoparticle at 70 ps results a very thin layer of the shell with large holes present. At 200 ps, the oxidation and contraction of the thin shell results the decrease of the nanoparticle’s radius and the increase of the overall shell thickness and the size of the holes. Visualization of the expansion process of this A6 system resembles that of the C9 system (see Fig. 4 (c)). But the shell at the expanding phase is uniform in the A6 system as compared to the non-uniform shape seen in the C9 system at 70 ps. This uniformity of the expanded shell is due to the initial uniform homogenous amorphous shell made for the ANP of A6 system. Because of the uniform expansion of the shell, holes on the shell are generated randomly, instead of the three-fold symmetry holes as seen in the C9 system. Figure 22 (c) shows the expansion of the A9 system at 70 and 200 ps. At 70 ps, the center shell ring is seen shattered by numerous uniformly distributed holes. The shell fragments seen at 70 ps become completely shattered into small clusters as seen at 200 ps. At 200 ps, the A9 system has no nanoparticle shape present at all. The center of the A9 system becomes a spherical region full of small oxygen and aluminum atom clusters. This complete shattering of the shell expose almost all the core aluminum atoms to the environmental oxygen and are ready for oxidation. Previously, we have discussed the reason of slow energy release in the C3 system. It was found that the shell of the nanoparticle in the C3 system is only slightly changed after the small expansion of the core. Atoms on the shell are mostly in a stable configuration. Thus the oxidation of the nanoparticle is mainly through the slow diffusion mechanism. 132 In the discussed C6, A3, C9 and A6 systems, shells of the nanoparticles have all experienced dramatic change in shape and structure as seen in Fig. 22 (a) and (b) and Fig 4 (b) and (c). In their expansion phase, large amount of aluminum atoms are emanated out of the nanoparticle to react with environmental oxygen. Destructive of the shell also causes the core Al atoms to react with under-coordinated oxygen atoms on the shell. These reactions consume the pure aluminum atoms and release large amount of heat. But in all these systems, large trunk of Al metal still exist in form of shell or the collapsed solid nanoparticle, which makes the complete oxidation of all the Al atoms of the nanoparticle difficult. Comparison of energy release among these systems shows that ANPs with weaker shell can release energy more easily (see Fig. 22 (c)). But present of holes on the shell allows further oxidation to core Al atoms, which can explain the exceeded energy release in C9 system to that in the A3 system. When holes are generated in both C9 and A6 systems, the amount of energy released for the system with weaker shell (system A6) is much larger and faster than that in a system with stronger shell (system C9), even the A6 system has a smaller initial driving force (initial core temperature of 6,000K as compared to the initial 9,000K core temperature in C9 system). Destruction and reconstruction of the shell as seen in C6, C9, A3 and A6 systems have resulted possible ballistic transportation of atoms for faster oxidation reaction. The energy release rate is dramatically enhanced as compared to the C3 system (see Fig. 21 above). But this energy release rate is still small as compared to the A9 system. In the expansion of the nanoparticle of the A9 system, it shows clearly different characteristics as compared to those in other systems. Complete shattering of the shell leaves no large aluminum oxide fragments. In addition to the complete exposure of Al atoms for 133 oxidation, the elimination of large aluminum oxide fragments makes the oxidation of all Al atoms much easier. By 200 ps, only about 15% Al atoms are not oxidized. Oxidation of large amount of Al atoms release superior amount of heat in the A9 system as compared to the heat released in other systems, where certain amount of aluminum atoms are either inside the core or attached to the shell and can not be easily oxidized. Figure 14-22. Snapshots of explosion status of the ANPAS systems at two time steps with different initial core temperature are shown by color coding the core aluminum with white, the shell aluminum with yellow, the shell oxygen with red and the environmental oxygen with blue. (a) expansion snapshots for A3 system at 70 and 200 ps; (b) expansion snapshots for A6 system at 70 and 200 ps; (c) expansion snapshots for A9 system at 70 and 200 ps. At 200ps, the A9 system shows dramatic difference in its configuration from the other two, indicating a different oxidation mechanism. 134 A plot of the number density distribution quantitatively shows the distribution of atoms at 200 ps for the A3, A6 and A9 systems, as seen in Fig. 23 (a). Figure 23 plots the radial distribution of the density of the number of atoms in all ANPAS systems. Fig. 23 (a) shows the A9 systems has almost uniform sparse density distribution over the whole system box. This small density may indicate a major gas phase of the systems. As compared to A9 system, the A6 system shows a density peak around 30 nm, which is the thick shell ring shown in Fig. 22 (b). The A3 system has a higher density peak at around 25 nm and a decreasing but non-zero density distribution inwards to the center. 0 0.02 0.04 0.06 0400 800 1200 r (Angstrom) 9000K 6000K 3000K (a) 0 2 4 6 0400 800 1200 r (Angstrom) 9000K 6000K 3000K (b) Figure 14-23. (a) Radial plot of the number density distribution for A3, A6 and A9 systems. (b) Radial plot of the composition distribution for A3, A6 and A9 systems. Fig. 23 (b) plots the ratio of number of Al atoms to Oxygen atoms radially from the center up to 120 nm. It can be seen from Fig. 23 (b) that even the A9 system has a less dense uniform distribution as seen from Fig. 23 (a), it still has a changing character from the center aluminum rich composition to the further away oxygen rich composition. Fig. 23 (b) also shows that inside the A3 and A6 systems, they are both aluminum rich, even the A6 system does not have a solid core. But the lower density of the shell region of the 135 A6 system shown in Fig. 23 (a) indicates an easier energy release situation than the A3 system. To quantitatively compare the different expansion dynamics among the ANPAS systems, we also plot the expansion of the nanoparticle’s radius with respect to time and the change of the particle shell’s minimum thickness in Fig. 24 (a) and (b) respectively. The radius of the nanoparticle is an average of the radii of the nanoparticle along the three Cartesian axes. Shell and the minimum thickness of the shell of the nanoparticles are identified using the same method discussed in section 4.1. In the case A9 system, the shell is later shattered into many small fragments, so the shell of it is identified as the composition of all the fragments having number of atoms larger than 500. This cutoff number for deciding the fragments that belong to the shell is arbitrary, and can cause error in the case of A9 system which shows no clear large fragment present seen from Fig. 22 (c). To reduce the ambiguity caused by the fracture of the shell and the non- uniform expansion of the shell at a later time, we will only show the plot for first 100ps for plots of ANPAS system only. To compare the expansion behavior quantitatively among all ANPAS and ANPCS systems, we have plotted the expansion of the nanoparticle’s radius and the minimum shell thickness for all the systems in Fig. 24 (c) and (d), respectively. Comparison of the expanded nanoparticle radius in the detonation process shows that the C3 system has expanded the least, and the largest radius it reaches is 25.9 nm, about 8% increase from the original 24 nm radius. The minimum shell thickness in the C3 systems goes back close to the original thickness after a first decrease to about 2 nm, a thickness of the shell that still restricts the opportunity of the oxidation of Al core atoms. The expanding of the 136 nanoparticle in C6, C9, A3 and A6 systems are similar to each other in a manner that they all expand firstly and shrinks afterwards. Shrinking of the nanoparticle shows the presence of the shell, where many Al atoms exist and oxidation of them are not easy. Fig. 24 (c) shows the increased expansion of the nanoparticle not only with the increasing initial core temperature, but also with the changing of shell from thicker and harder alumina crystalline to thinner and weaker amorphous phase. Fig. 24 (d) shows that only C9, A6 and A9 system have holes present at the end. Presence of holes enhances the oxidation opportunity of the Al atoms and increases the energy release rate in the detonation process. 200 400 600 800 0 20406080 100 time (ps) 9000K 6000K 3000K (a) 0 5 10 15 20 25 30 0 20406080 100 time (ps) 9000K 6000K 3000K (b) Figure 14-24. (a) Plot of the outside radius of the nanoparticle in A3, A6 and A9 systems along with time. (b) Change of shell’s minimum thickness of the three systems along with time. 137 200 400 600 800 1000 1200 050 100 150 200 A3 A6 A9 C3 C6 C9 time (ps) (c) 0 10 20 30 40 50 050 100 150 200 A3 A6 A9 C3 C6 C9 time (ps) (d) Figure 4-24: Continued (c) Comparison of nanoparticle radius expansion among all ANPAS and ANPCS system. (d) Comparison of minimum shell thickness among all ANPAS and ANPCS systems. 4.2.2.3 Fragment analysis In the previous section, it is seen that in the process of ANP’s detonation the shell is destructed and fragments of atoms are generated. Energy release process seen in Fig. 21 (c) shows three different regimes among the systems. But the magnitude does not resemble the Al consumption shown in Fig. 21 (d). The energy release is complicated due to the complex oxidation process of Al atoms, during which kinds of intermediate oxide products are generated. Experiments of pulsed-laser evaporating Al atoms have shown evidence of various species of aluminum oxides, such as Al 2 O 3 , AlO, AlO 2 , Al 2 O 4 , Al 2 O 2 , Al 2 O, etc 4,67 . Among these products, Al 2 O 3 is usually the final and the most stable product. AlO and AlO 2 are the two major intermediate products to form the final product Al 2 O 3 . In our simulation, the limitation of MD time step size makes it extremely difficult to simulation the whole burning process of the aluminum nanoparticle till its complete oxidation to form the final Al 2 O 3 products for all Al atoms. At the last time step of our simulation 200 ps, the products in the systems are mainly the intermediate products. By 138 performing the fragment analysis, we can check the progress of the aluminum nanoparticle’s oxidation by examining the fragment species, i.e. the oxygen rich fragments has had more reaction happened than the aluminum rich fragments. Also, since smaller fragments (fragments with less atom size) are easier for oxidation, fragment analysis of fragment size distribution can shed some light on the continuous oxidation progress of the nanoparticle. It’s also worth noting that, since we have not calibrated the oxidation process by our potential with either the experimental or the quantum calculations, the intermediate fragment products should only qualitatively but not quantitatively represent the real detonation process. In section 4.1, we have plotted the fragment size distribution for the C3, C6 and C9 systems. We found that, by increasing the initial core temperature, the number of fragment increases dramatically. And the majority parts of the fragments are the small sized fragments. We have found the same trend in the case of A3, A6 and A9 systems again. So, the result is not plotted separately again. Instead, we plotted the total number of fragments of each system at 200 ps together in Fig. 25 for the purpose of comparison among all the systems. We found from Fig. 8 that the total number of fragments in A9 system is superior to all other systems. The following system is the A6 system. Systems C9, A3 and C6 are the next three systems that have similar number of fragments following system A6. System C3 has the least number of fragments among all the systems. 139 0 0.2 0.4 0.6 15 16 17 Systems C3 A3 C6 A6 C9 A9 Figure 14-25. Comparison of the total number of covalently bonded fragments among all the ANPAS and ANPCS systems. Figure 25 shows that the order of number of fragments among the six systems resembles the magnitude hierarchy of the released energy among the systems (see Fig. 21 (c)). This can probably give us some hint about the energy release in the detonation process of aluminum nanoparticles. One possible explanation for the relationship between energy release and number of fragments could be as following: Originally, the system is in a stable structure, and aluminum atoms that can be oxidized are isolated by one piece of large fragment, the solid shell. With the destruction of shell due to the explosion, aluminum atoms are mixed with environmental oxygen atoms by jetting out of the shell or attracting oxygen atoms into the core. This procedure of mixing produces large number of small fragments. For a given number of aluminum atoms, they are more easily oxidized when divided into small fragments than clustered into a large piece, because more aluminum atoms in the first case can be exposed for oxidation. Similarly, for the aluminum oxide fragments in our simulation, each small fragment has aluminum atoms that can be oxidized easily as compared to when they are either inside the shell or 140 isolated into the core, thus more heat are produced when more fragments are generated during the oxidation process. We have also plotted the distribution for fragments with number of atoms less than 6 for the A3, A6 and A9 systems, and show them together with the results of C3, C6 and C9 together. In Fig. 26, 27 and 28, we compare the fragment distribution between systems of A3 and C3, A6 and C6, A9 and C9, respectively. 0 1 2 3 4 5 01 2 3 4 Fragment ratio n Al /n O AlO Al 2 O AlO 2 Al 3 O Amorphous (a) 0 1 2 3 4 5 0 123 4 Fragment ratio n Al /n O AlO Al 2 O Crystalline Al 3 O Al 4 O (b) Figure 14-26. (a) Statistics of fragment composition (ratio of number of aluminum atoms to number of oxygen atoms in a fragment) for fragments with number of atoms less than 6 in A3 system at 200ps, where majority of the small fragments are AlO and AlO 2 . (b) Statistics of fragment composition for fragments with number of atoms less than 6 in C3 system at 200ps, where majority of the small fragments are found to be aluminum rich fragments. 141 0 5 10 15 20 25 30 01 2 3 4 Fragment ratio n Al /n O AlO Al 2 O AlO 2 Amorphous (a) 0 5 10 15 20 25 30 01 2 3 4 Fragment ratio n Al /n O AlO Al 2 O AlO 2 Crystalline (b) Al 3 O Figure 14-27. (a) Statistics of size distribution of fragments in the 6K system at 200ps shows that the majority of the fragments are small size fragment with total number of atoms in the fragment less than 6 and the number of small fragments has as much 5 times of that in 3K system (b) Statistics of fragment composition for fragments with number of atoms less than 6 in 6K system at 200ps. As in 3K system, most of the small fragments are AlO or aluminum rich fragments, but there are also many oxygen rich fragments. 0 20 40 60 80 100 0 123 4 Fragment ratio n Al /n O AlO Al 2 O AlO 2 Amorphous (a) Al 2 O 3 0 5 10 15 20 0 123 4 Fragment ratio n Al /n O AlO Al 2 O AlO 2 Crystalline (b) Figure 14-28. (a) Statistics of size distribution of fragments in the 9K system at 200ps shows that most of the fragments are small sized fragments and the number is almost tripled as compared to the number in 6K system at 200ps. (b) Statistics of fragment composition for fragments with number of atoms less than 6 in 9K system at 200ps. The majority of the fragments are not aluminum rich fragments as contrary to the result in 3K and 6K system, and significant amount of the small fragments are oxygen rich fragments (AlO 2 and AlO 3 ). Figure 26 shows that, for the same initial 3,000K core temperature, when the shell of the aluminum nanoparticle is changed from alpha-alumina to amorphous alumina and 142 the thickness of the shell is reduced, the number of small size fragment is increased by several times (from 32250 to 137755). Other than the increase of the number of small fragments, the composition of the dominant small fragment is also changed from Al 2 O to AlO. Presence of more small fragments and the decreasing of the aluminum rich oxide fragments show that the oxidation process for the A3 system is faster than the C3 systems. Similarly, the number of small fragments also increases dramatically in the case of C6 to A6 system and C9 to A9 system (see Fig. 27 and 28). And the composition of the dominant fragment is changed from AlO to AlO 2 when shell is changed from alpha- alumina to amorphous alumina and the initial core temperature is 6,000K. When the initial core temperature is 9,000K, both the C9 and A9 system have AlO 2 as their dominant small fragments. But the number of small fragments is increased by almost an order of magnitude. And the number of small fragments of Al 2 O 3 in A9 system is also more obvious. It is interesting to note that the dominant fragment feature between systems A6 and C9 and between system A3 and C6 are same. If we sort the systems firstly according to the appearance of oxygen rich fragments as the dominant fragment and then by the magnitude of the number of fragments, the system that have undergone most oxidation reaction is the system A9 and followed by A6, C9, A3, C6 and C3. This order is exactly the same as the hierarchy order of the magnitude of released energy shown in Fig. 21 (c) at the point of 200 ps. This in no way should be the full explanation of the energy release difference among the systems, because the change of the shell structure is a combined factor of its thickness and structure and the current result is just an intermediate step of the detonation process, but it certainly provides us some hint on it. 143 4.2.2.4 Change of shell structure Figure 22 has shown the detonation process of all three ANPAS systems visually. And System A9 shows distinctive characters of detonation as compared to other systems. The energy release rate analysis and fragment analysis show the connection between the amount of energy released and the number of fragments produced in the detonation process. Fragments are mainly produced when the Al atoms jet out of the nanoparticle. Thus, in the process of detonation, destruction and reconstruction of the shell plays important roles for enhancing the energy release rate. Figure 23 (a) has quantitatively shown the system A9 has become almost uniform less dense system. To quantitatively check the dynamics of the shell structure’s change, we first plot the change of the total number of atoms on the largest fragment in the three ANPAS systems (see Fig. 29 below). In the case of A3 and A6 systems, the largest fragment is the shell as can be seen from Fig. 22 (a) and (b). However, there is no obvious large fragment in the A9 system, but only several not-small fragments are visible (see Fig 22 (c)). Plot of the number of atoms on the largest fragment in the system can give us some idea about how the shell is changed in case of A3 and A6 systems and how fast the A9 system breaks into pieces. 144 0 0.5 1 1.5 2 0 50 100 150 200 time (ps) 9000K 6000K 3000K Figure 14-29. Chang of the total number of atoms on the largest fragments of the three ANPAS system. Figure 29 shows that the total number of atoms on the largest fragment (the shell) of A3 system increases continuously. The number of atoms on the shell of A6 system first increases, then decreases significantly and then increases again. The decrease process shows the thinning of the shell as seen in Fig. 22 (b). But system A9 shows a quick drop at around 30 ps, which indicates the break of the shell at the early stage. In both A3 and A6 cases, the largest fragment is the shell and the increase of number of atoms on the shell after 100 ps shows the oxidation of the shell. On the contrary, Fig. 29 quantitatively shows that the largest fragment of the A9 system has number of atoms significantly less than those of A3 and A6 systems, indicating the completely shattering of the shell of the A9 system. To see the details of the shell’s destruction and oxidation, we have examined several events that happened to the shell in the process of detonation as we did in section 4.1. These invents include: ) ( e s Al →, ) ( c s Al →, ) ( s c Al → , and ) ( s e O → . Exact meanings of these symbols have been explained in section 4.1, and these events are plotted in Fig. 30 and 31 below, respectively. Both Fig. 17 and 31 plot the data for the 145 ANPAS system only, instead of comparing them directly with those of ANPCS systems, because the two set of systems originally have different number of atoms on the shell. Also, because the shell breaks completely in case of A9 system, the largest fragment alone does not represent the shell well at the later stage. So instead of using only the largest fragment to identify the shell, a set of non-small fragments (number of atoms on the fragment larger than 500) are used to identify the shell. Also to reduce the ambiguity caused by this definition, these plots show the data for A9 system only up to 100 ps. Figure 30 (a) and (b) plots the events of ) ( e s Al → and ) ( c s Al → to check the deterioration of the shell for ANPAS systems. Similar to the result in ANPCS systems, Fig. 30 shows that increasing of the initial core temperature increases the degree of degradation to the shell. In Fig. 30 (a) the red line representing the A9 system shows a saturation trend after 40 ps due to the fast expansion of the shell and the limitation of the MD box size. The decreasing of the number of shell Al atoms inside the core for A6 systems indicating that when the nanoparticle starts to shrink at around 60 ps (see Fig. 24 (a) or (c)), the shell Al atoms left inside the core during the expansion phase starts to react with the shell oxygen atoms. The reason of this reaction could be due to the shrinking of the core volume which makes the intimate contact for the Al atoms inside the core with the shell. 146 0 5 10 15 20 0 204060 80 100 time (ps) 9000K 6000K 3000K (a) n Al (s J e) 0 5 10 15 20 0 204060 80 100 time (ps) 9000K 6000K 3000K (b) n Al (s J c) Figure 14-30. Deterioration of the shells in A3, A6 and A9 nanoparticle system is quantitatively measured by (a) the number of shell aluminum atoms that jet out of the nanoparticle n Al (s → e) and (b) that are left inside the nanoparticle n Al (s → c). Figure 31 (a) and (b) plots the events of ) ( s c Al → and ) ( s e O → respectively to check the oxidation to the shell for ANPAS systems. Again, similar to the result in ANPCS systems, we find that increasing the initial core temperature can increase the easiness of oxidation to the shell. Exception is seen in Fig. 31 (b) for A9 system, where the decrease of the number of original core atoms on the shell could be caused by the shattering of the shell. 0 10 20 30 40 50 0 204060 80 100 time (ps) 9000K 6000K 3000K (a) n O (e J s) 0 40 80 120 160 0 204060 80 100 time (ps) 9000K 6000K 3000K (b) n Al (c J s) Figure 14-31. Oxidation of the shells in A3, A6 and A9 systems are measured quantitatively by (a) counting the number of environmental oxygen atoms that oxidize the shell n O (e → s) and (b) the number of core aluminum atoms that react with oxygen atoms and deposit on the shell ) ( s c n Al → . 147 The above plots have used either the largest oxide fragment or oxide fragments larger than a threshold for the identification of the shell. But the shell is a complicate structure mixed by oxide fragments and pure aluminum metal fragments. An example of the shell structure of the A6 system is shown in Fig. 32. Figure 32 (a) shows the overall picture of the shell of the aluminum nanoparticle in A6 systems. Figure 32 (b) shows the atomic detail of a piece of the shell, where aluminum oxide fragment is clearly mixed with cluster of pure Al atoms. So, to decide the number of atoms either on the shell or inside or outside the shell, we first use either the largest aluminum oxide fragment or those fragments larger than the threshold to identify the shell and the average radius of the shell and the minimum and maximum thickness of the shell. Then a virtual shell is constructed using the average radius of the fragments as the radius of the virtual shell. The maximum thickness is taken as the thickness of the virtual shell. Figure 14-32. (a) Overall view of the shell structure of the A6 system at 200 ps. Color coding is as following: red represents the oxygen atoms from shell; blue is the oxygen from environment; yellow is the Al atoms from shell; white is the Al atoms from core. (b) Atomistic view of a part of the shell from the left picture (a) shows the complexity of the shell structure is consist of aluminum oxide fragments mixed with pure aluminum metal clusters. Bond length cutoff is 2.5 Å and the color of red represent all oxygen atoms and color of white represent all aluminum atoms. 148 4.2.2.5 Direct oxidation by crossing the shell As seen from section 4.1, the direct oxidation through ballistic transportation of atoms crossing the shell is the main mechanism causing the higher energy release rate in the ANPCS systems. We have also seen from Fig. 21 the same trend for the ANPAS systems, except that the A9 system shows extraordinary enhancement of the energy release rate when compared to other systems. Figure 22 (c) shows the complete shattering of the shell in A9 system, which certainly could improve the opportunity of direct oxidation of Al atoms. To check the number of atoms that cross the shell for oxidation reaction in the ANPAS systems, we plotted the number of core Al atoms and environmental oxygen atoms that cross the shell in the process of detonation shown in Fig. 33 (a) and (b), respectively. Figure 33 (c) and (d) plot the two events again but together with the result of ANPCS systems. 1 10 10 2 10 3 10 4 10 5 10 6 10 7 0 20 406080 100 time (ps) 9000K 6000K 3000K (a) n Al (c J e) 1 10 10 2 10 3 10 4 10 5 10 6 10 7 0 20 406080 100 time (ps) 9000K 6000K 3000K (b) n O (e J c) Figure 14-33. Semi-log plot of the mixing of core aluminum atoms and environment oxygen atoms in A3, A6 and A9 system is measured quantitatively by counting (a) the number of core aluminum atoms that jet into the out-particle environment n Al (c → e) and (b) the number of environmental oxygen atoms that fly into the core n O (e → c). 149 10 0 10 1 10 2 10 3 10 4 10 5 10 6 10 7 0 50 100 150 200 A3 A6 A9 C3 C6 C9 time (ps) (c) n Al (c J e) 10 0 10 1 10 2 10 3 10 4 10 5 10 6 10 7 050 100 150 200 A3 A6 A9 C3 C6 C9 time (ps) (d) n O (e J c) Figure 4-33: Continued. (c) Plot of n Al (c → e) for both ANPCS and ANPAS systems. (d) Plot of n O (e → c) for both ANPCS and ANPAS systems. Comparisons shown in Fig. 33 (c) and (d) illustrate the effect of the shell on the direct oxidation process. In both ANPAS and ANPCS systems, the aluminum core and the outside oxygen environment are similar but the shells are different. Figure 33 (c) shows that for a given initial core temperature, the number of core Al atoms out of the nanoparticle in the ANPAS systems is much more than that of the corresponding nanoparticle in the ANPCS systems. Similarly, Fig. 33 (d) shows the dramatic enhancement for the environmental oxygen to fly into the core. Thus, reducing of the shell thickness and strength from the alpha alumina shell in ANPCS systems to the amorphous alumina shell in ANPAS systems has greatly increased the opportunity for the oxidation reaction between core Al atoms and environmental oxygen atoms. Two special cases among the six systems are the C3 and A9 systems. The C3 system shows no activity of atoms crossing the shell is observed clearly, and the reason as discussed in section 4.1 is contributed to the stable structure of the shell that allows only diffusion of the atoms in the shell instead of direct ballistic transportation. The A9 system shows continuous increasing of the number of environmental oxygen atoms that flies inside the 150 core, as compared to the convergence in other cases. The reason of this is the complete shattering of the shell in the A9 system as also seen visually in Fig. 22 (c). Also, due to the continuous expansion of the large virtual shell identified the layers of non-small aluminum oxide fragments (see Fig. 24 (c)), the number of core Al atoms jet out of the nanoparticle reaches saturation at an very early stage of the nanoparticle’s detonation. The direct transportation of atoms crossing the shell is provided by the paths (holes or weak points) on the shell. Study in section 4.1 has visually shown the core Al atoms jet out of the nanoparticle from the weak points on the shell of both C6 and C9 systems. But the number of weak points on the shell of C6 and C9 systems is limited. And the opening of holes at the weak points on the shell certainly reduces the opening of holes elsewhere on the shell. Nanoparticle with amorphous shells, on the contrary, has a uniform weak shell. In the expansion phase, pores can easily open almost everywhere on the shell. This explains the fact that even if the initial driving force (initial core temperature) of the A6 system is less than that of the C9 system, the number of atoms crossing the shell in A6 system is still much more than that of the C9 system. To compare the way of the holes opened up differently on the shells in A9 and C9 systems, we visualize the core aluminum atoms and the shell by excluding the environmental oxygen gas in Fig. 34 for these two systems. Figure 34 (a) shows the shell of the C9 system color coded by number density of the shell, and Figure 34 (b) shows the shell of the C9 system again with the same color code but together with the yellow colored core aluminum atoms at 100ps. In Fig. 34 (c) and 34 (d), we visualize the same features of the A9 system at 100 ps. Here we can see that the C9 system has six-fold symmetry weak points due to the initial shell structure of crystalline alpha alumina. But 151 the holes on the shell of the A9 systems are everywhere due to the initial uniform shell structure of amorphous alumina. Comparison of Fig 34 (b) and (d) also shows that at a given time (100 ps), the A9 system with a weaker shell has much more core Al atoms jet out of the nanoparticle. To quantitatively compare the large area of the holes generated on the shells of the C9 and A9 systems, we have calculated the percentage of the surface area that the holes occupy. A9 system has holes that occupy 81.55% of the virtual shell surface at 100 ps, but only 24.87% for C9 system at the same time. Figure 14-34. Snapshots of shell morphology of C9 system (a) and A9 system (c) at 100ps with colors representing the number density of the shell as shown by the color bar. Snapshots of the shell morphology together with yellow color core aluminum atoms are shown in (b) for C9 system and (d) for A9 system. Low High 152 4.2.2.6 Formation of hot spot In the laser flash heat induced MIC detonation, sea of hot spots are observed in the process of detonation. Formation of hot spots is considered necessary for the self- sustained detonation process. In section 4.1, we have examined the formation of the hot spot when the nanoparticle is exploded. Here in Fig. 35, we have visualized the temperature profiles of the three ANPAS systems at 200ps. Similar to what we found in Part I, for a given nanoparticle, a larger hot spot is generated when the core is initially heated to a higher temperature. Figure 14-35. Snapshots of the cross section of the center for A3 (a), A6 (b) and A9 (c) system with colors representing the temperature profile as shown by the color bar. With the increase of the initial core temperature, a larger and hotter spot is found. We have also quantitatively checked the temperature profiles for all the three ANPAS systems at 200 ps as shown in Fig. 36 (a). Figure 36 (b) shows the temperature profile comparison between all the six systems. The temperature on y-axis of Fig. 36 is an averaged temperature of the atoms in a shell of 2 Å thickness and by a distance of r Å away from the center of the nanoparticle. Figure 36 (a) shows that for the ANPAS systems, increasing the initial core temperature also results larger and hotter hot spot as 0K 6,000K 153 seen for the ANPCS systems. Figure 36 (b) compares the temperature profiles of both ANPCS and ANPAS systems together and it shows that for a given initial core temperature, when the shell is weakened by changing from the ANPCS systems’ thicker alpha alumina to ANPAS systems’ thinner amorphous alumina, the result hot spot at a given time is much larger and hotter. Figure 36 (b) also shows that the temperature of the shells in C6, A3, C9 and A6 are all around 6,000K, indicating a similar molten state for all the shells of them. But the temperature of shell of system C3 is less than the melting temperature 2,500K of alumina, indicating a solid state of the shell. System A9 does not have a shell, but only high temperature aluminum oxide fragments. 0 5 10 15 20 0400 800 r (Angstrom) 9000K 6000K 3000K (a) 0 5 10 15 20 0400 800 A3 A6 A9 C3 C6 C9 r (Angstrom) (b) Figure 14-36. (a) Temperature profile of the A3, A6 and A9 systems at 200ps. (b) Temperature profile comparison for ANPCS and ANPAS system. Temperature is averaged over a 2A thickness shell from the center of the nanoparticle. The vertical bar is where the outside surface of the nanoparticle located, and they are used merely to guide the viewer. 154 Chapter 5 Conclustion Thermal properties of ceramic SiC have been studied by multi-million-atom MD Simulations. We have calculated the thermal properties of both crystalline and amorphous SiC. The specific heat calculated from the density of states shows that there is little difference in the amount of energy required to increase the temperature of either 3C- SiC or a-SiC by a certain interval. Thermal conductivities of both 3C-SiC and a-SiC materials are simulated using non-equilibrium MD method. Simulation results show that the behavior of heat transportation in 3C-SiC and a-SiC is very different, as shown by the difference in the plots of the thermal conductivities as a function of temperature. The simulated absolute values of thermal conductivities of 3C-SiC at various temperatures are normally smaller than those from experiments. The reason of that could be caused by the difference of thermal expansion behavior between our MD system and the experiments. Nevertheless, our model still successfully calculated the changing trend of thermal conductivity with respect to increasing temperatures for both 3C-SiC and a-SiC. In particular, the results of thermal conductivity of amorphous SiC conform to Slack’s minimum thermal conductivity model at high temperatures. Phase difference of a ceramic has not only effects on the thermal conductivities, but also effects on the combustion behavior of some energetic materials. In this thesis, we have first studied the combustion mechanism of ANP by varying the initial core temperatures of the ANPCS systems, and then the effect of shell structure to the 155 combustion behavior of ANP when the alumina shell is changed from the thicker alpha- alumina state to the thinner amorphous one. We have performed millions-atom MD simulations of thermal explosion of ANPCS systems with the metal core preheated to 3000K, 6000K and 9000K, respectively. We firstly examined the energy release process for the three systems and found the amount and speed of energy release is positively related to the initial core temperatures, which is consistent with experiments. We have also observed a similar two-stage energy release process as that in experiments but in a faster manner. In our simulations, the first fast stage of energy release is attributed mainly to interaction of core aluminum and shell oxygen atoms, diffusion of which are enhanced by the mechanical expansion or broken of the shell. The slower energy release stage is attributed to the later continuous oxidation of the intermediate explosion products. We found for the system starting with a lower initial core temperature (C3 system), the energy release is mostly constrained to the nanoparticle inside by the diffusion mechanism. That is, the oxygen atoms near the inside of the alumina shell firstly diffuse into the core and react with the core aluminum atoms. Then the shell becomes aluminum rich gradually for the later oxidation by environmental oxygen atoms. But the diffusion of shell atoms to the core is slow and by the time of 200ps, the shell is still stable, and only few environmental oxygen atoms are seen to oxidize the shell. For the system with a medium initial core temperature (C6 system), the diffusion and oxidization of shell oxygen atoms to the core aluminum atoms are enhanced by the expansion and shrinking of the nanoparticle. Expansion of the nanoparticle generates less dense part in the shell for the core aluminum atoms to jet out of the nanoparticle and react with environmental oxygen atoms. These reacted atoms 156 form small oxygen rich fragments outside the nanoparticle, and they will be continuously oxidized until the final stable products Al 2 O 3 are formed. But the number of these jetted out aluminum atoms is limited. Meanwhile, the diffusion of shell oxygen atoms to core is greatly enhanced to at least 100m/s due to the expansion and shrinking of the nanoparticle, which quickly depletes the shell oxygen atoms and causes the shell of the nanoparticle become aluminum rich. This aluminum rich shell easily attracts more environmental oxygen atoms to oxidize it, thus the oxidization of the nanoparticle in this case becomes much faster than the C3 system. But due to the shrinking of the nanoparticle, a solid nanoparticle is formed again, and this newly formed shell will later prevent the complete oxidation of the nanoparticle. In other words, the burning of the nanoparticle will be incomplete, thus less efficient. For the system with a high initial core temperature, the expansion of the nanoparticle results broken of the shell. Majority of the small fragments are oxygen rich fragments. Oxidation of the core aluminum atoms is not limited by the diffusion mechanism now. Large holes on the shell open up free paths for the oxidation of core aluminum atoms. Core aluminum atoms are seen jet out of the nanoparticle to react with environmental oxygen, and the speed of these aluminum atoms are in the order of km/s. Also, the unclosed holes ensure the complete oxidation of the core aluminum atoms, which makes the burning of this nanoparticle more efficient than the other two systems. We then studied the effects of shell structure on the combustion behavior of ANPs by running three simulations of ANPAS systems. These three simulations are similar to those of ANPCS systems, except that the shell of ANP is changed from crystalline to amorphous. The results show that, for the ANPAS systems, the temperature effect on the 157 energy release rate, change of the shell structure, intermediate aluminum oxide products are all similar to those of ANPCS systems, i.e. increasing of the initial core temperature results faster energy release process and more violent destruction to the nanoparticle. For a given initial core temperature, comparison between the detonation process of ANPCS and ANPAS systems shows that, when the shell of the aluminum nanoparticle is changed from the hard thicker alpha alumina to weak thinner amorphous alumina, the energy release process becomes much faster. More oxygen rich aluminum oxide intermediate products are generated at a given time, indicating a faster oxidation process of the aluminum nanoparticle toward the final most stable products Al 2 O 3 . Comparison of the energy release rate between all the simulated systems shows that the detonation process can be divided into three groups. The first group is the system C3, which release energy in the process of detonation through diffusion mechanism. The shell of the nanoparticle in system C3 is still solid as indicated by the shell temperature less than the melting temperature of alumina. Shell composition is checked and found close to the stable structure of Al 2 O 3 . Oxidation of the core aluminum atoms is then controlled by the diffusion speed of them through the shell. The second group consists of C6, C9, A3 and A6 systems, which release energy in the process of detonation through mechanism of ballistic transportation of Al atoms to the oxygen environment or vice versa. Nanoparticles of systems in this group have a common character that the shells of these nanoparticles have experienced mechanical destruction due to the large impact force from the core. During the expansion and shrinking phase of the shells, either weak points or holes are generated on the shells to provide ballistic transportation of atoms for the oxidation reaction. In these systems, the shells are destructed, but are still present at the 158 end of the simulation, which leaves certain amount of Al atoms not fully oxidized either inside the shell enclosed core or inside the not-fully oxidized shell itself. The process of full oxidation to all Al atoms needs longer time, because their oxidations are again controlled by the diffusion mechanism. The third group consists of system A9, which has about 85% of its Al atoms either partially or fully oxidized within the 200 ps of simulation time. So, the energy release is dramatically enhanced as compared to the systems in the first and second groups. We have found this fast energy release is due to the complete shattering of its shell in the process of detonation, which leaves almost all the Al atoms either exposed for oxidation or inside small size aluminum oxide fragments. Most of Al atoms are in a cluster consisting of only tens of atoms. Nanochemistry of these superficial Al atom clusters results completely different reactivity as compared to the bulk aluminum. Oxidation of the these Al atoms is much easier resulting the fastest and the most complete oxidation of all Al atoms in the nanoparticle. 159 References: 1 M.P. Allen and D.J. Tildesley, Computer Simulation of Liquids (Oxford Univ. Press, Oxford, 1987). 2 H. C. Andersen, Journal of Chemical Physics 72, 2384-2393 (1980). 3 J. E. Angelo and M. I. Baskes, Interface Science 4, 47-63 (1996). 4 E. F. Archibong and A. St-Amant, Journal of Physical Chemistry A 102, 6877- 6882 (1998). 5 B. W. Asay, S. E. Son, J. R. Busse, and D. M. Oschwald, Propellants Explosives Pyrotechnics 29, 216-219 (2004). 6 Neil W. Ashcroft and N. David Mermin, Solid State Physics, College ed. (Harcourt College Publishers, 1976). 7 P. A. Atanasov, E. D. Eugenieva, and N. N. Nedialkov, Journal of Applied Physics 89, 2013-2016 (2001). 8 C. E. Aumann, G. L. Skofronick, and J. A. Martin, Journal of Vacuum Science & Technology B 13, 1178-1183 (1995). 9 M.I. Baskes, (ASME, New York, 1994). 10 T. Borca-Tasciuc, D. Achimov, W. L. Liu, G. Chen, H. W. Ren, C. H. Lin, and S. S. Pei, Microscale Thermophysical Engineering 5, 225-231 (2001). 11 M. Born and Von Karman, Physik. Z. 13, 297-309 (1912). 12 Mustafa Boyukata and Ziya B. Guvenc, Brazilian Journal of Physics 36, 720-724 (2006). 13 D. G. Cahill, W. K. Ford, K. E. Goodson, G. D. Mahan, A. Majumdar, H. J. Maris, R. Merlin, and S. R. Phillpot, Journal of Applied Physics 93, 793-818 (2003). 14 D. G. Cahill and R. O. Pohl, Annual Review of Physical Chemistry 39, 93-121 (1988). 15 D. G. Cahill and R. O. Pohl, Physical Review B 35, 4067-4073 (1987). 160 16 D. Chaussende, F. Mercier, A. Boulle, F. Conchon, M. Soueidan, G. Ferro, A. Mantzari, A. Andreadou, E. K. Polychroniadis, C. Balloud, S. Juillaguet, J. Camassel, and M. Pons, Journal of Crystal Growth 310, 976-981 (2008). 17 Wai-Kai Chen, The VLSI Handbook Second Edition, 2 ed. (CRC Press, 2006). 18 B.J. Clapsaddle, L. Zhao, A.E. Gash, J.H. Satcher Jr., K.J. Shea, M.L. Pantoya, and R.L. Simpson, in Synthesis and Characterization of Mixed Metal Oxide Nanocomposite Energetic Materials, Warrendale, PA, 2003 (Materials Res. Soc.), p. 91-96. 19 P. J. Daivis and D. J. Evans, Molecular Physics 81, 1289-1295 (1994). 20 M. S. Daw and M. I. Baskes, Physical Review B 29, 6443-6453 (1984). 21 Jonathan T. Desena and kenneth K. Kuo, Journal of Propulsion and Power 15, 794-800 (1999). 22 D. D. Dlott, Materials Science and Technology 22, 463-473 (2006). 23 F. Ercolessi and J. B. Adams, Europhysics Letters 26, 583-588 (1994). 24 Denis J. Evans and Gary P. Morriss, Statistical Mechanics of NonEquilibrium Liquids (Academic Press, London, 1990). 25 M. W. Finnis and J. E. Sinclair, Philosophical Magazine a-Physics of Condensed Matter Structure Defects and Mechanical Properties 50, 45-55 (1984). 26 F. Finocchi, G. Galli, M. Parrinello, and C. M. Bertoni, Physical Review Letters 68, 3044-3047 (1992). 27 Y. V. Frolov, P. F. Pokhil, and V. S. Logachev, Combustion Explosion and Shock Waves 8, 168-187 (1972). 28 D.L. Frost, J.H.S. Lee, and G. Ciccarelli, Shock Waves 1, 99-110 (1991). 29 V. Y. Gertsman and Q. S. M. Kwok, Microscopy and Microanalysis 11, 410-420 (2005). 30 L. Gil, M. A. Ramos, A. Bringer, and U. Buchenau, Physical Review Letters 70, 182-185 (1993). 31 H. Goldschmidt, Zeitschrift Fur Elektrochemie Und Angewandte Physikalische Chemie 14, 558-564 (1908). 161 32 J. J. Granier and M. L. Pantoya, Combustion and Flame 138, 373-383 (2004). 33 J. Guo, D. E. Ellis, and D. J. Lam, Physical Review B 45, 13647-13656 (1992). 34 S. E. Gustafsson and E. Karawacki, Review of Scientific Instruments 54, 744-747 (1983). 35 G. Gutierrez, A. B. Belonoshko, R. Ahuja, and B. Johansson, Physical Review E 61, 2723-2729 (2000). 36 V. Heine, C. Cheng, G.E. Engle, and R.J. Needs, in MRS Symp. Proc., Pittsburgh, PA, 1992 (Materials Research Society). 37 W. G. Hoover, Physical Review A 31, 1695-1697 (1985). 38 http://www.anter.com/TN65.htm. 39 http://www.anter.com/TN67.htm. 40 K. C. Huang and S. H. Ehrman, Langmuir 23, 1419-1426 (2007). 41 W. C. Huang and Y. C. Chen, Journal of Nanoparticle Research 10, 697-702 (2008). 42 G. V. Ivanov and F. Tepper, in Challenges in Propellants and Combustion 100 Years After Nobel, edited by k. K. Kuo (Begell House, New York, 1997), p. 636- 645. 43 P. Jund and R. Jullien, Physical Review B 59, 13707-13711 (1999). 44 A. K. Kapila, R. Menikoff, J. B. Bdzil, S. F. Son, and D. S. Stewart, Physics of Fluids 13, 3002-3024 (2001). 45 P. C. Kelires, Physical Review B 46, 10048-10061 (1992). 46 J. H. Kim, A. Feldman, and D. Novotny, Journal of Applied Physics 86, 3959- 3963 (1999). 47 M. W. Klein, Physical Review B 40, 1918-1925 (1989). 48 O Kordina, (1994). 49 kenneth K. Kuo, Grant A. Risha, Brian J. Evans, and Eric Boyer, in Potential Usage of Energetic Nano-sized Powders for Combustion and Rocket Propulsion, 2004, p. AA1.1.1-12. 162 50 Kenneth K. Kuo and Martin Summerfield, Fundamentals of solid-propellant combustion, Vol. 90 (American Institute of Aeronautics and Astronautics, New York, 1984). 51 J. E. Lennard-Jones, Proceedings of the Physical Society 43, 461-482 (1931). 52 V. I. Levitas, B. W. Asay, S. F. Son, and M. Pantoya, Journal of Applied Physics 101, 20 (2007). 53 J. Li, L. Porter, and S. Yip, Journal of Nuclear Materials 255, 139-152 (1998). 54 Z. Li and R. C. Bradt, Journal of Materials Science 21, 4366-4368 (1986). 55 Z. Lu, K. Nomura, A. Sharma, W. Q. Wang, C. Zhang, A. Nakano, R. Kalia, P. Vashishta, E. Bouchaud, and C. Rountree, Physical Review Letters 95, - (2005). 56 A. Maiti, G. D. Mahan, and S. T. Pantelides, Solid State Communications 102, 517-521 (1997). 57 A. Majumdar, Journal of Heat Transfer-Transactions of the Asme 115, 7-16 (1993). 58 I. Manassidis, A. Devita, and M. J. Gillan, Surface Science 285, L517-L521 (1993). 59 G. J. Martyna, M. L. Klein, and M. Tuckerman, Journal of Chemical Physics 97, 2635-2643 (1992). 60 M. Mench, C. Yeh, and kenneth K. Kuo, in Propellant Burning Rate Enhancement and Thermal Behavior of Ultra-Fine Aluminum Powders (ALEX), Fraunhofer Institute Chemische Technologie, Karlsruhe, Germany, 1998, p. 30-1- 30-15. 61 S. G. Muller, R. Eckstein, J. Fricke, D. Hofmann, R. Hofmann, R. Horn, H. Mehling, and O. Nilsson, Silicon Carbide, Iii-Nitrides and Related Materials, Pts 1 and 2 264-2, 623-626 (1998). 62 F. MullerPlathe, Journal of Chemical Physics 106, 6082-6085 (1997). 63 A. Nakano, R. K. Kalia, and P. Vashishta, Computer Physics Communications 83, 197-214 (1994). 64 A. Nakano, R. K. Kalia, and P. Vashishta, Computing in Science & Engineering 1, 39-47 (1999). 163 65 S. Nose, Molecular Physics 52, 255-268 (1984). 66 R. Orbach, Science 231, 814-819 (1986). 67 M. V. Pak and M. S. Gordon, Journal of Chemical Physics 118, 4471-4476 (2003). 68 Michelle Pantoya, Notable Accomplishments from CTS Awards (2006). 69 M. Parrinello and A. Rahman, Journal of Chemical Physics 76, 2662-2666 (1982). 70 E. Pearson, T. Takai, T. Halicioglu, and W. A. Tiller, Journal of Crystal Growth 70, 33-40 (1984). 71 A. Pechenik, R.K Kalia, and P. Vashishta, (Oxford Univ. Press, Oxford, 1999). 72 Peter Politzer, Pat Lane, and M. Edward Grice, Journal of Physical Chemistry A 105, 7473-7480 (2001). 73 E. K. Polychroniadis, A. Andreadou, and A. Mantazari, Journal of Optoelectronics and Advanced Materials 6, 47-52 (2004). 74 A. Rahman, Physical Review a-General Physics 136, A405-& (1964). 75 A. L. Ramaswamy and P. Kaste, Journal of Energetic Materials 23, 1-25 (2005). 76 A. L. Ramaswamy, P. Kaste, and S. F. Trevino, Journal of Energetic Materials 22, 1-24 (2004). 77 J. H. Rose, J. R. Smith, F. Guinea, and J. Ferrante, Physical Review B 29, 2963- 2969 (1984). 78 C. Rossi, K. Zhang, D. Esteve, P. Alphonse, P. Tailhades, and C. Vahlas, Journal of Microelectromechanical Systems 16, 919-931 (2007). 79 M. A. San Miguel, J. F. Sanz, L. J. Alvarez, and J. A. Odriozola, Physical Review B 58, 2369-2371 (1998). 80 P. K. Schelling, S. R. Phillpot, and P. Keblinski, Physical Review B 65, - (2002). 81 P. K. Schelling, S. R. Phillpot, and P. Keblinski, Physical Review B 65, 144306 (2002). 82 P. K. Schelling, S. R. Phillpot, and P. Keblinski, Applied Physics Letters 80, 2484-2486 (2002). 164 83 U. Schmid, M. Eickhoff, C. Richter, G. Krotz, and D. Schmitt-Landsiedel, Sensors and Actuators a-Physical 94, 87-94 (2001). 84 A. Shimizu, J. Saitou, and Y. J. Hao, Solid State Ionics 38, 261-269 (1990). 85 F. Shimojo, I. Ebbsjo, R. K. Kalia, A. Nakano, J. P. Rino, and P. Vashishta, Physical Review Letters 84, 3338-3341 (2000). 86 G.A. Slack, Solid State Physics, Vol. 34 (Academic, New York, 1979). 87 M. Soueidan and G. Ferro, Advanced Functional Materials 16, 975-979 (2006). 88 G. R. Stewart, Review of Scientific Instruments 54, 1-11 (1983). 89 W. B. Streett, D. J. Tildesley, and G. Saville, Molecular Physics 35, 639-648 (1978). 90 Z. H. Sun, X. X. Wang, A. K. Soh, and H. A. Wu, Modelling and Simulation in Materials Science and Engineering 14, 423-431 (2006). 91 I. Szlufarska, R. K. Kalia, A. Nakano, and P. Vashishta, Journal of Applied Physics 102, - (2007). 92 R.E. Taylor, H. Groot, and J. Ferrier, “Thermophysical Properties of CVD SiC, TRPL 1336,” Report No. TRPL 1336 (1993). 93 A. Tenenbaum, G. Ciccotti, and R. Gallico, Physical Review A 25, 2778-2787 (1982). 94 J. Tersoff, Physical Review B 49, 16349-16352 (1994). 95 T.M. Tillotson, A. E. Gash, R. L. Simpson, L. W. Hrubesh, and J. H. Satcher Jr., Journal of Non-Crystalline Solids 285, 338-345 (2001). 96 K. V. Tretiakov and S. Scandolo, Journal of Chemical Physics 120, 3765-3769 (2004). 97 D. H. Tsai, Journal of Chemical Physics 70, 1375-1382 (1979). 98 A. C. T. van Duin, S. Dasgupta, F. Lorant, and W. A. Goddard, Journal of Physical Chemistry A 105, 9396-9409 (2001). 99 P. Vashishta, R. K. Kalia, A. Nakano, and J. P. Rino, (to be submitted). 165 100 P. Vashishta, R. K. Kalia, A. Nakano, and J. P. Rino, Journal of Applied Physics 101, - (2007). 101 P. Vashishta, R. K. Kalia, A. Nakano, and J. P. Rino, Journal of Applied Physics 103, 083504-1-13 (2008). 102 Arthur F. Voter, in Intermetallic Compounds; Vol. 1, edited by J. H. Westbrook and R. L. Fleischer (John Wiley & Sons Ltd, 1994), p. 77-90. 103 Arthur F. Voter and Shao Ping Chen, in Accurate interatomic potentials for Ni, Al and Ni 3 Al, 1987, p. 175-180. 104 L. L. Wang, Z. A. Munir, and Y. M. Maximov, Journal of Materials Science 28, 3693-3708 (1993). 105 Shufeng Wang, Yanqiang Yang, Zhaoyong Sun, and Dana D. Dlott, Chemical Physics Letters 368, 189-194 (2003). 106 Shufeng Wang, Yanqiang Yang, Hyunung Yu, and Dana D. Dlott, Propellants, Explosives, Pyrotechnics 30, 148-155 (2005). 107 W. J. Weber, L. M. Wang, N. Yu, and N. J. Hess, Materials Science and Engineering a-Structural Materials Properties Microstructure and Processing 253, 62-70 (1998). 108 F.A. Williams, Combustion theory: the fundamental theory of chemically reacting flow systems, 2nd ed. (Benjamin/Cummings Pub. Co., Menlo Park, Calif, 1985). 109 K. Yamada and M. Mohri, Silicon Carbide Ceramics - 1 (Elsevier, London, 1991). 110 Goldberg Yu, M.E. Levinshtein, and Sergey L. Rumyantsev, Properties of Advanced SemiconductorMaterials GaN, AlN, SiC, BN, SiC, SiGe (John Wiley & Sons, Inc., New York, 2001). 111 Hyunung Yu, Selezion A. Hambir, and Dana D. Dlott, in Ultrafast Dynamics of Nanotechnology Energetic Materials, 2006, p. 9.

Abstract (if available)

Abstract Silicon Carbide (SiC) and Alumina (Al2O3) are two important engineering materials for their outstanding electronic, thermal and physical properties. Applications of SiC to high-temperature and igh-power devices require knowledge of its thermal properties. We have studied the specific heat and thermal conductivities of both cubic SiC and amorphous SiC as functions of temperature. Results show different heat transport behaviors for crystalline and amorphous SiC. A layer of about 4 nm Al2O3 is naturally formed on the surface of aluminum nanoparticles (ANPs). Presence of this layer affects the combustion behavior of ANPs. Using multi-million-atom molecular dynamics (MD) simulations, we have investigated the combustion behavior of a single ANP. Transition of the combustion mechanism from diffusion to ballistic transportation of atoms is found when the temperature of the preheated aluminum core of the ANP is changed from 3000K to 9000K. Higher initial temperature of the aluminum core results the mechanical breakdown of the alumina shell in the expansion phase of the ANP.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Heat-initiated oxidation of aluminum nanoparticles

PDF

Nanoindentation of silicon carbide and sulfur-induced embrittlement of nickel

PDF

Silicon carbide ceramic membranes

PDF

Nanomaterials under extreme environments: a study of structural and dynamic properties using reactive molecular dynamics simulations

PDF

Hypervelocity impact damage in alumina

PDF

A process-based molecular model of nano-porous silicon carbide membranes

PDF

Exploring properties of silicon-carbide nanotubes and their composites with polymers

PDF

Dynamics of wing cracks and nanoscale damage in silica glass

PDF

Large-scale molecular dynamics simulations of nano-structured materials

PDF

Atomistic modeling of the mechanical properties of metallic glasses, nanoglasses, and their composites

PDF

Synthesis, characterization, and mechanical properties of nanoporous foams

PDF

Structure and behavior of nano metallic multilayers under thermal and mechanical loading

PDF

Low-dimensional asymmetric crystals: fundamental properties and applications

PDF

The effect of lattice structure and porosity on thermal conductivity of additively-manufactured porous materials

PDF

Processing and properties of phenylethynyl-terminated PMDA-type asymmetric polyimide and composites

PDF

An experimental investigation by optical methods of the physics and chemistry of transient plasma ignition

PDF

The cathode plasma simulation

PDF

Development of a novel heterogeneous flow reactor: soot formation and nanoparticle catalysis

PDF

Molecular-scale studies of mechanical phenomena at the interface between two solid surfaces: from high performance friction to superlubricity and flash heating

PDF

Reactive and quantum molecular dynamics study of materials: from oxidation and amorphization to interfacial charge transfer

Asset Metadata

Creator Wang, Weiqiang (author)

Core Title Thermal properties of silicon carbide and combustion mechanisms of aluminum nanoparticle

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Materials Science

Publication Date 11/21/2008

Defense Date 06/13/2008

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag Al nanoparticle,molecular dynamics,OAI-PMH Harvest,SiC thermal conductivity

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Vashishta, Priya (committee chair), Haas, Stephan (committee member), Kalia, Rajiv K. (committee member), Mansfeld, Florian B. (committee member), Nutt, Steven R. (committee member)

Creator Email wang.weiqiang@gmail.com,wangweiq@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-m1822

Unique identifier UC1155928

Identifier etd-wang-2228 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-135493 (legacy record id),usctheses-m1822 (legacy record id)

Legacy Identifier etd-wang-2228.pdf

Dmrecord 135493

Document Type Dissertation

Rights Wang, Weiqiang

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Repository Name Libraries, University of Southern California

Repository Location Los Angeles, California

Repository Email cisadmin@lib.usc.edu