UFDC Home  myUFDC Home  Help 



Full Text  
PAGE 1 1 THERMAL TRANSPORT IN URAN IUM DIOXIDE AND DIAMOND BY ATOMIC LEVEL SIMULATIONS By TAKU WATANABE A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLOR IDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2008 PAGE 2 2 2008 Taku Watanabe PAGE 3 3 To my family PAGE 4 4 ACKNOWLEDGMENTS Since the time I arrived at the University of Florida, I have been fortunate to meet many people for whom I am deeply thankful. These people have made my experience at the University of Florida extraordinary. I cannot sufficiently thank supervisory committee chair, Prof. Phillpot, for his enthusiasm of science and education and patience to guide his student. I respect him as a scientist, educator, and individual He has been always willing to discuss about almost anything with a great sens e of humor. He had a positive e ffect on my character. I also thank to Prof. Schelling in the University of Cent ral Florida for helping me carrying out some of the key work I have done during my Ph.D program. He was always willing to help and I have learned so much from his deep insight of the subj ect. I am also grateful to Prof. Sinnott, Prof. Tulenko, Prof. Nino, Prof. Jones, and Prof. Anghai e for their guidance and willingness to help students. I was fortunate to have the experience of visiting Los Alamos National Laboratory in summer of 2006. The people in MST8 impressed me as scientists, mentors, and friends. Drs. Srivilliputhr, Uberuaga, and Stan ek were especially helpful a nd made my time there productive and enjoyable. It has been a great pleasure working with th e people in our com putational group in the Materials Science Department in the University of Florida. The group has been very active and supportive of each other, and it had positive impact on my research and life outside of the school. Finally, I thank to my family in Japan for th eir support and giving me the opportunity for the education and experience. Part of this work was supported by DOENERI Award DEFC0705ID14649. PAGE 5 5 TABLE OF CONTENTS page ACKNOWLEDGMENTS...............................................................................................................4 LIST OF TABLES................................................................................................................. ..........8 LIST OF FIGURES................................................................................................................ .........9 ABSTRACT....................................................................................................................... ............14 CHAPTER 1 GENERAL INTRODUCTION..............................................................................................15 Challenges in the Materials Science of Heat Transfer............................................................15 Defects in Crystalline Solids.................................................................................................. 17 Uranium Dioxide................................................................................................................ ....20 Diamond........................................................................................................................ .........23 Objective and Scope of Study.................................................................................................24 2 INTRODUCTION TO THERMAL TRANSPORT IN MATERIALS..................................34 Basic Concepts of Thermal Transport....................................................................................34 Three Modes of Heat Transfer Conduction, Radiation, Convection....................................35 Phonon Heat Transfer in Electronic Insulators.......................................................................37 Point Imperfections............................................................................................................ .....43 GBs and Surfaces............................................................................................................... .....45 Thermal Transport in Disordered Solids................................................................................46 3 SIMULATION METHODS...................................................................................................52 Introduction................................................................................................................... ..........52 Numerical Integration.......................................................................................................... ...53 Periodic Boundary Condition.................................................................................................53 Interatomic Potentials......................................................................................................... ....54 Longrange Interaction....................................................................................................54 Ewald sum................................................................................................................55 Direct summation.....................................................................................................56 Shortrange Interactions..................................................................................................57 Thermostat..................................................................................................................... .........58 Constant Pressure Algorithm..................................................................................................59 Nonequilibrium Thermal Conductivity Simulations..............................................................61 4 THERMAL TRANSPORT IN SINGLE CRYSTAL UO2.....................................................68 Molecular Dynamics Simulation of UO2................................................................................68 PAGE 6 6 Potential Models............................................................................................................... ......68 Thermal Expansion.............................................................................................................. ...70 Thermal Conductivity of Single Crystal UO2.........................................................................72 Discussion..................................................................................................................... ..........75 5 THERMAL TRANSPORT IN POLYCRYSTALLINE UO2.................................................84 Introduction................................................................................................................... ..........84 Structure of Model Polycrystal...............................................................................................84 Thermal Conductivity of Polycrystalline UO2.......................................................................86 Discussion..................................................................................................................... ..........89 6 THERMAL TRANSPORT IN UO2x.....................................................................................98 Introduction................................................................................................................... ..........98 Point Defects in UO2.......................................................................................................98 Thermal Transport in UO2 with Point Defects..............................................................101 Molecular Dynamics Simulations.........................................................................................103 Preparation of the Structures.........................................................................................104 Chemical Expansion of UO2x.......................................................................................105 Thermal Expansions of UO2x.......................................................................................106 Thermal Conductivity of UO2+ x.....................................................................................107 Thermal Conductivity of UO2x.....................................................................................109 Isotopic Effect...............................................................................................................1 09 Lattice Dynamics Calculations......................................................................................110 Discussion..................................................................................................................... ........112 7 METHODOLOGY FOR THE SIMULATION OF RADIATION DAMAGE IN UO2.......138 Introduction................................................................................................................... ........138 Background..................................................................................................................... ......139 Type of Radiation..........................................................................................................139 Fission and Radiation....................................................................................................140 Fission Products and th e Type of Stopping...................................................................141 Previous Experimental Studies on Radiation Damage in UO2......................................142 Previous Simulation Studies on Radiation Damage in UO2.................................................144 Implementation of the Molecular Dynamics Simulation Code............................................145 Universal ZieglerBiersakLittmark Potential...............................................................145 Variable Time Step Algorithm......................................................................................149 Thermostat..................................................................................................................... 153 Defect Identification Method........................................................................................156 Radiation Damage Simulations............................................................................................157 Time Step Size...............................................................................................................15 7 Berendsen Time Constant..............................................................................................159 Statistics..................................................................................................................... ....160 Active Volume Effect....................................................................................................161 Discussion..................................................................................................................... ........162 PAGE 7 7 8 INTERFACIAL THERMAL TRANSPORT IN DIAMOND..............................................187 Introduction................................................................................................................... ........187 Single Crystal Diamond........................................................................................................1 87 Interfacial Conductance of (001) Twist GBs........................................................................190 Interfacial Thermal Conductance: Diamond vs. Silicon.......................................................196 Extended ReadShockley Model For Interfacial Thermal Conductance..............................197 Discussion..................................................................................................................... ........200 9 SIMULATIONS OF THERMAL TRANSPORT IN NANOCRYSTALLINE DIAMOND AND COMPARISON WITH EXPERIMENTAL RESULTS.........................222 Introduction................................................................................................................... ........222 Experimental Approach........................................................................................................22 4 Experimental Results........................................................................................................... .225 Simulation Results............................................................................................................. ...227 Discussion..................................................................................................................... ........229 10 CONCLUSIONS................................................................................................................. .236 APPENDIX: LATTICE DYNAMICS.........................................................................................240 LIST OF REFERENCES............................................................................................................. 243 BIOGRAPHICAL SKETCH.......................................................................................................269 PAGE 8 8 LIST OF TABLES Table page 21. Values of the Debye temperature and the maximum in the experimental thermal conductivities of selected materials...................................................................................49 41. Parameters of interatomic potentials..................................................................................... 77 42. Lattice parameter and el astic constants at 300 K..................................................................78 61. Experimental and DFT values of the formation energies....................................................115 62. Potential parameters for th e shortrange interactions in UO2 x...........................................116 71. Spline parameters and connection point s for UU, OO, and UO interactions..................164 72. The ZBL screening function parameters.............................................................................165 81. Summary of the diamond (001) GBs properties..................................................................201 82. Comparison of the GB energies and coordination of the atoms of 29(001) =43.6 GB in diamond and Si............................................................................................................20 2 PAGE 9 9 LIST OF FIGURES Figure page 11. Point defects in crystalline solid....................................................................................... .....25 12. Pure relaxed single twist grain b oundary (GB) between two cubic lattices..........................26 13. An asymmetric tilt GB between two cubic lattices...............................................................27 14. GB energy as a function of rotation angle in <100> (a,b) and <110> (c,d) symmetric tilt GBs of aluminum...............................................................................................................2 8 15. Illustration of th e lattice structure in 5 (001) =36.5 GB. ...............................................29 16. Fluorite structure of UO2 unit cell. Blue and red denote uranium and oxygen atoms, respectively................................................................................................................... .....30 17. Phase diagram of UO system.............................................................................................. .31 18. Fission product yield as a function of at omic number in a pressurized water reactor (PWR) after 2.9 % burnup.................................................................................................32 19. Unit cell structure of diamond........................................................................................... ....33 21. Temperature dependence of the ther mal conductivity of selected materials.........................50 22. Temperature dependence of the thermal conductivity of Ge................................................51 31. The root mean square of the total ener gy as a function of time for velocity Verlet (circles), 4th order Gear (squares), 5th order Gear (triangles), and 6th order Gear (diamond) algorithms.........................................................................................................63 32. Two dimensional periodic system in simulation...................................................................64 33. A system of point charge s as a sum of screened point charges and smeared charge density with opposite sign..................................................................................................65 34. Electrostatic energy by the direct summation method...........................................................66 35. Simulation cell setup for the direct method for simulating thermal conductivity.................67 41. Normalized lattice parameter as a function of temperature from experiment 3, Busker and Yamada potentials.......................................................................................................79 42. Specific heat of UO2 from GULP, MD and experiment........................................................80 43. Size dependence of thermal conductivity for both the Busker (a) and Yamada (b) potentials..................................................................................................................... .......81 PAGE 10 10 44. Thermal conductivity of UO2 from the compilation of expe rimental data by Fink and simulations using the Busker and Yamada potentials (a) before and (b) after the anharmonic correction.......................................................................................................82 45. The same thermal conductivity data as in Figure. 44(b), but in loglog plot.......................83 51. Pair distribution function of the polycrystalline structure ...................................................91 52. The final polycrystalline structure used in the thermal conductiv ity calculations ..............92 53. Thermal conductivity of 3.8 nm grain polyc rystalline UO2 from simulation for the two potentials..................................................................................................................... .......93 54. Average thermal conductance for the [001] tilt GBs from 3.8 nm grain polycrystalline UO2............................................................................................................................... .....94 55. Kapitza length from 3.8 nm grain polycrystalline UO2 simulations.....................................95 56. Grain size dependence of the th ermal conductivity of polycrystalline UO2 using Yamadas potential............................................................................................................9 6 57. Interfacial conductances of selected interfaces as functio n of temperature..........................97 61. Normalized Frenkel and antiFrenkel pair formation energies by empirical potentials......117 62. The 2:2:2 defect cluster in UO2...........................................................................................118 63. Pair distribution fu nction between UO in UO2.25...............................................................119 64. Pair distribution fu nction between UO in UO1.80...............................................................120 65. Normalized lattice parameter of UO2+x as a function of defect concentration....................121 66. Normalized lattice parameter of UO2x as a function of defect concentration.....................122 67. Thermal expansion of UO2+x...............................................................................................123 68. Thermal expansion coefficient of UO2+x as a function of x ................................................124 69. Thermal expansion of UO2x................................................................................................125 610. Inverse thermal conductivity of UO2+x with various concentrations of oxygen interstitials as a function of simulation cell length at 800 K............................................126 611. Inverse thermal conductivity of UO2+x with various concentrations of oxygen interstitials as a function of simulation cell length at 1600 K..........................................127 612. Bulk thermal conductivity of UO2+x as a function of degree of offstoichiometry without anharmonic correction at 800 and 1600 K..........................................................128 PAGE 11 11 613. Bulk thermal conductivity of UO2+x as a function of degree of offstoichiometry with anharmonic correction at 800 and 1600 K.......................................................................129 614. Inverse thermal conductivity of UO2x with various concentrations of oxygen vacancies as a function of si mulation cell length at 800 K..............................................130 615. Inverse thermal conductivity of UO2x with various concentrations of oxygen vacancies as a function of si mulation cell length at 1600 K............................................131 616. Bulk thermal conductivity of UO2x as a function of degree of offstoichiometry at 800 K and 1600 K................................................................................................................... 132 617. Thermal conductivity of UO2x as a function of degree of offstoichiometry at 800 K (a) and 1600 K (b) for x 0.025........................................................................................133 618. Phonon Density of states of UO2 x for 0.020 x 0.125....................................................134 619. Participation ratios of UO2 x for 0.020 x 0.125..............................................................135 620. Projections of the polarization vectors of vibrational modes in UO2x at 1.5 THz.........136 621. Projections of the polarization vectors of vibrational modes in UO2 x at 10 THz............137 71. Potential kinetic, and total energy per at om from radiation damage simulations with various fixed time step size..............................................................................................166 72. Potential energy of OO interactions as a function of the interatomic distance. ...............167 73. Trajectory of the O atom placed at 0.1 away from another O atom at the origin............168 74. The calculation of the PKA trajectory in tw o ways: one large step and two half steps......169 75. Flow chart of the determination of the variable time step algorithm..................................170 76. Trajectory of 1 keV U PKA atom along [ 100] direction in 101010 single crystal UO2.171 77. Trajectory of 80 keV U PKA along [100] direction in 1010 10 single crystal UO2........172 78. Evolution of the time step size in the 80 keV U PKA simulation.......................................173 79. The configuration of the MD simulation cell used in ou r radiation damage simulations...174 710. Temperature evolution of 101010 unit cells single crystal UO2 during equilibration with the Berendsen thermostat with =400 fs..................................................................175 711. Total energy per atom of the equilibrated single crystalline UO2 at 300 K......................176 712. Standard deviation of the kinetic, potential and total energies per atom as a function of time constant, from the equilibrated single crystalline UO2 at 300 K..........................177 PAGE 12 12 713. Defect identification met hod in crystalline structures.......................................................178 714. Effect of the time step size on radi ation damage simulation in NVE ensemble...............179 715. Number of VU, UI, VO, and OI produced after the 1 keV U PKA radiation damage in a single crystal UO2............................................................................................................180 716. Effect of time step size in ra diation damage simulation in pNVE....................................181 717. Temperature of the active region duri ng the radiation damage simulations with different Berendsen time constants..................................................................................182 718. Total energy per atom during the radi ation damage simulations with different Berendsen time constant..................................................................................................183 719. Number of uranium (top) and oxygen ( bottom) vacancies produced in radiation damage simulations as a function of the Berendsen time constant..................................184 720. Number of uranium (top) and oxygen ( bottom) vacancies produced in radiation damage simulations with e quivalent PKA directions......................................................185 721. The number of U and O vacancies produ ced by the radiation damage simulations in 202020 and 303030 unit cells active volume..........................................................186 81. Lattice parameter of a single crysta l diamond as a function of temperature.......................203 82. Temperature profile of 44400 unit cells single crystal diamond after 4 ps.....................204 83. Time evolution of the temperature at the edge ( Lz=0) and center ( z = Lz/2) of the simulation cell................................................................................................................ ..205 84. Slopes of the linear region of the te mperature profile shown in Figure 82........................206 85. Slopes of the linear region of the te mperature profile shown in Figure 84........................207 86. Thermal conductivity versus temperat ure in both linear a nd loglog plots.........................208 87. Bonding of atoms in diamond structure along (a) [111], (b) [100] and (c) [110] directions..................................................................................................................... .....209 88. Thermal conductivity as a func tion of heat current density................................................210 89. GB energy as a function of the number of atoms removed.................................................211 810. A typical temperature profile of a system with GBs, in this case (001) =43.60 29 symmetric twist boundaries.............................................................................................212 811. Linear temperature gradient from the temperature profile of the 29 (001) =43.60 GB............................................................................................................................. .......214 PAGE 13 13 812. Size dependence of the thermal conductance of the 29 =43.60 GB fitted to a shifted exponential...........................................................................................................2 15 813. Temperature dependence of th e thermal conductance of the (001) 29 =43.60 GB.....216 814. Thermal conductance as a function of tw ist angle for diamond (001) symmetric twist GBs............................................................................................................................ ......217 815. Kapitza length in units of their respec tive lattice parameters as a function of twist angle for both diamond and Si at 1000 K........................................................................218 816. Cross section of 29(001) GB in diamond (a) and Si (b)................................................219 817. GB energy (symbols) and the fit to the extended ReadShockley model (solid line).......220 91. Compilation of the diamond ther mal conductivity versus grain size d ...............................232 92. UNCD thermal conductivity as a function of film thickness..............................................233 93. Microstructure and temperatur e profile of polycrystallin diamond.....................................234 94. Temperature profile through a diamond 29 =43.60 twist GB determined by simulation..................................................................................................................... ....235 PAGE 14 14 Abstract of Dissertation Pres ented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy THERMAL TRANSPORT IN URAN IUM DIOXIDE AND DIAMOND BY ATOMIC LEVEL SIMULATIONS By Taku Watanabe May 2008 Chair: Simon R. Phillpot Major: Materials Science and Engineering Thermal transport properties of UO2 and diamond are investig ated using atomic level simulation techniques. In partic ular, the effects of point defect s and GBs are investigated. A nonequilibrium molecular dynamics t echnique is used to calculat e the thermal c onductivities of UO2 and diamond. The correct treatment of anharm onicity is found to be important for accurate prediction of thermal transport in crystallin e solid. The thermal tr ansport properties of polycrystalline UO2 are investigated. Three effective medium approaches are compared and used to derive the interfacial conductance from our simulation results. The temperature and grain size dependence of the interfacial conductance is also characterized. Defect concentration dependence of the thermal conductivity of UO2 x is obtained. In the effort to understand the nature of vibrational modes in UO2 with point defects, the theory of highly disordered solid is applied. To a ddress the defect pro duction and annihilation mechanisms in UO2, a radiationdamage simulation code has been developed and validated. Thermal transport properties of diamond GB s are investigated and characterized. Connections to the experimental wo rk or collaborators are made. PAGE 15 15 CHAPTER 1 GENERAL INTRODUCTION Challenges in the Materials Science of Heat Transfer Heat transfer is a funda mental mode of energy transport. It is so fundamental that it can be observed at the atomic and molecular levels and al so at the astronomical sc ale. Controlling heat transfer has been technologically important for mo re than a century. Thus, relentless efforts have been made to tailor the thermal transport properties of materials for specific applications. This is often very challenging because of the complex ity of the transport mechanisms. Heat can propagate by means of conduction, convection, and radiation, and any co mbinations of these three mechanisms. In addition, the range of ther mal conductivities of materi als is quite narrow. For instance, the electrical c onductivity accessible in material s is roughly 20 orders of magnitude; however, the thermal conductivity of materials ranges only over about 5 orders of magnitude.1, 2 Despite the challenge, the dema nd for better management of the heat transfer has been continuously intensifying. In the effort to gain better control, one must look at the materials on both ends of the spectrum; namely, very hi gh and low thermal conductivity materials. UO2 is a good example of a low thermal conductiv ity ceramic material. It is al so the standard nuclear fuel for light water reactors for power generation. Its the thermal conductivity is roughly 27 W/mK depending on its temperature and microstructure3. In a nuclear reactor, the energy released from the fission reactions is converted to the thermal energy. For efficient transfer of the thermal energy, it is desirable for the fuel to have a high thermal conductivity. Improving the thermal conductivity of UO2 is a significant challenge; for example, see the literatures on cermet materials.4, 5 PAGE 16 16 An example of the other end of the spectru m is diamond. Diamond has the highest thermal conductivity among the natura lly occurring materials at room temperature.2, 6 There are a number of applications which require very hi gh thermal conductivity, incl uding heat spreader in microelectronics, coating for high wear machin ing devices, and high power IR laser windows. Because of this and its other unique propert ies, diamond, is being used in numerous applications.6, 7 Both UO2 and diamond are electronic insulators, with band gaps of ~2 eV8, 9 and ~5.5 eV1, 10,respectively. Thus, the main m echanism of the transport of h eat is by the vibration of the crystalline lattice. The quan tized lattice vibrations are call ed phonons and they are the main carriers of thermal energy in these materials. Thermal conductivity in an electronic insulator is determined by the scattering of phonons with other phonons and with defects. Thus the means to control phonon mediated heat transf er is to control point defects, grains, and addition of other materials. The detailed scattering mechanis ms of phonon with other phonons and defects are quite complicated. Although there are some anal ytical models available, there are severe approximations made and the re sults are often in poor agreement with the observed data. Difficulties are not only in the theoretical understandi ng, but also in experiments. In experiments, numerous effects coexist in a single sample, as a result not straightfo rward to separate the mechanisms in clear manner. In order to provide the necessa ry background to the work discu ssed in this thesis, in the following sections, the fundamental s of the defects in crystalline solids are introduced. In addition, the basic properties of UO2 and diamond are presented. In particular, point defects and GBs present in these two materials are introduced. PAGE 17 17 Defects in Crystalline Solids Solid materials with crystalline structure can contain a number of structural defects including point defects, disloc ations, GBs, pores, second phases, and surfaces. In the work discussed here, we describe the effect of point defects and GBs on the thermal transport properties. Thus in this secti on, the basic features of point de fects and GBs are introduced. Point defects in crystalline so lid can be divided into two ma jor categories: intrinsic and extrinsic point defects. Intrinsic point defects are defects of the original matrix material, including vacancies, interstitial s, Frenkel pairs, and Schottky defects: see Figure 11 for an illustration of these defects. The probability that a vacancy will form in thermal equilibration is proportional to the Boltzmann factor, exp(Ev/ kBT ), where Ev is the vacancy formation energy, kB is the Boltzmann constant, and T is the temperature. Thus the equilibrium number, nv, of vacancies in a crystal of N atoms is given by, T k E n N nB v v v/ exp (11) A typical value of Ev is of the order of ~1 eV. Assuming T =1000 K and Ev=1 eV, this yields nv/ N 105. The equilibrium number of interstitials, Frenkel pairs, and Schottky defects can be obtained in a similar fashion. Extrinsic point defects are atom species differe nt from the perfect solid itself and can take the form of interstitials and substitutional defects. Intrinsic and extrinsic defects often coexist in the same solid and they are extremely useful in many applications, such as ion transport in doped oxides. Intr insic defects such as vacancies and interstitials can be created by doping of extrinsic defect. For example, in UO2, doping Gd2O3 can create vacancies in the crystal: PAGE 18 18 O O U UOV O G O Gd 3 222 3 2. (12) Gd is commonly doped in UO2 to control the activity of the fuel in boiling water reactors11, 12. The vacancies created by the Gd2O3 doping are known to increase th e electrical conductivity of the UO2 fuel13, 14, but the reduce thermal conductivity by enhanced phonondefect scatterings 15, 16. The theory of phonon scattering mech anisms is discussed in Chapter 2. There are a number of models proposed to describe the st ructure and prope rties of GBs.17 We will cover discuss only rather simple grain boundary (GB) structures, namely, twist and tilt GBs in some detail since they are the most relevant to the rest of this thesis. However, a grain boundary in general can also be a combination of tilt or twist. (For more details, see reference 18.) Twist GBs can be viewed simply as follows. Take a block of single crystal and cut into half. Twist one half of the block by about an axis perpendicular to the interface; the bring the two half blocks back into contact. The resul ting interface is a twist GB; see Figure 12 for an illustration. Twist GBs can also be viewed as a crystal with a number of screw dislocations aligned in an atomic plane. The illustration in Figure 12 is a relaxed twist GB between two simple cubic crystals. The plane of th e interface is (100) in this case. Tilt GBs can be constructed by the relative rotation of the two blocks forming the GB about a rotation axis along the plan e of the interface. If the atomic planes parallel to the GB on either side are the same, then the GB is symmetric. If the atomic planes are not the same, it is an asymmetric GB. Figure 13 is an example of an as ymmetric tilt GB. Just as twist GBs consist of screw dislocations, tilt GBs can be vi ewed as a series of edge disloc ations situated in the plane of GB. PAGE 19 19 A useful physical quantity associated with a GB is the GB energy. The GB energy can be defined as the energy difference between the singl e crystal and crystal with the GB per area: A E E EGB 0 (13) E ( ) is the energy of the structure with the GB, E0 is the energy of a single crystal with the same number of atoms, and A is the area of the GB. The GB ener gy is a function of the rotation angle about the rotation axis. Figure 14 is an example of the calculate d and measured energies of tilt GBs in aluminum.17 Some of the notable features are th e rapid increase and decrease with the rotation angle near 0 degree and 90 degree. Ther e are also cusps at ce rtain angles, corresponding to particular small unit cells (the values on each figure). Twist GBs generally have much shallower cusps than those of tilt GBs. Read and Shockley19 first proposed a model of EGB( ) for the low angle GBs based on the dislocation models. Later Wolf18 extended the model to the higher angles using phenomenological model. Th is latter model is called the extended ReadShockley model and will be used in Chapter 8 to analyze our results on heat transport in diamond GBs. One of the useful models to describe the GB st ructures is the coincident site lattice (CSL) model. The CSL model utilizes the fact that the two overlapping two dimensional lattices form a periodic pattern of the coincide nt lattice points. See Figure 15 for an example. The open and closed symbols in the figure are the atoms at simp le cubic lattice site in two adjacent planes. The squares are used to show the two dimensional unit cells in each lattice. Imagine that the lattice made of the filled symbol is above the plane consisting of the open symbols. When the top plane is rotated by certain angles relative to the bottom plane, some atoms in the top and bottom plane coincide. In Figure 15 the rotation angle is 36.5 and the plan e is (100). This interface is called the 5 (100) GB: denotes the fact the area of the tw o dimensional repeat pattern has an PAGE 20 20 area that is 5 times larger than the original unit ce ll in the simple cubic lat tice. This model is a very useful way to catalog different types of GBs. Tilt GBs can be cataloged in a similar manner. More detailed explanations ar e widely available in the literature.17, 18 Uranium Dioxide Uranium dioxide has been a technologically important material for more than half a century,11 particularly as a nuclear fuel in fissio n reactors for electrical power generation. Uranium is a rather common element and is generally found in the form of U3O8 ore Major producers of uranium ore are Australia, Canada, K azakhstan, the United States, South Africa, in order of decreasing amount of reserves.20 Perhaps the single most important reason that uranium is used for the nuclear fission react ors is the natural occurrence of 235U. 235U is the only naturally occurring stable fissile elements. After physic al and chemical proces sing of the ore, the UO2 form used in reactors, is obtained.11 The main reasons why uranium is used in the form of UO2 for fission reactor fuel are: i) its high melting point ~3000 K, ii) its chemical st ability with common cladding materials and coolant, iii) its fluorite structure, which can relatively easily accomm odate fission products.11 UO2 of course is not without disadvantages in nuclear fuel applica tions: i) low thermal conductivity (~5 W/mK at 800 K), ii) lower density of fissile elements than its metallic form, and iii) the relative difficulty of processing compared to its metallic form.11 Other phases of uranium oxides were also used for colorants for glasses and glazes until the 1940s.21 The crystalline structure of UO2 is the cubic fluorite structure (space group of m Fm 3) with the lattice parameter of 5.4698 at 300 K.2225 Figure 16 is the unit cell of UO2. Uranium atoms sit at the face centered cubic (FCC) sites (4a Wycoff positions), while oxygen atoms occupy the tetrahedral positions (8c sites) in the FCC U lattice. The ionic radii of U4+ and O2are r (U4+)=0.95 and r (O2)= 1.35 .26 Therefore the radius ratio, r (U4+)/ r (O2), is 0.704. The PAGE 21 21 structure is consistent with Paulings rule27 since the oxygen sublattice forms a simple cubic structure around the cations. The elec tronegativities of U and O are 1.7 and 3.5,1 thus the difference is 1.8. According to th e Paulings law of electronegativity28, this means that the bonding in UO2 is ~55 % ionic. Figure 17 is the binary phas e diagram of UO system obt ained from the thermodynamic calculation by Lewis et al., based on available experimental data.29 There is only a single phase in stoichiometric UO2 from room temperature to the melting point.2931 The UO2+ x phase is fairly wide especially at high te mperature. With a high concentration of oxygen in UO2+ x, the U4O9 phase can form, but the exact phase boundary is still under debate.29 By contrast, UO2x phase is very limited. UO2 can accommodate a number of point defects, especially when it is used in the nuclear reactor fuel. Intrinsic defects are VU, UI, VO, and OI with various charge states. U itself can take a wide variety of oxidation states: U2+ U3+, U4+, U5+, and U6+.32, 33 U2+, however, is very rare. The predominant intrinsic point defects in UO2 are the antiFrenkel pair.3, 3436 In addition to the intrinsic defects, a large number of extrinsic defects exists in th e nuclear fuel, and many of them are produced during th e nuclear reactions.11 Figure 18 gives the yields of elements produced as extrinsic defects during the nuclear fission react ions in a typical pressurized water reactor (PWR).37 The effects of these extrinsic defects on properties depend on the type and on their concentrations. The study of the effects of fission products themselves is a large area of research and is beyond the scope of this work. Point defects in UO2 are known to affect its thermophys ical properties. For instance, Ronchi and Hyland35 analyzed the available heat capacity data of UO2. They found that below 1000 K the heat capacity is determined by the harmonic lattice vibrations and that the Debye PAGE 22 22 model applies very well. Between 1000 and 1500 K, the anharmonic interactions play an important role and as a result the heat capacity co ntinues to increase with increasing temperature. From 1500 to 2760 K, defect formation (mainly Fre nkel defects or is it antiFrenkel?) increase the rate of rise in the heat capacity. Finally above 2760 K, othe r point defects such as Schottky defects also become important and the heat capac ity increases further. Thus different types of defects operate in different range of temperat ure and influence the thermodynamic property. In addition to the point defects in UO2, electronic defects can be important. Electrons in dielectric materials like UO2, in which the ions are highly polarizable, can propagate when strongly coupled with lattice vibrations.38 The quasiparticles asso ciated with this phenomenon are called polarons. In particular, socalled small polarons, are common in UO2 and are involved in the charge transfer at high temperature. In the small polaron model, electrons hop from ion to ion, unlike the elec trons in the conduction band. Sma ll polarons also contribute to the thermal conductivity of UO2 at high temperature. However, the contribution is rather small (at 2000 K, the thermal conductivity by polaron is the 0.31 W/mK while by phonon is about 1.87 W/mK)39 and will not be considered further. GBs in UO2 also play very important roles in nuc lear fuel. As mentioned earlier, GBs contribute to the thermal resistance of crystallin e solids. However, there has been no systematic study of the grain size effect on the thermal transport of UO2. The reason may be that the grain size in UO2 fuel at beginning of life is fa irly large (size in a typical UO2 fuel pellet is of the order of 10 m)3, 40 and the interfacial resistance is not a si gnificant factor in thermal transport. However, it is know that the grain subdivision occurs in high burnup stru ctures, and the grain size can decrease to ~150 nm on average.40, 41 At this grain size, the interfacial thermal resistance may be significantly enhanced. Experi mentally it is difficulty to determine the effect PAGE 23 23 of grain size because the high burnup structures also contain a la rge number of point defects and pores which also decrease the thermal conductivity.40, 41 Other than the thermal transport, GBs also ha ve a number of other effects. One such example is the segregation of fission products. Pores formed during the usage in the nuclear reactor are generally found at the GBs.42, 43 Pores smaller than 2 m are removed in fine grain (<6 micron) sample, but pores are stable if the grain size is large (>20 m).4446 The amount of fission gas release and the rate of swelli ng can be reduced by increasing grain size.47 Mechanical properties such as hardness47, elastic44, and creep48 properties are also affected by the presence and the size of the grains in the sample. Diamond Diamond is an allotrope of carbon and has a number of unique prope rties. It is the hardest of known materials (hardness of roughly ~137167 GPa).49 It has the highest thermal conductivity of ~2200 W/mK at 300 K am ong naturally occurring materials6. The refractive index of diamond is 2.417 making it id eal for some optical applications7, and the brilliant refraction is appreciated for a gem stone. In diamond, carbon atoms are covalently bonded in the diamond structure with the space group of m Fd 3. Figure 19 is the unit cell of the diam ond structure. There are 8 carbon atoms in a unit cell at 8a sites, and each atom ha s 4 nearest neighbors. The lattice parameter determined by Xray diffraction is 3.5673 .50 Naturally occurring diamond contains less than 1 % nitrogen impurities and possibly much lower concentrations of hydrogen and boron. Thes e diamonds are categorized into typeIa, Ib, IIa, IIb and other rare types, based on the concen tration and the types of the impurities. Type IIa diamond is often used in t echnological applications,6 for which most important and the naturally occurring impurity is nitrogen. Nitrogen atoms resi de in diamond as substitutional defects, strain PAGE 24 24 the lattice and also modify the el ectronic structure. In addition to the impurities, natural diamond contains isotopic defects; in particular 13C appear at concentration of ~1.1 %.6, 51 Synthetic methods such as chemical vapor deposition (CVD) are often used to prepare diamonds of high purity.52, 53 Beside the single crysta ls, polycrystalline diamond is technologically important because of its unique physical properties. In pa rticular, ultrananocrystal line diamond (UNCD) is a very promising multifunctional material.54 Objective and Scope of Study The objective of this thesis is to elucidate the mechanisms by which atomic level defect structures influence therma l transport properties in UO2 and diamond. Atomic level simulation techniques are used to accomplish this tas k. A methodology for the simulation of radiation damage in UO2 is also developed. PAGE 25 25 Figure 11. Point defects in cr ystalline solid. Black and grey symbols represent atoms of different types; repres ents a vacant site. Vacancy Interstitial Schottky Defect Frenkel pair PAGE 26 26 Figure 12. Pure relaxed single tw ist grain boundary (GB) between tw o cubic lattices. The GB is in the plane of the page. Region A and B are the region of coherent and incoherent interfaces, respectively. [Adapted from the reference 17 (Page 317, Figure 13.10).] PAGE 27 27 Figure 13. Asymmetric tilt GB between two cubi c lattices. The GB is shown by the dashed line. Symbol is used to indicate the core of th e edge dislocation. [Adapted from the reference 17 (Page 316, Figure 13.9).] PAGE 28 28 Figure 14. The GB energy as a function of ro tation angle in <100> (a,b) and <110> (c,d) symmetric tilt GBs of aluminum. (a,c) are calculated values and (b,d) are measured data. [Adapted from reference 17 (Page 315, Figure 13.8).] PAGE 29 29 Figure 15. Lattice structure in 5 (001) =36.5 GB. Dotted circles ( ) and open circles ( ) represent atoms in two adjacent layers. Black dot ( ) indicates the coin cident sites of the two atomic layers. Black solid and dash ed lines are used to indicate the cubic unit cell in each crystal. The dash dot line shows the unit cell of 5 GB. [Adapted from reference 17 (Page 334, Figure 13.19 ) with minor modification.] PAGE 30 30 Figure 16. Fluorite structure of UO2 unit cell. Blue and red denote uranium and oxygen atoms, respectively. PAGE 31 31 Figure 17. Phase diagram of UO system. S and L in subscripts are used to indicate the solid and liquid phase. The numbers indicted in each phase are used in their reference for models of partial pressure of oxygen. [Adapted from reference 29 (Figure 4).] PAGE 32 32 Figure 18. Fission product yield as a function of atomic number in a pressurized water reactor (PWR) after 2.9 % burnup. [Adapted from reference 37 (Figure 1).] PAGE 33 33 Figure 19. Unit cell structure of diamond. Th e bars are used to guide for the eyes. PAGE 34 34 CHAPTER 2 INTRODUCTION TO THERMAL TRANSPORT IN MATERIALS Basic Concepts of Thermal Transport Heat transfer is a funda mental mode of energy transport in na ture. It is so fundamental that it exists everywhere from the interactions of mo lecules to the radiation of the cosmic microwave background of the universe. In todays worl d of technology, understand ing the mechanism of heat transfer is more important than ever. One example of the application of thermal ener gy management is in electronic devices. The size of logic circuits continue s to shrink, and the leng th scale of a transi stor is reaching ~45 nm today.55, 56 Accordingly, the energy de nsity of the device is also increasing. For example, chip level power densities as of Fall 2006 we re estimated to be the order of 100 W/cm2.57 When the temperature of the transistor rises, the performance of the chip degrades: mobility and threshold voltage lowers and the leakage current increases.58, 59 Therefore, such high power density electronic devices require be tter heat removal technologies. Another example is in the gene ration of electrical power. Both scientific and technological advances today require an increasing supply of energy. In addition, the need for minimizing the dependence on foreign sources for energy production is naturally diversifying the field of power generation technologies. As a result nuclear en ergy is again receiving considerable attention. Nuclear power technology has the advantage of high power density per weight of the fuel, and no emission of the greenhouse gases. However, the concerns for the safety of radioactive waste have been a major issue. Therefore the conti nuous improvement of the technologies involved in the nuclear materials is imperative. Both fuel and structural materials used in nuclear power generation are under extreme conditions of radiation and temperature. For example, as discussed in Chapter 1, UO2 is the standard nuclear fuel material today. A typical ~1cm diameter fuel pin PAGE 35 35 in a PWR under operating conditions experiences a temperature of ~800 K at the surface and nearly 1600 K. Ideally, all the he at generated by the nuclear fissi on reactions is quickly released to the coolant so that the thermal energy can be converted to create the steam which drives the turbine and generates electricity. However, beca use of the high temperature of the material, the fuel pins experience physical and chemical ch anges in their propertie s that often result in undesirable effects such as deformation and cr acking of fuel pellet, and possible chemical interaction with coolant.11 Therefore, understanding the heat transfer mechanisms in nuclear fuel materials is particularly important. In this chapter a general discussion on the type s of heat transfer and mechanisms in solids is given, with particular emphasis on phononmediated thermal transport. Three Modes of Heat Transfer Co nduction, Radiation, Convection Heat transfer in materials can occur th rough three distinct mechanisms: conduction, radiation and convection. Conducti on occurs by the motion of the energy carriers such as electrons, lattice vibrations; ra diation takes place via electromagnetic waves. By contrast, convection involves the mass transport of the constitu ent atoms. In solid materials, convection is therefore not relevant and will not be discussed further. In electronic conductors, electrons are the ma jor thermal energy carrier s. The electronic contribution to thermal transport depends on the el ectrical conductivity and temperature. This is captured in the relation in WiedemannFranz law; LT (21) and are the thermal and elec trical conductivities, and T is the temperature. L on the right hand side is a constant calle d the Lorenz number, (3/2)( kB/ e )2=1.11108 W /K2, where kB is the PAGE 36 36 Boltzmann constant and e is the electron charge. This la w holds to a good approximation if the collisions of the electrons are elastic.60 For intrinsic semiconductors and insulators lattice vibrations (phonons) are the major carrier of thermal energy. This mechanism will be discussed in detail. It should be noted that the conduction through electrons and phonons can c ouple and operate simultaneously; they are not, however, relevant for the materials of interest here. Whatever the carrier may be, once the system reaches steady state and the temperature distribution across the system is known, one can apply Fouriers law: T J (22)J is the heat current, is the thermal conductiv ity of the medium, and T is the temperature. Notice that J and T are vectors; therefore, is a secondrank tensor. The negative sign indicates the fact that the heat flows from higher to lower temperature regions. The assumption here is that the medium is a continuum such th at temperature can be de fined every point within the volume of interest, and that the cha nge of temperature is linear in space. Radiation is another form of heat transfer involving the transmission of electromagnetic waves. Every material emits radiation at finite temperature because of the thermal oscillation of ions and electrons.61 The oscillating charge d particles emit electromagnetic waves of the corresponding energy.62 The spectrum of the radiation follows the Planck distribution of a black body radiation.61 The total emissive power is63, 4T Js s (23) s is the StefanBoltzmann constant, 5.670108 W/m2K4, and s is the emissivity of the material. This mechanism is important at high temperature because of the T4 dependence. A good example is the thermal transport in thermal barrier coating. Yttria stabi lized zirconia (YSZ) is a PAGE 37 37 commonly used for thermal barrier co ating of gas turb ine engine blade 64, 65. At the surface of YSZ layer, the temperature is ~2000 K or so, a nd at the inner interface w ith the bond coating it is above 1300 K depending on the way the coating is applied. Although heat transfer in YSZ is dominated by the conduction mechanism, the radi ative component cannot be ignored at these high temperatures. For example, for a single la yered 8 wt.% YSZ the radiation heat current density can be almost 25 % of the total heat flux for YSZ between the surface temperature of ~1700 K and bond coating temperature of 1400 K with 250 micron thick coating.66 (The previous sentence is too complicated please simplify) In the case of UO2, Hyland67 estimated the radiation contribution to the thermal conductiv ity at very high temper ature (>2900 K). His estimate suggests that the radi ation component in the thermal conductivity is much less than other contributions. Bakker et al .68 and Hayes et al.69 also made estimates of the contribution of radiation to heat transfer in porous UO2 using finite element methods. Their calculations gave ~1 % of the total thermal conductivity at 3000 K is due to radiation. Because of this reason, the radiative thermal transport in UO2 is normally excluded from the analysis 3. In many situations, heat is transferred not only by a single mechanism but by more than one simultaneously. Effects of the individual m odes must be taken into account to understand the entire processes. Phonon Heat Transfer in Electronic Insulators As described above, in electronic insulators, phonons are the major ca rriers of the thermal energy. A model of phonon mediated heat transfer can be derived from the scattering of a phonon gas. The following discussion is based on standard discussions in solid state physics texts.38, 60, 70 To describe phononmediated heat transfer, the thermal conductivity of a solid can be considered as follows. Imagine phonons as partic les behaving like gas molecules. Assume that PAGE 38 38 there is a container having a num ber of phonons with concentration, n Suppose the system is in steady state between a temperature drop, T If c is the specific heat per phonon, a phonon will lose energy c T by moving from the hot end to the cold end is, c T If the local temperature gradient is dT/dx and l is the mean free path of the particles, then the temperature drop over distance l in direction i is T =( dT/dxi) li (The Einstein summation convention is assumed on repeated indices from here on unless explicitly stated otherwise.) Then, the net flux of the energy carried by the phonons is, j j i V j j i ix T v v c l v T c nv J (24) where vi is the mean velocity of the particles, li is the mean free path, cv is the volumetric specific heat, and is the average time between two collis ion evens of phonons, and is known as the relaxation time. By Fouriers law, th e thermal conductivity is the given by, j i V ijl v c (25) If the system is isotropic, th e mean velocity of the phonons is 3 v vi Therefore, the thermal conductivity is given by, l v cV3 1, (26) As discussed above, phonons, i.e. lattice vibratio ns, can be the major carrier of thermal energy. Despite the simplicity of the expression, th is kinetic formula is very useful. Indeed, it will be the foundation for our phenomenological anal yses of phononmediated thermal transport. To understand the mechanisms of the thermal resistance to heat fl ow, two types of phonon scattering processes must be introduced. Phonons pr opagating in a crystall ine lattice have wave vectors which are constrained by the periodic nature of the crystallin e structure. A structure that is crystalline in a real space exhibits periodicity in the reciprocal (or momentum) space. This PAGE 39 39 periodicity in the momentum space is given by the reciprocal lattice vector, K For example, an FCC structure in real space has a BCC structure in re ciprocal space, and vice ve rsa. Just as a real space lattice structure has a primitive cell gi ven by the WignerSeitz cell construction, a reciprocal space has a corresponding primitive cell given by the first Brillouin zone (BZ). The reciprocal lattice structure can be reproduced by replicating the first BZ and applying translational symmetry. The physical importance of th e lattice structure in th e reciprocal space is that it places constraints on the momentum conservation of phonons during the scattering processes. Consider a th reephonon scattering process: K p p p 3 2 1. (27) This can be seen as two phonons with momenta,1p and 2p colliding to form a single phonon with momentum3p. Because of this conservation law of the crystal momentum, there are two possible processes: normal (N) and umklapp (U) pr ocesses. An Nprocess is a phonon scattering event in which the total phonon mo mentum is conserved, i.e., K =0; this is analogous to the scattering processes of particles such as gas molecules. A Uprocess, however, is a scattering process in which total phonon momentum chan ges by a reciprocal lattice vector, i.e., K 0. Since the 1st BZ is defined by the recipr ocal lattice vector, the sum of all the phonon momenta for a Uprocess must be outside of the first BZ. Uprocesse s play distinct roles in the thermal transport. In Nprocesses, the total crys tal momentum is conserved. This means that, if phonons are created on one end, the same number of phonons must disappear on the other end to conserve the net flow. Since phononphonon scattering cons erves phonon energy, there cannot be any thermal resistance. Therefore, Nprocesses cannot contribute to the thermal resistance. On the other hand, in Uprocesses, the crystal momentum is not conserved, so that there can be a net heat flow and energy can dissipate. PAGE 40 40 As the thermal conductivity expression, Equation 26 include the specific heat, it is helpful to introduce the commonly used approximati on to gain insight to its behavior. The internal energy from the vibration of the atoms for e ach polarization is given by, 0 n D d U, (28) where D is the density of st ate (DOS) of phonons and n is the ensemble averaged number density of phonons, which follo ws BoseEinstein statistics: 1 exp 1 T k nB (29) The internal energy can be used to calculat e a number of thermodynami c quantities through the Maxwells relations. However, this is not a st raightforward task; the main complication being that D is a rather complicated function of and is singular. In orde r to gain physical insight through simpler analytic expressi on, there are two commonly used approximations: the Einstein approximation and Debye approximation. The Einstein model assumes ( k )= 0; that is the atoms in the lattice all vibrate independently of each other with the same frequency. Thus D ( )= N ( 0), where ( ) is the Dirac delta function. Because of the simple expression of D ( ), Equation 28 can be directly integrated. From the simple thermodynamics relation cv=( U / T)V, we get, 2 21 3 x x B Ve e x Nk c. (210) Here we introduced dummy variable x = / kBT It is well known that the Einstein model underestimate the specific heat at low temperatures ( T << D).38, 60 (Shouldnt there be an integral here the Einstein model gives cv vs T as the final result ). At high temperature, ( T >> D), it predict the DulongPetit law, as it must to be consistent wi th classical thermodynamics. PAGE 41 41 The Debye approximation assumes =vk where v is the speed of phonon assumed to be independent of k Under the Debye approximation, th e density of state is simply, D Dv V D 0 23 2 2, (211) where V is the volume of the solid. The maximum frequency, D of phonons is limited by the size of the primitive cell, which is determined by the number density of atoms, n : n vD 3 22 (212) Now the internal energy can be written as, Dx x D Be x dx T T k U0 3 31 9, (213) where, n v k kB B D D3 22 (214) is the Debye temperature. Finally the vol ume specific heat can be given by the, Dx x x D B Ve e x dx T Nk c0 2 4 31 9. (215) The integral is called the Debye integral of 4th order, and in general it must be computed numerically. However, it is instructive to consider the limiting cases of low and high temperature relative to the Debye temperature: D B D D B VT Nk T T Nk c 3 5 123 4. (216) The high temperature value is the same as the DulongPetit law. Although both of these models ar e rather oversimplified, they are useful approximations that can be used to gain insight of the act ual behavior of phonons of different modes. PAGE 42 42 Using the above results for the specific he at, the temperature dependence of thermal conductivity in solid can be described based on the phonon model. At high temperature ( D < PAGE 43 43 ~ T3 under the Debye approximation. As a result, th e thermal conductivity thus also has the same T3 dependence. This behavior agrees w ith a number of experimental observations 60. In the intermediate range of temperature ( T ~ D), the situation is rather complicated and we have to resort to extending the behavior from the high and low temperature limits. When Uprocesses are the only source of the thermal resistivity, rises exponentially with decreasing T ~exp( D/ T ), because of the quantization of the energy levels 70. Since the specific heat rises as T3 at low temperature, is commonly expressed as, bT T AD n D exp. (217) A, b, n are to be determined by fit. This function has the maximum at T = D/ nb where the value of nb is typically about 1535. (Table 21.) For summary, the temperature dependence the heat capacity of a pure crystalline solid is given by the Debye approximation. The temperat ure dependence of thermal conductivity varies from ~ T3 to Tnexp( D/ bT ), then ~1/ T Figure 22 illustrates this temperature dependence using data on Ge. Point Imperfections A crystalline solid can contain a variety of point defects such as interstitials, vacancies, isotopic defects, and substitutional defects. These defects can differ from the ions in solid with regards to their mass, size, and interactions with other atoms. In order to understand the contribution to the thermal resistance from point de fect scattering, we must return to the original expression of the thermal conductivity in Equation 26: l v cV3 1, PAGE 44 44 Here, it is useful to think the conductivity as a sum of the c ontribution from each mode of wave vector k j j j j V j j j j Vk d k k v k c k d k l k v k c ) ( ) ( ) ( 3 1 ) ( ) ( ) ( 3 12 (218) where ) ( kj is the relaxation time. The main pr oblem one faces in calculation of is to calculate the relaxation time. This may be done by so lving the Boltzmann transport equation. Exact solutions are generally unattainable due to ma thematical complexity, and approximations must be used. Klemens74 applied secondorder perturbation theo ry to the scatteri ng of phonons by point imperfections. The point defect in this analysis is assumed to be a generic atomic defect with mass Mi, force constant fi, and radius Ri such that the difference from the corresponding values of crystalline host material are given by, Mi= MiM0, fi= fif0 and Ri= RiR0, where M0, f0 and R0, are the mean values. The derivation of the formal expression is rather le ngthy and the details are not central for this study. The basic procedur e is following. Write the transition rate of phonon level by the phonondefect scatteri ng using the Fermis golden rule75. Then give Hamiltonian in terms of the defect properties such as the differen ces in masses, radii, or force constants. Then the rate of the number of phonons created by the scattering ca n be obtained by integration the transition rate in momentum space. Finally the scattering rate can be wr itten in terms of the change in the number of phonons. In the end, the scattering rate by the static poin t imperfections in the long wave length limit at low temperature ( T << D) is given by, i i i i i j j jR R f f M M x k v k k2 2 3 43 2 6 6 12 1 ) ( ) ( ~ ) ( 1 (219) PAGE 45 45 The overall numerical constants ar e suppressed since they are not re levant to the discussion. The summation is over the atomic species. xi is the fraction of each defect species. Notice the mass, force constant, and radius difference comes into th e scattering rate. Another important point is the dependence of th e scattering rate on 4. Since = vk and the original formal expression for the mean free path is proportional to the ratio of 4th order and 8th order Debye integrals, the mean free path is proportional to T4 at T << D. This in turn gives ~1/ T dependence since the specific heat contribution gives T3 at T << D. Equation 219 also suggests that the thermal resistance is linear in the concentration of de fects, but this dependence is ra ther artificial and introduced purely for the sake of mathematical convenien ce. A similar analysis shows that, at high temperature, the thermal resistance from point imperfections is independent of temperature70, 76. GBs and Surfaces GBs and surfaces are also sources of thermal resistance for phonon mediated heat transfer. In order to quantify the thermal transport property of an interface, a phenomenological relationship analogous to Fourie rs law is generally used:77, 78 T G JK (220) GK is the interfacial thermal conductance, also known as Kapitza conductance. T is the temperature drop at the interface. There are two widely used phenomenological models for studying the interfacial resistance. The first is the acoustic mismatch model (AMM). In this model, the crystal is treated as a continuum while the interfacial resistance arises due to the acoustic impedance mismatch between the two solids on both side of the interface79. A key weakness of this approach is that this model does not take into the discrete lattice structure of ions consisting of solid. Thus, this approach cannot be justified if the phonon wavelength is comparable to the lattice parameter of the crysta l. Furthermore, this model does not convey any PAGE 46 46 information of the effect of the interfacial structure on the therma l transport. Because of these limitations, AMM predicts zero interfacial resist ance in GB structures, a result that is in qualitative disagreement with experiment. Th e second model is the diffuse mismatch model (DMM). DMM assumes independent scattering of phonons with the interf ace. In addition, all phonons are assumed to be randomly scattered at th e interface such that they do not have any memory of their origin79, 80. Thus the fraction of energy transf erred across the interface is given by the ratio of the phonon density of states (DOS) on either side of the GB. The DMM has been shown to work in some cases80, but predicts that the phonon transm ission coefficient for all twist grain boundaries should be 50%, a result th at disagrees with simulation results.77 There is also an analytical approach by Klemens for GB scattering76. However, this approach is severely limited by the approximati ons which are needed to make the calculation tractable. Klemenss model assumes that the GB consists of multiple edge dislocations. He has shown that the for low angle scattering of l ong wavelength phonon at low temperature, the scattering rate is proportional to the square of the phonon frequency, 1/ ~ 2. He argued that the thermal conductivity of a solid with GB s and surfaces should increase with T3 below D following the specific heat but should be constant at higher temperatures76. Based on the Klemens work, Stoner and Maris gave the maxi mum limit to the Kapitza conductance in the form of the Debye integral81. However, their experimental data on heterogeneous interfaces were in poor agreement with their calculation. Because of the lack of a satis factory theoretical model, exte nsive work in the theory of interfacial thermal transport is still needed. Thermal Transport in Disordered Solids In crystalline solids, the thermal transport is limit ed either by phononphonon scattering or phonondefect (point defects, dislocations, GBs, su rfaces) scattering. In highly disordered solids, PAGE 47 47 including amorphous materials, th e thermal transport mechanism is fundamentally different from that in crystalline solids. In noncrystalline, the atomic arrangem ent is so disordered that the mean free paths of phonons are extremely s hort and are independent of temperature82, 83. However, the specific heat is almost lin ear in temperature at low temperature38. Above the Debye temperature, cV increases slowly with T and eventually becomes almost temperature independent. Therefore, assuming the form of Equation (26), the thermal conductivity of disordered solid follows that of its specific heat. Slack83 calculated the ther mal conductivity of the Debye solid, assuming that the mean free path of the phonon is half of its wavelength for each mode. Under this assumption, he defined the minimum thermal conductivity of glasses. Cahill and Pohl84 made the same assumption and gave the minimum thermal conductivity of amorphous solids in terms of a Debye integr al. Although these models give reasonable agreement for the thermal conductivity of highly disordered solids, the physical basis of describing thermal transport in an amorphous solid using phonon gas model is questionable because the phonon mean free path in these materials is so short (on the or der of the interatomic distance) that the concept of phonons itself seem inappropriate. Allen et al.85, 86 took a different approach to this problem, and developed a model for the thermal conductivity of amorphous materials at hi gh temperatures. Their model is based on a harmonic approximation analysis of the vibrationa l modes of the system. They demonstrated that the majority of the thermal energy is transferred through di ffusive modes termed diffusons. In their simulations of amorphous Si, approximately 93 % of the vibrational modes correspond to the diffuson modes. Other mode s in the higher frequency bands are called locons due to the fact that they are truly sp atially localized. At th e other extreme of the phonon spectrum, the low frequency modes are si milar to phonons, having well defined wave PAGE 48 48 vectors and polarization vectors. They are called propagons and constitute only about ~3 % of the entire vibrational modes for aSi. Although th e propagons are very effi cient in transfer of thermal energy, their contribution to thermal conductivity is small because of the small population. Within this approach, the temper ature dependence of the thermal conductivity mostly comes from the temperature dependence of specific heat. This in turn agrees with the temperature independent behavi or at high temperature. Th e approach developed by Allen et al. is useful not only in analyzing amorphous solids, but in analyzing the thermal conductivity in heavily doped crystalline materials;87 we will use this framework to understand the thermal conductivity of UO2 x. PAGE 49 49 Table 21. Values of the Debye temperature and the maximum in the experimental thermal conductivities of selected materials. D (K) Tmax (K) D/ Tmax NaF 49260 2038 24.6 LiF 73060 2060 36.5 C (diamond) 22006 702 31.4 Si 62560 4088 15.6 Ge 36060 1589 24.0 MgO 94090 302 63 31.1 Al2O3 90091 3070 30.0 Tmax is an estimate from the plot of ( T ). PAGE 50 50 Figure 21. Temperature dependence of the thermal conductivity of selected materials.92 PAGE 51 51 Figure 22. Temperature dependence of the thermal conductivity of Ge.38 Temperature dependence given by the Debye approximation is shown in the figure. The two data sets are used to compare the effect of is otopic defects. The enriched Ge has 96 % 74Ge. Normal Ge has 37 % 74Ge and the rest consists of 72Ge (27 %), 70Ge (20 %), 73Ge(8 %), and 76Ge(8 %). T / 1 ~ bT TD nexp ~ 3 ~ T PAGE 52 52 CHAPTER 3 SIMULATION METHODS Introduction There are a number of ways to simulate mate rials in condensed phases. There are several types of simulations which invol ve electronic structure calculati ons such as density functional theory (DFT) and HartreeFock method.93 Molecular dynamics (MD) and MonteCarlo (MC) treat atoms as the fundamental entities. Finiteelement94 and phase field95 methods are continuum approaches that can handle >m in size. The simulation method must be chosen based on the process of interest. For example, if one is interested in und erstanding the electronic structure of materials, DFT or HF may be appropri ate. If equilibrium behavior is required, MC is a common choice. MD is particularly useful in the study of nonequilibrium dynamic processes at the atomic level. MD treats atoms as point particles with no internal st ructure. For each atom, mass and charge are assigned. All atoms are treated clas sically and their dynamics are obtained by solving the Newtons equation of motion. Details of how to integrate the equation of motion are described in the following section. Since the elec tronic degrees of freedom of atoms are ignored, MD is capable of handling a larg e number of atoms and molecules. With todays computational power, MD can handle approximately 107 atoms on parallelized sma ll workstation computers. With a super computer, simulation of 109 atoms is possible. The onl y physical input required by the MD simulation is the interato mic potential. This has a signifi cant implication. It means that if the quality of intera tomic potential is good, the cal culated properties of the material as a whole should be reproduced to a reasonable extent. Un fortunately, while a given interatomic potential may describe well some aspects of the physical be havior of a material, no interatomic potential has been found to capture all the properties of a material in realistic manner. PAGE 53 53 Numerical Integration As mentioned in above, MD simulation requi res integration of Newtons equation of motion, which is a second order ordinary differen tial equation. There are a variety of methods available to integrate such equation.96, 97 Here, we focus only on the Gear predictorcorrector method since it is the algorithm used in our inhouse MD code. The 5th order Gear predictorcorrector method96 uses a high order Taylor expansion of atomic position about the current position, x ( t ), up to the 5th derivative. Using the current velocity, v ( t ), acceleration, a ( t ), and higher order derivatives of x ( t ), the next position of the atoms are predicted: 5 5 5 4 4 4 3 3 3 2) ( 5 1 ) ( 4 1 ) ( 3 1 ) ( 2 1 ) ( ) ( ) ( dt t dt x d dt t dt x d dt t dt x d dt t a dt t v t x dt t xp (31) Once the next positions of the atoms are pr edicted, the force on each atom is calculated using the predicted positions. With the new predicted force, the atomic positions are corrected by the corrector step: )} ( ) ( { ) ( ) () 2 ( 0 ) 2 ( ) ( ) (dt t x dt t x c dt t x dt t xp i i p i c (32) The superscript ( i ) is used to indicate the order of derivative of the trajectory. Thus, ) ( 1 ) () (t dt x d n t xn p n i p (33) is implied. ci is an array of numerical constants: c0=3/20, c1=251/360, c2=1, c3=11/18, c4=1/6, and c5=1/60.96 ) () 2 ( 0dt t x is obtained by the plugging ) () 2 (dt t xp into the equation of motion. ) ( 1 ) () 2 ( ) 2 (dt t x F m dt t xp o (34) Periodic Boundary Condition Although MD simulation technique can handle a large number of atoms, the number of atoms is quite small compared to that in real physical systems, which contain of the order of 1024 PAGE 54 54 atoms Therefore in order to simulate a la rge bulk material, periodi c boundary condition (PBCs) are used. Figure 32 illustrates PBCs in 2 dimens ions. The actual simulation cell is the central shaded cell of size L L Suppose, for simplicity, that th e simulation cell contains only 5 molecules as shown in the figure. When mol ecule 1 crosses the boundary of the simulation box, its periodic image reenters the cel l from the other side. Clearly th is is not what actually happens in a real materials; however but if there are a large number of atoms in the simulation cell, the environment any atom experience will resemble to that of an infinitely large system. Another important effect of the PBCs is in th e collective vibrations of atoms and molecules and pressure waves. If the PBC is not used, the simulation cell boundaries are terminated by the surfaces and only standing vibrational waves are allowed. Having PBCs actually allow traveling waves in the system allowing the material to be have like a large medium. It should be noted, however, that even if traveling waves are a llowed in the simulation cell, their maximum wavelength is limited by the size of the cell, with longest wave length being twice the length of the simulation cell size. As we shall see, thes e constraints are very important for thermal conductivity calculations. Interatomic Potentials The interatomic potential is perhaps the mo st important piece in the MD simulation technique. As mentioned in the introduction to this chapter, the interatomic potentials are the only model input for MD simulations, and all the phys ical properties are determined by it. It is quite intriguing that even the rather crude models of interatomic interactions typically used can reproduce physical properties w ith reasonable accuracies. Longrange Interaction Evaluation of the longrange 1/ r electrostatic interactions in ionic systems is one of the trickiest problems in MD simulations. We wi sh to compute the electrostatic energy by, PAGE 55 55 j i j i ESr q q E04 1, (35) where qi and qj are the charges of ions i and j separated by a distance r, and 0 is the electric permittivity of vacuum. The crux of the problem is that this sum is conditionally convergent. There have been several approaches proposed to get around this problem.96, 98 In following, two methods used in hour inhouse MD code are de scribed. The discussion focuses mainly on the idea behind the methodology rather than the details of derivation. Efficiency of the methods will also be discussed for the later purpose Ewald sum The idea of Ewald sum is to rewrite the su mmation in Equation (32) in terms of the electric charge density and expr ess these densities by a sum of de lta functions. If we assume every point charge qi is surrounded by a diffused charge of opposite sign, qi cancels exactly (Figure 33). The potentia l energy due to charge qi is the fraction of charge qi that is not cancelled, and it goes to 0 rapidly at large distances. Using th is approach, the conditionally convergent sum may be determined as: ij N j i ij j i N i i k ESr r q q q k k k V E erfc 2 1 4 exp 4 2 11 2 0 2 2 (36) The first term is the potential en ergy due the charge distribution ( r ) given in reciprocal space. The second term is the co rrection for selfinteract ion. The third term is the electrostatic energy of the point charges screened by oppositely charged Gaussians given in real space. Equation (36) is clearly compli cated and the implementation of this in an MD code is not straightforward. Furthermore, for the applicat ion of this method to non3D periodic system, further modifications are required. PAGE 56 56 In terms of efficiency, the time required for th e evaluation of electrostatic interaction by Ewald summation scales as N3/2 where N is the number of atoms in the system. This is better than a straight forward evaluation of N2 scaling (Equation (35)); th ere are alternative methods which scale much better with the system size.96, 98, 99 Direct summation The 1/ r truncated summation method developed by Wolf et al. is a powerful alternative method to the Ewald summation.99 From here on this method is re ferred to as di rect summation. The idea is based on the two key observations: 1. The electrostatic energy by the double sum of Equation 35 depends strongly on the environment of the atoms. 2. The electrostatic energy of an ionic cr ystal (Madelung energy) is shortranged. Consider the case of NaCl crystal, for example. Imagine a sphere of radius Rc in the crystal. If one calculates the electrostatic energy using Equation 35, the values depends significantly on the choice of Rc. (Figure 34) Assuming that the net charge of the ions in the sphere is q there is a clear linear depende nce in the electrostatic energy on q The Madelung energy of the crystal corresponds to the electrostatic ener gy for a charge neutral system. Thus, if it is possible to find the volume of the sphere such that q =0, a reasonably good value of the Madelung energy can be obtained. This can al ways be accomplished by smoothly compensating each charge in the sphere of Rc by equal and opposite charges at the surface of the volume. However, this itself is not sufficient for the di rect application because of the slow oscillatory decay of the potential as a function of Rc. Additional damping of the pair potential is required for the method to be practically useful. With all these in the consideration the electr ostatic energy of the i onic crystal of N atoms can be written as, PAGE 57 57 N i i c c N i R r i j ij ij j i R ij ij j i c ESq R R r r q q r r q q R Ec ij c ij1 2 1) ( erfc ) ( erfc lim ) ( erfc 2 1 ) ( (37) where is the damping constant, whose value is chosen to produce the most efficient, yet accurate, results.. This method is also shown to work with highly disordered condensed phase such as liquid or amorphous solid. The major advantage of this method is the computational efficiency. Because of the spherical truncation and fast convergence, the co mputation of the energy and forces is order N This high efficiency is necessary for simulations of very large systems. The direction summation method is used in all of the simula tions described in this thesis. Shortrange Interactions The shortrange interaction is determined from the type of material. Two commonly used shortrange interatomic potentials for ioni c and covalent solid s are the Buckingham100 and Tersoff101, 102 potentials. They are used he re for the description of UO2 and diamond respectively. The simplest form of the Buckingham potential is given by, 6) / exp( ) (ij ij ij ij ij ijr c r A r V (38) Subscript i and j are used to indicate the type of interaction between atom i and j Aij, ij, cij are fitting parameters. The first term represents the repulsive potential from the overlap of the electron clouds of the two atoms. The second te rm is the Van der Waals attraction. There are a large number of the set of parameters for i onic solids; see the work of Lewis and Catlow103, 104, for instance. We will use a variant of this potential in Chapters 47. PAGE 58 58 The Tersoff potential was develope d to describe covalent materials, taking into account the fact that the interatomic bonds depend on the environment of the atom s. The total energy E of a covalent solid is given by, ) exp( ) ( ) ( cos exp ) ( ) ( ) / exp( ) exp( ) exp( ) ( ) ( ) ( 2 12 1 0 2 1 ij ij c ij j i k ijk n ij ik ij ij ij ij ij ij ij ij c ij j i ijr r f r w d c r w r w z b z B B r B r A r f r V r V E (39) Here, Aij, B0, b c d n 1, and 2 are parameters fitted to some bulk properties of the solid. ijk is the angle between the bond ij and jk fc( rij) is given by, D R r D R r D R D R r D R r r fij ij c, 0 ) ( 2 sin 2 1 2 1 1 ) (. (310) R and D are the fitting parameters. It is not easy to give physical meaning to all these terms except that Aij and Bij terms in V ( rij) are repulsive and attractive te rms. The potential parameters for carbon are, Aij=1393.6 eV, B0=346.74 eV, b =, c =38049, d =4.3484, 1= 3.4879 n =0.72751, 2=2.2119 R =1.95 and D =0.15 .101 This potential is used to simulate diamond in Chapters 8 and 9. Thermostat In order to simulate materials properties under realistic condition, it is necessary to control the temperature of the material. There are a several commonly used techniques for controlling the temperature. The simplest of all is the velo city rescaling thermostat. The idea of this method is to simply modify the velocity of all or a re stricted number of atoms in the system such that PAGE 59 59 their mean temperature matches the target temper ature. The kinetic energy of a solid consisting of N atoms in the system is related to the temperature by, T Nk v mB i i i2 3 2 12, (311) where mi and vi are the mass and velocity of atom i T is the current temperature. To set the temperature of the thermostat to target temperature Tt, the velocity of the atom is scaled by, ) ( ) ( t v T T t vi t i (312) A discussion of some of the issues associated with velocity rescaling thermostat is given in Chapter 7 in the context of our development of methodologies for simu lations of radiation damage, along with a discussi on of alternative methods. Constant Pressure Algorithm Materials generally expand on heating. This implies that the size and shape of the simulation cell needs to be adjusted in real istic manner depending on the temperature and external pressure. In our MD simulations, this is done using a constant pressure algorithms.96, 98 There are several available methods depending on the constraint on the pressure. One of the methods called method of extended system first proposed by Andersen.105 To understand this method, imagine a system of volume V to which a piston of mass Q is attached. Thus, with the motion of the piston, the volume of the system changes. If the volume of the simulation cell is V then kinetic and potential energies of the piston are, PV E P dt dV Q E K 2 1 .2. (313) Here P is the pressure of the sy stem. Now define variable, si( t ), for atom i in terms of atomic position xi( t ): PAGE 60 60 ) ( ) (3 / 1t s V t xi i (314) Thus, si( t ) is reduced coor dinate of atom i Then the kinetic and poten tial energies of all the atoms in the simulation cell are, ) ( 2 1 .3 / 1 2 3 / 2i i i is V U E P dt s d m V E K (315) Now using either Hamiltonian or Lagrangian forma lism of classical mechanics, one can arrive at the two couple equations of motion: Q P p dt V d dt dV dt s d V mV t F dt s di i i 2 2 3 / 1 2 23 2 ) ( (316) Here Fi is the force on atom i p is the instantaneous pressure and it can be calculated from, i i i i i iF r V v m p 3 1 3 12. (317) Notice that in this scheme the pressure is isotropic; it means that the pressure is hydrostatic. However, it is easy to allow all the pressure co mponents to evolve indepe ndently. To do so, we need to define matrix, H such that, ) ( ) ( t s H t xb i ab a i, (318) instead of Equation (314). a and b here indicates the Cartesian directions. This expression implies that V =det( H ). Now the equation of motion, Equatio n (316), can be rewritten in terms of H instead of V This will give the constant stress simulation method developed by Parinello and Raman.106 This method has no restriction in the e volution of the compone nts of the pressure or the shape of the simulation cell. This A ndersen (and ParinellRahman) methods are used throughout the work presented in this thesis. PAGE 61 61 Nonequilibrium Thermal Conductivity Simulations The thermal conductivity of a material can be obtained from the nonequilibrium MD simulation. This approach used here is called direct method and was firs t developed by Jund and Jullien.107 The simulation cell used in this me thod is a square cylinder, long in the z direction, and narrow in the x and y directions: see Figure 35. Prior to the thermaltransport simulation, the system is heated to the temperature of interest using a constant temperature, constant pressure (NPT) simulation algorithm for thermal and st rain equilibration. After the system has equilibrated, two regions in the simulation cell are designated for the heat source and heat sink. The simulation cell dimensions are then fixed, and the heat source and heat sink are turned on to create a steady state heat flow. Once the system reaches steadystate, the temperature gradient is determined from an average over some duration of simulation time. The thermal conductivity is then calculated from Fouriers law, J = dT/dx where J is the heat current, is thermal conductivity tensor, and dT/dx is the temperature gradient. Some specific setup and condition such as length of the simulation, heat current, simulation cell size, etc. may depend on the material of interest. In a previous study on silicon, it was shown that the calculated thermal conductivity depends weakly on the cross sectional area 108. The same study on Si also showed that the calculated thermal conductiv ity depends strongly on the length of the simulation box. In Chapter 2, using the elementary kinetic theory for a phonon gas, it was shown that, (Equation 26) l v cv3 1 where cv is the volumetric specific heat, v is the mean sound velocity, and l is the effective mean free path of phonons. The appropriate specific heat above the Debye temperature is given by the DulongPetit value for a noninteracting gas: PAGE 62 62 n k cB v) 2 3 ( (319) where kB is the Boltzmann constant, and n is the number density of the ions. For a cubic lattice with Nc atoms in the unit cell, n = Nc/ a3. vs can be estimated by v =( vL+2 vT)/3, where vL and vT are the longitudinal and transverse sound ve locities. These valu es are given by 11/ C vL and 44/ C vT for a cubic lattice. The effective mean free path can be written as, using Matthiessens rule for the relaxation time 108 as, 1 1 1 BC effl l l, (320) where 1 l and 1BCl are the mean free path of phononphonon s cattering in an infinite media and phononboundary scattering, respecti vely. If the length of the simulation cell is Lz, BCl is approximated by Lz/4 since the heat source and sink are separated by Lz/2. Finally the thermal conductivity can be writ ten as a function of system length as: z C BL l v N k a 4 1 2 13. (321) This means that the best estimate for the infi nite size bulk thermal conductivity can be obtained from MD simulations by the extrapolat ion of the conductivity values to 1/ Lz=0. This approach was previously applied to Si and diamond single crystals 108, and will be used here. PAGE 63 63 Figure 31. The root mean square of the total en ergy as a function of time for velocity Verlet (circles), 4th order Gear (squares), 5th order Gear (triangles), and 6th order Gear (diamond) algorithms. The steep slopes of th e Gear method indicate that it has high accuracy is small time step size is used. [Adapted from reference 96 (Page 83, Figure 3.3).] PAGE 64 64 Figure 32. Two dimensional periodic system in simulation. The actual simulation cell is the central shaded cell of the size L L All the other cells labe led AH are its periodic image. The arrows indicate the molecule 1 is moving out of the simulation cell and coming in from the other size because of its periodic image. [Adapted from reference 96 (Page 24, Figure 1.9).] PAGE 65 65 Figure 33. A system of point charges as a sum of screened point charges and smeared charge density with opposite sign. [Reproduction fr om reference 96 (Page 30, Figure 1).] = + PAGE 66 66 Figure 34. (a) Electrostatic energy per ion as a function of cutoff, Rc. Dash line is used to indicate the Madelung energy of NaCl. Electrostatic energy show no systematic convergence even at large Rc. However, there are a few value of Rc that gives near the Madelung energy. (b) Electrostatic ener gy per ion as a function of the net charge, q in the sphere. Approximate linear relati on is clear between th e energy and the net charge of the sphere. The Madelung energy coincide with point, q =0. [Adapted from reference 99 (Figure 1).]. PAGE 67 67 Figure 35. Simulation cell setup for the direct method for simulating thermal conductivity. The heat source and sink are located at a quarter of the cell leng th away from the center of the cell. The same amount of energy, is added to the atoms in the heat source, and removed from those in the heat sink. Th is setup results in two equivalent heat currents J in opposite directions along the z axis. PAGE 68 68 CHAPTER 4 THERMAL TRANSPORT IN SINGLE CRYSTAL UO2 Molecular Dynamics Simulation of UO2 As discussed in Chapter 1, the thermal c onductivity is one of the most important performance metrics for a nuc lear fuel material. UO2 is used as the primary fuel in nuclear power reactors, and its thermal conductivity has been determined in a number of experiments. Therefore, it is an ideal material to test theore tical models. Once the validity of the model is established, simulations with point defects and GBs are possible. Potential Models The interatomic interactions (Chapter 3) cons ist of a longrange elec trostatic component and shortrange interactions which describe the mate rials specific largelyrepulsive component. In this work, the results obtained from two shortrange interaction models: those due to Yamada 109 and the Busker 110, 111. The Busker potential 110, 111 is a traditional Buckingham potential (Equation 38), 6) / exp( ) (ij ij ij ij ij ijr c r A r V The values of the parameters are given in Table I. One of the advantages of Busker model is the transferability of the potential between a variety of elements 110, 111 it also has parameters for U3+ and U5+, thereby allowing the effects of offstoichiometry to be studied. (Chapter 6) The Busker potential was originally devel oped to study defect properties in UO2 and ZrO2. By contrast, the Yamada model is a BushingIda type potential 112 given by 109: exp 2 2 exp ) / exp( ) (6r r r r D r c r A r Vij ij ij ij ij ij ij ij ij (41) The first two terms are the same as in the Buckingham potential. The last two terms are the Morse term, which provides a covalent component. However, this term is not strictly PAGE 69 69 covalent in the sense that it does not have di rectionality. Due to the assumption of partial covalency, the charges of the ions are given by nonformal values intended to represent the partial charge transfer between ions (Table I). Unlike the Busker potential, the Yamada potential was fit specifically to the UPuO system. To avoid the prohibitive computational expens e associated with the calculation of the electrostatic interaction through the Ewald method, particularly fo r systems with large numbers of ions, the electrostati c interactions are calculated using the direct summation method. (Chapter 3) Table II summarizes the structural parameters and elastic properties determined using the two potentials. By construction, both give good values for the la ttice parameters. The General Utility Lattice Program (GULP)113, 114 is used to determine the el astic properties at 300K. GULP uses a static method based on the quasiharmonic approximation113, 114 and thus provides a slightly lower estimated value for the elastic constants1 than does the direct MD simulation, which includes the dynamical motion of the ions. We can see that while the Busker potential reproduces the value of C12 rather well, it overestimates C11 and C44. By contrast, the Yamada potential gives good estimates of C11 and C44, but severely underestimates C12. As a result one overestimates and the other underestimates the bulk modulus B = ( C11 + 2 C12)/3. There are numerous other interatomic potentials for UO2 in the literature. Recently, Govers et al. undertook an extensive comparison of a number of empirical potentials for UO2, including the two used in this study 115. Their assessment include d both rigidion and shell models, and they calculated the cohesive energy, la ttice parameter, elastic constants, dielectric constants, point phonon frequencies, and defect form ation energies. While some potentials seemed to give a better physical description than others, their results showed that no single PAGE 70 70 potential faithfully reproduces al l of the physical properties of UO2. The two potentials used in here are thus repres entative of other UO2 potentials with regards to their materials fidelity. Thermal Expansion Another metric of the fidelity of a potential is to compare the thermal expansion of the system with experiment. The simulations to determine the thermal expansion are performed using Andersens constant pressure scheme 105. The simulation cell contains 666 unit cells (2592 atoms), and the temperature of the system is controlled by velocity re scaling. In order to reach equilibrium, the simulations are run for 7.5 ps with 0.25 fs time steps. The temperature is varied from 0 K to 2000 K at 100 K intervals. The lattice parameter at each temperature is obtained by the simulation cell vo lume average of the last 5.5 ps of the simulation. Figure 41 shows the normalized lattice parameter as a f unction of temperature; the experimental values are taken from Finks critical assessment of the experimental data 3. Yamadas potential shows good agreement with th e experimental thermal expansion up to about 1000 K, above which it is systematically lower. The Busker potential gives a systematically lower thermal expansion at all temperatures. At low temperatures, the thermal expansions are essentially temperature independent with Yamada = 8.4 106 K1 and Busker=6.2 106 K1. The discrepancy between our calculated value for th e Yamada potential and th e previously published value of 10.1106 K1 109 is at least in part due to the WignerKirkwood correction 116 to the free energy used in analyzing their simulations. Th e thermal expansion calc ulated with the two potentials are both smaller than the experimental value of Expt=11.8 106 K1. The thermal expansion coefficient is a result of the anharmonicity of the interactions in the material. This is encapsulated in the Grneisen relation: PAGE 71 71 B cv3 (42) where cv and B are the specific heat and bulk modulus respectively, and is the Grneisen parameter, which measures the dependence of the phonon frequencies on system volume, and is thus a drect measure of the anharmonicity of the interactions in the the system 60. If the interatomic interactions were purely harmonic, then the Grneisen parameter would be zero and there would be no thermal expansion. Thus th e higher thermal expansion of the Yamada model can be interpreted as a result of higher anharmonicity in the po tential compared to the Busker potential. An estimate of the lattice thermal conductivity in terms of the Grneisen parameter was first given by Leibfried and Schloemann 117, and refined by Klemens 118, T Mv h kB 3 3 24 10 24 ~ (43) Here kB is the Boltzmann constant, h is the Planck constant, v and M are the volume and the mass per atom. The only two materials consta nts that enter into Equation 43 are the Debye temperature, and the frequencyaveraged Grneisen parameter. From this relation, it is clear that the thermal conductivity decrease s with increasing anharmonicity. Through their dependences on the Grneisen pa rameter, we can use Equations 42 and 43 to give a simple relationship between the thermal conductivity a nd thermal expansion, in terms of the Debye temperature, the bul k modulus and the specific heat: B cv 2 3 2 (44) The constant subsumes all of the nonmaterials consta nts in Equation 42 and 43. In classical simulations at temperatures above the Debye temperature (370K for UO2), the specific heat is PAGE 72 72 essentially equal to the DulongPetit value of 3 kB. Also the Debye temperatures for the two potentials are assumed to be the same. Hence we find: B 2' (45) where = 3Cv 2. Using the values of and B determined above for the two pot entials, we then predict the following: 10 1 8 ~ 10 1 10 ~2 Yamada 2 Busker (46) Using the experimental values of B and Equation (45) gives Expt ~3.4 102 which is considerably smaller than the predictions for the Busker and Yamada potentials. Thus, based on this very nave analysis, we expect the direct si mulations with the two potentials to give thermal conductivities for the Busker and Yamada potentials that are 3.0 and 2.4 times larger respectively than the experimental values; as we shall see in the next section, these estimates are quite accurate.. Thermal Conductivity of Single Crystal UO2 The temperature dependence of the th ermal conductivity of single crystal UO2 is calculated using the direct method (Chapter 2). In this approach a heat current is set up; the resulting temperature gradient is identified from which the th ermal conductivity is calculated from Fouriers Law. Using the description given in Chapter 3, th e system is heated to the temperature of interest using a constant temperature, constant volume simulation algorithm before the thermaltransport simulation. We determined that 12.5 ps. is sufficient for thermal equilibration. The heat source and heat sink are then turned on, so as to create a steady st ate heat flow. Once the PAGE 73 73 system reaches steadystate, which takes about ~ 500 ps, the temperature gradient is determined from an average over 250 ps. Thermal conductivity depends weakly on th e cross sectional area (Chapter 2).108 Therefore the crosssection area is set to be 44 unit cells, which is the smallest size that can be simulated with the cutoff values used in these potentials of Rc=1.98 a where a is the lattice parameter. Since our simulation cell size is limited by the number of atoms that the available computational resources can reas onably handle, an analysis of the finite size length of the simulation cell must be taken into account to ob tain the bulk thermal conductivity (Chapter 3). The thermal conductivity of single crystal UO2 between 300 K and 2000 K is shown for the two potentials in Figure 43 as 1/ vs. 1/ LZ plots, in accordance with E quation 321. The results for the two potentials are similar at all temperatures and system sizes. The thermal conductivity for infin ite system size is determined by linear fits to the data in Figure 43. These infinite size limit thermal conduc tivities, which are the be st estimates of the intrinsic thermal conductivity of UO2 described by these potentials, are shown for the two potentials in Figure 44(a). as a function of temperature. Alth ough the error bars overlap, it does appear that the Yamada potenti al gives a somewhat higher esti mated thermal conductivity at low T than does the Busker potentia l, but has a somewhat stronger temperature dependence. This result is reasonably consistent with the only small difference in thermal conductivity predicted by our simple anharmonicity analysis. The higher values for the estimated ther mal conductivity derived using the Yamada potential compared to those previously published for this potential 109 mainly arise from the finite size effects which are particul arly important at low temper atures when the phonon meanfree path is long, and which were not addresses in the earlie r simulations. A secondary effect is that PAGE 74 74 in this study we have not used the WignerKir kwood correction to the te mperature, which has a significant effect below the Debye temperature. As Figure 44(a) indicates, both poten tials give significantly higher thermal conductivities at low temperatures than the e xperimental values. Fi gure 44(a) shows the thermal conductivity corrected according to the anha rmonicity analysis in the previous section; as we can see the match between the simulation a nd experiment is now much better especially above 750K, the temperature range of interest for nuclearfuel applications. It is worth stressing that this agreement is not the result of fitting th e simulation results to the experiments, but comes from a physical analysis of the anharmonicity in the experimental and simulated systems, via quantities that are both computati onally easy to calculate and gene rally experimentally available, even for materials in which thermal transport data are lacking. It thus attr ibutes the majority of the difference between the experimental and simula tion values as arising from the differences in the bulk modulus, which measures the harmonic properties of the syst em, and the thermal expansion, which measures the anha rmonic properties of the system. As discussed in Chapter 3, the anharmonic Umklapp processes lead to the temperature dependence of the thermal conduc tivity. Debye showed that ~ Tn, with n ~ 12.119 Figure 44. is a loglog plot of the data in 44(b). In each case, the thermal conductivity shows powerlaw behavior with temperature. The e xperimental results are fitted by Expt~0.79, while the simulations yield nYamada~1.14 and nBusker~1.30 respectively, which are consistent with the Debye analysis. The systems simulated are structurally much si mpler than the experimental systems; this also contributes to the discrepancy between the simulation and experi mental results. In particular, in the simulations there are no is otopic defects, no offstoichiometry, and no PAGE 75 75 microstructural defects (GBs, dislocations, second phases etc.) Characterization of the effects of polycrystalline microstructures is given in the next chapter; the effects of point defects are analyzed in Chapter 6. Discussion The simulation approaches used here are well suited to the characte rization of the thermal transport properties of electricall y insulating materials such as UO2 in which heat is transported by atomic vibrations. However, even in a relatively poor thermal conductor such as UO2, finite system size effects can lead to a significant underestimate of the thermal conductivity: the systematic variation in system size coupled with a finitesize scaling an alysis does appear to offer a viable method for taki ng these effects into account. The power low behavior of thermal conductivity from experimental data givesan exponent less than unity; this is in contradiction with the Debye model. The origin of this discrepancy is not clear. One possible origin of the disagreeme nt is the effects of point defects or grain boundaries, although grain size in typical reactor fuel may be too large to show significant influence. It is well known that defects can reduce the thermal conduc tivity and temperature dependence. We will return to this issue in later chapters. Neither potential can qua ntitatively match the experimental thermal conductivity without the consideration of the anharmonic effects. Based on the analysis of Govers et al. the two potentials used in the study are of a similar level of materials fidelity as others in the literature.115 The simple relationship, Equation (45), relati ng the elastic properties, thermal expansion and elastic constants, suggests that a potential with correct elastic properties and thermal expansion, should well reproduce the thermaltr ansport properties. A truly ge neral purpose potential would also need to well reproduce the po int defect properties. Govers et al.115 showed that none of the PAGE 76 76 twentyone potentials they exam ined could satisfactorily reprodu ce the formation and migration energies. There is thus considerable n eed for potentials whic h better describe UO2. PAGE 77 77 Table 41. Parameters of interatomic potentials Busker Yamada O2O2U3+O2U4+O2U5+O2UU OO UO Aij (eV) 9547.96 1165.65 1761.7752386.42 442.208 2346.149 1018.571 ij () 0.2192 0.3786 0.35643 0.3411 0.32 0.32 0.32 Cij (eV6) 32 0 0 4.1462 0 Dij (eV) 0 0 0.7810 ij (1/) 0 0 1.25 R*ij () 0 0 2.369 ZU (e) +4 +2.4 ZO (e) 2 1.2 PAGE 78 78 Table 42. Lattice parameter a nd elastic constants at 300 K. GULP MD Busker Yamada Busker Yamada Experiment3, 23, 120122 a () 5.481 5.482 5.479 5.481 5.478 C11 (GPa) 526 409 547 418 389396 C12 (GPa) 118 55.0 119121 C44 (GPa) 118 53.4 59.764.1 B (GPa) 257 174 209213 PAGE 79 79 0.995 1.000 1.005 1.010 1.015 1.020 1.025 1.030 0500100015002000 T [K]a(T)/a (0) [] Experiment Busker Yamada Figure 41. Normalized lattice parameter as a function of temper ature from experiment 3, Busker and Yamada potentials. PAGE 80 80 0 20 40 60 80 100 120 0500100015002000 T [K]cv [J/molK] GULP MD Experiment DulongPetit Figure 42. Specific heat of UO2 from GULP, MD and experiment.3 The DulongPetit value is shown as a reference. PAGE 81 81 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.000.020.040.06 1/ Lz [1/nm]1/ [mK/W] 300 K 500 K 750 K 1000 K 1500 K 2000 K 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.000.020.040.06 1/ Lz [1/nm]1/ [mK/W] 300 K 500 K 750 K 1000 K 1500 K 2000 K Figure 43. Size dependence of thermal conductiv ity for both the Busker (a) and Yamada (b) potentials. a b PAGE 82 82 0 5 10 15 20 25 30 35 40 45 0500100015002000 T [K] [W/mK] Experiment Busker Yamada 0 5 10 15 20 25 30 35 40 45 0500100015002000 T [K] [W/mK] Experiment Busker Yamada Figure 44. Thermal conductivity of UO2 from the compilation of experimental data by Fink and simulations using the Busker and Yamada potentials (a) before and (b) after the anharmonic correction PAGE 83 83 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 2.402.602.803.003.203.40 log( T )log( ) Experiment Busker Yamada Figure 45. The same thermal conductivity data as in Figure. 44(b) in loglog scales. PAGE 84 84 CHAPTER 5 THERMAL TRANSPORT IN POLYCRYSTALLINE UO2 Introduction Grain boundaries offer a significant obstacle to the transport of heat in phonon conductors as discussed in Chapter 3. Grai n boundaries are present in typica l fission reactor fuel material. (Chapter 1) Although the averag e grain size in such materials may be too large (order of 10 m)3, 40 to significantly influence thermal transport in UO2, finer grains (~150 nm on average) are also present in high burnup structures.40, 41 However, there has been no systematic study of the grain size dependence of th ermal transport properties in UO2. In this chapter, the thermal conductivity of a model fineg rained polycrystal of UO2 is determined. Predictions for the grainsize dependence of its thermal conductivity are also made. Structure of Model Polycrystal Experimental grain sizes of th e range of tens or hundred of microns are not accessible to MD simulation, since each grain would contain ~1013 1016 ions. In this study considerably smaller systems with grain sizes from 3.8 6.5 nm are considered; these small sizes maximize the area of the GB in the system, there by amplifying the interfacial effects. The polycrystalline structures used in the simulations consisted of 24 hexagonal columnar grains. When constructing the polyc rystalline structure, identical closepacked hexagons are arranged to form a completely peri odic structure. Each hexagon is filled with singlecrystal UO2 oriented with [001] along the columnar direction. The inplane orientations are chosen in such a way that the GBs between the grains are high energy tilt GBs, which ensures that the microstructure is stable agains t coarsening during the simulation. Because of the way each grain is constructed, initially there are always a small number of atoms in the GBs which are extremely close to each other: See Figure 51. To address this issue, if any two ions PAGE 85 85 are closer than 1.5 (66% of th e nearest neighbor distance between U4+ and O2, 2.29 ), one of the atoms is removed. This removal of atoms is carried out with care to ensure the charge neutrality of the entire system. Once the structur e is created, it is quenched at 0 K to equilibrate all of the atom positions and to eliminate any inplane stress on the system. It is found that no ions have anomalously high ener gies, indicating that the bonding in the sy stem is physically reasonable. The system is then annealed w ith a constantpressure, constanttemperature simulation at 2000 K and slowly relaxed to 0 K to ensure that the structure is equilibrated. Figure 52 shows the final relaxed polycrystalline UO2 structure for a grain size of 3.8 nm. The columnar microstructure used here allo ws the simulation cell to be thin along the columnar direction. The cutoff to the potentials is 10.4 which would allow the thickness to be as small as 4 unit cells: in our simulations, 5 unit cells thickness is used so as to minimize any effects of the system size in that direction. Since all the grains are equiaxial and of equa l size, it is easy to calculate the GB area and volume fraction. For a grain size of 3.8 nm, the area of the GB is 439.10 nm2. Assuming that the GBs have a thickness of one unit cell, the volume which the GB region occupies is approximately 30 % of the entire volume. The st ructural disorder at the GBs leads to a total volume expansion of 5.5 % and 5.3 % for Busker and Yamada potentials respectively. The corresponding average GB energies are 2.73 J/m2 and 1.89 J/m2, which are consistent with the GB energies of other ceramic materials.123, 124 After the preparation of the equilibrated stru cture, the thermal expansion of these models of polycrystalline UO2 were determined at 300 K. The values obtained were 7.57 106 K1 with the Busker potential and 8.92 106 K1 for the Yamada potential, which are almost PAGE 86 86 indistinguishable from the correspondi ng bulk single crystal values of 7.50 106 K1 and 8.83 106 K1. Thermal Conductivity of Polycrystalline UO2 The thermal conductivity of polycrystalline UO2 is calculated using the direct method described in Chapter 2. Figure 53 show s the temperature dependence of for a polycrystal with a grain size of 3.8 nm, as described by both potentials. These calculated thermal conductivities are considerably lower than those of the corresponding the perfec t crystals, attesting to the significant resistance of GBs to the flow of heat through the system. Unlike single crystalline UO2, the effects of the finite size of the simulatio n cell itself are small in polycrystals because the GB scattering dominates over the phononphonon scat tering; simulations ha ve directly shown that the effect of the size of the simulation cell in polycrystalline MgO is weak125,. Since UO2 has a significantly lower singlecrys tal thermal conductivity than MgO (UO2:7 W/mK and MgO 40 W/mK at 300 K)3, 126 the effect is expected to be even less important in UO2. The thermal conductivity for the polycrystal is considerably higher for the Busker potential than for the Yamada potential. Since the GB energy is a measure of the structural disorder at the interfaces, we would expect that the higher the energy asso ciated with the GBs, the higher the interfacial thermal resistance would be, and th e lower the thermal conductivity of the polycrystal would be. The GB energies obtained using the Bu sker and Yamada potentials are 2.73 J/m2 and 1.89 J/m2, which appears to be inconsistent with this argument. However, the Busker potential is a full charge model, while th e Yamada potential is a partial charge model, resulting in cohesive en ergies of 104.482 eV/UO2 and 45.54 eV/UO2 respectively. Thus, when normalized to the bulk cohesive en ergies, the GB energies are 0.00262 and 0.0016 2 for Yamada and Busker, respectively. That is, when described by the Busker potential, the GBs PAGE 87 87 actually offer less of an obstacle to heat transport than for the Yamada potential, which is consistent with the higher ther mal conductivity for the Busker potential than for the Yamada potential. The ensembleaveraged interfacial (Kapitza) re sistance of the GBs in the polycrystal can be extracted from the thermal conductivity using simple models. There have several effective medium models proposed, includi ng those of Nan and Birringer 127, Yang et al. 78, and Amrit.128 The model by Nan and Birringer assumes that the grains are iden tical ellipsoids and embedded in a continuum medium. The orientations of the grains can be random or aligned. The ratio of the single crystal thermal conductivity to that of the polycrystalline solid of identical columnar grains is, d lK 10 (51) d is the grain size and lK is called Kapitza length given by the ratio of the single crystal conductivity over the interfaci al conductance. The physical interpretation of lK is as the thickness of perfect crystal that would offer the same thermal resistance as the interface; thus a long Kapitza length corresponds to a high interfacial resistance. The model by Yang et al. assumes a one dimensional superlattice type structure. They assume that there is a temperature drop at the interface. The ratio, 0/ is then given by, d lK10, (52) where is a correction term, which accounting for th e conditions where the wavelength of the dominant phonon mode approaches the grain size. This can be taken safely to zero for high temperatures. PAGE 88 88 Amrits model assumes the polycrystalline struct ure consists of identic al square grains. This model assumes that the GBs along the dir ection of heat flow do not contribute to the thermal resistance. 0/ in this model is, d l n nK1 10 (53) Here n is the number of grains in the region of interest. Although these three approaches assume different models for the polycrystalline structures, the resulting expressi ons for the interfacial conductance are rather similar and none of them is obviously better than the others. We will come back to this point below. Most importantly the qualitative trend of the interfacial conductance will not be affected by the choice of the model. Thus, for this analysis, the model by Yang et al 78 is adopted for its simplicity. In, Yangs model (Equation 53), th e interfacial conductance, GK, is then given by: 0 01 d GK, (54) since lK= 0/ GK: The resulting values of GK are given in Figure 54.. Both potentials show moderate increases of interfacial conductan ce with temperature. This temp erature increase is consistent with the results of simulations of other interfacial systems 129 and, more importantly, with trends in experimental data for various systems.130 The physical origin of the increase in conductance with temperature can be understood in terms of the properties of the GB s. As the temperature increases, the anharmonicity of the interactions among the atoms is probed more strongly. While in the perfect crystal, the re sulting scattering lowers the th ermal conductivity, this anharmonic scattering more strongly couples modes across the interfaces, leading to be tter interfacial thermal transport. PAGE 89 89 The interfacial conductance can be recast in terms of the Kapitza length. As shown in Figure 55, the Kapitza le ngth decreases strongly with increasing temperature; this is a result of the decrease in the thermal c onductivity and the incr ease in the interfaci al conductance with temperature. For both potentials, the Kapitza length is signi ficantly larger than the grain diameter, particularly in the low temperature region. This is an indication that the thermal transport in our model system is dominated by the GBs. It also suggests that the funda mental assumption of separable and grainsize bulk a nd interfacial thermal properties used for the analysis may be violated at these small grain sizes There is thus a clear need fo r a model that treats the GB and grain interiors in an integrated manner. The effects of grain size on thermal conductivity have been investigated for grain sizes up to 6.5 nm grains. The inset to Figure 56 shows the increas e in thermal conductivity with increasing grain size determined from the simulati ons; this increase arises from the decrease in the relative volume of GBs in the system. Our data for polycrystals simulated with the Yamada potential are fit to the model of Yang et al. thereby allowing the thermal conductivity for large grain sizes to be estimated. Taking the bulk single crystal co nductivity of 15.2 W/mK at 300 K, the fit gives the Kapitza conductance of 0.15 GW/m2K, which is close to the value previously determined for the 3.8 nm polycrystal (see Figure 55). For the sake of comparison, the same fit to the grain size dependence are done using the models by NanBirringer and Amrit. lK are 8.4 nm for NanBirringer model and 12.6 nm for Amrit model, respectively. This compares well to 7.4 nm from the Yangs model. Discussion The analysis of the polycrystal to obtain the interfacial conductance is not unambiguous, since the calculated Kaptiza lengt hs are larger than the grain sizes. However, the fact that the PAGE 90 90 calculated thermal conductiv ities of polycrystals of different grains sizes, albeit over the narrow range of 3.8 6.5 nm, can be fit to the Yang et al model, suggest that their use is not unreasonable. Moreover, since the same analysis is used throughout, the trends of interfacial conductance with temperature are reasonable. Comparison among the three effective medium m odels showed fairly similar results. Although the Amrits model predicts a somewhat highe r Kapitza length, it is clear that the cause is due to the way his model includes the number of interfaces in between the heat source and sink. More importantly, however, this does not ch ange the qualitative trend of the grain size dependence although the absolute value of the cal culated Kapitza length, i.e. the interfacial conductance may be different. The interfacial conductance obtained from our simulations in fact compares well with those of the other heterogene ous interfaces: see Figure 57., indicating the validity of the methodology. PAGE 91 91 0 5 10 15 20 25 30 35 40 45 50 02468 r ()g(r) Original Removed Relaxed Figure 51. Pair distribu tion function of the polycrystalline stru cture. The original structure is the unrelaxed reference structure. Origi nal: All atoms are in perfect crystal like positions except the ones at the grain boundaries. Removed: The original structures has a number of atoms in the grain boundary region, which are so close to each other that they must be removed. Relaxed the structure after annealing and thermal equilibration. 0.0 0.2 0.4 0.6 0.8 1.0 01234 r () PAGE 92 92 Figure.52. The final polycrysta lline structure used in the th ermal conductivity calculations. Filled and open circles indicate uranium and oxygen ions, respectively. The view is along [001] in the fluorite cr ystal structure. The PBCs are applied in all three directions. PAGE 93 93 0.0 0.2 0.4 0.6 0.8 1.0 1.2 0500100015002000 T [K] [W/mK] Busker Yamada Figure 53. Thermal conductivity of 3.8 nm grain polycrystalline UO2 from simulation for the two potentials. PAGE 94 94 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0500100015002000 T [K]GK [GW/m2K] Busker Yamada Figure 54. Average thermal conductance for the [001] tilt GBs from 3.8 nm grain polycrystalline UO2. PAGE 95 95 0 10 20 30 40 50 60 70 80 90 0500100015002000 T [K]l [nm] Busker Yamada Figure 55. Kapitza length from 3.8 nm grain polycrystalline UO2 simulations. PAGE 96 96 0 2 4 6 8 10 12 14 16 02004006008001000 d [nm] [W/mK] Figure 56. Grain si ze dependence of the thermal conductivity of polycrystalline UO2 using Yamadas potential. The solid lin e is the fit of the model by Yang et al.78. Inset is the same plot for nanometer scale indicating the fit with our simulation data. 0.0 0.2 0.4 0.6 0.8 1.0 1.2 3.04.05.06.07.0 d [nm] [W/mK] PAGE 97 97 Figure 57. Interfacial conducta nces of selected interfaces as function of temperature. UO2 data point is from our MD simulation. DMM a nd LD are the theoretical predictions for the Al/Al2O3 interfaces. The right axis is the equivalent film thickness assuming the bulk conductivity of 1 W/mK. [Adapted from reference 130 (Figure 1).] UO2 PAGE 98 98 CHAPTER 6 THERMAL TRANSPORT IN UO2X Introduction Just as GBs do, point defects also act as scat tering centers in crystalline solids and reduce the mean free path of phonons. This results in a reduction of the ther mal conductivity of the solid. The thermal conductivity in a single crystal insu lator can be very well described by the dynamics of phonons as discussed in de tail in Chapters 2 and 4. In the case of highly disordered solids, such as highly doped materials or amor phous solids, the situation is quite different (Chapter 2 for general discussion). In practical applications, UO2 often contains a large number of point defects, dislocations, and GBs 11. Although understanding the mechanisms of how these defects couple to affect thermal transport is currently prohibitively complicated, understanding the effect of individual types of defect is possible. Therefore, in this chapter, the thermal transport properties of nonstoichiometric UO2 are characterized by simulation, and the results analyzed in terms of the theory of thermal transport in a disordered solid. While the main focus is on the anion defects, which are the dominant defects in UO2, the effects of other point defects are also discussed. The effects of is otopic defects are also briefly analyzed. Point Defects in UO2 UO2 contains a wide variety of point defects when it is used for the nuclear fuel application. The extreme conditions of temperat ure and radiation condition generates intrinsic, extrinsic (fission products), and electronic defects. In pure UO2, the possible in trinsic point defects are VU, UI, VO, OI, the Frenkel pair, the antiFrenkel pair, and the Schottky trio. The defect formation energies for intrinsic defect determined by experiments and density functional theory (DFT) calculations are s hown in Table 61. Agreement between the experimental and DFT values are reasonable except for the Schottky def ect. It can be deduced from this data that PAGE 99 99 the thermodynamically most relevant intrin sic defects are the antiFrenkel defects.3, 3436 Although the Schottky trio formation energy is hi gher than antiFrenkel pa ir, it also plays an important role at high temperature 3, 35. The most recent DFT calculations, by Freyss et al. 131 and by Crocombette et al. 132, determined the defect forma tion energies of individual OI, VO, UI, and VU, which are very difficult to determine fr om experiments. They found a negative formation energy for OI, which implies that UO2 readily incorporates oxygen ions under oxidizing condition. So far the DFT and experi ments show qualitativel y consistent results despite some quantitative discrepancies. Empirical potentials have been also used to evaluate the defect formation energies of intrinsic defects. Govers et al.115 calculate the Frenkel, antiFrenkel, Shottkey defect formation energies using 18 existing empirical potentials. Figure 61 gives a comparison of the normalized Frenkel and antiFrenkel pair formation energies from Govers work. An additional data point for the the rigidion Busker potential (used in our simulations of singlecrystal and polycrystalline UO2) is included and referred to as Busker1 in the figure. The error bars for the reference point indicates the scatter in the ava ilable DFT and experimental data. Most potentials do a fairly good job in reproducing the oxygen anti Frenkel pair formation energies, except for a few outliers. However, all the empirical models overestimate the uranium Frenkel pair formation energies. Intrinsic defects in UO2 often have rather complex structures. UO2 has a broad hyperstoichiometric phase at high temperature (See Figure 17.) because of the low defect formation energy of OI. Willis24, 25, 133135 and later Hutchings22 characterized the defect structures in UO2+ x, by neutron diffraction. Willis neutron diffraction data on UO2.13 showed that the excess anions were not at the body center interst itial site: rather, th e unit cell structures PAGE 100 100 were distorted and formed a complex of defects including oxygen vacancies25, 133. He argued that the dominant structural defect is the ani on Frenkel defect. A similar observation was made in fluorite (CaF2) by Cheetham et al.136, and the result was interprete d as an 2:2:2 interstitial cluster (Figure 62.). The 2:2:2 notation indicat es that the cluster consists of 2 vacancies, 2 atoms shifted in along <110>, and 2 shifted along < 111>. This defect cluster structure is often refereed to as the Willis cluste r in the literature. Hubbard et al. 137 performed optical absorption spectroscopy experiments and explai ned the change in the charge st ate of cation sublattice as due to the formation of the defect clusters. Allen et al.138, 139 also observed the def ect clusters in their diffraction data and proposed orde red defect substructures for UO2+x phases. There are several types of defect clusters in UO2+ x reported such as 2:1:2 and 4:3: 2, which can combine to form larger defect clusters.137, 138 A UO2x phase also can occur at high temper ature and low oxygen partial pressure.36 However, this UO2x phase occur more commonly by the formation of oxygen vacancies through extrinsic defects.34, 140 Catlow undertook the first general theoretica l survey of defects properties in UO2.34 He used an empirical model derived from experimental data of UO2 bulk properties and HartreeFock calculations. The defects considered were OV, IO UV IU, OV, IO Frenkel and antiFrenkel pairs, and the Schottky trio. His work on formation energies of these defects established that the antiFrenkel pair of a doubly charged interstitial and a v acancy is the dominant intrinsic atomic disorder in UO2. This result is in agreement with the neutron diffraction studies by Willis24, 25, 134, 135 although the charge state of the defects is not known in the experiment. Catlow also found that the lattice is distorted in the hypers toichiometric structure a nd that the interstitial is slightly shifted in <110> at the lowest energy state. Inte rstitials shifted along <100> and PAGE 101 101 <111> are also lower in energy compared to th e perfectly symmetric body centered position, but none of these shift are significant enough to account for the 2:2:2 defect cluster structure by Willis. In addition to the atomic defects, the elec tronic defects were also considered and turned out to be much more favorable than such atomic disorder.34, 141, 142 When used as a fuel in a LWR, UO2 also contains a certain leve l of isotopes of uranium. 235U is the fissile elements in UO2 and produces the neutrons necessary for the chain reactions. The dominant isotope is 238U; however it is not fissile. The concentration of 235U is controlled and in a typical fuel is about 34 wt.%. This isotope acts as a mass defect in UO2 and can scatter phonons in the crystal. The isotope effect is im portant in the thermal transport properties of materials such as diamond51 and Si143, 144. However, isotopic effects in the UO2 fuel thermal transport properties have not been investigated. In addition to the intrinsic defects, the effect of extrinsic defects such as fission products is critical.11 As discussed in Chapter 1, fission produc ts consist of a wide variety of elements having different masses and sizes.37 Some of them are even radioactive.11 They can reside in the UO2 matrix as gases42, solid solutions36, 145, or different phases36, 146. Although they can have significant effects on the UO2 properties including thermal transport40, the study of the effects of fission products is beyond the focus of this work. Thermal Transport in UO2 with Point Defects One of the earliest studies on the effect of stoichiometry on ther mal transport of UO2+x was undertaken by Hobson et al.147. They measured the thermal diffusivities from 550 to 2500 K by laser flash on samples with between 93 and 96 % of the theoretical density. The degree of offstoichiometry studied was x =0.006, 0.030, and 0.060. The thermal conductivity was determined from the measured thermal diffusivity and mass density of the sample, with the heat capacity data taken from the work by Affortit and Macron.148 In general, the ther mal conductivities of all PAGE 102 102 the samples showed monotonic decreases with increasing temperature, but there was a sudden change in the slope around 700 K for high concentration of oxygen interstitials at x =0.060. This change was later ascribed to the formation of U4O9 phases in the sample. Lucuta et al.149 determined the effect of thermal conductivity of UO2+ x between 298 and 1773 K. Thermal conductivities were determined by measurements of thermal diffusivities, heat capacity and mass density. All the samples used were 98 % of the theoretical density, but the conductivity values for 100 % dens e material were also estimated from the data. The thermal conductivity of UO2+ x was found to decrease with increasi ng degree of offstoichiometry, with reductions for x =0.007, 0.035, and 0.084 of 13, 37, and 56 % at 827 K, and 11, 23, and 33 % at 1773 K. Just as Hobson did, Lucuta et al. also observed a change in the behavior of the thermal conductivity for x =0.035 and 0.084. By Xray diffraction obs ervation, they determined that this effect as due to the precipitation of U4O9 phase. Even with careful pr eparation, the formation of the secondary phase in their samp les could not be avoided. Lewis and his collaborators undertook both e xperimental and nume rical studies on the oxidation behavior of UO2 in a defective fuel rod, and on th e effect of oxidation on the thermal conductivity of the fuel.150 They employed a phenomenological model for thermal conductivity as a function of degree of offstoichiometry. They showed that depends significantly on the value of x within the range of 0 x 0.2 below 2000 K. The contri butions from the radiation (Chapter 7) and polarons (Chapter 1) hardly made any difference in this temperature range and were not affected by the oxygen interstitia ls introduced in the material. There is also a recent MD simulation study performed by Yamasaki et al. 151 on the thermal conductivity of UO2+x. They developed a rigidion t ype interatomic potential based on the thermal expansion due to the incorporation of OI using the experimental values of Grnvold PAGE 103 103 152. They employed the GreenKubo (GK) method to calculate the ther mal conductivity by an equilibrium MD simulation. Their value of th e thermal conductivity was in good agreement with experimentally observed conductivity data. In particular, their data shows a decrease of conductivity with increasing concen tration of defects, the same trend as observed in the experiments 149, 153. However, the direct integration of the GK formula is known to result in ambiguous results87, 154156. In the case of si ngle crystal Si, Shelling et al.87 showed that the simulation must be run for longer th an 1 ns, and the time integral must be performed for at least 200 ps. Although UO2+ x has lower thermal conductivity than Si (experimental value of the isotopically enriched Si is ~200 W/mK143), it is not clear how long the necessary simulation time and integration interval must be for UO2+ x. Furthermore, Yamasaki et al constructed the structures based on the 2:2:2 clusters, such de fect structures are not consistent with the assumption of homogeneity in the GK method. Molecular Dynamics Simulations The theoretical stand point of view, there has been no reliable simulation work on the effect of stoichiometry on the thermal conductivity of UO2. In the following, the approach taken to perform thermal conductiv ity calculations of UO2 x using MD simulations is described. The range of offstoichiometry investigated was between x =0.02 and 0.25. This range was chosen based on the phase diagram of UO2 x between 800 and 1600 K (See Figure 17). Although this range of x is beyond the range of offstoichiometry studi ed in any of the previous experimental and simulation works, this broad range enhances th e effect of the defects to thermal transport and clearly indicates the concentrati on dependence of the thermal conduc tivity. Implications of the outcome of the simulations are discussed in te rms of the thermal conduc tivity of a disordered solid. PAGE 104 104 Preparation of the Structures In order to perform MD simulations with oxyge n defects, interatomic potentials must be modified to retain the charge neutrality. Sta ndard MD simulations do not allow the charges of the atoms to vary during the simulation. Therefore, in order to include O2and OV, interatomic interaction for U3+ and U5+ must be added The potential para meters are given in Table 41. These parameters were provided from Prof. Grimes in Imperial College London. Note that there is no shortrange interacti on between the cations. The UO2+x structures are prepared in the follo wing fashion. First a defect free UO2 structure is prepared in the desired size. Then a number of excess oxygen atoms corresponding to the desired concentration are in serted in randomly selected octahe dral (4b) sites. In addition, randomly selected U4+ atoms are replaced with U5+ atoms such that overall charge neutrality is maintained. Once the structure is created, the simu lation cell is heated to 3000 K for 10 ps. This allows the oxygen atoms to diffuse and minimizes the possibility of any sort of artificial configuration. After the high temp erature equilibration is complete the system is slowly cooled to 0 K in a step wise fashion at an average cool ing rate of 30 K/ps. Finally the structure is quenched by steepest decent al gorithm to ensure the minimu m energy configuration. A similar procedure is performed to create UO2x structures. From the single crystal UO2, an appropriate number of the randomly selected oxygen atoms are removed to create a desired concentration of oxygen vacancies. At the same time, randomly selected U4+ atoms are replaced with U3+ atoms for charge neutrality. After the st ructure is created, the same annealing and relaxation procedure is performed to achie ve the minimum energy configuration. Characterization of the oxygen excess and defi cient structures is performed by the pair distribution function analysis. Figure 63 shows the pair di stribution functions (PDF) of PAGE 105 105 UO2+x.between uranium and oxygen ions. The PDF of the defected st ructure before annealing is shown as a reference for the ideal single crys tal separation of uranium and oxygen ions. The PDF before annealing shows secondary peaks by the primary peaks. The primary peaks are the O2ions at 8c positions in the fluorite structure. These secondary peaks are the locations of oxygen interstitials at the octahedral (4b) sites. We also note that before the structure is relaxed, all O2ions are at 8c sites such that there is no difference between the PDF between U4+O2and U5+O2. After annealing and relaxing at 0 K, the peaks are broadened and their heights decreased. The differences in PDF between U4+O2and U5+O2appear. The secondary peaks corresponding to U5+O2are merged with the primary peaks. Furt hermore, the first peak of the PDF of U5+O2is shifted inward compared to the same peak of U4+O2by ~0.2 This shift is expected from the stronger electrostatic attraction between U5+O2than U4+O2. The higher order peaks are also shifted, but the differences ar e so small that they cannot be seen in the figure. In the UO2x, similar observations can be made. After relaxation, the peaks are also broadened and the heights are reduced. Although the U4+O2peak does not move, the U3+O2peak is shifted to larger distan ces. Just as in the case of UO2+ x, this is expected by the reduction of the electrostatic attraction between U3+O2. Chemical Expansion of UO2x Using the method described in the previous section, UO2 x structures are prepared, and are relaxed at 0 K. The range of offst oichiometry for this study is from UO1.98 to UO2.25. This range is chosen based on the phase diagram of UO system (Figure 17.). Although, in the actual reactor fuel, such a wide range of stoichiometry does not occur, this broad range of defect concentrations will enhance the in fluence on the thermophysical prop erties in our simulation. These structures show a change in lattice pa rameter due to the excess or deficiency of oxygen ions in the system. This is effect is know as the chemical expansion or contraction. PAGE 106 106 Figure 65 shows that the norma lized lattice parameter of UO2+ x at 0 K decreases with increasing concentration of oxygen interstitial s. The lattice parameter at th e highest defect concentration ( x =0.250) is ~0.5 % smaller than the la ttice parameter of stoichiometric UO2. This is the result of the increased electrostatic attractions by the competition between the contraction due the replacement of some of the U4+ ions by U5+, and the expansion from the addition of extra O2in the octahedral sites. This indi cates that the expansion of the la ttice due to the strain of the oxygen interstitials at 8c site s is overwhelmed by the strong electrostatic attraction. Grnvold152 performed Xray studies to dete rmine the thermal expansion of UO2+ x. His data indicates no significant diffe rence in the lattice parameter from that of stoichiometric UO2 over the range of x = 0.0~0.2, i.e, no chemical expansion. However, the phase diagram of UO2+x is rather complex and is stil l an active area of research.157 Furthermore, there has been no systematic study on the chemical expansion of UO2+ x after Grnvold, and more recent work is needed to draw a conclusion as to the pres ence or absence of a chemical effect. The UO2x phase is analyzed in the same fashion. Figure 66 gives th e normalized lattice parameter of UO2x at 0 K after relaxation determined fr om simulation. In contrast to UO2+ x, the lattice parameter increases with increasing concentration of def ects. The change in lattice parameter is linear because of the narrow range of defect concentrations. There is no experimental data on the chemical expansion of UO2x available due to experimental difficulties and due to the limited relevance to the applications. Thermal Expansions of UO2x Figure 67 gives the thermal expansion of UO2+ x obtained from the MD simulations. For each composition, the data is normalized by the 0 K lattice parameter, such that the effect of chemical expansion does not play a role. The temperature dependen ce of the lattice parameter is almost linear for all the con centrations of defects studied. The thermal expansion, indicated PAGE 107 107 by the slope of the curves, decreases from 6.8106 K1 to 5.4106 K1 as the defect concentration increases. Figur 68 is the plot of vs. x extracted from the data in Figure 67. in UO2+ x. It shows the gradual linear decrease of with increasing x Again, there is no systematic study on the effect of defect concentration on except for the Grnvolds work, which showed no effect of composition. The thermal expansion of UO2x is shown in Figure 69. The data indicates that the defect concentration has a very minor effect on the ther mal expansion of the lattice. This is the case because the low concentration of the defect introduced. Thermal Conductivity of UO2+ x Using the technique described in the beginning of this section, UO2+ x structures of 21.921.9 2 (44 unit cells) cross sections are prepared This is the same cross section size as the ones used in Chapter 4 for the single crystal UO2 study. The effect of the simulation cell length is still important and must to be inves tigated to obtain the infinite size bulk thermal conductivity. Therefore, the si mulation cell lengths of 262.5, 350.0, 437.4, and 524.9 were used. The simulations were run at 800 and 1600 K, which correspond approximately to the surface and centerline temperatures of a fuel pellet during the operation of a typical PWR. The rest of the simulation conditions were the same as in the single crystal UO2 (Chapter 3). The results are shown in Figure 610 and 611 as 1/ vs. 1/ Lz. The data points of the each defect concentration were fit to the straight lines expected from the kinetic theory of phonon gas (Equation 310). At 800 K, the data points fit well with the straight line for the defect concentrations. In particular, the quality of th e fit is better at the low concentrations ( x 0.75) with the correlation of fit being better than 0.95. For x>0.75, the correlation varies between 0.60 and 0.78. This is reasonable since Equation 310 assumes independent phonon scattering which PAGE 108 108 is most likely true at low concentrations of defect s. However, at higher concentrations of defects ( x 0.0125) at the same temperature, the data points do not follow a linear fit as we ll as that of lower concentrations. For x 0.175, the conductivity values are st atistically identical and error bars overlaps. At 1600 K (See Figure 611), the quality of fit is somewhat poorer ( r2 0.92 for x 0.75 and 0.62< r2 0.73 for x >0.75) than at 800 K, and the errors are slightly larger. At this temperature, the inverse conductivities are statis tically indistinguishable from each other for x 0.125. Slopes of the linear regressions for 800 K and 1600 K are rather si milar, indicating that the size dependence is not strongly dependent on the temperatures between 800 K and 1600 K. From the linear fit to the data in Figs. 610 and 611, the thermal conductivities of infinitely large system at 800 and 1600 K are obtai ned by extrapolation. The results are shown in the top of Figure 612. The data for bot h 800 and 1600 K show monotonic decreases of the conductivity with increasing concentration of oxygen interstitials. Interest ingly, the two curves converges to the same value of th e conductivity of ~3.6 W/mK above x =0.125. Besides the trend of th e thermal conductivity on OI concentration, there is an apparent discrepancy in the values of thermal conductivity rela tive to the experimental data by Lucuta 149. The thermal conductivity from the MD simulation is much higher than that of experiment by a factor of roughly 23. This difference is mo st likely due to the quality of the interatomic potential (Chapter 3 for the di scussion in single crystal UO2). Using the same scaling factor ,0.34, from the ratio of the square of the thermal expansion and bulk modulus between the experimental and simulation, we obtain the ther mal conductivity curves shown in Figure 613. The resulting rescaled conductiviti es agree well with experimental This suggests that this treatment is appropriate for the simulation of the effects of point defects on the thermal conductivity. PAGE 109 109 Thermal Conductivity of UO2x The thermal conductivity of UO2x was calculated using exactly the same methods as for UO2+ x. Figures 614 and 615 ar e the inverse thermal conduc tivity at 800 and 1600 K as a function of simulation cell length for th e oxygen vacancy concentrations of 0.00 x 0.02. Just as in the case of UO2+x, the simulation size effect is investigated The quality of the fits is fair at both temperatures ( r2>0.91 at 800 K, and r2>0.82 for 1600 K), and there is no significant difference between the two data sets. The predicted bulk thermal conduc tivity obtained from the linear fits is shown in Figure 616. The thermal conductivity values at both 800 and 1600 K decrease al most linearly with concentration. This is reasonable considering th e low concentration of defects. The linear fits give the slope of 122 W /mK at 800 K and 53 W /mK at 1600 K. (Note that the unit of concentration is hidden for simplicity.) The sm aller slope for 1600 K is an indication of the weaker concentration depe ndence because of the signi ficant phononphonon scattering. Unfortunately there is no reliable thermal c onductivity data from experiments to compare the data from the MD simulations. Finally the thermal conductivity at lo w defect concentrations in both UO2+ x and UO2x are characterized. Figure 617 show s the thermal conductivity for UO1.98.UO2.02 at 800 K and 1600 K. Although UO2+ x data is linear, UO2+ x data show some curvature. For analysis the data set from UO2+ x is fit to tangent line at x =0.000. The slopes of the tangent for 800 K and 1600 K data are 141 W/mK and 64 W/mK, respectiv ely, similar to the slopes of UO2x. Isotopic Effect The effect of uranium isotopes in UO2 was also studied. As desc ribed in Chapter 1 and the beginning of this chapter, UO2 in a typical LWR is isotopically enriched to increase the fissile 235U. The concentration of 235U may vary depending on the specifica tion of the reactor, but it is PAGE 110 110 typically in the range of 34 wt. %.11 In the MD simulation, 235U was introduced as a simple mass defect and the concentration is set to 5 at. %, which is comparable to the concentration of the isotopes in the actual fuel. obtained from the simulation at 300 K was 5. 8 W/mK, which is within the error bar of the value of the defect free UO2, 6.0 W/mK. Thus, there is no obvious influence of the mass defects to the ther mal conductivity. Lattice Dynamics Calculations Although the thermal conductivities of UO2 x have been quantitatively determined by the nonequilibrium MD simulations, further understanding in the mechanism of h eat transfer in this material can be gained by lattice dynamics (L D) calculations. LD is a method which allows analysis of vibrational modes in solid. The fo rmalism of LD is described in Appendix I. As discussed in Chapter 2, we fo llowed the approach of Allen et al.85 and performed LD calculations on 666 unit cells of UO2 x for this analysis. The cal culations are performed only at the zone center ( k =0) because of the large repeat unit. In principle, using a larger structure is better since the resolution of frequency is proportional to the nu mber of atoms in the repeat unit, N However, increasing the number of atoms increas es the size of the dynamical matrix which is proportional to N2, and diagonalization of a large matr ix is very demanding. Thus the computational load rapidly increases with the nu mber of atoms in the system. For example, using a 2.2 GHz 64bit AMD Opte ron processor with 4 GB of RAM, a single LD calculation on the 6x6x6 unitcell system takes approximately 48 hours to complete. As the first step in characterizing the the structures, phonon densities of states (DOS), g ( ), of UO1.98UO2.25 were calculated. Figures 615 (a) (f) are the phonon DOS as a function of concentration. In the defect free structure, UO2.00, the vibrational freque ncies are well defined resulting in very sharp peaks over the entire frequency range. The Debye model like behavior PAGE 111 111 (g( ) ~ 2) continues up to ~4 THz. As the defect c oncentration is increased either in the hyperor hypostoichiometric directions, the DOS smear out and the minor peaks disappear. There is no significant change in the DOS for 0.020 x 0.020 indicating that th e concentrations of defects are indeed low enough that they do not in fluence the dynamics of the vibrational modes in UO2. As discussed in Chapter 2, ma terials containing def ects can have vibrational modes that are spatially localized. These modes are called loco ns and they cannot carry thermal energy through the system;85 therefore they can act as a source of thermal resistance. The degree of localization of each mode can be quantified by the participation ratio, p, as given in Equation I5. 2 1 i i ie e N p In an ideal crystalline structure, p is of order of unity for some atoms and small for others, indicating that specific atoms are associated wit any given mode ; this is equivalent to the mode having a welldefined wave vector. 158 If the localization occurs, p approaches ~1/ N Figure 619 is the participati on of 666 single crystal UO2 x. The participation ratio of the defect free UO2 is shown in Figure 619 (c). p is widely scattered around 0.5, which is a characteristic of a crystalline soli d. Once the defects are introduced, p displays much narrower range of distribution. This means that the nearby vibrations in fr equency space in the UO2 x are coupled to each other, and should thus contribute approximately equa lly to the thermal transport. Even for the lowest concentrations ( x =0.010) of defects, this is still true. The propagon and diffuson modes can be di stinguished by using the polarization of vibrational modes in UO2 x. (Chapter 2) Figure 620 is a co mparison of the polarization vectors at 1.5 THz in UO2 x. Since the polarization vectors lies on the sphere of unit radius, two projection methods, i.e. direct and equalarea pr ojection, are used for comparison. The direct PAGE 112 112 projection is the simple projecti on of the vector components in th e Cartesian coordinates. The equalarea projection is achieved by scaling the two coordinates acco rding to their polar angle. For instance, if the equalarea projection is a bout the zaxis, x and y co mponents of the vectors are scaled by 2sin( /2), where is the angle between zaxis and the vector. The frequency of 1.5 THz was chosen since it is at the low end of the frequency ra nge which is expected to be dominated by the propagon modes. The pol arization vectors in the defect free UO2 are well defined and take the value of 1. For the low concentrations of defects (0.020 x 0.020), there are areas which data points are concentrated on both projection me thods. This indicating that there are still a fair number of modes with well defined polarization vectors. At x =0.125, the distribution of the data points is uniform and there is no well defined polarization vector. The story is different at 10 THz (Figure 621). Th e polarization vectors are homogeneous and isotropic at all the concentrations of defects in both projections. There is no obvious feature to distinguish the materials from the polarization plots. Thus, there is no well defined polarization vector, which is the feature of diffuson modes. Discussion Our analysis showed that the UO2+ x displays a chemical contraction compared to pure UO2. Further analysis showed that this contracti on is due to the electrost atic attraction by the OI and U5+. A similar, but opposite, effect takes place in UO2x. The thermal expansion of the UO2+ x decreased with increasing concentration of defects, while over the small range of stoichiometry considered, UO2x structures showed no such change. The reason for this difference is due to the range of concentrations considered. The small concentration of oxygen vacancies in UO2x also resulted in negligible change in the thermal expansion. If only the range of x 0.025 is considered, UO2+ x also show no discernible change in the thermal expansion. PAGE 113 113 The concentration dependence of the thermal conductivities of UO2 x was determined by the MD simulations using the direct me thod. At low con centrations, both UO2+ x and UO2x show approximate linear dependence. These results are consistent with the KlemensCallaway (KC) theory of phonondefect scattering which suggests th at the scattering rate to be proportional to the defect concentration.74, 159 The concentration dependence in the low concentration limit is rather similar in UO2+ x and UO2x despite the different defects in the structures. The ratios of the slopes of the concentration dependence are 2.3 and 2.2 when the absolute temperature is doubled from 800 K to 1600 K. This suggests that at these low concentra tions, the defect free bulk thermal conductivity temperature dependence of 1/ T is still valid; phononphonon scattering still dominates over the phonondefect sca ttering. At higher concentrations ( x >0.025), however, there is no useful model available and a mo re sophisticated theory is needed. Further analysis was performed using the LD calculations. The resulting picture was consistent with the AllenFeldman theory of amorphous materials, and work on YSZ by Schelling et al.160. From the particip ation ratio analysis, most of the vibrational modes are found to be delocalized. The only case where th ere was some localization was observed was x =0.125 at around 2728 THz. However, the fraction of modes localized is sti ll a small portion of the entire vibrational modes. Pola rization analysis did not show a clear transition frequency between propagon and diffuson modes, but it is determ ined to be in between 1.5 and 10 THz. This study showed that the UO2 with anion vacancies and in terstitials do not hinder the heat transfer by the localization of vibrationa l modes. Although the thermal conductivity is reduced by ~60 % for high concentration of oxyge n interstitials, the reduction comes from the low diffusivity of the diffuson modes. At low concentrations of  x <0.020, decreases with PAGE 114 114 increasing temperature just like a single crystalline bulk UO2; therefore, the phonon description appears to remain valid. As a final remark, it is worthwhile noting the e ffect of extrinsic point defects. There are several forms of fission products present in nuclea r fuel material: solid solutions, dispersions and fission gases. Fission products forming solid solutions lower the th ermal conductivity of UO2 by becoming mass defects or straini ng the host crystal, thereby ac ting as scattering centers for phonons.36, 145 The effect of formation of a second pha se consisting of fiss ion products is known to either increase or decreas e the thermal conductivity depend ing on the type of the phase.36, 146 Fission products are often metallic precipitates and their beha vior depends on the effect of burn up.40 Near the center line region, the relatively high temperature allows formation of 0.05 1 m metallic particles, and thei r electrical conductivity impr oves the heat transfer in UO2.146, 149 In a high burnup structure, the precipitates are of ten quite tiny (order of nanometers) and they contribute to phonon scattering and lower the thermal conductivity of the fuel. The effect of porosity has also been studied quite extensively. Again, there are diffe rent effects depending on the size of the pore. Microbubbl es often found in the crystal la ttice act as obstacles to the phonons. Larger porosities are commonly found at GBs and impede heat transfer. A very nice review on the topic of the effect of fission products on thermal trans port is included in the work by Lucuta et al.36. PAGE 115 115 Table 61. Experimental and DFT va lues of the formation energies. EExp (eV) EDFT (eV) OI 2.92.5131, 132 VO 6.16.7131, 132 UI 7.07.3131, 132 VU 3.34.8131, 132 O Frenkel 3.15.422, 140, 161165 3.63.9131, 132, 166, 167 U Frenkel 8.59.6140 9.512.6131, 132 Schottky 67140 4.96131, 132 PAGE 116 116 Table 62. Potential parameters for the shortrange interactions in UO2 x. O2O2U4+O2U5+O2U3+O2Aij (eV) 9547.96 1761.775 2386.42 1165.65 ij () 0.2192 0.35643 0.3411 0.3786 Cij (eV6) 32 0 0 0 PAGE 117 117 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 0.60.81.01.21.4 O E/EExpU E/EDFT Busker1 Basak Karakasidis Lewis a Morelon Tharmalingam 1 Walker Yamada Busker2 Catlow 1 Catlow 2 Grimes1 Jackson 1 Jackson 2 Lewis b Lewis c Meis 1 Meis 2 Sindzingre Reference Figure 61. Normalized Frenkel and antiF renkel pair formation energies by empirical potentials. Filled symbols are the rigid i on models and open symbols are for the shell models. Uranium FP and oxygen antiFP form ation energies are normalized by the DFT and experimental values, respectively. All data except Busker1 are taken from the work of Govers.115 PAGE 118 118 Figure 62. The 2:2:2 defect cluster in UO2. Blue and red are uranium and oxygen ions; only the interstitial oxygens are shown. Faded re d is used to indicate the oxygen vacancies. The lines are used to clarify the re lative positions of the atoms. PAGE 119 119 0 5 10 15 20 25 30 0246810 r ()pdf Before annealing After annealing (U(4+)O) After annealing (U(5+)O) Figure 63. Pair distributi on function between UO in UO2.25. Blue curve is the UO PDF for both U4+O2and U5+O2. There is no difference before relaxation. The PDF for U4+O2and U5+O2after annealing are shown by red a nd green curves, respectively. PAGE 120 120 0 20 40 60 80 100 120 0246810 r ()pdf Before annealing After annealing (U(4+)O) After annealing (U(3+)O) Figure 64. Pair distributi on function between UO in UO1.80. Blue curve is the UO PDFfor both U4+O2and U3+O2. There is no difference before relaxation. The PDF for U4+O2and U3+O2after annealing are shown by red a nd green curves, respectively. PAGE 121 121 0.994 0.995 0.996 0.997 0.998 0.999 1.000 1.001 0.000.050.100.150.200.250.30 xa/a0 Figure 65. Normalized lattice parameter of UO2+x as a function of defect concentration. The lattice parameters are given at 0 K after the high temperature annealing. a0 is the lattice parameter of UO2.000. PAGE 122 122 1.0000 1.0005 1.0010 1.0015 1.0020 1.0025 1.0030 0.0000.0050.0100.0150.0200.025 x a(T)/a0 Figure 66. Normalized lattice parameter of UO2x as a function of defect concentration. The lattice parameters are given at 0 K after the high temperature annealing. a0 is the lattice parameter of UO2.000. PAGE 123 123 1.000 1.002 1.004 1.006 1.008 1.010 1.012 0500100015002000 T (K)a(T)/a0 x=0.000 x=0.010 x=0.025 x=0.075 x=0.125 x=0.175 x=0.200 x=0.250 Figure 67. Thermal expansion of UO2+x. Lattice parameters are normalized by the 0 K lattice parameter obtained from the chemical expansion. PAGE 124 124 0 1 2 3 4 5 6 7 8 0.000.050.100.150.200.250.30 x [106 1/K] Figure 68. Thermal expa nsion coefficient of UO2+x as a function of x is obtained by a linear fit to the data in Figure 67. PAGE 125 125 1.000 1.002 1.004 1.006 1.008 1.010 1.012 0500100015002000 T (K)a(T)/a0 x=0.005 x=0.010 x=0.015 x=0.020 Figure 69. Thermal expansion of UO2x. Lattice parameters are normalized by the 0 K lattice parameter obtained from the chemical expansion. PAGE 126 126 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.000.010.020.030.040.050.06 1/Lz [1/nm]1/ [mK/W] x=0.000 x=0.010 x=0.025 x=0.075 x=0.125 x=0.175 x=0.200 x=0.250 Figure 610. Inverse thermal conductivity of UO2+x with various concentrations of oxygen interstitials as a f unction of simulation cell length at 800 K. PAGE 127 127 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.000.010.020.030.040.050.06 1/Lz (1/nm)1/ (mK/W) x=0.000 x=0.010 x=0.025 x=0.075 x=0.125 x=0.175 x=0.200 x=0.250 Figure 611. Inverse thermal conductivity of UO2+x with various concentrations of oxygen interstitials as a f unction of simulation cell length at 1600 K. PAGE 128 128 0 2 4 6 8 10 12 14 16 0.000.050.100.150.200.250.30 x (W/mK) MD (800 K) MD (1600 K) Lucuta (773 K) Lucuta (1673 K) Figure 612. Bulk thermal conductivity of UO2+x as a function of degr ee of offstoichiometry without the anharmonic corre ction at 800 and 1600 K. Th e experimental data are obtained from Lucuta et al. 149. PAGE 129 129 0 1 2 3 4 5 6 7 0.000.050.100.150.200.250.30 x (W/mK) MD (800 K) MD (1600 K) Lucuta (773 K) Lucuta (1673 K) Figure 613. Bulk thermal conductivity of UO2+x as a function of degr ee of offstoichiometry with anharmonic correction at 800 and 1600 K. The thermal conductivity from MD simulations were rescaled by the anharmonicity correction. PAGE 130 130 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 0.000.010.020.030.040.050.06 1/Lz (1/nm)1/ (mK/W) x=0.000 x=0.005 x=0.010 x=0.015 x=0.020 Figure 614. Inverse thermal conductivity of UO2x with various concentrations of oxygen vacancies as a function of simu lation cell length at 800 K. PAGE 131 131 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.000.010.020.030.040.050.06 1/Lz (1/nm)1/ (mK/W) x=0.000 x=0.005 x=0.010 x=0.015 x=0.020 Figure 615. Inverse thermal conductivity of UO2x with various concentrations of oxygen vacancies as a function of simu lation cell length at 1600 K. PAGE 132 132 0 1 2 3 4 5 0.0000.0050.0100.0150.0200.025 x (W/mK) 800 K 1600 K Figure 616. Bulk thermal conductivity of UO2x as a function of degree of offstoichiometry at 800 K and 1600 K. PAGE 133 133 0.0 1.0 2.0 3.0 4.0 5.0 6.0 0.0000.0050.0100.0150.0200.0250.030 x (W/mK) UO2x UO2+x 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 0.0000.0050.0100.0150.0200.0250.030 x (W/mK) UO2x UO2+x Figure 617. Thermal conductivity of UO2x as a function of degree of offstoichiometry at 800 K (a) and 1600 K (b) for x 0.025. a b PAGE 134 134 0.00 0.01 0.02 0.03 0.04 0.05 051015202530 f (THz)DOS (Arb. Unit) 0.00 0.01 0.02 0.03 0.04 0.05 051015202530 f (THz)DOS (Arb. Unit) 0.00 0.01 0.02 0.03 0.04 0.05 051015202530 f (THz)DOS (Arb. Unit) 0.00 0.01 0.02 0.03 0.04 0.05 051015202530 f (THz)DOS (Arb. Unit) 0.00 0.01 0.02 0.03 0.04 0.05 051015202530 f (THz)DOS (Arb. Unit) 0.00 0.01 0.02 0.03 0.04 0.05 051015202530 f (THz)DOS (Arb. Unit) Figure 618. Phonon density of states of UO2 x for 0.020 x 0.125. System size is 666 unit cells. The frequency bin size is 0.2 THz. ( a ) x =0.020 (b) x =0.010 (c) x =0.000 (d) x =0.010 (e) x =+0.020 (f) x =+0.125 PAGE 135 135 Figure 619. Participation ratios of UO2 x for 0.020 x 0.125. System size is 666 unit cells. The frequency bin size is 0.2 THz. ( a ) x =0.020 (b) x =0.010 (c) x =0.000 (d) x =+0.010 (e) x =+0.020 (f) x =+0.125 PAGE 136 136 Figure 620. Projections of the polarization vectors of vibrational modes in UO2 x at 1.5 THz. Figures correspond to the de fect concentrations of x =0.020 (a,b), 0.010 (c,d), 0.010 (e,f), 0.020 (g,h), and 0.125(i,j). Left fi gures are direct proj ection, and the right figures are the equalarea projections. a b c d e f g h i j PAGE 137 137 Figure 621. Projections of the polarization vectors of vibrational modes in UO2 x at 10 THz. Figures correspond to the concentrations of x =0.020 (a,b), 0.020 (c,d), and 0.125(e,f). Left figures are direct projection, and the ri ght figures are the equalarea projections. c d a b e f PAGE 138 138 CHAPTER 7 METHODOLOGY FOR THE SIMULATION OF RADIATION DAMAGE IN UO2 Introduction In nuclear fuel materials such as UO2, it is extremely important to consider the effects of radiation when assessing the performance in bo th operating and storage conditions. However, the complexity of the processes occurring under irradiation makes it very difficult to decouple each phenomenon and to understand the mechanisms of certain processes. For instance, high energy radiation causes the atoms to be knocked o ff their lattice sites. Because of the high energy deposited by the radiation in to the solid, extreme local temperature and pressure gradients can be created, with atoms diffusing accordingly. The local thermal transport properties will also be altered. At the same time, radiations can also cause electronic exc itations and can initiate chemical reactions. All these processes are enta ngled and the complexity is overwhelming. For this reason, simulations are a valuable tool to ch aracterize the details of such processes, albeit with certain limitations depending on the methodology. MD simulation is one of the tools which allow investigation of highly nonequilibrium phe nomenon such as radiation damage. However, the conventional MD algorithm does not work in such extreme conditions and some modifications are required. This chapter, ther efore, establishes a methodology for the simulation of the radiation damage studies in UO2 and presents the implemen tation and validation of a radiation damage simulation code. At least one of these methods, the form of the variable time step algorithm is completely new. While the ot her methods have been used previously, their applicability and their limitations have not been systematically established. PAGE 139 139 Background Type of Radiation Naturally occurring radioactive materials su ch as uranium, radium, and radon, are everywhere on the earth. Radiation also come s from space and constantly penetrates or is absorbed by the atmosphere. The natural backgr ound radiation is roughly 100 times greater than any of the human caused background radiations.61 Manmade radiation sources are used for a number of applications includi ng electrical power generation, medical imaging, radiation therapy, food processing, and weapons.61 All the materials used in these applications must be designed to take account of the e ffects of radiation. In orde r to understand the effects of radiation, one must understand the different types of radiations th at exist in nature. There are four types of radiations typical ly considered in nuclear appli cations: neutrons, alpha, beta, and gamma radiations. Neutron radiation consists of a ray of energe tic neutrons. The neut ron is one of the two nucleons, the other being the proton, which are the building blocks of the nuclei of atoms. It has no electric charge, and has mass a mass of 938.27 MeV/c2 (1.67261027 kg).168 Because of the absence of electric charge, neutrons often tr avel significant distan ce depending on the energy without interacting w ith other atoms. Alpha radiation consists of the nucleus of a helium atom, i.e., two neutrons and two protons. It has charge of +2 and a mass of 3.7274 GeV/c2 (6.64471027 kg).61 It is produced by radioactive decay processes and changes the at omic number of the atoms involved. Alpha radiation does not propagate a l ong distance because of the relative ly large charge and mass. Beta radiation is electrons or positrons released during certai n types of radioactive decays, called beta decays. The decay processes involve ionization of the decay products because of the loss or gain of electrons. PAGE 140 140 Gamma ray is the radiation of high ener gy electromagnetic waves (photons) produced during radioactive decay. In contrast to the ot her types of radiation, this radiation involves transition of the energy levels of atoms. Quite often atoms produced through alpha or beta decays are in excited states. When the excited at oms transition from their energy levels to lower energy states, the excess energy is released in the form of ener getic photons. A simple example is 60Co decaying to 60Ni through beta decay with th e half life of 5.27 years.169 The gamma radiation from this decay process is commonly us ed in medical and industrial applications. All these different types of radiation occur in nuclear materials and can change the materials properties. Fission and Radiation Nuclear power reactors for el ectrical power generation today involve fission reactions. Fission reactions are created by neut ron radiation to fissile elemen ts. When a neutron collides with the nucleus of a fissile element, the nucleus can split into nuclei of lighter elements and produces more neutrons. A good example is the fission of 235U. 235U is the only naturally occurring stable fissile element on earth. This is the reason that the uranium is used for nuclear fuel. A representative fission reaction of 235U is, 235U + n 92Kr + 141Ba+ 2n + 200 MeV. (71) In this reaction, with one incident neutron, 2 neutrons are produced. There are a number of fission reactions which result in producing more th an one neutron from a single incident neutron. Consequently, the number of neutrons can increase like an avalanche or chain reaction to sustain the reactions. The energy released during the fiss ion reactions dissipates as heat and can be used to produce electrical power. Ma ny of the fission products (FP) are unstable and spontaneously PAGE 141 141 decay into lighter elements including radioactiv e ones. Subsequent decays of these FPs produce n, and radiations. Fission Products and the Type of Stopping Before going into the specifics of the radiati on effects in UO2, the general concept of energy loss process of energetic particles in solid must be introduced. Energetic particles in a solid can lose their energies by either collision wi th nuclei or electrons. A projectile ion (often referred to as primary knockon atom, or PKA) in a solid can collide with th e ions in the solid. The energy transfer between the ions depends on the mass, charge, and momentum of both PKA and the target ion. While the PKA passes through the solid, the matrix atoms recoil and absorb energy, and the PKA is deflected. This process is described by the nuclearnuclear collisions and the collision is elastic. In addition to collid ing with the nuclei, the atoms can collide with electrons. The nature of electr onic energy loss is rather comple x and the following factors must be considered: Direct energy transfer by electronelectron collisions Excitation or ionization of target at oms (promotion of localized electrons) Excitation of transition in the band structur e (promotion of weakly localized or free electrons) Excitation, ionization, a nd electron capture of th e projectile ion itself There are several wellknown el ectronic energy loss models with different assumptions on the charge distribution,170172 however, the details of these mode ls are beyond the scope of this work. PAGE 142 142 Previous Experimental Studies on Radiation Damage in UO2 The effect of the radiation is determined by its type. radiation causes ionization of the surrounding matrix and is responsible for chemi cal reactions, bond rupture, valence changes, formation of oxygen bubbles, etc. Their effect is mainly chemical in nature and produces very little structural damage. Regard ing damage to the structure, decays are most detrimental in the fissile fuel matrix. A typical decay process produces a high energy particle of ~5 MeV, with the daughter atom receiving a recoil energy of ~0.1 MeV 173, 174. An particle loses its energy mostly by electronic stopping, thus causing ionization; it loses onl y a small portion of the energy for atomic displacement. By contrast, heavy recoil atoms from decays lose almost all of their energy by nuclear stopping and gene rate roughly 90 % of the lattice damage in the fuel matrix 173, 174. Because of the differences in the energies and masses, particles and the recoil atoms have different range s of propagation. particles often travels from 20 to 30 m, but the recoil atoms travel only about 25 nm.173, 174 This implies that the damage from the recoil atoms is spatially localized, but the damage from particles is dispersed. The effects of such structural damage on nuclear reactor fuel significantly a ffect its physical and ch emical properties. However, the details of such effects have never be en characterized in detail. Furthermore, there is no satisfactory model for the defect production and annihilation mechanisms. Neutrons also occasionally collide with the lat tice atoms, but they are not as effective as charged particles especially at the energy range t ypical in typical fission reactors. Therefore the direct collision of neutrons to nonfissile elements is rare and ha rdly damages the structure. The only way neutrons can contribu te to structural change is through the fission reactions 173. In fact, the energy dissipated by the slowi ng down of the fission products is the major source of thermal energy for the electrical power genera tion in the nuclear fission reactors.11 PAGE 143 143 In UO2, the high energy fission products lose most of their energy by el ectronic stopping. Therefore, a significant thermal effect is expected because UO2 is an electrical insulator; the energy deposited by the electronic st opping cannot dissipate like that of free electrons in metals. Roughly 30 % of the energy is lost in instantaneous Xray radiat ion and elastic waves, and the rest is transformed into heat which is referred to as the thermal spike.173, 174 The thermal spike has several effects, including generating a signif icant local temperature and pressure gradient, causing surrounding atoms to diffuse. This extrem e temperature and pressure gradient is known to cause the formation of the U4O9 phase in UO2.173, 174 The thermal spike also allows the recovery of defects formed by the impact of hi ghly energetic atoms. The recovery process is more significant in UO2 than in UN or UC, both of which ar e electrically conduc tive materials. Figure 71 is an Arrhenius plot of the cation diffusi on coefficient in (U,Pu)O2, (U,Pu)N, and (U,Pu)C. In the all three materials, the curve is strongly temperature de pendent for temperatures above ~1100 C indicating that the diffu sion process is thermally activated.173 However, below 1100 C, the selfdiffusion coefficien ts are independent of temperat ure and much higher than that of thermal diffusion. Detail of how the diffusion of cations takes place in this regime remains unclear. It has been postulated that because of the extremely high temperature and pressure of the fission spike uranium atom diffuse thr ough the core of the displacement cascade.173 It should be noted that the order of diffusion coefficients is D(UO2)>D(UN)>D(UC), which is opposite to the order of their thermal conductivities. Some experimental studies have been perfor med on simulated fuel and ion implantations rather than the actual fuel samples.175178 These studies show a consiste nt picture of the effect of radiation damage on the stability of the crystal structure, and in produci ng significant changes in the lattice parameter. Recently Ronchi et al .178 reported the effect of fission fragment on UO2 PAGE 144 144 single crystal film. Despite the extremely hi gh energy (~80 MeV) of the fission products, no significant structural change was observed; however, the U4O9 phase was present in the sample after the irradiat ion. Ronchi et al .178 attributed the very small st ructural damage to the shock wave produced in UO2, which is unloaded on the free surfaces. The above experimental studies show the unresolved issues in the radiation damage of UO2. It is not clear how the athermal diffusion of U atoms occur in the radiation environment. How does the recovery process take place to heal the crystal at such high energy irradiation? What is the role of oxygen ions in such processe s? Are these processes affected by the presence of defects such as point defect s, dislocations, and GBs? Previous Simulation Studies on Radiation Damage in UO2 A number of simulations have been performed to investigate radiation damage in solids, especially metals and oxide materials. Simu lations of radiation damage were pioneered by Gibson et al. .179 They studied the effects of PKA ener gy and direction on singl e crystal copper. Later extensive development on the simu lation method was done by Harrison and his collaborators180183 for the application of s puttering studies on Si and metal surfaces. During the time of Gibson and Harrison, the availability of the computational power was limited and their studies were restricted to small systems (<1000 atoms) with simple interatomic potentials. Consequently, the energy of PKA had to be low (~ a few 10s of eV). As more computational power became available, the area flourished. To day atomistic simulations of radiation damage are routinely performed on metals, oxides and se miconductors using the stateofart interatomic potentials.184189 The first MD simulation of radiation damage in UO2 was done by Van Brutzel et al .190 They developed their own UO2 potential based on the defect formation energies of U and O defects.191 The potential was a Bucki ngham type rigidion potential consisting of three radial PAGE 145 145 sections. They performed the radiation damage simulation in single crystalline UO2 from a few keV190 up to 80 keV192. They investigated the defect st ructure and the number of displaced atoms. There was no correlation between the numb er of defects produced and the directions of PKA. The number of Frenkel pairs produced as a function of PKA energy showed that the exponent is ~0.9 rather than the Norg ettRobinsonTorrens linear dependence.193 They have also performed simulations to investig ate the effect of cascade overla p. They found that the number of Frenkel pairs formed saturated after 7th PKA irradiation. The PKA dose rate of 0.79 keV/ps perhaps is unrealistically high since the defects recombination may take place on a much longer time scale. Furthermore their defect identi fication method is diffe rent in different publications,190, 192 and a number of technical issues associated with the thermostat and variable time step algorithm are not addressed. Implementation of the Molecula r Dynamics Simulation Code As discussed in the previous se ction, there a number of issues to be resolved in the effects of radiation damage on UO2. Although it is impossible to capt ure all the physical and chemical properties well by any simulations, simulation approaches have the advantage of having complete control over the material. In particular, MD simulations allow us to follow the highly nonlinear dynamic process of th e defect production and annih ilation by the radiation damage which occurs in very short time (order of pico seconds) with atomic reso lution. Since standard MD simulation codes are not capable of simulatin g extreme radiation cond ition, some significant modifications are needed. The algorithm a nd implementation of such modifications are described in the following. Universal ZieglerBiersakLittmark Potential During high energy collisions of atoms whic h occur in a radiation environment, two atoms can get extremely close (typically ~1 or le ss) compared to the normal distance given by PAGE 146 146 the thermodynamic equilibrium. Standard interato mic potentials are not designed to adequately describe the interatomic interactions at such short distances. For example, the commonly used Buckingham type potential has the functi onal form given in Equation 38: 6) / exp( ) (ij ij ij ij ij ijr c r A r V The second term, which represents the van der Waals interaction, goes to negative infinity when r goes to zero although the first exponential term remains finite. This means that the atoms collapse to a singularity if they are brought very close to each other. This is illustrated in Figure 72, which shows the plot of the Bu ckingham potential for OO interaction with the Yamada model. It can be seen in Figure 72 th at without the ZBL potential, the height of the potential energy barrier of the Buckingham potential is only about 260 eV or so, and a rather low energy PKA can overcome this barrier. This low ba rrier is an artifact of the functional form of the interatomic potential since we know that the repulsion of th e electron clouds should prevent this. Another example is the electrostatic inte raction: the attractive electrostatic interaction between two oppositely charged atoms al so goes to negative infinity as r approaches 0. In typical MD simulations, however, these singularities does not cause an problems due to small atomic displacements; however, once the force on an atom becomes very large, the atom will jump unrealistically large distance in a single step. (An exception is when the time step size is allowed to vary according to the force on the atom. If the time step is adjusted according to the force, the two atoms will keep approaching to the singularity indefinitely. The method of variable time step is discussed below.) Even if such catastr ophic failure does not happen, it is possible for the potential to fail in a rather subtle manner by giving unrealistic potential energies. Empirical potentials in general are fit to the bulk propert ies of the materials to reproduce the static and dynamic behavior of the materials near equilibrium. These potentia ls are therefore not meant to PAGE 147 147 reproduce the realistic potential energy surface in extreme conditions such as high energy collisions. A more realistic model of the very shortranged interacti on was given by Ziegler, Biersak and Littmark (ZBL).194 They found from two body collision experiments of a large number of elements that the interatomic interaction can be simplified in a universal form, which may be applied to any element in the periodic tabl e. This potential, referred to as the Universal ZBL potential, is based on the mode ls of ThomasFermi and Linhard.195 The ZBL potential is in fact a screened electrostatic potential of the nucleonnucleon interaction: ) ( ) (2 1r r e Z Z r V. (71) Z1 and Z2 are the atomic numbers of the atoms in collision, e is the electron charge, and r is the interatomic distance. ( r ) is the screening function, 4 1) / exp( ) (i u i ia r b A r. (72) Ai and bi are predefined parameters, whose values are given in Table 72. au is the length scale given by, Bohr ua Z Z a23 0 2 23 0 18854 0 (73) And aBohr is the Bohr radius, 0.539177 Details of the derivation of this potential is given in the literature.194 Since the ZBL potential is intended to opera te only at very short distances and the Buckingham potential works well for nearequilibrium separations it is importa nt to know how to connect these two potentials Unfortunately, there is no uni que way to define where the transition between the two potentials should occur, although in reality there is a single potential energy surface for a pair of atoms. The standa rd approach is to connect the Buckingham and ZBL potentials by a spline of the ex ponential of a polynomial of the form, PAGE 148 148 5 5 4 4 3 3 2 2 1 0exp ) ( r f r f r f r f r f f r V (74) The fi are constants determined by ma tching the V(r), dV(r)/dr, and d2V(r)/dr2 at the connection points to the ZBL and Buckingham potential. However, the choi ce of the spline region is arbitrary and depends very much on the atoms to be modeled. In principle, there is no need to use all the 6 parameters since the continuity of the 2nd derivatives are not required. However, in ionic solids, the gap between the ZBL and Buck ingham can be quite large, unlike metals, and having the extra parameters can be he lpful. (See the work by Fikar et al.196 for the spline with 4 parameters for metals.) The va lues of the spline parameters a nd connection points that we have chosen for are given in Table 71. How all the potentials are combined to form th e MD potential is illu strated in Figure 72 for the case of OO interaction. In order to show how the ZBL pot ential works, a simple simulation of an OO collision is performed. In th is demonstration, an oxygen atom is fixed at center of the simulation cell. Another O is simp ly placed 0.1 away from the center where the ZBL potential is active. Figure 73 shows the comparison of the tr ajectory of the O atom. With the ZBL potential, the atom simply is simply repelle d by the potential of th e O atom at the origin and rolls down the potential energy surface. By contrast, without the ZB L potential, the O atom is trapped in the potential well and violently rattles. In this simulation, the time step was fixed at 0.005 fs, but the interaction is so violent and th e trajectory cannot be a ccurately calculated and eventually the atom escapes from the well. If the time step is small enough, the O atom will approach the origin and will fall into the si ngularity. Thus, the ZBL potential must be implemented in the radiation damage simulation. PAGE 149 149 Variable Time Step Algorithm Even when the ZBL potential is used, the inte ratomic forces can still be extremely large when two atoms are very close to each other. This still poses a major problem if a normal MD time step (~1 fs) is used. Such a large force on an atom can cause the atom to move a very large distance in a single step, resulting in very inaccurate trajectory. One solution to this problem is to use an extremely short time step size. Howeve r, for such a case, the simulation will take an unreasonably long time to complete. A better approach is to allow the time step to change during the course of the simulation. When an atom is highly energetic, the time step is automatically reduced. Once the energy is dissipated, the time step size is allowed to incr ease. This will allow accurate atom trajectories to be determined, while maximizing the time simulated. The common algorithm to do this in radiation damage simulation is to set a restriction on the maximum distance any atom can travel in a single MD tim e step from which the necessary time step size can be calculated 197. The typical limit on the maximum distance for an atom to travel in a single step is <1 Although this algo rithm has been extensively used,190, 192, 195, 198200 the accuracy of the atomic trajectories has not been clearly established. We have therefore developed a new variable time step algorithm specifically designed to ensure that the trajectory of th e energetic atoms is accurate. The original idea is taken from the orbit calculation of an sa tellite under a gravitational field de scribed in the astrophysics literature 201. Orbit calculations are know n to require rather high accura cy; otherwise, the energy conservation deteriorates rapidly and the trajec tory will be incorrect. This is because the differential equation to be solved is of a t ype called stiff problem in mathematics whose solution indicates rapid change in certain regions of the pa rameter space. The high energy collision of atoms also falls into this stiff pr oblem category. The orig inal algorithm of the adaptive time step method was applied to a pur ely two body interaction and the RungeKutta PAGE 150 150 method was applied to integr ate the equation of motion.201 In our inhouse MD code, the scheme is modified to work with the 5th order Gear predictorcorrector integrator which is already incorporated for solving th e equation of motion. The idea is to calculate the trajectory of the most energetic atom in the system in two ways: one consisting of a single time step the other consisting of two half time steps; see Figure 74 for the illustration of this process. The two halfstep trajectory is in principle more accurate than the one large step. The two trajectories therefore sh ould agree within the accuracy of the integration method used. If the discrepancy in the trajectory is smaller than the intrinsic accuracy of the integration method or of a user defined accuracy, th en the time step used is adequate. However, if the discrepancy in the trajectory is too large, then the next trial time step size is determined based on the size of the discrepancy between th e two trajectories. The whole process repeats until an acceptable time step size is determined. By performing this analysis for the most energetic (i.e., fastest moving) atom in the syst em, accurate trajectories for all other atoms are guaranteed. Moreover, since only the trajectory of the PKA is used in the analysis, the method is both fast and natura lly parallelizable. The implementation of this algorithm naturally follows the process of the Gear predictorcorrector integration scheme. First the most ener getic atom in the system is identified by the predictor step, which is nothing but a calculation of the trajectory using the Taylor expansion: ... ) ( 2 1 ) ( ) ( ) (2 dt t a dt t v t x dt t x (75) where x ( t ), v ( t ), a ( t ) are the position, velocity, accel eration of the atom at time t This step is very fast, and the most energetic atom will give the largest displacement. Then the force calculation is performed using th e predicted positions of the atoms as in standard predictorcorrector step, but the calculation is done only on the most energe tic atom. Then the trajectory PAGE 151 151 of the energetic atom is corrected by the correc tor step. Once the new position of the atom is determined using time step, dt the code proceeds to calculate the same trajectory by two half steps. The entire process of the determination of the time step is shown in Figure 75 as a flow chart. The most expensive part of this scheme is the first half of the two small steps since this step requires the force calculation on all the atoms. In order to redu ce the calculation, this step is done only on the near neighbor atoms. When the computation is done using multiple processors, this part of the calculation can be done on a single processor. Another key step is the estimation of the new tim e step size. First the estimate of the time step is given by 6 / 1 0 x x dt dtcurrent est (76) where x0 and x are the user defined error of the trajec tory and the calculat ed difference in the trajectories from the large a nd half steps calculations. x0/ x defines the error ratio of the trajectory. dtest and dtcurrent are the new estimate and current time step size. The power 1/6 comes from the local truncation error being proportional to dt6 for the 5th order Gear predictorcorrector algorithm. The dtest is merely an estimate and it is possible that this value becomes either extremely large or small. In order to avoi d drastic change in a single iteration, a safety net is also applied: otherwise / if / if1 2 1 2 2 1 2 est current est current current est current newdt S S dt dt S S dt dt S dt S dt S dt. (77) S1 and S2 are the safety factors set to 0.9 and 4.0, but the exact value is not critical. The important point is that S1 should be slightly less than 1.0 and S2 should be somewhat larger than PAGE 152 152 1.0. The more these numbers deviate from 1.0, the more aggressive the estimate of the new time step is. In order to establish the validity of the variable time step scheme, the results of short test runs are compared with a fixed time step simulation with a small time step size. Figure 76 shows the trajectory of the PKA in a 101010 unit cells single crystal UO2. The PKA is 1 keV U in [100]; the statistical ensemble used was mi crocanonical, i.e., consta nt energy and constant volume. Each color in the plot corresponds to a specific error ratio as in dicated in the legend. The trajectory from the fixed time step of dt =0.10 fs simulation is indicated by the solid line. In order to obtain the base case scenario with the fixed time step simulation, a number of fixed time step runs were performed with different time st ep size. The time step size of 0.10 fs was determined to be sufficiently small for the genera tion of accurate trajectories in this case. Figure 76 clearly indicates that the variable time step ca n faithfully reproduce the base case trajectory if the error ratio is sufficiently small ( 104) In this simulation, the PKA traveled along the [100] direction and hit the adjacent U atom, transferri ng most of its momentum while small amount of the energy was dissipated to the rest of the system This process continued during the entire 2 fs simulation. This process is appa rent when the variable time step algorithm is used since the time step size is not uniform. When the collision o ccurred, the time step size decreased and thereby the density of the data points increased. To explore the power of this variable time step algorithm, simulations with 80 keV U PKA were performed using the same simple model sy stem. 80 keV is roughly the upper limit that we can handle in a more realistic simulation set up, because of the curr ent limitations on the computational resources. Figure 77 shows the tr ajectories of the PKA obtained from the runs with three different error ratios and one fi xed time step. Again, a number of fixed time PAGE 153 153 simulations were run prior to this comparison so as to determine the maximum time step size needed for the accurate calculation of the PKA tr ajectory. The result was that for 80 keV U PKA, the time step has to be at most 1104 fs. The variable time step algorithm was able to successfully follow the fixed time step trajectory when the error ratio was 1105. To analyze the efficiency of the variable time step scheme, the time evolutions of the time step size in the same 80 keV simulations were m onitored. (See Figure 78.) With the variable time step, the time step size decreases when collis ions occur. The decrease in the time step depends on the error ratio used in the simulation: the smaller the error ratio, the smaller the time step size at the collisions. The time step size reaches down to 1104 fs, only when the error ratio is 1105. In between the collision events, the time step size can go up 3 orders of magnitude from 104 fs to 101 fs This difference gives the gain in the efficiency of the variable time step scheme. In terms of the overall efficiency, the fixed time step simulation took approximately 208 min to complete 2 fs of the run, but the vari able time step took only about 10 min: more than a magnitude of improvement in efficiency, with no loss in the fidelity of the trajectory. Thermostat The last piece of modification needed in the MD code is in the control of the temperature of the system via the thermostat. During radi ation damage, a large amount of energy is converted from the ballistic energy of the PKA to the thermal vibrations of the atoms in the entire system. In reality, when the PKA collides with the ions in the mate rial, the ballistic energy is distributed over ~1024 atoms, and the energy dissipates at the rate determined by the materials thermal conductivity and the temperature gradient. However, MD simulations are limited by the number of atoms (typically less than a few millions of atoms for ionic materials). For such a small simulation cell this thermal energy can be overwhelming if the heat generated is not appropriately treated, even lead ing to melting of the system. Th erefore to mimic the thermally PAGE 154 154 absorbing effects of a very large system, a thermostat is applied around the region of the simulation cell where all the interesting dynamics takes place (the called active region). Figure 79 shows the configuration of th e MD simulation cell. Since the simulation cell configuration is a hybrid of NVE and NVT, this ensemble was dubbed the pseudoNVE (or pNVE) by Corrales et al.200, and we will follow the same nomenclature. The default thermostat in our inhouse MD c ode is the velocity rescaling thermostat. However, this thermostat is rather primitive and is known to have some problems,98, 202 the most prominent of which is that the statistical nature of the system is not canonical. In this thermostat, since the velocities of the atoms in the thermo stat region are changed simply by rescaling the ratio of the target and current temperature, T T0, the equation of motion is modified. Thus no energy and momentum conservation laws are obeyed It also means that the Liouville theorem no longer holds and the phase space distri bution function of the system is broken.203, 204 Moreover, the ergodic theorem will not hold and phys ical observables cannot be calculated from the ensemble average.203, 204 In short, it means that the dynamics of atoms in the thermostat region are not realistic. There are also other subt le problems as well. In particular, since the dynamics of the thermostat region is not correct the atoms adjacent to the thermostat will not behave in a realistic manner; it is not clear how far this effect ex tends. Another problem is that since the velocity rescaling will take place every time step during the simulation, the rate of change in the temperature may depend on the time step size. Because of these problems, one should be very cautious when performing simulati ons using velocity resc aling thermostat. There are other thermostat scheme s which are known to work better.96, 98 Among the available thermostats, the Berendsen thermostat was chosen for this study.205 The Berendsen thermostat is also a type of velocity rescali ng, but is both physically more reasonable and is PAGE 155 155 better founded in statistical mechanics. The Be rendsen thermostat is derived from the same formalism as the Langevin thermostat 205, 206. The idea of the Langevin thermostat is to couple the system with an imagin ary heat bath in appropriate manner such that the system remains in the canonical ensemble.206 The details of the derivation are given in Berendsens original paper.205 The Berendsen thermostat can be viewed as a first order approximation to the full Langevin approach. The way Berendsen thermostat adjust s the velocity is thr ough a scaling factor, 1 10T T dt (76) where dt is the current time step size, and is the userdefined time constant for thermal equilibration. T0 and T are the target and current temperature. From its form, it is clear that the Berendsen rescaling can be much gentler than the simple velocity rescaling. To illustrate the range of the strength of Berendsen thermost at, consider the two extreme cases. When is infinity, =1 and no rescaling is done; th is corresponds to the micro canonical ensemble. When = dt then it is identical to the simple velocity re scaling. The beauty of the Berendsen thermostat is that since it is a simplifica tion to the Langevin approach, it is approximately canonical and its strength can be controlle d to minimize the disturbance to the system unlike the simple velocity rescaling scheme. Thus the dynamics of the system can be kept close to the realistic phase space distribution. The Berendsen thermostat was implemented in the MD simulation code and several test runs were performed. Figure 710 shows the time evolution of the 101010 unit cells single crystal UO2 using the Berendsen thermostat. The entire simulation cell was set to be in the thermostat region. =400 fs was used, based on the recommendation to use =100400 fs in the original work by Berendsen on simulations of liquid water.205 We see that the thermal PAGE 156 156 equilibrium is clearly reached by 5 ps with very sm all fluctuation at the end. This indicates that the thermostat is correctly implemented. In order to identify th e appropriate value of for the solid UO2, NPT simulations with a wide range of were performed on the thermally equilibr ated single crystal st ructure at 300 K. The total energy of the system was monitored as a function of time (See Figure 718). The two extreme case of NVE and velocity rescaling simulations with dt =0.1 fs are also included for the purpose of comparison. It is a pparent that the velocity rescaling significantly interrupts the system and the the total energy violently fluctuates. With a small value of the Berendsen thermostat closely follows this behavior. As increases, the behavior approaches that of NVE simulation. Since the system is thermally equilibrated to star t with, the total energy should be constant. The variation in the total energy is qua ntified by the standard deviation and it is plotted as a function of in Figure 712. Note that the abscissa is in log scale. The solid lines are used to indicate the trend and are exte nded to the limiting values of the NVE ensemble. The standard deviation in the kinetic ener gy increases with increasing and that of the potential energy decreases as increases. The fluctuations in the ki netic and potential energy are actually correlated and out of phase. Thus these fluctuati ons cancel each other, a nd the the net fluctuation in the total energy is the difference of the fluctuat ions in the kinetic and potential energies. The variation in the total energy smoot hly decreases with increasing and the drop is minimal after ~1000 fs. This trend is very similar to that gi ven in the Berendsens original work on liquid water. Thus 100 500 fs is also applicable for the solid UO2. Defect Identification Method Before going into the analysis of the radiation damage simu lations, the method of defect identifications must be described. The c ode used was developed by our collaborator Dr. Srinivasan G. Srivilliputhur in MST8 at Lo s Alamos National Labor atory. Some minor PAGE 157 157 modifications were made to work for UO2 and to seamlessly integrate with our inhouse MD code output files. The identification method is to find defects in a structure by comparing the atom positions relative to the defect free structure. Figure 713 shows an illustration of this method. In this section, the defect free structure is structure A, and the structure w ith defects is structure B. To identify an interstitial, the analysis code looks through stru cture B and tries to find the corresponding atoms in structure A. Two atoms in the different structures are considered to be matched when they are within radius r0 of each other. If some at oms in structure B do not have any corresponding atom in structure A, then the at om in structure B is an interstitial. To find vacancies, the exact opposite is done; the code goes through structure A and look for neighbors of each atom in structure B. If there is no neighbor for some atoms, they are considered vacancies. Antisite defects can be found in a similar manner. The neighbor radius r0 is the critical parame ter in this method. The results of analysis depend on the value of r0, and this cannot be avoided if this algorithm is used. The dependence of the number of defects on r0 is characterized and discussed later. Radiation Damage Simulations The ZBL potential, variable time step algor ithm, and Berendsen thermostat have been successfully implemented and te sted individually. Now all the methods are combined and applied to radiation damage simulations. Time Step Size The first step is to characterize the base cas e scenario using a fixed time step. The result from this simulation will be the foundation to va lidate the successful implementation of all the modifications together. PAGE 158 158 The radiation damage simulations with 1 ke V U PKA were performe d with varying fixed time step sizes. A single crystal 101010 unit cells UO2 was used, and the simulation was done with NVE ensemble: no thermostat was applied. This small simple setup will provide an extreme case where the maximum time step determin ed here will be sufficient for 1 keV PKA in larger systems with thermostat. Figure 714 shows the vacancies of uranium and oxygen atoms produced during the simulation. The number of both U and O vacancies rapidly increased until 200 fs and then they heat plateau for ~300 fs. The maximum numbers of U and O vacancies are roughly 20 and 80, respectively. Around 500 fs, the number of both vacancies decreased at varying rates. The numbers of U and O vacanci es are exactly identical at each time step for dt 0.1 fs. Thus, 0.1 fs is sufficient for simulations with 1 keV U PKA. The result from a variable time step simulati on with the same system is also shown in Figure 714. The erro r ratio used is 104 as determined from the pr eliminary runs (see Figure 76). The data is taken at the time interval of 10 fs. The result is identical to those of dt <0.1 fs, and this proves that the variable time step algorithm works. After the maximum time step size was determined by the NVE simulations, pNVE runs were performed to show that the same time step can be used in pNVE. The system used was a single crystal 22 2222 unit cells UO2 with 1 keV U PKA. The system size was chosen following van Brutzels work190. A 1 unit cell thick Berendsen thermostat with =400 fs was applied around the simulation cell. The initia l location of the PKA was chosen to be ( x y z )=(6,0,0), based on some trial and error to ensure th at the defect cascade does not interact with the thermostat and boundary of the simulation cell. As mentioned above, the number of defects count ed in the analysis de pends on the value of r0. Figure 715 is the number of defe cts produced in the simulation with dt =0.1 fs as a function PAGE 159 159 of r0. VU, UI, VO, and OI all decrease with incr easing the size of radius r0 since the larger r0, the higher the likelihood of finding neighbor in the corresponding structure. In addition, we see that the number of the vacancies a nd interstitials of the same species do not agree when r0 is greater than roughly 1.2 The value of r0 is typically chosen to be slight ly smaller than the half of the nearest neighbor distance in th e perfect crystalline structure.186, 190 In UO2, the nearest neighbor distance is the distance between U and O atoms in 4a and 8c sites, which is ~2.37 The value of 1 seems to be appropriate. Thus, from hereon, r0=1 is used throughout. Figure 716 shows the result of 4 different runs with dt =0.100, 0.050, 0.001, and 0.0001 fs. It turns out that the number of defects produced was not identical even with dt<0.1 fs although they closely follow each other. This is because the time step size affects the way the thermostat draw energy out of the system. At least it suggests that the dt =0.1 is sufficient of producing the same number of defects as the smaller time step size. Berendsen Time Constant Now the time step size of 0.1 fs is used to ch aracterize the effect of the thermostat. The same single crystal 222222 unit cells UO2 with 1 keV U PKA was used to perform several radiation damage simulations in pNVE ensemble The value of the time constant was changed from 1 fs to 104 fs. The results from an NVE simulati on are also included for comparison. Figure 717 is the plot of how the temperatur e in the active region evolves with varying Initially the temperature starts out at approximately 380 K whic h corresponds to the energy of the 1 keV PKA in the thermally equilibrated sy stem at 300 K. Then the temperature drops rapidly to ~330 K in 370 fs. After that point, the behavior depends on the value of In the NVE ensemble, the temperature remain fairly constant around 330 K. As decreases, the rate of cooling becomes faster. For smaller than 100 fs, the temperature evolutions are almost identical and approach the targ et temperature of 300 K. PAGE 160 160 The change in the total energy is shown in Figur e 718. Curves in this plot indicate how fast the energy is removed by th e thermostat as a function of In the case of NVE, the energy remains constant, as it should. In real UO2, the cooling rate is dictated by the thermal conductivity of the material, and the heat dissipates to the 1024 atoms which are in thermal equilibrium at infinity. The cooling rate cannot go beyond the materials natural ability to transfer heat while keeping the thermostat at 300 K. Figure 718 clearly shows that the cooling rate is very similar for 100 fs. On taking a closer look, fluctu ations in the total energy is clear for 10 fs. Thus =100 fs is a good choice in terms of the fast cooling rate and minimal disturbance to the dynamics of the system. Finally we compare the number of defects in the same set of the simulations. Figure 719 is the result. Reassuringly, the number of defect s produced as a function of time is very similar in all the cases regardless of what value of is used. This means that even if the thermostat is not applied, the number of the defects is not affe cted despite the higher final temperature. The final temperature ranges from 300 K to 330 K. This 30 K difference makes hardly any difference in the defect pr oduction and annihilation. Statistics Since the defect production and annihilation process is a highly nonlinear stochastic process, simulations with slightly different cond itions can result in somewhat different result. The difference in the initial condition can be minute differences in the atom positions and velocities. A good analog of this may be the bil liard balls: a tiny difference in the relative positions and the amount of push you give can give quite different results. A solid crystalline well equilibrated structure is much more ordered than billiard balls and the variation in the results will not be completely chaotic. Nevert heless, it is necessary to characterize how much variation is expected wh en equivalent simulations were performed. PAGE 161 161 In order to characterize the variati on, a single crystal 222222 unit cells UO2 was used to perform 6 independent simulations. The PKA was 1 keV U atom but in each simulation, different U atoms were chosen. Their locations were (6, 0, 0), (0, 6, 0), (0, 0, 6) and in all cases the PKAs were pointing toward the center of the simulation cell. Figure 720 is the plot of the number of U and O vacancies produced during th e 6 simulations. The error bars are used to indicate the standard de viation from all 6 simula tions. The number of VU with the standard deviation is approximately 7 3. The same number for VO is 7020. Therefore, there is a significant variation in the da ta even though the PKA directio ns are crystallographically equivalent. We speculate that this is the effect of small differences in the atomic positions and velocities. Although, the system is thermally equilibrated, every PKA and its surrounding atoms in the 6 simulations have different relative pos itions and velocities. This can cause some difference in the number of defects formed. Active Volume Effect The size of the active region of the simulation ce ll can affect the defects produced even if the defects do not directly inter act with the thermostat region. This is because the energy density of the system is controlled by the amount of the PKA energy and volume of the active region. The effect of the active system size was inve stigated by comparison of two system sizes: 202020 and 303030 unit cells active region volum e. The volume of the thermostat was kept almost the same, ~2650 atoms since this volume determines the amount of the energy to be removed. The entire systems, therefore, ar e 222222 and 313131 unit cells. Figure 721 is the results of these two simulations in terms of the number of U and O vacancies produced as a function of time. The difference between the two runs are reasonably sim ilar and well within the variation expected from the statis tical uncertainty determined prev iously. This implies that the PAGE 162 162 202020 unit cells for the active region volume is already sufficiently large to perform the radiation damage simulations with 1 keV U PKA. Discussion Necessary modifications to perform radiation damage simulations in UO2 have been implemented in our inhouse MD code. These include the ZBL potential, a new and powerful variable time step algorithm, and the Berendsen th ermostat. Each modifications made is tested individually and also in the combined situation of the actual radiation damage simulation setup. Some of the basic parameters n eeded to perform the radiation damage simulation with 1 keV U PKA were characterized. There are other parameters yet to be explored In the actual simulation of the radiation damage in UO2, the effects of the PKA energy and incident direction need to be understood. The estimated PKA energy in a typical react or condition is a few hundreds of MeV173, 174 which is unattainable in the current MD simulation capabi lity because of the computational requirement for extraordinarily large system sizes. However, in repository condition, the required energy is much lower and is about 80 keV.173, 174 As mentioned before, va n Brutzel group has already performed the simulation at this energy although there were some inconsistencies with their own study with the low energy PKA. The direction of the PKA may also manifest some difference in the number and structure of the defects produced. However, in reality, the PKA travels in arbitrary directions; therefore, the statistics of the vari ation in the results is n eeded to be understood. Another important parameter to explore is th e temperature of the system. Under reactor condition, the temperature of the fuel pellet can rise up from 800 K to 1600 K. At these high temperatures, the thermal diffusion of the defects is important and the number and structure of the radiationinduced defects formed will most likely be affected. PAGE 163 163 The effect of the thermostat thickness was not studied, but it is expect ed not to affect the result least for 1 keV U PKA case. This is beca use essentially the same effect was studied when characterizing the Berendsen time constant. Simulations in single crystal structure are useful but in a real fuel pellet there are a number of defects such as point defects, dislocations, GB s, and surfaces. The eff ects of these defects are particularly interesting since they may indicate the possibility of c ontrolling the radiation properties of UO2. Finally the characterization of the number and structure of the defects and their clusters is interesting but challenging. Experimentally UO2 is known to be very radiation tolerant.207 Our simulations in fact showed that a significant numbe r of the defects were a nnihilated within a very short time (<2 fs) and the number remain steady from thereon. Thus the cr ystalline structure was still retained. However, this may depend on how high the PKA energy is. van Brutzel et al. saw no amorphization observed in their simulation,192 but their defect identification method used is questionable. Our interatomic poten tial is different from van Brutzel et al. and the result will depend on the interaction model. Complex defect structures, defect clusters su ch as 2:2:2 cluster discussed in Chapter 6 cannot be formed in our MD simulation because of the values of the fixed charges of the U and O. Although it is possible to fit a potential to allow the pr oduction of a particular defect cluster, it will not be universal. It is very likely that not all the defect clusters can be allowed using a fixed charge model. Nevertheless, further i nvestigation of the radiation damage in MD simulations will be able to address several issues such as the interaction of defect cascade and point defects and GBs, also influence of the defect cascade on the thermal transport. PAGE 164 164 Table 71. Spline parameters and connection po ints for UU, OO, and UO interactions. r1 () r2 () f0 () f1 () f2 (2) f3 (3) f4 (4) f5 (5) UU 2.0009 1.2003 452.3712 1502.56 1995.418 1301.83 415.5172 51.9609 OO 1.0005 0.1999 9.8149 15.5372 13.5114 12.2027 25.8868 10.4897 UO 1.2003 0.6002 31.64577 143.702 341.1256 408.432 236.589 52.9624 r1 and r2 are the connection between Buckinghamsp line and ZBLspline, respectively. PAGE 165 165 Table 72. The ZBL screening function parameters. i Ai bi 1 0.1818 3.2000 2 0.5099 0.9423 3 0.2802 0.4028 4 0.02817 0.2016 All the values are unitless. PAGE 166 166 Figure 71. Potential kinetic, and total energy per atom from ra diation damage simulations with various fixed time step size.173 PAGE 167 167 0 200 400 600 800 1000 1200 1400 1600 0.00.51.01.5 Distance ()Energy (eV) Buckingham Spline ZBL Combined Figure 72. Potential energy of OO interactions as a function of the interatomic distance. The connections between ZBLSpline and SplineBuckingham are made at 0.2 and 1.0 respectively. The combined potential is used in MD simulations for radiation damage. PAGE 168 168 0.4 0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 0.00.10.20.30.40.5 Time (fs)position () With ZBL Without ZBL Figure 73. Trajectory of the O atom placed at 0.1 away from another O atom at the origin. PAGE 169 169 Figure 74. The calculation of the PKA trajectory in two ways: one large step and two half steps. The final positions from the twoway calcula tions are compared to obtain the error in the calculation and used to es timate the next time step. x ( t ) x ( t+dt ) dt dt /2 dt /2 PAGE 170 170 Figure 75. Flow chart of the determina tion of the variable time step algorithm. Initial configuration Predictor step for all atoms (Identify the most energetic atom) Force calculation on PKA Corrector step on PKA Predictor step Force calculation Corrector step Compare xL( t ) and xS( t ) dx PAGE 171 171 1.6 1.5 1.4 1.3 1.2 1.1 1.0 0204060 Time (fs)Trajectory (a) 0.0001 0.0005 0.001 0.005 0.01 dt=0.1fs (fixed) Figure 76. Trajectory of 1 keV U PKA atom along [100] direct ion in 101010 single crystal UO2. The results from the variable time step and fixed time step simulations are compared. The values in the legend are the values of the error ratio. The data from the fixed time step run with dt=0.1 fs are indicated by the solid line. PAGE 172 172 8 7 6 5 4 3 2 1 0.00.51.01.52.0 Time (fs)Position () 1E04 1E04 5E04 dt=1.E04 (fixed) Figure 77. Trajectory of 80 keV U PKA along [100] directi on in 101010 single crystal UO2. The legend indicates the values of the erro r ratio, and the data from the fixed time step run with dt =1104 fs. PAGE 173 173 1.E05 1.E04 1.E03 1.E02 1.E01 1.E+00 0.00.51.01.52.02.5 Time (fs)Time Step Size (fs) 1E05 1E04 5E04 Figure 78. Evolution of the time step size in the 80 keV U PKA simu lation. The vertical axis is in log scale. PAGE 174 174 Figure 79. The configurati on of the MD simulation cell used in our radiation damage simulations. The central region is called active region where all the interesting dynamic processes take place. The surrounding shaded area is the thermostat region to control the temperature of the system. Thermostat Re g ion ActiveRe g ion PAGE 175 175 0 50 100 150 200 250 300 350 400 010002000300040005000 Time (fs)T (K) Figure 710. Temperature evolution of 101010 unit cells single crystal UO2 during equilibration with the Berendsen thermostat with =400 fs. PAGE 176 176 15.0995 15.0990 15.0985 15.0980 15.0975 15.0970 15.0965 15.0960 15.0955 15.0950 02004006008001000 Time (fs)Total Energy (eV/atom) NVE 1E3 fs 5E2 fs 1E2 fs 5E1 fs 1E1 fs 5E0 fs 1E0 fs Vel. Res. Figure 711. Total energy per atom of the equilibrated single crystalline UO2 at 300 K. The warmer color indicate the larger Berendsen time constant [what does this mean??]. Data from the NVE and velocity resca ling simulations are also included for comparison. PAGE 177 177 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 110100100010000 (fs)Standard Deviation (meV/atom ) Kinetic Potential Total Figure 712. Standard deviation of the kinetic, potential and total energies per atom as a function of time constant, from the equilibrated single crystalline UO2 at 300 K. PAGE 178 178 Figure 713. Defect identif ication method in crystalline structur es. Structure A is the defect free perfect crystal, and structure B is the structure with defect s. Black and grey are used to indicate two different types of atoms. The defects in structure B are numbered as (1) interstitial, (2) vacancy, and (3) antisite The arrows are used to indicate the same coordinates in both A and B. The dot lined circles indicates the sphere for finding the neighbors. structure A structure B 1 3 PAGE 179 179 0 5 10 15 20 25 02004006008001000 Time (fs)Number of VU 0.001 fs 0.01 fs 0.10 fs 0.25 fs 0.50 fs VTS 0 20 40 60 80 100 120 02004006008001000 Time (fs)Number of VO 0.001 fs 0.01 fs 0.10 fs 0.25 fs 0.50 fs VTS Figure 714. Effect of the time step size on radi ation damage simulation in NVE ensemble. The system is 101010 single crystal UO2 with 1 keV U PKA in [100]. The variable time step (VTS) was done w ith the error ratio of 104. PAGE 180 180 0 5 10 15 20 25 30 35 40 45 50 0.00.51.01.52.02.53.0 r0 ()Number of defects VU UI 0 20 40 60 80 100 120 140 160 0.00.51.01.52.02.53.0 r0 ()Number of defects VO OI Figure 715. Number of VU, UI, VO, and OI produced after the 1 keV U PKA radiation damage in a single crystal UO2. The simulation was run for 10 ps. PAGE 181 181 0 5 10 15 20 25 30 35 40 050010001500200025003000 Time (fs)Number of VU 0.10 fs 0.05 fs 0.01 fs 0.005 fs 0 20 40 60 80 100 120 140 160 180 050010001500200025003000 Time (fs)Number of VO 0.10 fs 0.05 fs 0.01 fs 0.005 fs Figure 716. Effect of time step size in radiat ion damage simulation in pNVE. The system is 222222 single crystal UO2 with 1 keV U PKA in [100]. PAGE 182 182 280 300 320 340 360 380 400 0200040006000800010000 Time (fs)Temperature (K) 1 fs 5 fs 10 fs 50 fs 100 fs 500 fs 1000 fs 10000 fs NVE Figure 717. Temperature of the active region during the radiation damage simulations with different Berendsen ti me constants. PAGE 183 183 15.098 15.097 15.096 15.095 15.094 15.093 15.092 15.091 15.090 15.089 0200040006000800010000 Time (fs)Total Energy (eV/atom) 1 fs 5 fs 10 fs 50 fs 100 fs 500 fs 1000 fs 10000 fs NVE Figure 718. Total energy per atom during the radiation damage simulations with different Berendsen time constant. The larger the ti me constant, the warmer the color. PAGE 184 184 0 5 10 15 20 25 30 35 40 0200040006000800010000 Time (fs)Number of VU 1 fs 5 fs 10 fs 50 fs 100 fs 1000 fs 10000 fs NVE 0 20 40 60 80 100 120 140 160 180 0200040006000800010000 Time (fs)Number of VO 1 fs 5 fs 10 fs 50 fs 100 fs 1000 fs 10000 fs NVE Figure 719. Number of uraniu m (top) and oxygen (bottom) v acancies produced in radiation damage simulations as a function of the Berendsen time constant PAGE 185 185 0 5 10 15 20 25 30 35 40 050010001500200025003000 Time (fs)Numver of V u 0 20 40 60 80 100 120 140 160 180 050010001500200025003000 Time (fs)Numver of V o Figure 720. Number of uraniu m (top) and oxygen (bottom) v acancies produced in radiation damage simulations with e quivalent PKA directions. PAGE 186 186 0 5 10 15 20 25 30 35 40 45 050010001500200025003000 Time (fs)Number of Vu 20x20x20 30x30x30 0 20 40 60 80 100 120 140 160 180 200 050010001500200025003000 Time (fs)Number of Vo 20x20x20 30x30x30 Figure 721. The number of U and O vacancies pr oduced by the radiation damage simulations in 202020 and 303030 unit cells active volume. PAGE 187 187 CHAPTER 8 INTERFACIAL THERMAL TRANSPORT IN DIAMOND Introduction Diamond is an ideal laboratory in which to characterize crysta llographic effects on interfacial conductance since it ha s the largest therma l conductivity of any material; thus the detrimental effects of the GBs can be expected to be particularly evident. In addition, diamond is of considerable technological interest.7, 208210 In this chapter, the interfaci al conductance of diamond GBs is determined as a function of temperature and twist angle. Furthermore, th e relation between GB energies and interfacial conductance is elucidated. From our simulations we find that the inte rfacial conductance is extremely high; however, when appr opriately normalized, it is less than that of the corresponding GBs in silicon, a result that we can unde rstand in terms of the bonding at the interfaces. Moreover, we find that the interfacial resistance is an approximately lin ear function of the GB energy. As a consequence, we find that it is possible to construct a ReadShockley type model for the thermal conductance of these boundaries. Single Crystal Diamond In order to test the fidelity of the Tersoff potential for thermaltransport applications, bulk single crystal properties have b een studied. First, energy mini mization of the diamond structure at 0 K yields a lattice parameter of 3.567 and the cohesive energy of 7.36819 eV. These values agrees well with expe rimental values of 3.567 1, 6, 211, 212 and 7.36 eV6, 209. Second, the thermal expansion of single crystal diamond is determined using the Andersen NPT ensemble (refer to Chapter 3). A 444 unitcell system of 512 atoms is equilibrated for 1.25 psec at 300, 500, 750, 1000, 1250, and 1500 K. The temperature dependence of the latti ce parameter is found to be linear ( r2=0.999) and the thermal expansion is determined to be 6.00106 1/K, independent PAGE 188 188 of temperature between 0 and 1500 K (see, Figure 81). Slack and Bartram 212 have compiled thermal expansion data of a wide variety of cer amics and metals includin g diamond. Their value of the linear thermal expansion is of course, temperature dependent ranging from nearly 0 below 100 K to 5.87106 1/K at 1600 K. The lattice parameter calculated from their thermal expansion data is also shown in Figure 81. The discrepancy between simu lation and experimental lattice parameter is at most 0.006 (0.16 %) over th is temperature range. This agreement can be considered as excellent. We can use thus use this potential with some confidence to characterize the thermal conductivity of single cr ystalline diamond using the direct method (Chapter 3). Figure 82 shows the temperature profile obtained from 4 4400 unit cells single crystal diamond after 400 psec. The zaxis is along [001], and the temper ature profile was obtained from the mean kinetic temperature of all of the atoms in each atomic pl ane averaged over 250 fsec of time. In the first 1 psec, the system is equilibrated at 1000 K, and then the heat source and sink are turned on to create steady state temperature gradient. Figur e 83 shows the evoluti on of temperature at z =0 and Lz/2. The planebyplane temperature is calculate d from the mean kinetic energy of atoms in each plane averaged over 250 fsec intervals. This data shows that after ~100 psec the temperature fluctuate within a pproximately 5 K of the nominal temperature of the system of 1000 K. The temperature profile shown in Figure 82 is determined from an average of the kinetic energy of atoms in each plane over the last 200 psec. Once the temperature profile is obtained, the linear temperature gradient is extracted to calculate the thermal conductivity according to Fouriers Law (Equation 22). In order to minimize the error from the noise, it is best to ta ke the linear temperatur e gradient over as long a distance as possible. However, the longer the se ction taken, the more lik ely it is to include PAGE 189 189 regions with a nonlinear temperatur e gradient, where Fouriers law is not valid. To determine, the region over which the temperature gradient is determined, two criteria are applied: (i) keep the center of the linear temperature gradient re gime half way between the heat source and heat sink; (ii) take as wide a re gion as possible without deviati ng from the linear behavior, as measured by a correlation of fit in excess of 0.95. The fit to the temperature gradient is shown in Figure 84. This careful selec tion of the linear regi on reduces the error fr om the noise, but not the error from the discrepancy of the temperature gradient between two regions to the left and the right of the heat source. The di fference of the slopes in Figure 84 is approximately 15 % of the average of the slopes. This is perhaps due to the long wavelength phonons having extremely long relaxation time, which cannot be easily controlle d in the simulation. Therefore, the thermal conductivity calculated from this method typically has a statistical error of the order of 10 15 %. As discussed in the context of UO2, the imposition of a heat current on the system leads to large temperature gradients. It is therefore important to establish th e effects of the size of the heat current and to compare with results from a method in which no heat current is imposed. Shelling et al. performed simulations to determine the thermal conductivity of Si using both equilibrium and nonequilibrium methods. They found that th e conductivities obtained from both methods agree when sufficiently large system is used. They also determined that the effect of heat current density and found that the energy density of ~2104 eV/nm2 per 0.55 fs gives least error in the final result. Since this value may depend on the material of inte rest, we performed seven test runs with heat current densities ranging from 1104 to 8104 eV/nm2 per 0.2 fsec. As apparent from Figure 88, the noise level increases as the current density decreases. The result from 1104 eV/nm2 is shown but excluded from since the maximum difference in the conductivity value PAGE 190 190 was 15 % of the average value and there was no syst ematic trend. A similar weak dependence in the heat flux was observed in previous work on Si GBs.213 Moreover, the pr evious study in Si for perfect crystals showed quantitative ag reement for the predic ted thermal conductivity between the imposed heat current method and a GreenKubo approach, in which there is no heat current at all.213 The heat current density for subse quent simulations was kept constant at 2.5103 eV/nm2fsec. The resulting temperature pr ofile in Figure 811 clearly shows the symmetry of the simulation setup. As discussed in Chapters 3 and 4, the size of the simulation cell can have a large effect on the calculated value for the bulk thermal conductivity. This is particularly true for materials with high thermal conductivities, such as diamond. A series of simulations were performed with different system lengths and temp eratures to investigate the size effect (Figure 85). The thermal conductivity extrapolated to 1/ Lz=0 is shown in Figure 86 as a function of temperature Three experimental data2, 51, 214 are also included in the plot for comparison. The values from Olson and Onn are from isotopically pure (13C<0.07 %) diamond and from Touloukian are for type IIa natural diamond (13C<1.01 %) (Chapter 1). Our simulati on results agree very well with the experimental data. The same da ta on loglog plot shows that the slope of experimental and simulation data to be approximately 1.35 and 1.16, respectively. Although the isotopically pure and natural diamonds give somewhat different values of the thermal conductivity, our results seem to lie right in between the experimental data. Interfacial Conductance of (001) Twist GBs The diamond crystal structure is character ized by four covalent bonds pointing along <111> directions. A GB on an arbitrary plane will thus have two covale nt bonds projecting back into the grain itself and two bonds across the GB (Figure 87).124 The two bonds across the interface can be expected to be significantly di fferent from the bonds in the bulk, due to the PAGE 191 191 structural disorder at the GB. If we examine the principal planes in the diamond lattice, GBs on the (111) plane have three bonds back into the crystal and one across th e interface; GBs on the (011) plane have one bond pointing back into the crystal, two para llel to the interface and one pointing across the interface. Neith er of these is similar to the typical GB described above. However, GBs on the (001) plane have two bonds poi nting back into the cr ystal and two pointing across the interface. A highangle twist boundary on the (001) plane may thus be considered as being, to a considerable extent, representative of the GBs in polycryst alline and nanocrystalline diamond.124 Therefore, to elucidate the structureproperty relationship fo r a generic GB in diamond, we determine how the thermal conductan ce of (001) GBs varies with twist angle. The diamond lattice has a twofold symmetry ax is in the <001> direction. Thus on forming symmetric twist boundaries on this pl ane, a rotation of one semicryst al with respect to the other by 180 will result in a single cr ystal. Therefore, to probe the full range of GBs, we need only consider twist angles from 0 to 90. We have simulated the thermal transport properties of 12 such GBs, with twist angles ranging from 8. 80 to 90.0. These GBs were formed by cutting a single crystal along an (001) plane and rotating one ha lf of the crystal with respect to the other. The unit cells for the GBs considered range in area l dimensions from being equal to that of the perfect crystal ( =90, =1) to being 101 times larger ( =11.42, =101). Lower angle GBs are excluded from this study because of their too large system sizes The 90 GB is a special case a nd is actually a symmetric twin GB. In the <001> direction, the stacking of the perfect crys tal can be written as AaBbA aBbAb , where the Aa and Bb layers come from the fcc Bravais lattice with a basis of two atom s at (0,0,0) and (0.25,0.25,0.25). The 90 twist GB then corresponds to a stack ing sequence of AaBbA bBaAb , in which we can see that the stacking sequence is symmetric about one of the A planes. Full crystallographic PAGE 192 192 details of all of the GBs considered and a summar y of the simulation result s are given in Table 81. To prepare the GBs for the th ermal conductivity simulations, th e unrelaxed GB structures were heated to 2000 K for 10 psec to minimize the st ress, then slowly cooled to near 0 K in a stepwise fashion; equilibrium T=0 K structures were then determined by steepest descents quenches to minimize the energies and to reduce the forces and stresses to negligible values. This process was determined to be appr opriate in earlier work on Si GBs.124 The coordination of atoms with in the typical highangle highenergy GB, 29(001) agree very well with the previous results (Table 61).124 To minimize the computer time, the relaxation process was performed on simulation cells that were only 16 unit cells long. After the re laxed structures were obtained, perfect crystal regions were inse rted to form structures long enough in the z direction for the thermal conductivity simulations. von Alfthan et al.215 investigated the zerotemperature st ructures of Si GB structures by removing atoms from the GB region and found GB structures with energies lower than those previously reported. We have explored this effect by performing the same procedure on a 29(001) GB in diamond, but we did not find any structures with energies lower than those formed by the process described above. Figure 89 shows the GB energies as a function of atoms removed. Initially a certain numb er of atoms are removed from unrelaxed 29(001) GB, and then the structures are rela xed by simulated annealing. Th e annealing and cooling processes are exactly the same as the process described ab ove. Atoms are systematically removed from (001) plane and relaxed to obt ain the GB energies. Since the diamond unit cell contains 4 distinct atomic planes in the [001] direction and a 4fold symmetry along [001], removal of two (004) planes is sufficient to inve stigate all of the possi bilities. Once a comp lete (004) plane is PAGE 193 193 removed, the grain is shifted by 1/ 4 of the unit cell in [001] direc tion to keep the density similar to the original 29(001) GB. This operation does not ma ke any difference in the final GB energies, as von Alfthan et al. s group also confirms. When th e atoms are removed the structure was relaxed by the same procedure described above to relax the GBs. Out data shows that the GB energies do not differ significantly from the relaxed 29(001) GB with the GB energy generally remaining below 8 J/m2. We have also observed a few cases where the GB energies are significantly higher reaching close to 12 J/m2. The 29(001) GB has one of the highest energies and is thus one of the most likely to undergo reconstruction, since the potential energy saving is the largest. Because, our simulations on this GB show no reconstruction, it is unlikely that the other, lower energy, GBs will further lower their energies by reconstruction. Before discussing the dependence of the therma ltransport properties on the twist angle, we examine the representative case of the 29 (001) =43.60 GB. Because, the simulation is set up as a 3D periodic system with crystallographically identical GBs at the cent er of the unit cell and at the edge, the heat source and heat sink are placed equidistant betw een the two GBs. The resulting temperature profile in Fi gure 811 clearly shows the symmet ry of the simulation setup. In our simulation, temperature in each 0.89 thick slice (i.e., a single (004) plane) through in the simulation cell is defined by the average kinetic energy of per atomic plane averaged over the last 560 psec. The temperatur e drop across a GB is quantified in the following manner: (i) the temperature profiles in the bulk regions surrounding the in terface are fit linearly, with the range of z values fit determined by th e minimizing the difference between the absolute values of the slopes. The linear fits are ma de over the region of between 12 nm and 25 nm, depending on the length of the simulation cell. The regions are considered bulk like when r2 of the linear fit is above 0.95. Because of the symm etry of the structure, all the slopes are found to PAGE 194 194 be the same within small error bars. (ii) Th e temperature jumps are determined from the equations to the linear fits. Analysis from the lin ear fits yield the values of 21.1 K and 23.8 K for the temperature drops at the tw o crystallographically identical GBs (see, Figure 811). The average temperature drop from these two GB s is thus 22.5 K; the corresponding thermal conductance is 8.80.6 GW/m2K, with the estimated uncerta inty in the conductance coming from the difference in the two temperature drops. As described above, the size of the simulati on cell has a considerable effect on the calculated thermal conductivity of perfect crystals.213, 216 The origin of this effect is the restriction on the maximum phonon mean free path set by the length of the periodic simulation cell. To explore the possibility of an analogous size effect for the in terface conductance, we have determined the temperature drop at the interface and the co rresponding interfacial conductance for the 29 GB for simulation cells ranging 14.3 nm to 287 nm in length. As Figure 812 shows, there is strong system si ze dependence; however, it appears that T reaches an asymptotic value at large system sizes. To confirm this, we have fit these data by a shifted exponential: ) exp(0Lz/L A T T (82) Figure 812 shows that this functional form fits the data well. The values of parameters A and L0 obtained from the fit are 60.9 K and 70.7 nm, re spectively. The formal description of the Kapitza conductance is given by Stoner et al.81 as an integration over the wave vector summed over all the branches. Due to the complexity of phonon density of states and dependence of transmission coefficients as a function of disper sion relation, extracting an analytical expression of the size dependence is not straightforward. At present, we do not have a theoretical argument PAGE 195 195 for this functional form, but simply justify its use a posteriori by the quality of the fit shown in Figure 812. The infinite limit temperature drop, T, is 17.6K. The calculated interfacial conductances show a corresponding increase w ith an asymptotic value of 12.3 GW/m2K. Interestingly, no such effect was observed in prev ious (albeit less extensive) simulations of this effect for Si GBs;77 however, the thermal conductivity of sing le crystal Si is much lower than that of singlecrystal diamond and the bonding of at oms at the interfaces are quite different, as is described below. Because the simulations for the longer systems are extremely computertime intensive, for our systematic studies of the e ffects of crystallography, we have used a unit ce ll of 142.4 nm (400 unit cells) in length. The va lues of the interfacial conducta nces discussed below are thus ~30 % lower than the asymptotic va lues we would predict for a simulation cell of infinite length. In addition to the effects of the length of the simulation cell, we have considered the effect of its cross section. In the previous study on Si single crys tal, the dependence on the width of simulation cell was found to be very weak beyond 33 unit cells for the cross section 213. The periodic planar unit cell is times larger for GBs than for th e single crystal. As a result the smallest cross sectional area is ~3.53.5 units (for the 29 GB), with most of the simulation cells have larger cross sectional areas. For the case of the 29 GB, we have also determined the temperature dependence of the interfacial conductance. As 813 shows, the interfacial conduct ance increases almost linearly with the temperature up to 1250 K, above which it de creases. This increase is in strong contrast to the bulk thermal conductivity which decreases strongly with temperature (Figure 86). This behavior is actually similar to that in UO2 (see Figure 54). We can understand this in a qualitative manner as arising from the significant obstacle that the GB offers to heat transport. As PAGE 196 196 such, there is poor coupling of the phonon mode s on opposite sides of the interface. As the temperature increases, the anharmonicity of the system increases, which increases the coupling of previously weakly coupled phonon modes acro ss the interfaces, thereby enhancing the thermal transport. This increase with temperature is consistent with experiments and simulation on a wide range of heterointerfaces.129, 130, 217 Detailed analysis shows that the drop in the conductance at 1500 K is due to the structural change at the interface, rather than a change in the ther maltransport mechanism itself. Interfacial Thermal Conductance: Diamond vs. Silicon To determine the effect of crystallogra phy on the interfacial c onductance, Figure 814 shows the calculated thermal c onductance at 1000 K of the twelve diamond GBs as a function of the twist angle for a simulation cell length of 142.4 nm. The interfacial conductance decreases from about 17 GW/m2K for low twist angles to 9 GW/m2K for = 43.6 before increasing again to 12 W/m2K for the 90 1 GB. The uncertainty in the calculated values is represented by the error bars, which are typically 1.3 GW/m2K. This significant structural dependence, whic h we analyze in detail below, is quite reasonable since the degree of structural disorder at the GB, for which the energy is a measure, is lowest at small and high twist angles, for whic h dislocation models of GBs are appropriate19 and highest at intermediate twist angles. In an earlier paper on struct ure and mechanical behavior 124, a comparison of the structure and properties of diamond and si licon GBs was found to be quite instructive. In particular, diamond readily forms sp2 graphitic bonding at GBs, while Si strongly favors the retention of sp3 bonding at the interface, even at the expense of the generatio n of a significant amount of structural disorder. Therefore, in addition to results for diam ond, Figure 814 also includes the PAGE 197 197 previously determined77 values for the interfacial conducta nce of two (001) twist GBs in Si. These values are close to an order of magnit ude lower than for the corresponding diamond GBs, which should not be too surprising since the room temperature thermal conductivity of bulk Si is also much less than that of di amond (150 vs. 2000 W/mK). This di fference in bulk properties can be taken into account th rough the Kapitza length, K KG l / which is the thickness of perfect crystal offering the same thermal resistance as the GB. Figure 815 shows the Kapitza length, given in terms of the lattice parameter of the respective materials, as a function of twist angle. We see that due to its lower si nglecrystal thermal conductivity and its larger lattice parameter, the Si GBs actually have somewhat lower Kapit za lengths, indicating that Si GBs actually allow more efficient heat flow th an the corresponding diamond GBs. To understand the origin in these differences in interfacial conductances, it is necessary to consider the structure of the GBs in more de tail. Figure 816 shows edgeon views of the 29 GB in both Si and C. It is evident from the fi gure, and quantified in Table II, that there is considerably less bonding across th e interface in diamond than in silicon. In particular, in diamond only 14 % of atoms at the GB have two bonds across the interface, while in Si 82% of atoms have two bonds across the interface. 124 Thus structurally, the Si GB is much better connected. This extra bonding across the Si inte rface naturally makes it easier for vibrational excitations, i.e. heat, to cross the interface also. Extended ReadShockley Model For Interfacial Thermal Conductance The ReadShockley (RS) model describes the energies of lowangle GBs in terms of dislocation cores and an associated strain field19. Wolf has made an empirical extension to the RS model the extended ReadShockley (ERS) model18 which fits the en ergies of GBs over PAGE 198 198 the full range of twist angle. Fo r the case of twofold rotational symmetry, as on the (001) plane of diamond, the ERS model is given by: 90 45 for / 2 180 sin ln 2 180 sin 45 0 for / )) 2 ln(sin ( 2 sin b E E E b E E Es c STGB s c gb (84) Ec and Es are the energies associated with di slocation core and strain field. E< and E> are used to indicate the parameters for angle lo wer than or greater than 45 twist angle. ESTGB is the energy of the symmetric twin GB ( =90). b and are the Burgers vector and twist angle, respectively. This functional form guarantees con tinuity in the energy and its first derivative at =45. Figurre 817 shows the calculated GB ener gy for a large number of (001) diamond twist GBs. (There are many more data points in this plot than in the plot fo r the Kapitza conductance because the GB energy calculations are not comput ationally very demanding.) The solid line in Figure 817 is the ERS fit to the data, which quite well reproduces the trend in the energies. Also shown in Figure 817, are the contributions to the energy from the dislocation core (dashed line) and the strain fiel d (dashdot line), as determined from the extended ReadShockley equation. At most twist angles the structural contribution, wh ich at low angles is due to dislocation cores, dominates the energy. In part icular, at 45 and 90, th ere is no strain field contribution. The only boundary condition that rest ricts the value of thes e coefficients is the continuity of the two functions at 45, which yields the constraint, STGB C CE b E b E GBE at 0 and d dEGB at 45 are zero by default. The para meters of the ERS model obtained from the fit are PAGE 199 199 2 2 2 2 2/ 75 1 / 42 5 / 19 4 / 24 5 / 94 5 m J E m J b E m J b E m J b E m J b ESTGB S C S C (85) Interestingly, the general trend in the Kapitza length shown in Figure 815 is similar to that of the GB energy shown in Figure 817. To s how this correlation more clearly, Figure 818 shows that the Kapitza length depends approximate ly linearly on the GB energy; the slope of the linear fit is 85.3 m3/J. This approximate linear relation indicates that the Kapitza length should also be describable by an ERS model: 90 45 for 2 180 sin ln 2 180 sin 45 0 for )) 2 ln(sin ( 2 sin s c STGB s c Kl l l l l l (86) The solid line in 815 is an ERS fit to the Kapitza length. In analogy with the energy of the dislocation core and the strain field energy, it is natural to identify the coefficients of this ERS model as the Kapitza length associated with the st ructural disorder in the GBs and with the strain field. From the fit obtained from the EGB, we determine the coefficients of lK to be the following: 0 0 0 0 00 121 1 71 4 13 3 105 3 134 a l a l a l a l a lSTGB S C S C (87) As in the case of the energy, the ERS form forces the condition C C STGBl l l by continuity. The strain contributions to th e Kapitza length differ by ~30% in contrast to the corresponding values for the energy that differ by only 5%. PAGE 200 200 Discussion The work presented here has shown that the in terfacial (Kapitza) resistance of (001) twist GBs in diamond depends systematic ally on the twist angle. Impor tantly, the Kap itza length is proportional to the energy of the GB, allowing us to fit it with an extended ReadShockley model. These results suggest that the energy of a GB, which is easily calculated, might be used as a surrogate quantity for the inte rfacial conductance, which is mo re difficult and computationally very expensive to determine. Further work on different GBs in diamond, and on different materials systems, is required to establish if this is a general behavior or if it is particular to this system. Regardless of the specific details however, this work does suggest that it should be possible to determine an empiri cal relationship for the interfaci al conductance as a function of the crystallography of the material. Such a relationship will be of considerable value as input to mesoscale simulations of thermal transport.218 PAGE 201 201 Table 81. Summary of the di amond (001) GBs properties. EGB V GK lK PAGE 202 202 Table 82. Comparison of the GB energi es and coordination of the atoms of 29(001) =43.6 GB in diamond and Si. 124 Diamond Si EGB (J/m2) 6.19 1.62 PAGE 203 203 3.565 3.570 3.575 3.580 3.585 3.590 3.595 3.600 3.605 0500100015002000 T [K]a(T) [A] Simulation Slack Figure 81. Lattice parameter of a single crys tal diamond as a function of temperature. Experimental data is taken from the comp ilation of experimental data by Slack and Bartram.212 PAGE 204 204 900 920 940 960 980 1000 1020 1040 1060 1080 1100 020406080100120140 z [nm]T [K] Figure 82. Temperature profile of 44400 uni t cells single crystal diamond after 4 ps. PAGE 205 205 940 950 960 970 980 990 1000 1010 1020 1030 1040 0100200300400 t [psec]T[K] z = 0 z = Lz/2 Figure 83. Time evolution of the temperature at the edge (Lz=0) and center (z=Lz/2) of the simulation cell. PAGE 206 206 T = 0.632z + 1040.3 950 960 970 980 990 1000 1010 1020 1030 1040 1050 405060708090100 z [nm]T [K] T = 0.7156z + 892.51 950 960 970 980 990 1000 1010 1020 1030 1040 1050 120130140150160170180 z [nm]T [K] Figure 84. Slopes of the linear region of the temperature profile shown in Figure 82. PAGE 207 207 0.000 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.010 0.0000.0100.0200.030 1/Lz [1/nm]1/ [mK/W] 100 K 300 K 500 K 1000 K 1500 K Figure 85. Slopes of the linear region of the temperature profile shown in Figure 84. PAGE 208 208 0 500 1000 1500 2000 2500 3000 3500 4000 0500100015002000 T [K]K [W/mK] Simulation Touloukian Olson Onn 2.0 2.5 3.0 3.5 4.0 4.5 1.52.02.53.03.5 log(T)log(K) Simulation Touloukian Olson Onn Figure 86. Thermal conductivity versus temperatur e in both linear and loglog plots from our simulation (filled) and the experimentals2 (open). PAGE 209 209 Figure 87. Bonding of atoms in diamond struct ure along (a) [111], (b) [100] and (c) [110] directions.124 All points are carbon atoms. Filled and empty points are used to distinguish the carbon atoms at FCC sites and tetrahedral inte rstitial sites for clarity. Interplanar distances ar e given by d1 and d2 in the unit of lattice parameter. [Adapted from reference 124 (Figure 6).] (a) (b) (c) PAGE 210 210 0 50 100 150 200 250 300 350 400 450 0.00.20.40.60.8 /A [x103 eV/nm2] [W/mK] Figure 88. Thermal conductivity as a function of heat current density. PAGE 211 211 0 2 4 6 8 10 12 14 0102030405060 Number of atoms removedEGB [J/m2] Figure 89. GB energy as a function of the number of atoms removed. PAGE 212 212 Figure 810. A typical temper ature profile of a system w ith GBs, in this case (001) =43.60 29 symmetric twist boundaries. The GBs are located at the center and the edge of the simulation cell. The average temper ature drop across the GBs is 22.5 K. PAGE 213 213 T = 0.9857z + 1080 1000 1005 1010 1015 1020 1025 1030 1035 1040 1045 1050 4050607080 z [nm]T [K] T = 0.8788z + 1052.1 950 955 960 965 970 975 980 985 990 995 1000 708090100110 z [nm]T [K] T = 0.8814z + 860.84 950 955 960 965 970 975 980 985 990 995 1000 110120130140150 z [nm]T [K] T = 0.8749z + 1010.7 1000 1005 1010 1015 1020 1025 1030 1035 1040 1045 1050 0102030 z [nm]T [K] PAGE 214 214 Figure 811. Linear temper ature gradient from the temperature profile of the 29 (001) =43.60 GB. PAGE 215 215 Figure 812. Size dependence of the thermal conductance of the 29 =43.60 GB fitted to a shifted exponential. T is the asymptotic temperature drop in the infinite size limit indicated by the dotted line. PAGE 216 216 0 1 2 3 4 5 6 7 8 9 10 0500100015002000 T [K]GK [GW/m2K] Figure 813. Temperature dependence of the thermal conductance of the (001) 29 =43.60 GB. The conductance increases almost lin early with the temperature up to about 1250 K. The drop at 1500 K is due to a change in structure. PAGE 217 217 Figure 814. Thermal conductance as a function of twist angle for diamond (001) symmetric twist GBs from our simulations (filled circles), and corresponding Si GBs (open circles)77. PAGE 218 218 Figure 815. Kapitza length in uni ts of their respective lattice parameters as a function of twist angle for both diamond and Si at 1000 K. Filled circles () are for diamond, and the open circles () are for Si. The solid line is th e extended ReadShoc kley fit to the diamond data. The dashed line is the disloc ation core contribution, while the dashdot line is the strain fi eld contribution. Both points for Si are about 30 a0 lower than for diamond. PAGE 219 219 Figure 816. Cross section of 29(001) GB in diamond (a) and Si (b). The freedom for C atoms to choose sp2 and sp3 hybridization results in the rather open structure within the GB. The dashed lines indicate the location of the GBs. a b PAGE 220 220 Figure 817. GB energy (symbols) and the fit to th e extended ReadShockley model (solid line). The dashed line is the dislocation core c ontribution to the GB energy. The dashdot line is the strain field contribution. PAGE 221 221 Figure 818. The Kapitza length is a reasonably linear function of the GB energy. PAGE 222 222 CHAPTER 9 SIMULATIONS OF THERMAL TRANSPORT IN NANOCRYSTALLINE DIAMOND AND COMPARISON WITH EXPERIMENTAL RESULTS1 Introduction As mentioned in Chapter 1, diamond has th e highest thermal c onductivity of any known materials (the thermal codunctivity of bulk si nglecrystal type IIa diamond is about 2200 W/mK). It also exhibits exceptional optical, thermal, tribological and electrochemical properties220. Unfortunately, the use of singlecrysta l diamond is not feasible for many applications because of the cost. Ultrananocrystalline diamond (UNCD)54 however, is an attractive alternative because it can be produced in thin film form at low processing temperatures with very little surface roughness. UNCD coatings consist of diam ond grains 35 nm in size with atomically abrupt high energy GBs.221 UNCD exhibits many exceptional materi als properties, incl uding high hardness and Youngs modulus209, low friction and low surface energies (as deposited) 222, tunable electrical conduction vi a doping with nitrogen223, hermeticity, bioinertness, and resistance to biofouling224. Given the extremely small grain size of UNCD, the thermal conductivity should not be nearly as high as natura l polycrystalline diamond with larg e grain sizes. In fact, if the interfacial resistance at UNC D GBs were similar to hete rophase interfaces (~150 MW/m2K), one would expect the room temperature thermal conduc tivity to be similar to, or less than, that of thermalbarrier coating materials (i.e., < 2 W/mK). It is thus an open question as to whether UNCD thin films can be a useful thermalmanagement material for nanoapplications. 1 In this chapter, I describe a collaboration on ther mal transport in ultrananoc rystalline diamond, including experiments done by the group at Argonne National Laboratory, and simulation work by Watanabe and Phillpot at UF and Bodapati and Keblinski group at RPI. This work was published in Ref. 219M. A. Angadi, T. Watanabe, A. Bodapati, X. C. Xiao, O. Au ciello, J. A. Carlisle, J. A. Eastman, P. Keblinski, P. K. Schelling,S. R. Phillpot, Thermal transport and grain boundary conductance in ultrananocrystalline diamond thin films Journal of Applied Physics 99 (2006). PAGE 223 223 There have been numerous experimental st udies of the thermal conductivity of diamond films prepared by chemical vapor de position. This data is collected in Figure 91. While there is a clear trend for the thermal conductivity to de crease with decreasing grain size, there is considerable scatter. The dashdot line is a fit to these experimental data using a simple thermal resistors model by Yang et al.78 which is the same effective medium model used in Chapter 5 for polycrystalline UO2. This model gives the experiment ally measured th ermal conductivity, in terms of the bulk, conductivity, 0, the grain size, d, and th e interfacial conductance, GK as: d GK 0 01 (91) The best fit to these experimental data yields a value of GK=325 MW/m2.K, a value that is only a little larger than typical values determined for he terophase interfaces 225. Essentially any structural imperfection will de crease the thermal conductivity of a sample. Thus, the significant scatter in these data for a si ngle nominal grain size ma y be attributed to the microstructural variations arisi ng from different film thickness growth conditions, and system geometry. Moreover, the high fraction of GBs, typically containing disordered sp2 bonded carbon, dangling bonds and other defects, can be expected to lower thermal conductivity for nanocrystalline diamond. (As described in the pr evious chapter, roughly 83 % of C atoms are 3 coordinated in MD simulation.) Interestingly, th ere are some data that give significantly higher values of the thermal conductivity than the trend in the rest of the data. Gi ven that any structural defects lower the thermal conductivity, it is natura l to interpret these va lues as being a better measure of the intrinsic thermal conductiv ity of polycrystalline diamond itself. PAGE 224 224 The consistent picture emerging from this stud ies is that the interf acial conductance of GBs in UNCD is extremely high, in excess of 3000 MW/m2K at room temperature, more than ten times that of any material in terface previously characterized. Experimental Approach A series of UNCD films with th ickness ranging from 0.8 7.5 m were synthesized by the Argonne group using microwave plasma chemic al vapor deposition using argonrich Ar/CH4 plasma chemistries.226 The films were deposited at 800C on silicon wafers (500 m thick) using a microwave plasma CVD system at 1275 W power and 200 mbar chamber pressure. To improve the seeding process and minimize porosit y at the UNCDsubstrate interface, for some samples, 10 nm thickness W layers were deposite d by magnetron sputterdeposition, resulting in much denser UNCD microstructure and much smoother surfaces. While, initially, hydrogen (210%) was added to the chamber to stabilize the plasma, after re aching a chamber pressure of 200 mbar the gas composition was maintained at 1.6% CH4 and 98.4% argon. The growth rate was approximately 0.5 m/h. For each seeding technique (i.e., growth on W/Si or on Si), visible (=632 nm) Raman spectroscopy indicated that sample s of different thicknesses had similar film qualities. The films exhibit nanostructures characteri zed by 35 nm diameter diamond grains separated by GBs that contain a mixture of sp3 and sp2 bonding.226 The hardness, Youngs modulus, refractive index, and friction propertie s are all essentially eq uivalent to natural diamond54; indeed several materials properties of UNCD have been shown to be superior those of natural diamond.227 The thermal conductivity wa s measured using the 3 method228 and analyzed using the offset model.229 According to this model, the film is treated as a thermal resistor, which results PAGE 225 225 in an experimental T vs. frequency curve that has the same shape as the curve calculated for uncoated silicon but is offset by an amount Tf given by b l t P Tf2 (92) where t is the film thickness, is the thermal conduc tivity of the film, l is the length of the metal line, b is the width of metal line, and P is the am plitude of the ac power generated in the metal line. Tf represents the difference between the thermal signal for the UNCD coated substrate and the bare silicon substrate. The UNCD film thic kness was determined by crosssection SEM. An important assumption in the analysis employed here is that the thermal transport is onedimensional (i.e., the rate of heat transfer in the direction normal to the film, into the substrate, is much larger than in the plane of the film). Si nce we observe that the thermal conductivity of the UNCD is small compared to that of silicon, this condition is met. Experimental Results Because the 3 experimental technique measures the conductivity thr ough the film, in addition to the effects of the UNCD GBs, it incl udes the effects of the interface between the substrate and the UNCD film. To specifically dete rmine the effects of th e heterophase interface, the Argonne group grew films both on a Si substr ate directly, and on a thin (10 nm) tungsten layer grown on the Si substrate. Thermal conductivities obtained by the Argonne group for a series of UNCD thin films grown on silicon substrates are shown as circles in Figure 92. We first note that for all film thicknesses, the measured thermal conductiv ity ranged from as little as 1.2 W/m.K to as high as 12 W/mK. Interestingly, as shown in Figure 91, this lower value, obtained for thinner films of UNCD grown directly on Si is quite consistent with the extrapolation of the best fit to the previous experimental results. An SEM analysis of the thin film microstructure of these films PAGE 226 226 (inset a to Figure 92) shows that there are large voids at the interface. It thus seems that for these thin films the intrinsic thermal properties of the diamond GBs are not being measured. By contrast, inset B shows that for UNCD grown on an intermediate W layer, there is no porosity at the interfaces, a consequence of the greatly enhanced initial nucleation density at the UNCD/W interface, thereby inhibiting the formation of voids seen at the UNCD/Si interface.230 For these films ~910 W/m.K is obtained; this higher value can be attributed to the absence of porosity. Moreover, the considerable increase in the th ermal conductivity of the UNCD/Si films over the thickness 3.5 5 mm can then be understood as th e result of a change in the interface structure from porous to nonporous. These films yield thermal conductivities as high as 12 W/m.K, which can be taken as an estimate of the lower bound to the true intrinsic thermal conductivity of UNCD. This best estimate is also shown in Fi gure 91, and is a considerable outlier compared with most of the previous experiments. Recall that Yangs model (Equation 54) gives the interfacial conductance by, 0 01 d GK. For an average grain diameter of 4 nm, a bulk conductivity of 2200 W/mK, and a measured thermal conductivity of 12 W/mK, GK is calculated to be 3000 MW/m2K at the measurement temperature of 310 K. This value is more than an orderofmagnitude larger than any previously reported roomtemperature GB or heterophase interface conductance 78, 225. It should be noted that GK was calculated from the measur ed thermal conductivity values under the assumption that Equation 54 is valid. An assumption imp licit in Equation 54 is that the UNCD grain interiors have th e conductivity of single crystal diamond. This assumption is supported by previous HREM obser vations that have indicated that UNCD grain interiors are indistinguishable in stru cture from crystalline.3,4 PAGE 227 227 It can be concluded from these experimental studies that UNCD allows an estimate of the intrinsic thermal properties of diamond GBs. However, the value obtained for GK=3000 W/m2.K is startlingly high compared with other materials systems. Simulation Results To provide additional credibility for this unusually high experimental value for Kapitza conductance obtained in the Argonne groups experiments, atomiclevel simulations of the thermal conductivity of both model diamond nanos tructures of a similar grain size to the experimentally investigated systems, and of an individual grain boundary, re presentative of those in nanocrystalline diamond have been performed. The value of the simula tions is that we can investigate systems without pores or impuritie s, thereby examining the intrinsic thermaltransport properties of diamond. We have constructed two model diamond nanocry stals with grain sizes of 2 nm and 3 nm respectively; recall that the grain size in UNCD is ~3.5 5nm. They are constructed in the same method as UO2 polycrystals (Chapter 5): 3 dimensional co lumnar grains rotated in the plane of the texture. The only difference from the UO2 polycrystals is that there is no restriction of charge neutrality of the system. In all of our simulations on both na nocrystals and individual GBs, diamond is described by the Tersoff manybody potential as described in Chapter 1. The basic bulk single crystal properties of diamond were established and described in the Chapter 7. The polycrystalline diamond structures were co nstructed by Watanabe while the thermal transport simulations were performed by Bodapati at RPI. To determine the thermal conductivity of our model nanostructures, the direct method was used (Chapter 3). The resulting temperature profile for the 2 nm system is shown in Figure 93. The thermal conductivity determined from this system, using Fouriers Law dx dT J is =9.2 PAGE 228 228 W/mK; for the 3 nm grain size, we obtain 13.8 W/mK. This 50% increase in thermal conductivity on increasing the grai n size by 50% is consistent w ith Equation 92 with the bulk conductivity being very high and the interface co nductivity being independent of grain size. These calculated values for the thermal conduc tivity are remarkably consistent with the experimental results for UNCD (at a grain size of 35 nm) of 12 W/m.K, decribed in the last section. Using Equation 54, these simulation da ta yield an average in terfacial conductance of GK= 4600 MW/m2K for the tilt boundaries in our model na nocrystal. We have used this estimate of the intrinsic interfaci al conductance obtained from simulati on to give the solid line in Figure 91. This gives a slightly higher predicted va lue of the thermal conduc tivity of UNCD than we found experimentally, and a higher value than almost all of the experimental values at any grain size. Because the simulation estimate of the interfacial conductance is based on the same rather simple analysis of the model nanocry stal as the experimental data (i.e., Equation 54), it is highly desirable to have a direct measur e of the interfacial resistance itself. While this is not possible experimentally, simulation can again help us (C hapter 8). Figure 94 shows the temperature profile through the GB at room temperature. The heat source a nd sink are at the left and right respectively. The temperature drop at the interface, T, is 25.4K. The smooth temperature decrease in the diamond perfect cr ystal is characteristic of th e bulk thermal conductivity. The imposed heat current was 217 GW/m2, allowing us to determine the Kapitza conductance to be 8500 MW/m2K. Conductance only varies by approximately a factor of two over a wide range of twist angles (Chapter 7), indicating that the sp ecific choice of the grainboundary twist angle is not important for this discussion. PAGE 229 229 This calculated value of the Ka pitza conductance of an isolat ed GB is somewhat larger than that deduced from the simulations of the mode l nanocrystals. This is attributable to (a) the bicrystal contains a si ngle (001) twist boundary, whereas there are a number of different (001) tilt boundaries in the nanocrystal, and (b) the mo re complicated microstructure in the model nanocrystals. Discussion The results from experiment and simulation pr ovide a reassuringly co nsistent picture of interfacial thermal transport in diamond, indicatin g that the interfacial co nductance is at least 3000 MW/m2K, more than tentimes that previous ly seen: specifically, in reviewing the literature, Cahill et al.225 noted that the interfacial conductance for a number of heterophase interfaces spanned th e range of 20200 MW/m2K. In the only previous study of a homophase system, Yang et al. found that the interfacial resistance of yttriastabilized zirconia lies at the high end of this range 78. When we put the diamond inte rfacial conductance into the ap propriate context, however, the interfacial conduc tance, though larger than any previous value, is not in fact anomalously large. As discussed in Chapter 7, materialsi ndependent measure of the interfacial conductance is provided by the Kapitza length lK=0/GK, which measures the thic kness of perfect crystal offering the same resistance to heat transport as the interface. Taking 0=2200W/mK and GK =3000 MW/m2K for UNCD, yields lK730nm. This is surprisingly si milar to the Kapitza length for the aluminumalumina interface 225: GK=170 MW/m2.K and = ( 0 Al+ 0 alumina) = (237+36)=150 W/mK, gives a Kapitza length of lK =800 nm.. The determination of the Kapitza length in diamond to be of the order of 1 m is consistent with thermal transport measurements of conve ntional diamond polycrystalline films which have PAGE 230 230 demonstrated that to reach bul klike conductivities grain size has to be larger than several microns. In this large grain size limit, th ermal resistance due to GBs becomes small by comparison with the intrinsic bulk resistance 53. It is interesting to note th at a UNCD conductiv ity of 12 W/mK is similar to the highend conductivity of what Morath et al. referred to as amorphous diamond 231, an sp3bonded glass structure. Both materials have much higher thermal conductivities th an either diamondlike carbon (a sp3bonded glass structure but with a large fraction of bonds terminated by hydrogen), or amorphous carbons c ontaining significant sp2 bonding component 232. This comparison indicates that while crystall ine grain interiors in UNCD e nhance thermal conductivity with respect to amorphous carbon stru ctures, the significant sp2bonding component in UNCD GBs offer higher thermal resistance than the co rresponding layer of an amorphous diamond sp3bonded structure. We have determined values of 12 W/mK, 3000 MW/m2K, and 750 nm for the room temperature thermal conductivity, in terfacial conductance, and Kap itza length, respectively of ultrananocrystalline diamond. The interfacial conductance is an order of magnitude higher than any previously reported, but th e relative interface conductance compared is not anomalously high, as indicated by a characteri stic Kapitza length similar to, or larger than, that of other materials. Although the thermal transport coefficient of UNCD is not as high as that of polycrystalline diamond with larger grains, the combination with a much smoother surface morphology makes these films very attractive fo r many applications. These include thermal spreaders for applications such as cooling of mi crochips in the next ge neration of high density CMOS devices, control of temperature in the he ads of magnetic storage media, and MEMS and PAGE 231 231 NEMS structures with combined optimum mechan ical and thermal stabi lity properties. Diamond structural layers with different grain sizes could be integrated together into a microsystem to effectively pipe heat to specific locations. PAGE 232 232 Figure 91. Compilation of diamond th ermal conductivity versus grain size d. (a)typical single crystal, (b) Ref: 233 (c) Ref: 234 (12) (d) Re f: 235(e) Ref: 236 (f) Ref: 237, (g) Ref: 238, (h) Ref: 239, (i) Ref: 240, (j) Ref: 52, (k) Ref: 241, (l) Ref: 242, (m) Ref: 239, (n) Ref: 243, (o) Ref: 244, (p) data from this study, (q) Equation (91) for GK=4600 MW/m2.K obtained from simulation of UNCD, (r) Fit of Equation 92 to data (a) (o) yields GK =325 MW/m2.K. Compilation and analysis by Watanabe. 219 PAGE 233 233 Figure 92. UNCD thermal conductivity as a function of film thickness. UNCD films grown on W/Si; UNCD films on Si. Insets show SEM cros ssections the film structure: (a) UNCD/Si with voids and, (b) UNCD/W/Si w ith a dense interface. Each data point corresponds to a single sample of the indi cated average grain size. Each error bar indicates the standard deviation from multiple measurements of that particular sample. This work is done by the Argonne group.219 PAGE 234 234 Figure 93. Microstructure (top ) and temperature prof ile (bottom) of polyc rystalline diamond. Model structure was created by Watanabe, the thermal conductivity simulation was performed by Bodapati at RPI.219 PAGE 235 235 270 280 290 300 310 320 330 525762677277828792z [nm]T [K] Figure 94. Temperature profile through a diamond 29 =43.60 twist GB determined by simulation. PAGE 236 236 CHAPTER 10 CONCLUSIONS A systematic investigation of the thermal transport properties in UO2 and diamond has been done using atomic level simulation technique. The main focus was to establish the basic physical properties of bulk single crystal materials, and to elucidate the effect of defects such as point defects and GBs on the thermal transport of UO2 and diamond. A number of contributions to our understanding of heat transfer in def ected crystalline solids have been made. The thermal transport properties of UO2 single crystal was es tablished using MD simulation technique. The theory of phonon medi ated heat transfer was used to relate anharmonicity and thermal conductivity of diel ectric solid. Results were compared with available experimental data and good agreement was obtained. The temperature and grain size dependence of the thermal conductivity of polycrystalline UO2 was investigated and characterized. An ef fective medium model for polycrystalline solids was used to calculate the interfaci al conductance. Analysis of th e Kapitza length showed that the interfacial resistance si gnificantly dominates the bulk therma l conductivity of very fine grain polycrystalline UO2. The dependence of the thermal conductiv ity of on offstoichiometry in UO2x was obtained at 800 K and 1600 K. The thermal conductivity of UO2+x was found to decrease with increasing concentration of the oxygen interstitials and eventually reac h a plateau when x>0.125. At low concentrations, the effect of oxygen vacancies and interstitials to thermal transport were found to be similar. Analysis using the theory of disordered solid supports this similarity: there were no significant difference in the vibraitonal modes at low concentrations of oxygen defects. A small fraction of localized modes was found only at very high concentration of interstitials. In PAGE 237 237 addition, it was found that the e ffect of different U isotopes on the thermal transport of UO2 was almost negligible. A molecular dynamics simulation code for the radiation damage simulations was developed. The inhouse MD code was modified to include ZBL potential, variable time step, and Berendsen thermostat. Each modifi cation was independently verified. A systematic study of the inte rfacial conductance of the diam ond GB was performed. The thermal conductance of the twist grain boundary was obtained as a function of twist angle. The grain boundary energy was found to follow the ex tended ReadSchockley (ERS) model. An ERS type model for the interfacial conductance was devloped. The relation between the grain boundary structure and interfacial conductance was characterized through the comparison of GBs in Si and diamond. Consistent values of the in terfacial conductance were obtai ned from both experiment and two simulations with distinct microstructures. The coordination number analysis also showed agreement between the experiment and simulation. Thermal transport in both UO2 and diamond are mediated by phonons as described in Chapter 1 and 2. Phonon mediated heat transfer is significantly affected by the presence of point defects and grain boundaries regard less of what the material was. The anharmonic analysis shown to work on UO2 should apply to any electrical insula tors because the theory is based on the phonon mediated thermal tran sport theory in bulk sing le crystalline solid. Both diamond and UO2 grain boundaries signif icantly reduced the th ermal conductivity of the solid. Comparing the ratios of the therma l conductivities of polycry stal of similar grain sizesd to the singlecrystal valu es, we find that the ratio is about ~12/500=1/80 for diamond and ~2/35=1/12.5 for UO2. Thus diamond grain boundaries of more significant obstacles to heat PAGE 238 238 flow than do GBs in UO2. The reason for this difference is most likely the difference in mean free path in these two materials. Since di amond naturally had larger mean free path, the reduction in thermal conductivity is more pronounced. A comparison of the effects of point def ects and grain boundaries is instructive; UO2.125 gives a minimum thermal conductivity of ~1.5 W/ mK at 300 K. By cont rast polycrystalline UO2 shows a thermal conductivity of ~0.7 W/mK at 300 K. This is inte resting since the effect of the oxygen interstitial has been satu rated and it seems that the thermal conductivity cannot be lowered with further increases in offstoich iometry. However, the 3.8 nm grain size polycrystalline UO2 has a thermal conductivity of almost half this value. This suggests that GBs are inhererntly stronger scatterers of phonons than oxygen interstitials. This work also suggests a number of directions for future work. One of the most important improvements one can make on the UO2 model is the interatomic potential. If the thermal expansion and bulk modulus are impr oved, the thermal conductivity of UO2 should be able to give values much closer to the experimental data. Since the effects point defect s and GBs have been identifi ed separately, it would be interesting to see what the effect of the combina tion of the two. How does the interaction of the defects affect the thermal transport? What is the effect of defects in GBs on the thermal conductivity? The same question applies to the radiation damage. Since real UO2 has point defects and GBs, the effects of the presence of point defect s and/or grain boundaries in radiation damage are of significant importance. Will these defects reduce or enhance the radiation resistance? In connection to both to the th ermal transport and radiation da mage, the effect of defect clusters will be intere sting to look at. Radiation damage in UO2 will create significant number of PAGE 239 239 defects with a variety of sizes a nd distributions of defect cluste rs. Can we categorize and make meaningful analysis of the effects of de fect cluster on the thermal transport? Finally, these studies show the power of at omiclevel simulations for elucidating the thermal transport properties of electronic insulators Also they offer the promise for helping the development of better thermal mana gement methods in solids. PAGE 240 240 APPENDIX LATTICE DYNAMICS Lattice dynamics (LD) is the theo ry of lattice vibrations in a solid. It allows detailed inspection of frequency of vibrational modes, a nd the indefication of the atoms associated with each mode. The following survey of the pertin ent aspects of lattice dynamics are based on discussions by Allen et al.85. The basis of LD lies in solving the equation of motion in the form of eigenvalue problem. Imagine a solid consisting of M atoms. If the structure is crystalline, the solid can be thought of as made of smaller set of N atoms. Thus N can be the number of basis atoms, or number of atoms in a unit cell. For example, a typical equati on of motion of point particles in a lattice can be cast into the form of, ) ( ) ( ) ( ) (, 2k e k D k e ki i (I1) where ) (2k and ) (,k ei are the eigenvalue and eigenvector. ) ( k D is a tensor of 2nd rank and known as the dynamical matrix, and given by, 'exp 1 ) (l j jl l j jr r k i l l j j m m k jj D (I2) Here j and j refer to the atoms within a repeat unit, and l and l are for the index of the repeat unit. mj is the mass of atom j, rjl is the coordinate of atom j in cell l. and indicate the direction in Cartesian coordi nates. Thus, the size of the dynamical matrix is 3N3N. is called the force constant matrix and is given by, ) ( ) ( '2l j u jl u V l l j j (I3) under the harmonic approximation. V is the total potential energy, and ) ( jl u is the displacement of atom j in cell l in direction. The size of fo rce constant matrix is 3M3M Note that when 0 k, the dynamical matrix is the force c onstant matrix divided by the mass. PAGE 241 241 The solution of the eigenvalue equation, solution of equation of motion (I1 ) is in the form of plane waves: t r k i e jl ujl j exp ) (, ,, (I4) Here ) (,k ej is the polarization vector which describes the direction in which atoms move. In (I4), is used explicitly to indi cate the vibra tional mode. The analysis on the localization of vibr ational modes can be quantified using the eigenvectors since the eigenvector of the vibrational mode contains all the spatial information. The participation ratio, p, can be calculated using the following relation: 2 1 i i ie e N p (I5)p takes values between 1/N and 1 depending on the degree of localization of the vibrational mode. If the vibration is localized p is small (~1/N), and if not, p will be of order 1. To further analyze the contribution of each element, it is helpful to decompose the contribution as following: A i i i Ae e p p ,, (I6) where A is the label for the species. To compare the polarizati on of particular modes, ie, should be appropriately normalized. The normalized pol arization can be given by i i i ie e e e, ,, (I7) This allows us to distinguish the nature of vibrat ional modes. As discussed in Chapter 2, in the case of highly disordered solid propagon modes are the only ones that have well defined PAGE 242 242 polarization vector. Thus, ie, can be used to distinguish propa gon and rest of the vibrational modes can be distinguished. PAGE 243 243 LIST OF REFERENCES 1D. R. Lide, Handbook of Chemistry and Physics (CRC Press, Cleaveland, OH., 2007). 2Y. S. Touloukian, Thermophysical Properties of Matter (IFI/Plenum, New York, 19701979). 3J. K. Fink, Thermophysical Properties of Uranium Dioxide, Journal of Nuclear Materials 279, 1 (2000). 4C. Haertling,R. J. Hanrahan, Literature Review of Ther mal and Radiation Performance Parameters for Hightemperat ure, Uranium Dioxide Fu eled Cermet Materials, Journal of Nuclear Materials 366, 317 (2007). 5J. H. Yang, K. W. Song, K. S. Kim,Y. H. Jung, A Fabrication Technique for a UO2 Pellet Consisting of UO2 Grains and a Continuous W Channel on the Grain Boundary, Journal of Nuclear Materials 353, 202 (2007). 6R. Berman, The Properties of Diamond (Academic Press, London, 1979). 7R. F. Davis, Diamond Films and Coatings (William Andrew P ublishing/Noyes, 1993). 8Y. Baer,J. Schoenes, ElectronicStructure and Coulom b CorrelationEnergy in UO2 SingleCrystal, Solid State Communications 33, 885 (1980). 9Y. Yun, H. Kim, H. Lim,K. Park, Electronic Structure of UO2 from the Density Functional Theory with onsite Coulomb Repulsion, Journal of The Korean Physical Society 50, 1285 (2007). 10C. D. Clark, P. J. Dean,P. V. Harris, Intrinsic Edge Absorption In Diamond, Proceedings of The Royal Society of London Series AMa thematical and Physical Sciences 277, 312 (1964). 11C. K. Gupta, Materials in Nuclear Energy Applications (CRC Press, Boca Raton, 1989). 12K. Kapoor, S. V. R. Rao, Sheela, T. Sanyal,A. Singh, Study on Solid Solubility of Gd in UO2 using Xray Diffraction, Journal of Nuclear Materials 321, 331 (2003). PAGE 244 244 13Z. A. Munir, The Electrical Conductivity of Doped and Undoped Uranium Oxide, International Journal of Thermophysics 2, 177 (1981). 14J. H. Yang, K. S. Kim, K. W. Kang, K. W. Song,Y. H. Jung, Electrical Conductivity and Nonstoichiometry in the (U,Gd)O2+/x System, Journal of Nuclear Materials 340, 171 (2005). 15M. Amaya, M. Hirai, T. Kubo,Y. Korei, Thermal Conductivity Measurements on 10 wt% Gd2O3 Doped UO2+x, Journal of Nuclear Materials 231, 29 (1996). 16M. Hirai,S. Ishimoto, Thermal Diffusivities And ThermalConductivities of UO2Gd2O3, Journal of Nuclear Science and Technology 28, 995 (1991). 17J. M. Howe, Interfaces in Materials (John Wiley & Sons, Inc, New York, 1997). 18D. Wolf,S. Yip, Materials Interfaces: Atomicl evel Structure and Properties (Springer, 1992). 19W. Shockley,W. T. Read, Quantitative Predictions From Dislocation Models Of Crystal Grain Boundaries, Physical Review 75, 692 (1949). 20W. I. Finch, Uranium Fuel for Nuclear Energy, U.S. Geological Survey, (2002), 1 21D. Strahan, Uranium in Glass, Glazes and Enamels: History, Identification and Handling, Studies in Conservation 46, 181 (2001). 22M. T. Hutchings, HighTemperature Studies of UO2 and ThO2 Using NeutronScattering Techniques, Journal of the Chemical Soci etyFaraday Transactions II 83, 1083 (1987). 23D. G. Martin, The ThermalExpansion of Solid Uo2 and (U, Pu) Mixed Oxides a Review and Recommendations, Journal of Nuclear Materials 152, 94 (1988). 24B. T. M. Willis, Neutron Diffraction Studies of Actinide Oxides.1. Uranium Dioxide and Thorium Dioxide at Room Temperature, Proceedings of the Royal Society of London Series aMathematical and Physical Sciences 274, 122 (1963). PAGE 245 245 25B. T. M. Willis, Neutron Diffraction Studies of Actinide Oxides.2. Thermal Motions of Atoms in Uranium Dioxide and Thorium Dioxide betw een Room Temperat ure and 1100 Degrees C, Proceedings of the Royal Society of London Seri es aMathematical and Physical Sciences 274, 134 (1963). 26R. D. Shannon, Revised Effective IonicRadi i and Systematic Studies of Interatomic Distances in Halides and Chalcogenides, Acta Crystallographica Section A 32, 751 (1976). 27L. Pauling, The Principles Determining the Structure of Complex Ionic Crystals, Journal of the American Chemical Society 51, 1010 (1929). 28L. Pauling, General Chemistry (Dover Publications, San Francisco, 1988). 29B. J. Lewis, W. T. Thompson, F. Akba ri, D. M. Thompson, C. Thurgood,J. Higgs, Thermodynamic and Kinetic Modell ing of Fuel Oxidation Behavi our in Operating Defective Fuel, Journal of Nuclear Materials 328, 180 (2004). 30P. Y. Chevalier, E. Fischer,B. Cheynet, Progress in the thermodynam ic modelling of the OU binary system, Journal of Nuclear Materials 303, 1 (2002). 31C. Gueneau, M. Baichi, D. Labroche, C. Chatillon,B. Sundman, Thermodynamic Assessment of the UraniumOxygen system, Journal of Nuclear Materials 304, 161 (2002). 32J. Schoenes, ElectronicTransitions, CrystalField Effects and Phonons in UO2, Physics ReportsReview Section of Physics Letters 63, 301 (1980). 33N. C. Pyper,I. P. Grant, Studies In Multiconfiguration Dira cFock Theory.3. Interpretation Of The ElectronicStructure Of Neutra l And Ionized States Of Uranium, Journal of The Chemical SocietyFaraday Transactions II 74, 1885 (1978). PAGE 246 246 34C. R. A. Catlow, PointDefect and Electronic Pr operties of UraniumDioxide, Proceedings of the Royal Society of London Seri es aMathematical Physical and Engineering Sciences 353, 533 (1977). 35C. Ronchi,G. J. Hyland, Analysis of Recent Measurements of the HeatCapacity of UraniumDioxide, Journal Of Alloys And Compounds 213, 159 (1994). 36P. G. Lucuta, H. Matzke,I. J. Hastings, A Pragmatic Approach to Modelling Thermal Conductivity of Irradiated UO2 Fuel: Review and Recommendations, Journal of Nuclear Materials 232, 166 (1996). 37R. W. Grimes,C. R. A. Catlow, The Stability of FissionPro ducts in UraniumDioxide, Philosophical Transactions of the Royal Society of London Series AMathematical Physical and Engineering Sciences 335, 609 (1991). 38C. Kittel, Introduction to Solid State Physics (John Wiley & Sons, Inc., 1996). 39C. Ronchi, M. Sheindlin, M. Musella,G. J. Hyland, Thermal Conductivity of Uranium Dioxide up to 2900 K from Simultaneous Measurement of the Heat Capacity and Thermal Diffusivity, Journal of Applied Physics 85, 776 (1999). 40C. T. Walker, D. Staicu, M. Sheindlin, D. Papaioannou, W. Goll,F. Sontheimer, On the Thermal Conductivity of UO2 Nuclear Fuel at a High Burnup of around 100 MWd/kgHM, Journal of Nuclear Materials 350, 19 (2006). 41J. Spino, K. Vennix,M. Coquerelle, Detailed Characterisation of the Rim Microstructure in PWR Fuels in the Burnup Range 4067 CWd/tM., Journal of Nuclear Materials 231, 179 (1996). 42P. Losonen, On the Behaviour of Intr agranular Fission Gas in UO2 Fuel, Journal of Nuclear Materials 280, 56 (2000). PAGE 247 247 43H. Matzke, H. Blank, M. Coquerelle, K. Lassma nn, I. L. F. Ray, C. Ronchi,C. T. Walker, Oxide Fuel Transients, Journal Of Nuclear Materials 166, 165 (1989). 44M. T. Aybers, A. A. Aksit, S. Akbal, S. Ekinci, A. Yayli, L. Colak, A. Van, B. Kopuz,N. Ates, in Euro Ceramics Viii, Pts 13 (Trans Tech Publications Lt d, ZurichUetikon, 2004), Vol. 264268, p. 985. 45K. C. Radford,J. M. Pope, UO2 Fuel Pellet Microstructure Modification Through Impurity Additions, Journal Of Nuclear Materials 116, 305 (1983). 46K. Une, M. Hirai, K. Nogita, T. Hosokawa, Y. Suzawa, S. Shimizu,Y. Etoh, Rim Structure Formation and High Burnup Fuel Behavior of Largegrained UO2 Fuels, Journal of Nuclear Materials 278, 54 (2000). 47K. Kurosaki, Y. Saito, H. Muta, M. Uno,S. Yamanaka, Nanoindentation Studies of UO2 and (U,Ce)O2, Journal of Alloys and Compounds 381, 240 (2004). 48F. Dherbey, F. Louchet, A. Mocellin,S. Leclercq, Elevated Temperature Creep of Polycrystalline Uranium Dioxide: From Micr oscopic Mechanisms to Macroscopic Behaviour, Acta Materialia 50, 1495 (2002). 49V. Blank, M. Popov, G. Pivovarov, N. Lvova, K. Gogolinsky,V. Reshetov, Ultrahard and Superhard Phases of Fuller ite C60: Comparison with Diamond on Hardness and Wear, Diamond and Related Materials 7, 427 (1998). 50W. Kaiser,W. L. Bond, Nitrogen, a Major Impurity in Common TypeI Diamond, Physical Review 1, 857 (1959). 51J. R. Olson, R. O. Pohl, J. W. Vandersande A. Zoltan, T. R. An thony,W. F. Banholzer, ThermalConductivity of Diamond betw een 170 and 1200K and the Isotope Effect, Physical Review B 47, 14850 (1993). PAGE 248 248 52K. E. Goodson, O. W. Kading,R. Zachai, in Proceedings of the Americ al Society of Mechanical Engineers Conference, Heat Transfer Division, edited by R. A. W. e. al., 1994), Vol. 292, p. 83. 53K. M. Leung, A. C. Cheung, B. C. Liu, H. K. Woo, C. Sun, X. Q. Shi,S. T. Lee, Measuring thermal conductivity of CVD diamond and di amondlike films on silicon substrates by holographic interferometry, Diamond and Related Materials 8, 1607 (1999). 54D. M. Gruen, Nanocrystalline diamond films, Annual Review of Materials Science 29, 211 (1999). 55International Technology Roadmap for Semiconductors, (2005) 56International Technology Roadm ap for Semiconductors: Update, (2006) 57E. Pop, S. Sinha,K. E. Goodson, Heat Generation and Transport in Nanometerscale Transistors, Proceedings of The IEEE 94, 1587 (2006). 58K. Oshima, S. Cristoloveanu, B. Guillaumot, H. Iwai,S. Deleonibus, Advanced SOI MOSFETs with Buried Alumina and Ground Plane: Selfheating and Shortchannel Effects, SolidState Electronics 48, 907 (2004). 59L. T. Su, J. E. Chung, D. A. Antoniadis, K. E. Goodson,M. I. Flik, Measurement And Modeling of SelfHeating in SOI NMOSFETs, IEEE Transactions On Electron Devices 41, 69 (1994). 60N. W. Ashcroft,D. N. Mermin, Solid State Physics (Saunders College Publishing, 1976). 61A. Beiser, The Concepts of Modern Physics (McGrawHill, Inc., New York, 1995). 62J. D. Jackson, Classical Electrodynamics (Wiley, New York, 1998). 63F. Kreith,M. S. Bohn, Principles of Heat Transfer (West Publinshing, St. Paul, 1993). 64D. R. Clarke,C. G. Levi, Materials Design for the Next Ge neration Thermal Barrier Coatings, Annual Review of Materials Research 33, 383 (2003). PAGE 249 249 65R. Siegel,C. M. Spuckler, Analysis of Thermal Radiation Ef fects on Temperatures in Turbine Engine Thermal Barrier Coatings, Materials Science and Engineer ing A Structural Materials Properties Microstructure and Processing 245, 150 (1998). 66X. Huang, D. M. Wang, P. Patnaik,J. Singh, Design and Computational Analysis of Highly Reflective Multiple Layered Thermal Barrier Coating Structure, Materials Science and Engineering A Structural Materials Pr operties Microstructure and Processing 460, 101 (2007). 67G. J. Hyland, ThermalConductivity of Solid UO2 Critique and Recommendation, Journal of Nuclear Materials 113, 125 (1983). 68K. Bakker, H. Kwast,E. H. P. Cordfunke, The Contribution of Ther malRadiation To The ThermalConductivity of Porous UO2, Journal of Nuclear Materials 223, 135 (1995). 69S. L. Hayes,K. L. Peddicord, Radiative HeatTransfer In Porous UraniumDioxide, Journal of Nuclear Materials 202, 87 (1993). 70J. M. Ziman, Electrons and Phonons (Oxford University Press, London, UK, 1960). 71R. E. Peierls, Quantum Theory of Solids (Oxford University Press, Oxford, 1955). 72C. Herring, Role of LowEnergy Phonons in Thermal Conduction, Physical Review 95, 954 (1954). 73H. B. G. Casimir, Note on the conduction of heat in crystals, Physica 5, 495 (1938). 74P. G. Klemens, The Scattering of LowFrequency Latt ice Waves by Static Imperfections, Proceedings of the Physical Society of London Section A 68, 1113 (1955). 75E. Merzbacher, Quantum Mechanics (Wiley, New York, 1997). 76P. G. Klemens, Thermal Conductivity and Lattice Vibrational Modes, Solid State PhysicsAdvances in Research and Applications 7, 1 (1958). PAGE 250 250 77P. K. Schelling, S. R. Phillpot,P. Keblinski, Kapitza Conductance and Phonon Scattering at Grain Goundaries by Simulation, Journal of Applied Physics 95, 6082 (2004). 78H. S. Yang, G. R. Bai, L. J. Thompson,J. A. Eastman, Interfacial thermal resistance in nanocrystalline yttriastabilized zirconia, Acta Materialia 50, 2309 (2002). 79E. T. Swartz,R. O. Pohl, ThermalBoundary Resistance, Reviews of Modern Physics 61, 605 (1989). 80R. M. Costescu, M. A. Wall,D. G. Cahill, Thermal Conductance of Epitaxial Interfaces, Physical Review B 67 (2003). 81R. J. Stoner,H. J. Maris, Kapitza Conductance and HeatFlow between Solids at Temperatures from 50 to 300 K, Physical Review B 48, 16373 (1993). 82C. Kittel, Interpretation of the Thermal Conductivity of Glasses, Physical Review 75, 972 (1949). 83G. A. Slack, Solid State Physics (Academic Press, New York, 1979). 84D. G. Cahill,R. O. Pohl, HeatFlow and Lattice Vibrations in Glasses, Solid State Communications 70, 927 (1989). 85P. B. Allen, J. L. Feldman, J. Fabian,F. Wooten, Diffusons, Locons and Propagons: Character of Atomic Vibrati ons in Amorphous Si, Philosophical Magazine BP hysics of Condensed Matter Statistical Mechanics Electronic Optical and Magnetic Properties 79, 1715 (1999). 86J. L. Feldman, M. D. Kluge P. B. Allen,F. Wooten, ThermalConductivity and Localization in Glasses Numerical Study of a Model of AmorphousSilicon, Physical Review B 48, 12589 (1993). PAGE 251 251 87P. K. Schelling,S. R. Phillpot, Mechanism of Thermal Transport in Zirconia and Yttriastabilized Zirconia by Mol ecularDynamics Simulation, Journal of The American Ceramic Society 84, 2997 (2001). 88C. J. Glassbrenner,G. A. Slack, Thermal Conductivity of Silicon + Germanium from 3 Degrees K to Melting Point, Physical Review A General Physics 134, 1058 (1964). 89J. A. Carruthers, T. H. Geballe H. M. Rosenberg,J. M. Ziman, The Thermal Conductivity of Germanium and Silicon between 2DegreesK and 300DegreesK, Proceedings of the Royal Society of London Series aMathem atical and Physical Sciences 238, 502 (1957). 90O. L. Anderson,I. Suzuki, Anharmonicity of 3 Minerals at HighTemperature Forsterite, Fayalite, and Periclase, Journal of Geophysical Research 88, 3549 (1983). 91M. Lucht, M. Lerche, H. C. Wille, Y. V. S hvyd'ko, H. D. Ruter, E. Gerdau,P. Becker, Precise measurement of the lattice param eters of alphaAl2O3 in the temp erature range 4.5250 K using the Mossbauer wavelength standard, Journal of Applied Crystallography 36, 1075 (2003). 92S. R. Phillpot,A. J. H. McGaughey, Introduction to the Thermal Transport, Materials Today, (2005), 18 93R. M. Martin, Electronic Structure Basic Theory and Practical Methods (Cambridge University Press, Cambridge, 2004). 94R. D. Cook, D. S. Malkus, M. E. Plesha,R. J. Witt, Concepts and Applications of Finite Element Analysis (Wiley, New York, 2001). 95L. Q. Chen, Phase Field Models for Mi crostructure Evolution, Annual Review of Materials Research 32, 113 (2002). 96M. P. Allen,D. J. Tildesley, Computer Simulation of Liquids (Oxford University Press, 1987). PAGE 252 252 97W. H. Press, B. P. Flannery, S. A. Teukolsky,W. T. Vetterling, Numerical Recipes in Fortran 77 (Cambridge University Press, Cambridge, 1992). 98D. Frenkel,B. Smit, Molecular Dynamics Simulation (Academic Press, 2002). 99D. Wolf, P. Keblinski, S. R. Phillpot,J. Eggebrecht, Exact Method for the Simulation of Coulombic Systems by Spheri cally Truncated, Pairwise r1 Summation, Journal of Chemical Physics 110, 8254 (1999). 100R. A. Buckingham, The Classical Equation of State of Gaseous Helium, Neon and Argon, Proceedings of the Royal Society of London Seri es aMathematical and Physical Sciences 168, 264 (1938). 101J. Tersoff, Empirical Interatomic Potential for Ca rbon, with Applications to AmorphousCarbon, Physical Review Letters 61, 2879 (1988). 102J. Tersoff, New EmpiricalModel for the Stru cturalProperties of Silicon, Physical Review Letters 56, 632 (1986). 103G. V. Lewis,C. R. A. Catlow, Potential Models For Ionic Oxides, Journal of Physics CSolid State Physics 18, 1149 (1985). 104C. R. A. Catlow, K. M. Diller,M. J. Norgett, Interionic Potentials For AlkaliHalides, Journal of Physics CSolid State Physics 10, 1395 (1977). 105H. C. Andersen, MolecularDynamics Simulations at C onstant Pressure andor Temperature, Journal of Chemical Physics 72, 2384 (1980). 106M. Parrinello,A. Rahman, Polymorphic Transitions in SingleCrystals a New MolecularDynamics Method, Journal of Applied Physics 52, 7182 (1981). 107P. Jund,R. Jullien, Molecular Dynamics Calculation of the Thermal Conductivity of Vitreous Silica, Physical Review B 59, 13707 (1999). PAGE 253 253 108P. K. Schelling, S. R. Phillpot,P. Keblinski, Comparison of atomiclev el simulation methods for computing thermal conductivity, Physical Review B 65 (2002). 109K. Yamada, K. Kurosaki, M. Uno,S. Yamanaka, Evaluation of thermal properties of uranium dioxide by molecular dynamics, Journal Of Alloys And Compounds 307, 10 (2000). 110M. Abramowski, R. W. Grimes,S. Owens, Morphology of UO2, Journal of Nuclear Materials 275, 12 (1999). 111G. Busker, in Department of Materials (Imperial College of Science, Technology, and Medicine, London, UK, 2002), Vol. Ph.D. 112Y. Ida, Interionic Repulsive Force A nd Compressibility Of Ions, Physics Of The Earth And Planetary Interiors 13, 97 (1976). 113J. D. Gale, GULP: A computer program for the sy mmetryadapted simulation of solids, Journal of the Chemical SocietyFaraday Transactions 93, 629 (1997). 114J. D. Gale,A. L. Rohl, The General Utility Lattice Program (GULP), Molecular Simulation 29, 291 (2003). 115K. Govers, S. Lemehov,M. Verwerft, Comparison of Interatomic Potentials for UO2. Part I: Static Calculations, Journal of Nuclear Materials 266, 161 (2007). 116E. Wigner, On the Quantum Correction fo r Thermodynamic Equilibrium, Physical Review 40, 0749 (1932). 117G. Leibfried,E. Schloemann, Nachr. Akad. Wiss. Gottingen, Math. Physik. K1 IIa, 71 (1954). 118P. G. Klemens, Thermal Conductivity and Lattice Vibrational Modes, 1954). 119G. H. Geiger,D. R. Poirier, Transport Phenomena in Metallurgy (AddisonWesley, 1973). 120J. B. Wachtman, M. L. Wheat, H. J. Anderson,J. L. Bates, Elastic Constants of Single Crystal Uo2 at 25 Degrees C, Journal of Nuclear Materials 16, 39 (1965). PAGE 254 254 121I. J. Fritz, Elastic Properties of Uo2 at HighPressure, Journal of Applied Physics 47, 4353 (1976). 122G. Dolling, R. A. Cowley,A. D. B. Woods, Crystal Dynamics of Uranium Dioxide, Canadian Journal of Physics 43, 1397 (1965). 123H. Yoshida, K. Yokoyama, N. Shibata, Y. Ikuhara,T. Sakuma, Hightemperature grain boundary and grain boundary energy in cubic sliding behavior zirconia bicrystals, Acta Materialia 52, 2349 (2004). 124P. Keblinski, D. Wolf, S. R. Phillpot,H. Gleiter, Role of bonding and coordination in the atomic structure and energy of diamond and silicon grain boundaries, Journal of Materials Research 13, 2077 (1998). 125P. Shukla, T. Watanabe,S. R. Phillpot, Thermal Transport Properties of MgO and Nd2Zr2O7 Pyrochlore by Molecular Dynamics Simulation, unpublished (2007). 126G. A. Slack, Thermal Conductivity of MgO, Al2O3, Mgal2O4, and Fe3O4 Crystals from 3 Degrees to 300 Degrees K, Physical Review 126, 427 (1962). 127C. W. Nan,R. Birringer, Determining the Kapitza resistanc e and the thermal conductivity of polycrystals: A simple model, Physical Review B 57, 8264 (1998). 128J. Amrit, Grain boundary Kapitza resistance and grainarrangement induced anisotropy in the thermal conductivity of polycryst alline niobium at low temperatures, Journal of Physics DApplied Physics 39, 4472 (2006). 129R. J. Stevens, L. V. Zhigilei,P. M. Norris, Effects of temperature and disorder on thermal boundary conductance at solidsolid interfa ces: Nonequilibrium molecular dynamics simulations, International Journal of Heat and Mass Transfer 50, 3977 (2007). PAGE 255 255 130D. G. Cahill, W. K. Ford, K. E. Goodson, G. D. Mahan, A. Majumdar, H. J. Maris, R. Merlin,S. R. Phillpot, Nanoscale thermal transport, Journal of Applied Physics 93, 793 (2003). 131M. Freyss, T. Petit,J. P. Crocombette, Point Defects in Uranium Dioxide: ab initio Pseudopotential Approach in the Ge neralized Gradient Approximation, Journal of Nuclear Materials 347, 44 (2005). 132J. P. Crocombette, F. Jollet, L. N. Nga,T. Petit, Planewave Pseudopotential Study of Point Defects in Uranium Dioxide, Physical Review B 6410 (2001). 133B. T. M. Willis, Crystallographic Studies of AnionExcess UraniumOxides, Journal of the Chemical SocietyFaraday Transactions II 83, 1073 (1987). 134B. T. M. Willis, Structures of UO2, UO2+x and U4O9 by Neutron Diffraction, Solid State Communications 2, 23 (1964). 135B. T. M. Willis, Structures of UO2, UO2+x and U4O9 by Neutron Diffraction, Journal de Physique 25, 431 (1964). 136A. K. Cheetham, B. E. F. Fender,M. J. Cooper, Defect Structure of Calcium Fluoride Containing Excess Anions.1. Bragg Scattering, Journal of Physics C Solid State Physics 4, 3107 (1971). 137H. V. S. Hubbard,T. R. Griffiths, An Investigation of Defect St ructures in SingleCrystal UO2+x by OpticalAbsorption Spectroscopy, Journal of the Chemical SocietyFaraday Transactions II 83, 1215 (1987). 138G. C. Allen,P. A. Tempest, Ordered Defects in th e Oxides of Uranium, Proceedings of the Royal Society of London Series aMathemati cal Physical and Engineering Sciences 406, 325 (1986). PAGE 256 256 139G. C. Allen, P. A. Tempest,J. W. Tyler, Coordination Model For The Defect Structure Of Hyperstoichiometric UO2+X And U4O9, Nature 295, 48 (1982). 140H. Matzke, Atomic TransportProperties in UO2 and Mixed Oxides (U, Pu)O2, Journal of the Chemical SocietyFaraday Transactions II 83, 1121 (1987). 141P. Ruello, K. D. Becker, K. Ullrich, L. Desgranges, C. Petot,G. PetotErvas, Thermal Variation of the Optical Absorption of UO2: Determination of the Small Polaron Selfenergy, Journal of Nuclear Materials 328, 46 (2004). 142F. Jollet, T. Petit, S. Gota, N. Thromat, M. GautierSoyer,A. Pasturel, The Electronic Structure of Uranium Dioxide: an Oxygen Kedge Xray Absorption Study, Journal of PhysicsCondensed Matter 9, 9393 (1997). 143W. S. Capinski, H. J. Maris, E. Bauser, I. Silier, M. AsenPalmer, T. Ruf, M. Cardona,E. Gmelin, Thermal Conductivity of Isotopically Enriched Si, Applied Physics Letters 71, 2109 (1997). 144D. G. Cahill, F. Watanabe, A. Rockett,C. B. Vining, Thermal Conductivity of Epitaxial Layers of Dilute SiGe Alloys, Physical Review B 71 (2005). 145J. Janeczek, R. C. Ewing, V. M. Oversby,L. O. Werme, Uraninite and UO2 in Spent Nuclear Fuel: a Comparison, Journal of Nuclear Materials 238, 121 (1996). 146P. G. Lucuta, H. Matzke,R. A. Verrall, Modeling of UO2Based Simfuel ThermalConductivity the Effect of the Burnup, Journal of Nuclear Materials 217, 279 (1994). 147I. C. Hobson, R. Taylor,Ainscoug.Jb, Effect of Porosity and Stoichiometry on ThermalConductivity of UraniumDioxide, Journal of Physics DApplied Physics 7, 1003 (1974). 148C. Affortit,J. P. Marcon, High TemperaturesHigh Pressures 7, 236 (1970). PAGE 257 257 149P. G. Lucuta, H. Matzke,R. A. Verrall, ThermalConductivity of Hy perstoichiometric Simfuel, Journal of Nuclear Materials 223, 51 (1995). 150B. J. Lewis, B. Szpunar,F. C. Iglesias, Fuel Oxidation and Therma l Conductivity Model for Operating Defective Fuel Rods, Journal of Nuclear Materials 306, 30 (2002). 151S. Yamasaki, T. Arima, K. Idemitsu,Y. Inagaki, Evaluation of thermal conductivity of hyperstoichiometric UO2+x by molecular dynamics simulation, Internationa l Journal of Thermophysics 28, 661 (2007). 152F. Gronvold, HighTemperature XRay Study of Uranium Oxides in the UO2U3O8 Region, Journal of Inorganic & Nuclear Chemistry 1, 357 (1955). 153M. Amaya, T. Kubo,Y. Korei, Thermal Conductivity Measurements on UO2+x from 300 to 1,400 K, Journal of Nuclear Science and Technology 33, 636 (1996). 154J. W. Che, T. Cagin, W. Q. Deng,W. A. Goddard, Thermal Conductivity of Diamond and Related Materials from Molecular Dynamics Simulations, Journal of Chemical Physics 113, 6888 (2000). 155J. Li, L. Porter,S. Yip, Atomistic Modeling of Finite Temp erature Properties of Crystalline BetaSiC II. Therma l Conductivity and Effect s of Point Defects, Journal of Nuclear Materials 255, 139 (1998). 156S. G. Volz,G. Chen, Molecular Dynamics Simulation of Thermal Conductivity of Silicon Crystals, Physical Review B 61, 2651 (2000). 157J. D. Higgs, W. T. Thompson, B. J. Lewis,S. C. Vogel, Kinetics of Precipitation of U4O9 from Hyperstoichiometric UO2+x, Journal of Nuclear Materials 366, 297 (2007). 158R. J. Bell, P. Dean,Hibbinsb.D.C., Localization of Normal Modes in Vitreous Silica, Germania and Beryllium Fluoride, Journal of Physics Part C Solid State Physics 3, 2111 (1970). PAGE 258 258 159J. Callaway, Model for Lattice Thermal Conductivity at Low Temperatures, Physical Review 113, 1046 (1959). 160P. K. Schelling, S. R. Phillpot,D. Wolf, Mechanism of the CubictoTetragonal Phase Transition in Zirconia and Yttriastabilized Zirconia by Molecular Dynamics Simulation, Journal of the American Ceramic Society 84, 1609 (2001). 161K. Clausen, W. Hayes, M. T. Hutchings, J. E. Macdonald, R. Osborn,P. Schnabel, Investigation of Oxygen Disorder, Thermal Parameters, LatticeVibrations and ElasticConstants of Uo2 and Tho2 at Temperatures up to 2930K, Revue de Physique Appliquee 19, 719 (1984). 162K. Clausen, W. Hayes, J. E. M acdonald, R. Osborn,M. T. Hutchings, Observation of Oxygen Frenkel Disorder in UraniumDioxide above 200 0K by Use of NeutronScattering Techniques, Physical Review Letters 52, 1238 (1984). 163K. N. Clausen, M. A. Hackett, W. Hayes, S. Hull, M. T. Hutchings, J. E. Macdonald, K. A. Mcewen, R. Osborn,U. Steigenberger, Coherent Diffuse NeutronScattering from UO2 and ThO2 at Temperatures above 2000 K, Physica B 156, 103 (1989). 164K. C. Kim,D. R. Olander, Oxygen Diffusion in UO2X, Transactions of the American Nuclear Society 39, 382 (1981). 165R. Szwarc, Defect Contribution to E xcess Enthalpy of Uranium DioxideCalculation of Frenkel Energy, Journal of Physics a nd Chemistry of Solids 30, 705 (1969). 166T. Petit, C. Lemaignan, F. Jollet, B. Bigot,A. Pasturel, Point Defects in Uranium Dioxide, Philosophical Magazine BPhysics of Condensed Ma tter Statistical Mechanics Electronic Optical and Magnetic Properties 77, 779 (1998). PAGE 259 259 167M. Freyss, N. Vergnet,T. Petit, ab initio Modeling of the Be havior of Helium and Xenon in Actinide Dioxide Nuclear Fuels, Journal of Nuclear Materials 352, 144 (2006). 168W. M. Yao,C. Amsler,D. Asner,R. M. Barnett, J. Beringer,P. R. Burchat,C. D. Carone,C. Caso,O. Dahl,G. D'Ambrosio,A. De Gouvea,M. Do ser,S. Eidelman,J. L. Feng,T. Gherghetta,M. Goodman,C. Grab,D. E. Groom,A. Gurtu,K. Ha giwara,K. G. Hayes,J. J. HernandezRey,K. Hikasa,H. Jawahery,C. Kolda,Y. Kwon,M. L. Mangano,A. V. Manohar,A. Masoni,R. Miquel,K. Monig,H. Murayama,K. Nakamura,S. Navas,K. A. Olive,L. Pape,C. Patrignani,A. Piepke,G. Punzi,G. Raffelt,J. G. Smith,M. Tanabashi,J. Tern ing,N. A. Tornqvist,T. G. Trippe,P. Vogel,T. Watari,C. G. Wohl,R. L. Workman,P. A. Zyla ,B. Armstrong,G. Harper,V. S. Lugovsky,P. S. Schaffner,M. Artuso,K. S. Babu,H. R. Band,E. Barberio,M. Battaglia,H. Bichel,O. Biebel,P. Bloch,E. Blucher,R. N. Cahn,D. Casper,A Cattai,A. Ceccucci,D. Chakraborty,R. S. Chivukula,G. Cowan,T. Damour,T. De Grand,K. Desler,M. A. Dobbs,M. Drees,A. Edwards,D. A. Edwards,V. D. Elvira,J. Erler,W. Fetscher ,B. D. Fields,B. Foster ,D. Froidevaux,T. K. Gaisser,L. Garren,H. J. Gerber,G. Gerbier, L. Gibbons,F. J. Gilman,G. F. Giudice,A. V. Gritsan,M. Grunewald,H. E. Haber,C. Hagmann,I. Hinchliffe,A. Hocker,P. IgoKemenes,J. D. Jackson,K. F. Johnson,D. Karlen,B. Kayser,D Kirkby,S. R. Klein, K. Kleinknecht,J. G. Knowles,R. V. Kowalewski,P. Kreitz,B. Krusch e,Y. V. Kuyanov,O. Lahav,P. Langacker,A. Liddle,Z. Ligeti,T. M. Liss,L. Littenberg,J. C. Liu,K. S. Lugovsky,S. B. Lugovsky,T. Mannel,W. J. Marciano,A. D. Martin,D. Milstead,M. Narain,P Nason,Y. Nir,J. A. Peacock,S. A. Prell,A. Quadt,S. Raby,B. N. Ratcliff,E. A. Razuv aev,B. Renk,P. R. Richardson,S. Roesler,G. Rolandi,M. T. Roman,L. J. Rosenberg,C. T. Sachrajda,Y. Sakai,S. Sarkar,M. Schmitt,O. Schneider,D. Scott,T. Sjostrand,G. F. Smoot,P. Sokolsky,S. Spanier,H. Spieler,A. Stahl,T. Staney,R. E. Streitmatter,T. Sumiyoshi,N. P. Tkachenko,G. H. Trilling,G. Valencia,K. van PAGE 260 260 Bibber,M. G. Vincter,D. R. Ward,B. R. Webbe r,J. D. Wells,M. Whalley,L. Wolfenstein,J. Womersley,C. L. Woody,A. Yamamoto,O V. Zenin,J. Zhang,R. Y. Zhu, Review of particle physics, Journal of Physics GNuclear and Particle Physics 33, 1 (2006). 169G. F. Knoll, Radiation Detection and Measurement (Wiley, New York, 1989). 170O. B. Firsov, A Qualitative Interpretation of the Mean Electron Excitation Energy in Atomic Collisions, Soviet Physics JETPUSSR 9, 1076 (1959). 171J. Lindhard,M. Scharff, Energy Dissipation by Ions in Kev Region, Physical Review 1, 128 (1961). 172O. S. Oen,M. T. Robinson, Computer Studies of Reflectio n of LightIons from Solids, Nuclear Instruments & Methods 132, 647 (1976). 173H. Matzke, RadiationDamage Effects in NuclearMaterials, Nuclear Instruments & Methods in Physics Research Section BBeam In teractions with Materials and Atoms 32, 455 (1988). 174H. Matzke, RadiationDamage in Crystalline Insula tors, Oxides and Ceramic NuclearFuels, Radiation Effects and Defects in Solids 64, 3 (1982). 175K. Hayashi, H. Kikuchi,K. Fukuda, RadiationDamage of UO2 Implanted with 100Mev Iodine Ions, Journal of Alloys and Compounds 213, 351 (1994). 176H. J. Matzke, Radiation Effects in Fuel Materials for Fission Reactors, Institute of Physics Conference Series, 423 (1983). 177K. Minato, T. Shiratori, H. Serizawa, K. Hayashi, K. Une, K. Nogita, M. Hirai,M. Amaya, Thermal Conductivities of Irradiated UO2 and (U,Gd)O2, Journal of Nuclear Materials 288, 57 (2001). 178C. Ronchi,T. Wiss, Fission Fragment Spikes in Uranium Dioxide, Journal of Applied Physics 92, 5837 (2002). PAGE 261 261 179J. B. Gibson, A. N. Goland, M. Milgram,G. H. Vineyard, Dynamics of Radiation Damage, Physical Review 120, 1229 (1960). 180D. E. Harrison, Application Of MolecularDynamics Simulations To The Study Of IonBombarded MetalSurfaces, Crc Critical Reviews In Solid State And Materials Sciences 14, S1 (1988). 181B. J. Garrison, N. Winograd, D. M. Deaven, C. T. Reimann, D. Y. Lo, T. A. Tombrello, D. E. Harrison,M. H. Shapiro, ManyBody EmbeddedAtom Potentia l for Describing the Energy and AngularDistributions of Rh Atoms De sorbed from IonBombarded Rh(111), Physical Review B 37, 7197 (1988). 182D. Y. Lo, T. A. Tombrello, M. H. Shapiro, B. J. Garrison, N. Winograd,D. E. Harrison, TheoreticalStudies of IonBom bardment ManyBody Interactions, Journal of Vacuum Science & Technology aVacuum Surfaces and Films 6, 708 (1988). 183M. H. Shapiro, T. A. Tombrello,D. E. Harrison, Simulation of Isotop ic Mass Effects in Sputtering.2., Nuclear Instruments & Methods in Physic s Research Section BBeam Interactions with Materials and Atoms 30, 152 (1988). 184K. Morishita, R. Sugano,B. D. Wirth, Thermal Stability of He liumVacancy Clusters and Bubble Formation Multiscale Modeling A pproach for Fusion Materials development, Fusion Science and Technology 44, 441 (2003). 185R. Smith, D. Bacorisen, B. P. Uberuaga, K. E. Sickafus, J. A. Ball,R. W. Grimes, Dynamical Simulations of Radiation Damage in Magnesium Alumi nate Spinel, MgAl2O4, Journal of PhysicsCondensed Matter 17, 875 (2005). PAGE 262 262 186B. P. Uberuaga, R. Smith, A. R. Cleave, G. Henkelman, R. W. Grimes, A. F. Voter,K. E. Sickafus, Dynamical Simulations of Radiation Damage and Defect Mobility in MgO, Physical Review B 71 (2005). 187R. Devanathan, W. J. Weber,T. D. de la Rubia, Computer Simulation of a 10 keV Si Displacement Cascade in SiC, Nuclear Instruments & Methods In Physics Research Section BBeam Interactions With Materials And Atoms 141, 118 (1998). 188T. D. delaRubia, Irradiationinduced defect pr oduction in elemental metals and semiconductors: A review of r ecent molecular dynamics studies, Annual Review Of Materials Science 26, 613 (1996). 189R. Smith,K. Beardmore, Molecular Dynamics Studies of Par ticle Impacts w ith Carbonbased Baterials, Thin Solid Films 272, 255 (1996). 190L. Van Brutzel, J. M. Delaye, D. Ghaleb,M. Rarivomanantsoa, Molecular Dynamics Studies of Displacement Cascades in the Uranium Dioxide Matrix, Philosophical Magazine 83, 4083 (2003). 191N. D. Morelon, D. Ghaleb, J. M. Delaye,L. Van Brutzel, A New Empirical Potential for Simulating the Formation of Defects and Their Mobility in Uranium Dioxide, Philosophical Magazine 83, 1533 (2003). 192L. Van Brutzel, M. Rarivomanantsoa,D. Ghaleb, Displacement Cascade Initiated with the Realistic Energy of the Recoil Nucleus in UO2 Matrix by Molecular Dynamics Simulation, Journal of Nuclear Materials 354, 28 (2006). 193M. J. Norgett, M. T. Robinson,I. M. Torrens, Proposed Method of Calculating Displacement DoseRates, Nuclear Engineering and Design 33, 50 (1975). PAGE 263 263 194J. F. Ziegler, J. P. Biersak,U. Littmark, The Stopping and Range of Ions in Matter (Pergomon, New York, 1985). 195R. Smith, Atomic & Ion Collisions in Solids and at Surfaces (Cambridge University Press, Cambridge, UK, 1997). 196J. Fikar,R. Schaeublin, Molecular Dynamics Simulation of Radiation Damage in bcc Tungsten, Nuclear Instruments & Methods in Physic s Research Section BBeam Interactions with Materials and Atoms 255, 27 (2007). 197J. M. Delaye,D. Ghaleb, Dynamic processes during displacemen t cascades in oxide glasses: A moleculardynamics study, Physical Review B 61, 14481 (2000). 198J. M. Delaye,D. Ghaleb, Dynamic Processes during Displacement Cascades in Oxide Glasses: a Molecular Dynamics Study, Physical Review B 61, 14481 (2000). 199L. Van Brutzel,M. Rarivomanantsoa, Molecular Dynamics Simulation Study of Primary Damage in UO2 Produced by Cascade Overlaps, Journal of Nu clear Materials 358, 209 (2006). 200L. R. Corrales, A. Chartier,R. Devanathan, Excess Kinetic Energy Di ssipation in Materials, Nuclear Instruments & Methods in Physics Re search Section BBeam Interactions with Materials and Atoms 228, 274 (2005). 201A. L. Garcia, Numerical Methods for Physics (PrenticeHall, Inc., 2000). 202Y. H. Hu,S. B. Sinnott, Constant Temperature Molecular Dy namics Simulations of Energetic Particlesolid Collisions: Comparison of Temperature control methods, Journal of Computational Physics 200, 251 (2004). 203F. Reif, Fundamentals of Statisti cal and Thermal Physics (McGrawHill, New York, 1965). 204K. Huang, Statistical Mechanics (John Wiley & Sons Inc., New York, 1987). PAGE 264 264 205H. J. C. Berendsen, J. P. M. Postma, W. F. Vangunsteren, A. Dinola,J. R. Haak, MolecularDynamics with Coupling to an External Bath, Journal of Chemical Physics 81, 3684 (1984). 206S. A. Adelman,J. D. Doll, Generalized Langevin Equation A pproach for AtomSolidSurface Scattering General Formulation for Cl assical Scattering Off Harmonic Solids, Journal of Chemical Physics 64, 2375 (1976). 207K. E. Sickafus, L. Minervini, R. W. Grimes J. A. Valdez, M. Ishimaru, F. Li, K. J. McClellan,T. Hartmann, Radiation Tolerance of Complex Oxides, Science 289, 748 (2000). 208H. O. Pierson, Handbook of Carbon, Graphite, Diam ond and Fullerenes Properties, Processing and Applications (William Andrew Publishing/Noyes, 1993). 209O. Auciello, J. Birrell, J. A. Carlisle, J. E. Gerbi, X. C. Xiao, B. Peng,H. D. Espinosa, Materials Science and Fabrication Proce sses for a New MEMS Technoloey Based on Ultrananocrystalline Diamond Thin Films, Journal af PhysicsCondensed Matter 16, R539 (2004). 210M. H. Nazare,A. J. Neves, Properties, Growth and A pplications of Diamond (Institution of Engineering and Technology, 2001). 211T. Sato, K. Ohashi, T. Sudoh, K. Haruna,H. Maeta, Thermal Expansion of a High Purity Synthetic Diamond Single Crys tal at Low Temperatures, Physical Review B 65 (2002). 212G. A. Slack,S. F. Bartram, ThermalExpansion of Som e DiamondLike Crystals, Journal of Applied Physics 46, 89 (1975). 213P. K. Schelling, S. R. Phillpot,P. Keblinski, Comparison of atomiclev el simulation methods for computing thermal conductivity, Physical Review B 65, 144306 (2002). PAGE 265 265 214D. G. Onn, A. Witek, Y. Z. Qi u, T. R. Anthony,W. F. Banholzer, Some Aspects of the ThermalConductivity of Isotopicall y Enriched Diamond SingleCrystals, Physical Review Letters 68, 2806 (1992). 215S. von Alfthan, P. D. Haynes, K. Kaski,A. P. Sutton, Are the structures of twist grain boundaries in silicon ordered at 0 K? Physical Review Letters 96, 055505 (2006). 216P. Shukla, T. Watanabe, J. C. Nino, J. S. Tulenko,S. R. Phillpot, Thermal Transport Properties of MgO and Nd2Zr2O7 Pyrochlore by Molecular Dynamics Simulation, Journal of Nuclear Materials (2007). 217H. K. Lyeo,D. G. Cahill, Thermal conductance of interfa ces between highly dissimilar materials, Physical Review B 73, 144301 (2006). 218T. Watanabe, B. Ni, S. R. Phillpot P. K. Schelling,P. Keblinski, Thermal Conductance Across Grain Boundaries in Diamond from Molecular Dynamics Simulation, Journal of Applied Physics 106, 063503 (2007). 219M. A. Angadi, T. Watanabe, A. Bodapati, X. C. Xiao, O. Auciello, J. A. Carlisle, J. A. Eastman, P. Keblinski, P. K. Schelling,S. R. Phillpot, Thermal transport and grain boundary conductance in ultrananocrystalline diamond thin films, Journal of Applied Physics 99 (2006). 220R. Berman, in The Properties of Diamond, edited by J. E. Field (Academic Press, London, 1979). 221J. Birrell, J. E. Gerbi, O. Auciello, J. M. Gibson, J. Johnson,J. A. Carlisle, Interpretation of the Raman spectra of ultrananocrystalline diamond, Diamond and Related Materials 14, 86 (2005). 222A. V. Sumant, A. R. Krauss, D. M. Gruen, O. Auciello, A. Erdemir, M. Williams, A. F. Artiles,W. Adams, Ultrananocrystalline diamond film as a wearresistant and protective coating for mechanical seal applications, Tribology Transactions 48, 24 (2005). PAGE 266 266 223S. Bhattacharyya, O. Auciello, J. Birrell, J. A. Ca rlisle, L. A. Curtiss, A. N. Goyette, D. M. Gruen, A. R. Krauss, J. Sc hlueter, A. Sumant,P. Zapol, Synthesis and Characte rization of Highly Conducting Nitrogendoped Ultrana nocrystalline Diamond Films, Applied Physics Letters 79, 1441 (2001). 224X. C. Xiao, J. Wang, L. Chao, J. A. Carlisle, B. Mech, R. Greenberg, D. Guven, R. Freda, M. Humayun, J. Weiland,O. Auciello, In Vitro and In Vivo Evaluation of Ultrananocrystalline Diamond for Coating of Impl antable Retinal Microchips, Journal of Biomedical Materials Research Part B: Applied Biomaterials 77, 273 (2006). 225D. G. Cahill, K. Goodson,A. Majumdar, Thermometry and thermal transport in micro/nanoscale solidstate devices and structures, Journal of Heat TransferTransactions of the Asme 124, 223 (2002). 226J. Birrell, J. E. Gerbi, O. Auciello, J. M. Gibson, D. M. Gruen,J. A. Carlisle, Bonding Structure in Nitrogen Doped Ultrananocrystalline Diamond, Journal of Applied Physics 93, 5606 (2003). 227S. Ertl, M. Adamschik, P. Schmid P. Gluche, A. Floter,E. Kohn, Surface micromachined diamond microswitch, Diamond and Related Materials 9, 970 (2000). 228D. G. Cahill,R. O. Pohl, ThermalConductivity of Amorphous Solids above the Plateau, Physical Review B 35, 4067 (1987). 229D. G. Cahill, M. Katiyar,J. R. Abelson, ThermalConductivity of AlphaSih ThinFilms, Physical Review B 50, 6077 (1994). 230A. V. Sumant, D. S. Grierson, J. E. Gerbi, J. Birrell, U. D. Lanke, O. Auciello, J. A. Carlisle,R. W. Carpick, Toward the ultimate tribological interface: Surface chemistry and nanotribology of ultr ananocrystalline diamond, Advanced Materials 17, 1039 (2005). PAGE 267 267 231C. J. Morath, H. J. Maris, J. J. Cuomo, D. L. Pappas, A. Grill, V. V. Patel, J. P. Doyle,K. L. Saenger, Picosecond Optical Studi es of Amorphous Diamond and DiamondLike Carbon ThermalConductivity and Longitudinal SoundVelocity, Journal of Applied Physics 76, 2636 (1994). 232A. J. Bullen, K. E. O'Hara, D. G. Cahill, O. Monteiro,A. von Keudell, Thermal conductivity of amorphous carbon thin films, Journal of Applied Physics 88, 6317 (2000). 233J. E. Graebner, S. Jin, G. W. Kamm lott, J. A. Herb,C. F. Gardinier, Unusually High ThermalConductivity in Diamond Films, Applied Physics Letters 60, 1576 (1992). 234D. T. Morelli, C. P. Beetz,T. A. Perry, ThermalConductivity of Synthetic Diamond Films, Journal of Applied Physics 64, 3063 (1988). 235J. Philip, P. Hess, T. Feygelson, J. E. Butle r, S. Chattopadhyay, K. H. Chen,L. C. Chen, Elastic, mechanical, and thermal proper ties of nanocrystalline diamond films, Journal of Applied Physics 93, 2164 (2003). 236V. B. Efimov,L. P. MezhovDeglin, Phonon scattering in diamond films, Physica B 263, 745 (1999). 237J. E. Graebner,S. Jin, J ournal of Applied Physics 76, 3 (1994). 238J. E. Graebner, M. E. Reiss, L. Seibles, T. M. Hartnett, R. P. Miller,C. J. Robinson, Physical Review B 50, 6 (1994). 239J. E. Graebner, V. G. Ralchenko, A. A. Smolin E. D. Obraztsova, K. G. Korotushenko,V. I. Konov, Diamond and Related Materials 5, 693 (1996). 240O. W. Kading, E. Matthias, R. Zachai, H.J. Fusser,P. Munzinger, Diamond and Related Materials 2, 1185 (1993). 241K. E. Goodson, O. W. Kading,R. Zachai, ASME HTD Proceedings 292, 83 (1994). PAGE 268 268 242K. Plamann, D. Fournier, E. Anger,A. Gicquel, Photothermal examina tion of the heat diffusion inhomgogeneity in diam ond films of submicron thickness, Diamond and Related Materials 3, 752 (1994). 243D. Fournier,K. Plamann, Di amond and Related Materials 4, 809 (1995). 244J. Hartmann, P. Voigt,M. Reichling, Measuring local thermal conduc tivity in po lycrystalline diamond with a high resolution photothermal microscope, Journal of Applied Physics 81, 2966 (1997). PAGE 269 269 BIOGRAPHICAL SKETCH Taku Watanabe was born in Seiro, Niigata, Ja pan. He received the B.E. degree from the Department of Physics, Stevens Institute of Technology in Hoboken, NJ in 2000 and the M.S. degree from the Department of Phys ics, University of Florida in Gainesville, FL, in 2004. He then began his Ph.D. study in the Department of Materials Science and Engineering with Prof. Simon R. Phillpot as his advisor. 