University of Southern California Dissertations and Theses

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

(USC Thesis Other)

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

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content NANOMATERIALS UNDER EXTREME ENVIRONMENTS: A STUDY OF STRUCTURAL AND DYNAMIC PROPERTIES USING REACTIVE MOLECULAR DYNAMICS SIMULATIONS By Adarsh Shekhar _____________________________________________________________________ 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) December 2013 Copyright 2013 Adarsh Shekhar ii (Lead me from Darkness to Light) Dedicated to my family, friends, Vidushak and Internet iii Acknowledgements Being a graduate student in the Collaboratory for Advanced Computing and Simulations (CACS) has been a great learning experience for me. I would like to express my heartfelt thanks to my advisor Prof. Priya Vashishta for guiding through this extremely important phase of my professional life. He has trusted me with important projects and taught me how to plan, execute, analyze and present a good research work. He has been a wonderful advisor and it was a pleasure working with him. I would also like to thank Prof. Rajiv Kalia and Prof. Aiichiro Nakano for helping me understand and explore deeper into my research projects. They have always asked the right questions and had the right suggestions to improve and polish our results. Next, I would like to thank Dr. Ken-ichi Nomura for teaching me the nuances of molecular dynamics simulations. He has patiently helped me through the learning phase of my research work, and has been a great sounding board for any kind of questions I might have had. I would also thank Dr. Richard Clark and Dr. Weiqiang Wang for getting me started and helping me through my first research project. I take this opportunity to thank my colleagues Dr. Zaoshi Yuan, Dr. Manaschai Kunaseth, Dr. Richard Seymour, Dr. Weiwei Mou, Amit Choubey, Ying Li and Pankaj Rajak for all their help and support during my graduate studies. Having such great colleagues made my time as a graduate student a lot of fun. I also thank Ms. Patricia Wong for taking care of all my academic requirements. She has always been very understanding and extremely helpful. I take this opportunity to thank my friends Dr. Joyita Dutta, Dr. Sushmita Allam, Dr. Rumi Ghosh, Dr. Anustup Choudhury and Krishnakali Dasgupta for the wonderful years iv we spent together as graduate students. They have always motivated me, supported me in tough times and have been a great company to be with. I also thank all the current and former members of Vidushak, the Indian improv theater group at USC, for providing and sharing a great platform to explore my creativity, having a lot of fun and helping me relax throughout my stay at USC. I also thank my undergraduate advisor Prof. N. Chakraborti at Indian Institute of Technology, Kharagpur for introducing me to computational materials science and teaching me the basics of academic research. He has been instrumental in my progression from an undergraduate to graduate student. I would like to thank my early school teachers Late Mr. D. Mukherjee, Mr. K. K. Panda, Mr. M. K. Chand, Mr. M. Hussain, and Mr. A. K. Yadav for taking me in their wings when I was a toddler and teaching me to learn, adapt and succeed. They helped me understand the importance of education and prepared me for the competitive world. Last, but not the least, I would to express my sincerest gratitude towards my parents Mr. Upendra Singh and Mrs. Deepti Singh, and my sister Priyanka Shekhar for always believing in me. It is because of their unconditional love and support that I am able to be here today. I owe all my success and accomplishments to them and will be forever indebted to them. v Table of Contents Acknowledgements ............................................................................................................ iii List of Figures ................................................................................................................... vii Abstract ............................................................................................................................... 1 Chapter 1. Introduction ....................................................................................................... 3 1.1 Oxidation of Aluminum Nanoparticles ..................................................................... 3 1.2 Shock-induced Collapse of Water Nanobubble near Silica surface ......................... 4 1.3 Structural and Dynamic Properties of Nanoconfined Water .................................... 6 Chapter 2. Molecular Dynamics ......................................................................................... 9 2.1 Overview ................................................................................................................... 9 2.2 Interatomic Potential ............................................................................................... 11 2.3 Discretization .......................................................................................................... 12 2.4 Linked List Cell Decomposition ............................................................................. 16 2.5 Parallel Molecular Dynamics .................................................................................. 19 2.6 Molecular Dynamics on BlueGene/P ..................................................................... 23 Chapter 3. Oxidation of Aluminum Nanoparticle Aggregate ........................................... 24 3.1 Background ............................................................................................................. 24 3.2 Simulation ............................................................................................................... 36 3.2.1 Interatomic Potential .................................................................................... 38 3.3 Oxidation of outer ANPs ........................................................................................ 41 3.3.1 Radial distribution of temperature ............................................................... 42 3.3.2 Rate of change of temperature ..................................................................... 44 3.3.3 Aluminum Oxide clusters ............................................................................ 45 3.4 Oxidation of central ANP ....................................................................................... 50 3.4.1 Temperature profile ..................................................................................... 53 3.4.2 Temperature gradient ................................................................................... 56 3.4.3 Reaction Heat ............................................................................................... 58 3.4.4 Penetration into central ANP ....................................................................... 60 3.5 Summary ................................................................................................................. 63 Chapter 4. Shock-Induced Collapse of Water Nanobubble .............................................. 65 4.1 Background ............................................................................................................. 65 4.2 Simulation ............................................................................................................... 76 4.2.1 Preparation of Amorphous Silica ................................................................. 78 4.2.1 Interatomic Potential .................................................................................... 79 4.3 Collapse of empty nanobubble ............................................................................... 85 4.3.1 High-speed water nanojet ............................................................................ 87 4.3.2 Pressure profile in water during collapse ..................................................... 90 vi 4.3.3 Cavitation damage on the silica slab ............................................................ 92 4.3.4 Pressure profile in silica slab ....................................................................... 96 4.4 Collapse of gas-filled nanobubble .......................................................................... 98 4.4.1 Lack of high-speed nanojet .......................................................................... 99 4.4.2 Reduced mechanical damage ..................................................................... 103 4.4.3 Pressure profile .......................................................................................... 105 4.4.4 Reduced chemical damage ......................................................................... 107 4.5 Summary ............................................................................................................... 108 Chapter 5. Water confined in Nanoporous Silica ........................................................... 109 5.1 Background ........................................................................................................... 109 5.2 Simulation ............................................................................................................. 111 5.3 Results and Discussion ......................................................................................... 113 5.4 Summary ............................................................................................................... 124 Chapter 6. Conclusion ..................................................................................................... 125 References ....................................................................................................................... 132 vii List of Figures Figure 2.1: (a) Schematic of a sample MD system (b) A typical trajectory of an MD system. ........................................................................................................................ 9 Figure 2.2: Schematic of (a) Periodic Boundary Conditions and (b) Minimum Image Convention ................................................................................................................ 14 Figure 2.3: Lennard-Jones Potential (blue) and its shifted form (red) .............................. 16 Figure 2.4: (a) Schematic showing linked list cell decomposition. (b) The connected components in a sample linked list. .......................................................................... 18 Figure 2.5: Schematic showing spatial distribution of system among multiple processors ................................................................................................................................... 19 Figure 2.7: (a) Schematic showing communication deadlock. (b) Two-phase message passing scheme to avoid deadlock. ........................................................................... 21 Figure 3.1: The system used in the simulation. (a) A single ANP, consisting of a spherical aluminum core, 400 Å in diameter, surrounded by a 30 Å thick spherical alumina shell. The aluminum core is shown in grey, the aluminum atoms in the shell are shown in red while the oxygen atoms in the shell are shown in yellow. (b) A slice of the system through the central place showing the linear arrangement of the three ANPs. The aluminum and oxygen atoms in the shell of the hot ANPs (left and right) are shown in red and yellow respectively, while the aluminum and oxygen atoms in the shell of the cold ANP (central) are shown in red and yellow, respectively. The aluminum atoms in the cores of the three ANPs are shown in grey. ........................ 37 Figure 3.2: Snapshots of the central slice of left ANP at (a) 0 ns, (b) 0.1 ns, (c) 0.2 ns, and (d) 0.3 ns. The aluminum core is shown in grey, the aluminum atoms in the shell are shown in red and the oxygen atoms in the core are shown in yellow. The reactions begin at the core-shell interface, resulting in the diffusion of shell atoms into the aluminum core. By 0.3 ns, the outer ANPs expand and undergo considerable deformation. .............................................................................................................. 41 Figure 3.3: Radial distribution of temperature inside left ANP. The ANP is divided into 5 Å thick spherical bins and the temperature of atoms falling in each bin is calculated and plotted against the distance from the center of the ANP. The red line represents the initial state. The distribution at 0.2 ns is shown in blue, while the same at 0.4 ns is shown in green. The double-sided arrows indicate the increase in temperature by 0.2 ns. By 0.4 ns, the temperature across the ANP reaches 4000 K. ........................ 43 viii Figure 3.4: Temperature and rate of increase in temperature, dT/dt, vs time in outer ANPs. The blue curve represents temperature vs time, with the temperature axis on the right. The rate of increase in temperature (dT/dt) vs time is shown in red, with dT/dt axis on the left. The green dash-dot lines demarcate the three stages of oxidation, while the purple dashed line represents the melting temperature of alumina, T = 2300K. ................................................................................................. 45 Figure 3.5: Distribution of clusters of compositions equivalent to (a) Al2O, (b) AlO, (c) AlO2, (d) Al2O3, (e) Al2O5, and (f) AlO3 inside the left ANP with time. Al-rich clusters (a) and (b) are observed from the beginning while O-rich clusters (c) – (f) appear after 0.3 ns. Al-rich clusters are considerably more in number than O-rich clusters ...................................................................................................................... 46 Figure 3.6: Radial distribution of composition in outer ANPs during the simulation. Composition is calculated as percentage of oxygen in 5 Å thick spherical bins. %Oxygen is plotted against distance from the center of ANP. ................................. 48 Figure 3.7: Snapshots of the system at (a) 0 ns, (b) 0.3 ns, (c) 0.6 ns, and (d) 1 ns. Only the central slice of shells of the three ANPs are being shown. Aluminum atoms in the outer (left and right) ANPs are shown in red, oxygen atoms in outer ANPs are in yellow; aluminum atoms in the central ANP are in cyan and oxygen atoms in the central ANP are shown in blue. Penetration into the central ANP is observed and the system forms an ellipsoidal aggregate by 1 ns. ......................................................... 51 Figure 3.8: a) Aluminum oxide clusters inside the central ANP are grouped based on their positions. Region I (red) represents a right circular cone of aperture = 10° along Z- axis, which is perpendicular to the zones of contact, while region II (blue) represents that along the X-axis, which go through the zones of contact. b) Histogram showing %increase in clusters in central ANP inside region I (red) and region II (blue) at different times. More clusters are observed in the zones of contact as compared to other regions inside the central ANP. ....................................................................... 52 Figure 3.9: Temperature across the central slice of the system at (a) 0 ns, (b) 0.3 ns, (c) 0.6 ns, and (d) 1 ns. The scale is shown on the right with blue representing temperatures below 500 K and red representing temperatures above 6000 K. ......... 54 Figure 3.10: Linear distribution of temperature along the X-axis. Heat flows from the outer ANPs to the central ANP through the zones of contact. The core of the central ANP is the last region to be heated. .......................................................................... 56 Figure 3.11: Temperature gradient (dT/dt) (blue) and the number of ejections from the alumina shell (red) in (a) left and (b) central ANPs. The temperature gradient curve corresponds to the axis on the left while the ejections curve corresponds to the axis on the right. The onset of stage III, in both the cases, coincides with the beginning of ejections of atoms form the shell. ............................................................................. 58 ix Figure 3.12: Reaction heat per particle generated inside the central ANP vs time. Initially, the rate of generation of heat is slow as the reactions are localized to the zones of contact. As the outer ANPs penetrate the central ANP and initiate more oxidation reactions, thereby increasing the reaction rate between 0.4 and 0.6 ns. After 0.6 ns, the rate declines again as ejections of atoms from central ANP into the background oxygen begin and the central ANP enters stage III of oxidation. ............................. 59 Figure 3.13: Radial distribution of outer atoms inside the central ANP. After 40.4 ns (brown), the outer atoms are still close to the alumina shell; by 0.6 ns, (green) the outer atoms have penetrated half the radius of the central ANP and by 0.8 ns (red), the penetration reaches 80 Å. By 1 ns (blue), the outer atoms can be seen to be present at the core. .................................................................................................... 60 Figure 3.14: Position of different layers of the penetration front inside the central ANP vs time. The atoms from the shell of the central ANP form the first layer (brown and magenta), the atoms from the shell of the outer ANPs from the second layer (red and blue) and atoms from the core of the outer ANPs form the third and last layer (green). The penetration fronts maintain their relative position with respect to the core of central ANP. ................................................................................................. 61 Figure 3.15: Maximum position of outer atoms inside the central ANP vs time. The average speed of penetration is 54 m/s. .................................................................... 63 Figure 4.1: Schematic of the system used for simulation ................................................. 77 Figure 4.2: Melt-Quench procedure for preparing amorphous silica ............................... 78 Fig. 4.3: The number of neighboring silicon/hydrogen atoms around an oxygen atom is a function of the distance between the Si/H and the O atoms.(source: Ref #[128]) .... 84 Figure 4.4: Snapshots of the system during collapse of empty nanobubble at (a) t = 5 ps, (b) t = 22 ps, (c) t = 35 ps, and (d) t = 50 ps. Silica slab is shown in yellow and orange. A cross section of water across the system is shown. Normal density water is shown in green while water under shock is shown in blue. ...................................... 86 Figure 4.5: Position of the tip of water nanojet during the simulation. The position of the shock front away from the bubble is shown in grey. ................................................ 87 Figure 4.6: (a) Water nanojet formed during shock-induced collapse of the nanobubble. The silica slab is shown in yellow. The central slice across the silica slab is color coded by the magnitude of the velocity of water molecules. The velocity field lines show a focused jet of water molecules. (b) Number of H 3 O + species along the central slice of the system. H 3 O + ion production increases significantly when the nanojet hits the distal side of the bubble at time t = 25 ps (blue), and t = 5 ps (red). ............ 89 Figure 4.7: Pressure profile in water at (a) t = 15 ps, (b) t = 25 ps, (c) t = 30 ps, (d) t = 46 ps, (e) t = 48 ps, and (f) t = 50 ps .............................................................................. 90 x Figure 4.8: (a) Cavitation damage due to water nanojet. (b) Side-view of the cross-section of the system at t = 50 ps. Silica slab is shown in yellow while the water is shown in cyan (hydrogen) and blue (oxygen). ......................................................................... 92 Figure 4.9: Maximum depth of cavitation pit with time during collapse of bubble of initial diameter of 40 nm (red) and 97.6 nm (blue). .................................................. 93 Figure 4.10: Volume of the cavitation pit in systems with nanobubbles of diameter 40 nm (red) and 97.6 nm (blue). .......................................................................................... 94 Figure 4.11: (a) Snapshot shows atoms at the pit edge of the silica slab in the system with 10 9 atoms. Silicon-oxygen bonds are yellow-red, and cyan spheres represent hydrogen atoms inside the silica slab. Blue spheres represent oxygen atoms in water molecules within 0.5 nm from the silica slab. (b) Silanol distribution in the silica slab as a function of the distance from the pit center for bubbles of diameters 40 nm (red) and 97.6 nm (blue). .......................................................................................... 95 Figure 4.12: Pressure profile across the silica slab during collapse of empty nanobubble at (a) t = 35 ps, (b) t = 40 ps, (c) t = 45 ps, and (d) t = 50 ps ........................................ 97 Figure 4.13: Snapshots of the system during collapse of gas-filled nanobubble at (a) t = 5 ps, (b) t = 25 ps, (c) t = 35 ps, and (d) t = 50 ps. Silica slab is shown in yellow and orange. A cross section of water across the system is shown. Normal density water is shown in green while water under shock is shown in blue. The gas particles are shown in magenta. .................................................................................................... 98 Figure 4.14: Comparison of jet position during shock-induced collapse of empty and gas- filled nanobubble. ..................................................................................................... 99 Figure 4.15: Vortex rings around the deformed gas bubble during shock-induced collapse. The inset shows the direction of local velocity vectors (black arrows). The color gradient represents the magnitude of the local velocity. ............................... 100 Figure 4.16: H 3 O + concentration in water around the periphery of the collapsing gas- filled nanobubble .................................................................................................... 102 Figure 4.17: Side-view of the cross-section of the system with (a) empty and (b) gas- filled nanobubble at t = 50 ps. Silica slab is shown in yellow while the water is shown in cyan (hydrogen) and blue (oxygen). Gas particles are shown in red. ..... 103 Figure 4.18: Comparison of depth of cavitation pit formed due to the collapse of empty and gas-filled nanobubbles of two initial diameters. .............................................. 104 Figure 4.19: Pressure profiles across the silica slab and inside the gas-filled bubble during shock-induced collapse at (a) t = 35 ps, (b) t = 40 ps, (c) t = 45 ps, and (d) t = 50 ps. ................................................................................................................................. 106 xi Figure 4.20: Comparison of chemical damage in the cavitation pit due to collapse of empty and gas-filled nanobubbles .......................................................................... 107 Figure 5.1: Preparation schedule for Nanoporous silica from β-cristobalite. ................. 112 Figure 5.2: (a) Structure of amorphous silica network with density 0.5 g/cc (red). The empty space is color coded by the distance of the nearest atom. (b) Close-up view of the silica network (red) around the largest pore, shown in cyan. c) Distribution of pores with respect to anisotropy. Perfect sphere has anisotropy = 0. d) Distribution of pores with respect to shape parameter. Negative shape parameter corresponds to oblate spheroid, while the pores with positive shape parameter resemble prolate spheroid. .................................................................................................................. 114 Figure 5.3: (a) Heterogeneous distribution of water in nanoporous silica. Silica network is shown in red, while yellow, green and blue represent the local density of water in the network. (b) Distribution of local density of water averaged over the entire system at various temperatures. .............................................................................................. 117 Figure 5.4: Local strucutre in LDW and HDW .............................................................. 119 Figure 5.5: Mean square displacement of a water molecule in nanoconfined water. The molecule is seen to break from one cage (green) and move into another cage (red) in 0.2 ns. ...................................................................................................................... 121 Figure 5.6: Cage correlation function vs time at (a) T = 150 K, (b) T = 200 K, and (c) T = 250 K. (d) Cage correlation function in vs normalized time for various temperatures. Stretched exponential curves for all temperatures collapse onto a single stretched exponential curve, with exponent β = d/(d+2). ....................................................... 122 1 Abstract Nanotechnology is becoming increasingly important with the continuing advances in experimental techniques. As researchers around the world are trying to expand the current understanding of the behavior of materials at the atomistic scale, the limited resolution of equipment, both in terms of time and space, act as roadblocks to a comprehensive study. Numerical methods, in general and molecular dynamics, in particular act as able compliment to the experiments in our quest for understanding material behavior. In this research work, large scale molecular dynamics simulations to gain insight into the mechano-chemical behavior under extreme conditions of a variety of systems with many real world applications. The body of this work is divided into three parts, each covering a particular system: 1) Aggregates of aluminum nanoparticles are good solid fuel due to high flame propagation rates. Multi-million atom molecular dynamics simulations reveal the mechanism underlying higher reaction rate in a chain of aluminum nanoparticles as compared to an isolated nanoparticle. This is due to the penetration of hot atoms from reacting nanoparticles to an adjacent, unreacted nanoparticle, which brings in external heat and initiates exothermic oxidation reactions. 2) Cavitation bubbles readily occur in fluids subjected to rapid changes in pressure. We use billion-atom reactive molecular dynamics simulations on a 163,840-processor BlueGene/P supercomputer to investigate chemical and mechanical damages caused by 2 shock-induced collapse of nanobubbles in water near amorphous silica. Collapse of an empty nanobubble generates high-speed nanojet, resulting in the formation of a pit on the surface. The pit contains a large number of silanol groups and its volume is found to be directly proportional to the volume of the nanobubble. The gas-filled bubbles undergo partial collapse and consequently the damage on the silica surface is mitigated. 3) The structure and dynamics of water confined in nanoporous silica are different from that of bulk water, and insight into the properties of confined water is important for our understanding of many geological and biological processes. Nanoporous silica has a wide range of technological applications because it is easy to tune the size of pores and their morphologies and to functionalize pore surfaces with a variety of molecular moieties. Nanoporous silica is used in catalysis, chromatography, anticorrosion coatings, desalination membranes, and as drug delivery vehicles. We use reactive molecular dynamics to study the structure and dynamics of nanoconfined water between 100 and 300 K. 3 Chapter 1. Introduction 1.1 Oxidation of Aluminum Nanoparticles Aluminum nanoparticles (AlNPs) are widely used in energetic applications such as explosives,[1, 2] rocket propulsions[3-5] and other pyrotechniques[6] due their high energy release rate as compared to micron-sized particles.[7-9] Although there have been many experimental studies on AlNPs,[10-15] focusing on combustion mechanism[16], energy release rates,[13] and size effects,[17] atomistic level of understanding the oxidation behavior is difficult due to small spatio-temporal scale. It has been proposed[14, 18] that while AlNPs oxidize via diffusion through growing oxide shell during slow heating, they undergo drastically different mechanism during fast heating. Under this mechanism, called melt dispersion mechanism, the aluminum core melts resulting in spallation of alumina shell and subsequent ejection of aluminum clusters out of the AlNP.[13, 14, 16, 18, 19] However, this mechanism considers AlNPs as isolated entities and does not describe the overall effect on the entire aggregate. Though experimental images show that the nanoparticles are very closely placed,[20, 21] it remains unexplained how the coupling between proximate AlNPs affects the burning behavior. With the advent of computing techniques and ever increasing computing resources, researchers have used numerical methods[22] such as molecular dynamics (MD) simulations[23, 24] to gain insight into properties like heat transfer,[25, 26] flow behavior[27] and size effects[17] in nanoparticles. So far, most of the numerical studies have involved only a few thousand atoms. To study such systems more comprehensively, large-scale simulations involving millions of atoms are required. In our previous studies, 4 we have used multi-million atom MD simulations to understand the mechanism behind the combustion of single AlNP when subjected to different heating methods.[23, 24, 28] In this work, we extend these methods to study the heat and mass transfer when two hot AlNPs are placed in close proximity of an identical but cold AlNP, in order to gain insight into the mode and mechanism behind the flame propagation in densely packed nanoparticles. 1.2 Shock-induced Collapse of Water Nanobubble near Silica surface Erosion due to bubble collapse is commonly observed in pipes, turbines, pumps and ship propellers,[29-39] Cavitation erosion damage has also been observed in spallation neutron sources,[40] where rapidly deposited heat energy from the proton beam generates pressure waves and cavitation bubbles in mercury.[41, 42] Collapse of these bubbles can cause severe damage to the vessel wall.[43] Quantum mechanical calculations have revealed[44] that the structure of water near a silica surface is different from that of bulk. Classical molecular dynamics simulations have also been carried out to study water- induced damage in self-assembled monolayers on amorphous silica.[45] Recently, cavitation bubbles have been used beneficially in a number of technological applications. This includes an environmentally friendly approach to generate nanoporous or mesoporous material surfaces with acoustic cavitation.[46-48] In this approach, high- intensity focused ultrasound produces nano- or micro-bubbles. When bubbles collapse, high-speed nanojets or microjets of water are generated in an environment akin to that of a microreactor,[49] i.e., extremely high temperatures and pressures in a localized volume 5 while the rest of the system is under ambient conditions. The impact of high-speed jets and oxidation due to water sonolysis can make the metal surface porous at the nanometer scale. By controlling the ultrasound intensity and frequency and by varying the solvent and sonication time, it is possible to control the roughness and porosity of sponges on the surfaces of various metals and metallic alloys. Mesoporous sponges are excellent stimuli- responsive systems to encapsulate and deliver a variety of active compounds including corrosion inhibitors, antibodies, DNA fragments, enzymes, and biocides.[46] Ultrasound- assisted cavitation is also used in medical applications such as dentistry and extracorporeal shock-wave treatment of renal stones.[50] Ultrasound, with in-vivo gas bubbles, is used to enhance the permeability of cell membranes[51] and thus improve the efficiency of protein, DNA and drug delivery. Cavitation bubbles are also used beneficially in an industrial process called water jet peening[52] in nuclear reactors. Pressurized water is injected into water, generating cavitation bubbles. The pressure released during the bubble collapse mitigates the residual tensile stress and prevents stress corrosion cracking in austenite-based stainless steel and nickel based alloys used in nuclear reactors. Collapse of a cavitation bubble near a solid has been captured with high-speed cinematography,[53-55] but with limited spatial resolution (~ a few micrometers). Philipp and Lauterborn[56] have observed experimentally that the shock generated by the collapse of acoustic cavitation bubbles is strongly attenuated with the distance, H 0 , between the bubble and the solid surface, and that for maximum damage, H 0 should be 6 less than twice the initial bubble radius R 0 . Bubble collapse has also been studied[32] with continuum simulation approaches[57] and the impact damage of bubble collapse on solid surfaces has been mapped out as a function of the stand-off parameter, s d = H 0 /R 0 . These simulations also show that the pressure generated on the wall due to shock-induced bubble collapse is maximum for sd = 2. It has also been found that shock[58] loading can initiate chemical decomposition of materials.[59, 60] Despite a great deal of research on cavitation bubbles, several important questions about bubble collapse near a solid surface remain unanswered: Is there turbulence at the nanoscale during bubble collapse and if there is, what are its characteristics? What is the extent of mechanical damage due to nanojet impact and how does it scale with the bubble size? What is the nature of chemical damage? What is the connection between chemical and mechanical damages? These questions are examined through billion-atom reactive molecular dynamics (MD) simulations of shock-induced collapse of a nanobubble near a silica slab. 1.3 Structural and Dynamic Properties of Nanoconfined Water Understanding the behavior of water confined in ultra small spaces[61] is important in widely varying research fields such as life sciences, geosciences, atmospheric sciences, energy and environment.[62] For example, in the rapidly emerging DNA sequencing technology based on nanopore sensing,[63] one of the key unresolved issues is how the structure and dynamics of water confined in nanopores sculpted in silica or silicon nitride membranes affect the DNA translocation speed. In geological processes films of nano- 7 confined water mediate crystallization, dissolution, and material transport along stressed surfaces in contact,[64-66] and mineralization reactions, such as CO 2 sequestration by carbonation, is limited by water transport through passivating silica-rich nano-porous layers formed at the reaction interfaces.[67, 68] Improved understanding of confined water in tight rocks, such as shales or carbonates, is also important to reduce environmental impacts of hydrocarbon production and to ensure permanent CO 2 storage in geological systems. Recent experimental studies indicate that nano-confined water may behave differently than bulk water.[69] It has been suggested that nanoconfinement may provide access to a metastable region, known as the “no man’s land (NML)”, in the phase diagram of bulk water.[70] At ambient pressure, the NML lies between the homogeneous nucleation temperature T h = 235 K and the glass transition temperature T g = 136 K. Below T g , bulk water is frozen either in the low-density amorphous (LDA) or high-density amorphous (HDA) phase depending on the pressure on the system.[71] The presence of LDA and HDA phases has led to an intriguing speculation, namely whether the NML is a mixture of low-density water (LDW) and high-density water (HDW) that are extensions of LDA and HDA, respectively. A couple of theoretical studies—one based on the assumption that there is a first order phase transition between two kinds of water terminating at a critical point[72, 73] and the other being a singularity-free scenario[74, 75]—support the existence of LDW and HDW in the supercooled region. An alternative theoretical scenario is spinodal decomposition,[76, 77] where response functions exhibit anomalous increases close to the spinodal boundary. 8 Recent experimental efforts[78-80] to examine the behavior of LDW and HDW have focused on water confined in nanometer size pores in nanoporous glasses and clays, in glass capillaries and thin films (a few nanometer in thickness) and as droplets suspended in another liquid, e.g., glycerin. Mallamace et al.[81] investigated the behavior of supercooled water confined in hexagonal cylindrical structures in a micelle-templated mesoporous silica (MCM-41-S) matrix with pores of diameter 1.4 nm. Zhang et al.[69] used neutron scattering to investigate liquid-liquid phase transition in deuterated water confined in the nanopores of MCM-41-S silica matrix. It has been suggested that dynamic heterogeneities also exist in supercooled water[78] and they decouple viscosity and self- diffusion coefficient below a certain temperature, causing the breakdown of Stokes- Einstein relation (SER).[82] NMR and quasi-elastic neutron scattering experiments indicate a breakdown of SER in supercooled water confined in a nanotube or a lysozyme protein layer and also in a methanol mixture. 9 Chapter 2. Molecular Dynamics 2.1 Overview Molecular dynamics[83-85] is a numerical technique to predict the spatio-temporal trajectory of a physical system represented by a system particles interacting under a given force field. The system consists of N particles, each representing a single atom, inside an MD box, as shown in Figure 2.1a. The trajectory of the system is the trace of the position of atoms in 3N-dimensional space. Each point on the trajectory is a 3N-element vector, r  N = (x 1 ,y 1 ,z 1 ,x 2 ,y 2 ,z 2 ,...,x N ,y N ,z N ), which represents the position of all the atoms in the system at any given instant (Figure 2.1b). Figure 2.1: (a) Schematic of a sample MD system (b) A typical trajectory of an MD system. 10 The velocity is given by the derivative of the position with respect to time. (2.1) where Δ is the infinitesimal increment in time. Acceleration is defined as the rate of change of velocity with respect to time, given by (2.2) The acceleration of a particle can be estimated from three consecutive positions on its trajectory separated by small time increments: (2.3) From Newton’s second law, we know that the acceleration of a particle is proportional to the force acting on the particle. (2.4) where m = mass of the particle. The same amount of force causes heavier particle to accelerate slower. Typically, MD is used for initial value problems. Given the positions and velocities of a system of N particles, it is used to predict those at later times. Thus, one can predict spatio-temporal evolution of a physical system represented by a system of N particles. 1 Molecular Dynamics Basics Basics of the molecular-dynamics (MD) method 1-3 are described, along with corresponding data structures in program, md.c. Newton’s Second Law of Motion TRAJECTORY, COORDINATE, AND ACCELERATION • Physical system = a set of atomic coordinates: € {  r i = (xi,yi,z i) | xi,yi,z i ∈ℜ,i = 0,...,N−1}, where ℜ is the set of real numbers (in the program, represented by a double precision variable) and we use a vector notation, r r i . (Data strcutures in md.c) int nAtom: N, the number of atoms. NMAX: Maximum number of atoms that can be handled by the program. double r[NMAX][3]: r[i][0], r[i][1], and r[i][2] are the x, y, and z coordinates of the i-th atom, where i = 0, .., N-1. x y z 1 2 N • Trajectory: A mapping from time to a point in the 3-dimensional space, € t∈ℜ  r i (t)∈ℜ 3 . In fact, a trajectory of an N-atom system is regarded as a curve in 3N-dimensional space. A point on the curve is then specified by a 3N-element vector, €  r N = (x 0 ,y 0 ,z 0 ,x 1 ,y 1 ,z 1 ,...,x N−1 ,y N−1 ,z N−1 ). r i (t=0) v i (t=0) r i (t=t 1 ) v i (t=t 1 ) • Velocity: Short-time limit of an average speed (how fast and in which direction the particle is moving), €  v i (t) =  ˙ r i (t) = d  r dt ≡ lim Δ→0  r i (t +Δ)−  r i (t) Δ . double rv[NMAX][3]: rv[i][0], rv[i][1], and rv[i][2] are the x, y, and z components of the velocity vector, r v i , of the i-th atom. • Acceleration: Rate at which a velocity changes (whether the particle is accelerating or decelerating), €  a i (t) =  ˙ ˙ r i (t) = d 2  r dt 2 = d  v i dt ≡ lim Δ→0  v i (t +Δ)−  v i (t) Δ . 1 Molecular Dynamics Basics Basics of the molecular-dynamics (MD) method 1-3 are described, along with corresponding data structures in program, md.c. Newton’s Second Law of Motion TRAJECTORY, COORDINATE, AND ACCELERATION • Physical system = a set of atomic coordinates: € {  r i = (x i ,y i ,z i ) | x i ,y i ,z i ∈ℜ,i = 0,...,N−1}, where ℜ is the set of real numbers (in the program, represented by a double precision variable) and we use a vector notation, r r i . (Data strcutures in md.c) int nAtom: N, the number of atoms. NMAX: Maximum number of atoms that can be handled by the program. double r[NMAX][3]: r[i][0], r[i][1], and r[i][2] are the x, y, and z coordinates of the i-th atom, where i = 0, .., N-1. x y z 1 2 N • Trajectory: A mapping from time to a point in the 3-dimensional space, € t∈ℜ  r i (t)∈ℜ 3 . In fact, a trajectory of an N-atom system is regarded as a curve in 3N-dimensional space. A point on the curve is then specified by a 3N-element vector, €  r N = (x 0 ,y 0 ,z 0 ,x 1 ,y 1 ,z 1 ,...,x N−1 ,y N−1 ,z N−1 ). r i (t=0) v i (t=0) r i (t=t 1 ) v i (t=t 1 ) • Velocity: Short-time limit of an average speed (how fast and in which direction the particle is moving), €  v i (t) =  ˙ r i (t) = d  r dt ≡ lim Δ→0  r i (t +Δ)−  r i (t) Δ . double rv[NMAX][3]: rv[i][0], rv[i][1], and rv[i][2] are the x, y, and z components of the velocity vector, r v i , of the i-th atom. • Acceleration: Rate at which a velocity changes (whether the particle is accelerating or decelerating), €  a i (t) =  ˙ ˙ r i (t) = d 2  r dt 2 = d  v i dt ≡ lim Δ→0  v i (t +Δ)−  v i (t) Δ . 2 double ra[NMAX][3]: ra[i][0], ra[i][1], and ra[i][2] are the x, y, and z components of the acceleration vector, r a i , of the i-th atom. v i (t) v i (t+Δ) v i (t) v i (t+Δ) a i (t)Δ The acceleration of a particle can be estimated from three consecutive positions on its trajectory separated by small time increments: €  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 NEWTON’S EQUATION Newton’s equation states that the acceleration of a particle is proportional to the force acting on the particle, € m  ˙ ˙ r i (t) =  F i (t) , where m is the mass of the particle. For a heavier particle, the same force causes less acceleration (or deceleration). • Initial value problem: Given initial particle positions and velocities, €  r i (0),  v i (0) ( ) i =1,...,N { } , obtain those at later times, €  r i (t),  v i (t) ( ) i =1,...,N;t≥ 0 { } . Note that both positions and velocities must be specified in order to predict future trajectories. Potential Energy We consider forces which are derived from a potential energy, V( r r N ), which is a function of all atomic positions. €  F k =− ∂ ∂  r k V  r N ( ) =− ∂V ∂x k , ∂V ∂y k , ∂V ∂z k       . Here partial derivative is defined as € ∂V ∂x k = lim h→0 V(x 0 ,y 0 ,z 0 ,...,x k +h,y k ,z k ,...x N−1 ,y N−1 ,z N−1 )−V(x 0 ,y 0 ,z 0 ,...,x k ,y k ,z k ,...x N−1 ,y N−1 ,z N−1 ) h , i.e., what is the rate of change in V when we slightly change one coordinate of an atom while keeping all the other coordinates fixed. LENNARD-JONES POTENTIAL To model certain materials such as neon and argon liquids, the Lennard-Jones potential is often used by scientists: 2 double ra[NMAX][3]: ra[i][0], ra[i][1], and ra[i][2] are the x, y, and z components of the acceleration vector, r a i , of the i-th atom. v i (t) v i (t+Δ) v i (t) v i (t+Δ) a i (t)Δ The acceleration of a particle can be estimated from three consecutive positions on its trajectory separated by small time increments: €  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 NEWTON’S EQUATION Newton’s equation states that the acceleration of a particle is proportional to the force acting on the particle, € m  ˙ ˙ r i (t) =  F i (t) , where m is the mass of the particle. For a heavier particle, the same force causes less acceleration (or deceleration). • Initial value problem: Given initial particle positions and velocities, €  r i (0),  v i (0) ( ) i =1,...,N { } , obtain those at later times, €  r i (t),  v i (t) ( ) i =1,...,N;t≥ 0 { } . Note that both positions and velocities must be specified in order to predict future trajectories. Potential Energy We consider forces which are derived from a potential energy, V( r r N ), which is a function of all atomic positions. €  F k =− ∂ ∂  r k V  r N ( ) =− ∂V ∂x k , ∂V ∂y k , ∂V ∂z k       . Here partial derivative is defined as € ∂V ∂x k = lim h→0 V(x 0 ,y 0 ,z 0 ,...,x k +h,y k ,z k ,...x N−1 ,y N−1 ,z N−1 )−V(x 0 ,y 0 ,z 0 ,...,x k ,y k ,z k ,...x N−1 ,y N−1 ,z N−1 ) h , i.e., what is the rate of change in V when we slightly change one coordinate of an atom while keeping all the other coordinates fixed. LENNARD-JONES POTENTIAL To model certain materials such as neon and argon liquids, the Lennard-Jones potential is often used by scientists: 11 2.2 Interatomic Potential To predict the future states, the forces acting on each particle need to be known. Generally, the forces are derived from a potential, which governs the motion of the atoms. The potential is usually a function of the 3N-dimensional position vector. V =V r  N ( ) ; r  N = (x 1 ,y 1 ,z 1 ,x 2 ,y 2 ,z 2 ,...,x N ,y N ,z N ) (2.5) The force on the k th atom is given by (2.6) Here the partial derivative is defined as (2.7) In other words, it is the rate of change in the potential, V, as one coordinate of an atom is incremented by an infinitesimal amount, h, with all the other coordinates fixed. Typically, the potential of the system is the sum of interatomic interactions, u(r ij ) , where r ij  = r i  −r j   is the relative position vector between i th and j th atoms. V = u(r ij ) i<j ∑ (2.8) A typical interatomic potential has multiple terms, each representing physical interactions such as steric repulsion, coulomb interactions, charge-dipole interactions, van der Waals interactions, etc. For a better prediction, the terms in the potential increase to represent as many physical forces as possible, but that also increases computational complexity. 2 double ra[NMAX][3]: ra[i][0], ra[i][1], and ra[i][2] are the x, y, and z components of the acceleration vector, r a i , of the i-th atom. v i (t) v i (t+Δ) v i (t) v i (t+Δ) a i (t)Δ The acceleration of a particle can be estimated from three consecutive positions on its trajectory separated by small time increments: €  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 NEWTON’S EQUATION Newton’s equation states that the acceleration of a particle is proportional to the force acting on the particle, € m  ˙ ˙ r i (t) =  F i (t) , where m is the mass of the particle. For a heavier particle, the same force causes less acceleration (or deceleration). • Initial value problem: Given initial particle positions and velocities, €  r i (0),  v i (0) ( ) i =1,...,N { } , obtain those at later times, €  r i (t),  v i (t) ( ) i =1,...,N;t≥ 0 { } . Note that both positions and velocities must be specified in order to predict future trajectories. Potential Energy We consider forces which are derived from a potential energy, V( r r N ), which is a function of all atomic positions. €  F k =− ∂ ∂  r k V  r N ( ) =− ∂V ∂x k , ∂V ∂y k , ∂V ∂z k       . Here partial derivative is defined as € ∂V ∂x k = lim h→0 V(x 0 ,y 0 ,z 0 ,...,x k +h,y k ,z k ,...x N−1 ,y N−1 ,z N−1 )−V(x 0 ,y 0 ,z 0 ,...,x k ,y k ,z k ,...x N−1 ,y N−1 ,z N−1 ) h , i.e., what is the rate of change in V when we slightly change one coordinate of an atom while keeping all the other coordinates fixed. LENNARD-JONES POTENTIAL To model certain materials such as neon and argon liquids, the Lennard-Jones potential is often used by scientists: 2 double ra[NMAX][3]: ra[i][0], ra[i][1], and ra[i][2] are the x, y, and z components of the acceleration vector, r a i , of the i-th atom. v i (t) v i (t+Δ) v i (t) v i (t+Δ) a i (t)Δ The acceleration of a particle can be estimated from three consecutive positions on its trajectory separated by small time increments: €  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 NEWTON’S EQUATION Newton’s equation states that the acceleration of a particle is proportional to the force acting on the particle, € m  ˙ ˙ r i (t) =  F i (t) , where m is the mass of the particle. For a heavier particle, the same force causes less acceleration (or deceleration). • Initial value problem: Given initial particle positions and velocities, €  r i (0),  v i (0) ( ) i =1,...,N { } , obtain those at later times, €  r i (t),  v i (t) ( ) i =1,...,N;t≥ 0 { } . Note that both positions and velocities must be specified in order to predict future trajectories. Potential Energy We consider forces which are derived from a potential energy, V( r r N ), which is a function of all atomic positions. €  F k =− ∂ ∂  r k V  r N ( ) =− ∂V ∂x k , ∂V ∂y k , ∂V ∂z k       . Here partial derivative is defined as € ∂V ∂x k = lim h→0 V(x 0 ,y 0 ,z 0 ,...,x k +h,y k ,z k ,...x N−1 ,y N−1 ,z N−1 )−V(x 0 ,y 0 ,z 0 ,...,x k ,y k ,z k ,...x N−1 ,y N−1 ,z N−1 ) h , i.e., what is the rate of change in V when we slightly change one coordinate of an atom while keeping all the other coordinates fixed. LENNARD-JONES POTENTIAL To model certain materials such as neon and argon liquids, the Lennard-Jones potential is often used by scientists: 12 The acceleration of k th atom can be derived from the interatomic potential by taking a partial derivative of the potential with respect to the instantaneous position of the atom. (2.9) resuling in (2.10) where δ ik is the Kronecker’s delta expression given by δ ik = 1 (i= k) 0 (i≠ k) " # $ (2.11) 2.3 Discretization In order to solve the equations of motion on a computer, we need to discretize the trajectories. Instead of progressing in continuous time, t, we progress in small timesteps, Δ. Thus the position and velocity vectors are calculated for a sequence of states. (2.12) Thus arises the need for an integration algorithm to predict the next state (r i  (t+Δ),v i  (t+Δ)) from (r i  (t),v i  (t)) . Verlet Discretization uses the position vector for current and previous states to predict the position and velocity for subsequent states: r i  (t+Δ)= 2r i  (t)−2r i  (t−Δ)+a i  (t)Δ 2 +O(Δ 4 ) (2.13) 4 In summary, the normalized Newton’s equation for atoms interacting with the Lennard-Jones potential is given below. (From now on, we always use normalized variables, and omit primes.) € d 2  r i dt 2 =− ∂V ∂  r i =  a i V(  r N ) = u(r ij ) i< j ∑ u(r) = 4 1 r 12 − 1 r 6       The figure in the right shows the normalized Lennard-Jones potential, u(r), as a function of interatomic distance, r. The distance, r 0 , at which the potential takes its minimum value, u(r 0 ), is obtained as follows. € du dr r=r 0 =0 ⇒ r 0 =2 1/6 ≈1.12 u(r 0 ) = 4 1 2 12/6 − 1 2 6/6       =−1 Figure: Lennard-Jones potential. ANALYTIC FORMULA FOR FORCES €  a k =− ∂ ∂  r k u(r ij ) i< j ∑ =− ∂r ij ∂  r k i< j ∑ du dr ij Let’s evaluate the two factors in the above expression: 1. € ∂r ij ∂  r k = ∂ ∂x k , ∂ ∂y k , ∂ ∂z k       x i − x j ( ) 2 + y i − y j ( ) 2 + z i − z j ( ) 2 = 2(x i − x j ),2(y i − y j ),2(z i − z j ) ( ) 2 x i − x j ( ) 2 + y i − y j ( ) 2 + z i − z j ( ) 2 δ ik −δ jk ( ) =  r ij r ij δ ik −δ jk ( )  d dx f (x) [ ] 1/2 = 1 2 f (x) [ ] −1/2 df dx       2. € du dr =4 − 12 r 13 + 6 r 7       =− 48 r 1 r 12 − 1 2r 6        d dr r −n =−nr −n−1       In the first result, € δ ik = 1 (i = k) 0 (i≠ k)    is called Kronecker’s delta expression. 5 Substituting these two results back into the expression for r a k , we obtain €  a k =  r ij − 1 r du dr       r=r ij i< j ∑ δ ik −δ jk ( ) where € − 1 r du dr = 48 r 2 1 r 12 − 1 2r 6       Summary of Molecular-Dynamics Equation System Given initial atomic positions and velocities, €  r i (0),  v i (0) ( ) i =1,...,N { } , obtain those at later times, €  r i (t),  v i (t) ( ) i =1,...,N;t≥ 0 { } , by integrating the following differential equation: €  ˙ ˙ r k (t) =  a k (t) =  r ij (t) − 1 r du dr       r=r ij (t) i< j ∑ δ ik −δ jk ( ) where € − 1 r du dr = 48 r 2 1 r 12 − 1 2r 6       The sum over i and j to evaluate the acceleration is implemented in a program as follows: for i = 1 to N,  a i = 0 , for i = 1 to N−1 for j = i+1 to N compute €  α =  r ij − 1 r du dr       r=  r ij €  a i + =  α €  a j − =  α Discretization We need to discretize the trajectories in order to solve the problem on a computer: Instead of considering €  r i (t),  v i (t) ( ) for t≥ 0 for continuous time, we consider a sequence of states €  r i (0),  v i (0) ( )   r i (Δ),  v i (Δ) ( )   r i (2Δ),  v i (2Δ) ( )  r i (0) r i (2Δ) r i (Δ) The question is: How to predict the next state, €  r i (t +Δ),  v i (t +Δ) ( ) , from the current state, €  r i (t),  v i (t) ( ) . double DeltaT: Δ in md.c. 5 Substituting these two results back into the expression for r a k , we obtain €  a k =  r ij − 1 r du dr       r=r ij i< j ∑ δ ik −δ jk ( ) where € − 1 r du dr = 48 r 2 1 r 12 − 1 2r 6       Summary of Molecular-Dynamics Equation System Given initial atomic positions and velocities, €  r i (0),  v i (0) ( ) i =1,...,N { } , obtain those at later times, €  r i (t),  v i (t) ( ) i =1,...,N;t≥ 0 { } , by integrating the following differential equation: €  ˙ ˙ r k (t) =  a k (t) =  r ij (t) − 1 r du dr       r=r ij (t) i< j ∑ δ ik −δ jk ( ) where € − 1 r du dr = 48 r 2 1 r 12 − 1 2r 6       The sum over i and j to evaluate the acceleration is implemented in a program as follows: for i = 1 to N,  a i = 0 , for i = 1 to N−1 for j = i+1 to N compute €  α =  r ij − 1 r du dr       r=  r ij €  a i + =  α €  a j − =  α Discretization We need to discretize the trajectories in order to solve the problem on a computer: Instead of considering €  r i (t),  v i (t) ( ) for t≥ 0 for continuous time, we consider a sequence of states €  r i (0),  v i (0) ( )   r i (Δ),  v i (Δ) ( )   r i (2Δ),  v i (2Δ) ( )  r i (0) r i (2Δ) r i (Δ) The question is: How to predict the next state, €  r i (t +Δ),  v i (t +Δ) ( ) , from the current state, €  r i (t),  v i (t) ( ) . double DeltaT: Δ in md.c. 13 v i  (t)= r i  (t+Δ)−r i  (t−Δ) 2Δ +O(Δ 2 ) (2.14) The Verlet algorithm, which uses equations 2.13 and 2.14, is shown below: Given r i  (t) and r i  (t−Δ) , 1. Compute a i  (t) from r i  (t) { } 2. r i  (t+Δ)= 2r i  (t)−2r i  (t−Δ)+a i  (t)Δ 2 3. v i  (t)= (r i  (t+Δ)−r i  (t−Δ))/(2Δ) Verlet algorithm requires the position at time t+Δ to calculate the velocity at time t. Therefore we use a variation of the Verlet algorithm called the Velocity-Verlet algorithm[86] which calculates the position and velocity at time t+Δ using those at time t. r i  (t+Δ)= r i  (t)+v i  (t)Δ+ 1 2 a i  (t)Δ 2 +O(Δ 4 ) (2.15) v i  (t+Δ)=v i  (t)+ a i  (t)−a i  (t+Δ) 2 Δ+O(Δ 3 ) (2.16) The Velocity-Verlet algorithm is as follows: Given (r i  (t),v i  (t)) , 1. Compute a i  (t) from r i  (t) { } 2. v i  (t+ Δ 2 )= v i  (t)+ Δ 2 a i  (t) 3. r i  (t+Δ)= r i  (t)+v i  (t+ Δ 2 )Δ 4. Compute a i  (t+Δ) from r i  (t+Δ) { } 14 5. v i  (t+Δ)= v i  (t+ Δ 2 )+ Δ 2 a i  (t+Δ) Figure 2.2: Schematic of (a) Periodic Boundary Conditions and (b) Minimum Image Convention To simulate bulk material using small number of atoms, we subject the MD box to periodic boundary conditions. The entire MD box is replicated in three dimensions and the atoms close to the edge of the box interact with the atoms in the replicated image. Figure 2.2a shows the schematic of replicated MD boxes for periodic boundary conditions. Displacement of an atom within the MD box is reflected in all the replicated images. Since all the images are just shifted copies of the original MD box, we need to store only the coordinates of the atoms in the original MD box. If an atom moves out of the MD, it moves into the replicated image of the MD box. This is equivalent to the atom entering the original box from the opposite side. At all times during the simulation, the following conditions are maintained: 0≤ x i ≤ L x ,0≤ y i ≤ L y ,0≤ z i ≤ L z ; L x ,L y ,L z = Length of MD box in x,y and z directions. Periodic boundary conditions are applied using the following equations: 15 x i = x i −L x , if x i > L x (2.17) x i = x i +L x , if x i <0 (2.18) Theoretically, an atom interacts with all the replicated images but due to computational limitations, the dimensions of the MD box are kept such that the contribution from the atoms at a distance, r > L x /2 can be neglected. Thus, an atom interacts with only the nearest image of another atom. Imagine a box of size L x x L y x L z centered at i th atom (Figure 2.2b), which interacts with all the atoms, j, in this box. Thus we have − 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.19) This formulation, called Minimum Image Convention, is implemented as x ij = x ij + L x 2 , if x ij ≤− L x 2 (2.20) x ij = x ij − L x 2 , if x ij ≥ L x 2 (2.21) Minimum Image Convention limits the cut-off distance for interatomic interactions at half of the length of the MD box. For small boxes, this truncation results in significant discontinuity in the potential function rendering the force indifferentiable. 16 Figure 2.3: Lennard-Jones Potential (blue) and its shifted form (red) To avoid this situation, we cut the potential at a pre-determined cut-off distance, , r c < L x /2. The potential is modified such that both u(r) and du/dr are continuous at r c . (2.22) Figure 2.3 shows a typical Lennard-Jones Potential in blue and its shifted version in red. Using a shifted potential allows us to limit the calculations by using a cut-off without incorporating boundary errors at the cut-off distance. 2.4 Linked List Cell Decomposition To calculate force on each atom using the above expressions, the simplest algorithm is as follows: 11 Or we can imagine a box of size L x × L y × L z centered at atom i, and this atom interacts only with other atoms, j, in this imaginary box. 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          during the computation of forces, where €  r ij =  r i −  r j = (x ij ,y ij ,z ij ). In the program, this is achieved by x ij ← x ij −SignR L x 2 ,x ij + L x 2       −SignR L x 2 ,x i − L x 2       y ij ← y ij −SignR L y 2 ,y ij + L y 2       −SignR L x 2 ,y ij − L x 2       z ij ← z ij −SignR L z 2 ,z ij + L z 2       −SignR L z 2 ,z ij − L z 2                Shifted Potential With the minimal image convention, interatomic interaction is cut-off at half the box length. If the box is not large enough, this truncation causes a significant discontinuity in the potential function and our derivation of forces breaks down (we cannot differentiate the potential). u(r) r L x /2 Potential energy jump To remedy this, we introduce a cut-off length r c < L x /2, etc., and modify the potential such that both u(r) and du/dr are continuous at r c . ′ u (r) = u(r)− u(r c )− (r− r c ) du dr r =r c r < r c 0 r≥ r c      For the Lennard-Jones potential, r c = 2.5σ is conventionally used. This causes some shift in the potential minimum but little shift in the minimum position. 17 As evident, this algorithm is of the order O(N 2 ), which is not scalable to large systems. In order to perform large-scale simulations, we use the method of linked list cell decomposition which is of the order O(N). We divide the MD box into a 3 dimensional grid (Figure 2.4a) with the length of each grid, lc, such that lc is slightly greater than the cut-off distance, r c . This is achieved by: mc = L box / r c ; L box = Length of the box, mc = number of grids lc = L box /mc We then loop over all the atoms and each atom is placed in its respective grid. We use a three dimensional array to store the current header of the grid and a one dimensional array to store the connectivity information. Both these arrays are initialized as NULL. As an atom is placed in a grid, the index of the atom becomes the header of the grid. When a new atom is placed in the same grid, the index of the new atom becomes the new header and it points to the previous atom in the grid, as shown in Figure 2.4b. Table 2.1 shows a sample connectivity list. To traverse through all the connected atoms, we start with the header of a grid and look for the value stored at that index in the connectivity list. This process is repeated till we reach NULL value, which means that there are no more connected atoms. This method of decomposition ensures that all the atoms within the cut- off distance lie within the nearest neighboring grids from the resident grid. Thus the 5 Substituting these two results back into the expression for r a k , we obtain €  a k =  r ij − 1 r du dr       r=r ij i< j ∑ δ ik −δ jk ( ) where € − 1 r du dr = 48 r 2 1 r 12 − 1 2r 6       Summary of Molecular-Dynamics Equation System Given initial atomic positions and velocities, €  r i (0),  v i (0) ( ) i =1,...,N { } , obtain those at later times, €  r i (t),  v i (t) ( ) i =1,...,N;t≥ 0 { } , by integrating the following differential equation: €  ˙ ˙ r k (t) =  a k (t) =  r ij (t) − 1 r du dr       r=r ij (t) i< j ∑ δ ik −δ jk ( ) where € − 1 r du dr = 48 r 2 1 r 12 − 1 2r 6       The sum over i and j to evaluate the acceleration is implemented in a program as follows: for i = 1 to N,  a i = 0 , for i = 1 to N−1 for j = i+1 to N compute €  α =  r ij − 1 r du dr       r=  r ij €  a i + =  α €  a j − =  α Discretization We need to discretize the trajectories in order to solve the problem on a computer: Instead of considering €  r i (t),  v i (t) ( ) for t≥ 0 for continuous time, we consider a sequence of states €  r i (0),  v i (0) ( )   r i (Δ),  v i (Δ) ( )   r i (2Δ),  v i (2Δ) ( )  r i (0) r i (2Δ) r i (Δ) The question is: How to predict the next state, €  r i (t +Δ),  v i (t +Δ) ( ) , from the current state, €  r i (t),  v i (t) ( ) . double DeltaT: Δ in md.c. 18 search space is reduced to 27 grids in 3 dimensional space, thereby making the algorithm of the order O(N). Having a search algorithm of O(N) enables us to perform very large- scale MD simulations involving upto billions of atoms, which would not be possible with O(N 2 ) algorithm. index 94 73 … 56 … 31 25 … 14 8 value 73 56 … 31 … 25 14 … 8 NULL Table 2.1 A sample linked list array representing the connectivity shown in Figure 2.2b. Figure 2.4: (a) Schematic showing linked list cell decomposition. (b) The connected components in a sample linked list. 19 2.5 Parallel Molecular Dynamics Due to memory and speed limitations, single processor results in a limit to the maximum size of the system that can be simulated. In order to perform large-scale MD simulations, we use multiple, inter-connected processors, which can perform parallel tasks. We partition the physical system to be simulated into sub-systems of equal volume, and assign each sub-system to a single processor.[87] This distribution is called Spatial Distribution. Processors in a parallel computer are logically arranged according to the topology of the physical system. We use a simple 3D mesh spatial decomposition. The subsystems are arranged in a 3D array of dimensions P x x P y x P z , where P x , P y , and P z are the number of processors in x, y and z directions, respectively. The processors follow the same logical arrangement. Figure 2.5 shows the schematic arrangement of processors in 3D. Figure 2.5: Schematic showing spatial distribution of system among multiple processors 20 Each processor is given a unique ID ranging between 0 and P -1, where P = P x x P y x P z is the total number of processors being used. To identify the processors in 3D, we also assign a vector ID, p, to each processor. The vector ID has three elements p x , p y , and p z , where p x = 0, …, P x -1; p y = 0, …, P y -1; and p z = 0, …, P z -1. The vector IDs are assigned as p x = p/(P y x P z ); p y = (p/P z ) mod P y ; and p z = p mod P z . Thus we have p = p x x P y x P z + p y x P z . Every subsystem has 26 neighbor subsystems. To compute interatomic interaction in a current (resident) processor with cut-off length, r c , the positions of the atoms near the boundary in its neighbor subsytems are required. The coordinates of the atoms located within r c from the subsystem boundary are copied to the resident processor. Figure 2.6 shows a schematic representation of this process. The resident processor is shown in brown, while the neighbors are shown in yellow. The green region, which is of thickness r c , represents the copied region. Figure 2.6: Schematic showing the copied region among neighboring subsystems 21 We use Message Passing Interface (MPI) to communicate between processors in a parallel framework through multiple send and receive commands. During sharing the copied information, there is a possibility of a deadlock if there arises a circular send and receive. In other words, due to finite size in the receiver’s communication system, a sender waits till the receiver’s buffer is clear. But in the situation shown in Figure 2.7a, none of the processors can start receiving until its send operation is completed. To avoid this situation, we perform the send and receive operations in two phases, as shown in Figure 2.7b. In phase 1, the even-tagged processors send information to the odd-tagged processors and the odd-tagged processors receive information from the even-tagged processors. In phase 2, we reverse the order and the odd-tagged processors send information to the even-tagged processors and the even-tagged processors receive information from the odd-tagged processors. In this way, after phase 2, all the processors have the necessary information, while avoiding deadlock. Figure 2.7: (a) Schematic showing communication deadlock. (b) Two-phase message passing scheme to avoid deadlock. 22 It is important to estimate the efficiency of a parallel program in order to properly plan and execute a large-scale simulation. Let us consider a system of size, W, running on P parallel processors. Let T(W,P) be the execution time of the this parallel program. Speed of the program is defined as S(W,P) = W/ T(W,P). Speedup is defined as the speed of P processors divided by that of 1 processor, giving S P (W,P) = S(W P ,P)/ S(W 1 ,1). We define parallel efficiency as E P = S P /P, since the ideal speedup on P processors is expected to be P. For constant problem-size scaling, the system size, W, is fixed and the workload per processor decreases with increasing number of processors. The speedup in this case is (2.23) and the parallel efficiency is given by (2.24) In case of isogranular scaling, the workload per processor, w, is fixed and the total system size increases with increase in the number of processors. The speedup is given by (2.25) and the isogranular parallel efficiency is given by (2.26) 11 Scalability Metrics for Parallel Molecular Dynamics Parallel Efficiency We define the efficiency of a parallel program running on P processors to solve a problem of size W. Let T(W, P) be the execution time of this parallel program. Speed of the program is then S(W, P) = W/T(W, P). Speedup, S P , on P processors is the speed of P processors divided by that of 1 processor, i.e., S P = S(W P , P)/S(W 1 , 1). To unambiguously define the speedup, we need to specify how the problem size, W P , scales as a function of the number of processors, P (which will be discussed in the next few paragraphs). The ideal speedup on P processors is expected to be P, and therefore we define the parallel efficiency, E P = S P /P. CONSTANT PROBLEM-SIZE SCALING In the constant problem-size scaling, the problem size W P = W is fixed constant. Therefore, the constant problem-size speedup is € S P = S(W,P) S(W,1) = W /T(W,P) W /T(W,1) = T(W,1) T(W,P) , and the parallel efficiency is € E P = S P P = T(W,1) P •T(W,P) . Amdahl’s law: Consider a case, in which a fraction, f (∈ [0, 1]), of the work W is inherently sequential and cannot be parallelized, but the rest, 1 − f, can be divided and executed in parallel on P processors. Then, T(W, P) = f•T(W, 1) + (1 − f)T(W, 1)/P. Therefore, the speedup is € S P = T(W,1) T(W,P) = 1 f + (1− f )/P . In the limit of large number of processors, P → ∞, the asymptotic speedup is S P → 1/f. The constant problem-size speedup is thus limited by the fraction of the sequential bottleneck of the parallel program. For example, if 1% of the work is sequential, the maximum achievable speedup is 0.01/1 = 100, and it does not make sense to use more than ~100 processors to run this program. ISOGRANULAR SCALING In the isogranular scaling, the problem size W P scales linearly with the number of processors: W P = P•w, where the granularity (or the work per processor), w, is constant. Therefore, the isogranular speedup is € S P = S(P •w,P) S(w,1) = P •w/T(P •w,P) w/T(w,1) = P •T(w,1) T(P •w,P) , and the corresponding isogranular parallel efficiency is € E P = S P P = T(w,1) T(P •w,P) . 11 Scalability Metrics for Parallel Molecular Dynamics Parallel Efficiency We define the efficiency of a parallel program running on P processors to solve a problem of size W. Let T(W, P) be the execution time of this parallel program. Speed of the program is then S(W, P) = W/T(W, P). Speedup, S P , on P processors is the speed of P processors divided by that of 1 processor, i.e., S P = S(W P , P)/S(W 1 , 1). To unambiguously define the speedup, we need to specify how the problem size, W P , scales as a function of the number of processors, P (which will be discussed in the next few paragraphs). The ideal speedup on P processors is expected to be P, and therefore we define the parallel efficiency, E P = S P /P. CONSTANT PROBLEM-SIZE SCALING In the constant problem-size scaling, the problem size W P = W is fixed constant. Therefore, the constant problem-size speedup is € S P = S(W,P) S(W,1) = W /T(W,P) W /T(W,1) = T(W,1) T(W,P) , and the parallel efficiency is € E P = S P P = T(W,1) P •T(W,P) . Amdahl’s law: Consider a case, in which a fraction, f (∈ [0, 1]), of the work W is inherently sequential and cannot be parallelized, but the rest, 1 − f, can be divided and executed in parallel on P processors. Then, T(W, P) = f•T(W, 1) + (1 − f)T(W, 1)/P. Therefore, the speedup is € S P = T(W,1) T(W,P) = 1 f + (1− f )/P . In the limit of large number of processors, P → ∞, the asymptotic speedup is S P → 1/f. The constant problem-size speedup is thus limited by the fraction of the sequential bottleneck of the parallel program. For example, if 1% of the work is sequential, the maximum achievable speedup is 0.01/1 = 100, and it does not make sense to use more than ~100 processors to run this program. ISOGRANULAR SCALING In the isogranular scaling, the problem size W P scales linearly with the number of processors: W P = P•w, where the granularity (or the work per processor), w, is constant. Therefore, the isogranular speedup is € S P = S(P •w,P) S(w,1) = P •w/T(P •w,P) w/T(w,1) = P •T(w,1) T(P •w,P) , and the corresponding isogranular parallel efficiency is € E P = S P P = T(w,1) T(P •w,P) . 11 Scalability Metrics for Parallel Molecular Dynamics Parallel Efficiency We define the efficiency of a parallel program running on P processors to solve a problem of size W. Let T(W, P) be the execution time of this parallel program. Speed of the program is then S(W, P) = W/T(W, P). Speedup, S P , on P processors is the speed of P processors divided by that of 1 processor, i.e., S P = S(W P , P)/S(W 1 , 1). To unambiguously define the speedup, we need to specify how the problem size, W P , scales as a function of the number of processors, P (which will be discussed in the next few paragraphs). The ideal speedup on P processors is expected to be P, and therefore we define the parallel efficiency, E P = S P /P. CONSTANT PROBLEM-SIZE SCALING In the constant problem-size scaling, the problem size W P = W is fixed constant. Therefore, the constant problem-size speedup is € S P = S(W,P) S(W,1) = W /T(W,P) W /T(W,1) = T(W,1) T(W,P) , and the parallel efficiency is € E P = S P P = T(W,1) P •T(W,P) . Amdahl’s law: Consider a case, in which a fraction, f (∈ [0, 1]), of the work W is inherently sequential and cannot be parallelized, but the rest, 1 − f, can be divided and executed in parallel on P processors. Then, T(W, P) = f•T(W, 1) + (1 − f)T(W, 1)/P. Therefore, the speedup is € S P = T(W,1) T(W,P) = 1 f + (1− f )/P . In the limit of large number of processors, P → ∞, the asymptotic speedup is S P → 1/f. The constant problem-size speedup is thus limited by the fraction of the sequential bottleneck of the parallel program. For example, if 1% of the work is sequential, the maximum achievable speedup is 0.01/1 = 100, and it does not make sense to use more than ~100 processors to run this program. ISOGRANULAR SCALING In the isogranular scaling, the problem size W P scales linearly with the number of processors: W P = P•w, where the granularity (or the work per processor), w, is constant. Therefore, the isogranular speedup is € S P = S(P •w,P) S(w,1) = P •w/T(P •w,P) w/T(w,1) = P •T(w,1) T(P •w,P) , and the corresponding isogranular parallel efficiency is € E P = S P P = T(w,1) T(P •w,P) . 11 Scalability Metrics for Parallel Molecular Dynamics Parallel Efficiency We define the efficiency of a parallel program running on P processors to solve a problem of size W. Let T(W, P) be the execution time of this parallel program. Speed of the program is then S(W, P) = W/T(W, P). Speedup, S P , on P processors is the speed of P processors divided by that of 1 processor, i.e., S P = S(W P , P)/S(W 1 , 1). To unambiguously define the speedup, we need to specify how the problem size, W P , scales as a function of the number of processors, P (which will be discussed in the next few paragraphs). The ideal speedup on P processors is expected to be P, and therefore we define the parallel efficiency, E P = S P /P. CONSTANT PROBLEM-SIZE SCALING In the constant problem-size scaling, the problem size W P = W is fixed constant. Therefore, the constant problem-size speedup is € S P = S(W,P) S(W,1) = W /T(W,P) W /T(W,1) = T(W,1) T(W,P) , and the parallel efficiency is € E P = S P P = T(W,1) P •T(W,P) . Amdahl’s law: Consider a case, in which a fraction, f (∈ [0, 1]), of the work W is inherently sequential and cannot be parallelized, but the rest, 1 − f, can be divided and executed in parallel on P processors. Then, T(W, P) = f•T(W, 1) + (1 − f)T(W, 1)/P. Therefore, the speedup is € S P = T(W,1) T(W,P) = 1 f + (1− f )/P . In the limit of large number of processors, P → ∞, the asymptotic speedup is S P → 1/f. The constant problem-size speedup is thus limited by the fraction of the sequential bottleneck of the parallel program. For example, if 1% of the work is sequential, the maximum achievable speedup is 0.01/1 = 100, and it does not make sense to use more than ~100 processors to run this program. ISOGRANULAR SCALING In the isogranular scaling, the problem size W P scales linearly with the number of processors: W P = P•w, where the granularity (or the work per processor), w, is constant. Therefore, the isogranular speedup is € S P = S(P •w,P) S(w,1) = P •w/T(P •w,P) w/T(w,1) = P •T(w,1) T(P •w,P) , and the corresponding isogranular parallel efficiency is € E P = S P P = T(w,1) T(P •w,P) . 23 2.6 Molecular Dynamics on BlueGene/P The IBM BlueGene/P supercomputer, called Intrepid, is located at the Argonne Leadership Computing Facility, Argonne National Laboratory in Chicago, IL. It consists of 40 racks with 1024 nodes per rack. Each node has 850 MHz quad-core processor and 2 GB RAM. There are 640 I/O nodes with 7.6 PB of storage. With a total of 163840 cores, 80 TB of RAM and a peak performance of of 557 Teraflops, Intrepid is a highly scalable parallel network suitable for a wide variety of science and engineering problems. While running an MD simulation on Intrepid, there are two important considerations: 1) Aggregated I/O To perform a wide range of post-simulation analyses, we need to store the system configuration at frequent intervals. It is not feasible to open 163840 files simultaneously to write and store intermediate MD data. We circumvent this problem by using a scheme of aggregated I/O. In this scheme, we write the MD output from a pre-specified number of processors into a single file. For example, FILE0 holds the MD output from processor ID 0 to 255; FILE1 holds the same for processor ID 256 to 511; and so on. Thus the number of simultaneously open files is reduced from 163840 to 640. 2) Binary format conversion We store the MD data in machine-readable binary format to save space. While HPCC file-system stores binary in ‘little-endian’ formant, Intrepid, does the same in ‘big- endian’ format. In order to read/write binary files for Intrepid, we use the intrinsic ‘convert=big_endian’ command from Intel’s Fortran compiler. By using this command in the read and write statements, we can use Intrepid binary files on HPCC for analyses. 24 Chapter 3. Oxidation of Aluminum Nanoparticle Aggregate 3.1 Background Aluminum nanoparticles are widely used in a wide range of energetic applications due to their high energy to volume ratio. Several researchers are investigating the behavior of aluminum nanoparticles in different conditions using experimental as well as numerical methods. Ermoline et al.[10] have observed that exothermic reactions leading to ignition begin at relatively low temperatures for nanocomposite materials prepared by Arrested Reactive Milling (ARM). The rate of such reactions is limited by diffusion through the oxide layer. Exposure of fresh metal surface to an oxidizing gas results in rapid oxidation till the thickness the oxide layers reached a temperature dependent limit. An increase in temperature results in further oxidation till the oxide thickness reaches the corresponding limit. Jayaraman et al.[88] have used exploding wire technique to characterize and study its combustion behavior of nanoaluminum. The nanoparticles this produced were spherical in shape and it was found that the addition of nanoaluminum doubles the burning rate as compared to micron-sized aluminum particles, which is attributed to near- complete combustion of nanoaluminum close to the burning surface. In their study of flame propagation rates in aluminum nanoparticles treated with perfluoroalkyl carboxylic acid, Dikici et al.[89] have found that the rate decreases by a factor of 22-95 in the absence of alumina passivation shell as compared to the nanoparticles with an oxide shell. The reaction efficiency increases as the particle size decreases due to the increase in corresponding specific surface area, resulting in higher flame speed. According to the authors, this might be the reason that aluminum powder mixtures with 20-50 nm sized 25 nanoparticles can react many times faster than mixtures with micron-sized particles. Huang et al.[11] performed experimental and analytical studies on the combustion of bimodal nano/micron-sized mixture of aluminum particles in air. They observed that low percentage of nanoparticles in fuel result in flame with wider regime and separated spiral structure, while the flame configurations overlapped for higher nanoparticle percentages, leading to the conclusion that the size of the particles influences the combustion mechanism. De Luca et al.[3] used Brunauer-Emmet-Teller method, electron microscopy, X-ray diffraction and X-ray photoelectron spectroscopy to study combustion mechanism of rocket propellants containing nanoaluminum. They found that the burning rates became steadier and more luminous as the mass fraction of nanoaluminum was increased. The presence of nanoaluminum reduced agglomeration and aggregation of metal particles. The authors speculate that due to the high reactivity of nanoaluminum, substantial combustion occurs within a narrow and very luminous layer close to the burning surface. Kwon et al.[12] measured temperature and radiation during combustion of superfine (~0.1 mm) aluminum powders produced by an electrically heated wire in static air at 1 atm. Using scanning electron microscopy, X-ray diffraction and chemical analysis, they identified two-stage self-propagating regime after ignition by local heating. The fist stage resulted in 70% by mass unreacted aluminum and amorphous oxides with traces of AlN. The AlN content exceeded 50% and the residual Al content decreased to ~10% after the second stage. At low temperatures, Al2O3, Al2O and AlN formed in addition to heat liberation whereas only Al2O3 was found to be exothermic at high temperatures. Yang and collaborators[90] studied the propagation of shock-induced chemical reactions over nanometer distances in energetic materials consisting of Al 26 nanoparticles in the range 30-110 nm in polymer oxidizers in nitrocellulose and Teflon. They used picosecond laser flash heating to vaporize the Al particles, which react with the surrounding oxidizer and generates a spherical shock wave with a rapidly dropping pressure that decomposes the oxidizer. The reaction volume is spherical with diameter in the range 50-1500 nm and is proportional to the laser fluence. To study the oxidation of aluminum powders at elevated temperatures, Trunov et al.[91] heated aluminum powders of different particle sizes and surface morphologies in oxygen upto 1500 °C at various heating rates. They analyzed partially oxidized samples from selected intermediate temperatures by X-ray diffraction. The reaction kinetics could be defined into several stages in the temperature range 300-1500 °C, which also included the structural transformation of alumina. In their experimental study of combustion of micron sized Al particles, Dreizin and collaborators[92] found that the particle temperature decreases continuously and the combustion terminates after the temperature reaches the melting point of Al 2 O 3 . It was also found that the combustion time of Al particles reduces under the influence of electric field. Bockmon et al.[7] have investigated the combustion speed of Al nanoparticles in the range 44-121 nm with MoO3 oxidizer. They observed that the combustion speed increased from 600 m/s to 1000 m/s when the size of the nanoparticle was decreased from 121 nm to 80 nm. The speed remained unchanged upon further reduction in size to 44 nm. They speculate that a larger fraction of oxide in smaller nanoparticles may act as a heat sink in the reaction process. Wang and collaborators[93] have studied the effect of oxide layer thickness on the combustion behavior of 60 nm Al nanoparticles embedded in nitrocellulose at different heating rates. It was found that the reaction propagation distance is larger for thinner oxide layer (2.5 nm) at heating rates < 27 1011 K/s and smaller for thicker oxide layer (6 nm). In contrast, the reaction propagation distance was larger in case of thicker oxide layer for heating rates > 1015 K/s. Similar effect was observed with nanoparticles in the range 30-215 nm with other oxidizing polymers. Zachariah and collaborators[94] studied the effect of the size of aluminum nanoparticle on the oxidation rate using single particle mass spectrometer. The reactivity of aluminum was found to increase with decrease in particle size. The smallest nanoparticle (~19 nm) was the most reactive with 68% oxidation at 900 °C, whereas nanoparticles greater than 50 nm underwent only 4% oxidation at 1100 °C. It was also found that as the particle size decreases, activation energy decreases and rate constant, limited by size-dependent diffusion, increases. Kwon et al.[70] used infrared and X-ray photoelectron spectroscopy, differential scanning calorimetry, high-resolution transmission electron microscopy as well as first principle calculations to study the role of interfacial chemistry on the reactivity of Al/CuO nanomaterial. They found that the alumina layers deposited via atomic layer deposition formed a stronger barrier against interlayer diffusion, even for ultrathin layers (~0.5 nm), as compared to the interface layer formed via physical deposition. They conclude that the kinetics of Al diffusion is controlled by the nature of the interface, instead of the thickness. Asay et al.[20] have used high-speed photography and pressure transducers to identify the mechanism of reaction propagation in a mixture of 44 nm Al particles and 15 nm MoO 3 particles. They have proposed convection as the most likely mechanism of combustion front propagation but the results are not conclusive. Bazyn et al.[95] compared ignition and combustion behavior of 10 mm, 2.7 mm and nano-Al particles, with Al:MoO3 ratios of 4:1, 8:1 and 16:1 respectively. Smaller Al particles were observed to have lower ignition temperature. 28 Spectroscopy results showed strong AlO features in all the samples. The ignition benefit increases at the expense of specific energy content with increasing amount of MoO 3 . Yan et al.[96] have investigated the thermal behavior of 30 – 85 nm aluminum nanoparticles reacting with propylene oxide. X-ray diffraction analysis of the products revealed different phases of alumina in different temperature regions. The emission temperature was determined to be 2705±150 K from the strength for Al-O spectrum. Pivkina et al.[97] have used electron and atomic force microscopy, X-ray powder and thermal analyses to study structural and combustion behavior of nanoaluminum (20-50 nm) mixed with ultrafine ammonium perchlorate powder (35 nm). They found that nanoaluminum sample gained 3.5 times more than that by micron-size aluminum. The activation energy of nanoaluminum was lower than that of micron-size aluminum. The burning rates were found to increase by an order of magnitude when micron-size aluminum was replaced by nanoaluminum in stoichiometric compositions. Morgan and collaborators[98] studied oxide-passivated micro-scale (100 mm) and nano-scale (30-250 nm) aluminum nanoparticles using pyrolysis combustion flow calorimetry. The temperature threshold for reactions was found to decrease from 600-650 °C for micron- size particles to 500-550 °C for 10 nm nanoparticles. Oxygen consumption was found to occur in two distinct steps in nanoparticles while it was a one-step process in case of micron-size particles. Brandstadt et al.[99] have investigated pre-ignition characteristics of nano- (<500 nm) and micro-scale (1-10 mm) Al particles in CO2 oxidizing system using simultaneous thermogravimetric analysis, differential scanning calorimetry, X-ray diffraction spectrometry and energy dispersion X-ray spectrometry. Nanoaluminum underwent highly exothermic reaction in the temperature range 500 – 775 °C. After the 29 temperature reached 1000 °C, the particles attained a hollow structure with Al2O3 exterior shells. The authors speculate that this structure is the result of outwards Al diffusion owing to small diffusion distances in nanoaluminum. Moore et al.[100] carried out experiments on ignition sensitivity and combustion velocity on a thermite composed of micron-sie Al particles (4 and 20 mm) and 80 nm Al nanparticle mixtures with MoO 3 . The ignition delay times for mixtures with 20% nanoaluminum were two orders of magnitude lower than those with micro-size particles. The combustion speed also increased when micron-size Al particles were replaced by nanoaluminum. In their study of the size effect on combustion on aluminum particle dust in a laminar air flow, Huang and collaborators[101] found that the flame speed increases and the combustion mechanism changes from diffusion-controlled to kinetic-controlled as the particle size decreases from micro to nanoscale. It was found that the flame speed followed d -0.92 -law for micron-size particles and d -0.52 -law for nano-sized particles. Gan and Qiao[102] conducted experiments on fuel droplets (n-decane and ethanol) with nano- (~ 80 nm) and micron- (~ 5 mm and ~ 25 mm) sized aluminum particles. Generally, nano-suspension is better than micron suspension, and ethanol based fuel achieved better suspensions than n- decane. Pre-heating and ignition, classical combustion, micro-explosion, surfactant flame and aluminum droplet flame were identified as the five stages of combustion. The formation of oxide shell on the Al nanoparticle surface resulted in longer and less complete combustion of large agglomerates of nano-suspensions. In contrast, micro- suspensions lead to much stronger explosion due to the rupture of the oxide shell. The authors conclude that nano-suspensions are dominated by random Brownian motion, while droplet surface regression, bubble formation, and internal circulation dominate the 30 micron suspensions. Firmansyah et al.[103] used High-Temperature X-Ray Diffraction analysis to examine the microstructural behavior of Al2O3 shell and Al core of 100 nm Al nanoparticles at various temperatures. The Al core was found to expand as the temperature increased from 25 to 800 °C. The amorphous alumina shell expanded from ~4.3 to ~10 nm and underwent phase transformations during heating. The inhomogeneous nature of the transformations resulted in shell fracture with sharp edges. The authors speculate that the regions in the shell with softer amorphous phase or amorphous-gamma interface provide buffer for the aluminum core to expand freely. Yan et al.[104] used X-ray diffraction, transmission electron microscopy, scanning electron microscopy and X-ray photoelectron spectroscopy to compare morphology, particle size and agglomeration process of micro- (4.5 µ) and nano-scale (75 nm) aluminum particles. The initial stage of sintering in nano-size powders was found to be dominated by grain boundary diffusion as opposed to volume diffusion in micro-size powders. Under the same shcok strength, it was found that the ignition times in micro- and nano-Al particles are 32 µs and 8 ns, respectively. Conner and Dlott[105] have compared laser flash heating induced combustion behavior of 62 nm boron and 50 nm aluminum nanoparticle in Teflon using time-resolved emission spectroscopy. In case of boron, the oxide shell boiled off and the metal core melted while in case of aluminum, oxide shell melted and the metal core vaporized. The energy release from vapor Al/Teflon reactions had a lifetime of ~100 ps while it was ~200 ps for liquid B/Teflon. Bocanegra et al.[106] compared the burning behavior of micron (6 µm) and nano (250 nm) sized aluminum particle clouds for various concentrations in air. Nanoparticle cloud exhibited faster (difference ~10 cm/s) flame propagation as compared to microparticle cloud. The flame 31 front had two stages: an initially irregular stage followed by a stable stage with uniform slope. The flame temperature reached a maximum of 3300 ± 100 K for cloud with micron size particles and 2900 ± 100 K for cloud with nano-size particles. Jones et al.[107] have carried out experiments to study the combustion behavior of nano-aluminum (n-Al, 50 nm) with amorphous alumina shell and nano-aluminum oxide (n-Al 2 O 3 , 36 nm) particles stably suspended in ethanol. Heat of combustion was found to increase almost linearly with n-Al concentration above 5%. There was no enhancement in the heat of combustion for lower concentrations. n-Al 2 O 3 and heavily passivated n-Al did not react during the experiment. Gesner et al.[108] have investigated the effect of oxide shell thickness on the flame speed in nanoaluminum. They used Transmission Electron Microscope to observe the damages in the oxide shell on Al nanoparticles with different M values (M = ratio of particle radius to shell thickness) during heating at 480 °C for up to 180 minutes in a pure oxygen environment. Their experiments suggest that the damage in the oxide shell reduces flame speed and there is a weak correlation between the flame speed and M values. Yetter et al.[109] discuss the combustion of metastable intermolecular composite (MIC), their application, and their synthesis and assembly as well as combustion of Al nanocomposites in their review article. Al particles typically have an oxide coating with a limiting thickness of about 3 nm at room temperature. With a 100-nm diameter particle having a 3-nm-thick coating, the energy loss per unit volume due to the presence of the oxide layer is about 10%. A particle with a diameter of 10 nm with the same oxide layer thickness would have an energy loss per unit volume of approximately 60%. Composite materials that limit the volume of non-energetic material in the powders have been under development in recent years. Coating Al nanoparticles with self-assembled monolayers 32 (SAMs) can be a possible strategy. The applications of Al nanoparticles are in the area of solid propellants, solid fuels, and thermites, and the manufacturing of this material with highest possible energy yield is promising with multifunctional and smart capabilities. Though there have been numerous experimental studies investigating the oxidation behavior of aluminum nanoparticles, a comprehensive understanding about the underlying mechanisms is still lacking. Researchers have conducted numerical studies and proposed models to understand the dynamics of the oxidation reactions at the atomistic level. Puri et al.[17] studied structural changes, stress development and phase transformations in Al nanoparticles using molecular dynamics simulations. The nanoparticles were in the range 5 – 10 nm with 1 – 25 thick crystalline or amorphous oxide layer. The simulations were carried out in NVE and NPH ensembles using Steritz- Mintmire potential[110] . It was found that the oxide shell melts at 1100 K, which is much lower than the experimentally observed value of 2327 K. The thickness of the oxide shell had weak influence on the melting temperature of the shell. Structural analysis showed that the stresses in the shell were caused by the volume dilation in the core during melting. The shell was found to crack if it melts at a slow rate. Alavi et al.[111] used Steritz-Mintmire electrostatic plus (ES+) potential[110], which takes into account the effect of charge transfer between Al and O using the electronegativity equalization principle, to simulate the oxidation of Al nanoparticles using molecular dynamics. They used NVT ensemble and the equations of motions were integrated using the leapfrog algorithm with a timestep of 2 fs. The nanoparticles were annealed by heating till 3000 K and subsequently cooling to 500, 1000 or 1500 K. The Al:O ratios in 33 the nanoparticles were 1:1 and 2:1. The 1:1 simulation results in the formation of oxide shell, which stabilizes the shape of the nanoparticle and retains the original structure. In case of 2:1 simulation, they observe partial oxidation due to insufficient oxygen and the nanoparticle attains irregular shapes and rough surfaces. Perron et al.[112] have performed variable charge molecular dynamics simulations using Steritz-Mintmire electrostatic plus (ES+) potential[110] to study the oxidation of multi-grain nanocrystalline (mean grain size = 5 nm) aluminum surfaces in the temperature range 300 – 750 K. The authors have identified two regimes: linear, governed by interfacial processes, and logarithmic, governed by diffusion. The linear constant was found to independent of temperature, while the logarithmic constant increased with the temperature. It was found that the kinetics follows a direct logarithmic law governed by diffusion process and tends to a limiting oxide thickness of ~3 nm. Chung et al.[113] presented a model describing the enthalpy of reaction for the oxidation of Al nanoparticle as a function of the diameter of the nanoparticle. Their model includes the size dependence of the cohesive energy of the reactant particles, the size dependence of the product lattice energy, extent of product agglomeration, and surface capping effects. As per their model, the maximum energy release should occur for sub 10-nm Al nanoparticles. Pantoya and Granier[9] have carried out combustion studies of different size Al particles – from micron scales to nano scales. They propose a model mechanism to account for the increased combustion rates for the nanoscale system. In their model, diffusive transport plays the important role in combustion at micron scale, whereas they claim that at nanoscale the mechanism is not diffusive, in stead the oxide shell of the nanoparticle explodes into fragments resulting in non-diffusive transport of Al from the 34 core into the oxygen environment. This leads to a dispersion of Al clusters in oxygen resulting in higher combustion rates. Levitas et al.[16] have proposed a melt dispersion mechanism to explain the oxidation of Al nanoparticles with an alumina shell. As per this model, the change in the volume of the molten Al increases pressure and results in spallation of the oxide shell. A later unloading wave generates high tensile pressure and disperses the molten Al core into small clusters, which are then ejected out at a high velocity. In another study, Levitas and collaborators[19] have suggested melt dispersion mechanism as a possible explanation of the observed decrease in burn time after fast heating rate. Tradition diffusion mechanism operates for oxidation at slow heating rates (103 °C/s) but at heating rates > 107 °C/s, melt dispersion mechanism comes into play. Fast heating causes fast melting of the core, which leads to high pressure in the melt and spallation of the entire shell instead of local cracks. Fast fracture of the oxide shell generates high tensile pressure in the unloading wave. Although there is no direct evidence for the proposed mechanism, the author claims that the observed heating rate effect strongly supports the melt-dispersion mechanism. Campbell et al.[114] performed large-scale reactive molecular dynamics simulations of the oxidation of Al nanocluster of diameter 20 nm. The interatomic interactions are described by the Steritz-Mintmire electrostatic plus (ES+) potential[110] and the simulations were carried out in NVE and NPH ensembles. In case of NVE ensemble, a continuous increase in the thickness and temperature of the oxide layer, followed by melting and explosion of the nanocluster. In simulations done in NPH ensemble, a 4 nm thick passivating amorphous oxide layer is seen to form and diffusion was determined to be the rate-limiting step. 35 Recently, our group has developed a different potential for the interaction between aluminum and oxygen atoms as well as an interpolation scheme[28] for smooth transition from this potential to EAM potential developed by Voter-Chen[115] for Al-Al interactions. In our previous works we have carried out large scale reactive molecular dynamics simulations based on this new potential combination to have atomistic level understanding of the oxidation behavior of aluminum nanoparticles. Wang et al.[23] have found a thermal to mechano-chemical transition of mechanism during the oxidation of a nanoparticle consisting of a laser flash heated Al core and an Al2O3 shell. The intermediate reaction products change from Al-rich Al2O to oxygen-rich AlO2 clusters as the mechanism shifts from thermal diffusion to mechanically enhanced diffusion to ballistic transport. High initial temperature of the core (9000 K) results in catastrophic failure of the shell, leading to direct oxidation of the core and faster energy release. In a follow-up study, Wang et al.[24] have compared the effect of amorphous and crystalline oxide shells on the oxidation of an Al nanoparticle. In case of amorphous shell, oxidized nanocluster fragments are formed due to shattering of the shell. This, combined with fragmentation and dispersion of the nanoparticle, results in faster oxidation reactions as compared to a nanoparticle with a crystalline shell. The crystalline shell expands and then contracts resulting in considerable amount of Al core covered by the shell. Using the same potential scheme, Clark[28] has studied the oxidation behavior of an Al 40 nm nanoparticle with a 3 nm amorphous alumina shell. The nanoparticle is heated to 1100 K and simulated in an NVE ensemble. Three distinct stages are identified during the oxidation process. In stage 1, the atoms from the shell penetrate the heated core and more heat is released from the exothermic reactions. The heat from the reactions at the core- 36 shell interface results in melting of the alumina shell initiating stage 2. The rate of change of temperature increases in stage 2 accompanied by deformation of the nanoparticle. Stage 3 begins when the temperature of the shell increases 3000 K. Clusters are seen to eject from the alumina shell, followed by those from the core into the surrounding oxygen. Temperature increases throughout the simulation as a result of the exothermic reactions inside the core. In this work, we have expanded this study to a chain of three identical nanoparticles arranged in a linear fashion. The left and right nanoparticles are heated and then the system is left to react in an NVE ensemble. 3.2 Simulation We simulate a linear chain of three ANPs, where each ANP consists of an aluminum core of 400 Å diameter covered with a 30 Å amorphous alumina (Al2O3) shell (Figure 3.1a). To create one ANP, a spherical core, 400 Å in diameter, is cut from an fcc lattice of aluminum. A slab of amorphous alumina is created from an initially crystalline slab using the melt-quench method.41 A 30 Å thick alumina shell, with inner diameter 400 Å, is then cut out of this slab of amorphous alumina. The aluminum core and the alumina shell are then combined together, and run under NVE ensemble until an interface is formed between the core and the shell. Three such ANPs are placed 5 Å apart in a linear fashion (Figure 3.1b). This system is then placed in an oxygen environment. The size of the MD box is 1938.4 × 969.2 × 969.2 Å3. Each ANP contains 1,973,159 atoms in the Al-core, 856,475 Al atoms and 1,284,281 O atoms in the alumina shell. The surrounding oxygen has 1,888,193 atoms, and the entire system has 14.23 million atoms. We use MPI-based parallel algorithm[87] to simulate this system on 512 processors. The outer ANPs (left 37 and right) are slowly heated to a temperature of 1200 K using velocity scaling. The central ANP and background oxygen are kept at 500 K. After heating, the system is left to react in NVE ensemble for 1 ns. The equations of motion are integrated using velocity- Verlet algorithm[86] with a timestep of 1 fs. Figure 3.1: The system used in the simulation. (a) A single ANP, consisting of a spherical aluminum core, 400 Å in diameter, surrounded by a 30 Å thick spherical alumina shell. The aluminum core is shown in grey, the aluminum atoms in the shell are shown in red while the oxygen atoms in the shell are shown in yellow. (b) A slice of the system through the central place showing the linear arrangement of the three ANPs. The aluminum and oxygen atoms in the shell of the hot ANPs (left and right) are shown in red and yellow respectively, while the aluminum and oxygen atoms in the shell of the cold ANP (central) are shown in red and yellow, respectively. The aluminum atoms in the cores of the three ANPs are shown in grey. 38 3.2.1 Interatomic Potential The interactions between aluminum atoms are described by the Voter-Chen EAM potential[115], where the total energy is given by: !"! = ! (3.1) The energy per atom is in two parts: ! = ! ! ∅ !! +( ! ! ) (3.2) ! = !" ! (3.3) Where the density function takes the form: !" = ! ( !!" +2 ! !!!" ) (3.4) The pairwise potential is a Morse potential: ∅ !" = ! 1− !! !(!!! !) ! − ! (3.5) ! = 3.7760 ; ! = 2.1176 Å; ! = 1.4859 Å !! ; ! = 5.555 Å; = 3.3232 Å !! The potential as well as density functions are shifted and smoothened at the cut-off distance (R c ). In this implementation, the density function, F(ρ) is determined numerically to match the “Rose Equation”: ! ∗ =− ! ∗ (1+ ∗ ) !! ∗ (3.6) ∗ = ( ! ! ! !!) ! !"! /!!" (3.7) Where, E U =cohesive energy of the scaled system, B=Bulk modulus = 79.38 GPa E 0 =E coh =equilibrium cohesive energy = -3.36 eV for Al, Ω= equilibirum atomic volume a=lattice constant, a 0 =equilibrium lattice constant=4.05Å for Al 39 The potential form used for alumina consists of a two-body term and a three-body term to incorporate various effects, including steric repulsion, Coulomb interaction, charge- dipole interaction, van der Waals interaction, and covalent bond bending and stretching. The total energy of the system is written as ∑ ∑ < < < + = k j i ik ij ijk j i ij ij r r V r V E ) , ( ) ( ) 3 ( ) 2 ( tot   , (3.8) where j i ij r r r    − = with i r  being the position of the i-th atom, and | | ij ij r r  = is the distance between atom i and j. In Eq. (3.8), the first, two-body, term sums over neighbor atoms, j, within a cutoff range r c from atom i, and has the form, V ij (2) (r)= H ij r η ij + Z i Z j r e −r r 1s − D ij 2r 4 e −r r 4s − w ij r 6 , (3.9) whereas the second, three-body, term, sums over atoms, j and k, that are covalently bonded with atom i within a distance r 0 , and the form of this three-body interaction is: ( ) 0 2 0 2 0 0 0 ) 3 ( , ) cos (cos 1 ) cos (cos exp ) , ( r r r C r r r r B r r V ik ij ijk ijk ijk ik ij ijk ik ij jik ≤ − + − ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − + − = θ θ θ θ ξ ξ   . (3.10) In Eq. (3.9), H ij is the steric repulsion strength, Z i is the effective charge in units of the electronic charge, D ij represents the charge-dipole strength, w ij is the strength of van der Waals attractions, η ij is the exponent of the steric repulsion, and r 1s and r 4s are the screening lengths for the Coulomb and charge-dipole interactions, respectively. In Eq. (3.10), B ijk is the strength of the three-body interaction, ξ and C ijk are constants, and θ ijk is the angle formed by  r ij and  r ik . For computational efficiency, the 2-body pair potential is truncated at r c and shifted for r < r c , so that the potential and its first derivative are continuous at r c : 40 [ ] ⎪ ⎩ ⎪ ⎨ ⎧ ≥ < − − − = = ) ( ) ( 0 / ) ( ) ( ) ( ) ( ) 2 ( ) 2 ( ) 2 ( ) shifted , 2 ( c c r r ij c c ij ij ij r r r r dr dV r r r V r V r V c , (3.11) The potential parameters for alumina have been determined to fit experimental values of bulk modulus, specific heat, and melting temperature of alumina. The parameters used in the simulation are given in Table 3.1 below. Al O Z i (e) 1.5237 -1.0158 r 1s = 5.0 Å r 4s = 3.75 Å r c = 6.0 Å e = 1.602×10 -19 C Two body Al-Al Al-O O-O η ij 7 9 7 H ij (eV Å η ) 12.4236 244.3038 553.3917 D ij (eV Å 4 ) 0 3.4374 1.52775 W ij (eV Å 6 ) 0 0 63.4894 Three body B ijk (eV) θ o (deg) C ijk ξ (Å) r 0 (Å) O-Al-O 12.4844 90.0 10 1.0 2.90 Al-O-Al 8.1149 109.47 10 1.0 2.90 TABLE 3.1: Parameters for two- and three-body part of our potential used for MD simulation of alumina. A detailed description of the interatomic potential used in this simulation and its validation can be found in Ref. #[28]. 41 3.3 Oxidation of outer ANPs Figure 3.2: Snapshots of the central slice of left ANP at (a) 0 ns, (b) 0.1 ns, (c) 0.2 ns, and (d) 0.3 ns. The aluminum core is shown in grey, the aluminum atoms in the shell are shown in red and the oxygen atoms in the core are shown in yellow. The reactions begin at the core-shell interface, resulting in the diffusion of shell atoms into the aluminum core. By 0.3 ns, the outer ANPs expand and undergo considerable deformation. Amorphous alumina was created by melt-quench method.[80, 116] The ANPs were then created by putting aluminum core inside alumina shell. They were then thermalized at 500 K till aluminum-alumina interface[117] was formed. Using velocity scaling, we heat 42 the outer ANPs to a temperature of 1200 K, keeping the central ANP and background oxygen at 500 K. At this stage, the aluminum core, with melting temperature of 900 K, is molten while the alumina shell is solid as it melts at 2300 K. As the system is left to react under NVE ensemble, the outer ANPs undergo oxidation reactions first. Figure 3.2 shows the snapshots of the left ANP at (a) t = 0 ns, (b) t = 0.1 ns, (c) t = 0.2 ns, and (d) t = 0.3 ns. Owing to dangling bonds at the aluminum-alumina interface, the oxidation reactions are localized at the core-shell interface during the first 0.1 ns (Figure 3.2b). As the reactions continue, heat is generated inside the outer ANPs and initiates further exothermic reactions between the atoms from the core and the shell. As a result, the shell begins to diffuse into the core (Figure 3.2c). As the shell atoms reach deeper into the core, they undergo oxidation reactions with the core atoms and continue to generate heat, resulting in expansion and deformation of the outer ANPs (Figure 3.2d). The average temperature of the outer ANPs increase from 1200 K at t = 0 ns to 3500 K at t = 0.3 ns. 3.3.1 Radial distribution of temperature The behavior of the outer ANPs is very close to that observed in a similar study of single ANPs.[28] To track the flow of heat in outer ANPs in the initial stages, we divide the outer ANPs in spherical bins of width 5 Å, starting from the respective center, and the average temperature of each bin at 0, 0.2 and 0.4 ns is plotted against the distance from the center in Figure 3.3. At t = 0 ns, the alumina shell is at a higher temperature than the aluminum core, with the temperature gradually decreasing as we move towards the center of the ANP (red line). During the first 0.2 ns (blue line), the increase in the region of core-shell interface is higher than that inside the core (shown by double-ended arrows), indicating that the reactions are localized to the core-shell interface. As the reaction 43 continues, the shell expands into the core and by 0.4 ns (green line) the shell reaches the center of the core, i.e. the shell merges into the core completely. This shows that the reactions begin at the core-shell interface and gradually proceed inwards until the entire ANP has been oxidized. Figure 3.3: Radial distribution of temperature inside left ANP. The ANP is divided into 5 Å thick spherical bins and the temperature of atoms falling in each bin is calculated and plotted against the distance from the center of the ANP. The red line represents the initial state. The distribution at 0.2 ns is shown in blue, while the same at 0.4 ns is shown in green. The double-sided arrows indicate the increase in temperature by 0.2 ns. By 0.4 ns, the temperature across the ANP reaches 4000 K. 44 3.3.2 Rate of change of temperature Figure 3.4 shows the temperature (blue line) and the rate of increase in the temperature, dT/dt, (red line) of the particles in outer ANPs versus time. The dT/dt curve can be divided into three sections, marked by green dash-dot line, representing the three stages of oxidation in outer ANPs. The same three stages were observed in the study of oxidation of a single ANP.41 In Stage I (from the beginning to the first inflexion point), we see a declining rate of increase in temperature. This is because during that period, the outer ANPs undergo localized heating at the core-shell interface. Stage II begins once the temperature of the alumina shell reaches the melting point (purple dashed line, representing T=2300K, around which alumina shell melts). We see an increasing rate of increase in temperature in this stage, as the molten alumina shell then expands into the core by reacting with the aluminum atoms in the core. Being an exothermic reaction, this generates heat, which initiates more such reactions thereby making the oxidation of the outer ANPs self-sustaining. After t = 0.3 ns, dT/dt begins to decrease again, marking the onset of Stage III. During this stage, we observe ejections from the outer ANPs into the background oxygen environment (shown later in Figure 3.7). As atoms eject out of the outer ANPs, the total ‘fuel’ remaining inside the ANP decreases. This reduces the rate of increase in temperature of the outer ANPs, helped further by the fact that the background oxygen environment is less dense than the alumina shell. It should be noted that while the rate of the change of temperature changes, the average temperature of the outer ANPs keep increasing throughout the simulation. 45 Figure 3.4: Temperature and rate of increase in temperature, dT/dt, vs time in outer ANPs. The blue curve represents temperature vs time, with the temperature axis on the right. The rate of increase in temperature (dT/dt) vs time is shown in red, with dT/dt axis on the left. The green dash-dot lines demarcate the three stages of oxidation, while the purple dashed line represents the melting temperature of alumina, T = 2300K. 3.3.3 Aluminum Oxide clusters Numerous aluminum oxide clusters, with varying compositions, form during the oxidation reactions. We list the nearest neighbors of each atom, subject to distance-based criteria corresponding to the respective bond lengths. We use breadth-first search on this neighbor list to identify independent clusters. 46 Figure 3.5: Distribution of clusters of compositions equivalent to (a) Al2O, (b) AlO, (c) AlO2, (d) Al2O3, (e) Al2O5, and (f) AlO3 inside the left ANP with time. Al-rich clusters (a) and (b) are observed from the beginning while O-rich clusters (c) – (f) appear after 0.3 ns. Al-rich clusters are considerably more in number than O-rich clusters 47 For each cluster, we calculate the composition and group them with the closest of the standard aluminum oxide compositions. Figure 3.5 shows the distribution of these clusters in the outer ANPs during oxidation. We see that Al-rich clusters (Al 2 O and AlO) outnumber O-rich clusters (AlO 2 , Al 2 O 3 , Al 2 O 5 and AlO 3 ). We observe a lot of Al 2 O clusters, along with small number of AlO clusters in the first 0.3 ns. It should also be noted that the number of each of these types of clusters increase with time. The fact that during the first 0.3 ns, only Al-rich clusters are formed indicates that oxygen from the shell is penetrating into the core i.e. the reaction is localized on the core-shell interface. O-rich clusters appear only after 0.3 ns. As shown in Figure 3.4, stage III in outer ANPs begins at t = 0.3 ns with the ejections of atoms from the outer ANPs into the background oxygen. This results in the sudden increase in the number of oxygen-rich clusters (AlO 2 , Al 2 O 3 , Al 2 O 5 and AlO 3 ). We also note that the number of such clusters increase with time. This is consistent with the continuous ejections from the outer ANPs as the pressure inside the ANPs, coupled with a weakened shell force the atoms out of the nanoparticles. Though the number of oxygen-rich clusters increases with time, they are still one order of magnitude less than the number of Al-rich clusters at any given point of time. This is because the background oxygen is less dense than the core and while more and more ejected atoms react with the background oxygen to form oxygen-rich cluster, the majority of the oxidation reactions occur inside the ANP itself. This is consistent with the declining rate of change of temperature in stage III, as seen in Figure 3.4. 48 Figure 3.6: Radial distribution of composition in outer ANPs during the simulation. Composition is calculated as percentage of oxygen in 5 Å thick spherical bins. %Oxygen is plotted against distance from the center of ANP. Figure 3.6 shows the radial distribution of composition inside the left ANP, measured in terms of %Oxygen. We divide the ANP into 5 Å thick spherical bins and calculate the percentage of oxygen atoms in each bin. This number is then plotted against the distance from the center. It can be seen that %Oxygen increases with time as one moves closer to the center of the ANP, indicating that oxygen from the shell is penetrating into the core. 49 At the end of the simulation (1 ns), the entire ANP attains a uniform composition of 40% oxygen. Table 3.2 shows the number of atoms in the largest cluster of different compositions, at different times during the oxidation process. We can see that the alumina shell, defined as the largest cluster, holds its composition in the first 0.2 ns, after which its composition changes to AlO. This implies that the shell gradually becomes Al-rich which can happen in two ways: either oxygen from the shell penetrates into the core, thus expanding the shell or oxygen from the shell ejects out of the shell, leaving the shell oxygen-deficient. The second possibility is discarded by the fact that Al-rich clusters were observed in the first 0.2 ns of reaction which can only happen if oxygen penetrated into the Al-core. Another noteworthy point is that Al-rich clusters are larger in size than the oxygen-rich ones. As the reaction continues, larger clusters are formed, especially Al-rich and Al 2 O 3 clusters. Type \ Time (ns) -> 0.1 0.2 0.3 0.4 0.5 0.6 Al 2 O 71 130 225 211 241 259 AlO 44 56 2001417 2495837 2483876 2434702 Al 2 O 3 2282993 2356584 18 50 68 77 AlO 2 6 11 9 15 25 28 Al2O 5 7 7 7 10 13 13 AlO 3 8 5 5 9 8 9 Table 3.2. Maximum size of clusters in Left ANP. The largest cluster is defined as the shell. Due to reactions between the shell and core atoms, the composition of the shell changes from that of Al2O3 to AlO after 0.2 ns. 50 3.4 Oxidation of central ANP Figure 3.7 shows snapshots of the atomic configuration during the simulation. The outer ANPs are the first ones to oxidize. The oxidation reactions in the central ANP are limited to the fused zones of contact during the first 0.3 ns (Figure 3.7a). Initially, the cores of the outer ANPs are molten while the shell is solid at 1200 K. The shell atoms of the outer ANPs react with their corresponding core atoms near the core-shell interface. Heat is generated inside the outer ANPs as well as the zones of contact due to the exothermic nature of such oxidation reactions. As the reactions continue, the shells of outer ANPs are seen to diffuse into the corresponding cores (Figure 3.7b). By 0.3 ns, the shells of the outer ANPs are completely molten (Figure 3.4). During the initial heating, the central ANP was kept at 500 K. At this temperature, both the aluminum core and the alumina shell are solid. Therefore, there are very limited reactions inside the central ANP in the regions away from the zones of contact with outer ANPs. After melting, the shells of outer ANPs begin to penetrate into central ANP through these zones of contact. Inside the central ANP, the shell region in the zones of contact melts first, which is then pushed into the core by the expanding outer ANPs. By 0.3 ns, half of the core of the central ANP is penetrated (Figure 3.7c). In addition to bringing external heat, the atoms from outer ANPs initiate exothermic oxidation reactions with the core of the central ANP. Hence, the central ANP gains heat due to a combination of two factors: 1) convective heat transfer, and 2) exothermic oxidation reactions involving atoms in the core and atoms from the shell and outer ANPs. By 1 ns, the three ANPs fuse together to form an ellipsoidal aggregate (Figure 3.7d). 51 Figure 3.7: Snapshots of the system at (a) 0 ns, (b) 0.3 ns, (c) 0.6 ns, and (d) 1 ns. Only the central slice of shells of the three ANPs are being shown. Aluminum atoms in the outer (left and right) ANPs are shown in red, oxygen atoms in outer ANPs are in yellow; aluminum atoms in the central ANP are in cyan and oxygen atoms in the central ANP are shown in blue. Penetration into the central ANP is observed and the system forms an ellipsoidal aggregate by 1 ns. 52 Figure 3.8: a) Aluminum oxide clusters inside the central ANP are grouped based on their positions. Region I (red) represents a right circular cone of aperture = 10° along Z-axis, which is perpendicular to the zones of contact, while region II (blue) represents that along the X-axis, which go through the zones of contact. b) Histogram showing %increase in clusters in central ANP inside region I (red) and region II (blue) at different times. More clusters are observed in the zones of contact as compared to other regions inside the central ANP. To identify the localized regions of reactions, the surface of the central ANP was divided into small bins, each having 5° as azimuthal and inclination angles (spherical 53 coordinates). The number of atoms in the central ANP, which have reacted and formed clusters, is calculated for each bin over time We compare such clusters in two regions, as shown in Figure 3.8a. Region I (red) is a right circular cone around z-axis, perpendicular to the axis of the linear chain, with origin at the center of the central ANP and aperture = 10°, whereas region II (blue) is covers the same volume along the axis of the chain, through the zones of contact. Inside region I, we find a 19% increase in the first 0.05 ns, 33% by 0.1 ns, and 53% by 0.2 ns (Figure 3.8b). Comparing this with region II, we notice a 21% increase in the first 0.05 ns, 43% by 0.1 ns, and 77% by 0.2 ns. This indicates that there are more reactions in the zones of contact as compared to the rest of the core-shell interface in the central ANP. 3.4.1 Temperature profile Figure 3.9 shows the temperature across the central slice of the system at different stages of simulation. We divide the system into cubic voxels of width 5 Å and calculate the temperature of atoms residing in each voxel. The value of each voxel is then averaged over the three neighboring voxels in each direction. The resulting value of temperature for each voxel is then displayed according to the color code shown on the right where blue represents temperature below 500 K and red represents temperature above 6000 K. It can be seen that initially the shells of the outer ANPs are at around 1400 K, while the corresponding cores are at around 1000 K, and the central ANP is at around 500 K (Figure 3.9a). As the reactions being and progress inside the outer ANP, we see that by 0.3 ns, the shell of the outer ANP reach temperatures higher than the melting point of 2300 K (Figure 3.9b). By this point, we also see increase in temperature of regions close to the zones of contact inside the central ANP. 54 Figure 3.9: Temperature across the central slice of the system at (a) 0 ns, (b) 0.3 ns, (c) 0.6 ns, and (d) 1 ns. The scale is shown on the right with blue representing temperatures below 500 K and red representing temperatures above 6000 K. 55 This is due to the reactions caused by the molten shells of outer ANPs in the fused zones of contact. By 0.6 ns, the outer ANPs attain temperature close to 5000 K (Figure 3.9c). Due to convective heat transfer into the central ANP, we see that the temperature inside the central ANP decreases as we go towards the center. By 1 ns, all the three ANPs fuse together and reach temperatures above 5500 K (Figure 3.9d). To further quantify the variation of temperature across the system, we divide the system into 5 Å thick slices along the X-axis, which is the axis along which the ANPs are arranged. The temperature of atoms inside each bin is calculated and plotted against the distance along the X-axis at different times during the simulation (Figure 3.10). Initially, we can clearly see the temperature difference between the outer ANPs and the central ANP (red curve). Since, the outer ANPs are the first ones to oxidize following the pre- heating, they gain much more heat as compared to the central ANP from 0 to 0.4 ns (blue curve). As the outer ANPs penetrate the central ANP, they bring in external heat through convective transfer. In addition, they initiate exothermic oxidation reactions inside the central ANP. As a result, the central ANP also begins to gain heat as reflected by the green curve representing the state at 0.6 ns. As seen earlier, all the ANPs fuse together into an ellipsoid aggregate by 1ns. The average temperature of this aggregate at the end is 5500 K, as seen in the brown curve representing the final state. 56 Figure 3.10: Linear distribution of temperature along the X-axis. Heat flows from the outer ANPs to the central ANP through the zones of contact. The core of the central ANP is the last region to be heated. 3.4.2 Temperature gradient To understand the mechanism of oxidation of the ANPs, we calculated the rate of increase in the temperature of the ANPs with respect to time (dT/dt). Temperature gradient is obtained as the difference in the temperature at two consecutive instances 5 ps apart divided by the time interval between the measurements i.e. 5 ps. For the outer ANPs and the central ANP, the value of dT/dt thus obtained is plotted against the instantaneous time in Figure 3.11a and 3.11b, respectively. As observed in the outer ANPs (Figure 3.4), the central ANP also exhibits three stages of oxidation. These stages, identified by the 57 inflexion points in the corresponding dT/dt curves, are denoted by I, II and III in Figure 11. Though both, the outer ANPs and the central ANP, undergo three stages during oxidation, the reaction dynamics in the central ANP is different from that of the outer ANPs. In case of outer ANPs (Figure 3.11a), dT/dt decreases with time as the oxidation reactions are localized to the core-shell interface. Melting of the alumina shell initiates stage II during which dT/dt increases with time. The outer ANPs expand and deform in this stage (Figure 3.2). This leads to weakening of the shell resulting in ejections of atoms from the outer ANPs into the background oxygen (red dashed curve in Figure 3.11a). The beginning of such ejections marks the onset of stage III during which we observe sharp decline in dT/dt with time. This is due to the fact that the background oxygen is less dense in comparison to the ANP, which results in decline in the rate of heating with time. It should be noted that the temperature of the ANPs always increases throughout the oxidation (Figure 3.4). In case of central ANP, stage I is observed during the first 0.1 ns which is also the case with outer ANPs. However, the rate of reaction during this period in the central ANP is much slower as compared to that in the outer ANPs. Initially the reactions in the central ANP are localized to the fused zones of contact (Figure 3.8b). When the alumina shells of the outer ANPs melt, they begin to heat the shell of the central ANP through the zones of contact at a faster rate. This initiates stage II in the central ANP during which we see that dT/dt increases at a much faster rate with time, as compared to outer ANPs. It is also observed that the stage II in central ANP has a larger time span and higher peak rate than those in outer ANPs. In stage II, atoms from the outer ANPs penetrate the central ANP through the zones of contact, pushing the shell of the central ANP towards the center of the central ANP in the process (Figure 3.7c). This 58 results in more oxidations reactions inside central ANP in stage II as compared to those in outer ANPs. At ~0.6 ns, dT/dt in central ANP begins to decrease with time, thus entering stage III. As seen in the case of outer ANPs, the onset of stage III in the central ANP coincides by the beginning the ejections of atoms out of the central ANP (red dashed curve in Figure 3.11b). Figure 3.11: Temperature gradient (dT/dt) (blue) and the number of ejections from the alumina shell (red) in (a) left and (b) central ANPs. The temperature gradient curve corresponds to the axis on the left while the ejections curve corresponds to the axis on the right. The onset of stage III, in both the cases, coincides with the beginning of ejections of atoms form the shell. 3.4.3 Reaction Heat To measure the amount of heat generated inside the central ANP, we track the atoms coming inside as well atoms going out of the central ANP in one timestep. We then calculate the reaction heat generated as E rxn = KE increase -KE new + KE out , where KE increase is the increase in kinetic energy of the central ANP, KE new is the kinetic energy of the atoms entering the central ANP at that instant, and KE out is the kinetic energy of the atoms exiting the central ANP at that instant. Figure 3.12 shows the progress of 59 calculated reaction heat with time. We observe that reaction heat generated within the central ANP increase with time during the entire oxidation process. During the initial 0.4 ns, the rate of generation of heat is slow. This is because initially the oxidation reactions inside the central ANP are localized to the fused zones of contact. As the outer ANPs penetrate into the central ANP, they initiate more oxidation reactions within the central ANP. This is reflected in the increase in the rate of generation of heat between 0.4 and 0.6 ns. After 0.6 ns, we see that the rate decreases again as the central ANP enters into stage III of oxidation. Figure 3.12: Reaction heat per particle generated inside the central ANP vs time. Initially, the rate of generation of heat is slow as the reactions are localized to the zones of contact. As the outer ANPs penetrate the central ANP and initiate more oxidation reactions, thereby increasing the reaction rate between 0.4 and 0.6 ns. After 0.6 ns, the rate declines again as ejections of atoms from central ANP into the background oxygen begin and the central ANP enters stage III of oxidation. 60 3.4.4 Penetration into central ANP Figure 3.13: Radial distribution of outer atoms inside the central ANP. After 40.4 ns (brown), the outer atoms are still close to the alumina shell; by 0.6 ns, (green) the outer atoms have penetrated half the radius of the central ANP and by 0.8 ns (red), the penetration reaches 80 Å. By 1 ns (blue), the outer atoms can be seen to be present at the core. As oxidation reactions progress inside the outer ANPs, they deform and expand (Figure 3.2). Once the shells of the outer ANPs melt, they begin to penetrate the central ANP through the zones of contact (Figure 3.7). To obtain the extent of penetration at any given instant during the simulation, we track the atoms, originally belonging to the outer ANPs, inside the central ANP. Based on their position, these outer atoms are distributed in 5 Å thick concentric spherical bins starting from the center of the central ANP and expanding outwards. The number of such atoms in each bin is plotted against the distance from the 61 center of the central ANP at different times is plotted in Figure 3.13. It can be seen that outer atoms are present at a minimum distance of 190 Å from the center of the central ANP at 0.4 ns (brown curve), which is very close to the zones of contact. By 0.6 ns (green curve), such atoms can be found at a distance of 150 Å. In the next 0.2 ns (red curve), the outer ANPs penetrate further and we can find outer atoms at a minimum distance of 50 Å from the center. By 1 ns (blue curve), they are present throughout the core, indicating complete penetration. Figure 3.14: Position of different layers of the penetration front inside the central ANP vs time. The atoms from the shell of the central ANP form the first layer (brown and magenta), the atoms from the shell of the outer ANPs from the second layer (red and blue) and atoms from the core of the outer ANPs form the third and last layer (green). The penetration fronts maintain their relative position with respect to the core of central ANP. 62 To identify the composition of the penetration front, we track the positions of atoms of different types (shell or core) and different origins (outer or central ANPs) and distribute them in 5 Å thick concentric spherical bins with the common center at the center of the central ANP. The position of innermost non-empty bin is plotted against time for each type and origin of atoms in Figure 3.14. It is evident that the different species maintain their relative positions during the penetration. The atoms from the shell of the central ANP constituted the first layer of the penetration front, followed by the atoms from the shells of the outer ANPs. The atoms from the core of the outer ANPs formed the third and last layer of the penetration front. This means that the shell of the central ANP does not break during penetration. The outer ANPs push the shell but not break through it as they penetrate the central ANP. As the atoms from the outer ANPs penetrate through the zones of contact, they form a convex front. To get the speed of penetration, we calculated the average position of the innermost atoms, belonging to tip of the penetration front. This position, which corresponds to the depth of penetration, is plotted against time in Figure 3.15. The slope of this curve gives the average speed of penetration as 54 m/s. In their studies on the combustion behavior of aluminum nanoparticles, Pantoya and Granier34 have reported a wide range of combustion velocities. It was found that the combustion velocity depends on the size and density of the aluminum nanoparticles. As the reaction pathways shift from solid oxidizer to ambient oxygen gas, the combustion velocity increases from 1 m/s to 1000 m/s. Our results fall within the observed range and may correspond to the case of high-density aluminum nanoparticles reacting with ambient oxygen. 63 Figure 3.15: Maximum position of outer atoms inside the central ANP vs time. The average speed of penetration is 54 m/s. 3.5 Summary In summary, we have used molecular dynamics to simulate the oxidation dynamics in a multi-million-atom system of three spherical ANPs. Each ANP consists of a spherical aluminum core, 400 Å in diameter, surrounded by a 30 Å thick amorphous alumina shell. The ANPs were arranged in a linear fashion and the left and the right (outer) ANPs were heated to a temperature of 1200 K. The central ANP and the background oxygen were kept at 500 K. Due to the pre-heating, the outer ANPs undergo exothermic oxidation reactions first, during which they expand and deform while gaining heat Initially, the oxidation reactions inside the central ANP are localized to the fused zones of contact. 64 Once the shells of the outer ANPs melt, they begin to penetrate the central ANP through the zones of contact. In addition to bringing in external heat via convective transfer, the outer atoms also initiate exothermic oxidation reactions inside the central ANP. Thus, although slow in the beginning, the central ANP oxidizes at a faster rate as compared to the outer ANPs. Both, the outer ANPs and the central ANP, show three stages of oxidation. Stage I is lasts 0.1 ns, when the rate of reaction decreases with time due to localized oxidation. Stage II is characterized by increasing rate of increase in temperature due to increased number of oxidation reactions. Stage III, with decreasing rate of increase in temperature, coincides with the ejections of atoms from the ANPs into the background oxygen. The penetration front consists of three layers: 1) atoms from the shell of the central ANP; 2) atoms from the shells of the outer ANPs; and 3) atoms from the cores of the outer ANPs. The layers of the penetration front maintain their relative position with respect to each other. This means that the outer ANPs push the shell of the central ANP towards the center of the central ANP but they do no break through it. The speed of penetration is calculated as 54 m/s, which is within range of experimentally observed flame propagation rates. By 1 ns, all the three ANPs fuse together into an ellipsoidal aggregate at an average temperature of 5500 K. From this study, we can say that oxidizing ANPs have a huge effect on the reaction dynamics of a neighboring ANP. By pre-heating one ANP, it is possible to initiate oxidation reactions in a neighboring ANP to a level of self-sustenance. Since proximity to other ANPs and their pre-condition can change the dynamics of oxidation of an ANP, these factors should be given appropriate consideration in any study, theoretical or experimental, to understand the mechanism of oxidation of ANPs. 65 Chapter 4. Shock-Induced Collapse of Water Nanobubble 4.1 Background Cavitation bubbles are generated when there is a rapid pressure drop inside hydrodynamic systems. Sudden change in pressure results in localized change in velocities inside the fluid, which lead to the formation of nanobubbles. Such bubbles are observed in propulsion systems, pumps, cooling systems in power plants, hydraulic mining machinery, ultrasound, etc. Researchers have been studying the cause and effects of cavitation for a long time. In a review article, Arndt[29] has identified the major factors involving the occurrence of cavitation, performance breakdown, and effects such as vibration, noise and damage. Cavitation is defined as the formation of vapor phase in a liquid. The formation of individual bubbles is related to reductions in pressure to some critical value. While generally considered as an unwanted side effect of fluid flow, cavitation also has useful applications, for example in homogenization of milk and industrial cleaning. The author has noted that dissolved gas can influence hydrodynamic loads in cavitating flows and should be properly investigated. Lin et al.[35] have carried out experiments to investigate cavitation erosion resistance of PHT electroless Ni-P coatings. The generation and collapse of bubbles at the surface of a component, thereby creating sudden changes of pressure, can lead to cavitation erosion. The pressure pulses, caused by the collapsing bubbles, result in the emission of shock waves and micro-jets. This leads to fatigue, fracture and material loss. They performed cavitation erosion test on stationary specimen, varying in micro- and nano-scale SiC content and extent of heat- treatment. It was found that the specimen incorporated with nano-SiC particles as well as 66 subjected to post-heat treatment offered the most resistance to cavitation erosion. Richman and McNaughton[37] have identified correlations between cavitation erosion behavior and mechanical properties of metals. As a result of cavitation erosion being a fatigue process, the predominant determinant of cavitation erosion resistance of a material is its fatigue resistance coefficient. The authors have combined fatigue resistance coefficient with cyclic strain-hardening exponent to further improve the correlation. Bourne and Field[53] studied the collapse of cylindrical cavities on solid surfaces using high-speed streak and framing photography. They induced asymmetric collapse using shock waves of varying amplitudes leading to formation of high-speed jets. The impact of such jets as well as strong compression waves formed on the rebound of the bubble were the most important features affecting the damage on the solid surface. Metals with low compressive but higher tensile strengths underwent plastic deformation resulting in a cavitation pit, where as brittle materials deformed by cracking. Lauterborn and Ohl[118] investigated the dynamics of optically and acoustically produced cavitation bubbles using high-speed photography with framing rates upto 20 million frames per second. They reported light emanating in form of filaments when the cavitation bubble field was observed in total darkness with dark-adapted eye, thus confirming sonoluminescence during cavitation bubble collapse. It was also found that asymmetric bubble collapse produces three shock waves: two jet-induced and one induced by bubble compression. Lindau and Lauterborn[55] have revisited the collapse and rebound of a cavitation bubble near a wall with modern experimental means. The bubble is generated by the optical breakdown of the liquid when a strong laser pulse is focused into water. Observations are made with high- speed cinematography; framing rates range between several thousand 67 and 100 million frames per second, and the spatial resolution is in the order of a few micrometres. After formation the bubble grows to a maximum size with a radius of 1.5mm at the pulse energy used, and in the subsequent collapse a liquid jet evolves on the side opposite the wall and penetrates through the bubble. Using a shadowgraph technique and high framing rates, the emission of shock waves, which is observed at minimum bubble size, is resolved in detail. For a range of stand-off distances between the bubble centre and the wall, a counterjet forms during rebound. The counterjet is clearly resolved to consist of cavitation micro-bubbles, and a quantitative measure of its height evolution is given. Its emergence might be caused by a shock wave, and a possible connection of the observed shock wave scenario with the counterjet formation is discussed. No counterjets are observed when the stand-off distance is less than the maximum bubble radius, and the bubble shape becomes toroidal after the jet hits the wall. The jet impact on the wall produces a pronounced splash, which moves radially outwards in the space between the bubble and the wall. The volume compression at minimum bubble size is found to depend strongly on the stand-off distance. In his study of erosion by water- hammer, Cook[119] has identified that the final act of erosion is due to water-hammer, with the concentration having the effect of increasing the velocity in the final stage of the collapse and making the water-hammer more intense. The water-hammer pressure generated by a nanojet is given as P = ρcv, where ρ = density of the jet, c = speed of sound in the jet and v = speed of the jet. Brujan et al.[30] have experimentally investigated the interaction of a laser-induced cavitation bubble with an elastic boundary and its dependence on the distance between the bubble and boundary. They used a transparent polyacrylamide gel with 80% water concentration and elastic modulus of 0.25 68 MPa. They observed bubble splitting, formation of liquid jets and jet-like ejection of boundary material into the liquid using high-speed photography with upto 5 million frames/s. The maximum jet speed was measured to be 960 m/s, which was found to penetrate the elastic boundary through a 0.35 mm thick water layer. Heathcock et al.[31] have studied cavitation erosion resistance of various kinds of steel. The mechanisms of resistance and modes of erosion were investigated using optical microscopy, scanning and transmission electron microscopy and X-ray diffractometry. It was found that martensitic steels have the highest erosion resistance. The erosion resistance of austenitic stainless and manganese steels depends on their work hardening rate. Ferritic stainless steels had very poor erosion resistance owing to the sensistivity of bcc structure to high strain rates. Johnson et al.[32] have investigated cavitation erosion resistance of various nickel aluminides. They performed cavitation erosion tests in distilled water at 20 ° C with an ultrasonic transducer coupled to a titanium-tipped exponential horn vibrating at 20 kHz. Alloys based on Ni 3 Al showed superior resistance to cavitation erosion compared to many commercial alloys. The volume loss rate from erosion is reduced with fine-grained material as compared to that with coarse-grained specimens. Cold working was also found to reduce the loss due to cavitation erosion. In their review article, Karimi and Martin[33] have examined cavitation erosion and the response of materials to cavitation. High-speed photography has been used to verify liquid jet formation and shock wave generation from the collapse of single bubbles. The shock wave is usually emitted by a high velocity jet of liwuid created during the end of the collapse period. The effective damage of collapse depends strongly on the distance between collapse center and the solid surface. Accumulation of shock-fatigue damage deforms the surface layers 69 plastically and produces cracks leading to material removal. Chemical and electrochemical factors can enhance the damage. Karunamurthy et al.[34] have experimentally investigated the mechanism of mechanical degradation during cavitation erosion in silicon nitride. Multiple intergranular and transgranular fractures were noted as the initial stages of erosion. Dislodging of grains from the surface led to the formation of small scale pits, resulting in a rough surface. The rough surface provided favorable condition for bubble nucleation, which accelerated the damage. Crack initiation and propagation were also noted showing that cavitation erosion is a surface fatigue process. The rate of erosion was found to largely depend on the microstructure of the material. In their review of ultrasonic cavitation at solid surfaces, Schuklin et al.[47] note that the effects of ultrasound derive primarily from cavitation, where bubble collapse results in an enormous concentration of energy from the conversion of the surface energy, kinetic energy of liquid motion into heat and chemical energy. The collapse creates locally extreme pressures and temperatures as well as shock waves and a liquid impinging on the surface. The difference between cavitation bubble at a surface and in bulk is that in contact with a solid also the surface energy becomes relevant and the liquid flow is asymmetric. The bubbles evolve as hemispherical caps during most of the lifetime of the doublet, except during the collapse phase where the bubbles tend to form jets towards each other and parallel to the wall. Wall-normal jets are well known for the collapse of a bubble close to a solid boundary which breaks the spherical symmetry of the inner flow during the collapse: one side of the bubble accelerates inward more rapidly than the opposite side leading to the formation of the re-entrant jets. Here the symmetry breaking is due to the neighboring bubble and the jets are wall parallel. Only those bubbles within 70 a distance from the solid equal to the bubble radius can produce any damage, and that if these collapse asymmetrically, cavitation must result from a combination of the impact by microjets and the shock wave generated by the bubble collapse. Power ultrasound was used to investigate physical and chemical effects of acoustic cavitation at the water-glass interface. Bubble collapse at the water-glass interface initiates not only mechanical erosion but also accelerates the (selective) leaching processes of the glass components. Erosion of the surface increases with sonication time and applied acoustic intensity for soda lime glass as well as for fused silica. Primary areas of erosion (cavities and microcracks) act by them- selves as nucleation sites, facilitating the formation of bubbles, and thus promoting the erosion of the surface. The cavitation erosion resistance of glass under acoustic cavitation is governed by the brittleness index. In other words, the harder is the glass, the less resistant under acoustic cavitation it is. Virot et al.[48] have used power ultrasound (20 kHz) to investigate acoustic cavitation at the water-glass interface. Physical and chemical effects resulting from ultrasound-induced erosion on glass surfaces were found to be highly dependent on the experimental settings such as distance, acoustic intensity, sonication time, hardness of surfaces, etc. Chemical analysis of water indicated that the bubble collapse at the water-glass interface accelerates the leaching process of the glass components, in addition to mechanical erosion. Cavitation erosion has also been observed in geological systems. Momber[120] has studied short-time cavitation erosion of rocks and cementitious composites in a laboratory flow cavitation chamber. It is reported that erosion mainly occurs in intergranular or intercrystalline mode and fracture toughness of rock and conventional 71 concrete is a good indication of cavitation erosion resistance. The author notes that materials with low degree of pre-existing flaws, e.g. glass and high-strength silica cement, show extraordinarily high erosion resistance. In another study, Momber[121] used Scanning Electron Microscope (SEM) to study the behavior of cement paste, mortar and concrete during erosion by high-speed water jets with speeds upto 470 m/s. Crack bridging, microcracking, crack branching and crack blunting were identified, thus verifying quasi-brittle fracture behavior of the materials during erosion. Cavitation erosion is also widely reported and studied in the cooling systems of power plants. Nagaya et al.[36] present a method of detecting cavitation externally, using accelerometers and microphones. Local rise in flow velocity due to decrease in area of the orifice leads to a drop in fluid pressure. An increase in flow area leads to decrease of fluid velocity, thus increasing fluid pressure and thereby collapsing cavitation bubbles generated earlier. Such a collapse generates impact pressure, inducing vibration, which can lead to erosion of pipes, thus damaging the piping system of power plants. To detect cavitation bubbles the authors used an accelerometer in high frequency range as well as directional and non-directional microphones. Aluminum was used to test erosion, as it is softer than different types of steel commonly used in pipes, hence can be used to evaluate the threshold of erosion. The authors have defined Cavitation number, σ as σ = 2(P – P v )/ ρV o 2 , where P = downstream pressure of the orifice [Pa], P v = saturated vapor pressure at experimental temperature [Ps], ρ = density of water (kg/m 3 ) and V o = average velocity in the orifice (m/s). In the experiments, the water flow velocity was between 15.0-15.4 m/s and the temperature of water was kept at 20°C. The authors observed that the number 72 of cavitation bubbles increased with decreasing cavitation number. Erosion pits were observed below σ = 0.8 at 50 and 75 mm from the orifice. At a distance of 100 mm from the orifice, erosion pits were only at σ = 0.7. The authors concluded that the cavitation number at incipient cavitation erosion was between 0.8 and 0.9. No cavitation erosion was observed above σ = 0.9. Lee et al.[43] investigated the cavitation erosion behavior of hardfacing alloys such as Co-based Stellite 6 and an alternate new Fe-based alloy. Cavitation erosion tests were performed using vibratory cavitation erosion testing machine. The test specimen was attached to the end of an ultrasonic horn and its weight loss was measured in a distilled water bath. For Co-based Stellite 6, the resistance to cavitation erosion was attributed to strain-induced phase transformation as well as continuous fine-scale twinning providing an efficient way to accommodate strain energy and strain-relaxation in blocked twins. The crack initiated in the carbide-matrix interface and was found to propagate mainly in the carbide region. In case of Fe-based alloy with austenitic structure, the low stacking fault energy prevents recombination of partial dislocations, which are necessary for cross-slip. Hence, planar slip is the primary mode of deformation in the initial stages of cavitation erosion. The authors conclude that the new Fe-based alloy was found to have similar cavitation erosion resistance as Co-based Stellite 6, which was attributed to the hardened matrix, thereby suppressing the propagation of cracks initiated at the carbide-matrix interface. Cavitation erosion damage is also being investigated in case of spallation neutron sources. Rapidly deposited heat energy from the proton beams generates pressure waves in mercury, resulting in the formation of cavitation bubbles. These bubbles, upon 73 collapse, can cause severe damage to the vessel wall. Haines et al.[40] have investigated cavitation erosion in Spallation Neutron Source (SNS) mercury targets. The authors performed a series of in-beam tests under varying conditions to quantify pitting damage. The results indicate that the damage is sensitive to beam intensity, surface treatment and gas injection. The authors found that reducing the beam intensity by a factor of two reduced the damage by an order of magnitude while the injection of a small amount of helium gas bubbles reduced the erosion by a factor of approximately four. Ida et al.[41] have carried out numerical studies to evaluate the effect of gas bubble injection into a liquid on the suppression of cavitation inception. It was found that the injected gas bubbles decrease the magnitude of the negative pressure in liquid mercury by radiating a positive-pressure wave in its early stage of expansion, which results in the suppression of the explosive expansion of cavitation near the injected bubbles. In their study of free gas bubbles in water, Ohl and Ikink[49] observed jet formation during the collapse of micron-size bubbles initiated by a shock front. The direction of the jet depended on the direction of wave propagation. The length of the jet tip increased with increasing bubble radius and the average jet speed increased from 20 to 150 m/s for bubbles ranging between 7 and 55 µm. Kodama and Tomita[122] studied the collapse of a single cavitation bubble near a gelatin surface using high-speed framing photography. The cavitation bubbles near the gelatin surface did not produce significant liquid jets directed at the surface, and tended to migrate away from the surface. The interaction of an air bubble with a shock wave yielded a liquid jet inside the bubble, penetrating the gelatin surface. The impact pressure was of the order of 10 MPa and the penetration 74 depth was ~ 1 mm and the speed of the jet was of the order of 10 m/s. Delius et al.[50] have carried out experiments to explore the effect of extracorporeal shock waves on hemoglobin release from red blood cells. It was found that shock waves applied under low static excess pressures led to considerably lower hemolysis. The authors identify shock wave-gas bubble interaction as the underlying mechanism of this effect. It is suggested that the size of the observable gas bubbles in the focal area after the passage of the shock wave decreases more at low repetition frequency and higher excess pressures reduce the effect of the subsequent shock waves. In their review of shock-bubble interactions (SBIs), Ranjan et al.[123] note that the most important aspect of SBIs is the vorticity deposition due to the misalignment of the pressure and the density gradients. As this field of shock waves (including the primary incident shock wave, along with the secondary refracted, reflected, diffracted, and focused waves) passes over the bubble, vorticity is produced in the flow. The interaction of a planar shock wave with a cylindrical or spherical density inhomogeneity produces three fundamental processes (compression/acceleration, shock reflection/refraction, and baroclinic vorticity generation), which are actually tightly coupled together. Their effects include the initial distortion of the bubble and the formation of primary and secondary vortex pairs or vortex rings. Johnsen and Colonius[124] have used a high-order accurate shock- and interface- capturing scheme to simulate the collapse of a gas bubble in water. In order to better understand the damage caused by collapsing bubbles, the dynamics of the shock-induced and Rayleigh collapse of a bubble near a planar rigid surface and in a free field are 75 analysed. Collapse times, bubble displacements, interfacial velocities and surface pressures are quantified as a function of the pressure ratio driving the collapse and of the initial bubble stand-off distance from the wall; these quantities are compared to the available theory and experiments and show good agreement with the data for both the bubble dynamics and the propagation of the shock emitted upon the collapse. Non- spherical collapse involves the formation of a re-entrant jet directed towards the wall or in the direction of propagation of the incoming shock. In shock-induced collapse, very high jet velocities can be achieved, and the finite time for shock propagation through the bubble may be non-negligible compared to the collapse time for the pressure ratios of interest. Several types of shock waves are generated during the collapse, including precursor and water-hammer shocks that arise from the re-entrant jet formation and its impact upon the distal side of the bubble, respectively. The water-hammer shock can generate very high pressures on the wall, far exceeding those from the incident shock. The potential damage to the neighbouring surface is quantified by measuring the wall pressure. The range of stand-off distances and the surface area for which amplification of the incident shock due to bubble collapse occurs is determined. While cavitation bubble collapse has interested researchers for a while now, there are a lot of unanswered questions: What are the characteristics of microturbulence generated during bubble collapse? Is there is a limiting factor to the damage due to nanojets? What is the extent of chemical damage? Are there chemical reactions leading to the mechanical damage? How does the extent of damage scale with size of the bubble? To have a comprehensive understanding of cavitation erosion damage, it is important to investigate 76 the factors affecting the chemical damage on a surface in addition to mechanical damage during cavitation erosion. Since silica has been an archetypal material for studying stress corrosion cracking, shock, indentation, and behavior under high-pressure, cavitation erosion damage at the silica-water interface has been a topic of recent experimental study. We thus use the silica-water system to study the factors affecting mechanical as well as chemical damage due to shock-induced collapse of nanobubble. 4.2 Simulation We examine above questions through billion-atom reactive molecular dynamics (MD) simulations of shock-induced collapse of a nanobubble near a silica slab. The simulations were performed on the full BlueGene/P supercomputer (163,840 processors) at the Argonne Leadership Computing Facility. In the billion-atom MD simulations, the system consists of a 60 nm thick silica slab placed in an MD box of dimensions 285 × 200 × 200 nm 3 in the x, y and z directions respectively; see Figure 4.1. The rest of the MD box contains water molecules and a spherical nanobubble of diameter 97.6 nm. The initial distance between the center of the bubble and the silica surface is 90 nm, and the stand- off parameter, s d = 1.84. We applied shock with a particle velocity of 2 km/s and found the shock velocity in water to be 5 km/s. This is in good agreement with experimental data of Rybakov and the MD simulation results of Vedadi et al.[125] We used periodic boundary conditions normal to the direction of shock-wave propagation (-x). The equations of motion were integrated with the velocity-Verlet algorithm using a time step of 0.5 fs. 77 System 10 8 atoms 10 9 atoms Thickness of silica slab 40 nm 60 nm Diameter of bubble 40 nm 97.6 nm #silica atoms 16.9 x 10 6 158.6 x 10 6 #water molecules 29.5 x 10 6 280.6 x 10 6 #gas particles 1.2 x 10 6 17.7 x 10 6 MD box 180 x 80 x 80 nm 3 285 x 200 x 200 nm 3 #processors 2880 (1 x 48 x 60) 163840 (2 x 320 x 256) Resources HPCC, USC BlueGene/P, ALCF Table 4.1 Table listing the details of the simulations done Figure 4.1: Schematic of the system used for simulation 78 4.2.1 Preparation of Amorphous Silica Starting with the β-cristobalite crystalline configuration at a density of 2.2 g/cc, we melt the system by scaling atomic velocities over a time period of 30 ps until the temperature of the system reaches 4000 K. After thermalizing for 20 ps at 4000 K, we cool the system to 1900 K over 10 ps and thermalize it for 10 ps. Repeating this cooling and relaxation cycle over a time period of 120 ps, the temperature of the system is reduced to 1400 K, 950 K, 600 K, 375 K and finally to 300 K. Figure 4.2 shows the temperature profile during the melt-quench schedule. The silica system at room temperature is amorphous, and its structure agrees very well with neutron scattering measurements[126]. To create a silica slab for shock simulations, we remove periodic boundary conditions and passivate the dangling bonds with silanols. The silanol density on the slab surface (4.6 nm -2 ) is in good agreement with experiments[127]. Figure 4.2: Melt-Quench procedure for preparing amorphous silica 79 4.2.1 Interatomic Potential We have designed and implemented a multi-body potential[128], called Environment Dependent Reactive (EDR) Potential, governing the interactions between silica and water atoms. We use FORTRAN 77 to program the simulation and FORTRAN90 for various analysis tools. Parallelization of the codes is done using Spatial Decomposition Method, for distributing the system among multiple processors and Message Passing Interface (MPI) for inter-processor communication. The potential form used for both silica and water consists of a two-body term and a three- body term to incorporate various effects, including steric repulsion, Coulomb interaction, charge-dipole interaction, van der Waals interaction, and covalent bond bending and stretching. The total energy of the system is written as , (4.1) where with being the position of the i-th atom, and is the distance between atom i and j. In Eq. (4.1), the first, two-body, term sums over neighbor atoms, j, within a cutoff range r c from atom i, and has the form, , (4.2) whereas the second, three-body, term, sums over atoms, j and k, that are covalently bonded with atom i within a distance r 0 , and the form of this three-body interaction is: . (4.3) In Eq. (4.2), is the steric repulsion strength, is the effective charge in units of the electronic charge, represents the charge-dipole strength, is the strength of van der € E tot = V ij (2) (r ij ) i< j ∑ + V ijk (3) (  r ij ,  r ik ) i< j<k ∑ €  r ij =  r i −  r j €  r i € r ij =|  r ij | € V ij (2) (r) = H ij r η ij + Z i Z j r e −r r 1s − D ij 2r 4 e −r r 4s − w ij r 6 € V jik (3) (  r ij ,  r ik ) = B ijk exp ξ r ij − r 0 + ξ r ik − r 0 $ % & & ' ( ) ) (cosθ ijk −cosθ 0 ) 2 1+ C ijk (cosθ ijk −cosθ 0 ) 2 r ij ,r ik ≤ r 0 ( ) € H ij € Z i € D ij € w ij 80 Waals attractions, is the exponent of the steric repulsion, and and are the screening lengths for the Coulomb and charge-dipole interactions, respectively. In Eq. (4.3), is the strength of the three-body interaction, and are constants, and is the angle formed by and . For computational efficiency, the 2-body pair potential is truncated at r c and shifted for r < r c , so that the potential and its first derivative are continuous at r c : , (4.4) The potential parameters for silica have been determined to fit experimental values of bulk modulus, equilibrium volume, and melting temperature of silica. Potential parameters for pure water have been determined to fit the force constants of symmetric bond stretching and bending of a single water molecule. The parameters for the silica and water potentials are given in Tables 4.2 and 4.3, respectively. Si O Z i (e) 1.2 -0.6 r 1s = 4.43 Å r 4s = 2.5 Å r c = 5.5 Å e = 1.602×10 -19 C Two body Si-Si Si-O O-O η ij 11 9 7 H ij (eV Å η ) 0.39246 78.3143 355.5263 D ij (eV Å 4 ) 0 3.456 1.728 W ij (eV Å 6 ) 0 0 0 Three body € η ij € r 1s € r 4s € B ijk € ξ € C ijk € θ ijk €  r ij €  r ik € V ij (2,shifted) (r) = V ij (2) (r)−V ij (2) (r c )−(r−r c ) dV ij (2) /dr [ ] r=r c 0 (r < r c ) (r≥r c ) $ % & ' & 81 B ijk (eV) θ o (deg) C ijk ξ (Å) r 0 (Å) O-Si-O 4.993 109.47 0 1.0 2.60 Si-O-Si 19.972 141.0 0 1.0 2.60 TABLE 4.2. Parameters for two- and three-body part of our potential used for MD simulation of silica. H O Z i (e) 0.32983 -0.66966 r 1s = 4.43 Å r c = 5.08 Å e = 1.602×10 -19 C Two body H-H H-O O-O r 4s (Å) 2.5 1.5113 2.5 η ij 9 9 9 H ij (eV Å η ) 0 0.61437 1965.88 D ij (eV Å 4 ) 0 0.2611 2.0887 W ij (eV Å 6 ) 0 0 10.0 Three body B ijk (eV) θ o (deg) C ijk ξ (Å) r 0 (Å) H-O-H 52.9333 97.9476 0 0.75 1.4 O-H-O 0 0 0 0 0 TABLE 4.3. Parameters for two- and three-body part of our potential used for MD simulation of water. Using the above two sets of parameters developed for pure silica and pure water using the same potential form, we now combine the two potentials to study the reaction between silica and water. The interaction between silicon and oxygen from water is taken as the 82 same as that between silicon and oxygen from silica. Likewise, the interaction between hydrogen and oxygen from silica is taken as the same as that in water. For the Si-H pair, we only account for the charge-charge interaction for simplicity. Screening length of the Si-H charge-charge interaction is taken to be the same as that of H-H charge-charge interaction. Additional three-body interaction is added for Si-O-H to describe the bond angles of silanol groups on silica surface, and to reproduce the correct interaction energy for reactions between silica and water. The three body bond cutoffs for Si-O and O-H pairs are taken to be the same as those in the corresponding silica and water 3-body potentials. The bond strength and equilibrium bond angle for Si-O-H interaction are taken to be 36.45265 eV and 110.0 degree, respectively.[128] The interaction between a pair of oxygen atoms is taken as an interpolation between the O-O interaction in silica and that in water by incorporating the many-body bond-order effect, i.e., the number of Si and H atoms surrounding the two oxygen atoms. This interpolation scheme ensures a smooth transition of the O-O interaction from silica potential to the water potential, when the environment of the oxygen changes from pure silica to pure water. In general, our interpolation scheme represents the interatomic potential energy E tot as a superposition of multiple effective force-field (EFF) models , (4.5) where and are the two- and three-body interatomic potential terms of the α- th EFF model, and the support functions satisfy the normalization conditions, . (4.6) € E tot = α ∑ f ij (2)α V ij (2)α i, j ∑ + α ∑ f ijk (3)α V ijk (3)α i, j.k ∑ € V i (2)α € V ij (3)α € f ij (2)α α ∑ = f ijk (3)α α ∑ =1 83 The interpolation scheme is formally identical to the divide-and-conquer (DC) approac to QM calculations based on DFT at both density and density-matrix levels, the partition-of- unity approach to the finite element (FE) method in continuum mechanics, and the graded hybridization scheme in the multiscale MD/FE simulation method. This scheme allows the reuse of the well-established silica potential, and instead to focus on the development of a water potential for describing water properties. Another advantage is that the program does not need to differentiate whether the oxygen belongs to water or silica to decide which parameters to apply, leading to an efficient implementation of the MD program. The specific expression for the two-body interaction between a pair of oxygen atoms is given by , (4.7) where and are the oxygen-oxygen two-body potential in pure silica and water, respectively. The weight factor is a smooth function of the number of Si and H atoms surrounding the two oxygen atoms and defined as , (4.8) where is the number of β-type atoms (β = silicon or hydrogen) around the i-th oxygen atom. Special care needs to be taken to handle zero number of Si and H atoms near both O atoms, i.e., in the case of oxygen gas, where an oxygen-oxygen interaction in the oxygen gas environment is needed. Since oxygen gas in the system we are interested in is not a major species, we model it using the O-O interaction in silica. The number is given as the sum of function that describes the fraction of the β-type atoms belonging to the i-th oxygen, € V ij (2) (r ij ) =V ij (2)silica (r ij )(1− f ij ) +V ij (2)water (r ij ) f ij € V ij (2)silica € V ij (2)water € f i, j € f ij (n i Si ,n i H ,n j Si ,n j H ) = n i H + n j H n i Si + n i H + n j Si + n j H € n i β € n i β € Θ 84 , (4.9) where the continuous function varies smoothly from 0 to 1 depending on the distance between the i-th oxygen and its k-th β-type neighbor: (4.10) A plot of is shown in Fig. 4. Values of parameters and are 2.0 and 0.3 Å, or 1.4 and 0.3 Å, for or , respectively. Fig. 4.3: The number of neighboring silicon/hydrogen atoms around an oxygen atom is a function of the distance between the Si/H and the O atoms. (source: Ref #[128]) The interatomic potential for silica was validated by neutron scattering data for static structure factor and T(r) (=r 2 g(r), g(r) is the radial distribution function).[126, 129, 130] Our molecular dynamics (MD) results for elastic moduli (Young’s modulus = 66.9 GPa, Bulk modulus = 39.2 GPa, Poisson’s ratio = 0.22) are in good agreement with experiments.[131] The MD result for fracture toughness, K 1C = 1 MPa•m 1/2 , is within the € n i β = Θ k∈β ∑ (r ik ) € Θ € Θ(r) = 1 (r≤R β −D β ) 1− r−R β +D β 2D β + sin(π(r−R β +D β ) D β ) 2π (R β −D β ≤r <R β +D β ) 0 (R β +D β <r) ' ( ) ) * ) ) € Θ(r) € R β € D β € β =Si € β =H 85 range of experimental values (0.8-1.2 MPa•m 1/2 ).[131] Recently we performed nanoindentation simulations and calculated the hardness of amorphous silica to be 10.6 GPa,[132] which is close to the experimental value of 10 GPa.[133] For water under ambient conditions, we find the O-H bond length is 0.97 Å and H-O-H bond angle to be 104°. These results are in good agreement with X-ray diffraction experiments.[134] Our MD simulation for water-silica interface shows that the silanol concentration lies between 4.3 and 5.4 nm -2 above 300 K, which is in good agreement with experimental results.[127] The calculated hydrolysis reaction energy between water and silica varies between 0.34 and 0.70 eV, depending on the reaction sites, which falls within the range of values obtained by ab initio calculations.[135] Also, the MD simulation results for the lowest energy structure of the gas phase orthosilicic acid are in good agreement with ab initio calculations.[136] The heat of immersion of a hydrated silica slab in water is calculated to be 0.18 J/m -2 , which is also in good agreement with experimental values.[137] 4.3 Collapse of empty nanobubble Figure 4.4 shows snapshots of the system during the collapse of empty nanobubble of initial diameter 97.6 nm. The system is given initial particle velocity of 2 km/s and the shock wave is generated using a momentum mirror. The shock wave arrives from the right (blue region in Figure 4.4a). The density of water in the shock region is 1.5 g/cc and the pressure is 10 GPa. The MD simulation reveals that the structure of water has changed significantly in the shock region, e.g., on average, the number of nearest 86 neighbors of a water molecule has increased from 4 to 8. The bubble almost completely collapses upon interacting with the shock wave by t = 22 ps (Figure 4.4b). A high-speed water jet is formed during the collapse of the empty nanobubble which proceeds to hit the slab at t = 35 ps (Figure 4.4c). A cavitation pit is formed on the silica slab as a result of erosion caused by the high-speed water nanojet (Figure 4.4d). The shock wave crosses the silica slab at t = 50 ps, at which point the simulation was stopped. Figure 4.4: Snapshots of the system during collapse of empty nanobubble at (a) t = 5 ps, (b) t = 22 ps, (c) t = 35 ps, and (d) t = 50 ps. Silica slab is shown in yellow and orange. A cross section of water across the system is shown. Normal density water is shown in green while water under shock is shown in blue. 87 4.3.1 High-speed water nanojet Figure 4.5: Position of the tip of water nanojet during the simulation. The position of the shock front away from the bubble is shown in grey. Figure 4.5 shows the progression of water nanojet during the collapse. The shock front hits the bubble at t = 10 ps. As the bubble collapses, a high-speed water nanojet forms, which hits the distal side of the bubble at t = 20 ps. The jet continues to accelerate till it hits the silica slab at t = 35 ps. Following the impact, the progression of the nanojet is suppressed at the erosion of the silica slab begins. When the shock front hits the nanobubble, water molecules at the periphery of the bubble rush in to form a focused nanojet in the direction of the shock-wave propagation. Figure 4.6a shows the nanojet velocity streamlines shortly before the bubble collapses completely. The nanojet reaches a speed of 7.8 km/s, which is significantly higher than the shock velocity (5 km/s). 88 Experimentally, Ohl and Ikink[49] have observed similar jetting phenomenon in the collapse of micron-size bubbles. They find the microjet length is approximately three times the initial bubble radius. Kodama and Tomita[122] find a similar scaling relation for bubbles of radii between 0.1 and 1 mm. Together with our previous simulation results[51, 125] on nanometer size bubbles, these experiments suggest that the scaling relation between the jet length and bubble radius holds for a wide range of bubble sizes. Chemical activity in water increases significantly with the onset of bubble collapse and nanojet formation. We observe an increase in the hydronium ion production with the nanojet growth. The blue curve in Figure 4.6b shows the spatial distribution of H 3 O + ions along the direction of nanojet propagation (i.e., along –x). The H 3 O + ions are counted in voxels of dimensions 1× 5 × 5 nm3 around the nanojet. At t = 25 ps, the H 3 O + population peaks at the location of the shock front inside the nanobubble. We observe an increase in the number of H 3 O + ions when the nanojet hits the distal side of the bubble at 202 nm. The peak in the red curve at 218 nm is due to a secondary shock wave generated by the nanojet impact on the distal side of the bubble. The secondary shock propagates in the direction opposite to the primary shock wave. Thus we observe significant chemical activity in the nanojet generated due to the shock-induced collapse of an empty nanobubble. The chemical composition of water changes with the pressure in the jet, with high-pressure regions resulting in increased concentration of H 3 O + ions. 89 Figure 4.6: (a) Water nanojet formed during shock-induced collapse of the nanobubble. The silica slab is shown in yellow. The central slice across the silica slab is color coded by the magnitude of the velocity of water molecules. The velocity field lines show a focused jet of water molecules. (b) Number of H 3 O + species along the central slice of the system. H 3 O + ion production increases significantly when the nanojet hits the distal side of the bubble at time t = 25 ps (blue), and t = 5 ps (red). 90 4.3.2 Pressure profile in water during collapse Figure 4.7: Pressure profile in water at (a) t = 15 ps, (b) t = 25 ps, (c) t = 30 ps, (d) t = 46 ps, (e) t = 48 ps, and (f) t = 50 ps We calculate pressure in the system using virial method. The entire system is divided into 1 nm thick cubic voxels. The average pressure in each voxel is then calculated from the contribution of each atom. For each voxel, a three-layer average is done to reduce the noise in the calculation. Figure 4.7 shows the spatial distributions of the pressure in the water at time t = 15, 25, 30, 46, 48, and 50 ps. The shock front hits the nanobubble at t = 10 ps, following which the nanojet begins to form (Figure 4.7a). The nanojet hits the distal side of the bubble at t = 20 ps leading to the formation of a region of high pressure 91 at the tip of the nanojet (Figure 4.7b). As the nanojet proceeds in the direction of the shock wave, a secondary shock wave is observed in the opposite direction (Figure 4.7c). The forward motion of the nanojet is halted when it hits the silica slab at t = 35 ps. The nanojet then begins to erode the silica surface, moving laterally along the slab. Figure 4.7d shows the reflected wave from the silica slab. By t = 48 ps (Figure 4.7e), the reflected wave splits from the forward wave and propagates along the direction opposite to that of the initial shock wave. The peak pressure (18 GPa) at t = 50 ps (Figure 4.7f) agrees very well with the estimate of water-hammer pressure, P = ρcv, where ρ is the mass density, c is the speed of sound and v is the jet speed.[53] Taking the sound velocity in water to be 1.5 km/s, and using the values ρ = 1.5 g/cc and v ≈ 7.8 km/s from our simulations, we estimate water-hammer pressure to be ~18 GPa. The pressure profiles reveal the dynamic nature of the shock-induced collapse of an empty nanobubble. The increased pressure in the high-speed nanojet results in increased concentration of hydronium ions (Figure 4.6b), thereby altering the chemical composition of water. 92 4.3.3 Cavitation damage on the silica slab Figure 4.8: (a) Cavitation damage due to water nanojet. (b) Side-view of the cross-section of the system at t = 50 ps. Silica slab is shown in yellow while the water is shown in cyan (hydrogen) and blue (oxygen). 93 The nanojet impact on the silica slab initiates pit damage at t = 35 ps (Figure 4.8a). At 50 ps the shock wave has already crossed the silica slab, leaving behind a pit of diameter 78 nm and depth 17.5 nm. The cavitation pit is spherical in shape as evident in Figure 4.8b, which shows the side view of the cross-section of the cavitation damage on the silica slab at t = 50 ps. Almost 52% of the silica surface is damaged due to erosion caused by the high-speed nanojet generated during the shock-induced collapse of empty nanobubble. Experimental results of Kodama and Takayama[138] indicate that the ratio between the size of the damaged pit on a gelatin surface to bubble radius is close to 0.1, whereas we find this ratio to be 0.3. Virot et al.[48] have observed similar erosion damage on fused silica and soda lime glass due to acoustic cavitation at the water-silica interface. Figure 4.9: Maximum depth of cavitation pit with time during collapse of bubble of initial diameter of 40 nm (red) and 97.6 nm (blue). 94 Figure 4.9 shows the maximum depth of the cavitation pit during the erosion process. The red curve represents the system with an initial bubble diameter of 40 nm while the blue curve represents that with an initial bubble diameter of 97.6 nm. The erosion process occurs at a constant rate as evident from the uniform slope of the cavitation depth vs time in both the systems. Figure 4.10: Volume of the cavitation pit in systems with nanobubbles of diameter 40 nm (red) and 97.6 nm (blue). To quantify mechanical damage in the system, we calculated the volume of the cavitation pit in the silica slab. For comparison, we performed another MD simulation on a system with 100 million atoms. The initial diameter of the bubble was 40 nm and the stand-off parameter was the same as in the case of the billion-atom MD simulation. Figure 4.10 shows the pits in silica resulting from nanojet impact in the billion-atom and 100 million- atom systems. These results indicate that the pit volume scales linearly with the bubble volume. 95 Figure 4.11: (a) Snapshot shows atoms at the pit edge of the silica slab in the system with 10 9 atoms. Silicon-oxygen bonds are yellow-red, and cyan spheres represent hydrogen atoms inside the silica slab. Blue spheres represent oxygen atoms in water molecules within 0.5 nm from the silica slab. (b) Silanol distribution in the silica slab as a function of the distance from the pit center for bubbles of diameters 40 nm (red) and 97.6 nm (blue). 96 Figure 4.11a shows an atomistic view of a 0.5 nm thick slice at the center of the damaged silica network after the nanojet impact at time t = 50 ps. We observe a large number of hydrogen ions forming bonds with Si and O in the damaged region of the silica slab. To quantify chemical damage, we show in Figure 4.11b the distribution of SiOH species in 1 nm deep annular disks centered at the damage pit. The radial distance between successive annular disks is 1 nm. The red and blue curves show SiOH distributions in the 100 million-atom and 1 billion-atom systems. Initially the number of SiOH groups on the surface was ~6 nm -2 . After the shock wave crosses the silica slab we observe a considerable increase in the population of SiOH, which peaks around the pit edge. Experimentally, Virot et al.[48] have observed chemical leeching in silica glass due to cavitation bubbles. The shape of the surface changes from spherical to planar at the edge of the pit resulting in discontinuity in the curve, marked by the circled regions in Figure 4.11b. 4.3.4 Pressure profile in silica slab As mentioned in section 4.3.2, we calculate pressure in the system using virial method. The entire system is divided into 1 nm thick cubic voxels. The average pressure in each voxel is then calculated from the contribution of each atom. For each voxel, a three-layer average is done to reduce the noise in the calculation. Figure 4.12 shows the pressure distribution across the silica slab during the erosion process. The nanojet hits the slab at t = 30 ps, initiating cavitation erosion. The pressure decreases radially outward from the center of the cross-section, which is the point of impact with the nanojet. The peak pressure obtained in the silica slab is 18 GPa. 97 Figure 4.12: Pressure profile across the silica slab during collapse of empty nanobubble at (a) t = 35 ps, (b) t = 40 ps, (c) t = 45 ps, and (d) t = 50 ps 98 4.4 Collapse of gas-filled nanobubble We also performed two sets of simulations in which the nanobubble contained an inert gas at a density of 0.08 g/cc. The gas-gas interatomic interaction was modeled using 12-6 Lennard-Jones Potential while the gas-water interaction had just the repulsive part of the gas-gas interaction. The systems contained 10 9 and 10 8 atoms and the initial bubble radii in these systems were again 97.6 nm and 40 nm, respectively. The stand-off parameter and the particle velocity were the same as before (s d = 1.84, v p = 2 km/s). Figure 4.13: Snapshots of the system during collapse of gas-filled nanobubble at (a) t = 5 ps, (b) t = 25 ps, (c) t = 35 ps, and (d) t = 50 ps. Silica slab is shown in yellow and orange. A cross section of water across the system is shown. Normal density water is shown in green while water under shock is shown in blue. The gas particles are shown in magenta. Figure 4.13 shows snapshots of the system during the collapse of gas-filled nanobubble of initial diameter 97.6 nm. The system is given initial particle velocity of 2 km/s and the shock wave is generated using a momentum mirror. The shock wave arrives from the 99 right (blue region in Figure 4.13a). The density of water in the shock region is 1.5 g/cc and the pressure is 10 GPa. The gas-filled bubble absorbs the shock wave (Figure 4.13b) and does not allow the formation of a high-speed nanojet. During the collapse, the bubble undergoes significant deformation (Figure 4.13c). By 50 ps, the bubble deforms into a torroid-shaped aggregate (Figure 4.13d) 4.4.1 Lack of high-speed nanojet Figure 4.14: Comparison of jet position during shock-induced collapse of empty and gas- filled nanobubble. Figure 4.14 shows the position of the tip of the jet formed during the collapse as a function of time. The red curve represents the collapse of empty nanobubble and the blue 100 curve represents the gas-filled nanobubble. It is evident that the gas particles in the bubble present significant resistance to the formation of a high-speed nanojet. The speed of the jet formed is only 2.4 km/s, as compared to 7.8 km/s in case of an empty nanobubble. The shock front, with a speed of 5 km/s, travels faster that the jet formed during the collapse. As will be shown later, the absence of high-speed nanojet results in significantly less damage on the silica slab. Haines et al.[40] have observed similar effect in case of mercury targets in Spallation Neutron Sources. They found that the presence of helium bubbles in the liquid mercury reduced the cavitation pitting damage by four times. Our simulations show that the effect translates very well to water and silica system. Figure 4.15: Vortex rings around the deformed gas bubble during shock-induced collapse. The inset shows the direction of local velocity vectors (black arrows). The color gradient represents the magnitude of the local velocity. 101 Following the collapse, the gas-filled bubble deforms into a toroid-shaped aggregate, as shown in Figure 4.15. Similar toroid-shaped deformation has also been seen in shock- wave experiments on spherical and cylindrical inhomogeneity[123], and in continuum simulations of shock-induced collapse of a bubble near a wall[124]. We find that the nanobubble shape change is due to vortices at the periphery of the gas bubble, shown in Figure 4.15 by the green lines. Here the velocity field is represented by streamlines around the deformed gas bubble and color coded by magnitude. The streamlines represent the path that a rigid particle would take if left at the origin of the streamline. These vortices are generated by the gradients in pressure and density at the periphery of the collapsing nanobubble. Experiments also reveal the formation of vortex rings from pressure and density gradients generated by the collapse of gas-filled bubbles in three dimensions.[123] The inset in Figure 4.15 shows the velocity vectors in the region of vortex rings. Since the density of the nanobubble is less than that of the medium, the vortex rings at the bubble periphery rotate towards the center of the bubble. This has also been observed in continuum simulations for bubble density less than that of the medium. These vortex rings reveal the highly dynamic flow pattern around the gas-filled nanobubble during shock-induced collapse. We are also able to identify hot spots inside the bubble during the collapse. The high-pressure regions inside the bubble coincide with the centers of vortex rings around the periphery of the bubble. Knowledge about such flow patterns can be useful in designing applications involving interaction of shock wave with spherical or cylindrical inhomogeneity. 102 Figure 4.16: H 3 O + concentration in water around the periphery of the collapsing gas- filled nanobubble To quantify chemical activity during the collapse, we divide the system into 1 nm 3 voxels and calculate the average number of H 3 O + ions in each voxel, as shown in Figure 4.16. We observe that the hydronium ion production is highly concentrated in the central region of the collapsed bubble and around the vortex rings. As the shock wave hits the gas-filled bubble, it deforms into a torroid-shaped aggregate and vortex rings are formed around the bubble. This highly dynamic flow pattern results in increased chemical activity in regions around the water-gas interface. It should be noted that the region beyond the reflected wave (yellow-green) has low concentration of H 3 O + ions, indicating that the concentration of hydronium ions depends on the pressure in the water. 103 4.4.2 Reduced mechanical damage Figure 4.17: Side-view of the cross-section of the system with (a) empty and (b) gas- filled nanobubble at t = 50 ps. Silica slab is shown in yellow while the water is shown in cyan (hydrogen) and blue (oxygen). Gas particles are shown in red. Figure 4.17 shows the side-view of the systems with empty and gas-filled bubbles at t = 50 ps. The silica slab is shown in yellow while water is shown in blue and cyan. The red region represents the gas-filled bubble. In the case of a gas-filled bubble, the mechanical damage is significantly less compared to the damage caused by the collapse of an empty bubble. The presence of gas particles in the nanobubble prevents the formation of high- speed nanojet (Figure 4.14). The gas bubble also absorbs much of the energy in the shock wave by deforming into a torroid-shape aggregate (Figure 4.15), thereby protecting the silica slab behind it from cavitation damage. This is consistent with the experimental observation of Ida et al.[41] that gas-filled nanobubbles suppress cavitation erosion damage in liquid mercury. This result can be used constructively as it provides a simple, yet effective method to reduce cavitation erosion damage in highly dynamics applications involving water. 104 Figure 4.18: Comparison of depth of cavitation pit formed due to the collapse of empty and gas-filled nanobubbles of two initial diameters. To quantify the mechanical deformation in the system, we calculate the depth of penetration inside the silica slab. For comparison, we have performed the same calculations for smaller, 100 million atom systems. The 100 million atom systems were simulated for both empty and gas-filled bubble. The initial diameter of the bubble as well as the initial distance of the bubble from the silica slab was 40 nm in the smaller systems, with the other simulation parameters remaining the same. Figure 4.18 shows the depth of the dent, caused by water nanojet, in the four simulated systems. For the same initial bubble diameter, the systems with empty bubble have a considerably deeper dent than the corresponding system with gas bubble. All the four systems simulated in this study have stand-off parameter, s d ~ 2. It is noteworthy that for the same standoff parameter, system 105 with larger initial bubble has a higher depth of dent as compared to the system with smaller initial bubble. The same trend is observed in the systems with gas bubbles. In summary, for the same standoff parameter, the larger bubble will cause more damage than a smaller bubble. However, the presence of inert gas particles significantly reduces the extent of damage on the silica slab. 4.4.3 Pressure profile In Figure 4.19, the silica slab and the gas-filled nanobubble are overlaid with the pressure distribution normal to the silica slab at (a) t = 35 ps, (b) t = 40 ps, (c) t = 45 ps, and (d) t = 50 ps. The pressure profile across the silica slab during the shock-induced collapse of gas-filled nanobubble is uniform across the cross-section. Comparing this with the case of an empty bubble, the pressure in the pit region decreases radially as one moves away from the pit center (Figure 4.12). The peak pressure in the silica slab is 9 GPa, which is half of the peak pressure in case of collapse of empty nanobubble. The uniform profile results from the planar shock wave as the gas particles in the nanobubble prevent the formation of high-speed nanojet. From the contrasting pressure profiles in the silica slab during the collapse of an empty nanobubble and a gas-filled nanobubble, it can be argued that the resultant silica surface is mechanically different and will respond differently to subsequent shock waves. In case of empty nanobubble, the silica slab has a eroded surface under very high pressure, rendering the eroded region more susceptible to further erosion. On the other hand, the collapse of a gas-filled nanobubble leaves a silica surface without any significant erosion, no specific hot spot, and half the peak pressure. 106 Looking at the pressure profiles in the gas-filled nanobubble, we observe that the pressure in the partially collapsed nanobubble is uniform around 7 GPa when the shock front reaches the silica slab at t = 35 ps (Figure 4.1a). The region of increased pressure in the right side of the bubble in Figure 4.19b indicates the presence of the reflected wave from the silica slab. As the shock wave crosses the silica slab, the nanobubble continues to collapse and by t = 50 ps it deforms into a torroid-shape at a uniform pressure of 9 GPa. Figure 4.19: Pressure profiles across the silica slab and inside the gas-filled bubble during shock-induced collapse at (a) t = 35 ps, (b) t = 40 ps, (c) t = 45 ps, and (d) t = 50 ps. 107 4.4.4 Reduced chemical damage Figure 4.20: Comparison of chemical damage in the cavitation pit due to collapse of empty and gas-filled nanobubbles To estimate the extent of this form of chemical damage on the silica slab, we plot the distribution of SiOH species along circular bins of width 0.5 nm, starting from the center of the slab and progressing outwards (Figure 4.20). Both systems (empty and gas) start with uniform distribution of SiOH across the slab at time t = 0 (green curves in Figure 4.20). We see an increase in the number of SiOH species across the slab, as indicated by the red and blue dashed lines representing the systems filled with gas for t = 25 and 50 ps, respectively, in Figure 4.20. In case of the systems with empty bubble (red and blue solid lines for t = 25 and 50 ps, respectively, the nanojet causes significant cavitation damage to the slab in form of a dent. We observe a shifted peak from the center (see the arrows in 108 Figure 4.20), which is a result of silica surface, eroded by the water nanojet, being flown outwards from the center by the vortex rings formed when the nanojet hits the slab. In case of empty nanobubble, we observe considerable increase in SiOH with time, as evidenced by the solid lines in Figure 4.20. However, in case of gas-filled nanobubble, the number of SiOH is much less compared to that in case of empty nanobubble (dashed lines in Figure 4.20). Thus it is clear that in addition to reduction in mechanical damage, the presence of the gas particles in the nanobubble reduces the chemical damage on the slab. 4.5 Summary In summary, large-scale reactive MD simulations of shock-induced nanobubble collapse reveal an increase in the number of hydronium ions when the nanobubble begins to collapse and water molecules from the bubble periphery rush in to form a focused jet. The nanojet impact on the distal side of the collapsing nanobubble generates a large number of hydronium ions and a secondary shock wave. When the nanojet hits the silica surface, it creates a reflected shock wave and a damage pit whose volume scales linearly with the nanobubble volume. We observe a large number of silanol and hydronium ions in the pit. We also performed reactive MD simulations with gas-filled nanobubbles of the same dimensions and for the same value of the stand-off parameter. These nanobubbles do not collapse completely nor do they cause much damage to the silica surface because the nanojets are much weaker than the nanojets from the collapse of empty bubbles. Pressure and density gradients around partially collapsed gas-filled nanobubbles give rise to vortex rings and enhanced hydronium ion production. 109 Chapter 5. Water confined in Nanoporous Silica 5.1 Background Understanding the behavior of water confined in ultra small spaces[61] is important in widely varying research fields such as life sciences, geosciences, atmospheric sciences, energy and environment.[62] For example, in the rapidly emerging DNA sequencing technology based on nanopore sensing[63], one of the key unresolved issues is how the structure and dynamics of water confined in nanopores sculpted in silica or silicon nitride membranes affect the DNA translocation speed. In geological processes films of nano- confined water mediate crystallization, dissolution, and material transport along stressed surfaces in contact [64-66] and mineralization reactions, such as CO 2 sequestration by carbonation, is limited by water transport through passivating silica-rich nano-porous layers formed at the reaction interfaces[67, 68] Improved understanding of confined water in tight rocks, such as shales or carbonates, is also important to reduce environmental impacts of hydrocarbon production and to ensure permanent CO 2 storage in geological systems. Recent experimental studies indicate that nano-confined water may behave differently than bulk water.[69] It has been suggested that nanoconfinement may provide access to a metastable region, known as the “no man’s land (NML)”, in the phase diagram of bulk water[70]. At ambient pressure, the NML lies between the homogeneous nucleation temperature T h = 235 K and the glass transition temperature T g = 136 K. Below T g , bulk water is frozen either in the low-density amorphous (LDA) or high-density amorphous 110 (HDA) phase depending on the pressure on the system[71]. The existence of LDA and HDA has led to an intriguing speculation, namely whether supercooled water in the NML is a mixture of low-density water (LDW) and high-density water (HDW) that are extensions of LDA and HDA, respectively. A couple of theoretical studies— one based on the assumption that there is a first order phase transition between two kinds of water terminating at a critical point[72, 73] and the other being a singularity-free scenario[74, 75]—support the existence of LDW and HDW in the supercooled region. An alternative theoretical scenario for the phase behavior of metastable water is spinodal decomposition[76, 77], where response functions exhibit anomalous increases close to the spinodal boundary. Recent experimental efforts[78-80] to examine the behavior of LDW and HDW have focused on water confined in nanometer size pores in nanoporous glasses and clays, in glass capillaries and thin films (a few nanometer in thickness) and as droplets suspended in another liquid, e.g., glycerin. Mallamace et al.[81] investigated the behavior of supercooled water confined in hexagonal cylindrical structures in a micelle-templated mesoporous silica (MCM-41-S) matrix with pores of diameter 1.4 nm. Using Fourier transform infrared spectroscopy, they were able to confirm the existence of supercooled LDW over a temperature range of 183 – 273 K. Zhang et al.[69] used neutron scattering to investigate liquid-liquid phase transition in deuterated water confined in the nanopores of MCM-41-S silica matrix. For temperatures ranging between 130 K and 300 K and pressure between 1 and 2,900 bars, they observed significant hysteresis in the density of heavy water, indicating a first-order phase transition between LDW and HDW. 111 Nanoconfined supercooled water also displays unusual transport properties. Under ambient conditions, the self-diffusion coefficient D and viscosity η obey the Stokes- Einstein relation (SER),[82] D ~ T/η. Upon cooling, the viscosity increases and the self- diffusion coefficient decreases dramatically. Below a crossover temperature T x , NMR and quasi-elastic neutron scattering experiments indicate a breakdown of the SER in supercooled water confined in a nanotube or a lysozyme protein layer and also in a methanol mixture. The temperature T x is associated with the fragile-to-strong[102] dynamic crossover: For T > T x the viscosity of water has a power-law dependence on temperature, whereas below T x the viscosity shows Arrhenius temperature dependence and decoupling from the self-diffusion coefficient or the structural relaxation time α. It has been suggested that the decoupling of viscosity and self-diffusion coefficient below T x arises from the presence of dynamic heterogeneities in supercooled water. These heterogeneities have been observed experimentally and in numerical simulations of other supercooled liquids as well.[78] Dynamic heterogeneities are also expected to affect various transport characteristics of nano-confined water, e.g., dynamic correlations in nano-confined water are expected to exhibit Kohlraush-Williams-Watts stretched exponential decay[139] with time: e − t /τ ( ) β . 5.2 Simulation In this work, we examine the structural and dynamical heterogeneities in supercooled water confined in nanopores of amorphous silica. We have performed molecular dynamics (MD) simulations based on reactive force fields for silica and water, which incorporate changes in the local atomic environment and allow bond breaking and bond 112 formation between atoms. The total interatomic potential energy of the system consists of two-body and three-body terms. The two-body terms include steric repulsion, charge transfer, charge-dipole, and van der Waals attractive interactions between atoms. The three-body potential takes into account bond bending and bond stretching between nearest-neighbor atoms. The force field parameters are fitted to experimental measurements and quantum mechanical calculations of structural, thermodynamic, dynamic, and mechanical properties of bulk amorphous silica (a-SiO 2 ) and water. The interatomic potential has been described in detail in section 4.2.1. Figure 5.1: Preparation schedule for Nanoporous silica from β-cristobalite. 113 In our MD simulation, nanoporous silica was prepared by the melt-quench method (Figure 4.2). Figure 5.1 shows the schedule for preparation of nanoporous silica. The starting atomic configuration of the MD system was the β-cristobalite lattice, containing 65836 atoms in a box of volume 99.7×99.7×99.7 Å 3 . The β-cristobalite configuration was chosen for the simple reason that it has the same mass density as bulk a-SiO 2 under ambient conditions. Silicon and oxygen atoms were assigned random initial velocities chosen from the Maxwell distribution and the system was heated to a temperature T = 4,500 K. Periodic boundary conditions were imposed and the equations of motion were integrated with the velocity-Verlet algorithm using a time step of 1 fs. After thermalizing the system for 250 ps, the MD box was expanded in increments of (5Å) 3 over a time period of 5 ps followed by thermalization for 10 ps. Following this expansion and thermalization schedule, we obtained several molten systems with mass densities ranging between 0.1 and 2.2 g/cc. During expansion and thermalization cycles, these systems were always coupled to a thermostat in order to maintain the temperature at 4,500K. Amorphous systems were obtained by quenching the molten systems from 4,500 K to 10 K over 250 ps. Each amorphous silica system thus obtained was heated to 300 K over 250 ps followed by thermalization for 500 ps. Rapid changes in internal pressure due to expansion created cavities in the molten systems and these cavities were quenched into nanopores in the amorphous systems. 5.3 Results and Discussion We focus on one of the nanoporous silica systems prepared by the melt-quench method. The mass density of the system is 0.5 g/cc. Figure 5.2a shows the network of pores and 114 interatomic bonds in the nanoporous system. The nanoporous silica network is shown in red, while the empty space has been color coded by the distance of the nearest atom. We found pores by dividing the system into voxels and then using the depth-first-search algorithm to connect nearest-neighbor voxels. The voxel size was chosen to be 2 Å, which corresponds to the first minimum in the Si-O radial distribution function in bulk amorphous silica (a-SiO 2 ). We found a wide range of nanopore shapes and sizes in the nanoporous silica system at 0.5 g/cc. Figure 5.2: (a) Structure of amorphous silica network with density 0.5 g/cc (red). The empty space is color coded by the distance of the nearest atom. (b) Close-up view of the silica network (red) around the largest pore, shown in cyan. c) Distribution of pores with respect to anisotropy. Perfect sphere has anisotropy = 0. d) Distribution of pores with respect to shape parameter. Negative shape parameter corresponds to oblate spheroid, while the pores with positive shape parameter resemble prolate spheroid. 115 We characterized the nanopore sizes by calculating their radii of gyration, R G , and nanopore shapes by computing the inertia tensor[140], (5.1) where R C α = α-coordinate of the position of the center of mass of i th voxel in a pore and N is the total number of voxels in the pore. From the inertia tensor we determined the anisotropy parameter, Δ, and shape parameter, S: (5.2) where . For a spherical pores, Δ = 0. The values of S range between -1/4 and 2; for oblate ellipsoidal pores S < 0 and for prolate ellipsoidal pores S > 0. Fig. 5.2b shows a close-up view of the largest nanopore and the surrounding Si-O bond network. The pore’s radius of gyration is 2.24 nm and its anisotropy parameter Δ = 0.32. The value of the shape parameters S = 0.33 indicates that this pore is a prolate ellipsoid.[140] Fig. 5.2c shows that nanopores have a wide distribution of Δ. Nearly 50% of the nanopores are spherical (Δ ≈ 0). The distribution of shape parameters in Fig. 5.2d shows that approximately 20% of the nanopores in the system are spherical, 20% are oblate ellipsoids, and the rest are prolate ellipsoids. rameters and l p we also describe the unusual structural characteristics of the ribosome, a large ribonucleoprotein complex. II. METHODS RNA structures. We computed several quantities to char- acterize the shapes of RNA using the atomic coordinates of their structures determined by x-ray crystallography, NMR, or cryo-electron-microscopy !cryo-EM". The coordinates for all RNA structures were obtained from the Protein Data Bank !PDB". 8 Our analysis is performed for over 1185 indi- vidual RNA chains with the number of nucleotides N!10 found in 642 RNArelated PDB files as of June 2005.Among these, 195 RNA chain structures are monomers, and the rest of the chains are part of oligomers or appear in complexes with other RNAmolecules or proteins. Structural features in the monomeric form can be different from those determined in an oligomer or complex because the intermolecular inter- action can affect the individual chain structure. Therefore, we analyzed the two groups of structures separately. For comparison,wehavealsocalculatedshapecharacteristicsfor a data set of proteins. The results for proteins enable us to assess certain unusual features of RNA-protein interactions especially in the ribosome. Size. The radius of gyration !R G " is an indicator of the overall size of RNA. The value of R G 2 , which can be mea- sured using small angle x-ray or neutron scattering, is calcu- lated using R G 2 = 1 2 # i M m i# j M m j # i M # j M m i m j !r i −r j " 2 , !1" where M is the number of atoms in the molecule and m i is the mass of the ith atom. In the calculation of R G 2 for RNA structure we used only the coordinates of the heavy atoms !C, N, O, P". Shape. The deviation from the spherical shape is charac- terized by the asphericity" and the shape parameter S, both of which are calculated from the inertia tensor, 11,12 T #$ = 1 2 # i M m i# j M m j # i M # j M m i m j !r i# −r j# "!r i$ −r j$ " = 1 # i M m i # i M m i !r i# −R C# "!r i$ −R C$ ", !2" where r i# is the#th Cartesian component of the position of atom i and R C =# i M m i r i /# i M m i is the center of mass. The square of the radius of gyration is R G 2 =trT. The eigenvalues % 1 , % 2 , and % 3 of the matrix T are the squares of the three principal radii of gyration. The extent of asphericity is char- acterized using"!0&"&1", "= 3 2 # i=1 3 !% i −% ¯ " 2 !trT" 2 , !3" where% ¯ =!% 1 +% 2 +% 3 "/3. For a perfect sphere"=0. Devia- tionfrom"=0 indicatestheextentofanisotropy.Theoverall shape of a molecule is assessed using S=27 $ i=1 3 !% i −% ¯ " !trT" 3 , !4" which satisfies the bound −1/4&S&2. Negative values of S correspond to oblate ellipsoids and S!0 are prolate ellip- soids. Most studies of packing in proteins and RNAs involve tessellation of space which always introduces certain amount of arbitrariness. 4,13 In contrast, the shape parameters" and S are directly computed using only the atomic coordinates. Knowledge of" and S is important in determining the over- all motion of RNA and their interaction with other biomol- ecules. Persistence length. A parameter that describes the flex- ibility of biomolecules is the persistence length l p , which is most clearly defined by assuming that RNAstructures can be described by a polymer model. Based on previous experi- mental studies it is suspected that the statistical properties of dsDNA, 10,14 ssDNA, 15 and RNA !Refs. 9 and 16" can be described using the wormlike chain !WLC" model. ForWLC models l p can be estimated provided that the distribution of the mean end-to-end distance R E or R G is known. Exact cal- culation of neither P E !R E " nor P!R G "is possible for WLC.A simple and accurate theoretical expression has been derived for P!R E " of wormlike chain using the mean field approximation. 17,18 The resulting distribution, which is in good agreement with computer simulations, 19 is P WLC !r E "= 4'Cr E 2 !1−r E 2 " 9/2 exp % − 3t 4!1−r E 2 " & , !5" where r E 'R E /L and t'L/l p . L is the contour length. For RNAmolecules, which from the perspective of polymers can be viewed as branched polyelectrolyte chains, the contour lengthisalsoanunknownparameter.Thenormalizationcon- stant C=1/(' 3/2 e −# # −3/2 !1+3# −1 +15/4# −2 ") with an # =3t/4. When l p is small P WLC !r E " reduces to a Gaussian chain, whereas for large l p P WLC !r E " approaches the rod limit as r E →1. Although direct measurements of P!R E " for biomol- ecules are not routinely performed it is conceivable that P!R E " may be obtained using single molecule fluorescence resonanceenergytransfer !FRET"experiments.However,the distance distribution function P!r" can be measured using SAXS experiments. 20–22 Based on general arguments, we ex- pect that the distribution functions P!r" and P!R E " should coincide provided that r(R G . Because *R E 2 +,*R G 2 +,Ll p for WLC provided that L is large it follows that P!r" should decay for large r as 194905-2 Hyeon, Dima, and Thirumalai J. Chem. Phys. 125,194905 !2006" Downloaded 19 Nov 2006 to 132.239.69.25. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp rameters and l p we also describe the unusual structural characteristics of the ribosome, a large ribonucleoprotein complex. II. METHODS RNA structures. We computed several quantities to char- acterize the shapes of RNA using the atomic coordinates of their structures determined by x-ray crystallography, NMR, or cryo-electron-microscopy !cryo-EM". The coordinates for all RNA structures were obtained from the Protein Data Bank !PDB". 8 Our analysis is performed for over 1185 indi- vidual RNA chains with the number of nucleotides N!10 found in 642 RNArelated PDB files as of June 2005.Among these, 195 RNA chain structures are monomers, and the rest of the chains are part of oligomers or appear in complexes with other RNAmolecules or proteins. Structural features in the monomeric form can be different from those determined in an oligomer or complex because the intermolecular inter- action can affect the individual chain structure. Therefore, we analyzed the two groups of structures separately. For comparison,wehavealsocalculatedshapecharacteristicsfor a data set of proteins. The results for proteins enable us to assess certain unusual features of RNA-protein interactions especially in the ribosome. Size. The radius of gyration !R G " is an indicator of the overall size of RNA. The value of R G 2 , which can be mea- sured using small angle x-ray or neutron scattering, is calcu- lated using R G 2 = 1 2 # i M m i# j M m j # i M # j M m i m j !r i −r j " 2 , !1" where M is the number of atoms in the molecule and m i is the mass of the ith atom. In the calculation of R G 2 for RNA structure we used only the coordinates of the heavy atoms !C, N, O, P". Shape. The deviation from the spherical shape is charac- terized by the asphericity" and the shape parameter S, both of which are calculated from the inertia tensor, 11,12 T #$ = 1 2 # i M m i# j M m j # i M # j M m i m j !r i# −r j# "!r i$ −r j$ " = 1 # i M m i # i M m i !r i# −R C# "!r i$ −R C$ ", !2" where r i# is the#th Cartesian component of the position of atom i and R C =# i M m i r i /# i M m i is the center of mass. The square of the radius of gyration is R G 2 =trT. The eigenvalues % 1 , % 2 , and % 3 of the matrix T are the squares of the three principal radii of gyration. The extent of asphericity is char- acterized using"!0&"&1", "= 3 2 # i=1 3 !% i −% ¯ " 2 !trT" 2 , !3" where% ¯ =!% 1 +% 2 +% 3 "/3. For a perfect sphere"=0. Devia- tionfrom"=0 indicatestheextentofanisotropy.Theoverall shape of a molecule is assessed using S=27 $ i=1 3 !% i −% ¯ " !trT" 3 , !4" which satisfies the bound −1/4&S&2. Negative values of S correspond to oblate ellipsoids and S!0 are prolate ellip- soids. Most studies of packing in proteins and RNAs involve tessellation of space which always introduces certain amount of arbitrariness. 4,13 In contrast, the shape parameters" and S are directly computed using only the atomic coordinates. Knowledge of" and S is important in determining the over- all motion of RNA and their interaction with other biomol- ecules. Persistence length. A parameter that describes the flex- ibility of biomolecules is the persistence length l p , which is most clearly defined by assuming that RNAstructures can be described by a polymer model. Based on previous experi- mental studies it is suspected that the statistical properties of dsDNA, 10,14 ssDNA, 15 and RNA !Refs. 9 and 16" can be described using the wormlike chain !WLC" model. ForWLC models l p can be estimated provided that the distribution of the mean end-to-end distance R E or R G is known. Exact cal- culation of neither P E !R E " nor P!R G "is possible for WLC.A simple and accurate theoretical expression has been derived for P!R E " of wormlike chain using the mean field approximation. 17,18 The resulting distribution, which is in good agreement with computer simulations, 19 is P WLC !r E "= 4'Cr E 2 !1−r E 2 " 9/2 exp % − 3t 4!1−r E 2 " & , !5" where r E 'R E /L and t'L/l p . L is the contour length. For RNAmolecules, which from the perspective of polymers can be viewed as branched polyelectrolyte chains, the contour lengthisalsoanunknownparameter.Thenormalizationcon- stant C=1/(' 3/2 e −# # −3/2 !1+3# −1 +15/4# −2 ") with an # =3t/4. When l p is small P WLC !r E " reduces to a Gaussian chain, whereas for large l p P WLC !r E " approaches the rod limit as r E →1. Although direct measurements of P!R E " for biomol- ecules are not routinely performed it is conceivable that P!R E " may be obtained using single molecule fluorescence resonanceenergytransfer !FRET"experiments.However,the distance distribution function P!r" can be measured using SAXS experiments. 20–22 Based on general arguments, we ex- pect that the distribution functions P!r" and P!R E " should coincide provided that r(R G . Because *R E 2 +,*R G 2 +,Ll p for WLC provided that L is large it follows that P!r" should decay for large r as 194905-2 Hyeon, Dima, and Thirumalai J. Chem. Phys. 125,194905 !2006" Downloaded 19 Nov 2006 to 132.239.69.25. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp λ = (λ 1 +λ 2 +λ 3 )/ 3 116 Water molecules were embedded in nanopores after passivating the under-coordinated silicon atoms with –OH and dangling oxygen atoms with +H groups. The passivation schedule is as follows: An oxygen atom was placed 2 Å away from an under-coordinated Si atom. Following this, a hydrogen atom was placed 2 Å away from any under- coordinated O atom in the system. Any extra H or O atom was removed to maintain zero net charge. The energy of the system was minimized using steepest descent. The system was then heated to 300 K resulting in passivated, nanoporous amorphous silica network. The passivated system was divided into 2 Å 3 voxels and a water molecule was placed in any empty voxel. The system was heated to 300 K and themalized for 250 ps to obtain a nanoporous amorphous silica network embedded with water. All in all, the MD box contained 404,379 atoms including 105,752 water molecules. Again periodic boundary conditions were imposed and the silica system with water in the nanopores was quenched by the steepest-descent method. Subsequently, MD simulations were carried out in the NVE ensemble to obtain a well-thermalized system at room temperature. The presence of hydrogen atoms made it necessary to reduce the time step to 0.5 fs. From the room temperature system, we obtained supercooled water confined in nanoporous silica at temperatures T = 250, 200, 150 and 100K. Fig. 5.3a is a snapshot of water molecules distributed in pores. Here the silica network is displayed in red and water molecules in green, yellow and blue. Green and blue represent water molecules with mass densities below 0.9 g/cc and above 1.1 g/cc, respectively; and water molecules with mass densities between 0.9 and 1.1 g/cc are shown here in yellow. We calculated the mass density of a water molecule by constructing a Voronoi cell[141] of nearest-neighbor water molecules and then dividing the mass of the molecule by the Voronoi cell volume. Fig. 5.3b shows 117 the mass density distribution of water molecules at various temperatures. The distributions sharpen but their shapes are hardly affected by cooling. Figure 5.3: (a) Heterogeneous distribution of water in nanoporous silica. Silica network is shown in red, while yellow, green and blue represent the local density of water in the network. (b) Distribution of local density of water averaged over the entire system at various temperatures. 118 Figure 5.4 shows spatial correlations in nanoconfined LDW and HDW in the supercooled state at T = 100 K. Here the mass density of LDW molecules ranges between 0.8 and 0.9 g/cc and the mass density of HDW molecules lies between 1.1 and 1.3 g/cc. In Fig. 5.4a we show the tetrahedral order parameters for LDW and HDW molecules confined in the nanopores of a-SiO 2 . The tetrahedral order parameter[108] is obtained from, (5.3) where ψ ijk represents the angle between the i th , j th and k th water molecules, with k th molecule at the center. The order parameters confirm that both LDW and HDW molecules are four-fold coordinated. Fig. 5.4b shows the oxygen-oxygen radial distribution functions, g OO (r), for supercooled LDW and HDW at T = 100 K. The first peaks in g OO (r) are located at r = 3.25 and 3.0 for LDW and HDW respectively. These peaks along with the tetrahedral order parameters in Fig. 5.4a indicate the existences of first solvation shells in both LDW and HDW. The g OO (r) for LDW has a broad second peak around 4.5 Å, which indicates the presence of a second solvation shell. In contrast, the second peak in the radial distribution function for HDW is shifted to 3.4 Å and there is a shallow second minimum around the second peak in the g OO (r) for LDW. This behavior suggests that the second solvation shell in HDW is broken. We found additional evidence for broken second solvation shell in the O-O-O bond-angle distribution for HDW. The presence of a second solvation shell would give rise to a peak at 70° in the O-O-O bond angle distribution, whereas our calculations show that the distribution for supercooled HDW is featureless. Tetrahedral order parameter Q k =1 3 8 3 X i=1 4 X j=i ✓ cos ijk + 1 3 ◆ Nanoconfined water has a high degree of tetrahedral order Unusual Properties of Water Confined in Nanoporous Silica Glasses C. N. Kirkemo 119 Figure 5.4: Local strucutre in LDW and HDW 120 Our results agree remarkably well with neutron scattering study of LDW and HDW by Soper and Ricci[80]. They examined the local molecular structure of water in a temperature-pressure regime where the anomalous physical properties of water are manifested most prominently. Using neutron diffraction with isotope substitution, they measured the partial structure factors of water as a function of pressure around the ice I/ice III triple point, namely at temperature T = 251 K and pressure P = 209 MPa. With an increase in pressure from 26 MPa to 400 MPa, they observed a continuous transformation from the LDW to HDW phase at T = 268 K. Assuming that the structure of water is a linear combination of the structures of LDW and HDW, they extracted from their data the mass densities of LDW and HDW to be 0.88 g/cc and 1.2 g/cc, respectively. Combining the measurements of oxygen-oxygen structure factors with simulations based on empirical pressure structure refinement, Soper and Ricci obtained the radial distribution functions for LDW and HDW. They found the first peaks in g OO (r) for LDW and HDW at ~2.9 Å and the second peaks around 4.5 Å and 3.4 Å, respectively. Furthermore, the radial distribution function for HDW had a shallow minimum around 4.5 Å. These experimental results imply that the first and second solvation shells in LDW are fully formed and that the second solvation shell has collapsed in HDW. In addition to spatial heterogeneities, we also find dynamical heterogeneities in supercooled nanoconfined water. Fig. 5.5 shows the mean square displacement as a function of time for a tagged water molecule confined in the largest nanopore in the nanoporous silica system. After the initial increase in the ballistic regime, the MSD exhibits plateuas and intervening diffusive regime. At a given temperature the durations 121 of plateau and diffusive regimes vary from particle to particle, although on average the plateaus last longer upon cooling. Analyzing the dynamics of a tagged water molecule and its neighbors, we find that the plateau region arises from the trapping of the tagged molecule in a cage formed by its nearest-neighbor water molecules.[107] When some of the nearest neighbors begin to diffuse, the cage breaks allowing the tagged particle to enter the diffusive regime. Figure 5.5: Mean square displacement of a water molecule in nanoconfined water. The molecule is seen to break from one cage (green) and move into another cage (red) in 0.2 ns. 122 The dynamics of water molecules is further analyzed in terms of a cage correlation function: (5.4) where the array l i (0) contains the IDs of water molecules which are nearest neighbors of the i th water molecule time t = 0 and l i (t) is an array of nearest neighbor IDs at a later time t. Figure 5.6: Cage correlation function vs time at (a) T = 150 K, (b) T = 200 K, and (c) T = 250 K. (d) Cage correlation function vs normalized time for various temperatures. Stretched exponential curves for all temperatures collapse onto a single stretched exponential curve, with exponent β = d/(d+2). Dynamics Cage correlation function: C(t)= < ` i (0) · ` i (t)> < ` i (0) 2 > Rabani, Gezelter, and Berne (1997) Water molecules di↵ use from cage to cage Unusual Properties of Water Confined in Nanoporous Silica Glasses C. N. Kirkemo 123 Fig. 5.6 shows the decay of the cage correlation function with time in supercooled nanoconfined water at temperatures T = 150, 200 and 250 K. These curves, along with the MD results for the cage correlation function at T = 100 K, collapse into a single universal curve presented in Fig. 5.6d. We find the stretched exponential function, e − t τ " # $ % & ' β , with β = 0.6 provides the best fit to cage correlation functions at all temperatures below T = 250 K. Note the value of β = 0.6 coincides with the simple relation β = d/(d+2), where d is the spatial dimensionality of the system. The MD result for the stretched exponent β in supercooled nano-confined water agrees with an exact result for a liquid diffusing in a medium with static random traps. In 1982 Grassberger and Proccacia[95] showed that the asymptotic decay of the density of diffusing particles in a medium with static random traps is slower than . Their analysis was based on the reasoning that trap-free spaces in the medium controlled the long-time diffusive behavior even though the probability of finding such spaces was small. A year later, Kayser and Hubbard[97] showed that in the case of a random distribution of static spherical traps the upper bound on the asymptotic decay of the density of diffusing particles is the same as the lower bound found by Grassberger and Proccacia[95]. Both these results follow from an exact asymptotic result proved by Donsker and Varadhan[142] in 1979. They showed that the survival probability P(t) averaged over all Brownian paths is , where denotes the set of points at a distance less than ε from the Brownian path at time t and <…> denotes the average over Brownian paths. In the limit , Donsker and Varadhan showed: exp −t /τ ( ) d d+2 ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ exp −n× S t ε { } ( ) S t ε t→∞ 124 (5.5) where γ is the lowest eigenvalue of the Laplacian operator in a sphere of unit volume with no boundary conditions. 5.4 Summary We investigate structural and dynamic behavior of water confined inside a nanoporous silica network with density of 0.5 g/cc. The nanoporous silica network consists of multiple pores of different size and shape. Approximately 20% of the nanopores in the system are spherical, 20% are oblate ellipsoids, and the rest are prolate ellipsoids. Nanoconfined water displayed a heterogeneous distribution of density, revealing localized Low Density Water (LDW) and High Density Water (HDW). Our results on the local structure in LDW and HDW agree very well with experiments. We also found that both LDW and HDW maintain good tetrahedral order. Under nanoconfinement, a water molecule is trapped in a cage formed by the neighboring water molecules. Mean square displacement calculations show that the trapped molecule breaks out the cage and then gets trapped in another cage, with a new set of neighbors. The decay of the cage is characterized by cage correlation function, which is found to follow exponential decay at all temperatures. However, the curves at different temperatures fall onto a single stretched exponential curve when plotted against normalized time. The coefficient of the stretched exponential decay agrees with an exact result for a liquid diffusing in a medium with static random traps. 1 t d/ d+2 ( ) ln exp −n×S t ε { } = −n 2/(d+2) d + 2 d ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ 2γ d ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ d/(d+2) Δ / 2 125 Chapter 6. Conclusion As technology advances further, understanding of material behavior at atomistic scale becomes increasingly important. Due to extremely small length and time scales involved, the current experimental techniques available cannot measure atomistic behavior in detail. Numerical methods are being used all around in order to bridge this gap between experimental limitations and a comprehensive understanding of atomistic behavior of nanomaterials. While any model has its own assumptions and limitations, increasingly accurate predictive models for a wide variety of materials and properties are constantly being proposed and existing models being improved upon over the years. A huge surge in computing power and resources available to researchers in recent years has played a large part in making numerical methods an important aspect of any research at the atomistic level. In this work, we have investigated structural, dynamic, mechanical, and chemical behavior of a range of nanomaterials under extreme environments using large-scale reactive molecular dynamics simulations. Aluminum nanoparticles (AlNPs) are widely used in energetic applications such as explosives, rocket propulsions and other pyrotechniques due their high energy release rate as compared to micron-sized particles. Although there have been many experimental studies on AlNPs, focusing on combustion mechanism, energy release rates, and size effects, atomistic level of understanding the oxidation behavior is difficult due to small spatio-temporal scale. We have used molecular dynamics to simulate the oxidation dynamics in a multi-million-atom system of three spherical ANPs. Each ANP consists of a spherical aluminum core, 400 Å in diameter, surrounded by a 30 Å thick amorphous 126 alumina shell. The ANPs were arranged in a linear fashion and the left and the right (outer) ANPs were heated to a temperature of 1200 K. The central ANP and the background oxygen were kept at 500 K. Oxidation reactions are initiated inside the outer ANPs due to the initial heat and release more energy as the reactions are exothermic in nature. As the reactions continue, the outer ANPs expand and deform. Initially, the oxidation reactions inside the central ANP are localized to the fused zones of contact. The shells of the outer ANPs melt first and begin to penetrate the central ANP through the zones of contact. In addition to bringing in external heat via convective transfer, the outer atoms also initiate exothermic oxidation reactions inside the central ANP. Thus, although slow in the beginning, the central ANP oxidizes at a faster rate as compared to the outer ANPs. As observed in previous studies of oxidation of single ANP, three stages of oxidation are identified in outer ANPs as well as the central ANP. Stage I, which lasts typically for 0.1 ns, is characterized by a decreasing rate of change of temperature with time. The behavior is attributed to the oxidation reactions being confined to the core-shell interface and the fused zones of contact. Stage II has an increasing rate of change of temperature with time due to increased number of oxidation reactions following the melting of the alumina shell. Stage II lasts longer and has a higher peak in case of the central ANP because it gains heat due to the penetration of hot atoms from outer ANPs as well as local exothermic oxidation reactions initiated by the atoms from the outer ANPs. As the reactions continue, the shell weakens resulting in the ejection of atoms from the ANPs into the background oxygen environment. Thus begins Stage III, where the rate of change of temperature decreases with time. The penetration front consists of three layers: 1) atoms from the shell of the central ANP; 2) atoms from the shells of the outer ANPs; 127 and 3) atoms from the cores of the outer ANPs. The layers of the penetration front maintain their relative position with respect to each other. This means that the outer ANPs push the shell of the central ANP towards the center of the central ANP but they do no break through it. The speed of penetration is calculated as 54 m/s, which is within range of experimentally observed flame propagation rates. By 1 ns, all the three ANPs fuse together into an ellipsoidal aggregate at an average temperature of 5500 K. Our simulation results indicate that the oxidizing ANPs have a huge effect on the reaction dynamics of a neighboring ANP. A pre-heated ANP can initiate oxidation reactions in a neighboring ANP to a level of self-sustenance. Proximity to other ANPs and the pre- heating schedule can dictate the oxidation behavior of an ANP. Hence, any theoretical or experimental study of the oxidation dynamics of aluminum nanoparticle aggregate should give due consideration to the above two factors for a comprehensive understanding. Continuing our studies of mechano-chemical behavior under extreme conditions at nanoscale, we investigate the effect of shock-induced collapse of a water nanobubble on a near-by silica slab. Cavitation bubbles are generated when there is a rapid pressure drop inside hydrodynamic systems. Sudden change in pressure results in localized change in velocities inside the fluid, which lead to the formation of nanobubbles. Such bubbles are observed in propulsion systems, pumps, cooling systems in power plants, hydraulic mining machinery, ultrasound, etc. Since silica has been an archetypal material for studying stress corrosion cracking, shock, indentation, and behavior under high-pressure, cavitation erosion damage at the silica-water interface has been a topic of recent 128 experimental study. We thus use the silica-water system to study the factors affecting mechanical as well as chemical damage due to shock-induced collapse of nanobubble. We performed billion-atom reactive molecular dynamics (MD) simulations of shock- induced collapse of a nanobubble near a silica slab. The simulations were performed on the full BlueGene/P supercomputer (163,840 processors) at the Argonne Leadership Computing Facility. In the billion-atom MD simulations, the system consists of a 60 nm thick silica slab placed in an MD box of dimensions 285 × 200 × 200 nm 3 in the x, y and z directions respectively. The rest of the MD box contains water molecules and a spherical nanobubble of diameter 97.6 nm. The initial distance between the center of the bubble and the silica surface is 90 nm, and the stand-off parameter, s d = 1.84. We applied shock with a particle velocity of 2 km/s and found the shock velocity in water to be 5 km/s. The density of water in the shock region is 1.5 g/cc and the pressure is 10 GPa. We find the structure of water to change drastically in the shock region with the average number of nearest neighbors increasing from 4 to 8. We perform two simulations: one with an empty bubble and the other with the bubble filled with gas-particles. As the empty nanobubble collapses, a high-speed water jet is formed during the collapse of the empty nanobubble, which proceeds to hit the slab resulting in the formation of a cavitation erosion pit. The nanojet reaches a speed of 7.8 km/s, which is significantly higher than the shock velocity (5 km/s). Chemical activity in water increases significantly with the onset of bubble collapse and nanojet formation. We observe an increase in the hydronium ion production with the nanojet growth. We estimate the water-hammer pressure as 18 GPa, which is consistent with theoretical predictions. The increased pressure in the high-speed nanojet results in increased concentration of hydronium ions, 129 thereby altering the chemical composition of water. The cavitation pit is spherical in shape and almost 52% of the silica surface is damaged due to erosion caused by the high- speed nanojet. Our results indicate that the cavitation pit volume scales linearly with the bubble volume. We also observe a large number of hydrogen ions forming bonds with Si and O in the damaged region of the silica slab. The pressure in the silica slab decreases radially outward from the center of the cross-section, which is the point of impact with the nanojet. The peak pressure obtained in the silica slab is 18 GPa. In case of gas-filled nanobubble, it deforms into a toroid-shaped aggregate following the collapse. Similar toroid-shaped deformation has also been seen in shock-wave experiments on spherical and cylindrical inhomogeneity, and in continuum simulations of shock-induced collapse of a bubble near a wall. Gas-filled nanobubbles do not collapse completely nor do they cause much damage to the silica surface because the nanojets are much weaker than the nanojets from the collapse of empty bubbles. Pressure and density gradients around partially collapsed gas-filled nanobubbles give rise to vortex rings and enhanced hydronium ion production. The pressure profile across the silica slab during the shock-induced collapse of gas-filled nanobubble is uniform across the cross-section. Comparing this with the case of an empty bubble, the pressure in the pit region decreases radially as one moves away from the pit center. In addition, we found that for the same standoff parameter, system with larger initial bubble has a higher depth of dent as compared to the system with smaller initial bubble. The same trend is observed in the systems with gas bubbles. In summary, for the same 130 standoff parameter, the larger bubble will cause more damage than a smaller bubble. However, the presence of inert gas particles significantly reduces the extent of mechanical as well as chemical damage on the silica slab. Thus, our results show that the gas-filled nanobubbles can be used in water to arrest an incoming shock wave and protect the surface behind it from cavitation erosion as well as chemical leeching. Exploring silica-water systems further, we examine the structural and dynamical heterogeneities in supercooled water confined in nanopores of amorphous silica. The structure and dynamics of water confined in nanoporous silica are different from that of bulk water, and insight into the properties of confined water is important for our understanding of many geological and biological processes. Nanoporous silica has a wide range of technological applications because it is easy to tune the size of pores and their morphologies and to functionalize pore surfaces with a variety of molecular moieties. In this work, we study nanoporous silica prepared by the melt-quench method, which is filled with water and has a mass density of 0.5 g/cc. The system contained 404,379 atoms including 105,752 water molecules. The nanoporous silica network has a lot of pores with a wide range in shape and size. The distribution of shape parameters shows that approximately 20% of the nanopores in the system are spherical, 20% are oblate ellipsoids, and the rest are prolate ellipsoids. From the room temperature system, we obtained supercooled water confined in nanoporous silica at temperatures T = 250, 200, 150 and 100K. Heterogeneous distribution of local density in water reveals the co- existence of Low Density Water (LDW) and High Density Water (HDW) under nanoconfinement. Our results agree remarkably well with neutron scattering study of 131 LDW and HDW. In addition to spatial heterogeneities, we also find dynamical heterogeneities in supercooled nanoconfined water. Under nanoconfinement, a water molecule is trapped in a cage formed by the neighboring water molecules. Mean square displacement calculations show that the trapped molecule breaks out the cage and then gets trapped in another cage, with a new set of neighbors. The decay of the cage is characterized by cage correlation function, which is found to follow exponential decay at all temperatures. However, the curves at different temperatures fall onto a single stretched exponential curve when plotted against normalized time. The coefficient of the stretched exponential decay agrees with an exact result for a liquid diffusing in a medium with static random traps. Thus we have performed large-scale atomistic simulations using reactive molecular dynamics on a variety of nanosystems with a wide range of applications such as energetics, rocket propulsion, cavitation erosion, shock wave lithotripsy, nanoporous membranes, drug delivery, etc. These simulations have provided useful insights into the mechanical, chemical, structural and dynamic behavior at atomistic levels under extreme environments. Our results can be helpful in designing future experiments to obtain a comprehensive understanding and underlying mechanisms of such highly dynamic processes. 132 References [1] A. Gromov et al., Propellants Explosives Pyrotechnics 31 (2006). [2] A. Ilyin et al., Propellants Explosives Pyrotechnics 27 (2002). [3] L. De Luca et al., Combustion, explosion, and shock waves 41 (2005). [4] L. Galfetti et al., Journal of Physics: Condensed Matter 18 (2006). [5] L. Meda et al., Materials Science and Engineering: C 27 (2007). [6] R. G. Sarawadekar, and J. P. Agrawal, Defence Science Journal 58 (2008). [7] B. S. Bockmon et al., Journal of Applied Physics 98 (2005). [8] D. D. Dlott, Materials Science and Technology 22 (2006). [9] M. L. Pantoya, and J. J. Granier, Propellants, Explosives, Pyrotechnics 30 (2005). [10] A. Ermoline, M. Schoenitz, and E. L. Dreizin, Combustion and Flame 158 (2011). [11] Y. Huang et al., Proceedings of the Combustion Institute 31 (2007). [12] Y. S. Kwon et al., Combustion and Flame 133 (2003). [13] V. I. Levitas, Combustion and Flame 156 (2009). [14] V. I. Levitas, M. L. Pantoya, and K. W. Watson, Appl Phys Lett 92 (2008). [15] A. Rai et al., Combustion Theory and Modelling 10 (2006). [16] V. I. Levitas et al., Appl Phys Lett 89 (2006). [17] P. Puri, and V. Yang, The Journal of Physical Chemistry C 111 (2007). [18] V. I. Levitas et al., Journal of Applied Physics 101 (2007). [19] V. I. Levitas, M. L. Pantoya, and B. Dikici, Appl Phys Lett 92 (2008). [20] B. W. Asay et al., Propellants, Explosives, Pyrotechnics 29 (2004). [21] K. W. Watson, M. L. Pantoya, and V. I. Levitas, Combustion and Flame 155 (2008). [22] G. Hall, and J. M. Watt, Modern numerical methods for ordinary differential equations (Clarendon Press, Oxford, 1976). [23] W. Wang et al., Appl Phys Lett 95 (2009). [24] W. Wang et al., Appl Phys Lett 96 (2010). [25] D. Eakins, and N. N. Thadhani, Journal of Applied Physics 101 (2007). [26] S. S. Mahajan, and G. Subbarayan, Physical Review E 76 (2007). [27] D. S. Wen, L. Zhang, and Y. R. He, Heat and Mass Transfer 45 (2009). [28] R. Clark, in Department of Physics and Astronomy (University of Southern California, Los Angeles, 2010). [29] R. E. A. Arndt, Annual Review of Fluid Mechanics 13 (1981). [30] E. A. Brujan et al., Journal of Fluid Mechanics 433 (2001). [31] C. J. Heathcock, B. E. Protheroe, and A. Ball, Wear 81 (1982). [32] M. Johnson et al., Wear 140 (1990). [33] A. Karimi, and J. L. Martin, International Materials Reviews 31 (1986). [34] B. Karunamurthy et al., Tribology International 43 (2010). [35] C. J. Lin, K. C. Chen, and J. L. He, Wear 261 (2006). [36] Y. Nagaya et al., Journal of Environment and Engineering 5 (2010). [37] R. H. Richman, and W. P. McNaughton, Wear 140 (1990). [38] R. H. Richman, A. S. Rao, and D. E. Hodgson, Wear 157 (1992). [39] A. Vogel, W. Lauterborn, and R. Timmi, Journal of fluid mechanics 206 (1989). [40] J. R. Haines et al., Journal of Nuclear Materials 343 (2005). [41] M. Ida, T. Naoe, and M. Futakawa, Physical Review E 76 (2007). 133 [42] H. Takahira, T. Matsuno, and K. Shuto, Fluid Dynamics Research 40 (2008). [43] M. Lee et al., Wear 255 (2003). [44] C. D. Lorenz et al., Journal of Computational and Theoretical Nanoscience 7 (2010). [45] J. M. D. Lane et al., Langmuir 24 (2008). [46] D. V. Andreeva et al., Small 8 (2012). [47] D. G. Shchukin et al., Advanced Materials 23 (2011). [48] M. Virot et al., The Journal of Physical Chemistry C 114 (2010). [49] C. D. Ohl, and R. Ikink, Phys Rev Lett 90 (2003). [50] M. Delius, F. Ueberle, and W. Eisenmenger, Ultrasound in Medicine & Biology 24 (1998). [51] A. Choubey et al., Appl Phys Lett 98 (2011). [52] Hitachi-GE, E-Journal of Advanced Maintenance 1. [53] N. K. Bourne, and J. E. Field, Journal of Applied Physics 78 (1995). [54] W. Lauterborn, and W. Hentschel, Ultrasonics 23 (1985). [55] O. Lindau, and W. Lauterborn, Journal of Fluid Mechanics, 479 (2003). [56] A. Philipp, Journal of fluid mechanics 361 (1998). [57] S. Succi, G. Smith, and E. Kaxiras, Journal of statistical physics 107 (2002). [58] D. D. Dlott, Annual Review of Physical Chemistry 50 (1999). [59] Z. A. Dreger et al., The Journal of Physical Chemistry B 106 (2002). [60] A. Strachan et al., Phys Rev Lett 91 (2003). [61] P. Gallo, M. Rovere, and S.-H. Chen, The Journal of Physical Chemistry Letters 1 (2010). [62] J. C. Rasaiah, S. Garde, and G. Hummer, Annu. Rev. Phys. Chem. 59 (2008). [63] P. Xie et al., Nature nanotechnology 7 (2011). [64] J. N. Israelachvili, Intermolecular and surface forces: revised third edition (Academic press, 2011). [65] H. C. Sorby, Journal of the Franklin Institute 77 (1864). [66] P. K. Weyl, Journal of Geophysical Research 64 (1959). [67] H. Béarat et al., Environmental science & technology 40 (2006). [68] H. E. King, O. Plümper, and A. Putnis, Environmental science & technology 44 (2010). [69] Y. Zhang et al., Proceedings of the National Academy of Sciences 108 (2011). [70] O. Mishima, and H. E. Stanley, Nature 396 (1998). [71] N. Giovambattista, H. E. Stanley, and F. Sciortino, Physical Review E 72 (2005). [72] C. Jeffery, and P. Austin, The Journal of Chemical Physics 110 (1999). [73] E. Ponyatovskii, V. Sinand, and T. Pozdnyakova, JETP LETTERS C/C OF PIS'MA V ZHURNAL EKSPERIMENTAL'NOI TEORETICHESKOI FIZIKI 60 (1994). [74] P. H. Poole et al., Phys Rev Lett 73 (1994). [75] S. Sastry, P. G. Debenedetti, and F. H. Stillinger, Nature 393 (1998). [76] S. K. Das et al., Physical Review E-Statistical, Nonlinear and Soft Matter Physics 73 (2006). [77] M. Laradji, S. Toxvaerd, and O. G. Mouritsen, arXiv preprint cond-mat/9609177 (1996). [78] M. Ediger, Annual Review of Physical Chemistry 51 (2000). [79] C. Huang et al., Proceedings of the National Academy of Sciences 106 (2009). 134 [80] A. K. Soper, and M. A. Ricci, Phys Rev Lett 84 (2000). [81] F. Mallamace et al., Proceedings of the National Academy of Sciences 104 (2007). [82] L. Xu et al., Nature Physics 5 (2009). [83] S. Yip, Handbook of materials modeling (Springer, 2005). [84] H. Stanley et al., Journal of Physics: Condensed Matter 12 (2000). [85] A. Rahman, Physical Review 136 (1964). [86] A. Angell, Nature Materials 11 (2012). [87] A. Nakano et al., AIP Conference Proceedings 583 (2001). [88] K. Jayaraman, Journal of propulsion and power 25 (2009). [89] B. Dikici et al., Energy & Fuels 23 (2009). [90] Y. Q. Yang et al., Journal of Applied Physics 95 (2004). [91] M. A. Trunov et al., Combustion and Flame 140 (2005). [92] E. L. Dreizin, Combustion and Flame 105 (1996). [93] L. Liu et al., Phys Rev Lett 95 (2005). [94] K. Park et al., The Journal of Physical Chemistry B 109 (2005). [95] P. Grassberger, and I. Procaccia, The Journal of Chemical Physics 77 (1982). [96] C. Angell, Journal of Physics: Condensed Matter 12 (2000). [97] R. Kayser, and J. Hubbard, Phys Rev Lett 51 (1983). [98] F. Delyon, and B. Souillard, Phys Rev Lett 51 (1983). [99] J. D. Gezelter, E. Rabani, and B. Berne, The Journal of Chemical Physics 107 (1997). [100] E. Rabani, J. D. Gezelter, and B. Berne, The Journal of Chemical Physics 107 (1997). [101] V. K. de Souza, and D. J. Wales, The Journal of Chemical Physics 129 (2008). [102] C. A. Angell, Science 319 (2008). [103] M. Potuzak, R. C. Welch, and J. C. Mauro, The Journal of Chemical Physics 135 (2011). [104] M. Ediger, C. Angell, and S. R. Nagel, The Journal of Physical Chemistry 100 (1996). [105] E. Rabani, J. D. Gezelter, and B. Berne, Phys Rev Lett 82 (1999). [106] E. B. Moore, and V. Molinero, The Journal of Chemical Physics 130 (2009). [107] E. R. Weeks, and D. Weitz, Phys Rev Lett 89 (2002). [108] P. Kumar, S. V. Buldyrev, and H. E. Stanley, Proceedings of the National Academy of Sciences 106 (2009). [109] R. A. Yetter, G. A. Risha, and S. F. Son, Proceedings of the Combustion Institute 32 (2009). [110] K. Wikfeldt, A. Nilsson, and L. G. Pettersson, Physical Chemistry Chemical Physics 13 (2011). [111] S. Alavi, J. W. Mintmire, and D. L. Thompson, The Journal of Physical Chemistry B 109 (2004). [112] S. H. Khan et al., Phys Rev Lett 105 (2010). [113] P.-L. Chau, and A. Hardwick, Molecular Physics 93 (1998). [114] T. Campbell et al., Phys Rev Lett 82 (1999). [115] A. F. Voter, and S. P. Chen, in Mater. Res. Soc. Proc (Materials Research Society, 1987), p. 175. 135 [116] P. Vashishta et al., Journal of Applied Physics 103 (2008). [117] A. Shekhar et al., to be published (2013). [118] W. Lauterborn, and C.-D. Ohl, Ultrasonics sonochemistry 4 (1997). [119] S. S. Cook, Proceedings of Royal Society of London Series A 119 (1928). [120] A. Momber, J Mater Sci 38 (2003). [121] A. Momber, Composites Part B: Engineering 34 (2003). [122] T. Kodama, and Y. Tomita, Applied Physics B: Lasers and Optics 70 (2000). [123] D. Ranjan, J. Oakley, and R. Bonazza, Annual Review of Fluid Mechanics 43 (2011). [124] E. Johnsen, and T. Colonius, Journal of Fluid Mechanics, 629 (2009). [125] M. Vedadi et al., Phys Rev Lett 105 (2010). [126] A. C. Wright, J Non-Cryst Solids 159 (1993). [127] L. T. Zhuravlev, Colloids and Surfaces A: Physicochemical and Engineering Aspects 173 (2000). [128] W. Wang et al., to be submitted. [129] S. Susman et al., Phys Rev B 43 (1991). [130] P. A. V. Johnson, A. C. Wright, and R. N. Sinclair, J Non-Cryst Solids 58 (1983). [131] Q. Wang et al., J Non-Cryst Solids 143 (1992). [132] K. Nomura et al., Appl Phys Lett 99 (2011). [133] K. Miyake, N. Satomi, and S. Sasaki, Appl Phys Lett 89 (2006). [134] A. H. Narten, and H. A. Levy, The Journal of Chemical Physics 55 (1971). [135] T. Bakos, S. N. Rashkeev, and S. T. Pantelides, Phys Rev Lett 88 (2002). [136] J. Sauer, The Journal of Physical Chemistry 91 (1987). [137] A. A. Hassanali, and S. J. Singer, The Journal of Physical Chemistry B 111 (2007). [138] T. Kodama, and K. Takayama, Ultrasound in Medicine & Biology 24 (1998). [139] O. Edholm, and C. Blomberg, Chemical Physics 252 (2000). [140] C. Hyeon, R. I. Dima, and D. Thirumalai, The Journal of Chemical Physics 125 (2006). [141] R. I. Saye, and J. A. Sethian, Proceedings of the National Academy of Sciences 108 (2011). [142] M. Donsker, and S. Varadhan, Communications on Pure and Applied Mathematics 32 (1979).

Abstract (if available)

Abstract Nanotechnology is becoming increasingly important with the continuing advances in experimental techniques. As researchers around the world are trying to expand the current understanding of the behavior of materials at the atomistic scale, the limited resolution of equipment, both in terms of time and space, act as roadblocks to a comprehensive study. Numerical methods, in general and molecular dynamics, in particular act as able compliment to the experiments in our quest for understanding material behavior. In this research work, large scale molecular dynamics simulations to gain insight into the mechano-chemical behavior under extreme conditions of a variety of systems with many real world applications. The body of this work is divided into three parts, each covering a particular system: ❧ 1) Aggregates of aluminum nanoparticles are good solid fuel due to high flame propagation rates. Multi-million atom molecular dynamics simulations reveal the mechanism underlying higher reaction rate in a chain of aluminum nanoparticles as compared to an isolated nanoparticle. This is due to the penetration of hot atoms from reacting nanoparticles to an adjacent, unreacted nanoparticle, which brings in external heat and initiates exothermic oxidation reactions. ❧ 2) Cavitation bubbles readily occur in fluids subjected to rapid changes in pressure. We use billion-atom reactive molecular dynamics simulations on a 163,840-processor BlueGene/P supercomputer to investigate chemical and mechanical damages caused by shock-induced collapse of nanobubbles in water near amorphous silica. Collapse of an empty nanobubble generates high-speed nanojet, resulting in the formation of a pit on the surface. The pit contains a large number of silanol groups and its volume is found to be directly proportional to the volume of the nanobubble. The gas-filled bubbles undergo partial collapse and consequently the damage on the silica surface is mitigated. ❧ 3) The structure and dynamics of water confined in nanoporous silica are different from that of bulk water, and insight into the properties of confined water is important for our understanding of many geological and biological processes. Nanoporous silica has a wide range of technological applications because it is easy to tune the size of pores and their morphologies and to functionalize pore surfaces with a variety of molecular moieties. Nanoporous silica is used in catalysis, chromatography, anticorrosion coatings, desalination membranes, and as drug delivery vehicles. We use reactive molecular dynamics to study the structure and dynamics of nanoconfined water between 100 and 300 K.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Thermal properties of silicon carbide and combustion mechanisms of aluminum nanoparticle

PDF

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

PDF

Dynamics of wing cracks and nanoscale damage in silica glass

PDF

Heat-initiated oxidation of aluminum nanoparticles

PDF

Molecular dynamics study of energetic materials

$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

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

PDF

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

PDF

Molecular dynamics simulations of nanostructures

PDF

Large-scale molecular dynamics simulations of nano-structured materials

PDF

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

PDF

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

PDF

Dynamics of water in nanotubes: liquid below freezing point and ice-like near boiling point

PDF

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

PDF

Multi-scale quantum dynamics and machine learning for next generation energy applications

PDF

Nanoindentation of silicon carbide and sulfur-induced embrittlement of nickel

PDF

Physics informed neural networks and electrostrictive cavitation in water

PDF

Interactions of planetary surfaces with space environments and their effects on volatile formation and transport: atomic scale simulations

PDF

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

PDF

Molecular dynamics studies of protein aggregation in unbounded and confined media

Asset Metadata

Creator Shekhar, Adarsh (author)

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

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Materials Science

Publication Date 08/19/2013

Defense Date 07/26/2013

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

Tag aluminum nanoparticles,high performance computing,large scale,molecular dynamics,nanomaterials,OAI-PMH Harvest,oxidation,shock,silica,simulations,Water

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Vashishta, Priya (committee chair), Kalia, Rajiv K. (committee member), Nakano, Aiichiro (committee member), Nutt, Steven R. (committee member)

Creator Email adarsh.shekhar@gmail.com,ashekhar@usc.edu

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

Unique identifier UC11287911

Identifier etd-ShekharAda-2011.pdf (filename),usctheses-c3-321802 (legacy record id)

Legacy Identifier etd-ShekharAda-2011.pdf

Dmrecord 321802

Document Type Dissertation

Rights Shekhar, Adarsh

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