Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Multimillion-to-billion atom molecular dynamics simulations of deformation, damage, nanoindentation, and fracture in silica glass and energetic materials
(USC Thesis Other)
Multimillion-to-billion atom molecular dynamics simulations of deformation, damage, nanoindentation, and fracture in silica glass and energetic materials
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
1 MULTIMILLION-TO-BILLION ATOM MOLECULAR DYNAMICS SIMULATIONS OF DEFORMATION, DAMAGE, NANOINDENTATION, AND FRACTURE IN SILICA GLASS AND ENERGETIC MATERIALS by Yi-Chun Chen A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (PHYSICS) August 2008 Copyright 2008 Yi-Chun Chen ii ACKNOWLEDGEMENTS I would first like to sincerely thank my Ph.D. advisor, Dr. Rajiv K. Kalia for his continuous guidance, support and encouragement. I greatly appreciate the valuable advice about both research and life given by Drs. Priya Vashishta and Aiichiro Nakano. I thank Drs. Nelson Eugene Bickers and Jia Grace Lu from the Department of Physics and Astronomy for serving on my thesis committee. I am grateful for many fruitful discussions and collaborations on studies of silica glass from Drs. Bouchaud and Rountree at Saclay and Dr. Yang at LLNL. I also would like to thank Drs. Goddard and van Duin at Caltech for providing the ReasFF potentials. It has been a wonderful and rewarding experience to be a member of CACS (Collaboratory for Advanced Computing and Simulations) for the past four years. I would like to thank my colleagues in CACS for their help and inspirational conversations. I deeply appreciate Dr. Ashish Sharma, Dr. Cheng Zhang, Dr. Zhen Lu, Dr. Kenichi Nomura, Richard Clark, Hsiu-Pin Chen, Richard Seymour, and Weiqiang Wang for their interactions and collaborations. I also would like to acknowledge Patricia Wong for her help in CACS. Finally, I would like to thank my dear family and friends. I would like to express my greatest appreciation to my parents, my brother and sisters for their love. Their support and encouragement have been an inestimable treasure in my life. iii TABLE OF CONTENTS ACKNOWLEDGEMENTS.............................................................................................ii TABLE OF CONTENTS ...............................................................................................iii LIST OF FIGURES.........................................................................................................v ABSTRACT ...................................................................................................................x CHAPTER 1: INTRODUCTION....................................................................................1 1.1 Glasses ...............................................................................................................1 1.2 Energetic materials .............................................................................................4 1.3 Deformation .......................................................................................................6 1.4 Nanoindentation .................................................................................................8 1.5 Fracture ............................................................................................................ 11 1.6 Multiscale simulations of deformation and fracture........................................... 15 CHAPTER 2: SIMULATION METHODOLOGY AND PARALLEL IMPLEMENTATION ................................................................................................... 20 2.1 Molecular dynamics simulation background ..................................................... 20 2.2 Interatomic potentials ....................................................................................... 24 2.2.1 Effective force field.......................................................................................... 25 2.2.2 Reactive force field........................................................................................... 27 2.2.2.1 Bond order and bond energy .......................................................................... 28 2.2.2.2 Lone pair energy............................................................................................ 30 2.2.2.2 Overcoordination and undercoordination energy............................................ 31 2.2.2.3 Valence angle energies .................................................................................. 31 2.2.2.4 Torsion angle energies ................................................................................... 33 2.2.2.5 Hydrogen bond energy................................................................................... 34 2.2.2.6 van der Waals and Coulomb energies............................................................. 34 2.3 Integration algorithms....................................................................................... 35 2.4 Physical properties............................................................................................ 36 2.4.1 Thermo-mechanical properties.......................................................................... 36 2.4.2 Structural properties.......................................................................................... 37 2.4.3 Dynamic properties........................................................................................... 40 2.5 Parallel algorithms............................................................................................ 41 2.6 Performance of parallel MD algorithm.............................................................. 44 CHAPTER 3: NANOCAVITATION DAMAGE IN SILICA GLASS........................... 48 3.1 Tensile and compressive fracture simulation..................................................... 49 3.2 Void coalescence under tension ........................................................................ 52 CHAPTER 4: SHEAR DEFORMATION IN SILICA AND RDX................................. 63 4.1 Shear deformation in silica ............................................................................... 63 iv 4.2 Shear deformation in RDX ............................................................................... 75 CHAPTER 5: NANOINDENTATION IN SILICA AND RDX..................................... 81 5.1 Nanoindentation of silica glass ......................................................................... 81 5.2 Nanoindentation of RDX .................................................................................. 91 CHAPTER 6: CONCLUSION ...................................................................................... 99 REFERENCES ........................................................................................................... 102 v LIST OF FIGURES Figure 1.1: Schematic diagram of temperature vs volume for crystalline and noncrystalline materials. Crystallization occurs at the melting temperature T m for crystalline materials. However, a characteristic of noncrystalline materials is the existence of glass transition temperature T g . ............................................................2 Figure 1.2(a): Schematic two-dimensional representation of an ordered SiO 2 network. The fourth bonds of Si points upward or downward from the plane of view. (b): Schematic two-dimensional representation of a random SiO 2 network. The fourth bonds of Si points upward or downward from the plane of view. .............................3 Figure 1.3: Structures of energetic materials having similar NO 2 groups, (a) RDX, (b) TATB, (c) TNT, and (d) PETN................................................................................5 Figure 1.4: An illustration of stress-strain behavior showing the elastic to plastic deformation. ............................................................................................................7 Figure 1.5: Modeled RDX structure when <100> shear is applied. The diagonal trace across the box is the {021} slip plane. .....................................................................8 Figure 1.6: The three modes of fracture. (a) Mode I, opening mode, tensile stress is applied normal to the direction of crack propagation. (b) Mode II, sliding mode, in- plane shear stress is applied parallel to the direction of crack propagation. (c) Mode III, tearing mode, out-of-plane shear stress is applied normal to the direction of crack propagation. ................................................................................................. 13 Figure 1.7: AFM images of silica surface, showing coalescence of cavities is shown..... 14 Figure 1.8: (a) AFM image of shocked RDX crystal. Cracks run across or along the boundaries. (b) AFM image of the TNT columnar grains in composite B containing river line patterns which are a few nanometers in depth. ........................................ 15 Figure 1.9: Snapshots from a fracture simulation of a pre-notched PMMA plate under unaxial tensile strain. ............................................................................................. 16 Figure 2.1: Experimental measurements of pair distribution functions of Si at different temperatures by Kimura et al.[66]. ........................................................................ 39 Figure 2.2: A schematic 2-D view of spatial decomposition used for parallel MD implementation. The internode communication is necessary for force calculation for atoms on other nodes and also for atom migration when the updated position is outside the local boundary. .................................................................................... 42 vi Figure 2.3: The benchmark result for silica on a Cray T3E cluster. The wall-clock time (open circles) and inter-processor communication time (open squares) per MD step with isogranular scaling are presented. The number of atoms per processor remains constant = 648,000. The result shows communication time is negligible compared to the computation time. ............................................................................................ 46 Figure 2.4: The isogranular scaling benchmark results for RDX on (A) BlueGene/L at Lawrence Livermore National Laboratory and (B) Altix 3000 at NASA-Ames. The grain sizes are 10,752 on the BlueGene/L and 36,288 on the Altix......................... 47 Figure 3.1: (a) Initial configuration of the system under uniaxial strain. The frozen layers are shown as dark gray regions. The strain is applied on the top and bottom layers by displacing boundary-layer atoms. (b) Snapshot shows cavities in front of the crack, and their coalescence leads to failure. .................................................................... 50 Figure 3.2: (a) Illustration of the wing crack simulation setup. The gray plate represents a rigid indenter and the blue region is the pre-cracked a-SiO 2 system. The indenter speed is 150 m/s. (b) Different nanocavities are shown in different colors. At 19 ps, nanocavities grow and coalesce into nanocrack columns. At 21 ps, nanocolumns continue to grow and form a wing crack. ............................................................... 51 Figure 3.3: (a) Snapshot of voids and nanocavities at a strain rate of 10 8 sec -1 . (b) Strain dependence of the average porosity per void in the two-void system, φ 2 , relative to the porosity of the single-void system φ 1 . (c) Snapshot of voids and the fractured inter-void ligament at ε = 8%. (d) The single void system at ε = 8% shows little change in the void size........................................................................................... 56 Figure 3.4: (a) Si-O pair distribution functions and (b) Si-O-Si bond angle distributions in the middle of the inter-void ligament in the unstrained and strained systems.......... 57 Figure 3.5: Bond switching mechanism involving non-bridging oxygens is shown in (a) and (b) by red and white dashed lines between blue and white atoms. In (c), the white dashed line indicates bond switching between white and green oxygens. (d) shows the Si-O-Si bond angle distribution for the rings involved in the bond- switching events. Red and yellow spheres represent Si and O atoms, respectively. 59 Figure 3.6: Voids in a billion atom system under a strain rate of 10 9 sec -1 . (a) shows a slice of the entire system; (b) shows cracks on void surfaces at ε=9%; (c) and (d) show void coalescence and inter-void ligament failure at ε=10.5% and ε=12%, respectively. .......................................................................................................... 61 Figure 3.7: Porosity vs. strain in billion-atom systems containing 500 identical voids distributed regularly and randomly. (a) and (b) show volumes of coalesced cavities in the latter system at ε=8% and 10%, respectively. ............................................... 62 vii Figure 4.1: (a) Shear deformation of a spherical void of diameter 10 nm into an ellipsoid after 15 ps; (b) time variation of the deformation parameter, D, for voids of initial diameters 3 nm (blue circles), 10 nm (red squares), and 50 nm (green triangles). The applied shear rate is 10 12 sec -1 ................................................................................ 67 Figure 4.2: Panels (a) and (b) show plastic deformations of voids of initial diameters 3 nm and 50 nm, respectively. Panel (c) shows a plastically deformed void which was initially a sphere of diameter 10 nm; and panel (d) shows its breakup at 40 ps. ...... 68 Figure 4.3: Snapshots of voids in elastic and plastic regimes. Panels (a) and (c) show deformation of voids under strain at 15 ps and 30 ps, respectively. We switch off the strain and relax the systems for 5 ps. Panels (b) and (d) show the deformation of those voids after 5ps. In panel (b) the nanocracks have healed, and in panel (d) they have not................................................................................................................. 69 Figure 4.4: Point defects responsible for shear deformation and flow, and crack nucleation in silica glass. (a) A 6-membered –Si-O-Si-O– ring (green) adjacent to a 9-membered ring (blue) with a fully coordinated O (magenta) and a 3-fold coordinated Si defect (blue). (b) At a strain of 8%, the magenta O atom becomes non-bridging and the blue Si atom remains under-coordinated. These point defects reside on a 13-membered ring (grey). (c) At a strain of 10%, the green Si and magenta O atoms become fully coordinated by bonding with each other, but the blue Si is still under-coordinated. The yellow region is inside a 14-membered ring. ..... 71 Figure 4.5: (a) Density of Si and O defects in the damage zone as a function of time; and (b) Si-O pair-distribution function for these defects. The initial void diameter was 3 nm......................................................................................................................... 73 Figure 4.6: (a) Mean square displacements for defects at strains of 15% and 55%. (b) Velocity auto-correlation functions (VAF) for defects at a strain of 15%. The VAF in the shear direction (Y axis) is close to the VAF in the plane perpendicular to the shear direction. ...................................................................................................... 73 Figure 4.7: Panels (a) and (b) are snapshots of voids under shear strains of 20% and 35%, respectively, at 1200 K. In panel (c), the void is under a strain of 35%. Here the strain rate is 2!10 12 sec -1 . In panel (d), the strain is 35% and the strain rate is 5!10 11 sec -1 . At a slower strain rate, the void breaks up into smaller voids. ....................... 74 Figure 4.8: (a) Setup of shear deformation simulations in RDX. (b) Time variation of the deformation parameter, D, for a void of initial diameter 3 nm. The applied shear rate is 10 12 sec -1 . .................................................................................................... 76 viii Figure 4.9: Snapshots of the center layer of RDX at different simulation times. The initial void diameter is 6 nm. (a) Initially, the molecules move with uniform displacement. (b) At a strain of 8%, the molecules from the top and bottom start to move more dramatically than those in the center. (c) The maximum displacement of molecules is found along 45° to the shear direction (strain = 10%). (d) At a strain of 12%, molecules around the void move inward. The color represents the molecular center- of-mass displacement. ........................................................................................... 77 Figure 4.10: Snapshots of temperature distribution around the void. (a) At a strain of 1%, the temperature of a few molecules around the void is about 400 K. (b), (c) and (d) show temperature distributions at strains of 7%, 13% and 19%, respectively. The void is almost filled with high-speed molecules at a strain of 13%......................... 79 Figure 4.11: The time variation of density of RDX molecules whose dipole moment differences are greater than 0.5au........................................................................... 80 Figure 5.1: (a) Pileup around the indenter at the maximum indentation depth. Indentation causes localized damage creating a high-temperature area underneath the indenter, which facilitates pile-up through plastic flow. The colors represent temperature. (b) Load versus displacement during loading and unloading phases of indentation in simulation 1........................................................................................................... 83 Figure 5.2: Density distribution at the maximum indentation-depth (a), and after unloading (b). For clarity, the indenter is not shown. Silica substrate deforms until the local density under the indenter reaches 2.6 g/cc, which prevents downward material flow. Significant pileup with relatively low density develops at the maximum indentation depth. The densified region relaxes after unloading............ 85 Figure 5.3 (a)-(c): Atomic configurations showing transport of oxygen defects; (a) the initial state, (b) the transition state, and (c) the final state. Atoms involved in the event are labeled as O1 (initial defect oxygen atom), Si2 (bond switching silicon) and O3 (final defect oxygen). The oxygen defect migrates by switching bonds from Si2-O3 to O1-Si2. (d)-(f): Atomic configurations showing an annihilation event for a pair of defect atoms; (d) the initial state, (e) the transition state, and (f) the final state. Atoms involved in the event are labeled as Si1 (initial silicon defect), O2 (bond switching oxygen), Si3 (bond switching silicon) and O4 (initial oxygen defect). Si1 and O4 defects annihilate by switching bonds from O2-Si3 to Si1-O2 and Si3-O4. 88 Figure 5.4: Snapshots of defect migration pathways ((a) and (b)) in a ball-and-stick representation. Yellow and red spheres represent silicon and oxygen atoms. Black dotted lines indicate covalent bonds before defect migration. Blue and red arrows represent initial and final positions of defect atoms, respectively. (a) Migration of a silicon defect. (b) Annihilation of a pair of silicon and oxygen defects. Panels (c) and (d) show the length of the string-like defect motion as a function of the number of migration events during loading and unloading, respectively. Mechanical loading by the indenter induces longer migration pathways................................................ 90 ix Figure 5.5: Temperature distribution of the RDX substrate at different indentation depths. (a) At 5.2 Å, no heating is observed. (b) At 10.45 Å, molecules around the indenter begin to heat up. (c) At 15.7 Å, the “temperature” of some of the molecules exceeds 500K. At indentation depths of (d) 20.95 Å, (e) 26.20 Å and (f) 31.6 Å, we observe high “temperature” of RDX molecules within ~ 10 Å of the indenter contact area. 93 Figure 5.6: Indentation size-dependence of hardness value. ........................................... 94 Figure 5.7: The center-of-mass displacement, com r ! , of RDX molecules at various indentation depths. For clarity, only half of the system is shown. Panels (a) and (b) are snapshots at indentation depths of 10.45 Å and 25.15 Å, respectively. ............. 95 Figure 5.8: Panel (a) shows that the indentation damage is concentrated within the yellow dashed region. Panel (b) is an illustration of how molecules in the yellow dashed region are displaced in panel (a). Each sphere represents an RDX molecule, where the molecules move in the direction of ! 1 20 [ ] . ....................................................... 96 Figure 5.9: (a) and (b) show fragmented RDX molecules on the indenter surface at time 14.5 and 19.65 ps. Red, dark blue, light blue and white spheres represent oxygen, nitrogen, carbon and hydrogen atoms, respectively. (c): The mean square displacement of a few fragmented RDX molecules as a function of time. The dynamics of these molecules can be decomposed into translational and rotational motions. ................................................................................................................ 98 x ABSTRACT Multimillion-to-billion molecular dynamics (MD) simulations are carried out to study atomistic mechanisms of deformation, damage and failure in silica glass and energetic materials. The simulations are based on experimentally validated interatomic potentials and employ highly efficiently algorithms for parallel architectures. The onset of void-void interaction is investigated by performing MD simulations of amorphous silica under hydrostatic tension. The simulations reveal that nanocavities in amorphous silica (a-SiO 2 ), which are linked to Si-O rings, play an important role in void-void coalescence and inter-void ligament failure. Nanocracks nucleated by the migration of three-fold coordinated Si and non- bridging O on —Si-O-Si-O— rings are observed in the multimillion MD simulations of a single void in amorphous silica subjected to a high shear rate. With the increase in shear strain, nanocracks appear on void surfaces and the voids deform into a threadlike structure. At a strain of 40%, the voids break into fragments. The results are similar to experimental and theoretical studies of bubble deformation and breakup under shear. Defects such as voids are known to be important in the detonation of energetic materials. To investigate deformation of a void in an RDX crystal under high shear rate, we have performed million-atom reactive force field (ReaxFF) MD simulations. Simulations reveal that without breaking a bond, the excess strain energy leads to translational and rotational motion of RDX molecules. At a strain of 13%, molecules with high kinetic energy collapse inward without affecting the rest of the system. xi MD simulations of nanoindentation in amorphous silica reveal migration of defects and their recombination in the densified plastic region under and the material pileup region around the indenter. The plastic flow of silica glass is related to the defect transport mechanism where a defect migrates a considerable distance via a chain of bond-switching events[44]. We obtained a hardness value of 7.2 GPa using a sharp indenter and 8.0 GPa for a slightly blunt indenter. We have also performed nanoindentation simulation on a (100) α-RDX crystal surface using ReaxFF. Simulation reveals localized melting and decomposition of RDX molecular fragments. We have found a distinct (210) plane boundary, where molecules above the ! 210 ( ) plane have displaced dramatically and molecules below the plane remain intact. Simulation also shows the fragmented RDX molecules diffuse from the substrate and walk on the indenter surface. 1 CHAPTER 1: INTRODUCTION This chapter serves as an introduction to the dissertation. The basic characteristics of glasses and energetic materials are discussed first, followed by brief descriptions of deformations, nanoindentation and fracture. This chapter concludes with an overview of various simulation models and their applications to glasses and energetic materials. 1.1 Glasses Glassy materials are ubiquitous in the modern society. Applications of glasses range from containers, lenses, and optical waveguides to electronic displays. Over the years, scientists have tried to answer the question “what is glass?” Tammann[117] suggested that: “In the glassy state, there are solid, uncrystallized materials.” A more precise definition was mentioned in 1945 at the American Society for Testing and Materials: “Glass is an inorganic product of fusion which has been cooled to a rigid condition without crystallizing.” However, the finding of organic glasses (such as polymenthyl methacrylate) and novel manufacturing methods, including vacuum evaporation, have made this definition no longer appropriate. Glassy, or noncrystalline, materials do not solidify in the same manner as crystalline materials. In the latter case, crystallization occurs at the freezing point, T m , with a discontinuous decrease in volume (see figure 1.1). As the temperature continues to drop, the volume decreases with a lower temperature coefficient. However, a glassy material becomes more viscous with decreasing temperature, thus the volume drops continuously upon cooling. At the glass transition temperature (also 2 called the transformation temperature), a slight deviation from the original slope is observed. Below this temperature, the liquid becomes a solid glass. Figure 1.1: Schematic diagram of temperature vs volume for crystalline and noncrystalline materials. Crystallization occurs at the melting temperature T m for crystalline materials. However, a characteristic of noncrystalline materials is the existence of glass transition temperature T g . Turning to the glass structure, Goldschmidt[51] first proposed the network hypothesis in 1926. Zachariasen[129] later pointed out that the energy differences between a glass and crystal of the same composition are very small, which implies the coordination number in glass must be approximately the same as in the crystal. For example, in crystalline SiO 2 , the silicon atom is surrounded by four oxygen atoms which form SiO 4 tetrahedra. Thus the structure in silica glass is expected to be composed of SiO 4 tetrahedra. Figure 1.2 (a) and (b) show how SiO 2 may form either a crystal or a glass with SiO 4 tetrahedra. The crystal has a perfect regular structure whereas the glass has distorted bond angles which give a random effect and lead to the loss of any long-range periodicity. Melt Supercooled melt Glass Crystallization Crystal T g T m Temperature Volume 3 Figure 1.2(a): Schematic two-dimensional representation of an ordered SiO 2 network. The fourth bonds of Si points upward or downward from the plane of view. (b): Schematic two-dimensional representation of a random SiO 2 network. The fourth bonds of Si points upward or downward from the plane of view. In addition, Zachariasen proposed four conditions for oxide glass formation: (1) An oxygen atom may not be linked to more than two cations. (2) The coordination number of the cation must be small. (3) The oxygen polyhedra share corners with each other, not edges or faces. (4) At least three corners in each oxygen polyhedron must be shared. Oxides of the type R 2 O 3 , RO 2 and R 2 O 5 satisfy these conditions, which are confirmed by the existence of B 2 O 3 , SiO 2 , GeO 2 , As 2 O 3 and P 2 O 5 . The formation of glass is also possible with systems of several components. Zachariasen also suggested that an oxide glass may be formed if (1) there is a dominant percentage of cations which are surrounded by oxygen triangles or tetrahedra, (2) which share only corners with each other, and (3) some oxygen atoms (b) (a) 4 are linked to only two such cations. Thus, an oxide glass must contain a large amount of the oxides of B, Si, P, Ge, As. Another important glass formation criterion is based on Sun’s finding[116], which relates to the concept of single bond strength. Single bond strength is defined as the dissociation energy (the energy needed to dissociate into gaseous atoms) of the cation with its oxygen neighbors, divided by the coordination number of the cation. He suggested that oxides with high bond strengths (>100 kcal) may form glass. Examples of high bond strength oxides of B, Si, Ge, P agrees with this observation. Some of the properties of glass have been researched thoroughly because of their commercial and industrial importance. Scientists have made numerous investigations of mechanical properties which include elastic moduli, strength and microhardness. Bartenev[21] has summarized the work in this field and these properties are discussed in several handbook publications (Uhlmann and Kreidl[123], Harper[57]). Further discussion can be found in the later sections of this chapter. 1.2 Energetic materials Materials that release a large amount of chemical energy are categorized as energetic materials. Pyrotechnic compositions (gunpowder, for example), explosives, fuels and propellants are categorized as energetic materials. These materials have played an important role in human history starting from the discovery of gunpowder to recent explosives and rocket fuels. Low explosives deflagrate and create a subsonic explosion (below 1000 m/s). On the other hand, high explosives undergo detonation with a speed between 1000 to 9000 m/s. Detonation is defined as an explosive effect with high-density mass flow. Conventionally, high explosives are divided into two 5 subgroups: primary and secondary explosives. Lead azide and lead styphnate are two commonly used primary explosives, as they are sensitive to mechanical load, friction, and heat. They are often used as a trigger to ignite secondary explosives. Secondary explosives are less sensitive, hence have more applications. Examples of secondary explosives are cyclotrimethylenetrintramine (RDX) and trinitrotoluence (TNT). Energetic materials consist of fuel and oxidizer that release a large amount of energy when ignited. The energy density depends on the ratio of fuel to oxidizer. Crystalline materials such as nitrates, nitro compounds, nitramines and perchlorates are used as oxidizers for the high bond energy characteristic. Nitrogen oxides, in particular, have applications ranging from C-nitro derivatives such as trinitrotoluence (TNT), to O-nitro compounds such as nitrocellulose, to N-nitro such as HMX and RDX. Figure 1.3 shows the molecular structures of some commonly used energetic materials. Figure 1.3: Structures of energetic materials having similar NO 2 groups, (a) RDX, (b) TATB, (c) TNT, and (d) PETN. In order to maximize the energy release and to avoid unexpected explosions from materials of this kind, recent research has focused on the structure and properties of energetic materials[77], dislocations[14, 15, 17], decomposition, combustion and (a) (b) (c) (d) 6 detonation chemistries[31]. Mechanical properties of some explosives will be discussed in later sections. 1.3 Deformation The change in a body’s shape and size when subjected to an applied force is referred to as deformation. When the stress σ is proportional to the strain ε, the deformation is elastic and the constant of proportionality E is the Young’s modulus, or modulus of elasticity The relationship is known as Hooke’s law, ! " E = . (1.1) Elastic deformation is not permanent. Upon releasing the applied load, the sample returns to its original undeformed state. If the strain along the direction of applied tensile stress is taken as ε z , and the strains from contraction in the directions perpendicular to the applied stress are ε x , ε y , we get Poisson’s ration ν from the relationship, z y z x - - ! ! ! ! " = = . (1.2) For isotropic materials, one can obtain shear modulus, G, through ( ) ! + = 1 2G E . (1.3) Magnesium, for instance, has a Young’s modulus of 45 GPa, shear modulus of 17 Gpa and Poisson’s ratio of 0.29. When the stress is no longer proportional to the strain, plastic deformation takes place. Figure 1.4 shows schematically stress versus strain curve. Here the plastic deformation occurs after the point P. At the nanometer scale, plastic deformation results from bond-breaking mechanism and formation of new bonds when molecules or atoms move away from their original positions. 7 Figure 1.4: An illustration of stress-strain behavior showing the elastic to plastic deformation. It is found that Hooke’s law is satisfied by glasses under normal conditions, which indicates that glassy materials are elastic and can be described by Young’s modulus, E. Fused silica has values of E = 70 GPa, G = 30.5 GPa, K = 37 GPa, where G is the shear modulus and K is bulk modulus. According to Deeg’s experiment[43], alkali oxide (Na 2 O and K 2 O) weakens the structure of silica glass. Thus Young’s modulus decreases as the content of alkali oxide increases. On the contrary, introducing oxide which has high field strength cation increases the Young’s modulus[20]. The composition dependence of shear modulus is similar to that of Young’s modulus. Only the values are lower, for most cases, between 20 to 40 GPa. The bulk modulus generally increases when other oxides are introduced into vitreous silica. However, it is shown by Weir and Shartsis[127] that when the cations themselves are deformable, with an increase in ionic radii, the bulk modulus decreases. The deformation mechanisms in energetic materials are also complex due to their molecular structures. In the presence of an applied load, the RDX molecules deform, Strain Stress Elastic Plastic P σ 8 shift and rotate. The chemical reactions take place in the heavily deformed regions. Figure 1.5 shows a slip plane {021}, which is found when the RDX crystal undergoes shear displacement[15]. The ab initio calculations by Kuklja et al.[67, 68] show molecular decomposition mechanism near the dislocation site. Reorientation of the molecular structure in an anthracen crystal can reduce the stacking fault energy as proposed by Ramdas et al.[104]. However, despite much research effort, the understanding of deformation at the atomistic level is limited. To study localized features in explosive materials, the second part of my thesis discusses the deformation of the secondary explosive RDX at the nanometer scale. Figure 1.5: Modeled RDX structure when <100> shear is applied. The diagonal trace across the box is the {021} slip plane. 1.4 Nanoindentation Indentation test is a commonly used technique to characterize the mechanical properties of materials. In an indentation test (nano- or micro-indentation), an indenter tip, such as diamond, is gradually pressed into the sample and removed afterward. The load on the indenter and the penetration depth are measured to obtain 9 the load-displacement curve, which yields information about plastic deformation and compliance. The resistance to plastic deformation due to mechanical load is hardness. Hardness is defined as the maximum load, F max , divided by the residual impression, A, namely, ! H = F max A . (1.4) Young’s modulus can be obtained from the relation developed by Pharr et al.[100], ! dF unload dh =S = 2 " E (1#$ 2 ) A , (1.5) where the slope of the unloading curve, dF unload /dh, is referred to as the stiffness of contact, E is the Young’s modulus, ν is the Poisson’s ratio, and A is the projected area. The penetration depth and the geometry of the indenter tip together give information about the contact area. The four-sided Vickers pyramid indenter has an angle of 136° between the opposite faces. The three-sided Berkovich indenter has a face angle of 65.3°, which is popular for its ease of manufacturing. The Knoop indenter has two different semi-angles, which is useful for testing brittle materials. Microhardness, the ability of a material to resist abrasion and scratching from mechanical loading, is often used to describe the strength of glassy materials. Research efforts in this area focus on better understanding of the underlying mechanisms. During indentation, elastic deformation occurs first followed by plastic deformation and damage which are measured after unloading. Fused silica has a 10 Vickers hardness of 6.18 Gpa. High-temperature nanoindentation tests of soda-lime glass and fused silica[22] show that these two glassy materials have different behavior. As the temperature increases, the modulus increases for fused silica and decreases for soda-lime glass; the degree of changes in hardness also differs. The Knoop microhardness measurements of silica glass in presence of different liquids[60] show a time-dependent load only in the presence of water, suggesting that water diffusion into the glass plays an important role during indentation. The indentation- induced fracture experiments[39] show different crack initiation points, where radial cracks form during the loading phase for crystalline materials and in the unloading phase of normal glasses (soda-lime, for example). The appearance of flow lines after loading was first described by Marsh[80] and Peter[99]. However, the flow mechanisms in the case of solid materials are linked to their crystal structures; glass, which has a disordered network, does not fall into this category. Plastic flow processes in glasses are not yet well understood. Indentation can also be used to obtain crack propagation data, by rapidly unloading the indenter. Cook and Liniger[40] found their perturbed indentation fracture system enabled the study of various crack velocities in one sample. For energetic materials, hardness test is a key experimental and modeling tool. The measured crystal surface energy γ from indentation experiments is often used to determine the cleavage susceptibility factor[16], 1/2 Gb ! " # $ % & ' , where G is the shear modulus and b is the Burgers vector. If the susceptibility is low, it indicates that the cleavage is relatively easy to form; on the other hand, high susceptibility factor 11 suggests that dislocation may nucleate easily. The reported values for cleavage susceptibilities are 0.066, 0.070 and 0.17 for RDX, PETN and MgO, respectively. Microindentation tests on RDX and MgO samples[13] support this hypothesis. Results show that RDX crystal has a hardness value of 245 Nmm -2 and is more crack- prone than MgO, but cracks do not propagate as easily in RDX as in MgO. Significant heating due to collapse of material pile-up is also reported. 1.5 Fracture The undesirable failure of engineering components makes fracture an important subject. There are two possible modes of fracture: ductile and brittle. In ductile fracture, substantial plastic deformation occurs before fracture. In contrast, brittle fracture exhibits little plastic deformation. Fracture process involves crack formation and propagation. Since ductile fracture has extensive plastic deformation, the crack length extends relatively slowly. In brittle fracture, crack grows relatively fast and crack propagation continues even without an increase in applied load when the load exceeds a critical value. Griffith developed one of the earliest criteria for crack propagation by performing an energy balance between the strain and surface energies[9, 53]. Consider a plate with a crack of length 2a subjected to tensile loading σ. The Griffith energy balance for a crack in equilibrium, can be described as 0 dA dW dA dU dA d s = + = ! , (1.6) where ψ is the total energy, U is the elastic energy, W s is the work required to form new surfaces, and A is the crack area. Applying Inglis’s stress analysis[62] gives 12 E a - U U 2 2 0 !" = , (1.7) where U 0 is the potential energy of a plate without a crack. Since W s can also be expressed in terms of surface energy, s s 4a W ! = . (1.8) Thus, crack propagates when the stress exceeds a critical value, σ c : 1/2 s c a E ! " # $ % & = ' ( ) 2 . (1.9) The stress analysis around the crack tip is another important topic. The stress distribution in front of a crack tip in an isotropic linear elastic material can be expressed as ( ) ( ) ! " ! # ) ( ) ( ) ( 0 2 , lim I ij I I ij r f r K r = $ where I is the mode of fracture, K is the stress intensity factor which depends on the geometry of the crack and mode of fracture, and f ij is a dimensionless function. There are three distinct modes by which a load can enable crack propagation. Mode I is an opening (tensile) mode. Modes II and III are sliding (in-plane shear) and tearing (out- of-plane) modes (see figure 1.6), respectively. 13 Figure 1.6: The three modes of fracture. (a) Mode I, opening mode, tensile stress is applied normal to the direction of crack propagation. (b) Mode II, sliding mode, in- plane shear stress is applied parallel to the direction of crack propagation. (c) Mode III, tearing mode, out-of-plane shear stress is applied normal to the direction of crack propagation. The resistance to breakage of glasses is very important. To manufacture stronger glasses, scientists have attempted to explain the cause[29, 69] for crack initiation and propagation in glasses. For decades, it was believed that under the influence of tensile stress, glass acts as a brittle solid with no observable plasticity. Thus, the brittleness of a glass was described in terms of its strength. The theoretical strength σ th of a glass is determined by the strength of the bonds between the individual components, 0 th a E! " = , (1.11) where γ is the surface energy and a 0 is the atomic spacing. For a typical silicate glass, a value of 35 GPa is obtained. However, most glasses have flaws, which make the practical strength several orders of magnitude smaller than the theoretical strength. Despite the efforts to understand brittleness of glasses, fundamental mechanisms of 14 fracture in glasses are still not well understand. My thesis will focus on fracture in amorphous silica to better understand atomistic mechanisms of damage observed during crack propagation and fracture. In a recent stress corrosion cracking experiment (SCC)[34], cavitation in front of a crack tip is shown to be important in the failure of fused silica. Figure 1.7 shows snapshots of AFM images of a glass sample under tensile loading. As time evolves, cavities grow and coalesce. The experiment provides an explanation of fracture mechanism, i.e., that crack propagation is mediated by cavities ahead of the crack tip. Figure 1.7: AFM images of silica surface, showing coalescence of cavities is shown. Sharma et al.[109] have used AFM as an imaging method to reveal planar and curvilinear cracks running throughout the shocked crystal RDX specimen, see figure 1.8(a). The submolecular steps are found with constant elevation ranging from 0.05 to 0.4nm. Hot spots with size greater than a molecule which can overcome this step height are thus necessary for initiation. Friction caused by the submolecular steps due to changes in the relative positions of molecules on the crack surfaces may initiate chemical reaction, which adds new insight to the requirement of decomposition processes of brittle molecular crystals like RDX. Lanzerotti et al.[72] applied a 15 similar technique to a composite B energetic material (59% of RDX, 40% of TNT and 1% of wax) which undergoes fracture across columnar grains of the TNT, as well as the river line patterns, see figure 1.8(b). Sharma et al. proposed the possibility of partial dislocations that penetrate the fracture surface and unit dislocation that is normal to the cleavage plane. They also suggested that the separation of molecules is related to the breaking of intramolecular bonds which causes the detonation of energetic materials. Figure 1.8: (a) AFM image of shocked RDX crystal. Cracks run across or along the boundaries. (b) AFM image of the TNT columnar grains in composite B containing river line patterns which are a few nanometers in depth. 1.6 Multiscale simulations of deformation and fracture Depending on the simulation model, system sizes and time scales have different constraints. Ab-initio quantum mechanical approaches, such as density-functional theory (DFT), are used to calculate the ground state of many-body systems. Significant research efforts are devoted to realize the O(N) scalable efficiency and extend the system size [28, 46, 47, 50, 110] in DFT based simulations. However, due 16 to the inherent complexity, the practical system size is generally limited to the order of a few hundred atoms[110] for picoseconds. The finite element method (FEM), on the contrary, deals with macroscopic systems and is often used to find an approximation for partial differential equations. Recent finite element simulation[11] on PMMA plate shows dynamic processes of fracture. The crack nucleates, propagates, branches and eventually the system fractures (see figure 1.9). The dynamic fracture experiment on a 229 mm×190 mm×9.5 mm sheet made of Homalite-100 and Loctite-384 are compared with simulation of the same setup and there is good agreement regarding the crack tip position/velocity history and crack initiation time. However, the resolution in FEM depends on the mesh size chosen for discretization. Figure 1.9: Snapshots from a fracture simulation of a pre-notched PMMA plate under unaxial tensile strain. Recently, nanoindentation has become a commonly used technique to study the mechanical properties of various materials. A multiscale nanoindentation simulation performed by Smith et al.[111] provides characteristics of phase transformation as well as load vs. displacement curves. They employ the extended local quasicontinuum model which combines either Stillinger-Weber empirical 17 potential[112] or tight-binding Hamiltonian[24] and finite element method. Close to the indenter, atomistic description is necessary due to the possible occurrence of structural phase transformations, dislocations or cavitations. On the other hand, coarser stress descriptions are used within the continuum model away from the indenter. The challenge of this multiscale technique rests on the smooth interphasing of different computing models at the boundaries. Smith et al. report reversible phase transformations and a reproducible hysteresis loop. These phases formed underneath the indenter are bct5, β-Sn, simple cubic and bcc. The load vs. displacement is close to the experimental values when the mesh size is small. Molecular dynamics (MD) simulation method has a wide variety of applications in materials science, biochemistry, physics, etc. The method has been used to study different materials under normal and extreme conditions. The method involves the solution of Newton’s equations of motion for interacting atoms/molecules. The atomic/molecular trajectories in the phase space are obtained by discretizing time and integrating Newon’s equations of motion. Comparison with experimental studies and first-principles quantum-mechanical calculations validates the interatomic/inter- molecular interaction potentials. Material properties of the system are calculated from atomic positions and momenta. There are limitations in the simulation time scale, e.g., material response such as fatigue is not practical for MD simulation. MD simulation has been used in this dissertation to study deformation, indentation and fracture in atomistic level. The first MD simulation was carried out by Rahman using 864 argon atoms to study structural and dynamical properties of liquid argon. The pair-correlation function and diffusion coefficient obtained from the 18 simulation agree with experimental data. The microcanonical ensemble used in Rahman’s simulation is also known as the NVE ensemble, which assumes constant number of particles, N, volume, V, and total energy, E. NVE ensemble represents a thermodynamically isolated system where the potential and kinetic energies may change but the total energy is conserved. Besides NVE ensemble, the canonical ensemble is also commonly used in MD simulations. It is sometimes referred to as constant temperature (NVT) ensemble. NVT ensemble describes an isolated system which exchanges energy with a heat bath (thermostat) to maintain a constant temperature. The isothermal-isobaric (NPT) ensemble is used to mimic experiments under ambient temperature and pressure. Other than thermostat, a barostat is also used to control the pressure by varying the system volume. In addition to equilibrium phenomena, MD simulation can be used to study non- equilibrium systems. Non-equilibrium MD (NEMD) algorithm is the response to an external perturbation. NEMD can provide properties such as viscosity and thermal conductivity. Recently, the NEMD method has become popular in indentation and fracture simulations. We have used this approach in our studies of indentation of energetic materials and amorphous silica. Recent advances have made massively parallel computers commercially available at a reasonable cost. Teraflop parallel machines along with the O(N) linear scaling algorithms make MD simulation sizes comparable to those used in experiments. It is now possible to perform a billion-atom MD simulation on a teraflops machine. 19 The calculation of forces at each time step is the most computationally intensive part of MD simulations. Straightforward force evaluation involves O(N 2 ) operations. The number of operation can be reduced to O(N) by using linked-cell lists approach for finite-range potentials and forces. In the linked-cell lists, the MD box is divided into cells whose length is slightly greater than the cut-off of the potential. In each of these cells, atoms are sorted and their corresponding information is stored in the linked list. To calculate the force on an atom, one simply needs to calculate the contributions from atoms in the same cell and atoms in the 26 nearest-neighbor cells. For long-range forces, such as Coulomb interactions, a divide-and-conquer technique such as the fast multipole method is used to achieve O(N) performance Multiple time steps (MTS)[115, 120] is often used to improve the efficiency of force calculation. MTS technique separates the force on an atom into a rapidly varying part and a slowly varying part. The rapidly varying part of the force is associated with the nearest neighbor atoms while the slowly varying part is for long- range interactions. The separation of force components leads to different calculation rates for these components. Rapidly varying part is computed every time step and the slowly varying part is calculated only after a few time steps. This thesis is arranged as follows: Chapter 2 contains a description of atomistic simulation methodology and implementation of MD simulations on parallel computers. In Chapter 3, nanocavitation damage in glass undergoing tensile and compressive fracture is discussed. Molecular dynamics simulations of shear deformation in silica glass and RDX are presented in Chapter 4. Chapter 5 describes nanoindentation simulations in silica and RDX followed by conclusions in Chapter 6. 20 CHAPTER 2: SIMULATION METHODOLOGY AND PARALLEL IMPLEMENTATION The first part of this chapter discusses the molecular dynamics (MD) simulation methodology. An overview of various statistical ensembles is presented and the interatomic potentials used in this thesis are discussed. Section 2.3 focuses on the integration algorithms used in our MD simulations. This is followed by a discussion of thermo-mechanical, structural, and dynamic properties. MD simulation of a system consisting of millions-to-billions of atoms is only feasible on a parallel computer. A large-scale system is divided into smaller subsystems which are mapped onto individual processors simultaneously. Section 2.5 focuses on a parallel implementation of an MD simulation code. The implementation is based on: (1) algorithms that provide high parallel efficiency, and (2) machine architectures that support high performance computing. This chapter concludes with the benchmark results for MD algorithms. 2.1 Molecular dynamics simulation background MD is a computational approach to simulate a system of interacting atoms by numerically solving Newton’s equation: i i i m F r = & & , (2.1) where ! m i is the mass, i r & & the acceleration, and i F is the total force acting on atom i. Atomic forces are derived from the total potential energy ! V , i N 1 i ) ,..., V( - r r r F ! ! = , (2.2) where N is the number of atoms. 21 In an MD simulation, the equation of motion are integrated and thus the atomic positions and momenta, i.e., the phase-space trajectories are obtained. MD simulations can be carried out in microcanonical and canonical ensembles. In the microcanonical ensemble, which is also referred to as the NVE ensemble, the number of particles, N, volume, V, and total energy, E, remain constant. In the canonical ensemble, N, V, and temperature, T, remain constant. In the NVE ensemble, the Hamilton’s equations of motion are, i i i i H H r p p r ! ! ! ! " = = & & , , (2.3) where p i are the particle momenta, and H is the Hamiltonian of the system, ! H = p i 2 2m i i=1 N " +V({r i }). (2.4) In the NVT ensemble first developed by Nosè[92, 93], the Hamiltonian includes an additional degree of freedom corresponding to a heat bath which keeps the temperature of the system constant: ! H = p i 2 2m i s 2 i=1 N " +V r i { } ( ) + p s 2 2Q +Gk B T ext ln(s). (2.5) The last two terms represent the effect of the heat bath. Here s is a “coordinate” for the heat bath, p s is the corresponding “momentum”, G=3N+1 is the number of degrees of freedom, k b is the Boltzmann constant, and Q is an effective “mass” of the heat bath. The variables s and Q rescale time and kinetic energy of the system to satisfy the constraint of canonical ensemble. Due to rescaling, the simulation time and the real time are related by ! dt'=dt/s. 22 Hoover[61] has modified Nosè’s scheme to eliminate the time rescaling problem. Martyna et al.[81] have further extended it to describe linked heat baths. The Hamiltonian for this extended Nosé-Hoover model with M linked heat baths is ! H = p i 2 2m i i=1 N " +V r i { } ( ) + p s i 2 2Q i i=1 M " +3Nk B T ext s 1 + k B i=2 M " Ts i , (2.6) which remains valid even for small and stiff systems. The derived dynamics can be expressed as, i i i m p r = & , ! ˙ p i ="# r i V r i { } ( ) " p i p s 1 Q 1 , i s i Q p s i = & , 2 1 2 2 1 1 3 Q p p T Nk m p p s s ext B N i i i s ! " # $ % & ' ! = ( = & , 1 1 2 1 1 + ! + ! ! " " # $ % % & ' ! = j s s ext B j s s Q p p T k Q p p j j j j & , ext B M s s T k Q p p M M ! = ! ! 1 2 1 & . (2.7) Andersen[8] developed an extended Lagrangian scheme to simulate systems under constant pressure. Parrinello and Rahman[95, 97] extended Andersen’s scheme to allow both the size and shape of the MD box to change when the system is subjected to an external stress. Parrinello and Rahman introduce an h-matrix, 23 ! H = h 1 ,h 2 ,h 3 ( ) = h 1x h 2x h 3x h 1y h 2y h 3y h 1z h 2z h 3z " # $ $ $ % & ' ' ' , (2.8) which is constructed from vectors of the three axes of the simulation box. The h- matrix relates the atomic coordinates, ! r i =(r ix ,r iy ,r iz ) , to the reduced coordinates, ! s i =(s ix ,s iy ,s iz ) , by, ! r i =s i1 h 1 +s i2 h 2 +s i3 h 3 =Hs i . (2.9) This technique can be combined with Nosè’s heat bath approach to simulate systems in the isothermal-isobaric (NPT) ensemble. Martyna et al.[82] proposed a flexible NPT system described by the Hamiltonian, ! H = p i 2 2m i i=1 N " +V r i { },H ( ) + p s 2 2Q + 1 2W g TrP #1 gP g [ ] + 3N +1 ( ) k B T req ln(s) +P ext detH [ ] , (2.10) where P g and W g are the box momentum and mass; s, p s , Q and T req are respectively the position, momentum, mass and temperature associated with the thermostat; and P ext is the external pressures. The equations of motions are, i g g i i i W m r P p r + = & i s i g g i g g i i Q p W Tr N W p p P p P F p ! ! ! = ] [ 3 1 & g g W H P H = & 24 g s N i i i ext g Q p m N P P V P I p I P ! " " # $ % % & ' + ! = ( =1 2 int 3 1 ) ( & Q p s s = & ! = + " + = N i B 2 g g g i i s T k d N Tr W m 1 2 ) 3 ( ] ~ [ 1 P P p p & , (2.11) where I is the identity matrix and P int is the internal pressure. Also, when d=1, the equations are reduced to the case of uniform scaling. 2.2 Interatomic potentials Interatomic potentials are the essential ingredients of materials modeling at the atomistic scale. The parameters in an interatomic potential are generally fitted to experimental data such as the bulk modulus, lattice constants, and melting temperature. The interatomic potential may be formulated as a sum of one-body, two- body, three-body terms etc., ! ! ! + + + = N k j i k j i N j i j i N i i V V V V , , 3 , 2 1 .... ) , , ( ! 3 1 ) , ( ! 2 1 ) ( r r r r r r (2.12) The one-body potential is the potential energy of atoms in an external field. Examples of two-body potential include Coulumb and van der Waals interactions. Three-body potential is due to the covalent effects. One of the most widely used models is the Lennard-Jones potential, V(r ij ) =4! " r ij # $ % % & ' ( ( 12 ) " r ij # $ % % & ' ( ( 6 * + , , - . / / , (2.13) 25 where the variables σ and ε are related to the length and energy scales intrinsic to the system. This potential is well-suited for systems composed of inert gas atoms such as He, Ne, Ar, Kr, Xe, etc.. The most widely used three-body potential is developed by Stillinger and Weber: 2 ) cos (cos exp ijk ijk ik ij ijk Weber Stillinger a r C a r C B V ! ! " # # $ % & & ' ( " + " = " (2.14) where a is the cutoff distance, B ijk and C are constants, and ijk ! is the optimal angle between atoms i, j and k. For system such as SiO 2 and GaAs, ijk ! =109 0 . The term ! cos" ijk -cos"ijk ( ) represents the effects of bond bending and the exponential terms represent bond stretching. To describe various types of chemical bond, Abell[2] proposed an analytic form for the binding energy. Tersoff[119] adopted this idea and developed bond-order potentials which relate the strength of the covalent bond to the coordination number. This model is suitable for covalently-bonded systems in which single, double and triple bonds occur. 2.2.1 Effective force field The interatomic potential used for silica in this thesis consists of two-body and three-body terms[126]: ! V = V ij (2) (r ij ) i< j " + V ijk (3) (r ij ,r jk ) i<k " j " . (2.15) The two-body potential incorporates steric repulsion between atoms, screened Coulomb interaction due to charge transfer between Si and O, and charge-dipole interactions. The functional form of the two-body potential is: 26 ( ) s ij s ij ij r r ij i j j i r r ij j i ij ij ij ij (2) e r Z Z e r Z Z r H r V 4 1 / 4 2 2 / 2 1 ) ( ! ! + ! + = " " # (2.16) where r ij is the separation between atoms i and j, H ij and η ij are the strengths and exponents of the steric repulsion, respectively, Z i is the effective charge and α i is the ionic polarizability. r 1s and r 4s are the decay lengths for the Coulomb and charge- dipole interactions. These parameters are determined by fitting to experimental data on structural and mechanical properties of amorphous silica. The three-body term has the form of Stillinger-Weber potential, see Eq. 2.14. Only Si-O-Si and O-Si-O three-body interactions are taken into account. The parameters used are given in Table 2.1. 27 Table 2.1: Parameters for SiO 2 potential. Z Si (e - ) 1.2 Z O (e - ) -0.6 α Si (Å 3 ) 0 α O (Å 3 ) 2.4 r 1s (Å) 4.43 r 4s (Å) 2.5 η Si-Si 11 η Si-O 9 η O-O 7 Η Si-Si 0.057 Η Si-O 11.387 Η O-O 51.692 B Si-O-Si (erg) 3.2×10 -11 B O-Si-O (erg) 0.8×10 -11 C Si-O-Si 1.0 C O-Si-O 1.0 ! " Si-O-Si 141.0 ! " O-Si-O 109.47 r 0 (Å) 2.6 r c (Å) 2.2.2 Reactive force field MD simulations of RDX crystal are carried out with a Reactive force field (ReaxFF[114, 124]) which describes large-scale reactive systems. Compared to MM3[5-7], UFF[105] and DREIDING[84] force field methods, ReaxFF can accurately describe chemical reactivity, such as bond formation and bond breaking, and dynamical interatomic charge transfer. Brenner’s potential[30] can also describe bond breaking, but it does not include van der Waals and Coulomb interactions. ReaxFF is much faster than any electronic structure calculation, which makes it 28 practical for simulations of large chemical systems. Nitramines potential is introduced into ReaxFF [114] using QM calculations designed for atomistic interactions under various environments. ReaxFF uses a general relationship between the bond distance and bond order and between bond order and bond energy to describe dissociation and formation of bonds between atoms. Valence angle and torsion angle terms in the force field are also defined in terms of bond orders and these dependencies ensure that all these terms vanish smoothly as bonds break. ReaxFF also includes non-bonded Coulomb and van der Waals interactions which are shielded at short range and vanish as the distance between a pair of atoms goes to zero. To avoid abrupt energy increase as interatomic distance goes to zero, the non-bonded interactions are smeared at short range. The energy has the following form: ! E = E bond + E lp + E over + E under + E val + E pen + E coa + E tors + E conj + E hbond + E vdWaals + E Coulomb (2.17) where E bond denotes the bond energy, E lp is the lone pair energy, E over is the overcoordination energy, E under is the undercoordination energy, E val is the valence angle energy, E pen is the penalty function, E coa is the three-body conjugation energy, E tors is the torsion angle energy, E conj is the four-body conjugation energy, and E hbond is the hydrogen bond interaction energy. E vdWaals and E Coulumb are the van der Waals and Coulumb energies, respectively. 2.2.2.1 Bond order and bond energy The basic quantity in ReaxFF, the bond order, BO ij , between atoms, depends on the interatomic distance r ij . Contributions to BO ij are separated into σ, π and double-π bonds. 29 ! ! " # $ $ % & ' ' ( ) * * + , + ! ! " # $ $ % & ' ' ( ) * * + , + ! ! " # $ $ % & ' ' ( ) * * + , = - + - + - = - bo6 o ij bo5 bo4 o ij bo3 bo2 o ij bo1 ij ij ij ij p r r p p r r p p r r p O B O B O B O B .. . / .. . / exp exp exp , (2.18) where ! r ij =|r ij |=|r i "r j | is the distance between atoms i and j. This uncorrected bond order, BO′ ij is used to obtain the uncorrected overcoordination term, Δ′ i , ! " i ' =#Val i + BO ij ' j$n(i) % . (2.19) One can construct the corrected bond order term, BO ij , using: ! BO ij " =BO ij '" # f 1 ($ i ' ,$ j ' )# f 4 ($ i ' ,BO ij ' )# f 5 ($ j ' ,BO ij ' ), ! BO ij " =BO ij '" # f 1 ($ i ' ,$ j ' )# f 1 ($ i ' ,$ j ' )# f 4 ($ i ' ,BO ij ' )# f 5 ($ j ' ,BO ij ' ), ! BO ij "" =BO ij '"" # f 1 ($ i ' ,$ j ' )# f 1 ($ i ' ,$ j ' )# f 4 ($ i ' ,BO ij ' )# f 5 ($ j ' ,BO ij ' ), ! BO ij =BO ij " +BO ij # +BO ij ## . (2.20) The functions f 1 , f 2 , f 3 , f 4 , f 5 , are defined as follows: ! f 1 (" i ," j ) = 1 2 Val i + f 2 (" i ' ," j ' ) Val i + f 2 (" i ' ," j ' ) + f 3 (" i ' ," j ' ) + Val j + f 2 (" i ' ," j ' ) Val j + f 2 (" i ' ," j ' ) + f 3 (" i ' ," j ' ) # $ % % & ' ( ( , ) exp( ) exp( ) , ( 2 j boc1 i boc1 j i p p f ! " + ! " = ! ! , ! f 3 (" i ," j ) =# 1 p boc2 ln 1 2 exp #p boc2 " i ( ) + exp #p boc2 " j ( ) [ ] $ % & ' ( ) , ) ) ( exp( 1 1 ) , ( 4 boc5 boc i ij ij boc4 boc3 ij i p BO BO p p BO f + ! " " + = ! , ) ) ( exp( 1 1 ) , ( 5 boc5 boc j ij ij boc4 boc3 ij j p BO BO p p BO f + ! " " + = ! . (2.21) boc i ! is the correction for atoms with lone electron pairs, 30 ! " + # = $ ) ( i n j ij boc i boc i BO Val . (2.22) where n(i) are the neighbor atoms of i that are within single covalent-bond cutoff distance. The valency of atom i can be calculated from the corrected bond order term, ! " i =#Val i + BO ij j$n(i) % . (2.23) The the bond energy E bond , which represents the strength of a covalent bond between atomic pairs, include the σ, π and double-π bond terms: ! E bond = "D e # $BO ij # $exp p be1 1" BO ij # ( ) p be2 ( ) % & ' ( ) * "D e + $BO ij + "D e ++ $BO ij ++ % & ' ( ) * j , i , .(2.24) 2.2.2.2 Lone pair energy The lone electron pairs affect the valence angles around the atoms which in turn influence the presence of overcoordinated and undercoordinated atoms, especially oxygens and nitrogens. The valency term, ! " i e , is expressed in terms of the number of outer shell electrons and sum of the full bond order, ! " i e =#Val i e + BO ij j$n(i) % . (2.25) The number of lone electron pairs for atom i, ! n lp,i , is described as, ! n lp,i =int " i e 2 # $ % & ' ( exp )p lp1 2 +" i e )2int " i e 2 # $ % & ' ( # $ % & ' ( 2 * + , , - . / / . (2.26) The penalty energy function for atoms with lone electron pairs exceeding the optimal number is ! E lp = p lp2 " i lp 1+exp#75" i lp ( ) i $ 31 where ! " i lp is the deviation from the optimal number of lone electron pair: ! " i lp =n lp,opt #n lp,i . (2.27) 2.2.2.2 Overcoordination and undercoordination energy The effect of overcoordination (Δ i > 0) and undercoordination (Δ i < 0) are described by a penalty term. The coordination difference Δ i decreases when a lone electron pair appears, ! " i lpcorr =" i # " i lp 1+ p ovun3 exp p ovun4 " j #" j lp ( ) (BO ij $ +BO ij $$ ) j%n(i) & ' ( ) ) * + , , . (2.28) Energies for the overcoordinated and undercoordinated atoms are calculated using the corrected ! " i lpcorr , ! E over = p ovun1 D e " BO ij j#n(i) $ % i lpcorr +Val i % i lpcorr 1 1+exp p ovun2 % i lpcorr ( ) & ' ( ( ) * + + i $ , ! E under = ("p ovun5 ) 1"exp p ovun6 # i lpcorr ( ) 1+exp("p ovun2 # i lpcorr ) $ i % 1 1+ p ovun7 exp p ovun8 # j "# j lp ( ) (BO ij & +BO ij && ) j'n(i) % ( ) * * + , - - . (2.29) 2.2.2.3 Valence angle energies The energy contribution arising from the deviation of valence angle Θ ijk from equilibrium is described in terms of E val . The equilibrium angle Θ 0 is obtained by, 32 ! SBO j = BO jm " +BO jm "" ( ) m#n(j) $ + 1% exp %BO jm 8 ( ) m#n( j) & ' ( ) ) * + , , %- j angle %p val8 n lp,j ( ) , ! " j angle =#Val j angle + BO jm m$n(i) % , ! SBO2 = 0 (SBO"0) SBO p val9 (0 < SBO <1) 2#(2# SBO) p val9 (1< SBO <2) 2 (SBO >2) $ % & & ' & & ! " 0 (SBO) =#$" 0,0 1$exp $p val10 2$SBO2 ( ) [ ] { } . (2.30) The function E val goes to zero as the bond order within an atomic triplet goes to zero: ! E val =" i " j " k f 7 (BO ij )f 7 (BO jk )f 8 (# j )$ p val1 %p val1 exp %p val2 & 0 SBO j ( ) %& ijk ( ) 2 ' ( ) * + , - . / 0 1 2 ! f 7 (BO ij ) =1"exp"p val3 BO ij p val4 ( ) ! f 8 (" j ) = p val5 # p val5 #1 ( ) 2 +exp p val6 " j angle ( ) 1+exp p val6 " j angle ( ) +exp #p val7 " j angle ( ) . (2.31) When two double bonds share an atom in a triplet, such as allene, an additional penalty term is introduced: ! E pen = p pen1 f 9 " j ( ) exp#p pen2 BO ij #2 ( ) 2 [ ] exp#p pen2 BO jk #2 ( ) 2 [ ] k $ j $ i $ , ! f 9 " j ( ) = 2 +exp#p pen3 " j ( ) 1+exp#p pen3 " j ( ) +exp p pen4 " j ( ) . (2.32) Three-body conjugation energy, E coa , sets in for NO 2 -group triplet, which has the following form: 33 ! E coa = p coa1 1 1+exp p coa2 " j val ( ) k # j # i # $exp %p coa3 %BO ij + BO im m&n(i) # ' ( ) ) * + , , 2 - . / / 0 1 2 2 exp %p coa3 %BO jk + BO km m&n(k) # ' ( ) ) * + , , 2 - . / / 0 1 2 2 $exp %p coa4 BO ij %1.5 ( ) 2 [ ] exp %p coa4 BO jk %1.5 ( ) 2 [ ] . (2.33) 2.2.2.4 Torsion angle energies The torsion angle ω ijkl is the angle between the two planes defined by two atom triplets (i, j, k) and (j, k, l). The torsion angle energy E tors decreases as the valence angles approaches π, or as the bond orders goes to zero, that is, ! E tors = 1 2 f 10 (BO ij ,BO jk ,BO kl )sin" ijk sin" jkl l # k # j # i # $ V 1 (1+cos% ijkl ) +V 2 exp p tor1 2&BO jk ' & f 11 (( j ,( k ) ( ) 2 { } 1&cos2% ijkl ( ) +V 3 1+cos3% ijkl ( ) ) * + + + + , - . . . . , ! f 10 (BO ij ,BO jk ,BO kl ) = 1"exp "p tor2 BO ij ( ) [ ] 1"exp "p tor2 BO jk ( ) [ ] 1"exp"p tor2 BO kl ( ) [ ] , ! f 11 (" j ," k ) = 2 +exp #p tor3 " j angle +" k angle ( ) [ ] 1+exp #p tor3 " j angle +" k angle ( ) [ ] +exp p tor4 " j angle +" k angle ( ) [ ] (2.34) The four-body conjugation energy E conj gives: ! E conj = p cot1 f 12 (BO ij ,BO jk ,BO kl )1+ cos 2 " ijkl #1 ( ) sin$ ijk sin$ jkl [ ] l % k % j % i % , ! f 12 (BO ij ,BO jk ,BO kl ) = exp"p cot2 BO ij "1.5 ( ) 2 [ ] exp "p cot2 BO jk "1.5 ( ) 2 [ ] exp "p cot2 BO kl "1.5 ( ) 2 [ ] . 34 2.2.2.5 Hydrogen bond energy The hydrogen bond interaction E hbond depends on bond order and the arrangement of atomic triplets. For the triplet i-j-k, where the atom i is a neighbor of hydrogen j and an atom k within the designated cutoff (~10 Å) from atom i, ! E hbond = p hb1 1"exp p hb2 BO ij ( ) [ ] exp p hb3 r hb o r jk + r jk r hb o "2 # $ % % & ' ( ( ) * + + , - . . sin 8 / ijk 2 # $ % & ' ( k 0 j 0 i 0 . (2.36) 2.2.2.6 van der Waals and Coulomb energies To describe Pauli principle orthogonalization for short interatomic distance and dispersion for long interatomic distance, van der Waals and Coulumb interactions are taken into account among all atom pairs. A distance-corrected Morse-potential is used for van der Waals interactions: ! E vdWaals = D ij Tap(r ij ) exp " ij 1# f 13 (r ij ) r vdW $ % & ' ( ) * + , - . / #2exp 1 2 " ij 1# f 13 (r ij ) r vdW $ % & ' ( ) * + , - . / 0 1 2 3 2 4 5 2 6 2 j 7 i 7 , ! f 13 (r ij ) = r ij p vdW1 + 1 " w # $ % & ' ( p vdW1 ) * + + , - . . 1/p vdW1 , where Tap has the following definition: 0 4 4 5 5 6 6 7 7 35 84 70 20 ) ( ij ij cut ij cut ij cut ij cut ij r r R - r R r R - r R r Tap + + = (2.37) R cut is the non-bonded cutoff radius. Coulomb potential is defined as: ! E Coulomb =C q i q j r ij 3 + 1 " ij # $ % % & ' ( ( 3 ) * + + , - . . 1/3 (2.38) 35 where γ ij is a parameter for the smeared Coulombic function. The atomic charges are updated dynamically using the electronegativity equalization method (EEM). 2.3 Integration algorithms Among the various numerical integration algorithms, integrators that allow a large time step and can still remain accurate and stable over long time are used for the calculation of forces in a MD simulation. In our MD simulations, the velocity-Verlet algorithm is used because it is time- reversible and symplectic. The velocity-Verlet algorithm is a modification of the Verlet algorithm where the position of atom i at time t t ! + is obtained using Taylor series expansion. The expansion around r(t) gives ! r(t +"t) =r(t) +"tv(t) + 1 2 "t 2 a(t) + 1 6 "t 3 b(t) +O("t 4 ) r(t#"t) =r(t)#"tv(t) + 1 2 "t 2 a(t)# 1 6 "t 3 b(t) +O("t 4 ) . (2.39) Adding the two equations, one gets: ) ( ) ( ) ( ) ( ) ( 2 4 t O t t t t t 2 t t ! + ! + ! " " = ! + a r r r , (2.40) and the velocity is obtained from, ! v i t ( ) = r i t +"t ( )#r i t#"t ( ) 2"t +O"t 3 ( ) . (2.41) Therefore, the rounding error in the velocity is higher than the error in the position. In simulations where velocity plays an essential role, the drawbacks of the Verlet algorithm can be eliminated by the velocity-Verlet algorithm. In this algorithm, the velocity at half the time step is given by, t t t t t ! + = ! + ) ( ) ( ) ( 2 1 2 1 a v v . (2.42) 36 The position at time t t ! + is 2 t t t t t t t ! + ! + = ! + ) ( ) ( ) ( ) ( 2 1 a v r r (2.43) From the updated positions, forces and accelerations are computed. The velocity at time t t ! + is given by, t t t t t t t ! ! + + ! + = ! + ) ( ) ( ) ( 2 1 2 1 a v v (2.44) where a(t + ! "t) are the accelerations of atoms after the position update. Recursive use of this algorithm updates positions, velocities and accelerations[83]. 2.4 Physical properties An MD simulation generates phase-space trajectories from which the physical properties of the system can be calculated and direct comparison can be made with experiments as well as more accurate calculations. Based on ergodic hypothesis, physical properties in an MD simulation can be calculated from, ! f =lim "#$ 1 " f(t)dt 0 " % = f(r,p)&(r,p)drdp % . (2.45) ! f represents the time average of a physical quantity f, which, according to the ergodic hypothesis, is equivalent to an ensemble average. r and p are the atomic positions and momenta, respectively. ! "(r,p) is the probability distribution function. 2.4.1 Thermo-mechanical properties Phase-space configurations are used in an MD simulation to calculate kinetic energy, temperature, stress tensor, heat capacity, etc. The average kinetic energy K is related to the temperature of the system, 37 ! K = 3 2 Nk B T =lim "#$ 1 2 m i v i t ( ) [ ] 2 dt i % 0 " & , (2.46) where N is the number of atoms, T is the temperature, k B is the Boltzmann constant, m i and ! v i are the mass and velocity of atom i, respectively. The stress tensor, σ, can be obtained from the Virial theorem, ! " #$ = 1 V m i v i # v i $ i % + 1 2 r ij # f ij $ j % i % (2.47) where α and β are Cartesian coordinate indices {x,y,z}, and ! r ij " and ! f ij " are the vector and force between atoms i and j. One can calculate pressure P from, ! P = 1 3 tr(" #$ ). (2.48) The fluctuations of kinetic energy in the NVE ensemble is related to the specific heat of the system via ! "K 2 = K # K ( ) 2 = 3 2 Nk B 2 T 2 1# 3Nk B 2C V $ % & ' ( ) , (2.49) where V C is the constant-volume heat capacity. 2.4.2 Structural properties Calculations of various correlation functions can provide the spatial distribution of atoms in a system. One such function is the pair distribution function, g, defined as, ! g "# (r 1 ,r 2 ) = V 2 N " N # $ r 1 %r i ( )$ r 2 %r j ( ) j& # { } i' j ( i& " { } ( , (2.50) 38 where V is the volume, ! N and ! N are the number of atoms of species! and ! , repectively, and <…> is the ensemble average. In a uniform system, the pair- distribution function depends on r = r 1 – r 2 : ! g "# (r) = V N " N # $ r%r ij ( ) j& # { } i' j ( i& " { } ( , (2.51) The expression can be simplified for an isotropic system, where the pair-distribution function depends only on r: ! g "# (r) = V 4$r 2 N " N # % r&r ij ( ) j' # { } i( j ) i' " { } ) . (2.52) The total pair-distribution function is a sum over the partial pair distribution-functions. ! g(r) = c " c # g "# (r) ",# $ , (2.53) where c α,β are the concentrations of species, N N c , , / ! " ! " = . Figure 2.1 is an example of pair-distribution of molten silicon at various temperatures. The atomic coordination number are obtained from, N !" (r) = 4#N " V g !" ( $ r ) $ r 2 d $ r 0 r % . (2.54) 39 Figure 2.1: Experimental measurements of pair distribution functions of Si at different temperatures by Kimura et al.[66]. The static structure factor, ! S "# (q), is the Fourier transform of ) (r g !" and can be measured by neutron scattering, ! S "# (q) =$ "# + c " c # N V g "# (r)e iqr dr % . (2.55) If the system is isotopic, the static structure factor becomes, ! S "# (q) =$ "# +4% c " c # N V g "# (r)&1 [ ] sin(qr) qr r 2 dr ' . (2.56) The neutron scattering static structure factor, ! S N (q) , can be obtained from, ! S N (q) = b " b # ",# $ c " c # S "# (q)%& "# + c " c # [ ] b " c " " $ ' ( ) * + , 2 , (2.57) 40 where ! b and ! b are the coherent neutron scattering lengths of species ! and ! , respectively. Bond-angle distribution is a three-body correlation that provides additional information about the spatial distribution of atoms. The nearest neighbor distance r c obtained from the first peak of ! g(r) is used to determine whether a bond between a pair of atoms is formed. If r ij < r c and r jk < r c , the angle of this triplet i-j-k is calculated from ! cos" ijk = v r ij # v r jk v r ij # v r jk and binned. The bond-angle distribution can be compared with NMR measurements[45]. 2.4.3 Dynamic properties The cross correlations between two quantities A and B can be expressed as: ! G AB = A- A B- B A 2 - A 2 ( ) B 2 - B 2 ( ) . (2.58) The correlation is dynamic when both A and B vary with time. If A and B are the same quantities, the function is called an auto-correlation function, ! G auto r,t ( ) = "A r +r o ,t +t o ( ) ( ) "A r o ,t o ( ) ( ) r o ,t o "A r o ,t o ( ) ( ) 2 r o ,t o , (2.59) where ! "A is the fluctuation of A around the mean value. The diffusion constant, D α , of specie α can be computed from the velocity auto- correlation function, ! v auto," t ( ) = v t +t o ( )#v t o ( ) ",t o v t o ( ) ( ) 2 ",t o , 41 ! D " = k B T m " v t ( ) v t o ( ) " v t o ( ) ( ) 2 " 0 # $ dt. (2.60) One can also calculate the diffusion constant from the mean square displacement, ! ! 2 t t t t D ) ( ) ( 6 1 lim 0 r r " = # $ . (2.61) The Fourier transform of the velocity auto-correlation function is proportional to the phonon density of states, ! D(")# e $i"t v auto (t)dt % . (2.62) Dynamic structure factor is another important quantity, which describes the dynamics of the system. The dynamic structure factor, S(q,ω), is obtained from, ! s(k,") = 1 2# e $ikr i (0) e ikr j (t) % e $i"t dt i,j & . (2.63) 2.5 Parallel algorithms In an MD simulation, the force acting on an atom i is determined by the positions of surrounding neighbor atoms within a finite interaction range. This characteristic has made spatial decomposition a commonly used parallel decomposition scheme for MD simulations. The spatial decomposition technique utilizes the divide-and-conquer strategy, where the simulation volume is divided into ! P =P x "P y "P z sub-domains of equal volume, which are mapped to processors in a parallel computer. Atom i at position r i = (r ix , r iy , r iz ) belongs to the spatial domain of processor p by: ! p = p x P y P z + p y P z + p z , ! " # # # # L P r p i / = . (2.64) 42 where α denotes the Cartesian coordinates index x, y, z and L α represent the lengths of the MD box. Information about positions, velocities and atom types of all atoms within the spatial sub-domain is assigned to the mapped processor p. The domain of every processor is extended to access positions of atoms within a cutoff distance from the boundaries of neighboring processors, and the subsequent force calculation is then local to the processor. Figure 2.2 shows the “extended” domain (red) of a processor p. Atom information near the boundary of the neighboring processors need to be cached into p. Figure 2.2: A schematic 2-D view of spatial decomposition used for parallel MD implementation. The internode communication is necessary for force calculation for atoms on other nodes and also for atom migration when the updated position is outside the local boundary. To minimize the communication time, the system maintains a large atom/processor ratio so that caching of atoms from neighboring processors is p 43 negligible compared to the local computation time. This internode communication involves 26 neighboring domains, where data are being sent and received. The “caching” process is done in six steps by sending attributes of the skin-layer atoms to left, right, top, bottom, back and front neighbors consecutively. When atoms move out of the domain of a processor after the position updates, data attributes of those atoms are transferred to the neighboring processor. Message passing interface (MPI) library is used to communicate between processors. Each processor is assigned a unique identification number. For a 3D simulation, the system with P processors is divided into an array of processors ! P x "P y "P z . The sequential index i and the vector indices i x , i y , i z of a processor are defined as follows: ! i x =i P y "P z ( ) ! i y = i P z ( )modP y ! i z =imodP z ! i =i x "P y "P z +i y "P z +i z . (2.65) The arrangement is associated with spatial decomposition, and the assigned information is acquired by the specific processor. Linked-cell list method and neighbor-list method are also used in our MD simulation. The linked-cell list method decomposes the subsystem into smaller cells whose dimensions are slightly larger than the potential cutoff, r c +! , where the “skin” δ is chosen so that atoms do not move more than δ/2 in time δt. In the linked- list data structure, an array of data is stored with each of the record pointing to the 44 next one. In each cell, atoms are sorted and their corresponding information is stored in the linked list. To calculate the force on an atom i, one simply needs to calculate the contributions from atoms in the same cell and from atoms in the 26 nearest- neighbor cells. The Verlet neighbor list method calculates only half of the nearest- neighbor cells to obtain the neighbor-list. Although the neighbor list requires additional memory, the force computation is reduced to <M>N steps, where <M> is the average number of neighbors of an atom. 2.6 Performance of parallel MD algorithm Shared-memory and distributed-memory systems are available for parallel computations. A shared-memory system represents machines where more than one processor accesses the same memory. The size of the private memory assigned to the processor is limited compared to the major memory. Multi-processor computer belongs to this category, where each processor runs the tasks independently of other processes. Processors residing in one node has the advantage of fast data communication speed. In a distributed-memory machine, each processor has its private memory and inter-processor communication is done via the interconnect. The data communication scheme in such a machine has to be implemented by end user using parallel programming libraries such as Message Passing Interface (MPI) and Parallel Virtual Machine (PVM). This architecture has the scalability advantage. A Beowulf cluster is an example of a distributed-memory system and has been adopted as a popular architecture due to its cost efficiency. 45 A hybrid distributed-shared memory system is a combination of the above two architectures. It consists of linked computing nodes where each node is a multi- processor machine. The 6,020-processor Linux cluster at the High Performance Computing Center (HPCC) and the 2,048-processor Linux clusters at the Collaboratory for Advanced Computing and Simulations (CACS) use this architecture. The efficiency, E, is an important measure of parallel algorithm. It is defined as, p) , T(W PW T(W,1) W p S E p 1 p p = = . (2.66) S p is the speedup of the program on p processors: ! S p = S(W,p) S(W,1) = W p T(W,1) W 1 T(W,p) , (2.67) where T(W,1) and T(W,p) are the total parallel program execution times on 1 and p processors, respectively. In an ideal situation, a simulation on p processors should be p times faster than on a single processor with the same work load. Therefore, an ideal algorithm has an efficiency of 1. 46 Figure 2.3: The benchmark result for silica on a Cray T3E cluster. The wall-clock time (open circles) and inter-processor communication time (open squares) per MD step with isogranular scaling are presented. The number of atoms per processor remains constant = 648,000. The result shows communication time is negligible compared to the computation time. Figure 2.3[64] shows the Cray T3E benchmark for the performance of the MD algorithm for silica. The overall efficiency is above 90% and the highest value is 97% on 1024-processor. We choose a large number of atoms per processor (648,000 atoms) to minimize the communication overhead. Figure 2.4 shows the benchmark results for the ReaxFF algorithm used in MD simulations of energetic materials (RDX)[91]. The grain size is chosen to be 10,752 on the BlueGene/L at Lawrence Livermore National Laboratory and 36,288 on Columbia at NASA-Ames. The parallel efficiency is close to 1 on these machines. 47 Figure 2.4: The isogranular scaling benchmark results for RDX on (A) BlueGene/L at Lawrence Livermore National Laboratory and (B) Altix 3000 at NASA-Ames. The grain sizes are 10,752 on the BlueGene/L and 36,288 on the Altix. 48 CHAPTER 3: NANOCAVITATION DAMAGE IN SILICA GLASS When a crack propagates under the influence of an external load, a part of the mechanical energy is dissipated ahead of the crack tip in a region called the process zone. The dissipation occurs through plastic deformation, damage, release of elastic waves, chemical processes, etc.[9, 128] In a brittle material, it is believed that fracture occurs solely by de-cohesion of atomic bonds. However, experimental studies of the morphology of fracture surfaces for a variety of materials reveal scaling features that are independent of the precise nature of the material under consideration or conditions of loading.[27] This suggests that damage may play a generic role within the process zone. A series of MD simulations have been carried out to study the nature of deformation and fracture in silica glass under different mechanical loading. Rountree et al.[107] reported nanocavitation in front of the crack tip in silica glass under Mode I tensile strain. Lu et al.[78] performed multiaxial compressive loading simulations which also revealed that nanocavitation plays a pivotal role in the release of mechanical stress ahead of the crack. These two simulations are summarized in the first part of this chapter. To study atomistic mechanisms of deformation, damage, and fracture, we have performed multimillion-to-billion atom MD simulations of amorphous silica under hydrostatic tension. We find that nanocavities in a-SiO 2 arise from the migration of non-bridging oxygen atoms in –Si-O-Si-O-Si– rings and that these nanocavities play an important role in void-void coalescence and inter-void ligament failure. Mott also 49 proposed this defect transport mechanism for flow in silica glass[89]. The simulation results are presented in the second part of this chapter. 3.1 Tensile and compressive fracture simulation Rountree et al.[107] performed molecular dynamics simulations to study dynamic fracture in a-SiO 2 under uniaxial tensile strain. The simulated system consisted of 110 million atoms in a cube of (120.6nm) 3 . Figure 3.1(a) shows the initial setup of the simulation where atoms at the top and bottom layers are kept frozen. Figure 3.1(b) shows cavities in front of the crack tip at 55ps. These cavities grow, coalesce among themselves and with the advancing crack front to cause material failure. The average crack velocity can be grouped into two regimes: (1) the crack propagation velocity is about 400m/s when the strain is between 2.2% to 3.6%; and (2) the crack velocity is about 1300 m/s at higher strains. The simulation shows that the system fractures when the crack velocity reaches about 50% of the speed of Rayleigh waves (3,000 m/s). 50 Figure 3.1: (a) Initial configuration of the system under uniaxial strain. The frozen layers are shown as dark gray regions. The strain is applied on the top and bottom layers by displacing boundary-layer atoms. (b) Snapshot shows cavities in front of the crack, and their coalescence leads to failure. Recent MD simulations of compressive fracture in a pre-notched amorphous silica system[78] reveals how the pre-crack kinks and forms a wing crack. The multimillion-atom system has the dimensions 120 ! "120 ! "15 nm 3 . Figure 3.2(a) shows the simulation setup and the direction of the applied strain. 51 Figure 3.2: (a) Illustration of the wing crack simulation setup. The gray plate represents a rigid indenter and the blue region is the pre-cracked a-SiO 2 system. The indenter speed is 150 m/s. (b) Different nanocavities are shown in different colors. At 19 ps, nanocavities grow and coalesce into nanocrack columns. At 21 ps, nanocolumns continue to grow and form a wing crack. Figure 3.2(b) shows the process by which the wing crack forms. At 16ps, nanocavities form at the pre-crack tip due to frictional sliding of pre-crack surfaces. The pre-crack kinks at an angle of 70° relative to the z-axis. This is in agreement with a theoretical prediction[41] and dynamic compression experiments[18, 19, 35, 75, 76]. The nanocavities continue to nucleate around the kink and coalesce with it to form nanocolumns. The maximum damage zone of nanocavities in front of the crack tip is about 300 Å. At 21 ps, the wing crack propagates with a velocity of 1000 m/s. 52 Subsequently, the wing crack heals at an average speed of 1300 m/s due to reflection of stress waves from the boundaries of the system. 3.2 Void coalescence under tension Cavitation is a ubiquitous form of damage in ductile fracture of metallic alloys[121]. Voids in metallic alloys nucleate at grain-boundary junctions and from decohesion or cracking of inclusions[23, 74]. Argon et al. proposed that cavities could form if the inclusion diameter is greater than 10nm and the interfacial stress exceeds a critical value[10]. The growth of an isolated void and various mechanisms by which voids join in metallic systems have been studied extensively for many decades[33, 37, 38, 85, 106]. Experimental studies reveal localized, microscopic plastic deformation in intervoid ligaments prior to coalescence and ligament failure by microcracking or shearing. In high strength AISI-4340 steel, ligament failure is observed to be due to the merging of a sheet of cavities that are an order-of- magnitude smaller than the voids around inclusions[42]. Brown and Embury[32] proposed a simple criterion for void coalescence in a perfectly plastic material under tension. According to them, a shear band forms and the ligament fractures when the void length equals the center-to-center separation between neighboring voids. Another criterion for the onset of void coalescence depends critically on porosity, which is regarded as a material constant [85]. Simulation studies of void growth and coalescence in metallic systems have been based mostly on continuum models[103, 122]. Recently, however, molecular dynamics (MD) simulations have also been used to investigate the onset of void-void coalescence in copper under hydrostatic tension [108]. 53 A long-held belief is that brittle materials fracture solely by de-cohesion of atomic bonds at the crack tip [73]. However, MD simulations under Mode I tension[107] and compressive loading[78] indicate that damage cavities release mechanical stresses ahead of the crack tip. Recent Atomic Force Microscopy (AFM) studies have generated considerable controversy about the nature of damage in stress corrosion cracking of brittle amorphous solids. The Saclay group observes the nucleation, growth and coalescence of nanometer scale cavities as the key mechanism for crack extension in silica and aluminosilicate glasses[26, 34], whereas the NIST group finds no experimental evidence for cavitation around a crack in silica glass[54]. We have performed multimillion-to-billion atom MD simulations of void coalescence in silica glass under hydrostatic tension. Our simulations show nucleation of nanometer scale cavities in inter-void ligaments, which is triggered by the coalescence of nearest-neighbor Si-O rings when a non-bridging oxygen atom binds with a silicon atom and concurrently the latter severs its bond with another oxygen atom. At higher strains, nanocracks appear on void surfaces and ligaments fracture due to the growth and coalescence of ligament nanocavities. Many of these features in silica glass are surprisingly similar to those observed in ductile metallic alloys. We have performed a series of MD simulations of void coalescence that contain between 1 million to 1 billion atoms. Two of the simulations reported here contain 500 voids in systems with a billion atoms each in a 319.5×296.7×179.7 nm 3 MD box. In addition, we have investigated growth and interaction between a pair of voids of varying sizes and separations in systems with 1 million and 15 million atoms in MD boxes of dimensions (25.6nm) 3 and (62.8nm) 3 . We impose periodic boundary 54 conditions and apply a dilatational strain at a constant rate of 10 8 or 10 9 sec -1 using the Parrinello-Rahman approach[96]. The equations of motion for the atoms are integrated with the velocity-Verlet algorithm[4]. Silica glasses are generated from an ideal β−crystobalite structure using the melt-quench method[90]. After preparation, the amorphous systems are well thermalized at room temperature and spherical voids are created by removing atoms. These systems with voids are relaxed with the conjugate gradient method and subsequently heated gradually back to room temperature where they are well thermalized again. Each system with a million atoms was run for 550 ps and each billion-atom system for 120 ps. Figure 3.3(a) shows the growth of two voids in a million-atom a-SiO 2 system. In this simulation the initial diameter of the voids is 3nm, the center-to-center void separation is 6 nm, and the applied dilatational strain rate is 10 8 sec -1 . With an increase in the applied strain, the void diameters grow slowly at first from their initial size of 3 nm to 3.5 nm and when the strain, ε, is about 4% nanometer size cavities nucleate in the inter-void ligament region. In order to delineate damage nucleation and evolution due to void-void interaction from that caused by independent void growth, we have also simulated the growth dynamics of one void under the same initial conditions (1 million-atom system with a void of diameter 3 nm) and strain rate (10 8 sec -1 ). Figure 3.3(b) shows the ratio of the average volume per void, φ 2 , in the system with 2 voids to the void volume in the single-void system, φ 1 . Up to a strain of about 7%, the average growth rate of two voids is nearly the same as that of a single void but there are no damage cavities in the latter system. Starting at a strain of 7%, nanocracks appear on the surfaces of the two voids and their average void 55 growth rate far exceeds the growth rate of a single void; see Fig 3.3(b). In the two- void case, inter-void ligament cavities grow and at ε = 8% the inter-void ligament fractures due to the coalescence of voids with the nanocavities, see Fig 3.3(c). As a comparison, figure 3.3(d) shows the single void system at ε = 8%. The void growth in the two-void system is dramatic. At coalescence the inter-void ligament distance is close to the radius of either void, which agrees qualitatively with the Brown-Embury criterion for ductile materials. 56 Figure 3.3: (a) Snapshot of voids and nanocavities at a strain rate of 10 8 sec -1 . (b) Strain dependence of the average porosity per void in the two-void system, φ 2 , relative to the porosity of the single-void system φ 1 . (c) Snapshot of voids and the fractured inter-void ligament at ε = 8%. (d) The single void system at ε = 8% shows little change in the void size. To determine the influence of the void-void interaction on structural changes in the intervoid ligament region prior to ligament fracture, we monitor the Si-O pair correlation function, ! g Si"O (r), and Si-O-Si bond angle distribution, ! P Si"O"Si (#), as a function of the applied strain, ε. Figure 3.4(a) shows the changes in the first peak in 57 ! g Si"O (r) at ε = 8% relative to that in the unstrained system. The height of the first peak drops significantly because of the decrease in the coordination of Si resulting from Si-O bond breaking. The height of the second peak in ! g Si"O (r) also drops, which indicates changes in the –Si-O-Si-O-Si– ring structure. These changes are also evident in the bond angle distribution ! P Si"O"Si (#) shown in Fig. 3.4(b): At ε =4%, ! P Si"O"Si (#) is narrower and shifted to higher angles compared to the unstrained system. At higher values of ε, Si-O bond breaking causes a partial shift in ! P Si"O"Si (#) back towards the unstrained case. Figure 3.4: (a) Si-O pair distribution functions and (b) Si-O-Si bond angle distributions in the middle of the inter-void ligament in the unstrained and strained systems. Detailed analysis of cavity nucleation reveals a novel mechanism involving strain-enhanced defect transport; see Fig. 3.5(a) - (c). In the unstrained a-SiO 2 , each Si atom (red) is connected to 4 O atoms (yellow) and these SiO 4 tetrahedra are linked into nanometer size –Si-O-Si-O-Si– rings through corner-sharing O atoms. In Fig. 58 3.5(a), blue, green and gray regions are inside 7-, 6- and 5-membered Si-O rings, respectively. The blue atom in Fig. 3.5(a) is a non-bridging O bonded to a Si in the 7- membered blue ring. The oxygen atoms, shown in green and white, play a pivotal role in the nucleation of a nanocavity. Figure 3.5(b) shows that at a strain of 1% the non-bridging O (blue) becomes fully coordinated by binding with a Si while the white O atom becomes under-coordinated by dissociating from that Si atom. As a result, the 5-membered gray ring becomes a 12-membered ring and the blue oxygen is now a part of a 3-membered ring. At a strain of 4%, the green O dissociates from a Si atom and the latter binds with the white O atom, thus creating an 11-membered ring (gray) adjacent to a 10-membered ring (green); see Fig. 3.5(c). The transport of non- bridging O is driven by stress gradients. This mechanism introduces a natural length scale, i.e., ring diameter (~ 1 nm), into the nanocavity nucleation problem. Figure 3.5(d) shows a considerable change in ! P Si"O"Si (#) for the rings with bond-switching events involving non-bridging O: instead of a single broad peak at 143º (Fig. 3.4a), the distribution has another peak at 160º due to strain-induced growth of rings. Void- void coalescence studies in a 15 million-atom SiO 2 glass at a hydrostatic strain of 10 9 sec -1 reveal the same mechanisms for the nucleation and growth of damage nanocavities and inter-void ligament failure. 59 Figure 3.5: Bond switching mechanism involving non-bridging oxygens is shown in (a) and (b) by red and white dashed lines between blue and white atoms. In (c), the white dashed line indicates bond switching between white and green oxygens. (d) shows the Si-O-Si bond angle distribution for the rings involved in the bond- switching events. Red and yellow spheres represent Si and O atoms, respectively. Figure 3.6 shows evolution and coalescence among multiple voids in a billion- atom system subjected to a dilatational strain of 10 9 sec -1 . The initial configuration is an array of 500 voids distributed at the center of the specimen, about 100 nm from the boundaries in the x and y directions and 60 nm from the boundaries in the z direction; 60 see Fig. 3.6(a). The initial diameter of each void is 3 nm and the smallest center-to- center void separation is 9 nm. As the strain is increased to 4.5%, sheets of nanocavities nucleate in inter-void ligaments and small cracks appear on void surfaces. Nanocavities form in ligaments between nearest neighbor as well as next nearest neighbor voids (diagonal ligaments). The voids and ligament cavities grow and merge to fracture the ligaments. The ligaments first fail at a strain of 8% around the corner voids. With a slight increase in strain (9%), ligaments between many different voids fracture (Figure. 3.6(b)). This releases the tensile stress on voids close to the fractured ligaments and as a result they shrink in size and cracks on their surfaces begin to heal. Figures 3.6(c) and 3.6(d) show inter-void ligament fracture and crack formation resulting from the coalescence of many voids at a strain of 12%. 61 Figure 3.6: Voids in a billion atom system under a strain rate of 10 9 sec -1 . (a) shows a slice of the entire system; (b) shows cracks on void surfaces at ε=9%; (c) and (d) show void coalescence and inter-void ligament failure at ε=10.5% and ε=12%, respectively. We have also investigated the effect of an applied dilatational strain (at a strain rate of 10 9 sec -1 ) on a random spatial distribution of 500 non-overlapping voids in a billion-atom a-SiO 2 system. The initial size of each void is 3 nm and the void centers are chosen from a uniform random distribution. Figure 3.7 shows the porosity of this system, and for comparison the porosity of the system with regularly distributed voids, as a function of the applied strain. As the strain is increased, the voids first 62 grow independently of one another without any significant damage around them. However, when the strain exceeds 8% small cracks appear on the surfaces of some voids along with the onset of void coalescence; see Fig. 3.7(a). The void-void interaction and eventual coalescence are again mediated by the nucleation and growth of inter-void ligament nanocavities. Figure 3.7(b) shows a mosaic of cracks resulting from the percolation of voids through the system at a strain of 10%. Figure 3.7: Porosity vs. strain in billion-atom systems containing 500 identical voids distributed regularly and randomly. (a) and (b) show volumes of coalesced cavities in the latter system at ε=8% and 10%, respectively. 63 CHAPTER 4: SHEAR DEFORMATION IN SILICA AND RDX This chapter focuses on shear deformations in silica and RDX. Our results reveal that shear-induced void deformation, damage and flow in silica glass have the same underlying mechanism involving Si-O bond breaking and the migration of three-fold coordinated silicon and non-bridging oxygen defects on —Si-O-Si-O— rings. This mechanism is similar to the one observed in void coalescence under hydrostatic tension, which was described in the previous chapter. The second part of this chapter deals with shear deformation simulation of an RDX crystal. Our results show molecular structure of RDX crystal plays a pivotal role in the deformation process. 4.1 Shear deformation in silica Bubble deformation and breakup in sheared viscous fluids is a forefront research area with applications ranging from food processing, ink-jet printing, glass manufacturing, blending of polymer melts to drug delivery systems and the rheology of magmas[3, 79, 102, 113]. Here we describe atomistic mechanisms of deformation and breakup of voids in silica glass subjected to a high shear rate in multimillion atom molecular dynamics simulations. With an increase in the shear strain, we observe two kinds of point defects—three-fold coordinated silicon and non-bridging oxygen atoms—as spherical voids deform elastically into ellipsoidal shapes. For shear strains ε > 15%, the number of Si and O defects increases as nanocracks appear on void surfaces and voids deform plastically into a threadlike structure before breakup. Nanocracks are nucleated by the migration of three-fold coordinated Si and non- bridging O on —Si-O-Si-O— rings. For ε > 40%, the voids break up into several 64 fragments. These atomistic simulation results in silica glass have much in common with macroscopic-level experimental observations and theoretical descriptions of bubble deformation and breakup in shearing fluids. Applications of drop deformation and breakup phenomena in fluids span a wide range of spatio-temporal scales. At the continuum level, drop deformation and breakup have been investigated quite extensively[3, 102, 113] as a function of the ratio of the drop viscosity (λµ) relative to that of the fluid (µ) and the capillary number ! C = Gaµ " , where G is the applied shear rate, and a and γ are the drop radius and interfacial tension, respectively. The capillary number represents a competition between shear- flow forces which tend to deform the drop and the interfacial tension which tends to maintain the spherical shape. For small Reynold’s number of fluid motions ( ! "au/ µ <<1, where ρ and u are the fluid density and velocity, respectively), a spherical drop deforms into an ellipsoid and the deformation parameter D (= (L - B)/(L + B), where L and B are the length and breadth of the drop, respectievly) is proportional to C for C << 1. As C increases, the drop becomes more elongated and turns into a threadlike structure before breaking up into smaller droplets above a critical value of C. Long, slender shapes are observed for an inviscid drop or a bubble[113] and in these instances of low viscosity ratios λ << 1, the drop behavior is well described by the slender-body theory[58]. Hinch and Acrivos predicted on the basis of this theory[59] that a spherical bubble under simple shear would deform into an S-shaped, long, thin bubble in the limit ! C"#. Shape measurements on air bubbles in viscous liquids support their prediction and rheology data on bubble shapes in lava flow at high capillary numbers (rhyolite) also reveal pointed ends on 65 the bubbles[79, 102, 113]. Joseph has argued that cavitation in liquids have much in common with crack nucleation in amorphous solids and that a glass is an ideal system to test this commonality[63]. We present molecular dynamics simulation results for a single void in amorphous silica subjected to a high shear rate between 10 11 and 10 12 sec -1 , where the strain is applied every 0.1 ps. The simulations involve up to 633 million atoms. The initial diameter of the void ranges between 3 nm and 50 nm, which covers nearly the entire range of void sizes observed in the damage zone in dynamic fracture simulations[107] and quasi-static stress corrosion cracking experiments[34] on silica glass. The simulations are based on an interatomic potential, which incorporates ionic and covalent effects through a combination of two-body and three-body terms[126]. The interatomic potential has been validated extensively by comparing the MD simulation results for structural and mechanical properties of a-SiO 2 with experimental measurements and quantum mechanical calculations based on density functional theory. Figure 4.1 shows the degree of deformation for a void in silica glass at room temperature. The initial diameter of the void is 10 nm. As time evolves, the void shape changes to ellipsoidal and damage in the form of nanocracks nucleates on the void surface at the two ends of the short axis of the ellipsoid; see the snapshot in Fig. 4.1(a). The same kind of deformation and damage is also observed for voids of initial diameters 3 nm and 50 nm. Stress calculations[36] reveal that nanocracks form along the direction of the maximum tensile stress. Joseph has proposed a cavitation criterion in a flowing liquid, according to which a cavity nucleates in the direction of 66 maximum tensile stress much like the way fracture occurs in amorphous solids[63]. In the case of a Newtonian fluid in pure shear, the cavitation direction is predicted to be 45° from the shear plane. This is indeed the direction of nanocrcak initiation on void surfaces in our simulations. Figure 4.1(b) shows the time variation of the deformation parameter, D, for initially spherical voids of diameters 3 nm (blue), 10 nm (red) and 50 nm (green). In all three cases D increases linearly with time up to 15 ps, which is an indication that the void deformation is elastic up to a shear strain of 15%. To distinguish between elastic and plastic void deformations, we switch off the shear strain at a certain value of the strain and let the system relax without shear. We find that if the shear strain rate is dropped to zero before the strain reaches 15% voids recover their initial spherical shapes and sizes. This confirms that the void deformation is elastic so long as the strain is below 15%. 67 Figure 4.1: (a) Shear deformation of a spherical void of diameter 10 nm into an ellipsoid after 15 ps; (b) time variation of the deformation parameter, D, for voids of initial diameters 3 nm (blue circles), 10 nm (red squares), and 50 nm (green triangles). The applied shear rate is 10 12 sec -1 . Plastic deformation in shearing silica glass appears after 15 ps (strain > 15%) as the deformation parameter begins to increase non-linearly with time. Figures 4.2(a) and 4.2(b) are snapshots showing long, thin voids with nanocracks on the surfaces near the middle of the long axis. Initially, these were spherical voids of diameters 3 nm and 50 nm. In both cases, the voids become threadlike before breakup, which occurs above a strain of 40%. Figure 4.2(c) shows plastic deformation of a void, which was initially a sphere of diameter 10 nm. The snapshot taken at 35 ps shows a long, thin void with ends deformed roughly like an S shape. There are nanocracks on the surface and near the ends of the void. As time evolves, nanocracks grow and the void becomes more elongated, and at 40 ps the void fragments; see Fig. 4.2(d). 68 Figure 4.2: Panels (a) and (b) show plastic deformations of voids of initial diameters 3 nm and 50 nm, respectively. Panel (c) shows a plastically deformed void which was initially a sphere of diameter 10 nm; and panel (d) shows its breakup at 40 ps. To confirm that the void deformation is indeed plastic in the non-linear deformation regime, we have performed several simulations in which we switch off the shear strain after 15 ps and let the system evolve without any external strain. Figure 4.3(a) shows the snapshot of the void at 15 ps, and the yellow arrow indicates the location of the nanocrack. When the shear strain is switched off, the nanocrack on the void surface heals, see Fig. 4.3(b). However, we find that neither the nanocracks nor the pointed ends heal completely if the strain exceeds 15% (Fig. 4.3(c) and (d)). This indicates that the system has undergone plastic deformation. When the shear is turned off in the case of the S-shaped void, the pointed end fragments in a manner akin to the “end-pinching” mechanism of drop breakup in fluids[52, 113, 118]. Note that Hinch and Acrivos[59] predicted an S-shaped drop in the limit λ << 1 and C >> 69 1, and experiments on bubbles in viscous liquids in simple shear flow support their theoretical analysis based on the slender-body theory[113]. Threadlike deformations of inviscid drops have also been observed experimentally in fluids under high shear rates[113, 118]. Figure 4.3: Snapshots of voids in elastic and plastic regimes. Panels (a) and (c) show deformation of voids under strain at 15 ps and 30 ps, respectively. We switch off the strain and relax the systems for 5 ps. Panels (b) and (d) show the deformation of those voids after 5ps. In panel (b) the nanocracks have healed, and in panel (d) they have not. Detailed analyses of elastic-to-plastic void deformation and crack initiation and growth reveal a novel mechanism involving strain-enhanced defect transport; see Fig. 4.4(a) - (c). In the unstrained a-SiO 2 , each Si atom (yellow) is connected to four O 70 atoms (red) in the form of a SiO 4 tetrahedron and these tetrahedra are linked into nanometer size –Si-O-Si-O– rings through corner-sharing O atoms. In Fig. 4.4(a), green and blue regions represent 6- and 9-membered rings, respectively, at a strain of 5%. The magenta atom is a bridging O and the blue atom is a 3-fold coordinated Si atom in the 9-membered ring. Figure 4.4(b) shows that at a strain of 8% the magenta O atom becomes non-bridging while the blue Si atom remains under-coordinated. The grey ring with these Si and O defects is a 13-membered ring. At a strain of 10%, the green Si and magenta O atoms become fully coordinated by bonding with each other. The blue Si is still under-coordinated, but now belongs to a 14-membered ring (yellow). Non-bridging oxygen and three-fold coordinated silicon atoms play pivotal roles in the nucleation of cracks through the enlargement of –Si-O-Si-O– rings, and in void deformation and shear flow in silica glass. 71 Figure 4.4: Point defects responsible for shear deformation and flow, and crack nucleation in silica glass. (a) A 6-membered –Si-O-Si-O– ring (green) adjacent to a 9-membered ring (blue) with a fully coordinated O (magenta) and a 3-fold coordinated Si defect (blue). (b) At a strain of 8%, the magenta O atom becomes non-bridging and the blue Si atom remains under-coordinated. These point defects reside on a 13-membered ring (grey). (c) At a strain of 10%, the green Si and magenta O atoms become fully coordinated by bonding with each other, but the blue Si is still under-coordinated. The yellow region is inside a 14-membered ring. 72 Figure 4.5(a) shows the time variation of the density of under-coordinated Si and O defects in a region that encloses the deformed void (initial diameter 10 nm) and the damage zone around it. We find that this region has approximately the same number of three-fold coordinated Si and non-bridging O atoms. The defect density remains small for about 15 ps and then increases over the next 25 ps before leveling off as steady plastic flow is established in the system. Voids break up as soon as the defect density saturates. Figure 4.5(b) shows the Si-O pair-correlation function, g Si-O (r), for defects when the void is elastically deformed (15 ps) and after it breaks up (55 ps). In both cases, there is only one prominent peak at 1.57 Å, which is shifted to the left relative to the first peak in the g Si-O (r) in bulk a-SiO 2 (1.62 Å) because the defects are under-coordinated. The peak height drops significantly after the system undergoes elastic-to-plastic transformation. We have also investigated the dynamics of these defects and find steady plastic flow along the shear direction. Perpendicular to the direction of shear flow, the mean square displacement for these defects oscillates in the elastic regime, but increases linearly with time after the plastic flow is established, see Fig. 4.6(a). The diffusion constant for defects in the plastic regime is ~ 10 -5 cm 2 /s, which is typical of liquids. Figure 4.6(b) shows the velocity auto-correlation function for defects at a strain of 15% which we have monitored for 5 ps. These results on defect characteristics are in accordance with Mott’s theoretical model[89] in which under-coordinated Si and non-bridging O point defects diffusing on —Si-O-Si-O— rings are believed to be responsible for plastic flow in silica glass. 73 Figure 4.5: (a) Density of Si and O defects in the damage zone as a function of time; and (b) Si-O pair-distribution function for these defects. The initial void diameter was 3 nm. Figure 4.6: (a) Mean square displacements for defects at strains of 15% and 55%. (b) Velocity auto-correlation functions (VAF) for defects at a strain of 15%. The VAF in the shear direction (Y axis) is close to the VAF in the plane perpendicular to the shear direction. (a) (b) 55% 15% 74 We have also investigated changes in the void shape and elastic-to-plastic transition in silica glass at 1200 K, see Fig. 4.7(a) and (b). In contrast to the case at room temperature, the voids are partially filled with atoms released from the void surfaces and the density of these atoms inside the voids increases with time. The voids still undergo spherical-to-ellipsoidal shape transformation, but the time variation of the deformation parameter is slightly less than that at room temperature. Nanocracks appear on void surfaces with the onset of plasticity, and voids deform into long ellipsoidal shapes with pointed ends and their breakup is again preceded by threadlike structures just as at room temperature. Figure 4.7: Panels (a) and (b) are snapshots of voids under shear strains of 20% and 35%, respectively, at 1200 K. In panel (c), the void is under a strain of 35%. Here the strain rate is 2!10 12 sec -1 . In panel (d), the strain is 35% and the strain rate is 5!10 11 sec -1 . At a slower strain rate, the void breaks up into smaller voids. (a) (b) (c) (d) 75 Simulations have also been carried out to study the effect of strain rate on shear- induced deformation of a void. Figure 4.7(c) shows the deformation of a void at a strain of 35% at a higher strain rate (2!10 12 sec -1 ). The void has a threadlike shape and nanocracks are found on the void surface. The void breaks around a strain of 40%. At a lower strain rate (5!10 11 sec -1 ), the void still undergoes shape transformation and breaks up into smaller cavities at 35%, see Fig. 4.7(d). Nanocracks on void surfaces in these two systems exhibit similar results as in the system at a strain rate of 10 12 sec -1 . Taken all together, MD simulations reveal that shear-induced void deformation, damage and flow in silica glass have the same underlying mechanism involving Si-O bond breaking and the migration of three-fold coordinated silicon and non-bridging oxygen defects on –Si-O-Si-O– rings. The density of these defects remains small and nearly constant in the elastic deformation regime. With the onset of plasticity at a shear strain of 15%, the defect density increases and reaches a plateau soon after the void breaks up. Despite the enormous differences in spatio-temporal scales, the observed shape changes and fragmentation of voids in silica glass are remarkably similar to the deformation and breakup of macroscopic inviscid drops at high shear rates. 4.2 Shear deformation in RDX Million-atom reactive force field (ReaxFF) MD simulations were carried out to investigate deformation of a void in an RDX crystal under a constant shear strain rate of 10 12 sec -1 . The simulation consists of 1.3 million atoms. Figure 4.8(a) shows the simulation setup. The dimensions of the system are 237×231×235 Å 3 . The initial void 76 diameters are 3 nm and 6 nm. The void was created in an RDX crystal by removing molecules around the center of the system. After inserting the void, the system was relaxed at 1 K for 2 ps. The x, y, and z axes in Fig. 4.8(a) are the ! [100], ! [010], and ! [001] crystallographic orientations of the RDX crystal. The shear strain is applied along [010]. The time step in the simulation is 0.1fs. Figure 4.8: (a) Setup of shear deformation simulations in RDX. (b) Time variation of the deformation parameter, D, for a void of initial diameter 3 nm. The applied shear rate is 10 12 sec -1 . Again, as in the case of a-SiO 2 , we calculate deformation parameter D (= (L - B)/(L + B), where L and B are the length and breadth of the void, respectievly). Figure 4.8(b) shows the time variation of the deformation parameter, D, for an initially spherical void of diameter 3 nm. D increases linearly with time up to a strain of 7% and then saturates at a strain of 8%. Subsequently, the void breaks up gradually when it is filled with RDX molecules. In a-SiO 2 , we observe non-linearity in D after the system enters the plastic regime. y x z (a) (b) 77 Figure 4.9: Snapshots of the center layer of RDX at different simulation times. The initial void diameter is 6 nm. (a) Initially, the molecules move with uniform displacement. (b) At a strain of 8%, the molecules from the top and bottom start to move more dramatically than those in the center. (c) The maximum displacement of molecules is found along 45° to the shear direction (strain = 10%). (d) At a strain of 12%, molecules around the void move inward. The color represents the molecular center-of-mass displacement. To analyze the dynamics of RDX molecules, we investigate the center-of-mass displacement of molecules from, ! r r "com = m i r r i (t)# m i r r i (0) i=1 21 $ i=1 21 $ m i i=1 21 $ , where m i and ! r r i (t) y z 78 are the mass and position vector of atom i in the molecule at time t, respectively, and ! r r i (0) is the initial position of atom i. Figure 4.9 shows the molecular displacement at various stages of the simulation. The initial void diameter is 6 nm. We find the displacement to be even in the early stages of the simulation, see Fig. 4.9(a). In Fig. 4.9(b), the yellow dashed box shows that the molecules in the center slab but away from the void move less dramatically than the rest of the molecules at a strain of 8%. The maximum displacement of molecules is along 45° relative to the shear direction (strain = 10%). Starting around 8% strain, the molecules around the void surface move into the void, and the void completely collapses at a strain of 12%. 79 Figure 4.10: Snapshots of temperature distribution around the void. (a) At a strain of 1%, the temperature of a few molecules around the void is about 400 K. (b), (c) and (d) show temperature distributions at strains of 7%, 13% and 19%, respectively. The void is almost filled with high-speed molecules at a strain of 13%. The temperature profile is monitored during the simulation and shown in figure 4.10. The shearing process heat up the molecules around the void. Without breaking a bond, the excess strain energy leads to translational and rotational motion of these molecules. Molecules with high kinetic energy collapse into the void. At a strain of 13%, we find the void is almost filled by these thermally excited RDX molecules whereas the rest of the system remains cool. The maximum temperature is as high as 400 K. 80 We have also investigated the dipole moments of the molecules. For an RDX molecule, the dipole moment, d r , is calculated from ! r d = q i ( r r i " r r com ) i=1 21 # , where q i is the charge of atom i, r i is the atomic position and ! r r com is the center of mass of the molecule. The negatively charged oxygen atoms and positively charged hydrogen atoms contribute to the orientation of the RDX dipole moment which is perpendicular to the C 3 N 3 -ring and pointing inward from the NO 2 -group. Figure 4.11 shows the density of molecules, whose dipole moment difference, (0) d - (t) d r r is greater than 0.5 au. Initially, the density increases linearly with the strain, but exhibits non-linear behavior above a strain of 8%. Figure 4.11: The time variation of density of RDX molecules whose dipole moment differences are greater than 0.5au. 0 1 2 0 5 10 15 20 Density (nm -3 ) Strain (%) 81 CHAPTER 5: NANOINDENTATION IN SILICA AND RDX This chapter focuses on nanoindentation simulations in silica glass and RDX. In the next section, MD simulations of nanoindentation in amorphous silica are presented. These simulations reveal that undercoordinated silicon and non-bridging oxygen defects in the amorphous silica network cause plastic deformation under the indenter and pileup around the indenter. These defects are akin to those postulated by Mott to theoretically explain shear flow in silica glass. As discussed in the previous chapter, these defects were also observed in MD simulations of (a) shear-induced void deformation, damage and flow in silica glass and (b) coalescence of nanovoids and fracture in amorphous silica under tension. The second part of this chapter presents nanoindentation simulation on an RDX crystal. In this simulation, the RDX molecules diffuse/decompose underneath and on all four sides of the diamond indenter. Hardness values vary with the indentation depth and range between 50 and 391 Nmm -2 . Mechanical loading induces localized heating of the RDX substrate and the excess stress is released by the diffusion of RDX molecules on the sides of the indenter. 5.1 Nanoindentation of silica glass In recent years, considerable progress has been made in estimating hardness from indentation testing by combining it with atomic force microscopy (AFM)[87, 94, 101, 128]. In indentation experiments, the hardness of a material is obtained from the ratio of the maximum applied load (P max ) to the contact area (A c ). Since the contact area in nanoindentation experiments is too small for direct observation by optical 82 microscopy, one estimates the projected area (A p ) from a functional relationship between A p and the measured contact depth (h c ). This can introduce considerable uncertainty in the value of hardness. AFM circumvents this problem by allowing direct observation of the impression created by the indenter and precise determination of the contact area[87]. Another advantage of AFM indentation tests over conventional nanoindentation is that indentation impressions can be made shallower by applying smaller normal loads in the range of nano to pico newtons. Miyake et al.[87] have used AFM indentation tests to measure the hardness of fused silica and their estimated value (10 GPa) is considerably lower than the one obtained by using the Oliver-Pharr relationship between the projected area and the contact depth. We have performed molecular dynamics (MD) simulations of nanoindentation in which the applied load and the size of the indentation impression are comparable to those in an AFM nanoindentation experiment[87]. The MD simulation, however, has an advantage over an AFM indentation experiment in that it can provide atomistic- level stress distribution as well as the structure and dynamics of defects associated with deformation and plasticity in the system. Figure 5.1(a) shows a snapshot of the system at the maximum indentation depth. The system was prepared by melting β-cristobalite and then quenching the melt incrementally to a very low temperature[125]. Subsequently, periodic boundary condition in the z direction is removed and atoms within 5.5 Å from the bottom x-y plane are frozen. The amorphous film thus created contains 40 million atoms and has dimensions 100×100×60 nm 3 . The film is heated gradually to 300K and then indented normal to the x-y plane. The indenter, a square-based rigid pyramid with an 83 apical angle of 90°, interacts with the substrate through a repulsive potential. We have used two kinds of indenters to study the effect of the indenter curvature: an atomistically sharp indenter and an indenter whose end is spherical with a radius of 10 nm. Henceforth, the corresponding results are referred to as simulation 1 (atomistically sharp indenter) and 2 (indenter with curvature). Figure 5.1: (a) Pileup around the indenter at the maximum indentation depth. Indentation causes localized damage creating a high-temperature area underneath the indenter, which facilitates pile-up through plastic flow. The colors represent temperature. (b) Load versus displacement during loading and unloading phases of indentation in simulation 1. Figure 5.1(b) displays the applied load in the MD simulation as a function of the penetration depth of the indenter. The load is calculated from the z component of the force on the indenter atoms. From the load and residual indentation impression, we obtain the hardness, H, of a-SiO 2 . In the case of an atomistically sharp indenter, H = 84 7.2 GPa, and for the slightly blunt indenter H = 8.0 GPa. The AFM indentation experiment on fused silica yields a hardness value of 10 GPa [87]. Figure 5.2 shows a snapshot of the indented silica film at an indentation depth of h = 32nm (panel (a)). For clarity, only a slice of the system without the indenter is shown here. At low indentation depths (e.g., h = 11nm), the material close to the bottom of the indenter densifies but there is no material pile-up around the indenter. At the maximum indentation depth of 32nm, we observe significant pileup and a 20% increase in the mass density of a-SiO 2 [25, 65] beneath the indenter relative to the mass density of bulk silica glass under ambient conditions (2.2 g/cc). In the pileup region, the density is 18% lower (1.8 g/cc) than in the bulk system. During the unloading phase of the simulation, the silica substrate exhibits slight elastic recovery and a decrease in the density under the indenter. Panel (b) is a snapshot of a slice of the system after the indenter is pulled out of the substrate. Here the density just beneath the indenter goes back to the normal bulk a-SiO 2 density (2.2 g/cc), but there is still a considerable pileup around the indenter. 85 Figure 5.2: Density distribution at the maximum indentation-depth (a), and after unloading (b). For clarity, the indenter is not shown. Silica substrate deforms until the local density under the indenter reaches 2.6 g/cc, which prevents downward material flow. Significant pileup with relatively low density develops at the maximum indentation depth. The densified region relaxes after unloading. Perriot et al. have used Raman microspectroscopy to characterize plastic behavior in an indentation experiment on fused silica[98]. They have combined Raman spectrum and density to obtain a spatial distribution of densification underneath the indenter. Fused silica is known to densify up to 20% under a b 86 mechanical loading. Perriot et al. have also observed pressure-induced hardening of fused silica in diamond-anvil cell experiments[98]. Detailed analyses of atomistic configurations reveal how plastic deformation and pileup occur in the nanoindentation of a-SiO 2 . In the normal SiO 2 glass, each Si atom is bonded to 4 O in the form of a tetrahedron and each O bonds with two Si atoms to connect neighboring SiO 4 tetrahedra[88]. These corner-sharing tetrahedra form closed –O–Si–O–Si-O– rings of various sizes, with a peak in the distribution around 6-membered rings i.e., 6 Si–O connected pairs. According to Mott, under- coordinated Si and non-bridging O point defects diffusing on –Si-O-Si-O– rings are responsible for plastic flow in silica glass[89]. These point defects in the silica network glass play a central role in the plastic deformation and pileup we observe in both nanoindentation simulations. The defects are involved in bond-switching events, where either a silicon or an oxygen atom severs its bond with one of its neighbors and shortly thereafter (~ 10 ps in simulation 1) forms a new bond with a different neighbor atom. Figures 5.3 (a)-(c) show how a bond-switching event occurs in the case of a non-bridging oxygen atom[1]. The three atoms involved in this event are labeled as O1 (initial oxygen defect), Si2 (bond switching silicon) and O3 (the final oxygen defect). Initially, O1 is a non-bridging oxygen whereas Si2 and O3 have normal coordinations (4 O neighbors for Si2 and 2 Si neighbors for O3). At the transition state, Si2-O3 bond breaks and these two atoms become undercoordinated. Subsequently, Si2 forms a new bond with O1 and they both become fully coordinated whereas the oxygen labeled O3 remains a defect. 87 Another bond-switching event we observe involves the annihilation of a non- bonded pair of undercoordinated silicon and oxygen atoms. Figures 5.3 (d)-(f) show this pair-annihilation event. The four atoms involved in this event are labeled as Si1 (an undercoordinated silicon), O2 (bond-switching oxygen), Si3 (bond-switching silicon) and O4 (a non-bridging oxygen). Initially, O2 is bonded with Si3 and they both have normal coordinations of 2 and 4, respectively. In the transition state, O2- Si3 bond breaks and they both become undercoordinated. In the final state, Si1 forms a bond with O2 and Si3 with O4 and now these four atoms have normal coordinations. 88 Figure 5.3 (a)-(c): Atomic configurations showing transport of oxygen defects; (a) the initial state, (b) the transition state, and (c) the final state. Atoms involved in the event are labeled as O1 (initial defect oxygen atom), Si2 (bond switching silicon) and O3 (final defect oxygen). The oxygen defect migrates by switching bonds from Si2-O3 to O1-Si2. (d)-(f): Atomic configurations showing an annihilation event for a pair of defect atoms; (d) the initial state, (e) the transition state, and (f) the final state. Atoms involved in the event are labeled as Si1 (initial silicon defect), O2 (bond switching oxygen), Si3 (bond switching silicon) and O4 (initial oxygen defect). Si1 and O4 defects annihilate by switching bonds from O2-Si3 to Si1-O2 and Si3-O4. Our simulations also reveal a defect transport mechanism in which a defect migrates a considerable distance via a chain of bond-switching events[44]. This is an important mechanism in the plastic flow of silica glass. Figures 5.4 (a) and (b) show two examples of a bond-switching chain event. In panel (a) an undercoordinated 89 silicon migrates via a bond-switching chain involving seven atoms and in panel (b) the annihilation of a defect pair is mediated by a bond-switching chain of six atoms. Here yellow and red spheres represent positions of silicon and oxygen atoms, respectively, after each multi-atom bond-switching event. In the transport event for a silicon defect, the defect effectively migrates about 1nm in 10 ps. The length of the chain of bond-switching events varies with mechanical loading. Figures 5.4 (c) and (d) show the dependence of the number of bond-switching events N chain on the chain length l at several indentation depths in simulation 1. In the loading phase, longer bond-switching chains appear as the indentation depth increases. Note that the slopes at the three indentation depths are almost constant in the semi-log plot, which implies exponential distributions N chain ∝ exp(-l/l 0 ) with nearly the same characteristic length l 0 . The slopes differ significantly before (h = 31nm) and after (h = 16 and 21nm) the indenter is pulled out of the substrate, which shortens the chain length. 90 Figure 5.4: Snapshots of defect migration pathways ((a) and (b)) in a ball-and-stick representation. Yellow and red spheres represent silicon and oxygen atoms. Black dotted lines indicate covalent bonds before defect migration. Blue and red arrows represent initial and final positions of defect atoms, respectively. (a) Migration of a silicon defect. (b) Annihilation of a pair of silicon and oxygen defects. Panels (c) and (d) show the length of the string-like defect motion as a function of the number of migration events during loading and unloading, respectively. Mechanical loading by the indenter induces longer migration pathways. 91 5.2 Nanoindentation of RDX Mechanical loading of energetic materials is studied to understand the initiation of fast chemical decomposition[12] and the importance of the structure characteristics in their explosive properties. Localized features in deformation and cracking behavior of RDX have been studied experimentally[12]. Vickers micro-indentation experiments have been carried out on RDX single crystals to determine the characteristics[55] of cleavage planes and fracture surface energy, which is regarded as one of the key factors in the performance of explosives. Similar micro-indentation experiments reveal dislocation pile-up activity and localized heating caused by mechanical forces[12, 13]. Infrared measurements reveal significant temperature enhancement in the early stages of plastic flow[86]. To further our microscopic understanding of deformation and thermal properties of RDX, it is necessary to employ atomistic simulations to study chemical reactions under mechanical loading. We have performed such simulations on a (100) α-RDX crystal surface at 150K with a reactive force field (ReaxFF[114, 124]). A unit cell of RDX contains eight covalently bonded molecules [CH 2 N(NO 2 )] 3 and has a {021} slip plane[48, 56] within the orthorhombic crystal structure. The dimensions of the RDX substrate are 127.26×128.52×99.08 Å 3 and it contains 133056 atoms. Periodic boundary conditions are applied in the plane of the substrate. The hydrogen- terminated indenter is constructed by removing atoms from (111) diamond. Dangling carbon atoms on the facets are passivated with hydrogen atoms to eliminate surface radicals. The pyramidal shape diamond indenter has dimensions 92 108.37×108.37×88.66 Å 3 and contains 43342 atoms. In this simulation, the (100) surface of RDX is indented since it has a lower energy than (010) or (001) surface. Both the indenter and the substrate are first relaxed for 11ps in order to bring them to minimum energy states. Atoms within 20 Å of the bottom layer of the RDX substrate are frozen and the rest of the substrate and the indenter are initially at 150K. The simulation time step is 0.15 fs and the indentation speed is 200 m/s, which is 6% of the sound speed in RDX (2.78 km/s [49]). For every 5Å indentation depth, the whole system is relaxed for 1.5 ps and the total indentation depth is 30 Å. The simulation reveals significant increase in temperature under the indenter, which causes localized melting and decomposition of RDX molecular fragments. The temperature distributions are shown in figure 5.5. At the initial stage of the simulation, the area under the indenter deforms without any significant increase of thermal energy. When the indentation depth increases beyond 10 Å, the temperature of the contact region rises dramatically due to plastic deformation and the energy released from deformation is transformed into heat. The area around and underneath the indenter gradually heats up (figure 5.5 (c) to (f)) and eventually melts. Temperature in the damage zone is more than 200K higher than in the rest of the sample. RDX molecules that have the excess energy are found to rotate and their dipole moments point away from the indenter. The hydrogen atoms on the surface of the indenter attract oxygen atoms in the RDX molecules, causing the oxygen atoms to orient toward the indenter surface. There is a limited amount of pile-up and some parts of the fragmented RDX molecules diffuse. The area outside the deformed region retains its original crystal structure. 93 Figure 5.5: Temperature distribution of the RDX substrate at different indentation depths. (a) At 5.2 Å, no heating is observed. (b) At 10.45 Å, molecules around the indenter begin to heat up. (c) At 15.7 Å, the “temperature” of some of the molecules exceeds 500K. At indentation depths of (d) 20.95 Å, (e) 26.20 Å and (f) 31.6 Å, we observe high “temperature” of RDX molecules within ~ 10 Å of the indenter contact area. Hardness is calculated from the ratio of the load to the projected area of the residual impression after the load removal. The size-dependent hardness value ranges from 50 to 391 Nmm -2 . Previous experiments[13, 55] reported Vickers hardness number between 236 and 290 Nmm -2 and Tukon microhardness value between 245 and 363 Nmm -2 . The significance of the size effect is shown in figure 5.6. As the indentation size increases, the hardness value drops dramatically due to the growth of the plastic deformation region. 94 100 200 300 400 5 10 15 20 25 Hardness (Nmm -2 ) Figure 5.6: Indentation size-dependence of hardness value. To study the dynamics of RDX molecules under mechanical loading, we calculated the center of mass displacement, com r ! , which is the difference in the center of mass calculated at time t and 0: ! ! ! = = = " # = 21 1 21 1 21 1 ) 0 ( ) ( i i i i i i i i com m r m t r m r . Figure 5.7 shows com r ! at various simulation times. We find little change in com r ! at the initial stages of the simulation (see figure 5.7(a)). However, com r ! increases dramatically with further indentation. Figure 5.7(b) shows the center-of-mass displacement at the maximum indentation depth of 25.15 Å. Large values of com r ! are found close to the indenter surface and the fragmented RDX molecules move along the indenter surface. The changes in the center-of-mass are concentrated in the ! 1 20 [ ] direction, see figure 5.8. At the maximum indentation depth, the slip becomes obvious. Figure 5.8(b) is a 95 schematic view of the displaced molecules. For clarity, only their center of mass is shown. Molecules above the ! 210 ( ) plane have displaced more than 5 Å from their original positions. However, despite the movements of other molecules, those below the ! 210 ( ) plane remain intact. Figure 5.7: The center-of-mass displacement, com r ! , of RDX molecules at various indentation depths. For clarity, only half of the system is shown. Panels (a) and (b) are snapshots at indentation depths of 10.45 Å and 25.15 Å, respectively. 96 Figure 5.8: Panel (a) shows that the indentation damage is concentrated within the yellow dashed region. Panel (b) is an illustration of how molecules in the yellow dashed region are displaced in panel (a). Each sphere represents an RDX molecule, where the molecules move in the direction of ! 1 20 [ ] . The simulation shows that in order to release the mechanical stress, the RDX crystal exhibits a nontrivial mechanism in addition to pile-up, dislocation emission, and slip. To release the excess stress, RDX molecules diffuse from the substrate. Those molecules are later found “walking” on the indenter surface, see figure 5.9(a) and (b). The dynamics of these fragmented molecules can be decomposed into translational and rotational motions. Figure 5.9(c) is the mean square displacement of such fragmented molecules. The NO 2 groups are attracted to the hydrogen-terminated 97 indenter while the ring is parallel to the indenter surface. This diffusion behavior is similar to 9,10-dithioanthracene (DTA) and 9-thioanthracene (TA) “walking” molecules on a copper surface where the sulfur atoms are attached to the copper substrate[70, 71]. TA molecules tend to rotate around their single sulfur substrate linker whereas DTA have a unidirectional motion in which their two sulfur atoms move alternatively. An RDX molecule has three NO 2 substrate linkers; the alternating motion of oxygen atoms produces the movement of the NO 2 legs and the combined activities of the three legs result in either directional motion towards the top of the indenter or the rotational motion around the molecule itself. The translational velocities of the diffused RDX molecules reach as high as 373 m/s. 98 Figure 5.9: (a) and (b) show fragmented RDX molecules on the indenter surface at time 14.5 and 19.65 ps. Red, dark blue, light blue and white spheres represent oxygen, nitrogen, carbon and hydrogen atoms, respectively. (c): The mean square displacement of a few fragmented RDX molecules as a function of time. The dynamics of these molecules can be decomposed into translational and rotational motions. Molecule A Molecule A Molecule A Molecule B Molecule B (a) (b) (c) 99 CHAPTER 6: CONCLUSION In this dissertation, I have presented results on atomistic mechanisms of deformation, damage, and fracture in silica glass and energetic materials. The results were obtained with multimillion-to-billion atom MD simulations based on an effective force field for silica and reactive force field for in energetic materials. The MD simulations employ highly efficient algorithms for massively parallel architectures. A series of MD simulations under different mechanical loading conditions have been carried out to study the nature of deformation, damage, and fracture in silica glass. Simulations reveal nucleation, growth, and coalescence of nanocavities. The system fractures as a result of coalescence of damage cavities with the crack front. Damage nanocavities also play an important role in the formation of wing cracks under multiaxial compressive loading[78, 107]. We have performed multimillion-to-billion atom MD simulations of amorphous silica under hydrostatic tension. We have observed the onset of void-void interaction. The nanocavities in a-SiO 2 which arise from the migration of non-bridging oxygen atoms in –Si-O-Si-O-Si– rings play an important role in void-void coalescence and inter-void ligament failure. At the onset of coalescence, the radii of voids are close to half the void center-to-center separation. Cavitation in the glass is linked to Si-O rings, which are smaller than cavities around inclusion in metallic systems. To study atomistic mechanisms of deformation and breakup of voids in silica glass under shear strain, we have performed multimillion molecular dynamics (MD) simulations of a single void in amorphous silica subjected to a high shear rate. We 100 have observed two kinds of point defects during the simulations: three-fold coordinated silicon and non-bridging oxygen atoms. The number of undercoordinated defects increases as nanocracks appear on void surfaces and the voids deform plastically into a threadlike structure with the increase in shear strain. We have also found that nanocracks are again nucleated by the migration of three-fold coordinated Si and non-bridging O on —Si-O-Si-O— rings. For ε > 40%, the voids break up into several fragments. We have performed million-atom reactive force field (ReaxFF) MD simulations to investigate deformation of a void in an RDX crystal under a shear strain rate of 10 12 sec -1 . Our results show the molecular structure of RDX crystal plays a pivotal role in the deformation process. The deformation parameter increases linearly first and saturates at a strain of 8%. The MD simulations reveal that the molecules in the center slab but away from the void move less dramatically than the rest of the molecules at a strain of 8%. The maximum displacement of molecules is found 45 o relative to the shear direction. Starting around 8%, the molecules with high kinetic energy around the void surface move into the void. We have found that without breaking a bond, the excess strain energy leads to translational and rotational motion of these molecules. At a strain of 13%, we find the void is almost filled by these thermally excited RDX molecules whereas the rest of the system remains cool. We have performed molecular dynamics (MD) simulations of nanoindentation in which the applied load and the size of the indentation impression are comparable to those in an AFM nanoindentation experiment[87]. The simulations reveal migration of defects and their recombination in the densified plastic region under and the 101 material pileup region around the indenter. The defect transport mechanism where a defect migrates a considerable distance via a chain of bond-switching events[44] is an important mechanism in the plastic flow of silica glass. The simulation result for hardness is close to the value estimated from atomic force microscopy indentation experiment on fused silica. To study the deformation and thermal properties of RDX at the atomistic scale, we have performed nanoindentation simulations on a (100) α-RDX crystal surface using reactive force field. We have observed localized melting and decomposition of RDX molecular fragments due to significant increase in temperature under the indenter. At the maximum indentation depth, we have found: (1) molecules above the ! 210 ( ) plane have displaced more than 5 Å from their original positions; (2) molecules with their positions below the ! 210 ( ) plane have remained intact. Simulation also shows that in order to release the excess stress, the fragmented RDX molecules diffuse from the substrate and walk on the indenter surface. 102 REFERENCES [1] Movie files of animations of the bond-switching events can be obtained from the EPAPS home page (www.aip.org/pubservs/epaps.html). [2] G. C. Abell, Physical Review B 31, 6184 (1985). [3] A. Acrivos, Annals of the New York Academy of Sciences 404, 1 (1983). [4] M. P. Allen, and D. J. Tildesley, Computer Simulation of Liquids ( Oxford University Press, Oxford, 1987). [5] N. L. Allinger, Y. H. Yuh, and J. H. Lii, Journal of the American Chemical Society 111, 8551 (1989). [6] N. L. Allinger, F. B. Li, and L. Q. Yan, Journal of Computational Chemistry 11, 848 (1990). [7] N. L. Allinger et al., Journal of Computational Chemistry 11, 868 (1990). [8] H. C. Andersen, Journal of Chemical Physics 72, 2384 (1980). [9] T. L. Anderson, Fracture Mechanics Fundamentals and Applications (CRC Press, Florida, 1995). [10] A. S. Argon, J. Im, and R. Safoglu, Metallurgical Transactions A 6, 825 (1975). [11] I. Arias et al., Computer Methods in Applied Mechanics and Engineering 196, 3833 (2007). [12] R. W. Armstrong, and W. L. Elban, in Proc. Seventh Symp. (Int.) on DetonationAnnapolis, MD, p. 976. [13] R. W. Armstrong, and W. L. Elban, Materials Science and Engineering a- Structural Materials Properties Microstructure and Processing 111, 35 (1989). [14] R. W. Armstrong, Journal De Physique Iv 5, 89 (1995). [15] R. W. Armstrong et al., Thermochimica Acta 384, 303 (2002). [16] R. W. Armstrong, and W. L. Elban, Materials Science and Technology 22, 381 (2006). [17] R. W. Armstrong, W. Arnold, and F. J. Zerilli, Metallurgical and Materials Transactions a-Physical Metallurgy and Materials Science 38A, 2605 (2007). 103 [18] M. F. Ashby, and S. D. Hallam, Acta Metallurgica 34, 497 (1986). [19] M. F. Ashby, and C. G. Sammis, Pure and Applied Geophysics 133, 489 (1990). [20] N. P. Bansal, and D. R.H, Handbook of glass properties (Academic Press, Orlando, 1986). [21] G. M. Bartenev, The structure and mechanical properties of inorganic glasses (Wolters-Nordhoff Publ., Groningen, 1970). [22] B. D. Beake, and J. F. Smith, Philosophical Magazine a-Physics of Condensed Matter Structure Defects and Mechanical Properties 82, 2179 (2002). [23] F. M. Beremin, Metall Trans A 12, 723 (1981). [24] N. Bernstein, and E. Kaxiras, Physical Review B 56, 10488 (1997). [25] N. Binggeli, and J. R. Chelikowsky, Nature 353, 344 (1991). [26] D. Bonamy et al., International Journal of Materials & Product Technology 26, 339 (2006). [27] E. Bouchaud, Surface Review and Letters 10, 797 (2003). [28] D. R. Bowler et al., Physica Status Solidi B-Basic Solid State Physics 243, 989 (2006). [29] R. C. Bradt, and A. G. Evans, Fracture mechanics of ceramics (Plenum Press, New York, London, 1974). [30] D. W. Brenner, Physical Review B 42, 9458 (1990). [31] T. B. Brill et al., Decomposition, Combustion, and Detonation Chemistry of Energetic Materials (Materials Research Society, Pittsburgh, 1996). [32] L. M. Brown, and J. D. Embury, in Third International Conference on Strength of Metals and Alloys (Institute of Metals, London, 1973), p. 164. [33] B. J. Budianski, J. W. Hutchinson, and S. Slutsky, in Mechanics of Solids. The Rodney Hill 60th Anniversary Volume (Pergamon Press, Oxford, 1982). [34] F. Celarie et al., Physical Review Letters 90 (2003). [35] W. N. Chen, and G. Ravichandran, Journal of the Mechanics and Physics of Solids 45, 1303 (1997). 104 [36] K. S. Cheung, and S. Yip, Journal of Applied Physics 70, 5688 (1991). [37] A. C. F. Cocks, and M. F. Ashby, Metal Science 14, 395 (1980). [38] A. C. F. Cocks, and M. F. Ashby, Progress in Materials Science 27, 189 (1982). [39] R. F. Cook, and G. M. Pharr, Journal of the American Ceramic Society 73, 787 (1990). [40] R. F. Cook, and E. G. Liniger, Journal of the American Ceramic Society 76, 1096 (1993). [41] B. Cotterell, and J. R. Rice, International Journal of Fracture 16, 155 (1980). [42] T. B. COX, and J. R. LOW, Metallurgical Transactions 5, 1457 (1974). [43] E. Deeg, Glastech. Ber. 31, 1 (1958). [44] C. Donati et al., Physical Review Letters 80, 2338 (1998). [45] E. Dupree, and R. F. Pettifer, Nature 308, 523 (1984). [46] J. L. Fattebert, and F. Gygi, Computer Physics Communications 162, 24 (2004). [47] J. L. Fattebert, and F. Gygi, Physical Review B 73 (2006). [48] H. G. Gallagher et al., Philosophical Transactions of the Royal Society of London Series a-Mathematical Physical and Engineering Sciences 339, 293 (1992). [49] T. R. Gibbs, and A. Popolato, LASL Explosive Property Data (1980). [50] S. Goedecker, Reviews of Modern Physics 71, 1085 (1999). [51] V. M. Goldschmidt, Geochemische Verteilungsgesetze der Elemente. 8, 137 (1926). [52] H. P. Grace, Chemical Engineering Communications 14, 225 (1982). [53] A. A. Griffith, Philosophical Transaction 221A, 163 (1920). [54] J. P. Guin, and S. M. Wiederhorn, International Journal of Fracture 140, 15 (2006). [55] J. T. Hagan, and M. M. Chaudhri, Journal of Materials Science 12, 1055 (1977). [56] P. J. Halfpenny, K. J. Roberts, and J. N. Sherwood, Journal of Crystal Growth 65, 524 (1983). 105 [57] C. A. Harper, Handbook of ceramics, glasses, and diamonds (McGraw-Hill Companies, New York, 2001). [58] E. J. Hinch, and A. Acrivos, Journal of Fluid Mechanics 91, 401 (1979). [59] E. J. Hinch, and A. Acrivos, Journal of Fluid Mechanics 98, 305 (1980). [60] K. Hirao, and M. Tomozawa, Journal of the American Ceramic Society 70, 497 (1987). [61] W. G. Hoover, Physical Review A 31, 1695 (1985). [62] C. E. Inglis, Transaction of the Institute of Naval Architects 55, 219 (1913). [63] D. Joseph, Journal of Fluid Mechanics 366, 367 (1998). [64] R. K. Kalia et al., Int J Fracture 121, 71 (2003). [65] B. B. Karki et al., Physical Review B 55, 3465 (1997). [66] H. Kimura et al., Applied Physics Letters 78, 604 (2001). [67] M. M. Kuklja, and A. B. Kunz, Journal of Physics and Chemistry of Solids 61, 35 (2000). [68] M. M. Kuklja, and A. B. Kunz, Journal of Applied Physics 89, 4962 (2001). [69] C. R. Kurkjian, Strength of inorganic glass (Plenum, New York, 1985). [70] K. Y. Kwon et al., Abstracts of Papers of the American Chemical Society 229, U746 (2005). [71] K. Y. Kwon et al., Physical Review Letters 95 (2005). [72] Y. D. Lanzerotti, J. Sharma, and R. W. Armstrong, Metallurgical and Materials Transactions a-Physical Metallurgy and Materials Science 35A, 2675 (2004). [73] B. Lawn, Fracture of Brittle Solids (Cambridge U. Press, New York, 1993). [74] B. J. Lee, and M. E. Mear, Journal of the Mechanics and Physics of Solids 47, 1301 (1999). [75] S. Lee, and G. Ravichandran, Optics and Lasers in Engineering 40, 341 (2003). [76] S. Lee, and G. Ravichandran, Engineering Fracture Mechanics 70, 1645 (2003). 106 [77] D. H. Liebenberg, R. W. Armstrong, and J. J. Gilman, Structure and Properties of Energetic Materials (Materials Research Society, Pittsburgh, 1993), Vol. 296. [78] Z. Lu et al., Physical Review Letters 95, 135501 (2005). [79] M. Manga et al., Journal of Volcanology and Geothermal Research 87, 15 (1998). [80] D. M. Marsh, Proceeding of Royal Society 282A, 33 (1964). [81] G. J. Martyna, M. L. Klein, and M. Tuckerman, Journal of Chemical Physics 97, 2635 (1992). [82] G. J. Martyna, D. J. Tobias, and M. L. Klein, Journal of Chemical Physics 101, 4177 (1994). [83] G. J. Martyna et al., Molecular Physics 87, 1117 (1996). [84] S. L. Mayo, B. D. Olafson, and W. A. Goddard, Journal of Physical Chemistry 94, 8897 (1990). [85] F. A. McClintock, Mechanical Engineering 90, 56 (1968). [86] P. J. Miller, C. S. Coffey, and V. F. Devost, Journal of Applied Physics 59, 913 (1986). [87] K. Miyake et al., Japanese Journal of Applied Physics 43, 4602 (2004). [88] S. C. Moss, and D. Price, Physics of Disordered Materials (Plenum, New York, 1985). [89] N. F. Mott, Philosophical Magazine B-Physics of Condensed Matter Statistical Mechanics Electronic Optical and Magnetic Properties 56, 257 (1987). [90] A. Nakano, R. K. Kalia, and P. Vashishta, Journal of Non-Crystalline Solids 171, 157 (1994). [91] K. I. Nomura et al., Comput Phys Commun 178, 73 (2008). [92] S. Nose, Journal of Chemical Physics 81, 511 (1984). [93] S. Nose, Molecular Physics 52, 255 (1984). [94] W. C. Oliver, and G. M. Pharr, Journal of Materials Research 7, 1564 (1992). [95] M. Parrinello, and A. Rahman, Journal of Applied Physics 52, 7182 (1981). 107 [96] M. Parrinello, and a. Rahman, Journal of Applied Physics 52, 7182 (1981). [97] M. Parrinello, and A. Rahman, Journal of Chemical Physics 76, 2662 (1982). [98] A. Perriot et al., Journal of the American Ceramic Society 89, 596 (2006). [99] K. W. Peter, Journal of Non-Crystalline Solids 5, 103 (1970). [100] G. M. Pharr, W. C. Oliver, and F. R. Brotzen, Journal of Materials Research 7, 613 (1992). [101] G. M. Pharr, Materials Science and Engineering a-Structural Materials Properties Microstructure and Processing 253, 151 (1998). [102] J. M. Rallison, Annual Review of Fluid Mechanics 16, 45 (1984). [103] S. Ramaswamy, and N. Aravas, Computer Methods in Applied Mechanics and Engineering 163, 11 (1998). [104] S. Ramdas, J. M. Thomas, and M. J. Goringe, Journal of the Chemical Society- Faraday Transactions Ii 73, 551 (1977). [105] A. K. Rappe et al., Journal of the American Chemical Society 114, 10024 (1992). [106] J. R. Rice, and D. M. Tracey, Journal of the Mechanics and Physics of Solids 17, 201 (1969). [107] C. L. Rountree et al., Annual Review of Materials Research 32, 377 (2002). [108] E. T. Seppala, J. Belak, and R. E. Rudd, Physical Review B 71 (2005). [109] J. Sharma et al., Applied Physics Letters 78, 457 (2001). [110] F. Shimojo et al., Computer Physics Communications 167, 151 (2005). [111] G. S. Smith et al., Acta Materialia 49, 4089 (2001). [112] F. H. Stillinger, and T. A. Weber, Physical Review B 31, 5262 (1985). [113] H. Stone, Annual Review of Fluid Mechanics 26, 65 (1994). [114] A. Strachan et al., Physical Review Letters 91 (2003). [115] W. B. Streett, D. J. Tildesley, and G. Saville, Molecular Physics 35, 639 (1978). [116] K. H. Sun, Journal of American Ceramic Society 30, 277 (1947). 108 [117] G. Tammann, Der Glaszustand (Voss, Leipzig, 1933). [118] G. I. Taylor, Proceedings of the Royal Society of London Series a-Containing Papers of a Mathematical and Physical Character 138, 41 (1932). [119] J. Tersoff, Physical Review B 37, 6991 (1988). [120] M. Tuckerman, B. J. Berne, and G. J. Martyna, Journal of Chemical Physics 97, 1990 (1992). [121] V. Tvergaard, Y. Huang, and J. W. Hutchinson, European Journal of Mechanics a-Solids 11, 215 (1992). [122] V. Tvergaard, and a. Needleman, International Journal of Solids and Structures 32, 1063 (1995). [123] D. R. Uhlmann, and N. J. Kreidl, Glass: Science and technology (Academic Press, New York, London. [124] A. C. T. van Duin et al., Journal of Physical Chemistry A 105, 9396 (2001). [125] P. Vashishta et al., Solid State Ionics 40-1, 175 (1990). [126] P. Vashishta et al., Physical Review B 41, 12197 (1990). [127] C. E. Weir, and L. Shartsis, Journal of American Ceramic Society 38, 299 (1955). [128] S. Yip, Handbook of Materials Modeling (Springer, Berlin, 2005). [129] W. H. Zachariasen, Journal of American Chemical Society 54, 3841 (1932).
Abstract (if available)
Abstract
Multimillion-to-billion molecular dynamics (MD) simulations are carried out to study atomistic mechanisms of deformation, damage and failure in silica glass and energetic materials. The simulations are based on experimentally validated interatomic potentials and employ highly efficiently algorithms for parallel architectures.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Shock sensitivity of energetic material and nanometric damage mechanisms in silica glass
PDF
Dynamics of wing cracks and nanoscale damage in silica glass
PDF
Nanomaterials under extreme environments: a study of structural and dynamic properties using reactive molecular dynamics simulations
PDF
Molecular dynamics simulations of nanostructures
PDF
Hypervelocity impact damage in alumina
PDF
Nanoindentation of silicon carbide and sulfur-induced embrittlement of nickel
PDF
Shock-induced poration, cholesterol flip-flop and small interfering RNA transfection in a phospholipid membrane: multimillion atom, microsecond molecular dynamics simulations
PDF
Molecular dynamics study of energetic materials
PDF
Large-scale molecular dynamics simulations of nano-structured materials
PDF
Simulation and machine learning at exascale
PDF
Multi-scale quantum dynamics and machine learning for next generation energy applications
PDF
Controlling electronic properties of two-dimensional quantum materials: simulation at the nexus of the classical and quantum computing eras
PDF
Neural network for molecular dynamics simulation and design of 2D materials
PDF
Physics informed neural networks and electrostrictive cavitation in water
PDF
Atomistic modeling of the mechanical properties of metallic glasses, nanoglasses, and their composites
PDF
High throughput computational framework for synthesis and accelerated discovery of dielectric polymer materials using polarizable reactive molecular dynamics and graph neural networks
PDF
Modeling of bicontinuous nanoporous materials and investigation of their mechanical and hydraulic properties
PDF
Theoretical studies of lipid bilayer electroporation using molecular dynamics simulations
PDF
Quantum molecular dynamics and machine learning for non-equilibrium processes in quantum materials
PDF
Dynamics of water in nanotubes: liquid below freezing point and ice-like near boiling point
Asset Metadata
Creator
Chen, Yi-Chun
(author)
Core Title
Multimillion-to-billion atom molecular dynamics simulations of deformation, damage, nanoindentation, and fracture in silica glass and energetic materials
School
College of Letters, Arts and Sciences
Degree
Doctor of Philosophy
Degree Program
Physics
Publication Date
07/31/2008
Defense Date
05/02/2008
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
damage,deformation,fracture,nonoidentation,OAI-PMH Harvest
Language
English
Advisor
Kalia, Rajiv K. (
committee chair
), Bickers, Nelson Eugene, Jr. (
committee member
), Lu, Jia Grace (
committee member
), Nakano, Aiichiro (
committee member
), Vashishta, Priya (
committee member
)
Creator Email
yichunc@usc.edu
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-m1464
Unique identifier
UC1128655
Identifier
etd-Chen-20080731 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-91543 (legacy record id),usctheses-m1464 (legacy record id)
Legacy Identifier
etd-Chen-20080731.pdf
Dmrecord
91543
Document Type
Dissertation
Rights
Chen, Yi-Chun
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Repository Name
Libraries, University of Southern California
Repository Location
Los Angeles, California
Repository Email
cisadmin@lib.usc.edu
Tags
deformation
fracture
nonoidentation