University of Southern California Dissertations and Theses

Molecular dynamics study of energetic materials

(USC Thesis Other)

Molecular dynamics study of energetic materials

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content MOLECULAR DYNAMICS STUDY OF ENERGETIC MATERIALS by Ying Li A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (MATERIALS SCIENCE) August 2014 Copyright 2014 Ying Li ii Acknowledgements 致给那些在我之前的人，致给那些在我之后的人。 Firstly, I would like to sincerely thank my Ph.D. advisor, Dr. Priya Vashishta, for his continuous guidance, support and warm encouragement. I feel lucky to have had him as my advisor, specially after encouraging me try out every possible resolution of doing research and facing life. The support from Dr. Vashishta would be making you firmly believe that there is still hope for certain kind of answer to what you desire, which you are very close to, but just need to try more out patiently and calmly. This is such great reflection of leadership, which will influence and empower me for lifetime long. I acquired extraordinarily treasure to have had the opportunity to work with Dr. Aiichiro Nakano. Dr. Nakano is an extremely hard working, modest and caring person that I ever know. I have learned lot from him about doing research in the frontier, including literature research, proposal writing, simulation conducting, data analysis, results extrapolation, manuscript constructing. His philosophy towards life “be ambitious and do something important”, “Go for broke” and “Yes, we can!” strongly laid marks in the bottom of my heart along with his unconditionally caring and comforting of me. Dr. Nakano made himself a real life role model to stimulate me pursuing further on the way of academia with curiosities of life and science. I greatly appreciate the valuable advice in research and life from him. I hope that the connection with him will be steady and progressive thorough continuously learning from him as my growing and that I would be able to pass the essential spiritual through educating more. iii I certainly admire the persistence towards physics from Dr. Rajiv K. Kalia. His genuine and fundamental questions in my qualifying exam made me realize that I had lots of imperfection to be improved. The lecture, such as how to correct the angular momentum from a system, he taught me from daily contact would benefit me as an educator forever. I thank Dr. Andrea M. Armani and Dr. Malancha Gupta from Department of Chemical Engineering and Materials Science at University of Southern California (USC) for serving on my disseration committee. In regards the study on ab inito molecular dynamics simulations, I thank Dr. Fuyuki Shimojo from Kumamoto University, Japan for his quantum code and original insights as well as continuous efforts on developing the program. It has been a wonderful and rewarding experience to be a member of CACS (Collaboratory for Advanced Computing and Simulations) for the past six years. I am grateful for the fruitful discussions and collaborations on studies of aluminum nanostructures from Dr. Weiqiang Wang, Dr. Richard Clark and Dr. Adarsh Shekar. I would like to specially thank Dr. Ken-ichi Nomura, who has been really a great support through my PhD life. I deeply appreciate his enthusiasm and his willingness to work alongside by me on debugging the programming code and special talent of visualization and amazing observation of simulation results. Without the help from him, my time spent on PhD career will be doubled. Other group memebers that I have had the pleasure to work with are Dr. Amit Choubey, Dr. Zaoshi Yuan, Dr. Manaschai Kunaseth, Pankaj Rajak, Dr. Hikmet Dursun, Dr. Richard Seymour, Dr. Liu Peng, Dr. Anne Hemeryck, iv Chunyang Sheng, Subodh Tiwari, and Evan Brown. I would like to thank all my colleagues in CACS for their help and inspirational conversations. I deeply appreciate for their interactions and collaborations. I also would like to acknowledge Ms. Patricia Wong for her help in CACS. I gratefully acknowledge the resources that support my research. The work was started with funding support from the Basic Research Program of Defense Threat Reduction Agency (DTRA), and was completed with funding support from the Office of Naval Research (ONR). Simulations were performed at the Center for High Performance Computing and Communications (HPCC) of USC. I thank so many encounters in the past six years life in USC. I cannot think of finer individuals than Ken-ichi Nomura and Amit Choubey, who have been exceptional friends to me both in and out of CACS. My brotherly beloved friend Dr. Liang Xiao had made my graduate years so fun and colorful. Finally, I would like to express my greatest appreciation to my parents for their love. They have been ineffable supportive in my life. This dissertation is dedicated to them. Ying Li University of Southern California Summer 2014 v Table of Contents Acknowledgements ................................................................................................. ii Lists of Figures ..................................................................................................... vii Lists of Tables ...................................................................................................... xiv Abstract ................................................................................................................. xv Chapter 1. Introduction ........................................................................................... 1 1.1 Energetic materials ..................................................................................... 1 1.1.1 Aluminum nanostructures ...................................................................... 1 1.1.2 High explosives .................................................................................... 14 1.2 Molecular dynamics simulation of energetic materials .............................. 22 1.3 Dissertation overview ................................................................................. 23 Chapter 2. Simulation Method: Molecular Dynamics .......................................... 25 2.1 Molecular dynamics simulation .................................................................. 26 2.2 Interatomic force fields ............................................................................... 27 2.2.1 Many-body potentials .......................................................................... 29 2.2.2 Reactive force field (ReaxFF) .............................................................. 34 2.3 Simulation algorithm .................................................................................. 42 2.3.1 Numerical integrator ............................................................................ 42 2.3.2 Linked-list cell decomposition ............................................................. 42 2.3.3 Periodic boundary conditions .............................................................. 44 2.3.4 Parallelization ...................................................................................... 46 2.4 Performance ................................................................................................ 49 Chapter 3. Size Effect on the Oxidation of Single Aluminum Nanoparticle ........ 53 3.1. Background and motivation ....................................................................... 53 3.2 System setup for three single Al-NPs MD simulations .............................. 56 3.3 Analysis of reactive molecular dynamics simulations ............................. 59 3.3.1 Oxidation dynamics of three Al-NPs ................................................... 59 3.3.2 Temperature profile of three Al-NP systems ....................................... 61 3.3.3 Local aluminum atoms concentration .................................................. 63 3.3.4 Formation of Al x O y clusters ................................................................. 69 3.3.5 Evolution of the shell structure ............................................................ 73 3.3.6 Reaction delay ...................................................................................... 76 Chapter 4. Oxidation Dynamics of Aluminum Nanoparticle Aggregates ............ 80 4.1 Background and motivation ........................................................................ 80 4.2 System setup for three Al-NPs aggregates MD simulations ....................... 82 vi 4.3 Analysis of reactive molecular dynamics simulations ................................ 84 4.3.1 Oxidation dynamics of three Al-NP aggregates .................................. 84 4.3.2 Agglomeration of different layers of Al-NPs aggregates .................... 85 4.3.3 Consumption of remaining core Al in Al-NP aggregates .................... 88 4.3.4 Reaction crossover of small intermediate oxidation fragments ........... 89 4.3.5 Evolution of largest oxide fragments in the aggregates ....................... 96 Chapter 5. Oxidation Dynamics of Aluminum Nanorods .................................... 98 5.1 Background and motivation ........................................................................ 98 5.2 Simulation setup for three Al nanorod systems .......................................... 99 5.3 Analysis of the reactive molecular dynamics ........................................... 102 5.3.1 Temperature evolution of the three Al-NRs ...................................... 102 5.3.2 Evolution of remaining core Al atoms for the three Al-NRs ............. 104 5.3.3 Reaction dynamics on the oxide shell of Al-NRs .............................. 105 5.3.4 Oxide shell as a reactive medium ...................................................... 108 5.3.5 Spatial propagation of oxidation reactions ........................................ 109 5.3.6 Size effect on the burning rate of Al-NRs ......................................... 111 5.3.7 Heat propagation in three Al-NRs ..................................................... 113 Chapter 6. Shock-Induced Detonation Reaction in High Explosives RDX ........ 115 6.1 Background and motivation ...................................................................... 115 6.2 Simulation setup for RDX crystal ............................................................. 116 6.3 Analysis of detonation results ................................................................... 118 6.3.1 Temperature and pressure evolution .................................................. 118 6.3.2 Fragment analysis: two-stage reaction times ..................................... 119 6.3.3 Confirmation of the two-stage reaction in an alternative simulation . 121 6.3.4 Reaction pathways analysis ............................................................... 122 6.3.5 Origin of the constituent atoms of N 2 and H 2 O products ................... 125 6.3.6 Reaction products: large clusters ....................................................... 128 6.3.7 Quantum mechanical validation ........................................................ 132 Chapter 7. Conclusion ......................................................................................... 133 References ........................................................................................................... 137 vii Lists of Figures Figure 1.1 Experimentally determined (a) ignition temperature and (b) burning time of aluminum as a function of the used sample size [17]. ........................................ 3 Figure 1.2 (a) High-resolution transmission electron microscope (HR-TEM) image of cross-section of single Al nanoparticle in an oxide shell [29]. (b) Schematic representation of the flame structures observed in aluminum combustion proposed by Bazyn [15]. ..................................................................................... 4 Figure 1.3 Schematic of the process occurring during fast heating of an aluminum nanoparticle covered by a strong oxide shell and representing a melt-dispersion mechanism [24]. .................................................................................................. 5 Figure 1.4 Exposure of Al-NPs to a flash in air and the melt-dispersion mechanism of Al oxidation. (a–e) Schematics and (f–i) TEM images illustrate the oxidation process of the Al-NPs. (a) Initial Al-NPs are covered by a Al 2 O 3 shell. (b and f) Al melts upon rapid heating, which pushes the shell outwardly. (c and g) The shell ruptures and the melted Al becomes in contact with air. (d and i) the large hollow sphere corresponds to the expanded Al-NPs, where most Al has flown away from the particle center and been oxidized at the particle surface. (e and j) the hollow sphere fractures into small clusters of 3–20 nm in sizes. These clusters are consisted of both Al 2 O 3 particles and partially oxidized Al particles. [32] ...................................................................................................................... 6 Figure 1.5 Photographs of (a) agglomerated Ultra Fine Al-NPs (UFAP) and (b) unagglomerated UFAP [34] ................................................................................ 7 Figure 1.6 Comparison of Al-NPs and Al/CuO thermite mixture before and after exposure to the camera flash. (a–c) Optical and scanning electron microscope (SEM) images of Al-NPs before and after the exposure. The original spherical Al-NPs break up into smaller clusters after burning. (d–f) Optical and SEM images of Al–CuO NPs before and after the exposure. The products of the thermite reaction agglomerated into much larger, micron sized particles. [32] .. 8 Table 1.1 Summary of measurements on the Al–CuO and Al–MoO 3 thermites [42]. ....... 9 Figure 1.7 Self-assembled nanoenergetic composite. (a) Schematic of self-assembly of Al-NPs around CuO nanorods using poly(4-viny pyridine). (b) TEM image of CuO nanorods after calcination at 450 o C for 3 hrs. (c) TEM image of Al-NPs assembly around CuO nanorods [46]. ............................................................... 10 Figure 1.8 SEM image of Al nanorods synthesized in a commercial membrane (a) top view and (b) side/titled view [54]. .................................................................... 11 Figure 1.9 HRTEM images of Al nanorods with deposited TiO 2 . (a) A homogenous TiO 2 layer (dark color) is covering the whole surface of the nanorod tip. (b) The TiO 2 layer is also covering the narrow spacing in between the nanorods. (c) There is viii a 4 nm thick Al 2 O 3 layer between the TiO 2 layer and the nanorod [54]. .......... 12 Figure 1.10 Typical SEM images of the Al products: (a) commercial Al powders (b) as- deposited Al nanorods [50]. .............................................................................. 12 Figure 1.11 SEM micrographs of Al nanodisks with (a) a mean diameter D = 61 nm and height h = 20 nm and (b) D = 115 nm and h = 20 nm. (c) Atomic force microscopy (AFM) scan of a single D = 288 nm Al nanodisks. The 3D image and the line scan reveal a rather rough surface with pronounced grainy microstructure [53]. ........................................................................................... 13 Figure 1.12 (a) Angle resolved X-ray photoelectron spectroscopy (XPS) analysis of 20 nm Al thin films, evaporated under the same conditions as in the nanodisk fabrication. The native oxide layer thickness is found to be of the order of 2.5-3 nm and forms within a few hours after air exposure. The red solid line represents a direct logarithmic fit to the XPS data. (b) Schematic depiction of the oxidation process in the supported Al nanodisks used in the experiment. It is assumed that a homogenous oxide layer is formed around the particle surface that are exposed to ambient environment, while no oxidation takes places at the glass-Al interface. The oxide thickness is increasing with time [53]. ............... 14 Figure 1.13 Nanodiamonds synthesized from nanostructured and microstructured explosives. (a–d) TEM micrographs showing nanodiamonds obtained from microstructured (a, b) and nanostructured (c, d) explosive charges. Inset in (a and c) Electronic diffraction patterns showing the diffraction ring of the (111), (220), and (311) diamond planes (d hkl = 0.207 nm, 0.125 nm, 0.107 nm). (e) Size distributions and cumulative frequencies for both samples [58]. .............. 16 Figure 1.14 Snapshot of a chemically sustained shock wave at 15 ps initiated by a four- layer flyer plate with an impact velocity of 6 km/s. The two types of atoms are depicted in black and white. The shock front is propagating from left to right [60]. ................................................................................................................... 17 Figure 1.15 Final structure of the largest carbon cluster formed in 10-fold expanded TATB (d = 0.188 g/cm 3 ) after cooling down from T = 3,000 to 500 K. The configuration was taken from the last step of the sudden expansion simulation [65]. ................................................................................................................... 18 Figure 1.16 Potential energy profile for (a) N-N bond cleavage and (b) the concerted HONO elimination from RDX and the following decomposition of the radical [72]. ................................................................................................................... 19 Figure 1.17 Conductivity profiles in two HMX charges of different grain sizes. Conductivity σ in S/cm, S (Siemens) is 1/Ω. Coordinate x is measured from the detonation front. The records are terminated at 0.8 mm (shot 151, fine-grained) and at 1.8 mm (shot 150, coarse-grained) because of slit shortening. The definition of the peak width Δ is shown in the upper right [77]. ....................... 20 Figure 1.18 Experimental configuration used for launching and probing shock waves in ix pure aluminum films deposited onto glass. Part of the clipped pump pulse (~ 270 ps duration, centered on 800 nm, spot diameter ~20 µm, rise ~12 ps) is absorbed by the aluminum and consequent plasma expansion drives a shock toward the free Al surface. The time history of the motion of the free surface (up to 250 ps) is obtained in one shot by interfering two probe pulses (~ 270 ps, 800 nm, ~200 µm, separated in time by 10 ps, pulse durations not to scale) in a grating spectrometer. Measurements are made separately at the free surface of the two steps (positions 1 and 2). Thicknesses h 1 and h 2 are, respectively, ~ 0.72 µm and 1.44 µm (not to scale with focal diameters) [80]. ................................ 20 Figure 1.19 Diagram of dissertation overview ................................................................. 24 Figure 2.1 (a) Schematic representation of a sample MD box containing N particles; (b) representation of trajectories and velocities of a single particle at different time; (c) representation of velocities and acceleration of a single particle with respect of time. .............................................................................................................. 26 Figure 2.2 Lennar-Jones potential (from http://www.wikipedia.org/) .............................. 29 Figure 2.3 Cell decomposition of a system for linked list algorithm. ............................... 43 Figure 2.4 Data structure for linked list algorithm, where “E” denotes empty. ............... 44 Figure 2.5 Periodically repeated images of an original simulation box (shaded) in 2 dimensions. Minimum image convention applied on primary atom i with imaginary atom j crossing the boundary. .......................................................... 46 Figure 2.6 Schematic view of assigning process ID for domains. .................................... 47 Figure 2.7 A schematic of a projection view of spatial decomposition of 3D system implemented in parallel MD program. (a) Atoms attributes in the extended-cell are cached from 26 neighboring domains. (b) Atoms attributes move out the primary domain to other proper processors. ...................................................... 49 Figure 2.8 Results of the isogranular scaling (2,044,416 atoms per processor) using finite- range MD simulation on silica systems. The weak-scaling parallel efficiency is 0.975 on 131,072 BlueGene/P processors. The weak-scaling parallel efficiency on 212,992 BlueGene/L processors is 0.985 based on the speedup over 4,096 processors [115]. ............................................................................................... 51 Figure 2.9 Execution time of ReaxFF MD on energetic materials (RDX) per time step as a function of the number of processors P using our parallel ReaxFF MD algorithm for energetic materials (RDX) on BlueGene/L (squares) and BlueGene/P (circles) computers, where the number of atoms is N = 16,128P. The efficiency for BlueGene/L is 0.996, which is close to the ideal efficiency, 1 [115]. ................................................................................................................. 52 Figure 3.1 Radial temperature profile of D = 26 nm Al-NP during (a) preparation of aluminum core and alumina shell system, and (b) preheating of the system to 1,100K. The dashed line indicates the interface between the Al-NP and x environmental oxygen. The arrows indicate the time sequence of the preparation and preheating. ............................................................................... 58 Figure 3.2 Snapshots of the central slice for D = 26 nm Al-NP at time = 0.0, 0.2, 0.4, 0.6, 0.8, and 1 ns. ...................................................................................................... 59 Figure 3.3 Snapshots of the central slice for D = 36 nm Al-NP at time = 0.0, 0.2, 0.4, 0.6, 0.8, and 1 ns. ...................................................................................................... 59 Figure 3.4 Snapshots of the central slice for D = 46 nm Al-NP at time = 0.0, 0.2, 0.4, 0.6, 0.8, and 1 ns. ...................................................................................................... 60 Figure 3.5 (a) Temperature of the system (i.e., Al-NP and the environmental oxygen) during the oxidation as a function of time for D = 26, 36 and 46 nm. (b) Radial temperature of the three Al-NPs at time 0.3 ns. Position of the alumina shell is indicated by arrows. .......................................................................................... 61 Figure 3.6 (a) Schematic of a quarter slice of the whole system. Dark blue, grey and yellow circles are aluminum atoms from alumina shell, aluminum atoms from aluminum core, and oxygen atoms, respectively. (b) Local concentration is calculated at a distance R by averaging the concentration within a spherical shell of thickness dR = 10 Å at radius R from the center of the nanoparticle. .. 63 Figure 3.7 Local concentration of Al atoms, C Al Eq. (3.1), at 0 ns, 0.3 ns and 0.6 ns for (a) D = 26 nm, (b) D = 36 nm and (c) D = 46 nm. ............................................ 65 Figure 3.8 Local concentration of Al atoms from the alumina shell, C Al shell Eq. (3.2), at 0 ns, 0.3 ns and 0.6 ns for (a) D = 26 nm, (b) D = 36 nm and (c) D = 46 nm. ..... 67 Figure 3.9 Local concentration of Al atoms from core, C Al core Eq. (3.3), at 0 ns, 0.3 ns and 0.6 ns for (a) D = 26 nm, (b) D = 36 nm and (c) D = 46 nm. ............................ 68 Figure 3.10 Number of (a) Al 2 O, (b) AlO and (c) AlO 2 clusters as a function of time in D = 26 nm respectively. ........................................................................................ 70 Figure 3.11 Number of (a) Al 2 O, (b) AlO and (c) AlO 2 clusters as a function of time in D = 36 nm respectively. ........................................................................................ 71 Figure 3.12 Number of (a) Al 2 O, (b) AlO and (c) AlO 2 clusters as a function of time in D = 46 nm respectively. ........................................................................................ 72 Figure 3.13 Different shell structure of D = 26, 36 and 46 nm systems are illustrated by (a) the radius of Al-NPs as a function of time, and (b) the ratio of oxygen and aluminum within the shell region. ..................................................................... 74 Figure 3.14 (a) Semi-log plot of the number of ejected Al shell atoms as a function of time for D = 26, 36 and 46 nm. (b) Oxide shell temperature during oxidation as a function of time for the three Al-NP systems. Ejection of Al atoms from the shell starts when the shell is in molten state at 2,900 K. ................................... 76 Figure 3.15 (a) The global temperature, (b) The rate of change in the global temperature xi (dT system /dt) and (c) Oxide shell temperature as a function of time for D = 26, 36 and 46 nm .......................................................................................................... 78 Figure 4.1 Transmission electron microscopy (TEM) image of Al nanopowders aggregates [161]. ............................................................................................... 81 Figure 4.2 (a) Core-shell structure of a single Al-NP. An Al core (blue) is covered by a 3 nm thick amorphous Al 2 O 3 shell (orange). (b) Schematic of an aggregate, where Al-NPs are close packed in the ABCBA stacking structure. The five layers are labeled by Roman numerals. ............................................................................. 84 Figure 4.3 Snapshots of Al-NP aggregates for D = 26 nm (a-c), 36 nm (d-f), and 46 nm (g-j) at times 0, 1 and 2 ns. The temperature is color-coded. ............................ 85 Figure 4.4 Time evolution of the center-of-mass positions of Al-NPs in layers I-V in the x direction. (a)-(c) are for Al-NP diameter D = 26, 36 and 46 nm, respectively. 86 Figure 4.5 Time evolution of Al-NP agglomeration as measured by the distance ΔX I-V (t) between the center-of-mass positions of layers I and V in the x direction. The normalized distance ΔX I-V (t)/ ΔX I-V (0) is plotted as a function of time t for D = 26, 36 and 46 nm. .............................................................................................. 87 Figure 4.6 Fraction of unreacted core Al atoms as a function of time for 26 (blue curve), 36 (green curve) and 46 nm (red curve) Al-NP aggregates. Five stages with different core Al consumption rates are shown as regions of different colors and are denoted with Roman numerals. The insets show snapshots of the core Al for the 46 nm Al-NP in layer IV at times t = 1.5 and 1.8 ns. .................................. 89 Figure 4.7 A snapshot of oxide fragments in the 26 nm Al-NP aggregate at t = 2 ns, showing the formation of Al 2 O 3 from AlO and AlO 2 . The largest oxide fragment is the shell shown in the right. ........................................................... 90 Figure 4.8 Number of Al 2 O (a), AlO (b) and AlO 2 (c) fragments as a function of time in the 26 nm Al-NP aggregate. .............................................................................. 91 Figure 4.9 Number of Al 2 O (a), AlO (b) and AlO 2 (c) fragments as a function of time in the 36 nm Al-NP aggregate. .............................................................................. 92 Figure 4.10 Number of Al 2 O (a), AlO (b) and AlO 2 (c) fragments as a function of time in the 46 nm Al-NP aggregate. .............................................................................. 93 Figure 4.11 Time evolution of the fraction of major small oxide fragments, (a) Al 2 O, (b) AlO, and (c) AlO 2 , in 26, 36 and 46 nm Al-NP aggregates, shown in red, green and blue lines, respectively. .............................................................................. 95 Figure 4.12 The ratio between the number of O atoms and that of Al atoms in the shell as function of time for D = 26 nm (red), 36 nm (green) and 46 nm (blue). The initial ratio for the amorphous Al 2 O 3 shell is 1.5. ............................................. 97 Figure 5.1. Schematic of an Al-NR with diameter D and length L = 3D. The Al core (shown in green) is covered with a 3 nm-thick alumina shell (orange), and is xii embedded in oxygen environment (blue). ....................................................... 100 Figure 5.2 Middle-xz slice snapshots of D = 26, 36 and 46 nm Al-NR systems at times 1 ps, 500 ps, and 1 ns. Only shown are Al atoms from the shell (red dots), O atoms from the shell (yellow dots), and environmental O atoms (black dots).102 Figure 5.3 (a) Temperatures of the three Al-NRs as a function of time. (b) The rates of temperature increase dT/dt of the three Al-NRs as a function of time. ........... 103 Figure 5.4 The fraction of remaining core Al atoms as a function of time for the diameter D = 26, 36 and 46 nm. ..................................................................................... 104 Figure 5.5 Snapshots of part of the middle xz slice of the 36 nm system, showing the largest oxide fragment at (a) 50 ps and (b) 200 ps. Original shell Al, shell O, core Al and environmental O atoms are shown in red, yellow, blue and black colors, respectively. ......................................................................................... 106 Figure 5.6 Time evolution of the numbers of atoms from different sources in the oxide shell. The composition of the shell is shown by stacked bar graphs in (a), (c) and (e) for D = 26, 36 and 46 nm Al-NRs, respectively. The numbers atoms from different sources are shown individually in (b), (d) and (f) for D = 26, 36 and 46 nm systems, respectively. Environmental O, core Al, shell O and shell Al atoms are shown in black, blue, yellow and red colors, respectively. ....... 107 Figure 5.7 The ratio of the number of O atoms to that of Al atoms in the oxide shell as a function of time for D = 26, 36 and 46 nm Al-NRs. ...................................... 109 Figure 5.8 (a) Snapshots of the upper half middle xz slice of the largest fragment in D = 36 nm Al-NR at different times. Al and O atoms originated from shell are shown in red and yellow. Al atoms originated from core are shown in blue. The red arrows show the position of the oxidation front. (b) Radial (R) and axial (X) temperature distribution of the same Al-NR. .................................................. 110 Figure 5.9 Oxidation-front (red) and melting-front (blue) positions as a function of time for D = (a) 26, (b) 36 and (c) 46 nm Al-NRs. ................................................. 112 Figure 5.10 Radial temperature distributions at yz cross sections in (a) left, (b) center and (c) right parts of D = 36 nm Al-NR at times 1 ps, 300 ps, 600 ps and 1 ns. The black lines delineate the boundary between core and shell in different sections as the oxidation progress. Also shown are the temperature distributions on the same yz cross sections in (d) left, (e) center and (f) right parts show for D = 36 nm Al-NR, temperature distribution along the radial on Y-Z plane for different sections at different times. ............................................................................... 114 Figure 6.1 (a) RDX molecule, where gray, cyan, blue and red spheres represent C, H, N and O atoms, respectively. (b) RDX unit cell. The 8-molecule (or 168-atom) unit cell has the lattice parameters of a = 13.182 Å, b = 11.574 Å, c = 10.709Å, and α = β = γ = 90˚. (c) Schematic diagram of shock-induced detonation and propagation. ..................................................................................................... 117 xiii Figure 6.2 Colored map of (a) temperature and (b) pressure as a function of the x position and time. .......................................................................................................... 119 Figure 6.3 Number of molecular fragments as a function of time. As the shock front begins to reflect, a rapid production of H 2 O, OH, and N 2 is observed. Shortly after the expansion phase begins, various chemical products such as CO, CO 2 , and NO are produced. ...................................................................................... 121 Figure 6.4 Number of molecular fragments as a function of time in an alternative simulation, where the expansion phase directly follows the compression phase. ......................................................................................................................... 122 Figure 6.5 (a) Initial configuration of RDX molecules in a unit cell, where gray, cyan, blue and red spheres represent C, H, N and O atoms, respectively. Atoms that later will be bonded to form N 2 and H 2 O molecules are highlighted with different colors, i.e., H, O, and N are in green, yellow and magenta, respectively. (b) Molecular fragments formed behind the reversed shock front traveling rightward, where N 2 and H 2 O molecules formed by the highlighted atoms in (a) are enclosed in circles. (c) Fraction of N 2 and H 2 O molecules formed by atoms from single RDX molecules as a function of time. ............. 124 Figure 6.6 Fraction of N 2 molecules that are formed by unimolecular pathways as a function of time. The red line shows the fraction of N 2 formed by N atoms from a single RDX molecule. The blue line shows N 2 originating from an N-N bond within a single RDX molecule. ....................................................................... 126 Figure 6.7 Fraction of H 2 O molecules formed by different reaction pathways as a function of time. The blue line shows H 2 O formed by atoms from a single RDX molecule, i.e., unimolecular pathways. The red line shows H 2 O formed by atoms from two RDX molecules, i.e., dimolecular pathways. ........................ 127 Figure 6.8 Time evolution of large clusters. (a) Fraction of the number of atoms in larger cluster N C vs. the total number of atoms in the system N T . (b) Stoichiometric composition of the large clusters as a function of time. The dashed lines show the normalized stoichiometric number for different elements in the RDX molecule, C 3 H 6 N 6 O 6 . The arrows indicate how the stoichiometric value changed at the final state. (c) Number of large clusters and large C- and O-rich clusters as a function of time. .......................................................................... 131 xiv Lists of Tables Table 1.1 Summary of measurements on the Al–CuO and Al–MoO 3 thermites [42]. ....... 9 Table 1.2 Burn rate of various Al/CuO composites [46]. ................................................... 9 Table 2.1 Parameters for two- and three-body part of our potential used for MD simulation of alumina ........................................................................................ 30 Table 3.1 Reactive molecular dynamics simulation setup for aluminum (core) and amorphous alumina (shell) nanoparticles in oxygen environment .................... 56 Table 3.2 Relationship between the delay time and the nanoparticle size ........................ 79 Table 5.1 Number of atoms of different components in D = 26, 36 and 46 nm aluminum nanorod (Al-NR) systems. ............................................................................... 100 xv Abstract This dissertation is focused on molecular dynamics (MD) simulation of energetic materials. Energetic materials are a class of materials with high amount of stored chemical energy that can be released. Propellants, explosives, and pyrotechnic are typical classes of energetic materials. In the past decade, energetic materials have been playing very important roles in the fields of defense and other technologies, since they are secure, convenient and energetically compact, which will be continuously functionalized for further improvements indispensably in the new era. The significance of how to utilize the novel and promising features more profitably and constructively will necessitate fundamental understanding of dynamics of the material releasing energy through reactions. Among all kinds of energetic materials, aluminum nanostructures and high explosives are the subjects studied in this dissertation. For aluminum nanostructures, we have performed multimillion atom reactive MD simulations to investigate the oxidation dynamics of aluminum nanoparticles (Al-NPs), Al-NP aggregates, and aluminum nanorods (Al-NRs). For high explosives, we have addressed the shock-induced detonation of cyclotrimethylenetrinitramine (RDX, an initialism for Research Department Explosive). The scientific interest in aluminum nanostructures research basically lies in the increased reactivity during oxidation as compared with conventional micro-sized structures. However, what is the size effect on the reactivities of the Al-NPs still puzzles researchers. Practically, Al-NPs commercially produced always have a thin oxide layer xvi encapsulating the pure metal particles. The oxide shell has been included in all our simulated aluminum nanostructures to accurately represent the reality as much as possible. We have systematically simulated different sizes of single Al-NPs (diameter, D = 26, 36 and 46 nm) to study their reactivities in oxygen environments as isolated individuals. We have also systematically extended our study to their assemblies – Al-NP aggregates with different sizes in oxygen environments. Aluminum nanorods (Al-NRs), a different type of morphology, also have been investigated due to their high contact areas with oxidizers. Our MD simulations reveal, within the range of our studied cases, the smaller sizes of aluminum nanostructures have higher reactivities (e.g., faster reaction rate, shorter reaction time and faster oxidation speed). Meanwhile, our atomistic analyses show interesting and unexpected roles of the interfaces of the oxide shell with aluminum and oxygen environment during the initial oxidation. In addition, the processes of final product formation show the evolution of different major intermediate species. Such atomistic understanding of the interface, reaction pathways and the size effect should have profound impacts on the design of the efficient burning of intermolecular products. The second topic of this dissertation is the detonation of a high explosive, RDX. Despite many years of intense research, atomistic mechanisms underlying the reaction time and intermediate reaction products of detonating high explosives have been elusive. Our reactive MD simulations reveal a novel two-stage reaction mechanism during the detonation of RDX crystal. Rapid production of N 2 and H 2 O is followed by delayed production of CO molecules. We have found that further decomposition towards the final products is inhibited by the formation of large, metastable carbon- and oxygen-rich xvii clusters. Such atomistic understanding may pave a way toward rational design for detonation synthesis of materials. 1 Chapter 1. Introduction 1.1 Energetic materials Energetic materials have drawn great attentions ever since the 9 th century in China [1]. However, the increased availability and investigation of the energetic materials did not emerge until the appearance of recent improvised devices [2, 3]. A wide variety of energetic materials can release energy; only a small number of them can be mass- produced safely, inexpensively, stably, and non-toxically. The common types of energetic materials are explosives, pyrotechnic compositions, and propellants. The primary focuses in this dissertation on the energetic materials are metallic additives (e.g., aluminum) [4, 5] and energetic compounds (e.g., high explosives) [6, 7]. 1.1.1 Aluminum nanostructures Aluminum (Al) is the most abundant metal element on the earth. It is a lightweight metal that is used to make airplanes, cars, bicycles, building, electronics, and many other things. In nature, aluminum is combined with other elements, since it is too reactive to present as free metal. The aluminum objects that we use daily, like soda cans, have a protective oxide coating that prevents this reaction from happening. The oxidation of Al is an extremely powerful exothermic reaction that can give off vast amount of heat, 83.8 kJ/cm 3 for Al to become aluminum oxide (Al 2 O 3 ) [8]. Combining the properties of thin layer of protective oxide and high volumetric energy density, it can be used as solid rocket propellants. Because of the high volumetric energy density [8], low production 2 cost, and the formation of protective oxide layer for safe storage [9, 10], aluminum composites have been widely employed as energetic additives in propellants and explosives [11-13]. Composites containing micron-size Al have been applied as suitable fuel in practical system, because it is inexpensive, is safe to handle, and contains more reactive pure Al metal in single entirety due to lower fraction of dead mass and volume of inert Al 2 O 3 shell. There have been experiments conducted on the Al micron-size composites about the oxidation mechanism [14-16]. However, as a result of the ability to routinely synthesize and characterize the aluminum nanostructures, there is growing interests in Al composites with the features reduced to the nanometer scale [17-19]. In particular, aluminum nanoparticles (Al-NP) have attracted much attention as promising nanoenergetic materials due to their high reactive power [20] and low melting temperature [17]. Al-NPs are widely used for propellants and explosives because of their high volumetric energy density [8], high heat-release rate [21], low ignition temperature [22], and rapid heat propagation [23]. A number of experiments were performed in an attempt to understand the combustion mechanisms based on a single Al-NP [24, 25]. For example, Yetter et al. [17] found that Al nanostructures actually have more advantages as metallic additives in the field of combustion. Ivanov et al. [26] reported that the addition of Al-NPs into propellant greatly enhances the burning rate. Aumann et al. [10] showed that the activation energy for the oxidation of aluminum nanopowder is lower than that for aluminum bulk. The higher surface-to-volume ratio in nanostructures leads to lower melting temperatures, lower sintering temperatures, and higher theoretical densities 3 comparing with Al microparticles. Trunov et al. [27] showed that the ignition temperature is much lower when the size of Al powder is reduced to the nanometer scale. In Figure 1.1, we see the lower ignition temperature and shorter burning time for smaller diameter of Al particles. Figure 1.1 Experimentally determined (a) ignition temperature and (b) burning time of aluminum as a function of the used sample size [17]. Theoretical study about the combustion of aluminum on a wide range of diameters by Huang et al. [28] claimed that there is a transition from a diffusion- controlled to a kinetically controlled model as the size of particle decreasing from micron to nano scale. Bazyn et al. (see Figure 1.2(b)) have concluded that there should be a surface process for the Al particles less than 10 µm at the beginning of combustion to explain the lowering combustion temperature for Al-NPs. They explained the phenomena as oxidizer diffuses to the nanoparticle surface and then through the aluminum oxide layer of the particle (see the right panel in Figure 1.2(b)). Accordingly, the combustion 4 temperature is highest inside the particle, not related to aluminum oxide dissociation temperature or ambient temperature. Figure 1.2 (a) High-resolution transmission electron microscope (HR-TEM) image of cross-section of single Al nanoparticle in an oxide shell [29]. (b) Schematic representation of the flame structures observed in aluminum combustion proposed by Bazyn [15]. Levitas et al. [24, 30, 31] proposed a melt-dispersion mechanism to explain the combustion of Al-NPs at fast heating rate. Figure 1.3 illustrates the schema of the combustion as firstly, rapid melting of the Al core that leads to shell spallation. Then, an unloading wave propagates to the center of Al molten core and generates tensile pressure, which disperses small Al clusters. However, since the combustion of Al-NP is a dynamic process at relative high temperatures (at least above 1,000 K, see Figure 1.1(a)), it is extremely difficult to capture the reaction mechanism through experiments. 5 Figure 1.3 Schematic of the process occurring during fast heating of an aluminum nanoparticle covered by a strong oxide shell and representing a melt-dispersion mechanism [24]. Ohkura et al. at Stanford University [32] found a direct experimental evidence of the melt-dispersion mechanism though TEM. Their analysis reveals that the oxidation mechanism of Al-NPs under the heating rate of 10 6 K/s or higher is the rupture of the amorphous Al 2 O 3 shell leading to the contact of melted Al with the air/oxidizer. The images from TEM show that, towards the end of the burning of Al-NPs, the finial products are smaller fragmented clusters of both Al 2 O 3 particles and partially oxidized Al particles. 6 Figure 1.4 Exposure of Al-NPs to a flash in air and the melt-dispersion mechanism of Al oxidation. (a–e) Schematics and (f–i) TEM images illustrate the oxidation process of the Al-NPs. (a) Initial Al-NPs are covered by a Al 2O 3 shell. (b and f) Al melts upon rapid heating, which pushes the shell outwardly. (c and g) The shell ruptures and the melted Al becomes in contact with air. (d and i) the large hollow sphere corresponds to the expanded Al-NPs, where most Al has flown away from the particle center and been oxidized at the particle surface. (e and j) the hollow sphere fractures into small clusters of 3–20 nm in sizes. These clusters are consisted of both Al 2O 3 particles and partially oxidized Al particles. [32] In real applications, Al-NPs always form into group/aggregations. During combustion, Al-NPs nanoparticles tend to melt and coalesce into large agglomerates [33, 34], which makes the combustion behavior highly complex [35, 36]. The dynamics of such collective combustion has never been explained. As shown in Figure 1.6, Ohkura et al. [32] discovered differences of final combustion products Al-NPs in air and with CuO NPs. Figure 1.6 (c) shows that the original spherical Al-NPs break up into smaller 7 clusters after burning. Figure 1.6(f) shows that the products of the thermite reaction agglomerated into much larger, micron sized particles. Even though the oxidation conditions for the two sets of experiments are different — one is Al-NPs reacting with air (gas phase oxidizer), while the other is Al-NPs reacting with CuO-NPs (solid phase oxidizer) — the fundamental issues are why are there different sizes of final products, and what is the size effect on the burning of Al-NP aggregates? Since none of the experiments can eliminate the effect of size distributions of Al-NPs and experiments confirmed a relationship between the primary size of Al-NPs and their reactivities [25, 28, 37], it is essential to know the effect of aggregated Al-NPs’ size on the oxidation dynamics. Since the oxidization processes of multiple Al-NPs are more complicated than single Al-NP, the fundamental questions are, e.g., how does the reaction dynamics proceed, what reaction pathways do assembled Al-NPs go through, whether do aggregated Al-NPs sustain the oxidation reaction or cease to react? Figure 1.5 Photographs of (a) agglomerated Ultra Fine Al-NPs (UFAP) and (b) unagglomerated UFAP [34] 8 Figure 1.6 Comparison of Al-NPs and Al/CuO thermite mixture before and after exposure to the camera flash. (a–c) Optical and scanning electron microscope (SEM) images of Al-NPs before and after the exposure. The original spherical Al-NPs break up into smaller clusters after burning. (d–f) Optical and SEM images of Al–CuO NPs before and after the exposure. The products of the thermite reaction agglomerated into much larger, micron sized particles. [32] Price et al. [38] proposed that slow burning of these agglomerates prevent continuous combustion. Since the agglomeration also results in uneven size and composition distributions of the mixed fuel and oxidizer particles, some researchers believe that aggregated Al-NPs are detrimental to the combustion performance through, e.g., incomplete Al combustion [39], instabilities [40], and slag formation [16]. Various nanomaterial architectures have been explored to achieve a better combustion performance by enhancing the homogeneous distribution of oxidizer and fuel [41, 42]. Table 1.1 summarizes the burning behavior of different aluminum and oxidizer allotments. This set of experiments justifies the notion that mixture of aluminum fuel and oxidizer is too complicated to analyze the results, since not each aspect of measurements 9 shows distinguishable influence on the oxidation dynamics. Table 1.1 Summary of measurements on the Al–CuO and Al–MoO 3 thermites [42]. Linear burning rate (m/s) Mass burning rate (kg/s) Energetic mass burning rate (kg/s) m avg per run (g) % Mass Al 2 O 3 T (k) Nano Al/nano CuO 980 3.8 3.1 0.31 17 2390 Micron Al/nano CuO 660 4.8 4.7 0.58 1.0 2540 Nano Al/micron CuO 200 1.3 1.1 0.53 17 2340 Micron Al/micron CuO 180 2.0 2.0 0.89 1.0 2550 Diluted micron Al/nano CuO 230 1.7 1.4 0.58 17 - Diluted micron Al/micron CuO 86 0.69 0.57 0.64 17 - Nano Al/nano MoO 3 680 2.0 1.4 0.23 26 - Micron Al/nano MoO 3 360 1.5 1.5 0.34 1.7 - Nano Al/micron MoO 3 150 0.45 0.33 0.24 26 - Micron Al/Micron MoO 3 47 0.52 0.51 0.89 1.7 - Among various nanomaterial architectures, mesoporous structures [43, 44] and nanorods based on self-assembly approaches [45, 46] were used instead of spherical oxidizers to achieve higher contact areas with fuel nanoparticles around them and lower resistance for overall diffusional processes. Figure 1.7 shows self-assembled nanoenergetic composite by Subramaniam et al. [46]. It has been shown that the burn rate of the self-assembled composite is significantly higher than that of composite by simple mixing, see Table 1.2. Table 1.2 Burn rate of various Al/CuO composites [46]. Composites Burn rate, m/s CuO nanoparticles mixed with Al nanoparticles 550-780 CuO nanorods mixed with Al nanoparticles 1500-1800 CuO nanorods self-assembled with Al nanoparticles 1800-2200 10 Figure 1.7 Self-assembled nanoenergetic composite. (a) Schematic of self-assembly of Al-NPs around CuO nanorods using poly(4-viny pyridine). (b) TEM image of CuO nanorods after calcination at 450 o C for 3 hrs. (c) TEM image of Al-NPs assembly around CuO nanorods [46]. The rod shape of nano oxidizer can enlarge the overall contact area during the oxidation reaction. Alternatively, the same goal may be achieved using non-spherical nanofuels. In particular, the technology of developing and synthesizing non-particle shaped Al nanostructures are already emerging for Al nanorods (Al-NRs) and Al nanodisks, including electrodeposition [47-49], vapor deposition [50], electromigration [51, 52], and colloidal lithography [53]. Figure 1.8 shows the morphology of Al-NRs 11 under scanning electron microscope (SEM) in top and side views, which shows the distinguished aspect ratio of the aluminum nanostructure. Figure 1.8 SEM image of Al nanorods synthesized in a commercial membrane (a) top view and (b) side/titled view [54]. Much higher resolution microscopy shows that the coating of TiO 2 can cover the outside of the Al-NR, since the structure is for nanobattery purposes [54]. Between the coating of TiO 2 and Al-NR, however, a 4 nm thick Al 2 O 3 layer is intrinsically formed due to the exposure of the Al-NR to air, as in other Al nanostructures. This oxide layer not only provides extra protection to the Al nanorods from reactions with the electrolyte in a battery but also facilitates the growth and adhesion of the TiO 2 layer on the rods. 12 Figure 1.9 HRTEM images of Al nanorods with deposited TiO 2 . (a) A homogenous TiO 2 layer (dark color) is covering the whole surface of the nanorod tip. (b) The TiO 2 layer is also covering the narrow spacing in between the nanorods. (c) There is a 4 nm thick Al 2 O 3 layer between the TiO 2 layer and the nanorod [54]. Figure 1.10 shows SEM images of Al nanostructure products. It has been shown that commercial Al-NRs have more uniform size distribution than Al-NPs. Al-NRs would be very good candidates for studying the oxidation dynamics of Al nanostructures, since they eliminate the size distribution among assembled Al-NRs. Figure 1.10 Typical SEM images of the Al products: (a) commercial Al powders (b) as- deposited Al nanorods [50]. 13 Besides Al-NRs, Al nanodisks also have been introduced to the world of Al nanostructures. Because of the distinguished aspect ratio, D/h (D is diameter of the nanodisk, h is the height of the naondisk), this metal nanostructure has shown significant sensitivity for the localized surface plasmon resonances (LSPR). As is other Al- nanostructures, there is an oxide coating on the surface of the Al nanodisks as well. As shown in Figure 1.12, the thickness of the native oxide layer is around 3 nm, where the thickness increases with time following a logarithmic law. Figure 1.11 SEM micrographs of Al nanodisks with (a) a mean diameter D = 61 nm and height h = 20 nm and (b) D = 115 nm and h = 20 nm. (c) Atomic force microscopy (AFM) scan of a single D = 288 nm Al nanodisks. The 3D image and the line scan reveal a rather rough surface with pronounced grainy microstructure [53]. 14 Figure 1.12 (a) Angle resolved X-ray photoelectron spectroscopy (XPS) analysis of 20 nm Al thin films, evaporated under the same conditions as in the nanodisk fabrication. The native oxide layer thickness is found to be of the order of 2.5-3 nm and forms within a few hours after air exposure. The red solid line represents a direct logarithmic fit to the XPS data. (b) Schematic depiction of the oxidation process in the supported Al nanodisks used in the experiment. It is assumed that a homogenous oxide layer is formed around the particle surface that are exposed to ambient environment, while no oxidation takes places at the glass-Al interface. The oxide thickness is increasing with time [53]. Having surveyed various Al nanostructures, we can conclude the following characteristics of Al fuel as energetic materials: (1) Al nanostructure has significantly higher burning efficiency than Al micron-structures; (2) among various Al nanostructures, Al-NPs have the most variation in the size distribution and oxide shell thickness; (3) Al-NRs have better distinguished aspect ratios to enhance the oxidation activities between oxidizer and fuels; and (4) it is inevitable to have an approximately 3 nm thick oxide coating layer on the surface of Al nanostructures in no matter what form. 1.1.2 High explosives High explosives (HE) have been considered as a type of energetic materials for a long time. In 2012 alone, the consumption of explosives was 3.38 million metric tons (Mt) for a broad range of applications including mining, demolition, automobile airbag and military applications [55]. Despite the extensive production and practical use for so many years, the microscopic properties during the detonation of those explosives have yet to be unraveled [56]. Many researches have invested the reaction dynamics in the 15 process of HE detonation. Reed et al. [56] studied a transient semimetallic layer in detonating nitromethane (CH 3 NO 2 ). They discovered that the wide-bandgap insulator nitromethane undergoes chemical decomposition and a transformation into a semimetallic state for a limited distance behind the detonation front. They found that this transformation is associated with the production of charged decomposition species and provides a mechanism to explain recent experimental observations. Wu et al. [57] studied the catalytic behavior of dens hot water produced from the detonation of high explosive pentaerythritol tetranitrate (PETN) by first principle atomistic simulations. They discovered that H 2 O (source), H (reducer) and OH (oxidizer) act in concert to transport oxygen between reaction centers. Their finding suggests that water may catalyze reactions in other explosives. Research by Pichot et al. [58] has found an unconventional synthesis route for building nanodiamond particles through firing trinitrotoluene (TNT) and hexogen (RDX) mixtures with a negative oxygen balance in a non-oxidative environment. They provided new understanding of nanodiamond formation-mechanisms: discontinuity of the explosive at the nanoscale level plays the key role in modifying the size of the diamond particle. Figure 1.13 shows the difference of nanodiamonds synthesized from different samples of explosives. 16 Figure 1.13 Nanodiamonds synthesized from nanostructured and microstructured explosives. (a–d) TEM micrographs showing nanodiamonds obtained from microstructured (a, b) and nanostructured (c, d) explosive charges. Inset in (a and c) Electronic diffraction patterns showing the diffraction ring of the (111), (220), and (311) diamond planes (d hkl = 0.207 nm, 0.125 nm, 0.107 nm). (e) Size distributions and cumulative frequencies for both samples [58]. Detonation of these HEs involves complex interplay between mechanical shock loading, thermo-mechanical response, and induced chemistry [59]. Brenner et al. [60] have developed an A-B model to describe the chemically sustained shock wave with properties that are consistent with experimental results and the classic continuum theory of planar detonation. Their results demonstrated that simulations provide probe to study the properties between continuum shock waves and atomic scale chemistry. 17 Figure 1.14 Snapshot of a chemically sustained shock wave at 15 ps initiated by a four- layer flyer plate with an impact velocity of 6 km/s. The two types of atoms are depicted in black and white. The shock front is propagating from left to right [60]. The process of detonation starts when the leading shock front compresses and heats up the unreacted material. For sufficiently strong shock, the high temperature and pressure trigger the chemical decomposition of the HE. Once initiated, the HE spontaneously reacts to release energy in supersonic detonation waves. Immediately behind the detonation wave front, the density and pressure increase, in the form of von Neumann spike, then decrease as the chemical reactions progress and the material expands into the reaction zone. The length of this reaction zone, as well the associated reaction time, is a fundamental quantity that rules out various important properties of HEs such as sensitivity [61]. The reaction time in turn is essentially determined by what intermediated reaction products are produced, i.e., the reaction pathways [61-63]. It has been suggested that the slow formation time of certain products play an essential role in determining the reaction time of some HEs [61, 64, 65]. However, the molecular processes that determine the reaction time are largely unknown [61]. For instance, Figure 1.15 shows the largest cluster formed after cooling of a sudden expansion simulation. The mechanism of this large cluster formation still remains elusive. 18 Figure 1.15 Final structure of the largest carbon cluster formed in 10-fold expanded TATB (d = 0.188 g/cm 3 ) after cooling down from T = 3,000 to 500 K. The configuration was taken from the last step of the sudden expansion simulation [65]. Cyclotrimethylenetrinitramine (C 3 H 6 N 6 O 6 or RDX) is one of the most powerful HEs, yet with good stability under ambient conditions. Consequently, this archetypal HE has been studied extensively both experimentally [66-69] and theoretically [70, 71]. For example, experimentally Raman spectroscopy was used to study shock wave-induced phase transition in RDX single crystals [67], and theoretically large scale molecular dynamics simulation has been performed on RDX crystal to reveal the transition from a diffuse shock front with well-ordered molecular dipoles behind it to a disordered dipole distribution behind a sharp front at the initial shock stage [71]. 19 Figure 1.16 Potential energy profile for (a) N-N bond cleavage and (b) the concerted HONO elimination from RDX and the following decomposition of the radical [72]. In particular, the reaction pathways of a single RDX molecule have been mapped out in detail [72] as shown in Figure 1.16. The final products of RDX detonation include carbon monoxide (CO), nitrogen (N 2 ), and steam (H 2 O) [73, 74]. The reaction time for RDX detonation has been inferred from electrical conductivity measurements to be on the order of ns [75-77]. Apart from the overall reaction time, however, little is known about the reaction pathways within ns. Thus, a serious knowledge gap exists on sub-ns reaction dynamics of RDX under high-pressure, high-temperature conditions [66, 78, 79]. 20 Figure 1.17 Conductivity profiles in two HMX charges of different grain sizes. Conductivity σ in S/cm, S (Siemens) is 1/Ω. Coordinate x is measured from the detonation front. The records are terminated at 0.8 mm (shot 151, fine-grained) and at 1.8 mm (shot 150, coarse-grained) because of slit shortening. The definition of the peak width Δ is shown in the upper right [77]. Figure 1.18 Experimental configuration used for launching and probing shock waves in pure aluminum films deposited onto glass. Part of the clipped pump pulse (~ 270 ps 21 duration, centered on 800 nm, spot diameter ~20 µm, rise ~12 ps) is absorbed by the aluminum and consequent plasma expansion drives a shock toward the free Al surface. The time history of the motion of the free surface (up to 250 ps) is obtained in one shot by interfering two probe pulses (~ 270 ps, 800 nm, ~200 µm, separated in time by 10 ps, pulse durations not to scale) in a grating spectrometer. Measurements are made separately at the free surface of the two steps (positions 1 and 2). Thicknesses h 1 and h 2 are, respectively, ~ 0.72 µm and 1.44 µm (not to scale with focal diameters) [80]. Experimental research on explosion (e.g. shock-induced phase transition and shock-induced chemistry) has been a challenge [81, 82]. There is no study on the detonation mechanism of HEs to systematically explain or illustrate the reaction dynamics, especially no direct data is available for the reaction time and reaction pathways. Not only because the reaction products of high explosives are too complex, but also because the reaction normally undergoes extreme high pressures or temperatures, the reaction dynamics is extremely difficult to measure in reality. Due to the advancement of experiments, reactions under shock can now be observed with extremely fine spatial and temporal resolutions [80, 83, 84], as shown in Figure 1.18. We can approximately infer the reaction time according to the value of detonation velocity from experiments. The characteristic of the detonation of HEs can be generalized as: (1) detonation undergoes very high temperature and high pressure exothermic reactions; (2) self-sustained detonation speed is higher than the sonic speed in HE material; and (3) detonation involving complicated chemical reactions can form various reaction products. In this dissertation, we focus on cyclotrimethylenetrinitramine, normally called as 22 RDX, which is an initialism for Research Department eXplosive [85]. Its structural formula is hexahydro-1,3,5-trinitro-1,3,5-triazine or (CH 2 -N-NO 2 ) 3 . RDX has major advantages of possessing greater explosive power than TNT, and requires no additional raw materials for its manufacturing [86]. At density of 1.76 g/cm 3 , the detonation velocity of RDX is 8.75 km/s [87]. At room temperature, it is very stable as a colorless solid. It starts to decompose at about 170 °C and melts at 204 °C. It burns rather than explodes and detonates only with a detonator, being unaffected even by small arms fire. (This is one of the properties that make it a useful military explosive.) It is less sensitive than pentaerythritol tetranitrate (PETN). RDX when exploded in air has about 1.5 times the explosive power of TNT per unit weight and about 2.0 times per unit volume [88, 89]. 1.2 Molecular dynamics simulation of energetic materials For studying oxidation/combustion processes, atomistic simulation is an ideal tool, since it can acquire the atomic/electronic behavior based on the physics modeling for well-defined systems, which can provide prediction and calculations of the properties of materials under extremely high temperatures and pressures. Simulation methods such as finite element method (FEM) are too approximate, and accordingly will lose the predictability [90]. On the other hand, ab initio quantum mechanical methods require too sophisticated computations, [91] which only can be applied to small systems. Molecular dynamics (MD) simulations can provide microscopic data, including positions, velocities, and partial charges of atoms [92]. The conversion of these microscopic data to 23 macroscopic information such as pressure, energy, temperatures, etc. is achieved on the basis of statistical mechanics. Since the size of Al nanostructures normally range between 1 and 100 nanometers, which is 10 to 1,000 times larger than the single atom size. The number of atoms inside nanostructures (in terms of dimensions that we intend to simulate) will not exceed a billion (10 9 ), which is well in the calculation capability of MD simulations. For studying detonation dynamics, the time resolution of reaction is normally within picosecond. [93, 94] MD simulation serves on atomistic modeling well under the extreme condition for duration of picosecond detonation dynamics. Given the time step of MD (on the order of femtosecond), it normally requires 100 ps or 10 5 of MD steps. This is well in the range of actual computational time [95]. With all being said, the sufficient precision and practical amount of time of simulating the energetic materials using MD is the most model method to enhance our understanding of reaction dynamics. The accuracy of calculation relies on the interatomic potentials. The MD simulation on reaction dynamics can complement the knowledge about the structural, dynamical, and mechanical processes for energetic materials by contributing reasonable and useful system design. 1.3 Dissertation overview Parallel MD simulation methods, including MD simulation algrithms and interatomic potentails used in MD simulation, are presented in Chapter 2. The scientific contribution of this dissertation is the investigation of reaction dynamics of energetic 24 materials. The discussion covers four topics of different materials: size effect on the oxidation of single aluminum nanoparticle (Chapter 3); oxidation dynamics of aluminum nanoparticle aggregates (Chapter 4); oxidation dynamics of aluminum nanorods (Chapter 5); and shock-induced detonation reaction in high explosives RDX (Chapter 6). To better guide the reading, Figure 1.19 illustrates the organization of this dissertation, which serves as a route map between chapters and methodologies that are used. Figure 1.19 Diagram of dissertation overview 25 Chapter 2. Simulation Method: Molecular Dynamics In molecular dynamics (MD), atoms are represented as single point masses, which interact with each other under the influence of an interatomic potential energy (or force field). An example is the van der Waals potential, which is typically represented by the 6- 12 Lennard-Jones potential. This means that they cannot move into each other below certain distance determined by a repulsive surface of “hard sphere”. In addition, atoms have partial charges that are calculated by quantum-mechanical methods to represent the distribution of electronic charges of those molecules. Bonds between atoms are typically represented as simple harmonic oscillators along with angular constrains. To start a simulation, a range of velocity is chosen to match the Boltzmann distributions for the specific temperature. Then, those velocities are assigned randomly to each atom. The simulation is run according to the Newton’s equations of motion, with the molecules transferring energy and momentum to one another via electrostatic and van der Waals interactions. The structural measurements from MD can be directly compared with experiments such as neutron scattering [96], fluorescence anisotropy measurements [97], nuclear magnetic resonance (NMR) experiments [98]. This allows the validation of MD by taking actual materials properties into account. The MD simulation approach toward the laws of nature to predict/make precise statements about materials is recognized as the third scientific discipline besides experiments and theory. 26 2.1 Molecular dynamics simulation Molecular dynamics [99, 100] is a computer simulation of physical movements of particles under a given interatomic force field. An N-body simulation deals with particle positions as an 3N-element vector, ! r N = (x 1 ,y 1 ,z 1 ,x 2 ,y 2 ,z 2 ,...,x N ,y N ,z N ) and the velocities of particles also as an 3N-element vector, ! v N = (v 1x ,v 1y ,v 1z ,v 2x ,v 2y ,v 2z ,...,v Nx ,v Ny ,v Nz ) . The trajectories of atoms and molecules are predicted by numerically solving the Newton’s equations of motion, where forces between particles and the associated potential energy are defined by molecular mechanics force fields. ! a i (t)= ! F i (t) m i (2.1) where m i is the mass of ith particle. Figure 2.1 (a) Schematic representation of a sample MD box containing N particles; (b) representation of trajectories and velocities of a single particle at different time; (c) representation of velocities and acceleration of a single particle with respect of time. 27 If we set Δ to be a small increment in time, the velocity of given particle can be inferred from the change of position with respect to time, ! v i (t)= ! " r i (t)= d ! r dt ≡ lim Δ→0 ! r i (t+Δ)− ! r i (t) Δ (2.2) The acceleration of that particle is the change of velocity with respect of time, ! a i (t)= ! "" r i (t)= ! " v i (t)= d ! v dt ≡ lim Δ→0 ! v i (t+Δ)− ! v i (t) Δ (2.3) In the velocity-Verlet algorithm, the velocity expression in the acceleration form is given by half the time step, ! a i (t)= lim Δ→0 ! v i (t+ Δ 2 )− ! v i (t− Δ 2 ) Δ = lim Δ→0 ! r i (t+Δ)− ! r i (t) Δ − ! r i (t)− ! r i (t−Δ) Δ Δ = lim Δ→0 ! r i (t+Δ)− 2 ! r i (t)+ ! r i (t−Δ) Δ 2 (2.4) Through Eq. (2.1) to Eq. (2.4), we can relate the positions and velocities of a system of N particles with given interatomic force fields. Thus as seen from Figure 2.1, the spatio-temporal evolution of a physical system represented by that system could be predicted with given initial values. 2.2 Interatomic force fields To predict the states of N particles system in the next steps, the force applying on each particle is essential. In MD, we consider atomic forces that are derived from the total potential energy, which govern the motion of the atoms. The potential energy can be written as a function of a 3N-dimensional position vector. 28 V =V( ! r 1 ,..., ! r N )=V(x 1 ,y 1 ,z 1 ,x 2 ,y 2 ,z 2 ...,x N ,y N ,z N ) (2.5) Force on the i th atom is represented as ! F i =− ∂V( ! r 1 ,..., ! r N ) ∂ ! r i =− ∂V ∂x i , ∂V ∂y i , ∂V ∂z i # $ % & ' ( (2.6) Here the partial derivative of the potential with respect of the x component for position is defined as, ∂V ∂x i = lim h→0 V(x 1 ,y 1 ,z 1 ,...,x i +h,y i ,z i ,...,x N ,y N ,z N )−V(x 1 ,y 1 ,z 1 ,...,x i ,y i ,z i ,...,x N ,y N ,z N ) h (2.7) In general, the total interatomic potential energy normally contains multiple terms. Here, the form of each term represents a certain physics effect. Lennard-Jones potential [101] has been used for over 80 years to describe the interaction between a pair of atoms. It is one of the most basic atomic interactions and consist of both short-range repulsion and long-range attraction between pair of i th and j th atoms as shown in Figure 2.2, where r ij = ! r ij = ! r i − ! r j is the relative distance between the atoms. V(r)= V(r ij ) i<j ∑ = 4ε 0 i< j ∑ σ 0 r ij " # $ $ % & ' ' 12 − σ 0 r ij " # $ $ % & ' ' 6 ) * + + , - . . (2.8) 29 Figure 2.2 Lennar-Jones potential (from http://www.wikipedia.org/) In this dissertation, we have used two kinds of interatomic force fields, a hybrid many-body potential form and reactive force field (ReaxFF) to describe the aluminum nanostructures and high explosives as real as possible, respectively. 2.2.1 Many-body potentials 2.2.1.1 Potential for Al 2 O 3 To simulate aluminum nanostructures, for aluminum oxide (Al 2 O 3 ), we use a many-body interatomic potential developed by Vashishta et al. [102], which incorporates ionic and covalent effects through a combination of two- and three-body terms, where N is the number of atoms. V = V ij (2) (r ij ) i<j N ∑ + V jik (3) (  r ij , i<j<k N ∑  r ik ) (2.9) 30 The two-body potential V (2) in Eq. (2.9) includes the effects of steric repulsion, Coulomb interaction, charge-dipole interaction and van der Waals interaction: V ij (2) (r)= H ij r η ij + Z i Z j r e −r λ − D ij r 4 e −r ξ − w ij r 6 (2.10) In Eq. (2.10), H ij is the strength of steric repulsion between atoms i and j, Z i is the effective charge of the i th atom, D ij represents the strength of charge-dipole interation, w ij is the strength of van der Waals attraction, η ij is the exponent of the steric repulsion, λ and ξ are the screening lengths for the Coulomb and charge-dipole interactions, respectively, and r is the interatomic distance. The three-body potential V (3) in Eq. (2.9) reflects bond-bending and bond- stretching effects: V jik (3) (r ij ,r ik )= R (3) (r ij ,r ik )⋅P (3) (θ jik ) (r ij ,r ik ≤r 0 ) (2.11) R (3) (r ij ,r ik )=B jik exp γ r ij −r 0 + γ r ik −r 0 " # $ $ % & ' ' Θ r 0 −r ij ( ) Θ r 0 −r ik ( ) P (3) (θ jik )= (cosθ jik −cosθ 0 ) 2 1+C jik (cosθ jik −cosθ 0 ) 2 (2.12) In Eq. (2.12), B ijk is the strength of the three-body interaction between atoms i, j and k. γ is constant and Θ(r 0 - r ij ) is the step function. θ jik is the angle formed by vectors  r ij and  r ik , and r 0 is the cutoff distance for the three-body interaction. θ 0 and C jik are constants. The parameters used in our simulations are listed in Table 2.1. Table 2.1 Parameters for two- and three-body part of our potential used for MD simulation of alumina Al O Z i (e) 1.5237 -1.0158 31 r 1s = 5.0 Å r 4s = 3.75 Å r c = 6.0 Å e = 1.602 x 10 -19 C Two-body Al-Al Al-O O-O η ij 7 9 7 H ij (ev Å η ) 12.4236 244.3038 533.3917 D ij (ev Å 4 ) 0 3.4374 1.52775 W ij (ev Å 6 ) 0 0 63.4894 Three-body B ijk (ev) θ 0 (deg) C ijk ξ (Å) r 0 (Å) O-Al-O 12.4844 90 10 1.0 2.90 Al-O-Al 8.1149 109.47 10 1.0 2.90 The Al 2 O 3 interatomic potential has been validated by comparing the MD results for structural and mechanical properties of both crystalline and amorphous alumina with quantum-mechanical (QM) calculation and experimental data [102]. These include the cohesive energy, elastic constants, and melting temperature of the crystalline alumina. For liquid Al 2 O 3 , the calculated neutron and X-ray structure factors are in good agreement with experimental results. 2.2.1.2 Potential for aluminum For aluminum in the core of the nanoparticles and nanorods, we use the embedded atom method (EAM) form of the potential, where E i is the energy of atom i: E= E i i N ∑ (2.13) E i = 1 2 φ(r ij ) j ∑ +F i (ρ i ) (2.14) The first term in Eq. (2.14) describes the pair interaction between atom i and its neighboring atoms j. The second term describes the attractive interaction due to metallic 32 bonding. This is modeled by the placement of a positively charged atom in the electron density due to the free valence sea of electrons created by host atoms. The embedding function F i (ρ i ) is a function of superimposed charge densities ρ i , due to the charge density distribution of neighboring atoms: ρ i = ρ(r ij ) j ∑ (2.15) where ρ(r ij ) is a pair-wise electronic density as a function of the distance r ij between atoms i and j. In the case of Voter-Chen EAM expression, the pair-wise interaction is taken to be a Morse potential:[103] φ(r)=D M {1−exp[α M (r−R M )]} 2 −D M (2.16) 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: ρ(r)=r 6 (e −βr +2 9 e −2βr ) (2.17) where β is an adjustable parameter. In order to implement the EAM in MD with O(N) scaling, both φ(r) and ρ(r) are truncated at a cutoff radius r cut as h(r)=h(r)−h(r cut )+ r cut m 1− r r cut " # $ % & ' m ( ) * * + , - - dh dr " # $ % & ' r=r cut (2.18) where h(r) represents either φ(r) or ρ(r), and m = 20. 2.2.1.3 Hybrid potential for interface 33 To describe the oxidation of Al, a bond-order based interpolation scheme is used to smoothly interpolate the Al 2 O 3 potential by Vashishta et al. and the EAM potential for Al. The total potential energy of the system is written as the sum of the potential energy of all aluminum atoms and that of all oxygen atoms [104]: E tot = E i = E i i(Al) ∑ i ∑ + E i i(O) ∑ (2.19) The potential energy of each Al atom i is a weighted average between Vashishta’s potential energy for alumina and the EAM energy for Al: E i = E Vashishta.i ⋅ f +E EAM,i ⋅(1− f ) = 1 2 V ij (2) j≠i ∑ % & ' ' ( ) * * ⋅ f + 1 2 φ(r i,j ) j≠i ∑ +F ρ(r i,j ) j≠i ∑ % & ' ' ( ) * * % & ' ' ( ) * * ⋅(1− f ) (2.20) where the weight f is a function of the oxidation degree n i o of atom i, f (n i o )= 0 n i o ≤ 0 1 2 sin((n i o − 0.5)⋅π )+1 ( ) 0< n i o <1 1 n i o ≥1 % & ' ' ( ' ' (2.21) The oxidation degree n i o of atom i is given by n i o = Ω(r ik ) k ∑ (2.22) where Ω(r) is a continuous function that counts the degree of oxidation according to the distance of the Al atom under consideration to a neighboring O atom: Ω(r)= 1 r≤R α −D α 1− r−R α +D α 2D α + sin(π (r−R α +D α ) /D α ) 2π R α −D α <r<R α +D α 0 r≥R α +D α % & ' ' ( ' ' (2.23) We use R a = 3.0 Å and D a = 0.5 Å. 34 The interpolation scheme has been validated by comparing the MD results for the structure of several Al x O y clusters with the corresponding quantum mechanical results [105]. We have also checked the structural and dynamical stability of the interface between the Al core and alumina shell, once they are put together to form an aluminum nanostructures [104]. 2.2.2 Reactive force field (ReaxFF) For detonation reaction of high explosives, we used a reactive force field (ReaxFF) [106, 107] to describe the hydrocarbon systems of RDX crystals. ReaxFF potential can describe chemical reactions and dynamic properties for large-scale reactive system, with both bonded and nonbonded interactions. For bonded terms, bond order is applied to take account of the formation and dissociation of the chemical bond. For nonbonded terms, ReaxFF considers shielded van der Waals and Coulomb interactions with a finite range. The total energy of the simulated system is composed of many different energy terms as the following: E=E bond +E lp +E over +E under +E val +E pen +E coa +E tors +E conj +E hbond +E vdWaals +E Coulomb (2.24) where E bond is the bond energy, E lp is the lone pair energy, E over and E under are the over- and under- coordination energies, E val is the valence angle energy, E pen is the penalty energy, E coa is the three-body conjugation energy, E tors is the torsion angle energy, E conj is the four-body conjugation energy, E hbond is the hydrogen bond interaction energy, E vdWaals and E Coulomb are the van der Waals and Coulomb energies, respectively. 35 2.2.2.1 Bond order and bond energy The bond order, BO ij , a principle concept of ReaxFF, is a function of the interatomic distance r ij between a pair of atoms. BOʹ′ ij is calculated from three exponent terms, which are from σ, π and double-π bonds. ! ! " # $ $ % & ' ' ( ) * * + , + ! ! " # $ $ % & ' ' ( ) * * + , + ! ! " # $ $ % & ' ' ( ) * * + , = - + - + - = - bo6 bo5 bo4 bo3 bo2 bo1 exp exp exp p r r p p r r p p r r p O B O B O B O B o ij o ij o ij ij ij ij ij ππ π σ ππ π σ (2.25) where r ij =| ! r ij |=| ! r i − ! r j | is the interatomic distance between atoms i and j, and 6 5 4 , 3 , 2 1 , , , bo bo bo bo bo bo P P P P P P , are fitted parameter on σ, π and double-π bonds. The above defined uncorrected bond order, BOʹ′ ij is used to obtain the uncorrected overcoordination term, Δʹ′ i , which is the degree of deviation of the sum of the uncorrected bond orders around an atomic center from its valency Val i , i i n j ij i Val BO − = Δ ∑ ∈ ) ( ' (2.26) The corrected bond order term, BO ij is then calculated as BO ij σ = BO ij 'σ ⋅ f 1 (Δ i ' ,Δ j ' )⋅ f 4 (Δ i ' ,BO ij ' )⋅ f 5 (Δ j ' ,BO ij ' ) BO ij π = BO ij 'π ⋅ f 1 (Δ i ' ,Δ j ' )⋅ f 1 (Δ i ' ,Δ j ' )⋅ f 4 (Δ i ' ,BO ij ' )⋅ f 5 (Δ j ' ,BO ij ' ) BO ij ππ = BO ij 'ππ ⋅ f 1 (Δ i ' ,Δ j ' )⋅ f 1 (Δ i ' ,Δ j ' )⋅ f 4 (Δ i ' ,BO ij ' )⋅ f 5 (Δ j ' ,BO ij ' ) BO ij = BO ij σ +BO ij π +BO ij ππ (2.27) Here, the functions f 1 , f 2 , f 3 , f 4 , f 5 , are defined as following: 36 f 1 (Δ i ,Δ j )= 1 2 Val i + f 2 (Δ i ' ,Δ j ' ) Val i + f 2 (Δ i ' ,Δ j ' )+ f 3 (Δ i ' ,Δ j ' ) + Val j + f 2 (Δ i ' ,Δ j ' ) Val j + f 2 (Δ i ' ,Δ j ' )+ f 3 (Δ i ' ,Δ j ' ) " # $ $ % & ' ' f 2 (Δ i ,Δ j )= exp(−p boc1 Δ i )+ exp(−p boc1 Δ j ) f 3 (Δ i ,Δ j )=− 1 p boc2 ln 1 2 exp −p boc2 Δ i ( ) + exp −p boc2 Δ j ( ) ) * + , - . / 0 1 2 f 4 (Δ i ,BO ij )= 1 1+ exp(−p boc3 (p boc4 BO ij BO ij −Δ i boc )+ p boc5 ) f 5 (Δ j ,BO ij )= 1 1+ exp(−p boc3 (p boc4 BO ij BO ij −Δ j boc )+ p boc5 ) (2.28) boc i Δ is a correction for atoms with lone electron pairs, ∑ ∈ + − = Δ ) ( boc boc i n j ij i i BO Val (2.29) where n(i) is the set of neighbor atoms of i that are within a single covalent-bond cutoff distance. The valency of atom i is calculated from the corrected bond order term, ∑ ∈ + − = Δ ) (i n j ij i i BO Val (2.30) The contribution of bond energy E bond , which represents the strength of a covalent bond between atomic pairs, includes the σ, π and double-π bond terms, is computed as € E bond = −D e σ ⋅BO ij σ ⋅exp p be1 1− BO ij σ ( ) p be2 ( ) % & ' ( ) * −D e π ⋅BO ij π −D e ππ ⋅BO ij ππ % & ' ( ) * j ∑ i ∑ (2.31) where e be be D p p , , 2 , 1 , are additional bond parameters. 2.2.2.2 Lone pair energy The valence angles around atoms are affected by the lone electron pairs, which in turn influence the presence of overcoordinated and undercoordinated atoms, especially 37 oxygens and nitrogens. The valency term, € Δ i e , is given by the number of outer shell electrons and a summation of the full bond order, € Δ i e =−Val i e + BO ij j∈n(i) ∑ (2.32) The number of lone electron pairs for atom i, € n lp,i , is as follows: n lp,i = int Δ i e 2 " # $ % & ' exp −p lp1 2+Δ i e −2int Δ i e 2 " # $ % & ' " # $ % & ' 2 ) * + + , - . . (2.33) Once the bond order and the number of lone electron pairs are known, the penalty energy function for atoms with lone electron pairs exceeding the optimal number is € E lp = p lp2 Δ i lp 1+exp −75Δ i lp ( ) i ∑ (2.34) where € Δ i lp is the deviation from the optimal number of lone electron pair: € Δ i lp = n lp,opt − n lp,i (2.35) 2.2.2.3 Overcoordination and undercoordination energy For the effect of overcoordination (Δ i > 0), a penalty term is imposed in energy contribution. For undercoordination (Δ i < 0) the energy contribution of resonance of the π-electron between attached undercoordinated atomic centers is taken into account. The coordination difference Δ i decrease when a lone electron pair appears. The corrected € Δ i lpcorr is € Δ i lpcorr =Δ i − Δ i lp 1+ p ovun3 exp p ovun4 Δ j −Δ j lp ( ) (BO ij π + BO ij ππ ) j∈n(i) ∑ ' ( ) ) * + , , (2.36) 38 Energies for the overcoordinated and undercoordinated atoms are calculated in terms of the corrected € Δ i lpcorr , E over = p ovun1 D e σ BO ij j∈n(i) ∑ Δ i lpcorr +Val i Δ i lpcorr 1 1+exp p ovun2 Δ i lpcorr ( ) $ % & & ' ( ) ) i ∑ E under = (−p ovun5 ) 1−exp p ovun6 Δ i lpcorr ( ) 1+exp(−p ovun2 Δ i lpcorr ) × i ∑ 1 1+p ovun7 exp p ovun8 Δ j −Δ j lp ( ) (BO ij π +BO ij ππ ) j∈n(i) ∑ $ % & & ' ( ) ) (2.37) 2.2.2.4 Valence angle energies The energy contribution arising from the deviation of valence angle Θ ijk from the equilibrium value is described by E val . The equilibrium angle Θ 0 is obtained by SBO j = BO jm π +BO jm ππ ( ) m∈n(j) ∑ + 1− exp −BO jm 8 ( ) m∈n(j) ∏ % & ' ' ( ) * * −Δ j angle −p val8 n lp,j ( ) Δ j angle =−Val j angle + BO jm m∈n(i) ∑ SBO2= 0 (SBO≤ 0) SBO p val9 (0<SBO<1) 2− (2−SBO) p val9 (1<SBO< 2) 2 (SBO> 2) - . / / 0 / / Θ 0 (SBO)=π −Θ 0,0 1− exp −p val10 2−SBO2 ( ) % & ( ) { } (2.38) The function E val goes to zero as the bond order within an atomic triplet goes to zero: 39 E val =∑ i ∑ j ∑ k f 7 (BO ij )f 7 (BO jk )f 8 (Δ j )× p val1 − p val1 exp −p val2 Θ 0 SBO j ( ) −Θ ijk ( ) 2 & ' ( ) * + { } f 7 (BO ij )=1−exp −p val3 BO ij p val4 ( ) f 8 (Δ j )= p val5 − p val5 −1 ( ) 2+exp p val6 Δ j angle ( ) 1+exp p val6 Δ j angle ( ) +exp −p val7 Δ j angle ( ) (2.39) When two double bonds share an atom in a triplet, such as allene, an additional penalty term is introduced: E pen = p pen1 f 9 Δ j ( ) exp −p pen2 BO ij −2 ( ) 2 # $ % & ' ( exp −p pen2 BO jk −2 ( ) 2 # $ % & ' ( k ∑ j ∑ i ∑ f 9 Δ j ( ) = 2+exp −p pen3 Δ j ( ) 1+exp −p pen3 Δ j ( ) +exp p pen4 Δ j ( ) (2.40) Three-body conjugation energy, E coa , has the following form: E coa = p coa1 1 1+exp p coa2 Δ j val ( ) k ∑ j ∑ i ∑ ×exp −p coa3 −BO ij + BO im m∈n(i) ∑ & ' ( ( ) * + + 2 , - . . / 0 1 1 exp −p coa3 −BO jk + BO km m∈n(k) ∑ & ' ( ( ) * + + 2 , - . . / 0 1 1 ×exp −p coa4 BO ij −1.5 ( ) 2 , - . / 0 1 exp −p coa4 BO jk −1.5 ( ) 2 , - . / 0 1 (2.41) 2.2.2.5 Torsion angle energies The torsion angle energy E tors is calculated from 4-body dihedral angles. The dihedral angle w ijkl is the angle between the two planes defined by two atom triplets (i, j, k) and (j, k, l). The torsion angle energy E tors decreases as the valence angles approaches π, or as the bond orders goes to zero, that is 40 E tors = 1 2 f 10 (BO ij ,BO jk ,BO kl )sinΘ ijk sinΘ jkl l ∑ k ∑ j ∑ i ∑ × V 1 (1+ cosω ijkl ) +V 2 exp p tor1 2−BO jk π − f 11 (Δ j ,Δ k ) ( ) 2 { } 1− cos2ω ijkl ( ) +V 3 1+ cos3ω ijkl ( ) & ' ( ( ( ( ( ) * + + + + + f 10 (BO ij ,BO jk ,BO kl )= 1− exp −p tor2 BO ij ( ) & ' ) * 1− exp −p tor2 BO jk ( ) & ' ) * 1− exp −p tor2 BO kl ( ) & ' ) * f 11 (Δ j ,Δ k )= 2+ exp −p tor3 Δ j angle +Δ k angle ( ) & ' ) * 1+ exp −p tor3 Δ j angle +Δ k angle ( ) & ' ) * + exp p tor4 Δ j angle +Δ k angle ( ) & ' ) * (2.42) The 4-body conjugation energy E conj is given by E conj = p cot1 f 12 (BO ij ,BO jk ,BO kl ) 1+ cos 2 ω ijkl −1 ( ) sinΘ ijk sinΘ jkl # $ % & l ∑ k ∑ j ∑ i ∑ f 12 (BO ij ,BO jk ,BO kl )= exp −p cot2 BO ij −1.5 ( ) 2 # $ ( % & ) exp −p cot2 BO jk −1.5 ( ) 2 # $ ( % & ) exp −p cot2 BO kl −1.5 ( ) 2 # $ % & (2.43) 2.2.2.6 Hydrogen bond energy The hydrogen bond energy E hbond is in a form of the bond order and ijk Θ of atomic triplets. E hbond is thus calculated from a triplet i-j-k, where atom j is hydrogen, atom i is chosen from j’s neighbor, and atom k is selected within the finite cutoff radius (~10 Å) from atom i: E hbond = p hb1 1−exp p hb2 BO ij ( ) " # $ % exp p hb3 r hb o r jk + r jk r hb o −2 & ' ( ( ) * + + " # , , $ % - - sin 8 Θ ijk 2 & ' ( ) * + k ∑ j ∑ i ∑ (2.44) 2.2.2.7 Van der Waals and Coulomb energies 41 In ReaxFF, van der Waals and Coulomb forces are calculated among all pairs of atoms to account for Pauli exclusion principle for short interatomic distances and dispersion for long interatomic distances. For van der Waals interaction, a shielded interaction is included in a distance-corrected Morse-potential as follows: E vdWaals = D ij Tap(r ij ) exp α ij 1− f 13 (r ij ) r vdW " # $ % & ' ( ) * + , -−2exp 1 2 α ij 1− f 13 (r ij ) r vdW " # $ % & ' ( ) * + , - . / 0 1 0 2 3 0 4 0 j ∑ i ∑ f 13 (r ij )= r ij p vdW1 + 1 γ w " # $ % & ' p vdW1 ( ) * * + , - - 1/p vdW1 (2.45) Here, Tap has the following definition: Tap(r ij )= 20 R cut 7 r ij 7 - 70 R cut 6 r ij 6 + 84 R cut 5 r ij 5 - 35 R cut 4 r ij 4 +r ij 0 (2.46) where R cut is the non-bonded cutoff radius. The Coulomb potential is expressed as: E Coulomb = 1 2 q i q j j ∑ i ∑ Tap(r ij ) r ij 3 + 1 γ ij " # $ $ % & ' ' 3 ( ) * * + , - - 1/3 (1−δ ij ) (2.47) where γ ij is a parameter for the smeared Coulombic function[108]. In ReaxFF, atomic charges are updated dynamically during force-field optimization procedure using the electronegativity equilibration method (EEM) [108]. The validation of ReaxFF potential on hydrocarbon system, especially on RDX system has been shown good agreement with quantum mechanic calculations [109, 110]. 42 2.3 Simulation algorithm 2.3.1 Numerical integrator In order to realize the dynamical simulation of computing atomic positions as a function of time, the velocity-Verlet algorithm is applied in this dissertation tasking to find € ! r i (t +Δ)and € ! v i (t +Δ), positions and velocities of atoms in the system at given time after an initial condition assigned to all atoms. In this algorithm, the velocity at half the time step is given by € ! v i (t + Δ 2 ) = ! v i (t) + ! a i (t) Δ 2 (2.48) The position at time € t +Δ is € ! r i (t +Δ) = ! r i (t) + ! v i (t + Δ 2 )Δ (2.49) From the updated positions € { ! r i (t +Δ)}, the corresponding forces € ! F i (t +Δ)and accelerations € ! a i (t +Δ) can be computed. The velocity at time t t Δ + is then given by € ! v i (t +Δ) = ! v i (t + Δ 2 ) + Δ 2 ! a i (t +Δ) (2.50) Here, given € ! r i (t), ! v i (t) ( ), velocity-Verlet algorithm predicts € ! r i (t +Δ), ! v i (t +Δ) ( ). 2.3.2 Linked-list cell decomposition As previously mentioned in the interatomic potential part, even for the most simple Lennard-Jones potential, MD simulation of a N atoms system will always have at least the computational complexity of O(N 2 ) due to the doubly nested loops for calculating forces between all pair of atoms. 43 In fact, with a finite cut-off length r c , an atom interacts with only a limited amount of atoms. So with precondition that potential form is shifted to eliminate the force discontinuity at the cut-off, implementation of a linked-list cell algorithm reduces the computational complexity of entire interaction to O(N) operations. The whole system is divided into cells with edges α c r , which are slightly larger than the cutoff length: α α α c c L L r / = , (2.51) ! " c c r L L / α α = , (2.52) Figure 2.3 Cell decomposition of a system for linked list algorithm. where α denotes Cartesian coordinates index x, y, z, and α L represents system length. The edge lengths of each cell r cα must be larger than r c . Atoms in each cell will only interact with atoms in the same cell and other 26 neighboring cells (as shown in Figure 2.3 for a 2D schematic). The central part of the linked-list algorithm is to setup the list of atoms in each decomposed cell as shown in Figure 2.4 [111]. Each cell is assigned a scalar ID c in the 44 range of 0 ~ N C -1 (N C is the total number of cells in the system). It requires two arrays as in form of data structures. The 1 st array is denoted as “head[]” in Figure 2.4 with size of N C elements, and each element stores the first atom ID in each cell. The 2 nd array is denoted as “lscl[]” in Figure 2.4 with size equal to the number of atoms in the system. The linked-list array “lscl[]” holds the atom ID to which the atom points. Figure 2.4 Data structure for linked list algorithm, where “E” denotes empty. 2.3.3 Periodic boundary conditions To understand the importance of boundary conditions for atomistic simulations, let us look at one mole of substance, which contains Avogadro’s number, N A =6.022×10 23 molecules. This is such a big number that even billion-atom simulations on massive parallel computers cannot reach. Correspondingly, statistical 45 behaviors of materials cannot be sufficiently represented in MD simulations, unless more number of atoms is calibrated through proper boundary conditions and without causing excess of the data. In this dissertation, since we are not focusing on surface/interface properties of materials, the usage of Periodic Boundary Condition (PBC) can allow atoms in the boundary layer to adjust their positions in response to the motion of inner atoms ([112]). Hence, there are no boundaries and no artificial surface effects when PBC is used. PBC can be applied along any direction of the system cell according to the need of simulation. In principle, atoms interact with other atoms (even with their own replicas except for themselves) in all periodic replicas. In practical, we make the simulation box larger enough so that atoms interact with only the nearest image of other atoms by applying the minimum image convention. This convention can avoid the ambiguity so that atoms interact only within an imaginary box L x x L y x L z . Therefore, − L x 2 ≤ x ij ≤ L x 2 − L y 2 ≤ y ij ≤ L y 2 − L z 2 ≤ z ij ≤ L z 2 # $ % % % & % % % (2.53) 46 Figure 2.5 Periodically repeated images of an original simulation box (shaded) in 2 dimensions. Minimum image convention applied on primary atom i with imaginary atom j crossing the boundary. 2.3.4 Parallelization Some of the system has the region of interest in the range of over several nanometers. The atomic interaction between the surrounding environment and the specific region plays very important role during the dynamic reaction. The artifact that brought by PBC cannot reflect the nature of physics inside materials. The large-scale simulation involving multi-millions to billions of atoms is critical to simulate this kind of system. However, due to the limitation of contemporary computability, the execution of such massive operation is not achievable on a single computer. Parallel computing on CPUs, which breaks down physical system into domains of equal volume and maps the computation inside and among the domains to parallel computers accordingly, as spatial decomposition [113, 114]. In an MD simulation, we 47 use a simple 3D mesh decomposition to arrange the processors that computes physical properties of a domain, since the interatomic force in MD only carrying on within a finite range. It makes the whole idea of parallel computing for MD simulation feasible. In this dissertation, message passing interface (MPI) library carries out the communication between processors. If there are P computing processors, we rationally organize them as P=P x ×P y ×P z , where P x ,P y ,P z is number of processors along the x, y, z direction, respectively. Each processor is given a unique processor ID, p, ranged from 0 to P-1 to compute a cell, as shown in Figure 2.6. Atom i at position ! r i = (x i ,y i ,z i ) will be assigned to processor p as p= p x P y P z + p y P z + p z (2.54) p α = r iα /P α L α ! " # $ (α= x,y,z) (2.55) Figure 2.6 Schematic view of assigning process ID for domains. 48 where Lα represents the length of the edge length of each domain, and pα indicates the vector index in the x, y and z directions. Each processor p deals with atoms types, positions, velocities and partial charges within the domain. In order to perform the O(N) computation in a N atoms MD simulation, while also keep the integrity of the simulation, it is important to calculate the all interactions between domains, even for the atoms that are not in the primary domain. To obtain the information at the boundaries of each domain, extended-subdomain (including a corner, an edge, and a face sharing neighbor domains by a cut-off-length-thick layer as shown in Figure 2.7) is considered through caching by inter-processor communications. This inter-processor communication involves updating information from the 26 neighboring domain processors. In atom caching, coherence is maintained by copying the latest neighbor surface atoms every time before atomic accelerations are computed. After the atomic positions are updated according to the velocity Verlet algorithm, some atoms from primary domain have moved out of the domain boundary. Those updated atoms have to be moved to proper processors with coordinated migration schema. 49 Figure 2.7 A schematic of a projection view of spatial decomposition of 3D system implemented in parallel MD program. (a) Atoms attributes in the extended-cell are cached from 26 neighboring domains. (b) Atoms attributes move out the primary domain to other proper processors. 2.4 Performance Parallel computing efficiency is mainly decided by the computing hardware (e.g., computer architecture). Computing hardware normally has three types of memory organizations: shared-memory system, distributed-memory system or hybrid distributed- shared memory system. In this dissertation, we use a distributed-memory system: each processor has access to its own local memory; buffer data is shared by inter-processor network. Message Passing Interface (MPI) is used widely as a language independent communication protocol for this type of shared-memory system. p 50 In terms of scalability of a parallel MD simulation, we define the efficiency of a parallel program running on P processors to solve a problem of size W in T(W,P) execution time. The speed of the program is given by S(W,P)= W T(W,P) . Speedup of the program S p is therefore S p = S(W,P) S(W,1) , and the efficiency E can be defined asE P = S P p . In a constant problem-size (or strong) scaling test, the speedup is S P = S(W,P) S(W,1) = W T(W,P) W T(W,1) = T(W,1) T(W,P) , the efficiency isE P = S P P = T(W,1) P•T(W,P) . Whereas in an isogranular (or weak) scaling test, the problem size W P scales linearly with the number of processors P: W P =P ⋅w , where the granularity w is constant (the work per processor). The speedup isS P = S(W P ,P) S(w,1) = W P T(W P ,P) w T(w,1) = P ⋅w T(P ⋅w,P) w T(w,1) = P ⋅T(w,1) T(P ⋅w,P) , the efficiency is E P = S P P = T(w,1) T(P ⋅w,P) . The following figures are the weak scalability tests done on finite range potential parallel MD algorithm for 2,044,416P atoms silica systems (Figure 2.8) and parallel ReaxFF MD algorithm for 16,128P atoms RDX systems (Figure 2.9) [115]. They show benchmark results using IBM BlueGene/L (open symbols) and BlueGene/P (solid symbols) computers. The wall-clock time (circles) and inter-processor communication time (squares) per MD step are shown. 51 Figure 2.8 Results of the isogranular scaling (2,044,416 atoms per processor) using finite- range MD simulation on silica systems. The weak-scaling parallel efficiency is 0.975 on 131,072 BlueGene/P processors. The weak-scaling parallel efficiency on 212,992 BlueGene/L processors is 0.985 based on the speedup over 4,096 processors [115]. 52 Figure 2.9 Execution time of ReaxFF MD on energetic materials (RDX) per time step as a function of the number of processors P using our parallel ReaxFF MD algorithm for energetic materials (RDX) on BlueGene/L (squares) and BlueGene/P (circles) computers, where the number of atoms is N = 16,128P. The efficiency for BlueGene/L is 0.996, which is close to the ideal efficiency, 1 [115]. 53 Chapter 3. Size Effect on the Oxidation of Single Aluminum Nanoparticle This chapter primarily discusses the size effect on the oxidation of aluminum nanoparticles (Al-NPs). Background and motivation of the simulation are discussed in section 3.1. In section 3.2, we describe the system setup of MD simulations to study single Al-NP with different sizes. The size effect of Al-NPs within the studied range is discussed in section 3.3 by treating the simulated system as an entirety (to calculate the temperature and radial distribution) and taking different parts (core, shell and outside clusters) into consideration. Section 3.4 concludes the size effects on the oxidation of Al- NPs. Within the cases that studied in this chapter, the smaller Al-NPs oxide faster in terms of shorter reaction time and shorter reaction delay. However, all size Al-NPs have the same temperature for shell Al atoms ejection onset. 3.1. Background and motivation Combustion of aluminum is a subject of great interest in science, engineering and technology. It has been studied since the last century due to its wide use in the field of solid propellants and explosives. For instance, combustion of micron size Al particles were studied by Drezin [116] and Jordan et al.[117] using different experimental approaches. The effect of temperature on oxidation[24, 27, 118-120] and the effect of oxide shell on the surface[121-128] have been studied extensively, especially for 54 nanometer-size Al particles. Many different mechanisms have been proposed to explain the combustion behavior of Al nanoparticles (Al-NPs).[129-132] By adding Al-NPs into micron-size particles, numerous experiments have observed the enhanced oxidation reactivity in terms of burning rate,[133, 134] flame speed,[28, 135] activation energy,[136] ignition sensitivity,[135, 137-139] combustion velocity,[137] and agglomeration.[118] Several investigators have specifically investigated the size effect on the combustion of Al-NPs.[25, 140] Most of them concluded that the smaller the size of Al-NPs is, the higher the reactivity becomes. However, the effect of particle size on reactivity is still controversial. For example, Gan et al.[141] have conducted experiments on fuel droplets with nano (~80 nm) and micron (~5 mm and ~25 mm) size Al particles. They found that the combustion is longer and less complete for large agglomerate of nano-suspensions due to the formation of oxide shell on the Al-NPs surface. In a review article, Yetter et al.[17] stated that for smaller Al particles, the energy loss per unit volume due to the presence of the same oxide layer thickness is significantly larger. Besides the above-mentioned experimental studies, there have been numerous simulations of oxidation dynamics of Al-NPs. Campbell et al.[142, 143] studied the oxidation of aluminum nanoclusters with a parallel molecular dynamics (MD) approach based on dynamic charge transfer among atoms in both micorcanonical and canonical ensembles under atomic and molecular oxygen environments. Alavi et al.[144] have simulated the oxidation of Al-NPs using MD with Streitz-Mintmire electrostatic plus (ES+) potential.[145, 146] Puri et al.[147] performed MD simulation to study the 55 thermo-mechanical behavior of nano Al particles coated with crystalline and amorphous oxide layers. Perron et al.[148] have studied the oxidation of multi-grain nanocrystalline (mean grain size = 5 nm) Al surfaces in the temperature range 300 – 750 K using variable charge molecular dynamics simulations. Structures of g- and amorphous Al 2 O 3 have been investigated using MD and density function theory (DFT) by Gutierrez et al.[149, 150] Vashishta et al.[151] have developed an interaction potential consisting of two- and three-body terms for alumina to simulate the amorphous and liquid phases. Combining Vashishta’s potential[151] and embedded atom method (EAM) potential [152], Wang et al.[153] have studied atomistic mechanisms of oxidation in a laser flash heated core (aluminum) – shell (alumina) nanoparticle. In another paper by Wang et al.[154], they studied effects of the crystalline and amorphous structures of alumina shells on the dynamics of oxidation of an Al-NP. Despite the extensive experimental and simulation research efforts mentioned above, several key questions remain unanswered. Most important of them are: the collective reaction behaviors among particles; and spatially and temporally resolved oxidation dynamics of individual particles at atomic resolution. It is essential to delineate these issues unambiguously and study each factor independently. Regarding the collective reaction, Shekhar et al.[155] have investigated the mechano-chemial reaction in the oxidation of three adjacent Al-NPs by initiating the oxidation in the middle nanoparticle and studying the nature and speed of the oxidation front. However, there is no systematic study on the atomistic mechanisms underlying size-dependent oxidation dynamics of individual Al-NPs. This chapter focuses on the size dependence of oxidation 56 dynamics rate and that of reaction temperature, as well as the role of oxide shell in Al- NPs oxidation. Here, we report the results of large-scale, parallel MD simulations of the size effect on single Al-NP of diameters D = 26, 36 and 46 nm. 3.2 System setup for three single Al-NPs MD simulations The initial setup for the oxidation simulation is as follows (see Table 3.1). The Al- NP diameters are chosen to be D = 26, 36 and 46 nm, with 3 nm thick amorphous alumina shell. The metallic Al core is cut from crystalline Al at 300 K. The amorphous shell is prepared by first heating the crystalline alumina, and then slowly quenching the molten alumina into amorphous alumina at 300 K, and then, from the amorphous alumina, removing the outer and inner concentric parts accordingly to retain an oxide layer of 3 nm thick. After combining the crystalline aluminum core and amorphous alumina shell structure together, the single Al-NP is relaxed and thermalized before it is put into an oxygen environment. Table 3.1 Reactive molecular dynamics simulation setup for aluminum (core) and amorphous alumina (shell) nanoparticles in oxygen environment Particle diameter (nm) Shell thickness (nm) No. of core Al atoms No. of shell Al 2 O 3 atoms No. of environmental oxygen atoms Total No. of atoms in system Simulation box size (nm) 26 3 252,288 440,106 237,264 929,658 52 36 3 851,423 897,290 802,498 2,551,211 72 46 3 1,979,315 1,529,708 1,889,162 5,398,185 92 The following relaxation procedure is performed for the Al-NP assembled after putting the core (aluminum) and shell (amorphous alumina) together. The reactive MD was run on the Al-NP for 30 ps, in microcanonical ensemble, to form and stabilize the aluminum/alumina interface. Figure 3.1(a) displays the radial temperature distribution 57 during preparation of 26 nm Al-NP system. We plot temperature within a 10 Å thick spherical shell as a function of radial distance R. The red curve in Figure 3.1(a) shows the temperature profile at 30 ps, which is the end of interface-formation process. The shell temperature is higher than core and the environmental oxygen because of the shell reconstruction. The Al-NP is then quenched to 5 K (blue curve in Figure 3.1(a)) and the whole system is gradually heated to 300 K (as shown in time sequence of temperature profiles from cyan to green to orange in Figure 3.1(a)). In order to start the oxidation dynamics in the Al-NP, velocity scaling is performed to preheat nanoparticle from 300 K to 1,100 K in steps of 100 K (see Figure 3.1(b)). In each step, the Al-NP is heated at a constant heating rate, followed by thermalization for 0.1 ns during each interval. During the Al-NP heating, the oxygen environment is kept at 300 K as shown in Figure 3.1(b). After the preheating, systems undergo oxidation for 1 ns in microcanonical ensemble. Exactly the same procedure is applied to D = 36 and 46 nm Al-NPs. 58 Figure 3.1 Radial temperature profile of D = 26 nm Al-NP during (a) preparation of aluminum core and alumina shell system, and (b) preheating of the system to 1,100K. The dashed line indicates the interface between the Al-NP and environmental oxygen. The arrows indicate the time sequence of the preparation and preheating. 59 3.3 Analysis of reactive molecular dynamics simulations 3.3.1 Oxidation dynamics of three Al-NPs Figure 3.2 Snapshots of the central slice for D = 26 nm Al-NP at time = 0.0, 0.2, 0.4, 0.6, 0.8, and 1 ns. Figure 3.3 Snapshots of the central slice for D = 36 nm Al-NP at time = 0.0, 0.2, 0.4, 0.6, 0.8, and 1 ns. 60 Figure 3.4 Snapshots of the central slice for D = 46 nm Al-NP at time = 0.0, 0.2, 0.4, 0.6, 0.8, and 1 ns. To visualize the reaction process of the three Al-NPs, Figure 3.2, Figure 3.3 and Figure 3.4 show snapshots of the central slice, respectively, for D = 26, 36 and 46 nm Al- NPs at 0.0, 0.2, 0.4, 0.6, 0.8 and 1.0 ns. Blue color represents environmental oxygen. Green and black colors indicate the oxygen and aluminum atoms, respectively. In the alumina shell, the red color represents the core aluminum in the initial configuration. The oxidation process starts after the initial conditions are set up (core and shell preheated to 1,100 K). As the oxidation begins, due to the molten core expansion at 1,100 K, the size of the Al-NP grows larger compared to its initial size. No shell breaking or shattering of the shell is observed in the oxidation process; only the deformation of Al-NP is observed. From the extent of the core and shell 61 deformation, reaction is apparently more intense for D = 26 nm Al-NP than D = 36 and 46 nm Al-NPs. 3.3.2 Temperature profile of three Al-NP systems Figure 3.5 (a) Temperature of the system (i.e., Al-NP and the environmental oxygen) during the oxidation as a function of time for D = 26, 36 and 46 nm. (b) Radial 62 temperature of the three Al-NPs at time 0.3 ns. Position of the alumina shell is indicated by arrows. To determine the size dependence on oxidation reactivity of Al-NPs, we study the heat release for all three systems by plotting temperature of whole systems as functions of time in Figure 3.5(a). Significant amount of heat is produced as a result of oxidation within the nanoparticles in all three systems. The radial temperature distribution of the three systems at 0.3 ns is shown in Figure 3.5(b). Temperature at the center of the core is higher than the shell temperature for the smallest system (D = 26 nm) as shown in Figure 3.5(b). However, the temperature at the center of the core is lower than the shell temperature at 0.3 ns for larger systems (D = 36 and 46 nm). Comparison between radial temperature distribution and global temperature of all systems (at time less than 0.3 ns) shows that the shell temperature is higher than the average temperature of the system, which indicates that oxidation at core-shell interface dominates the heat release of the whole system at early stage. 63 3.3.3 Local aluminum atoms concentration Figure 3.6 (a) Schematic of a quarter slice of the whole system. Dark blue, grey and yellow circles are aluminum atoms from alumina shell, aluminum atoms from aluminum core, and oxygen atoms, respectively. (b) Local concentration is calculated at a distance R by averaging the concentration within a spherical shell of thickness dR = 10 Å at radius R from the center of the nanoparticle. 64 In order to study the distribution of aluminum atoms inside Al-NPs, we calculate local concentration of aluminum atoms from their two different origins, i.e., those from the alumina shell and from the aluminum core (see Figure 3.6(a)). Local concentration at distance R from the center of Al-NP is calculated by averaging the concentration within a spherical shell of thickness dR = 10 Å. We compute three concentrations as follows: C Al = n Al n total (3.1) C Al shell = n Al shell n total (3.2) C Al core = n Al core n total (3.3) C Al is the local concentration of Al, C Al shell is the local concentration of Al from alumina shell, and C Al core is the local concentration of Al from aluminum core. Heren total is the total number of atoms, n Al is the total number of Al atoms, n Al shell is the number of Al atoms from alumina shell, n Al core is the number of Al atoms from aluminum core, respectively, at distance R from the center of Al-NP within a spherical shell of thickness dR = 10 Å. Figure 3.7 shows C Al as function of radius R for time t = 0, 0.3 and 0.6 ns for D = 26, 36 and 46 nm. As shown in blue lines in Figure 3.7 (a), (b) and (c) for time 0 ns, there are shoulders around R = 10~13, 15~18, and 20~23 nm, respectively. These shoulders reflect the core-shell structure of Al-NPs. Namely, the center part, up to R =10, 15 and 20nm, respectively, for the three systems, is pure Al. The shell part is Al 2 O 3 extends for 3 nm beyond the core. There are no Al atoms beyond the shell radius of the three systems. 65 As time progresses, these shoulders become less pronounced indicating the intermixing of atoms between shell and core. Figure 3.7 Local concentration of Al atoms, C Al Eq. (3.1), at 0 ns, 0.3 ns and 0.6 ns for (a) D = 26 nm, (b) D = 36 nm and (c) D = 46 nm. 66 To distinguish origin of Al atoms, Figure 3.8 and Figure 3.9 show C Al shell and C Al core , respectively. In Figure 3.8 , the movement of shell Al atoms is observed from the radial distribution at different times. At 0 ns, blue lines, the local concentration is 0.4, reflecting the composition of Al 2 O 3 . At 0.3 ns, green lines, the alumina shells have moved outwards and become wider due to the expansion of the molten core. For D = 26 nm, more Al atoms spread inwards, which is not observed in D = 46 nm until 0.6 ns. At 0.6 ns, red lines, more shell Al atoms spread outwards than inwards for D = 26 and 36 nm. At times beyond 0.6 ns, the same behavior is seen in D = 46 nm. In Figure 3.9, at 0 ps, blue lines, the center part is pure (metallic) Al, while there is no core aluminum in alumina shell or oxygen environment. At 0.6 ns, red lines, the core Al spreads outward. The core Al atoms in D = 26 nm moves outwards earlier in time than other systems at the same instant. 67 Figure 3.8 Local concentration of Al atoms from the alumina shell, C Al shell Eq. (3.2), at 0 ns, 0.3 ns and 0.6 ns for (a) D = 26 nm, (b) D = 36 nm and (c) D = 46 nm. 68 Figure 3.9 Local concentration of Al atoms from core, C Al core Eq. (3.3), at 0 ns, 0.3 ns and 0.6 ns for (a) D = 26 nm, (b) D = 36 nm and (c) D = 46 nm. 69 3.3.4 Formation of Al x O y clusters To understand the oxidation in the core-shell region for three different sizes of Al- NPs, we next investigate the chemical composition of oxide clusters using fragment analysis. As the oxidation starts, the aluminum oxide clusters with a variety of aluminum- oxygen ratios have been found both in simulations[143, 156] and experiments.[157-159] Here we take the fragment as the atoms that are covalently bonded. Except the large fragments (those having over a thousand atoms, which correspond to the oxide shell), in the early stage of the reaction, the majority of fragments in all systems are small clusters, having less than six atoms. Figure 3.10, Figure 3.11, and Figure 3.12 show the number of clusters Al 2 O, AlO and AlO 2 in D = 26, 36 and 46 nm, respectively, in the reaction process. We found that, in all three systems, the number of Al 2 O clusters (blue) decreases after an initial increase. The AlO (red) and AlO 2 (green) clusters become more numerous after certain times. From the figures, we observe that the majority of the clusters in all three systems change from Al-rich to oxygen-rich after different time durations. For D = 26 nm, the population of AlO exceeds that of Al 2 O around 0.3 ns, which is earlier than 0.5 ns for D = 36 nm and 0.7 ns for D = 46 nm. 70 Figure 3.10 Number of (a) Al 2 O, (b) AlO and (c) AlO 2 clusters as a function of time in D = 26 nm respectively. 71 Figure 3.11 Number of (a) Al 2 O, (b) AlO and (c) AlO 2 clusters as a function of time in D = 36 nm respectively. 72 Figure 3.12 Number of (a) Al 2 O, (b) AlO and (c) AlO 2 clusters as a function of time in D = 46 nm respectively. 73 3.3.5 Evolution of the shell structure During the oxidation, the shell structure of Al-NPs changes, in width due to core expansion and also in its chemical composition. Figure 3.13(a) shows the outer radius of the Al-NPs. Within 1 ns, D = 26 and 36 nm systems show expansion followed by contraction, whereas 46 nm system only shows expansion. The contraction phase for the D = 46 nm system is beyond 1 ns. The extent of expansion (the increment percentage of the radius from the initial value to the maximum value) of the D = 26, 36 and 46 nm systems are 18.89%, 22,33% and 24.18%, respectively. The oxide shell radius is calculated by averaging distances from the center of the system to all atoms contained within the largest oxide fragment. Figure 3.13(b) shows the ratio between oxygen atoms and aluminum atoms in the largest oxide fragments (i.e., the oxide shell). Initially, the O/Al ratio of 1.5 stands for the amorphous Al 2 O 3 shell. The formation of the shell structure is a dynamic process involving Al atoms coming into the shell from the core and Al atoms leaving the shell into the oxygen environment. The overall tendency of the ratio is declining. At early times, the O/Al ratio in the shell region decreases with time for all the three systems. For D = 26 nm, oxygen incorporation from the environment increases the ratio until it peaks around 1.45 at 0.5 ns and then decreases. A similar behavior is observed for D = 36 and 46 nm, but at later times and with different ratio values. By 1 ns, the ratio for all systems are under 1.4, which means the largest fragment in the system changes from perfect amorphous oxide shell to an Al-rich shell. 74 Figure 3.13 Different shell structure of D = 26, 36 and 46 nm systems are illustrated by (a) the radius of Al-NPs as a function of time, and (b) the ratio of oxygen and aluminum within the shell region. In the process of oxidation, we observe some shell Al atoms eject out into the surrounding oxygen. In order to quantitatively determine the ejection process, after the 75 oxide shell structure has been defined (the largest oxide fragment), we count the number of shell Al atoms that are located beyond the oxide shell radius. Figure 3.14(a) plots the number of ejected shell Al atoms as a function of time. The onset time for Al ejection from the shell are t 0 = 0.18 ns, 0.28 ns and 0.43 ns for D = 26, 36 and 46 nm, respectively. For D = 26 nm, the shell aluminum atoms eject earlier than D = 36 and 46 nm, resulting in a faster heat release. This is reflected as higher global temperature in D = 26 nm than D = 36 and 46 nm after the same amount of time as shown in Figure 3.14(b), which plots the temperature of oxide shell during oxidation as a function of time. Comparison of Figure 3.14 (a) and (b) shows that for all three Al-NP systems the shell temperature at the onset Al ejection time is the same, i.e., ~2900 K. The dashed arrows specify the onset time of Al ejections as obtained from Figure 3.14(a). The black dashed double arrows mark the shell temperature when Al ejections begin. Thus, we conclude that only after the shell temperature reaches 2,900 K, the shell Al atoms will start their ejections into the oxygen environment. At 2,900 K, the alumina shell is in molten state in all three systems. 76 Figure 3.14 (a) Semi-log plot of the number of ejected Al shell atoms as a function of time for D = 26, 36 and 46 nm. (b) Oxide shell temperature during oxidation as a function of time for the three Al-NP systems. Ejection of Al atoms from the shell starts when the shell is in molten state at 2,900 K. 3.3.6 Reaction delay In order to explain the time dependence of reaction for the three systems, we plot temperature of Al-NPs including the surrounding oxygen at different time in Figure 77 3.15(a) and the rate of change in temperature (dT system /dt) for the three systems in Figure 3.15(b). Comparing Figure 3.15(a) and (c) (same as Figure 3.14(b)) the temperature of the global systems and the oxide shells have the same increasing trend. This restates the fact that the temperature of the shell is higher than the average temperature of the whole system, which indicates that reaction in the core-shell region dominates the heat release as mentioned in 3.3.1 Oxidation dynamics of three Al-NPs. Figure 3.15(b) shows that the largest rate of change in temperature (dT system /dt) occurs at t 1 = 0.2 ns, 0.32 ns and 0.51 ns for D = 26, 36 and 46 nm, respectively. At these times, the reaction is most intensive. Once passing this peak (i.e., the largest rate of change in temperature (dT system /dt) max ), temperature still increases but at a slower rate. In Figure 3.15(a), the temperature at the culmination point is given for each Al-NP (2,500 K, 2,700 K and 2,800 K for D = 26, 36 and 46 nm, respectively). 78 Figure 3.15 (a) The global temperature, (b) The rate of change in the global temperature (dT system /dt) and (c) Oxide shell temperature as a function of time for D = 26, 36 and 46 nm 79 It is worth noting that the delay of the culmination points from the onset of shell- Al ejection is different for the three Al-NPs systems. The delay for D = 26, 36 and 46nm are (t 1 - t 0 ) = 0.02, 0.04 and 0.09 ns respectively. Namely, smaller system has shorter delay. Table 3.2 shows the detailed comparison of Figure 3.15(a) and Figure 3.15(c). Table 3.2 Relationship between the delay time and the nanoparticle size D Diameter of the Nanoparticle t 0 Time for the onset of Al ejection from the shell Temperature of the shell at the onset Al ejection t 1 Peak time of (dT system /dt) (time derivative of the temperature of the whole system) Temperature of the shell at peak of (dT system /dt) Temperature of whole system (including environmental oxygen) at the Peak of (dT system /dt) t 1 - t 0 Reaction Delay 26 nm 0.18 ns 2,900 K 0.20 ns 3,100 K 2,500 K 0.02 ns 36 nm 0.28 ns 2,900 K 0.32 ns 3,300 K 2,700 K 0.04 ns 46 nm 0.42 ns 2,900 K 0.51 ns 3,500 K 2,800 K 0.09 ns Finally, we would like to point out the possible effects of oxygen density on the reaction dynamics. In this paper, we used a higher oxygen density so that the entire oxidation process can be observed during the simulation time. We expect that the earlier stage of the reaction is less sensitive to the oxygen density because the reaction is occurring inside the nanoparticle at the metal core and oxide shell. However, the reaction at the later stage involves the environmental oxygen, thus it should depend on the oxygen density.[160] This can be seen in the D = 26 nm Al-NP, where the depletion of environmental oxygen slows down the oxidation reaction. From Figure 3.15(b), the oxygen from the environment for 26 nm system is almost consumed after 0.6 ns, which leads to a decrease of the reaction rate (dT system /dt), as shown in Figure 3.15(c). 80 Chapter 4. Oxidation Dynamics of Aluminum Nanoparticle Aggregates In this chapter, the collective oxidation behaviors of Al nanoparticle aggregates are discussed. Uniformly sized Al-NP aggregates are considered systematically changing the particle size. Section 4.1 provides the background and motivation of studying Al-NPs aggregates. The systematical simulation setup for uniformly distributed Al-NP aggregates is described in section 4.2. Section 4.3 presents the oxidation dynamics of different Al- NP diameters (D = 26, 36 and 46 nm) aggregates from multiple aspects, including the agglomeration degree, consumption of remaining core Al, reaction crossover of small intermediate fragments (Al 2 O, AlO and AlO 2 ), and evolution of large oxide fragments. Section 4.4 concludes the significance of the agglomeration of Al-NPs and the production of small intermediate aluminum-oxide fragments. Comparison of oxidation of different diameters Al-NPs aggregates shows faster reactions of smaller diameter Al-NP aggregates. We find crossovers of the major intermediate products from Al 2 O to AlO 2 then to AlO. Also found is the production of the final product, Al 2 O 3 , directly from the AlO and AlO 2 intermediates. Our simulation suggests an unexpectedly active role of the oxide shell as a nanoreactor. 4.1 Background and motivation While the combustion of Al-NPs are highly useful in the field of rocket fuel when combustion, in reality, Al-NPs do not burn independently from each other. Transmission electron microscope (TEM) images show Al-NPs to be close packed as loosely 81 compacted aggregates or fused into agglomerates [20, 32, 161] as shown in Figure 4.1. Ohkura et al. [32] observed different morphologies of the reaction products from Al-NPs after flash heating under different conditions. Al-NP aggregates were found to form small clusters after oxidation in air, while they formed micron-size agglomerations using CuO as oxidizer. Rashkovskii [33] found a potential relationship between the normal acceleration (g-load) and the size distribution of Al-NP aggregations, i.e., a nonmonotonic dependence of the mean-mass size of agglomerates on the magnitude of normal acceleration. Agglomerated ultrafine Al powders (with size less than 100 nm) and unagglomerated ultrafine Al powders were studied in air to examine the their reactivity [34]. Under normal pressure, lower concentration of bound nitrogen was found in agglomerated powders. Figure 4.1 Transmission electron microscopy (TEM) image of Al nanopowders aggregates [161]. 82 Besides these experimental studies, several simulation studies were conducted to study the properties of Al-NPs. Alavi et al. [144] studied the structure and charge distribution of Al-NPs with different Al-to-O ratios. Puri et al. [147] studied the thermo- mechanical behavior of oxide-coated Al-NPs of different sizes during melting. However, most simulations were restricted to the Al-NP size less than 10 nm, which is smaller than typical particle sizes found in experiments, tens of nanometers (e.g. 30 nm) [162]. Previously, there are reactive molecular dynamics (RMD) simulations on massively parallel computers [163] to study various aspects of oxidation dynamics involving larger Al-NPs, such as different heating methods [164], different oxide-shell structures [165], and size effect [166]. Recently, Shekhar et al. [167] also studied the collective oxidation behavior of a chain of three Al-NPs. However, the oxidation dynamics of Al-NP aggregates has never been studied by reactive MD simulations, primarily due to the enormous computation that is required. Key scientific questions are: Do Al-NPs agglomerate or fragment during oxidation; what is the size effect on the reaction time; and what are the reaction pathways (or intermediate reaction products)? Here, we answer these questions using multimillion-atom reactive MD simulations of the oxidation of close-pack aggregates of uniformly sized Al-NPs. To study the size effect, Al-NP aggregates of different diameters are compared. 4.2 System setup for three Al-NPs aggregates MD simulations We simulate aggregates of uniformly sized spherical Al-NPs of diameter D. Each Al-NP consists of an Al core of diameter Dʹ′ (= D−6 nm) covered by a 3 nm thick 83 amorphous alumina shell [168]. Figure 4.2(a) shows the core-shell structure of a single Al-NP. To study the size effect, three aggregates with different diameters (D = 26, 36 and 46 nm) are considered. In each aggregate, Al-NPs are placed 0.5 nm apart from each other and are stacked along the x direction in 5 layers as labeled in Roman numerals I to V in Figure 4.2(b). In the close-packed ABCBA stacking structure, an odd layer (I, III, or V) contains three Al-NPs, whereas an even layer (II or IV) contains a single Al-NP. For each aggregate, three Al-NPs in layer I (colored red in Figure 4.2(b)) are preheated to 1,100 K [168], while the rest of the Al-NPs are kept at 300 K. The aggregate is then placed in a 300 K oxygen environment within a box of lengths 5.3D × 2.9D × 3D in the x, y and z directions. The closest distance between the boundary of the box and the Al-NPs is 0.5D. Periodic boundary conditions are applied in the y and z directions, whereas the moment mirrors are placed at both ends in the x direction. We refer the three systems as D = 26, 36 and 46 nm Al-NP aggregates. For each aggregate, simulation is performed in the microcanonical ensemble for 2.0 ns. The equations of motion are integrated using the velocity Verlet algorithm with a time step of 1 fs. The D = 26, 36 and 46 nm Al-NP aggregate systems contain 9.3 million, 23.4 million and 48.7 million atoms, respectively, which are distributed among 256, 512 and 768 processors in a parallel computer. 84 Al2O3 shell Core Al D 3 nm (a) (b) Figure 4.2 (a) Core-shell structure of a single Al-NP. An Al core (blue) is covered by a 3 nm thick amorphous Al 2 O 3 shell (orange). (b) Schematic of an aggregate, where Al-NPs are close packed in the ABCBA stacking structure. The five layers are labeled by Roman numerals. 4.3 Analysis of reactive molecular dynamics simulations 4.3.1 Oxidation dynamics of three Al-NP aggregates To visualize the reaction processes of the 26, 36 and 46 nm Al-NP aggregates, Figure 4.3 shows their snapshots at 0, 1 and 2 ns, where the temperature is color-coded. At time t = 0, the spherical Al-NPs are closely packed (Figure 4.3 (a), (d) and (g)). Here, the amorphous shells of the cooler Al-NPs (light blue) in layers II-V protect the Al cores from reacting with environmental oxygen. Immediately after the reaction begins, Al-NPs start to deform and fuse together in a layer-by-layer fashion. The reaction proceeds with different speeds for different systems. At t = 1 ns, for instance, layer IV is still intact for D = 46 nm (Figure 4.3(i)), while it has partially reacted for D = 26 nm (Figure 4.3(b)) and 36 nm (Figure 4.3(e)). As the reaction progresses, at t = 2 ns, all the Al-NPs for D = 26 85 nm fuse into an ellipsoidal agglomerate with temperature as high as 6,000 K (seen red in Figure 4.3(c)). Similar agglomerates are observed at t = 2 ns for D = 36 nm (Figure 4.3(f)) and 46 nm (Figure 4.3(j)), though layer V still remains unreacted in the latter. When the left layer of Al-NPs starts to fuse with the right layers, the boundaries between core and shell structures of the right Al-NPs become blurred. Also, core Al atoms of these Al-NPs start to react with oxygen. We see vast shell and core atoms are ejected into the environment backwards the direction of combustion. Figure 4.3 Snapshots of Al-NP aggregates for D = 26 nm (a-c), 36 nm (d-f), and 46 nm (g-j) at times 0, 1 and 2 ns. The temperature is color-coded. 4.3.2 Agglomeration of different layers of Al-NPs aggregates As shown in Figure 4.3, a common phenomenon for all three aggregates is the agglomeration of Al-NPs. In order to qualify the process of agglomeration, we first calculate the center-of-mass position of each layer Xα(t) (α = I, II, III, IV, V) in the x direction as a function of time t. From Figure 4.4, as the agglomeration progresses for 86 larger t, we observe the decrease of the distances between the layers. As a measure of this decrease, we calculate the distance between layers I and V: ΔX I-V (t) = X V (t) − X I (t). Smaller ΔX I-V (t) values signify more extensive agglomeration. Figure 4.4 Time evolution of the center-of-mass positions of Al-NPs in layers I-V in the x direction. (a)-(c) are for Al-NP diameter D = 26, 36 and 46 nm, respectively. 87 Figure 4.5 plots the normalized distance ΔX I-V (t)/ ΔX I-V (0) as a function of time. We see more rapid decreases of the distance for smaller diameters. At t = 2 ns, ΔX I- V (t)/ ΔX I-V (0) = 0.39, 0.49 and 0.71 for D = 26, 36 and 46 nm, respectively. This result indicates faster agglomeration of Al-NP aggregates with smaller diameters. This is understandable considering the larger surface-to-volume ratio of smaller Al-NPs. As demonstrated in our previous study, the penetration of hot atoms from reacted Al-NPs into adjacent unreacted Al-NPs is the main cause of the combustion propagation [167]. The larger surface-to-volume ratios of smaller Al-NPs accelerate the oxidation and heat propagation between the particles. This in turn leads to an increase of the extent of agglomeration. Figure 4.5 Time evolution of Al-NP agglomeration as measured by the distance ΔX I-V (t) between the center-of-mass positions of layers I and V in the x direction. The normalized distance ΔX I-V (t)/ ΔX I-V (0) is plotted as a function of time t for D = 26, 36 and 46 nm. 88 4.3.3 Consumption of remaining core Al in Al-NP aggregates To quantify the different burning rates of the three systems, we measure the rate of core Al consumption. This is done by calculating the fraction of unreacted core Al atoms as a function of time for D = 26, 36 and 46 nm as shown in Figure 4.6. The slope of each curve in Figure 4.6 represents the consumption rate of core Al for the corresponding system. We see that the consumption rates are different not only for different systems, but also for different layers of the close pack stacking. Each curve in Figure 4.6 exhibits five distinct stages of reaction, which are shown as regions with different colors. We have found that different stages correspond to the reaction of different layers of Al-NPs. This is evidenced in the insets of Figure 4.6, which shows snapshots of the core Al for the 46 nm Al-NP in layer IV at two different times within reaction stage IV. At the beginning of stage IV (t = 1.5 ns), the Al-NP in layer IV is unreacted, whereas later in stage IV (t = 1.8 ns) Al core atoms are jetting out due to intense reaction. The different reaction stages are accordingly labeled with layer numbers I to V in Figure 4.6. An important observation is that the core Al consumption rate is higher for layers I and III (red and green regions) compared with layers II and IV (yellow and grey regions). Note that layers I and III have three Al-NPs each, while each of layers II and IV only has one Al-NP (Figure 4.2(b)). Namely, the more is the number of Al-NPs within a layer, the higher is the reaction rate. The aggregation of Al-NPs thus accelerates the consumption of core Al. The combustion of larger diameter Al-NP aggregates takes 89 longer and remains less complete due to the smaller surface to volume ratio. A similar size dependence of reaction completeness was obtained in an experiment[141] on the combustion of fuel droplets. I II III IV V Figure 4.6 Fraction of unreacted core Al atoms as a function of time for 26 (blue curve), 36 (green curve) and 46 nm (red curve) Al-NP aggregates. Five stages with different core Al consumption rates are shown as regions of different colors and are denoted with Roman numerals. The insets show snapshots of the core Al for the 46 nm Al-NP in layer IV at times t = 1.5 and 1.8 ns. 4.3.4 Reaction crossover of small intermediate oxidation fragments To investigate oxidation mechanisms, we next examine intermediate oxide products using fragment analysis. An oxide fragment is identified as a set of atoms that are connected by Al-O bonds. A large number of oxide fragments are ejected out from the aggregates. In all three systems, the majority of ejected fragments are found to be 90 small clusters having less than six atoms each. Figure 4.7 shows an example of oxide fragments for 26 nm Al-NP aggregates at t = 2 ns. In the figure, Al 2 O, AlO and AlO 2 fragments are highlighted. A notable observation is the formation of an Al 2 O 3 fragment from AlO and AlO 2 . Al 2 O 3 is the final product of the Al oxidation, and its production from two primary intermediates, AlO and AlO 2 , was previously proposed [157]. Figure 4.7 A snapshot of oxide fragments in the 26 nm Al-NP aggregate at t = 2 ns, showing the formation of Al 2 O 3 from AlO and AlO 2 . The largest oxide fragment is the shell shown in the right. Figure 4.8, Figure 4.9 and Figure 4.10 show the number of Al 2 O, AlO and AlO 2 fragments as a function of time in D = 26, 36 and 46 nm Al-NP aggregates, respectively. In all three systems, the Al 2 O production exhibits peaks at times coincident with the high 91 core Al-consumption rates in stages I, III and V (see Figure 4.6). This suggests that Al 2 O is the direct product from core Al consumption. Such features are not visible in the AlO and AlO 2 curves, suggesting that they are secondary products. Figure 4.8 Number of Al 2 O (a), AlO (b) and AlO 2 (c) fragments as a function of time in the 26 nm Al-NP aggregate. 92 Figure 4.9 Number of Al 2 O (a), AlO (b) and AlO 2 (c) fragments as a function of time in the 36 nm Al-NP aggregate. 93 Figure 4.10 Number of Al 2 O (a), AlO (b) and AlO 2 (c) fragments as a function of time in the 46 nm Al-NP aggregate. To study the time evolution of the primary intermediate products of the reaction, Figure 4.11 plots the fraction of small oxide fragments[62] fα (α = Al 2 O, AlO, AlO 2 ) for 26, 36 and 46 nm Al-NP compared to the total number of fragments. For all three 94 systems, initially produced fragments are predominantly Al 2 O (Figure 4.11(a)). The fraction of Al 2 O drops sharply, accompanied by rapid increase of the fraction of AlO 2 (Figure 4.11(b)). The crossover of the major fragment species from Al-rich Al 2 O to O- rich AlO 2 occurs at t = 0.3, 0.4 and 0.5 ns for 26, 36 and 46 nm Al-NP aggregates, respectively. The crossover occurs the fastest in the 26 nm system, indicating a faster oxidation for smaller Al-NPs. Subsequently, the fraction of AlO 2 starts to decrease, while that of AlO increases gradually. This leads to a second crossover of the major intermediate product from AlO 2 to AlO at t = 2.0, 1.6 and 1.5 ns for 26, 36 and 46 nm Al- NP aggregates, respectively. The crossover occurs slower for smaller diameter Al-NP aggregates. The two crossovers of the major intermediate oxide products may be understood as follows. Initially, a large amount of Al atoms react with environmental oxygen to form Al-rich clusters according to reactions such as 4Al + O 2 → 2Al 2 O. At around 0.4 ns, a large amount of oxygen flows from the right (O-rich) cold region into the left (O- depleted) hot region, which leads to locally oxygen-rich environment. This results in rapid production of AlO 2 fragments with reactions such as Al + O 2 → AlO 2 . This intense oxidation depletes the environmental oxygen around Al-NPs, and consequently less oxygen becomes available for further oxidation of Al. This oxygen-poor environment preferentially drives reactions such as 2Al + O 2 → 2AlO and 2AlO 2 → 2AlO + O 2 . The formation of Al 2 O 3 fragments through the reaction, AlO + AlO 2 → Al 2 O 3 , was observed and seems a preferred pathway over others, but its scarcity implies a rather long reaction time. 95 Figure 4.11 Time evolution of the fraction of major small oxide fragments, (a) Al 2 O, (b) AlO, and (c) AlO 2 , in 26, 36 and 46 nm Al-NP aggregates, shown in red, green and blue lines, respectively. 96 4.3.5 Evolution of largest oxide fragments in the aggregates Between the Al core consumed by the reaction (Figure 4.6) and the intermediate oxide fragments ejected into the environment (Figure 4.11) lies the oxide shell. We identify the shell as the largest bonded aluminum-oxide fragment. The oxidation dynamics involves a large number of atoms moving into (e.g. core Al consumption) and out of (e.g. intermediate oxide fragments ejection) the shell. Figure 4.12 shows the time evolution of the composition of the largest oxide fragment, i.e., the ratio between the number of O atoms and that of Al atoms in the shell N O shell /N Al shell [169]. In the 26 nm system, the ratio decreases until it reaches the minimum at ~ 1.3 ns, and starts to increase. The minimum coincides with the end of the core Al consumption. This is understandable because more Al atoms are incorporated into the shell due to the core Al consumption. The similar minimum is observed for 36 nm system, but at much later time ~ 1.8 ns. In contrast, for 46 nm system, the ratio is a decreasing function of time, because the core Al atom is not consumed up at the end of the simulation time. 97 Figure 4.12 The ratio between the number of O atoms and that of Al atoms in the shell as function of time for D = 26 nm (red), 36 nm (green) and 46 nm (blue). The initial ratio for the amorphous Al 2 O 3 shell is 1.5. 98 Chapter 5. Oxidation Dynamics of Aluminum Nanorods In this chapter, we extend our study on Al nanopowders from the spherical shape of particles to a different morphology, i.e., rods, since manufacturing technology is ripe to synthesize rod/wire-like Al nanostructures, along with advantages of distinctive reaction behaviors of non-spherical structures. Section 5.1 provides the background and motivation of studying oxidation of Al nanorods (Al-NRs). The simulation methods and results are described in Sections 5.2 and 5.3, respectively, followed by the conclusion in Section 5.4. 5.1 Background and motivation In order to overcome various issues [16, 33, 34, 38-40] during the combustion of conventional Al nanopowders mainly composed of spherical Al nanoparticles, new morphologies of Al powders are highly demanded. State-of-the-art technologies are now able to synthesize Al nanostructures with varying aspect ratios (nanorods [45, 46] and nanodisks [53]) to enhance the oxidation efficiency. Many different nanostructures have been used to in the field of nanocomposites, such as mesoporous structures [43, 44] and nanorods based self-assembly access [45, 46] to increase the reactivities between the fuel and oxidizer. Among all non-spherical shapes of Al nanostructures, Al nanorod (Al-NR) is the most accessible and efficient one that will be an ideal candidate for future nanofuel. With more advanced technologies [50-53] to produce Al-NRs, it is essential to understand fundamental oxidation mechanisms that differentiate the reaction dynamics of Al-NRs 99 from that of Al-NPs. An important question is what is the size effect on the burning speed? Since a native oxide coating exists in all Al-NRs synthesized by the above- mentioned methods, also essential to answer is what is the role of the oxide shell in the oxidation dynamics? 5.2 Simulation setup for three Al nanorod systems We have performed reactive MD simulations of three oxide-coated Al-NRs with diameters D = 26, 36 and 46 nm. Each cylindrical Al-NR has a length of 3D and consists of a cylindrical Al core of diameter Dʹ′ (= D−6 nm) covered by a 3 nm amorphous alumina (Al 2 O 3 ) shell (see Figure 5.1). The cylindrical core is cut from face-centered cubic Al crystal. The amorphous shell is prepared by quenching molten alumina to a temperature of 300 K, then removing the concentric innermost and outermost cylindrical parts to have a 3 nm-thick oxide layer. The combined Al core and alumina shell is then embedded in oxygen environment at 300 K. The rectangular simulation box has side lengths of 5D, 3D and 3D in the x, y and z directions, where the NR axis is in the x direction. The system is first thermalized for 30 ps with the environmental oxygen atoms intact, in order to form an interface between the Al core and alumina shell. The total number of atoms is 4,518,228, 11,585,041 and 23,566,179 in the 26, 36 and 46 nm Al- NR systems, respectively. Breakdown of the number of atoms into the core, shell and environment subsytems is provided in Table 5.1. The interatomic potentials used for Al 2 O 3 [170], metallic Al [171], and Al oxidation are identical to those used in our previous RMD simulations [166, 172, 173]. 100 The equations of motion are integrated using the velocity Verlet algorithm with a time step of 1 fs. D" D" L=3D" Oxygen' D" D" D’=D'6'nm' Al' Al 2 O 3' D" L’=3D'6'nm' x" z" y" Figure 5.1. Schematic of an Al-NR with diameter D and length L = 3D. The Al core (shown in green) is covered with a 3 nm-thick alumina shell (orange), and is embedded in oxygen environment (blue). Table 5.1 Number of atoms of different components in D = 26, 36 and 46 nm aluminum nanorod (Al-NR) systems. D = 26 nm D = 36 nm D = 46 nm Number of core Al atoms 1,360,903 4,343,761 9,989,618 Number of environmental O atoms 1,360,360 3,631,504 7,526,316 Number of shell Al atoms 719,199 1,444,090 2,419,850 Number of shell O atoms 1,077,765 2,165,687 3,630,395 The oxidation dynamics is studied by heating the left one-third of the Al-NR to a temperature of 1,100K within 9 ps, while the atomic positions of the rest of the Al-NR 101 and the environmental oxygen are fixed at those in the 300 K simulation. The radial distribution function shows that the Al core in the heated left end is liquid at 1,100 K. Subsequently, the constraints are removed and the system undergoes oxidation dynamics in the microcanonical ensemble for 1 ns. Figure 5.2 shows middle-xz slice snapshots of the D = 26 ((a)-(c)), 36 ((d)-(f)) and 46 nm ((g)-(i)) systems at times t = 1 ps, 500 ps, and 1 ns. Initial heating of the left one- third of the Al-NR causes slight volume expansion of the left end as shown in Figure 5.2(a), (d) and (g) for D = 26, 36 and 46 nm, respectively. Subsequently, the oxidation front propagates from left to right along the nanorod. At the same time, each Al-NR shortens and its left end deforms into an ellipsoid as shown in Figure 5.2(b), (e) and (h) for D = 26, 36 and 46 nm, respectively. By 1 ns, the entire 26 nm rod has become an ellipsoidal aggregate (Figure 5.2(c)), and the majority of the 36 nm rod except its right end has become an ellipsoidal aggregate (Figure 5.2(f)). In contrast, only half of the 46 nm rod aggregates as shown in Figure 5.2(i). This indicates complete burning of D = 26 nm Al-NR by the end of the simulation, while D = 36 and 46 nm Al-NRs remain incompletely burned. 102 Figure 5.2 Middle-xz slice snapshots of D = 26, 36 and 46 nm Al-NR systems at times 1 ps, 500 ps, and 1 ns. Only shown are Al atoms from the shell (red dots), O atoms from the shell (yellow dots), and environmental O atoms (black dots). 5.3 Analysis of the reactive molecular dynamics 5.3.1 Temperature evolution of the three Al-NRs To quantify the size dependence of the reactivity, Figure 5.3(a) plots the temperature as a function of time for D = 26, 36 and 46 nm systems. For all Al-NRs, the temperature increases monotonically as soon as the oxidation reaction starts. Also note that the temperature is always higher for thinner Al-NRs. To study heat release in the three systems, Figure 5.3(b) plots the rate of temperature increase dT/dt as a function of time. While dT/dt > 0 at all times in all three systems, the dT/dt exhibits a valley, corresponding to the lowest heat release, at t = 100 ps for all Al-NRs. The minimum 103 dT/dt values are different for different diameters, with D = 26 nm Al-NR having the highest value. After passing the valley, the rates start increasing, signifying the accelerated heat release. The dT/dt value is always larger, corresponding to faster heat release, for smaller Al-NRs. Figure 5.3 (a) Temperatures of the three Al-NRs as a function of time. (b) The rates of temperature increase dT/dt of the three Al-NRs as a function of time. 104 5.3.2 Evolution of remaining core Al atoms for the three Al-NRs The decrease of dT/dt for D = 26 nm Al-NR after 800 ps is due to the completion of burning. To confirm this interpretation, Figure 5.4 shows the time evolution of the fraction of core Al atoms remained non-bonded to O atoms to quantify the consumption rate of the metallic core Al for the three Al-NRs. Smaller Al-NRs exhibit faster consumption rates. The Al core of D = 26 nm Al-NR is nearly completely consumed in 800 ps, while that of D = 46 nm Al-NR is only half way consumed by the end of the simulation. Figure 5.4 The fraction of remaining core Al atoms as a function of time for the diameter D = 26, 36 and 46 nm. 105 5.3.3 Reaction dynamics on the oxide shell of Al-NRs In order to understand the cause of the crossover from deceleration to acceleration of heat release in Figure 5.3(b), Figure 5.5 shows part of the middle xz cross section of the oxide shell in D = 36 nm Al-NR before (a) and after (b) the crossover. Here, the oxide shell is identified as the largest fragment interconnected by Al-O bonds, which may contains some of the Al (red) and O (yellow) atoms from the original shell from the beginning, Al atoms from the original Al core (blue), and environmental O atoms (black). In the early stage, some O atoms from the shell start to react with core Al atoms (see Figure 5.5(a) at 50 ps). Namely, we observe de-oxidation of the amorphous alumina shell involving reactions between liquid Al(l) and solid Al 2 O 3 (s) at the core-shell interface. A similar high-temperature de-oxidation of alumina at Al 2 O 3 (s)-Al(l) interfaces was observed by Oh et al. using high-resolution transmission electron microscopy during vapor-liquid-solid growth of alumina nanowires [174]. After the initial dissolution of O atoms from the shell, steady reaction between Al and O atoms generates heat, which assists continuing reaction. Subsequently, environmental O atoms start to get absorbed onto the outer surface of the shell, while more reactions occur at the core-shell interface (see Figure 5.5(b) at 200 ps). Namely, the de-oxidized amorphous shell re-oxidizes, involving reactions between gas-phase O(g) and Al 2 O 3 (s). These results indicate an unexpected role of the oxide shell as a reactive medium. Rather than being an inert separator between the inner Al fuel and outside oxidizer, the oxide shell continuously exchanges its constituent atoms with core Al and environmental O atoms, thus promoting interfacial reactions and heat release. 106 !"#$ %$&'$ ()$&'! *+,-$./$"0+'1! 23$41$ !"#$ %$&'$ (&)*+,&'-&./0$1$/.,'2! 344$52$ Figure 5.5 Snapshots of part of the middle xz slice of the 36 nm system, showing the largest oxide fragment at (a) 50 ps and (b) 200 ps. Original shell Al, shell O, core Al and environmental O atoms are shown in red, yellow, blue and black colors, respectively. To study the dynamics of the oxide-shell growth, Figure 5.6 plots the time evolution of the composition of its constituent atoms from different sources, i.e., the numbers of Al and O atoms from the original shell, original core Al atoms, and original environmental O atoms. Here, the oxide shell is identified as the largest oxide fragment interconnected by Al-O bonds. The composition changes continuously reflecting reactions. The compositions of the oxide shell are shown by stacked bar graphs in Figure 5.6 (a), (c) and (e) for D = 26, 36 and 46 nm systems, respectively. The numbers atoms from different sources are shown individually in Figure 5.6 (b), (d) and (f) for D = 26, 36 and 46 nm systems, respectively. For all three Al-NRs, core Al atoms are continuously incorporated into the oxide shell, while some of the O atoms are dissolved into the Al core, starting immediately at time 0. The rapid initiation of the core-shell reaction is likely due to the high pressure and temperature of the liquefied Al core. This is similar to the sub-oxide formation mechanism at alumina-liquid Al interfaces suggested by Laurent et al.[175] On the other hand, environmental O atoms only appear after 200 ps in the 107 oxide shell. This onset time roughly corresponds to the accelerated heat release in Figure 5.5(b). We thus attribute the accelerated heat release after 200 ps to the onset of O(g)- Al 2 O 3 (s) reactions. Figure 5.6 Time evolution of the numbers of atoms from different sources in the oxide shell. The composition of the shell is shown by stacked bar graphs in (a), (c) and (e) for D = 26, 36 and 46 nm Al-NRs, respectively. The numbers atoms from different sources are shown individually in (b), (d) and (f) for D = 26, 36 and 46 nm systems, respectively. 108 Environmental O, core Al, shell O and shell Al atoms are shown in black, blue, yellow and red colors, respectively. 5.3.4 Oxide shell as a reactive medium To study the stoichiometry of the oxidation product, Figure 5.7 plots the ratio of the number of O atoms to that of Al atoms in the oxide shell as a function of time for the three Al-NRs. The oxide shell, identified as the largest fragment interconnected by Al-O bonds, changes from the original amorphous Al 2 O 3 with stoichiometry 1.5 to more Al- rich for all three Al-NRs. Once oxidation dynamics starts, Al atoms from metallic core at the left heated side react with shell O atoms rapidly. Only for D = 26 nm Al-NR, the stoichiometry of the largest oxide fragment ascends around 800 ps. This is primarily because the metallic Al atoms in the system are completely consumed. But there are still environmental O atoms absorbed onto the surface of the largest oxide fragment after 800 ps, so the stoichiometry continues to increase until then. For D = 36 and 46 nm Al-NRs, the ratios keep descending until 1 ns. The stoichiometry of D = 46 nm Al-NR is the lowest by 1 ns, since the reaction area between the metallic core Al and the oxide shell is much more significant than smaller Al-NRs, which leads to a greater change of the stoichiometry. The time evolution of the stoichiometry in D = 26 nm Al-NR directly reflects the different stages of the oxidation mentioned above, i.e., de-oxidization at the core-shell interface, followed by re-oxidization of the Al-rich shell at the environment- shell interface. In this system, the complete consumption of core Al at 800 ps (see Figure 5.4) leads to the diminishing partial pressure of Al, and the rest of reaction is between the hot Al-rich aggregate (Figure 5.1(c)) and environmental oxygen. For the larger (D = 36 109 and 46 nm) Al-NRs, the second reaction stage is not reached before 1 ns (Figure 5.7) due to the slower reaction speeds (Figure 5.4). This confirms the active role of the oxide shell as a reactive medium, instead of an inert dead weight. Figure 5.7 The ratio of the number of O atoms to that of Al atoms in the oxide shell as a function of time for D = 26, 36 and 46 nm Al-NRs. 5.3.5 Spatial propagation of oxidation reactions Next, we study the spatial propagation of the oxidation reaction along the Al-NR. To do so, Figure 5.8(a) shows snapshots of the upper half middle xz slice of the oxide shell (i.e., the largest oxide fragment) in D = 36 nm Al-NR at different times. The positions of the oxidation front (see red arrows that identify the rightmost position of the left side of the oxide shell at the axial center) at different times shows that the oxidation reaction propagates from left to right along the Al-NR axis. The temperature distribution along the radial (R) and axial (X) directions of the same Al-NR is shown in Figure 5.8 (b). The hot spots shown in red color are located on the outer left surface of the shell. 110 Oxidation reaction on the shell generates heat, which promotes further reactions and pushes the oxidation front rightward along the Al-NR axis. Figure 5.8 (a) Snapshots of the upper half middle xz slice of the largest fragment in D = 36 nm Al-NR at different times. Al and O atoms originated from shell are shown in red and yellow. Al atoms originated from core are shown in blue. The red arrows show the 111 position of the oxidation front. (b) Radial (R) and axial (X) temperature distribution of the same Al-NR. 5.3.6 Size effect on the burning rate of Al-NRs In order to quantify the size effect on the combustion speed, we calculate the oxidation front position as the average position of the rightmost 100 core Al atoms in the largest oxide fragment, as well as the melting front position in core Al as the rightmost position where the temperature of the un-reacted core reaches the melting point, 933.52 K. Figure 5 shows the oxidation and melting fronts as a function of time for D = (a) 26, (b) 36 and (c) 46 nm Al-NRs, respectively. Linear regression yields the average oxidation- front speed in D = 26, 36 and 46 nm Al-NRs to be v o = 88, 70 and 60 m/s, respectively. In comparison, the average melting front in core Al propagates at v m = 91 m/s, 72 m/s and 66 m/s for D = 26, 36 and 46 nm Al-NRs, respectively. The average v m is larger than v o in all systems, which signifies that the former dictates the burning of Al-NRs. The faster v o ’s of the smaller Al-NRs are likely caused by the larger surface-to-volume ratios. This is understandable in the light of the active roles of the core-shell and environment-shell interfaces as reaction spots and heat sources to support continuous melting of the metallic core Al. Here, we should note that both the melting and oxidation fronts in D = 26 nm Al-NR reach the right end before 1 ns, with the former reaching earlier than the latter (Figure 5.9 (a)). In D = 46 nm Al-NR, v o is slower than v m in the first 600 ps, exceeds the latter around 600 ps, and then slows down below v m toward 1 ns. This can be understood as a result of oxidation in the radial (in addition to axial) direction that contributes 112 significantly in larger Al-NRs. In D = 46 nm Al-NR, the axial and radial oxidation fronts merge abruptly around 600 ps to cause a sudden advance of the oxidation front. Figure 5.9 Oxidation-front (red) and melting-front (blue) positions as a function of time for D = (a) 26, (b) 36 and (c) 46 nm Al-NRs. 113 5.3.7 Heat propagation in three Al-NRs For a closer look at the heat propagation in Al-NRs, Figure 5.10 shows the radial temperature distribution at yz cross-sections in the left one-third, middle one-third and right one third parts of D = 36 nm Al-NR at different times. The distributions at (a) 300 ps (green) in the left section, (b) 600 ps (yellow) in the middle section, and (c) 1 ns (red) in the right section clearly show that the temperature near the shell is higher than that in the core. Figure 5.10(d), (e) and (f) respectively show the temperature distribution of the same cross sections in the left, middle and right parts of the 36 nm Al-NR at different times. The temperature distribution from left to right at different time shows how the heat generates and propagates along the nanorod on the shell. We see that the hottest area always is the shell region. This establishes the role of the shell as a heat source to sustain the oxidation reaction. 114 Figure 5.10 Radial temperature distributions at yz cross sections in (a) left, (b) center and (c) right parts of D = 36 nm Al-NR at times 1 ps, 300 ps, 600 ps and 1 ns. The black lines delineate the boundary between core and shell in different sections as the oxidation progress. Also shown are the temperature distributions on the same yz cross sections in (d) left, (e) center and (f) right parts show for D = 36 nm Al-NR, temperature distribution along the radial on Y-Z plane for different sections at different times. 115 Chapter 6. Shock-Induced Detonation Reaction in High Explosives RDX In this chapter, we switch our research topic to high explosives, explicitly on RDX solid. Section 6.1 provides the background and motivation of studying detonation of RDX solid. The simulation methods and results are described in Sections 6.2 and 6.3, respectively, followed by the conclusion in Section 6.4. 6.1 Background and motivation To study high explosives, we focus on MD simulation of RDX crystals because of their practical usage and the availability of reliable force field. The key scientific questions for detonation of RDX solid are: What is the reaction time and what are the intermediate reaction products that determine it? It has been realized that the chemical dynamics behind the shock front in energetic materials located in a ~100 nm thin layer and on the ~100 ps time scale is critical to understand microscopic details of detonation [84, 93]. Unfortunately, no MD simulation capable of describing chemical reactions has been performed to encompass the large length (~100 nm) and time (~100 ps) scales [176]. In order to overcome this computational challenge, we have developed a scalable parallel implementation of reactive force-field (ReaxFF) MD simulation based on spatial decomposition and message passing [177]. To describe chemical reactions with moderate computational costs, the first principles-based ReaxFF allows bond breaking and bond formation through reactive bond orders and dynamical charges by employing an 116 electronegativity-equalization scheme [106, 178] In this chapter, we present ReaxFF MD simulations validated by quantum molecular dynamics (QMD) simulations to study shock-induced detonation of RDX. 6.2 Simulation setup for RDX crystal The initial system setup is composed of 168 × 5 × 5 unit cells of RDX crystal in an MD simulation box of dimensions 222.8 × 5.787 × 5.354 nm 3 at room temperature. Real HEs are defective (e.g. plasticizer bonded, wax added, and porous) [69], and these heterogeneities lead to complex multidimensional flows even in macroscopically unidirectional detonation [179]. In particular, voids cause localized hot spots, which essentially drive the decomposition reaction [180]. To mimic this effect, a void of size 3 × 3 × 3 unit cells is randomly inserted in every 5 × 5 × 6 unit cells of the RDX crystal. Figure 6.1(a) and (b) show the atomic structure of an RDX molecule and the configuration of the RDX crystalline unit cell, respectively. The number of RDX molecules inside the system is 27,552, amounting to the total number of atoms to be 578,592. Starting from this initial configuration, we perform MD simulations with a time step of 0.1 fs up to 70 ps. Figure 6.1(c) shows a schematic diagram of shock loading. To model the shock wave and subsequent detonation and expansion, we employ a planar impact loading by a rigid-wall piston with the speed of v p = 6 km/s from the right end onto the system. When the shock front hits the rigid wall at the left end of the simulation box, it is reflected back to the right. It takes about 11 ps for the shock wave to travel for a 117 distance of 110 nm, indicating the shock speed of ~ 10 km/s. This is above the RDX detonation speed of 8.75 km/s, suggesting an overdriven detonation [181]. After the shock front bounces back towards right, the piston is removed to allow the detonated RDX to expand freely. The expansion continues until the volume of the system becomes three times the original value. Figure 6.1 (a) RDX molecule, where gray, cyan, blue and red spheres represent C, H, N and O atoms, respectively. (b) RDX unit cell. The 8-molecule (or 168-atom) unit cell has the lattice parameters of a = 13.182 Å, b = 11.574 Å, c = 10.709Å, and α = β = γ = 90˚. (c) Schematic diagram of shock-induced detonation and propagation. 118 6.3 Analysis of detonation results 6.3.1 Temperature and pressure evolution To demonstrate the time evolution of temperature distribution, Figure 6.2(a) shows a colored map of the temperature as a function of the x position and time. The initial length of RDX is X 0 = 222.8 nm. The highest temperature of the system reaches 3,000 K near the shock front (shown as diamond-shaped red spots in Figure 6.1(d)) before 20 ps, and during compression at 20-30 ps. To quantify the time evolution of pressure distribution, we calculate the virial pressure from the atomic positions and velocities, and average it over a slab of 1 nm thickness along the x direction around each atom. Figure 6.2(b) shows a colored map of pressure as a function of the x position and time. While the peak value of the pressure at the detonation front is as high as 80 GPa, the final overall pressure of the system is about 26 GPa, which is close to a reported experimental value of 27.2 GPa [75]. 119 Figure 6.2 Colored map of (a) temperature and (b) pressure as a function of the x position and time. 6.3.2 Fragment analysis: two-stage reaction times To study the chemical reaction pathways, we have performed fragment analysis, where a cluster of covalently bonded atoms is counted as a molecule. Here, a pair of atoms are considered connected if the bond order of the pair is greater than 0.3 [182] and their distance is less than a critical value slightly larger than the corresponding covalent bond length. Figure 6.3 shows the number of molecular products as a function of time. These include the final products (H 2 O, N 2 , and CO) of the overall reaction of RDX decomposition: 120 C 3 H 6 N 6 O 6 → 3 H 2 O + 3 N 2 + 3 CO . (6.1) We observe rapid production of N 2 and H 2 O during the initial expansion phase followed by plateaus in the later stage. We also observe a very slow production of CO. To estimate the time constant τα for the production of the α-th final product (α = H 2 O, N 2 , CO), and associated error bar, we first calculate the normalized number of molecular fragments by dividing the number of fragments in Figure 6.3 by the number of initial RDX molecules (i.e. 27,552). We fit the yield of the corresponding molecular fragments as a function of time as η α (t)=1−exp(−t /τ α ) . (6.2) Here, ηα(t) is defined as the ratio of the observed number of the α-th fragments at time t to the asymptotic number for t → ∞ expected from Eq. (6.1). Namely, the expected number of molecules produced from one RDX molecule is 3 for H 2 O, N 2 , or CO. The calculated yield ηα(t) as a function of time (α = H 2 O, N 2 , or CO) is then fitted to Eq. (6.2) (the origin of time in Eq. (6.2) is the start of the expansion phase). Least square fit produces both the expected time constant and its error bar: τ N 2 = 10 ± 2 ps, τ H 2 O = 30 ± 4 ps, and τ CO = 900 ± 400 ps. Namely, the detonation reaction proceeds in two stages: Rapid production of N 2 and H 2 O, followed by much slower production of CO. The rapid production of N 2 and H 2 O at an early stage of detonation is consistent with earlier simulation by Strachan et al. [183]. In their simulation for 5 ps, the dominant products were N 2 and H 2 O. Also, slightly different multistage reactions were observed experimentally by Anisichkin [184] for a related HE, RDX/TNT mixture, using an 121 isotope tracer method. Figure 6.3 Number of molecular fragments as a function of time. As the shock front begins to reflect, a rapid production of H 2 O, OH, and N 2 is observed. Shortly after the expansion phase begins, various chemical products such as CO, CO 2 , and NO are produced. 6.3.3 Confirmation of the two-stage reaction in an alternative simulation In order to verify that the two-stage reaction is not an artifact of the simulation schedule, we have performed another simulation using a different schedule. The alternative simulation only has a compression phase (instead of the compression and holding phases in the original simulation) before the expansion phase. Figure 6.4 shows two distinct stages for the production of the final products using the new schedule. Namely, a large number of N 2 and H 2 O molecules are produced rapidly well below 100 ps, whereas not many CO molecules are produced in the same time frame, signifying 122 their much delayed production. This conclusion is identical to that from the original simulation, suggesting that the two-stage reaction we have found is an intrinsic property of RDX. Figure 6.4 Number of molecular fragments as a function of time in an alternative simulation, where the expansion phase directly follows the compression phase. 6.3.4 Reaction pathways analysis In order to identify the reaction pathways for the rapid production of N 2 and H 2 O, we backtrack where the atoms composing individual N 2 and H 2 O molecules originate in the initial configuration of the simulation. For some of the N 2 and H 2 O products, all the atoms that constitute a molecule originate from a single RDX molecule, i.e., unimolecular pathways. For others, atoms from different RDX molecules form N 2 and H 2 O molecules, i.e., intermolecular pathways. Examples are shown in Figure 6.5, where the H, O and N atoms that constitute the circled N 2 and H 2 O products in Figure 6.5(b) are highlighted, 123 respectively, with green, yellow and magenta colors in the initial configuration in Figure 6.5(a). In Figure 6.5(a), one green H atom belongs to a –CH 2 group of one RDX molecule, while another green H atom belongs to a –CH 2 group of another RDX molecule and one yellow O atom belongs to an –NO 2 group of the latter RDX molecule. Immediately after the reversed shock front traveling rightward passes, cleavages of C and H atoms from –CH 2 groups and N and O atoms from –NO 2 group release the highlighted H and O atoms, which together form the H 2 O molecule circled in Figure 6.5(b), signifying an intermolecular pathway. In contrast, for the N 2 molecule circled in Figure 6.5(b), both magenta N atoms composing the N-N bond originate from one RDX molecule in the initial configuration in Figure 6.5(a), i.e., a unimolecular pathway. We have examined the origin of the other N 2 molecules as well, as shown in Figure 6.6. Interestingly, most of them are from N-N bonds within single RDX molecules, even though –NO 2 functional groups cleave first when the shock front passes RDX, as was seen in a previous simulation [70]. Experimental results also indicate that NO 2 is a direct product of N-N hemolysis in the initial reaction stage under shock [185]. 124 O H H N N H H O N N (b) (a) (c) Figure 6.5 (a) Initial configuration of RDX molecules in a unit cell, where gray, cyan, blue and red spheres represent C, H, N and O atoms, respectively. Atoms that later will be bonded to form N 2 and H 2 O molecules are highlighted with different colors, i.e., H, O, and N are in green, yellow and magenta, respectively. (b) Molecular fragments formed behind the reversed shock front traveling rightward, where N 2 and H 2 O molecules formed by the highlighted atoms in (a) are enclosed in circles. (c) Fraction of N 2 and H 2 O molecules formed by atoms from single RDX molecules as a function of time. 125 To better identify the source of the most abundant fragments (N 2 and H 2 O), Figure 6.5(c) plots the fraction of N 2 and H 2 O molecules that are produced by unimolecular pathways, i.e., all constituent atoms of each molecule originate from a single RDX molecule in the initial configuration: f α intra (α = N 2 , H 2 O). We see that N 2 production mechanism is predominantly unimolecular. Namely, 75% of the N 2 products are from N-N bonds within single RDXs. In contrast, only 15% of H 2 O products are formed by H and O atoms from single RDXs, i.e., H 2 O production mechanism is intermolecular. Figure 6.7 provides detailed analysis of the origin of the H and O atom of H 2 O. This implies that N atoms react locally, while H and O atoms move more actively to react with those from further non-adjacent RDX molecules. This is consistent with simulation results by Wu et al. [57], who found a catalytic behavior of water in the detonation of HE. They found that H 2 O actively participates in reactions by transporting oxygen between different fragments, instead of being a mere stable final product. The role of hydrogen as long-ranged reaction participants is understandable because of their lightest mass. Accordingly, H atoms can travel farther to bond with O atoms from other reactant molecules to form H 2 O. 6.3.5 Origin of the constituent atoms of N 2 and H 2 O products To study the reaction pathways for N 2 and H 2 O production, we determine where the constituent atoms of the N 2 and H 2 O products originate from, as shown in Figure 6.6 and Figure 6.7, respectively. In Figure 6.6, the red line shows the fraction of N 2 molecules produced from N atoms within a single RDX molecule, i.e., unimolecular 126 pathways. Also shown by the blue line is the fraction of N 2 molecules that are formed by an N-N bond within a single RDX. We see that most unimolecular pathways in fact involve an existing N-N bond in a RDX molecule. This highlights the highly localized nature of reaction pathways for N 2 production. Namely, they are mostly localized in the region of N-N bonds within a single RDX molecule. These direct reaction pathways lead to the most rapid production of N 2 molecules as shown in the main text. Figure 6.6 Fraction of N 2 molecules that are formed by unimolecular pathways as a function of time. The red line shows the fraction of N 2 formed by N atoms from a single RDX molecule. The blue line shows N 2 originating from an N-N bond within a single RDX molecule. 127 Figure 6.7 Fraction of H 2 O molecules formed by different reaction pathways as a function of time. The blue line shows H 2 O formed by atoms from a single RDX molecule, i.e., unimolecular pathways. The red line shows H 2 O formed by atoms from two RDX molecules, i.e., dimolecular pathways. In Figure 6.7, the blue line shows the fraction of H 2 O products that are formed by H and O atoms that originate from a single RDX molecule, i.e., unimolecular pathways. Also shown by the red line is the fraction of H 2 O molecules that are composed of atoms originating from two RDX molecules, i.e., dimolecular pathways. We see that atoms from two adjacent RDX molecules form more H 2 O products than from a single RDX molecule. In RDX crystal, RDX molecules are spatially arranged so that a number of O and H atoms from adjacent RDX molecules are located in spatial proximity. Thus, it is not difficult for these proximate intermolecular H-O pairs to react rapidly. Furthermore, H atoms can react with O atoms located in remote RDX molecules due to their large 128 diffusivity. As a result, the fraction of dimolecular pathways is larger than that of unimolecular pathways in Figure 6.7, since more combinations of H and O atoms exist for the former to provide broader reaction channels. 6.3.6 Reaction products: large clusters To understand the reaction pathways involving the rest of the elements (particularly carbon) and the reason behind the slow production of CO molecules, we study the time evolution of the population of large clusters remained in the system. Here, we define a large cluster as that containing more number of atoms than the 21-atom RDX molecule. Figure 6.8(a) shows the fraction N C /N T of the number of atoms in larger clusters, N C = C(i)i i>21 ∑ , vs. the total number of atoms in the system, N T = C(i)i i=1 ∞ ∑ , as a function of time. Here, C(i) is the number of molecular fragments each consisting of i atoms. In the compression phase, the number of coagulated atoms, which are from the compressed part of RDX, increases linearly with time. Once no further shock loading is applied on the system, the large clusters start to decompose. Especially in the initial expansion phase, the sudden release from the right opening side leads to massive decomposition of the large clusters, which produce vast amount of small fragments such as N 2 and H 2 O. However, in the later stage, lower temperature and density cannot sustain further chemical decomposition reactions. As a result, more than 10% of the atoms remain in large clusters even at the end of the total simulated time of 70 ps. The characteristic decay time of the large clusters is estimated by the exponential fitting as τ decay = 800 ± 500 ps, which is close to the reaction time of CO: τ CO = 900 ± 400 ps. 129 Therefore, we can identify the large metastable clusters to be the major retardant of the CO-production reaction. We next perform stoichiometric analysis of the large clusters. Figure 6.8(b) shows the averaged atomic compositions, c(α) (α = C, H, O, N), of the large clusters as a function of time, which is normalized as c α =1 α ∑ . In the figure, the dashed lines show the composition for different elements in the RDX molecule: c 0 (C) = 3/21 while c 0 (H) = c 0 (O) = c 0 (N) = 6/21. Comparison of the compositions c(α) of the large clusters with those of the RDX reactant c 0 (α) shows that the large clusters at the end of the simulation are not only carbon-rich but also oxygen-rich. Namely, the average stoichiometry of the large clusters at 70 ps is summarized as C 4.94 H 4.76 O 6.99 N 4.31 (where the numbers are normalized to add up to 21—the number of atoms in an RDX molecule), which contains more C and O atoms than the RDX reactant, C 3 H 6 O 6 N 6 . Since these large C- and O-rich clusters remain metastable, it inhibits the formation of final CO products. This mechanism of large C- and O-rich clusters as reaction retardants is akin to a recently proposed mechanism in another HE crystal, TATB. Manaa et al. [61] found that N-rich clusters impede the formation of the final product of N 2 and solid carbon. Similarly, in our simulation of RDX the persistent formation of large C- and O-rich clusters prohibit the formation of CO, which occur at different stages during the simulation. Carbon-rich clusters form immediately once the simulation starts. On the other hand, oxygen-rich clusters only form after the expansion phase; see Figure 6.8(b). It can be seen in Figure 6.8(b) that originally N-rich clusters also become N-poor during the expansion of the system. The change in composition slows down after 40 ps, consistent 130 with the slow decay time (~ 1 ns) of the large clusters in Figure 6.8(a). In order to study the configuration of large retardant clusters, which inhibit the production of the final products, Figure 6.8(c) compare the number of large C- and O-rich clusters (red line) with that of all large clusters (blue line) as a function of time. As in the main text, a large cluster is defined to have the number of constituent atoms, n > 21. Though the stoichiometry varies from cluster to cluster, we see that 60% of the large clusters are C- and O-rich. 131 Figure 6.8 Time evolution of large clusters. (a) Fraction of the number of atoms in larger cluster N C vs. the total number of atoms in the system N T . (b) Stoichiometric composition of the large clusters as a function of time. The dashed lines show the normalized stoichiometric number for different elements in the RDX molecule, C 3 H 6 N 6 O 6 . The 132 arrows indicate how the stoichiometric value changed at the final state. (c) Number of large clusters and large C- and O-rich clusters as a function of time. 6.3.7 Quantum mechanical validation In order to validate the metastability of the large C- and O-rich clusters in the high-temperature detonation condition, we perform QMD simulations, in which interatomic forces are computed quantum mechanically in the framework of density functional theory. We pick one of the carbon-rich clusters (C 12 H 20 O 13 N 13 ) from the population of large clusters in the ReaxFF-MD simulation, and perform QMD simulations on them using the VASP software package [186] at a temperature of 1,300 K in the canonical ensemble for 0.5 ps. In the simulations, electronic states are calculated using the projector-augmented-wave method [187], which is an all-electron electronic- structure-calculation method within the frozen-core approximation. The generalized gradient approximation [188] is used for the exchange-correlation energy. The plan-wave cutoff energy is set as 400 eV. Within the time period of the QMD simulation, the selected fragment remains stable. 133 Chapter 7. Conclusion Reactive MD simulations results have provided quantitative information on the combustion chemistry and the atomistic mechanism of detonation for aluminum nanostructures and high explosives, respectively. For Al-NPs with different sizes, D = 26, 36 and 46nm in which the Al core diameter is varied from 20, 30 and 40 nm and the alumina shell thickness is held constant at 3 nm. In the initial state of the system the nanoparticles are uniformly heated to 1,100 K in ambient oxygen environment. We have found that the oxidation dynamics starts within the Al-NP at the core-shell interface when molten aluminum starts reacting with alumina shell by taking some oxygen from expanded alumina shell to form aluminum rich Al x O y (x>y) clusters. This reaction at the aluminum/alumina interface generates heat first. By analyzing the structure of the oxide shell of all three Al-NPs in the process of this reaction, we found that the shell of Al-NPs does not break or shatter, but only deforms, where the deformations differ slightly depending on the size of the Al-NP. For smaller Al-NP, which has larger surface-to-volume ratio and lower ratio of core-radius- to-shell-thickness, the reaction is faster. In early stage of the reaction, majority of Al x O y clusters formed are small, having less than six atoms. Most of the clusters in this stage are AlO and Al 2 O as shown in Figure 3.10, Figure 3.11 and Figure 3.12. The induced heat is mainly as a result of the oxidation dynamics at the aluminum/alumina interface. When the temperature reaches ~ 2,900 K, the alumina shell starts to boil. This shell boiling is coincident with the start of aluminum ejections from the Al-NP into oxygen environment. 134 The aluminum ejections start earlier in smaller Al-NP – at t 0 = 0.18, 0.28 and 0.42 ns for D = 26, 36 and 46 nm. We would like to emphasize here that in all the three systems the shell temperature at t 0 = 0.18, 0.28 and 0.42 ns is 2,900 K at which time the shell is in molten state. After the Al ejection starts from Al-NP, the oxidation also starts at the alumina/oxygen interface. The initial Al x O y clusters found in this outer region are oxygen rich (y > x). The start time of AlO 2 clusters production in the outer regions of the Al-NP is after the shell temperature has reached ~2,900K, as can be inferred from Figure 3.10 (c), Figure 3.11(c), and Figure 3.12(c). We have also examined the rate of heat production in the three Al-NP systems. As the oxidation proceeds, the total system temperature (including the environmental oxygen) increases monotonically, however the time derivative of the total system temperature, (dT system /dt), has a peak that occurs at t 1 = 0.20, 0.32 and 0.51 ns for D = 26, 36 and 46 nm. At this peak in the time derivative of the total temperature, in (dT system /dt) , the shell temperature for the three Al-NPs are 3,100 K, 3,300 K, and 3,500 K, respectively, as given in the Table 3.2. However, because all the oxygen is not uniformly heated due to continuing oxidation dynamics inside and outside the Al-NP, the full system temperature is lower than the shell temperature. At the peak in (dT system /dt) the system temperature for the three Al-NP systems are 2,500 K, 2,700 K, and 2,800 K, respectively. The time lag between when the shell temperature reaches ~2,900 K and the peak in time derivative of the total system temperature, (dT system /dt) specifies maximum rate of heat production are t 1 – t 0 = 0.02, 0.04 and 0.09 ns for D = 26, 36 and 46 nm, clearly indicating the Al-NP size effect. 135 For Al-NP aggregates, multimillion-atom reactive MD simulations revealed that the oxidation dynamics is governed by the interplay between the agglomeration of Al- NPs and the production of small intermediate aluminum-oxide fragments such as Al 2 O, AlO, and AlO 2 . Smaller diameter Al-NP aggregates were found to oxidize faster than larger diameter ones, resulting in more rapid agglomeration, core Al consuming and shell oxidation. Major fragments at the end of simulation were AlO and AlO 2 , which can directly form the final product, Al 2 O 3 , by the reaction, AlO + AlO 2 → Al 2 O 3 . Our simulation results show that the oxide shell of Al-NP is not inert container of Al fuel, but rather plays an unexpectedly active role in the oxidation mechanism. Such atomistic understanding of the reaction pathways and the size effect should have profound impacts on the rational design for the efficient burning of metastable intermolecular composites. For Al-NRs, MD simulations of this oxide-coated aluminum nanostructures with different diameters reveal the atomistic mechanisms behind the faster oxidation reactions of the smaller Al-NRs. When one end of an Al-NR is heated to 1,100 K and the rest of the nanorod is kept at 300 K in oxygen environment, the nanorod burns from the heated end to the other end along the nanorod axis. The reaction initiates at the interface between the alumina shell and Al core, where shell O atoms react with core Al at high temperature to form bonded interfaces. Subsequently, the fuel-rich oxide shell starts absorbing environmental oxygen for further oxidation. This in turn accelerates the heat release. The smaller nanorods have the faster oxidation-front speeds. The fastest oxidation-front speed in our simulations is 88 m/s for D = 26 nm Al-NR. In front of the oxidation, the Al core melts with heat provided by the oxidation reaction. The melting speed has been found to 136 be faster than the oxidation speed for every studied nanorod. For D = 26 nm Al-NR, for instance, the melting-front speed is 91 m/s. The faster oxidation-front and melting-front speeds are due to the larger surface-to-volume ratios. Namely, heat produced from the shell interfaces promotes the melting and oxidation of the core Al volume. Thus enlarging the surface-to-volume ratio is critical for enhancing the burning efficiency. These findings might help understand the combustion behaviors of energetic compounds containing Al-NRs. Last but not least, ReaxFF-MD simulations results have provided the first evidence that the detonation of RDX proceeds in two stages: First, N 2 and H 2 O (and OH) are formed rapidly within ~ 10 ps, resulting from uni- and inter-molecular reaction pathways, respectively. Subsequently, CO starts to form at a very low rate (i.e., the reaction time ~ 1 ns). We found that the CO production is retarded by the formation of large C-rich and O-rich clusters. Such atomistic understanding of the reaction time and intermediate products provides valuable insight into broad technologies involving HEs. An example is rational design of insensitive energetic materials and detonation synthesis of materials such as nanodiamond [58]. Our simulations on energetic materials have provided valuable quantitative information on mechanical, chemical, structural and dynamics reactions at atomistic levels. Our results can be useful in designing better experiments to obtain comprehensive understanding and underlying mechanisms for energetic materials. 137 References 1. Kelly, J., Gunpowder: Alchemy, Bombards, and Pyrotechnics: The History of the Explosive that Changed the World. 2005: p. 2-5. 2. Londoño, E., U.S. Troops in Iraq Face A Powerful New Weapon: Use of Rocket- Propelled Bombs Spreads. Washington Post., 2008: p. A01. 3. Dillow, C., How it Would Work: Destroying an Incoming Killer Asteroid With a Nuclear Blast. Popular Science (Bonnier). 2013. 4. Timothee L. Pourpoint, T.D.W., Mark A. Pfeil, John Tsohas, and Steven F. Son, Feasibility Study and Demonstration of an Aluminum and Ice Solid Propellant. International Journal of Aerospace Engineering, 2012. 2012: p. 11. 5. Beloni, E., V.K. Hoffmann, and E.L. Dreizin, Combustion of Decane-Based Slurries with Metallic Fuel Additives. Journal of Propulsion and Power, 2008. 24(6): p. 1403-1411. 6. Davis, W.C., High explosives the interaction of chemistry and mechanics. p. 48- 75. 7. Bdzil, J.B., Understanding the effects of a finite-length reaction zone. Los Alamos Science, 2003. 28(High-Explosives Performance): p. 96-110. 8. Fischer, S.H. and M.C. Grubelich, Theoretical energy release of thermites, intermetallics, and combustible metals. Twenty-Fourth International Pyrotechnics Seminar, 1998: p. 231-286. 9. Rossi, C., et al., Nanoenergetic materials for MEMS: A review. Journal of Microelectromechanical Systems, 2007. 16(4): p. 919-931. 10. Aumann, C.E., G.L. Skofronick, and J.A. Martin, Oxidation Behavior of Aluminum Nanopowders. Journal of Vacuum Science & Technology B, 1995. 13(3): p. 1178-1183. 11. Dlott, D.D., Thinking big (and small) about energetic materials. Materials Science and Technology, 2006. 22(4): p. 463-473. 12. Dreizin, E.L., Metal-based reactive nanomaterials. Progress in Energy and Combustion Science, 2009. 35(2): p. 141-167. 13. Armstrong, R.W., et al., Enhanced propellant combustion with nanoparticles. Nano Letters, 2003. 3(2): p. 253-255. 14. Olsen, S.E. and M.W. Beckstead, Burn time measurements of single aluminum particles in steam and CO2 mixtures. Journal of Propulsion and Power, 1996. 12(4): p. 662-671. 15. Bazyn, T., H. Krier, and N. Glumac, Evidence for the transition from the diffusion-limit in aluminum particle combustion. Proceedings of the Combustion Institute, 2007. 31: p. 2021-2028. 138 16. Marion, M., C. Chauveau, and I. Gokalp, Studies on the ignition and burning of levitated aluminum particles. Combustion Science and Technology, 1996. 115(4- 6): p. 369-&. 17. Yetter, R.A., G.A. Risha, and S.F. Son, Metal particle combustion and nanotechnology. Proceedings of the Combustion Institute, 2009. 32: p. 1819- 1838. 18. Mench, M.M., et al., Comparison of thermal behavior of regular and ultra-fine aluminum powders (Alex) made from plasma explosion process. Combustion Science and Technology, 1998. 135(1-6): p. 269-292. 19. Ohkura, Y., et al., Reducing minimum flash ignition energy of Al microparticles by addition of WO3 nanoparticles. Applied Physics Letters, 2013. 102(4): p. 043108. 20. Rai, A., et al., Understanding the mechanism of aluminium nanoparticle oxidation. Combustion Theory and Modelling, 2006. 10(5): p. 843-859. 21. Son, S.F., Performance and characterization of nanoenergetic materials at Los Alamos. Synthesis, Characterization and Properties of Energetic/Reactive Nanomaterials, 2004. 800: p. 161-172. 22. Wu, H.C., et al., Explosion Characteristics of Aluminum Nanopowders. Aerosol and Air Quality Research, 2010. 10(1): p. 38-42. 23. Bockmon, B.S., et al., Combustion velocities and propagation mechanisms of metastable interstitial composites. Journal of Applied Physics, 2005. 98(6): p. 064903. 24. Levitas, V.I., Burn time of aluminum nanoparticles: Strong effect of the heating rate and melt-dispersion mechanism. Combustion and Flame, 2009. 156(2): p. 543-546. 25. Park, K., et al., Size-resolved kinetic measurements of aluminum nanoparticle oxidation with single particle mass spectrometry. Journal of Physical Chemistry B, 2005. 109(15): p. 7290-7299. 26. Gennady V. Ivanov, F.T., 'Activated' aluminum as a stored energy source for propellants International Journal of Energetic Materials and Chemical Propulsion, 1997. 4(1-6): p. 10. 27. Trunov, M.A., et al., Effect of polymorphic phase transformations in Al2O3 film on oxidation kinetics of aluminum powders. Combustion and Flame, 2005. 140(4): p. 310-318. 28. Huang, Y., et al., Effect of particle size on combustion of aluminum particle dust in air. Combustion and Flame, 2009. 156(1): p. 5-13. 29. Yu. A. Kotov, I.V.B., A. I. Medvedev, and O. R. Timoshenkova, Synthesizing Aluminum Nanoparticles in an Oxide Shell. Nanotechnologies in Russia, 2009. 4(5-6): p. 354-358. 139 30. Levitas, V.I., M.L. Pantoya, and B. Dikici, Melt dispersion versus diffusive oxidation mechanism for aluminum nanoparticles: Critical experiments and controlling parameters. Applied Physics Letters, 2008. 92(1): p. 011921. 31. Levitas, V.I., et al., Melt dispersion mechanism for fast reaction of nanothermites. Applied Physics Letters, 2006. 89(7): p. 071909. 32. Ohkura, Y., P.M. Rao, and X.L. Zheng, Flash ignition of Al nanoparticles: Mechanism and applications. Combustion and Flame, 2011. 158(12): p. 2544- 2548. 33. Rashkovskii, S.A., Effect of acceleration on agglomeration of aluminum particles during combustion of composite solid propellants. Combustion Explosion and Shock Waves, 2007. 43(6): p. 654-663. 34. Il'in, A.P., et al., Combustion of agglomerated ultrafine aluminum powders in air. Combustion Explosion and Shock Waves, 2002. 38(6): p. 665-669. 35. Grigorev, V.G., V.E. Zarko, and K.P. Kutsenogii, Experimental Investigation of the Agglomeration of Aluminum Particles in Burning Condensed Systems. Combustion Explosion and Shock Waves, 1981. 17(3): p. 245-251. 36. Karasev, V.V., et al., Formation of charged aggregates of Al2O3 nanoparticles by combustion of aluminum droplets in air. Combustion and Flame, 2004. 138(1- 2): p. 40-54. 37. J. Eckert, J.C.H., C.C. Ahn, Z. Fu, W.L. Johnson, Melting behavior of nanocrystalline aluminum powders. Nanostructured Materials, 1993. 2(4): p. 407- 413. 38. Price, E.W., Sigman, R. K., Sambamurthi, J. K. and Park, C. J., Behavior of Aluminum in Solid Propellant Combustion. Scientific Report, 1982: p. 91. 39. Mullen, J.C. and M.Q. Brewster, Reduced Agglomeration of Aluminum in Wide- Distribution Composite Propellants. Journal of Propulsion and Power, 2011. 27(3): p. 650-661. 40. Malchi, J.Y., et al., Nano-aluminum flame spread with fingering combustion instabilities. Proceedings of the Combustion Institute, 2007. 31: p. 2617-2624. 41. Pivkina, A., et al., Nanomaterials for heterogeneous combustion. Propellants Explosives Pyrotechnics, 2004. 29(1): p. 39-48. 42. Weismiller, M.R., et al., Effects of fuel and oxidizer particle dimensions on the propagation of aluminum containing thermites. Proceedings of the Combustion Institute, 2011. 33: p. 1989-1996. 43. Bezmelnitsyn, A., et al., Modified Nanoenergetic Composites with Tunable Combustion Characteristics for Propellant Applications. Propellants Explosives Pyrotechnics, 2010. 35(4): p. 384-394. 44. Prakash, A., A.V. McCormick, and M.R. Zachariah, Aero-sol-gel synthesis of nanoporous iron-oxide particles: A potential oxidizer for nanoenergetic materials. Chemistry of Materials, 2004. 16(8): p. 1466-1471. 140 45. Shende, R., et al., Nanoenergetic composites of CuO nanorods, nanowires, and Al-nanoparticles. Propellants Explosives Pyrotechnics, 2008. 33(2): p. 122-130. 46. Subramaniam, S., et al., Self-assembled nanoenergetic composite. Multifunctional Energetic Materials, 2006. 896: p. 9-14. 47. Perre, E., et al., Direct electrodeposition of aluminium nano-rods. Electrochemistry Communications, 2008. 10(10): p. 1467-1470. 48. Pang, Y.T., et al., Synthesis of ordered Al nanowire arrays. Solid State Sciences, 2003. 5(7): p. 1063-1067. 49. Pomfret, M.B., et al., Electrochemical Template Deposition of Aluminum Nanorods Using Ionic Liquids. Chemistry of Materials, 2008. 20(19): p. 5945- 5947. 50. Li, C.S., et al., Metallic aluminum nanorods: Synthesis via vapor-deposition and applications in al/air batteries. Chemistry of Materials, 2007. 19(24): p. 5812- 5814. 51. Saka, M. and R. Nakanishi, Fabrication of Al thin wire by utilizing controlled accumulation of atoms due to electromigration. Materials Letters, 2006. 60(17- 18): p. 2129-2131. 52. Lu, Y.B., et al., Enhancement of Al thin wire fabrication by using electromigration in relation to the discharge resistance of the atoms. Optoelectronics and Advanced Materials-Rapid Communications, 2011. 5(11): p. 1219-1222. 53. Langhammer, C., et al., Localized surface plasmon resonances in aluminum nanodisks. Nano Letters, 2008. 8(5): p. 1461-1471. 54. Cheah, S.K., et al., Self-Supported Three-Dimensional Nanoelectrodes for Microbattery Applications. Nano Letters, 2009. 9(9): p. 3230-3233. 55. Apodaca, L.E., http://minerals.usgs.gov/minerals/pubs/commodity/explosives/myb1-2012- explo.pdf. United States Geological Survey 2005 Minerals Yearbook, 2011. 56. Reed, E.J., et al., A transient semimetallic layer in detonating nitromethane. Nature Physics, 2008. 4(1): p. 72-76. 57. Wu, C.J., et al., Catalytic behaviour of dense hot water. Nature Chemistry, 2009. 1(1): p. 57-62. 58. Pichot, V., et al., Understanding ultrafine nanodiamond formation using nanostructured explosives. Scientific Reports, 2013. 3: p. 2159. 59. Rice, B.M., et al., Molecular-dynamics study of detonation .1. A comparison with hydrodynamic predictions. Physical Review E, 1996. 53(1): p. 611-622. 60. Brenner, D.W., et al., Detonations at Nanometer Resolution Using Molecular- Dynamics. Physical Review Letters, 1993. 70(14): p. 2174-2177. 141 61. Manaa, M.R., et al., Nitrogen-rich heterocycles as reactivity retardants in shocked insensitive explosives. Journal of the American Chemical Society, 2009. 131(15): p. 5483-5487. 62. Losada, M. and S. Chaudhuri, Theoretical study of elementary steps in the reactions between aluminum and teflon fragments under combustive environments. Journal of Physical Chemistry A, 2009. 113(20): p. 5933-5941. 63. Losada, M. and S. Chaudhuri, Finite size effects on aluminum/Teflon reaction channels under combustive environment: A Rice-Ramsperger-Kassel-Marcus and transition state theory study of fluorination. Journal of Chemical Physics, 2010. 133(13): p. 134305. 64. Shaw, M.S. and J.D. Johnson, Carbon clustering in detonations. Journal of Applied Physics, 1987. 62(5): p. 2080-2085. 65. Zhang, L.Z., et al., Carbon cluster formation during thermal decomposition of octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine and 1,3,5-triamino-2,4,6- trinitrobenzene high explosives from ReaxFF reactive molecular dynamics simulations. Journal of Physical Chemistry A, 2009. 113(40): p. 10619-10640. 66. Yetter, R.A., et al., Development of Gas-Phase Reaction-Mechanisms for Nitramine Combustion. Journal of Propulsion and Power, 1995. 11(4): p. 683- 697. 67. Patterson, J.E., Z.A. Dreger, and Y.M. Gupta, Shock wave-induced phase transition in RDX single crystals. Journal of Physical Chemistry B, 2007. 111(37): p. 10897-10904. 68. Hudson, R.J., P. Zioupos, and P.P. Gill, Investigating the Mechanical Properties of RDX Crystals Using Nano-Indentation. Propellants Explosives Pyrotechnics, 2012. 37(2): p. 191-197. 69. Li, T., C. Hua, and Q. Li, Shock Sensitivity of Pressed RDX-Based Plastic Bonded Explosives under Short-Duration and High-Pressure Impact Tests. Propellants Explosives Pyrotechnics, 2013. 35: p. 1-5. 70. Strachan, A., et al., Shock waves in high-energy materials: The initial chemical events in nitramine RDX. Physical Review Letters, 2003. 91(9): p. 098301. 71. Nomura, K.I., et al., Dynamic transition in the structure of an energetic crystal during chemical reactions at shock front prior to detonation. Physical Review Letters, 2007. 99(14): p. 148303. 72. Chakraborty, D., et al., The mechanism for unimolecular decomposition of RDX (1,3,5-trinitro-1,3,5-triazine), an ab initio study. Journal of Physical Chemistry A, 2000. 104(11): p. 2261-2272. 73. Zel’dovich, Y., Theory of Detonation. Theory of Detonation, 1960(Academic Press): p. 208-210. 74. Robertson, A.J.B., The Thermal Decomposition of Explosives .2. Cyclotrimethylenetrinitramine and Cyclotetramethylenetetranitramine. Transactions of the Faraday Society, 1949. 45(1): p. 85-93. 142 75. Duff, R.E. and E. Houston, Measurement of the Chapman-Jouguet pressure and reaction zone length in a detonating high explosive. Journal of Chemical Physics, 1955. 23(7): p. 1268-1271. 76. Lubyatinsky, S.N.L., B. G., Reaction zone measurements in detonating aluminized explosives. AIP Conference Proceedings, 1996. 370: p. 779. 77. Ershov, A.P. and N.P. Satonkina, Electrical conductivity distributions in detonating low-density explosives - Grain size effect. Combustion and Flame, 2010. 157(5): p. 1022-1026. 78. Ermolin, N.E., et al., Processes in Hexogene Flames. Combustion Explosion and Shock Waves, 1988. 24(4): p. 400-406. 79. Davidson, J.E. and M.W. Beckstead, Improvements to steady-state combustion modeling of cyclotrimethylenetrinitramine. Journal of Propulsion and Power, 1997. 13(3): p. 375-383. 80. Crowhurst, J.C., et al., Invariance of the Dissipative Action at Ultrahigh Strain Rates Above the Strong Shock Threshold. Physical Review Letters, 2011. 107(14). 81. Graham, R.A., Solids under high-pressure shock compression. Solids under high- pressure shock compression, 1993. 82. Hare, D.E., J. Franken, and D.D. Dlott, A New Method for Studying Picosecond Dynamics of Shocked Solids - Application to Crystalline Energetic Materials. Chemical Physics Letters, 1995. 244(3-4): p. 224-230. 83. Dang, N.C., et al., Femtosecond Stimulated Raman Scattering Picosecond Molecular Thermometry in Condensed Phases. Physical Review Letters, 2011. 107(4). 84. Tokmakoff, A., M.D. Fayer, and D.D. Dlott, Chemical-Reaction Initiation and Hot-Spot Formation in Shocked Energetic Molecular Materials. Journal of Physical Chemistry, 1993. 97(9): p. 1901-1913. 85. Partnership, M.a.M., Final Properties Report: Newport Army Ammunition Plant. National Park Service, AD-A, August 1984. 175(818). 86. Urbański, T., Laverton, Silvia, ed., Chemistry and Technology of Explosives III, translated by Marian Jureck (First English ed.). Warszawa: PWN - Polish Scientific Publishers and Pergamon Press. 87. http://en.wikipedia.org/wiki/Table_of_explosive_detonation_velocities. 88. Elderfield, R.C., Werner Emanual Bachmann: 1901–1951. Washington DC: National Academy of Sciences, 1960. 89. Baxter III, J.P., Scientists Against Time. Cambridge, MA: MIT Press, 1968. 90. Codina, R., Comparison of some finite element methods for solving the diffusion- convection-reaction equation. Computer Methods in Applied Mechanics and Engineering, 1998. 156(1-4): p. 185-210. 143 91. Galli, G., Linear scaling methods for electronic structure calculations and quantum molecular dynamics simulations. Current Opinion in Solid State & Materials Science, 1996. 1(6): p. 864-874. 92. Streitz, F.H. and J.W. Mintmire, Energetics of aluminum vacancies in gamma alumina. Physical Review B, 1999. 60(2): p. 773-777. 93. Dlott, D.D., Picosecond Dynamics Behind the Shock Front. Journal De Physique Iv, 1995. 5(C4): p. 337-343. 94. Strachan, A., A.C.T. van Duin, and W.A. Goddard, Initial chemical events in the energetic material RDX under shock loading: Role of defects. Shock Compression of Condensed Matter - 2003, Pts 1 and 2, Proceedings, 2004. 706: p. 895-898. 95. Dror, R.O., et al., Exploring atomic resolution physiology on a femtosecond to millisecond timescale using molecular dynamics simulations. Journal of General Physiology, 2010. 135(6): p. 555-562. 96. Stillinger, F.H., and A. Rahman, Molecular dynamics calculation of neutron inelastic scattering from water. Molecular Motions in Liquids, Springer Netherlands, 1974: p. 479-494. 97. Axelsen, P.H., E. Gratton, and F.G. Prendergast, Experimentally Verifying Molecular-Dynamics Simulations through Fluorescence Anisotropy Measurements. Biochemistry, 1991. 30(5): p. 1173-1179. 98. Mantsch, H.H., H. Saito, and I.C.P. Smith, Deuterium Magnetic-Resonance, Applications in Chemistry, Physics and Biology. Progress in Nuclear Magnetic Resonance Spectroscopy, 1977. 11: p. 211-271. 99. Alder, B.J. and T.E. Wainwright, Studies in Molecular Dynamics .1. General Method. Journal of Chemical Physics, 1959. 31(2): p. 459-466. 100. Rahman, A., Correlations in Motion of Atoms in Liquid Argon. Physical Review a-General Physics, 1964. 136(2A): p. A405-&. 101. Jones, J.E., On the determination of molecular fields III - From crystal measurements and kinetic theory data. Proceedings of the Royal Society of London Series a-Containing Papers of a Mathematical and Physical Character, 1924. 106(740): p. 709-718. 102. Vashishta, P., et al., Interaction potentials for alumina and molecular dynamics simulations of amorphous and liquid alumina. Journal of Applied Physics, 2008. 103(8): p. 083504. 103. Voter, A.F.C., S. P., Accurate Interatomic Potentials for Ni, Al and Ni3Al. Mat Res Soc Porc, 1987. 82: p. 175-180. 104. Wang, W., Thermal properties of silicon carbide and combustion mechanisms of aluminum nanoparticle. Department of Chemical Engineering and Materials Science, 2008. Univeristy of Southern California. 105. Politzer, P., P. Lane, and M.E. Grice, Energetics of aluminum combustion. Journal of Physical Chemistry A, 2001. 105(31): p. 7473-7480. 144 106. van Duin, A.C.T., et al., ReaxFF: A reactive force field for hydrocarbons. Journal of Physical Chemistry A, 2001. 105(41): p. 9396-9409. 107. Strachan, A., et al., Shock waves in high-energy materials: The initial chemical events in nitramine RDX. Physical Review Letters, 2003. 91(9). 108. Nomura, K., et al., A scalable parallel algorithm for large-scale reactive force- field molecular dynamics simulations. Computer Physics Communications, 2008. 178(2): p. 73-87. 109. Chenoweth, K., A.C.T. van Duin, and W.A. Goddard, ReaxFF reactive force field for molecular dynamics simulations of hydrocarbon oxidation. Journal of Physical Chemistry A, 2008. 112(5): p. 1040-1053. 110. Warrier, M., P. Pahari, and S. Chaturvedi, Validation Of A Reactive Force Field Included With An Open Source, Massively Parallel Code For Molecular Dynamics Simulations Of RDX. International Conference on Physics of Emerging Functional Materials (Pefm-2010), 2010. 1313: p. 278-280. 111. Allen, M.P., D.J. Tildesley, and J.R. Banavar, Computer simulation of liquids. Physics Today, 1989. 42: p. 71-80. 112. Sinclair, J., et al., Flexible boundary conditions and nonlinear geometric effects in atomic dislocation modeling. Journal of Applied physics, 1978. 49(7): p. 3890- 3897. 113. Nakano, A., R.K. Kalia, and P. Vashishta, Multiresolution molecular dynamics algorithm for realistic materials modeling on parallel computers. Computer Physics Communications, 1994. 83(2): p. 197-214. 114. Nakano, A., R.K. Kalia, and P. Vashishta, Scalable molecular-dynamics, visualization, and data-management algorithms for materials simulations. Computing in Science & Engineering, 1999. 1(5): p. 39-47. 115. Nomura, K., et al. A metascalable computing framework for large spatiotemporal-scale atomistic simulations. in Parallel & Distributed Processing, 2009. IPDPS 2009. IEEE International Symposium on. 2009. IEEE. 116. Dreizin, E.L., Experimental study of stages in aluminum particle combustion in air. Combustion and Flame, 1996. 105(4): p. 541-556. 117. Jordan, J.L., et al., Equation of state of aluminum-iron oxide-epoxy composite. Journal of Applied Physics, 2007. 101(9). 118. Brandstadt, K., D.L. Frost, and J.A. Kozinski, Preignition characteristics of nano- and micrometer-scale aluminum particles in Al-CO(2) oxidation systems. Proceedings of the Combustion Institute, 2009. 32: p. 1913-1919. 119. Ermoline, A., M. Schoenitz, and E.L. Dreizin, Reactions leading to ignition in fully dense nanocomposite Al-oxide systems. Combustion and Flame, 2011. 158(6): p. 1076-1083. 120. Stratakis, E., et al., Generation of Al nanoparticles via ablation of bulk Al in liquids with short laser pulses. Optics Express, 2009. 17(15): p. 12650-12659. 145 121. Bockmon, B.S., et al., Combustion velocities and propagation mechanisms of metastable interstitial composites. Journal of Applied Physics, 2005. 98(6). 122. Yan, Z.X., et al., Shock-induced thermal behavior of aluminum nanoparticles in propylene oxide. Journal of Applied Physics, 2007. 101(2). 123. Wang, S.F., et al., Dynamical effects of the oxide layer in aluminum nanoenergetic materials. Propellants Explosives Pyrotechnics, 2005. 30(2): p. 148-155. 124. Jones, M., et al., Experimental study of combustion characteristics of nanoscale metal and metal oxide additives in biofuel (ethanol). Nanoscale Research Letters, 2011. 6. 125. Gesner, J., M.L. Pantoya, and V.I. Levitas, Effect of oxide shell growth on nano- aluminum thermite propagation rates. Combustion and Flame, 2012. 159(11): p. 3448-3453. 126. Firmansyah, D.A., et al., Microstructural Behavior of the Alumina Shell and Aluminum Core Before and After Melting of Aluminum Nanoparticles. Journal of Physical Chemistry C, 2012. 116(1): p. 404-411. 127. Conner, R.W. and D.D. Dlott, Comparing Boron and Aluminum Nanoparticle Combustion in Teflon Using Ultrafast Emission Spectroscopy. Journal of Physical Chemistry C, 2012. 116(4): p. 2751-2760. 128. Lewis, W.K., et al., Comparison of post-detonation combustion in explosives incorporating aluminum nanoparticles: Influence of the passivation layer. Journal of Applied Physics, 2013. 113(4). 129. Huang, Y., et al., Combustion of bimodal nano/micron-sized aluminum particle dust in air. Proceedings of the Combustion Institute, 2007. 31: p. 2001-2009. 130. Kwon, Y.S., et al., The mechanism of combustion of superfine aluminum powders. Combustion and Flame, 2003. 133(4): p. 385-391. 131. Levitas, V.I., et al., Melt dispersion mechanism for fast reaction of nanothermites. Applied Physics Letters, 2006. 89(7). 132. Dikici, B., et al., Influence of Aluminum Passivation on the Reaction Mechanism: Flame Propagation Studies. Energy & Fuels, 2009. 23: p. 4231-4235. 133. De Luca, L.T., et al., Burning of nano-aluminized composite rocket propellants. Combustion Explosion and Shock Waves, 2005. 41(6): p. 680-692. 134. Jayaraman, K., et al., Production, Characterization, and Combustion of Nanoaluminum in Composite Solid Propellants. Journal of Propulsion and Power, 2009. 25(2): p. 471-481. 135. Morgan, A.B., et al., Heat release measurements on micron and nano-scale aluminum powders. Thermochimica Acta, 2009. 488(1-2): p. 1-9. 136. Pivkina, A.N., Y.V. Frolov, and D.A. Ivanov, Nanosized components of energetic systems: Structure, thermal behavior, and combustion. Combustion Explosion and Shock Waves, 2007. 43(1): p. 51-55. 146 137. Moore, K., M.L. Pantoya, and S.F. Son, Combustion behaviors resulting from bimodal aluminum size distributions in thermites. Journal of Propulsion and Power, 2007. 23(1): p. 181-185. 138. Yan, Z.X., J. Deng, and Z.M. Luo, A comparison study of the agglomeration mechanism of nano- and micrometer aluminum particles. Materials Characterization, 2010. 61(2): p. 198-205. 139. Bazyn, T., et al., Combustion Measurements of Fuel-Rich Aluminum and Molybdenum Oxide Nano-Composite Mixtures. Propellants Explosives Pyrotechnics, 2010. 35(2): p. 93-99. 140. Chung, S.W., et al., Size-dependent nanoparticle reaction enthalpy: Oxidation of aluminum nanoparticles. Journal of Physics and Chemistry of Solids, 2011. 72(6): p. 719-724. 141. Gan, Y.A. and L. Qiao, Combustion characteristics of fuel droplets with addition of nano and micron-sized aluminum particles. Combustion and Flame, 2011. 158(2): p. 354-368. 142. Campbell, T., et al., Dynamics of oxidation of aluminum nanoclusters using variable charge molecular-dynamics simulations on parallel computers. Physical Review Letters, 1999. 82(24): p. 4866-4869. 143. Campbell, T.J., et al., Oxidation of aluminum nanoclusters. Physical Review B, 2005. 71(20). 144. Alavi, S., J.W. Mintmire, and D.L. Thompson, Molecular dynamics simulations of the oxidation of aluminum nanoparticles. Journal of Physical Chemistry B, 2005. 109(1): p. 209-214. 145. Streitz, F.H. and J.W. Mintmire, Electrostatic Potentials for Metal-Oxide Surfaces and Interfaces. Interface Control of Electrical, Chemical, and Mechanical Properties, 1994. 318: p. 679-684. 146. Keffer, D.J. and J.W. Mintmire, Efficient parallel algorithms for molecular dynamics simulations using variable charge transfer electrostatic potentials. International Journal of Quantum Chemistry, 2000. 80(4-5): p. 733-742. 147. Puri, P. and V. Yang, Thermo-mechanical behavior of nano aluminum particles with oxide layers during melting. Journal of Nanoparticle Research, 2010. 12(8): p. 2989-3002. 148. Perron, A., et al., Oxidation of nanocrystalline aluminum by variable charge molecular dynamics. Journal of Physics and Chemistry of Solids, 2010. 71(2): p. 119-124. 149. Gutierrez, G., A. Taga, and B. Johansson, Theoretical structure determination of gamma-Al2O3. Physical Review B, 2002. 65(1). 150. Gutierrez, G. and B. Johansson, Molecular dynamics study of structural properties of amorphous Al2O3. Physical Review B, 2002. 65(10). 147 151. Vashishta, P., et al., Interaction potentials for alumina and molecular dynamics simulations of amorphous and liquid alumina. Journal of Applied Physics, 2008. 103(8). 152. Chen, S.P., A.F. Voter, and D.J. Srolovitz, A Simulation Study of Interfaces in Ni, Al and Ni3al with and without Boron. Journal De Physique, 1988. 49(C-5): p. 157-163. 153. Wang, W.Q., et al., Fast reaction mechanism of a core(Al)-shell (Al2O3) nanoparticle in oxygen. Applied Physics Letters, 2009. 95(26). 154. Wang, W.Q., et al., Effects of oxide-shell structures on the dynamics of oxidation of Al nanoparticles. Applied Physics Letters, 2010. 96(18). 155. Shekhar, A., et al., Collective oxidation behavior of aluminum nanoparticle aggregate. Applied Physics Letters, 2013. 102(22). 156. Hasnaoui, A., et al., Molecular dynamics simulations of the nano-scale room- temperature oxidation of aluminum single crystals. Surface Science, 2005. 579(1): p. 47-57. 157. Piehler, T.N., et al., Temporal evolution of the laser-induced breakdown spectroscopy spectrum of aluminum metal in different bath gases. Applied Optics, 2005. 44(18): p. 3654-3660. 158. Wang, S.F., et al., Fast spectroscopy of energy release in nanometric explosives. Chemical Physics Letters, 2003. 368(1-2): p. 189-194. 159. Pangilinan, G.I. and T.P. Russell, Role of Al-O-2 chemistry in the laser-induced vaporization of Al films in air. Journal of Chemical Physics, 1999. 111(2): p. 445- 448. 160. Prentice, J.L., Combustion of Laser-Ignited Aluminum Droplets - Rates and Mechanisms in Wet and Dry Oxidizers. Journal of the Electrochemical Society, 1975. 122(8): p. C260-C260. 161. Trunov, M.A., et al., Oxidation and melting of aluminum nanopowders. Journal of Physical Chemistry B, 2006. 110(26): p. 13094-13099. 162. Ananthapadmanabhan, P.V., et al., Formation of nano-sized alumina by in-flight oxidation of aluminium powder in a thermal plasma reactor. Scripta Materialia, 2004. 50(1): p. 143-147. 163. Campbell, T.J., et al., Oxidation of aluminum nanoclusters. Physical Review B, 2005. 71(20): p. 205413. 164. Wang, W., et al., Fast reaction mechanism of a core(Al)-shell (Al2O3) nanoparticle in oxygen. Applied Physics Letters, 2009. 95(26): p. 261901. 165. Wang, W., et al., Effects of oxide-shell structures on the dynamics of oxidation of Al nanoparticles. Applied Physics Letters, 2010. 96(18): p. 181906. 166. Li, Y., et al., Size effect on the oxidation of aluminum nanoparticle: Multimillion- atom reactive molecular dynamics simulations. Journal of Applied Physics, 2013. 114(13): p. 134312. 148 167. Shekhar, A., et al., Collective oxidation behavior of aluminum nanoparticle aggregate. Applied Physics Letters, 2013. 102(22): p. 221904. 168. Li, Y., et al., Size effect on the oxidation of aluminum nanoparticle: Multimillion- atom reactive molecular dynamics simulations. Journal of Applied Physics, 2013. 114(13). 169. See supplementary material for detailed simulation data. 170. Vashishta, P., et al., Crack propagation and fracture in ceramic films - Million atom molecular dynamics simulations on parallel computers. Materials Science and Engineering B-Solid State Materials for Advanced Technology, 1996. 37(1- 3): p. 56-71. 171. Cohen, J.M. and A.F. Voter, Convergence of Surface-Diffusion Parameters with Model Crystal Size. Surface Science, 1994. 313(3): p. 439-447. 172. Wang, W.Q., et al., Fast reaction mechanism of a core(Al)-shell (Al2O3) nanoparticle in oxygen. Applied Physics Letters, 2009. 95(26): p. 261901. 173. Wang, W.Q., et al., Effects of oxide-shell structures on the dynamics of oxidation of Al nanoparticles. Applied Physics Letters, 2010. 96(18): p. 181906. 174. Oh, S.H., et al., Oscillatory Mass Transport in Vapor-Liquid-Solid Growth of Sapphire Nanowires. Science, 2010. 330(6003): p. 489-493. 175. Laurent, V., et al., Wettability of Monocrystalline Alumina by Aluminum between Its Melting-Point and 1273-K. Acta Metallurgica, 1988. 36(7): p. 1797-1803. 176. Germann, T.K., K.; Lomdahl, P., 25 Tflop/s Multibillion-Atom Molecular Dynamics Simulations and Visualization/Analysis on BlueGene/L. Proceedings of IEEE/ACM Supercomputing ’05, 2005. 061: p. 59593. 177. Nomura, K.I., et al., A scalable parallel algorithm for large-scale reactive force- field molecular dynamics simulations. Computer Physics Communications, 2008. 178(2): p. 73-87. 178. Zhou, T.T., et al., Anisotropic shock sensitivity for beta-octahydro-1,3,5,7- tetranitro-1,3,5,7-tetrazocine energetic material under compressive-shear loading from ReaxFF-lg reactive dynamics simulations. Journal of Applied Physics, 2012. 111(12). 179. Borne, L., J. Mory, and F. Schlesser, Reduced sensitivity RDX (RS-RDX) in pressed formulations: Respective effects of intra-granular pores, extra-granular pores and pore sizes. Propellants Explosives Pyrotechnics, 2008. 33(1): p. 37-43. 180. Bourne, N.K., On the collapse of cavities. Shock Waves, 2002. 11(6): p. 447-455. 181. Liu, Z., S. Nagano, and S. Itoh, Overdriven detonation phenomenon in high explosive. Shock Compression of Condensed Matter-1999, Pts 1 and 2, 2000. 505: p. 227-230. 182. van Duin, A.C.T., et al., Atomistic-scale simulations of the initial chemical events in the thermal initiation of triacetonetriperoxide. Journal of the American Chemical Society, 2005. 127(31): p. 11053-11062. 149 183. Strachan, A., et al., Thermal decomposition of RDX from reactive molecular dynamics. Journal of Chemical Physics, 2005. 122(5). 184. Anisichkin, V.F., Mechanism of the detonation of TNT-RDX mixtures as studied by the isotope tracer method. Russian Journal of Physical Chemistry B, 2011. 5(2): p. 304-307. 185. Patterson, J.E., et al., Shock wave induced decomposition of RDX: Time-resolved spectroscopy. Journal of Physical Chemistry A, 2008. 112(32): p. 7374-7382. 186. Kresse, G. and J. Furthmuller, Efficient iterative schemes for ab initio total- energy calculations using a plane-wave basis set. Physical Review B, 1996. 54(16): p. 11169-11186. 187. Blochl, P.E., Projector Augmented-Wave Method. Physical Review B, 1994. 50(24): p. 17953-17979. 188. Perdew, J.P., K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple. Physical Review Letters, 1996. 77(18): p. 3865-3868.

Abstract (if available)

Abstract This dissertation is focused on molecular dynamics (MD) simulation of energetic materials. Energetic materials are a class of materials with high amount of stored chemical energy that can be released. Propellants, explosives, and pyrotechnic are typical classes of energetic materials. In the past decade, energetic materials have been playing very important roles in the fields of defense and other technologies, since they are secure, convenient and energetically compact, which will be continuously functionalized for further improvements indispensably in the new era. The significance of how to utilize the novel and promising features more profitably and constructively will necessitate fundamental understanding of dynamics of the material releasing energy through reactions. ❧ Among all kinds of energetic materials, aluminum nanostructures and high explosives are the subjects studied in this dissertation. For aluminum nanostructures, we have performed multimillion atom reactive MD simulations to investigate the oxidation dynamics of aluminum nanoparticles (Al-NPs), Al-NP aggregates, and aluminum nanorods (Al-NRs). For high explosives, we have addressed the shock‐induced detonation of cyclotrimethylenetrinitramine (RDX, an initialism for Research Department Explosive). ❧ The scientific interest in aluminum nanostructures research basically lies in the increased reactivity during oxidation as compared with conventional micro‐sized structures. However, what is the size effect on the reactivities of the Al-NPs still puzzles researchers. Practically, Al-NPs commercially produced always have a thin oxide layer encapsulating the pure metal particles. The oxide shell has been included in all our simulated aluminum nanostructures to accurately represent the reality as much as possible. We have systematically simulated different sizes of single Al-NPs (diameter, D = 26, 36 and 46 nm) to study their reactivities in oxygen environments as isolated individuals. We have also systematically extended our study to their assemblies—Al-NP aggregates with different sizes in oxygen environments. Aluminum nanorods (Al-NRs), a different type of morphology, also have been investigated due to their high contact areas with oxidizers. Our MD simulations reveal, within the range of our studied cases, the smaller sizes of aluminum nanostructures have higher reactivities (e.g., faster reaction rate, shorter reaction time and faster oxidation speed). Meanwhile, our atomistic analyses show interesting and unexpected roles of the interfaces of the oxide shell with aluminum and oxygen environment during the initial oxidation. In addition, the processes of final product formation show the evolution of different major intermediate species. Such atomistic understanding of the interface, reaction pathways and the size effect should have profound impacts on the design of the efficient burning of intermolecular products. ❧ The second topic of this dissertation is the detonation of a high explosive, RDX. Despite many years of intense research, atomistic mechanisms underlying the reaction time and intermediate reaction products of detonating high explosives have been elusive. Our reactive MD simulations reveal a novel two‐stage reaction mechanism during the detonation of RDX crystal. Rapid production of N₂ and H₂O is followed by delayed production of CO molecules. We have found that further decomposition towards the final products is inhibited by the formation of large, metastable carbon‐ and oxygen‐rich clusters. Such atomistic understanding may pave a way toward rational design for detonation synthesis of materials.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

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

PDF

Quantum molecular dynamics and machine learning for non-equilibrium processes in quantum materials

PDF

Molecular dynamics simulations of nanostructures

PDF

Large-scale molecular dynamics simulations of nano-structured materials

PDF

Shock sensitivity of energetic material and nanometric damage mechanisms in silica glass

$Multimillion-to-billion atom molecular dynamics simulations of deformation, damage, nanoindentation, and fracture in silica glass and energetic materials$

PDF

Multimillion-to-billion atom molecular dynamics simulations of deformation, damage, nanoindentation, and fracture in silica glass and energetic materials

PDF

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

PDF

Modeling of bicontinuous nanoporous materials and investigation of their mechanical and hydraulic properties

PDF

Towards groundbreaking green energetic materials

PDF

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

PDF

Ice: an impedance spectroscopy and atomistic simulation study

PDF

Controlling electronic properties of two-dimensional quantum materials: simulation at the nexus of the classical and quantum computing eras

PDF

Simulation and machine learning at exascale

PDF

Dynamics of wing cracks and nanoscale damage in silica glass

PDF

Shock-induced poration, cholesterol flip-flop and small interfering RNA transfection in a phospholipid membrane: multimillion atom, microsecond molecular dynamics simulations

PDF

Theoretical studies of lipid bilayer electroporation using molecular dynamics simulations

PDF

Neural network for molecular dynamics simulation and design of 2D materials

PDF

High throughput computational framework for synthesis and accelerated discovery of dielectric polymer materials using polarizable reactive molecular dynamics and graph neural networks

PDF

Quantum molecular dynamics simulations of hydrogen production and solar cells

PDF

Metascalable hybrid message-passing and multithreading algorithms for n-tuple computation

Asset Metadata

Creator Li, Ying (author)

Core Title Molecular dynamics study of energetic materials

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Materials Science

Publication Date 08/08/2016

Defense Date 06/03/2014

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

Tag computational materials science,Computer Science,energetic materials,OAI-PMH Harvest,Physics

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Vashishta, Priya (committee chair), Armani, Andrea M. (committee member), Gupta, Malancha (committee member), Nakano, Aiichiro (committee member)

Creator Email 9yingli9@gmail.com,yli12@usc.edu

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

Unique identifier UC11286892

Identifier etd-LiYing-2800.pdf (filename),usctheses-c3-458055 (legacy record id)

Legacy Identifier etd-LiYing-2800.pdf

Dmrecord 458055

Document Type Dissertation

Format application/pdf (imt)

Rights Li, Ying

Type texts

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

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA