Thermal Transport at Dislocations and Interfaces by Atomistic Simulation

Material Information

Thermal Transport at Dislocations and Interfaces by Atomistic Simulation
Deng, Bowen
Place of Publication:
[Gainesville, Fla.]
University of Florida
Publication Date:
Physical Description:
1 online resource (119 p.)

Thesis/Dissertation Information

Doctorate ( Ph.D.)
Degree Grantor:
University of Florida
Degree Disciplines:
Materials Science and Engineering
Committee Chair:
Committee Co-Chair:
Committee Members:
Graduation Date:


Subjects / Keywords:
Atoms ( jstor )
Conceptual lattices ( jstor )
Grain boundaries ( jstor )
Kapitza resistance ( jstor )
Phonons ( jstor )
Silicon ( jstor )
Simulations ( jstor )
Temperature dependence ( jstor )
Thermal conductivity ( jstor )
Wave packets ( jstor )
Materials Science and Engineering -- Dissertations, Academic -- UF
defects -- interface -- phonon -- thermal -- transport -- uo2
bibliography ( marcgt )
theses ( marcgt )
government publication (state, provincial, terriorial, dependent) ( marcgt )
born-digital ( sobekcm )
Electronic Thesis or Dissertation
Materials Science and Engineering thesis, Ph.D.


The object of this thesis is to using atomic simulation technique to investigate thermal transport in selected crystalline materials, especially to understand how defect structures affect thermal conduction. This task is accomplished primarily using the tools of molecular dynamics (MD) simulations and the phonon wave packet dynamics (PWD) method. Non-equilibrium molecular-dynamics is used to characterize the effects of dislocations on the thermal transport properties of UO2. The temperature dependence of the effect is obtained. The results are also compared with the classic Klemens-Callaway analysis. The effect of dislocation density is also quantified. The simulation results are fit to the pertinent part of the empirical formula for the thermal conductivity used in the FRAPCON fuel-performance code, which gives the overall effects of temperature and dislocation effects on thermal conductivity. The fitted results can be well-described within this formalism, indicating that the results of molecular-dynamics simulations can be used as a reliable source of parameters for models at longer length scales. A phonon wave packet dynamics method is used to in addition to the MD method to characterize the Kapitza resistance of twist grain boundaries in UO2 and the Si/SiO2 interface in a Si/SiO2/Si heterostructure. The wave-packet dynamics method provides additional detailed insights on how phonon scatters at the grain boundary and interface. The energy transmission coefficients for individual phonon mode is quantified. Also, an elementary chain model is used to simulate and illustrate the phonon transport at a grain boundary. Connections to the results of phonon wave packet dynamics method are produced. ( en )
General Note:
In the series University of Florida Digital Collections.
General Note:
Includes vita.
Includes bibliographical references.
Source of Description:
Description based on online resource; title from PDF title page.
Source of Description:
This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
Thesis (Ph.D.)--University of Florida, 2014.
Co-adviser: SINNOTT,SUSAN B.
Statement of Responsibility:
by Bowen Deng.

Record Information

Source Institution:
Rights Management:
Applicable rights reserved.
LD1780 2014 ( lcc )


This item has the following downloads:

Full Text




© 2014 Bowen Deng


To my wife and parents


4 ACKNOWLEDGMENTS research and life. I am always inspired by his ideas, and all of which are always accompanied by his sense of humor. Besides, his intelligence a nd vision in work always has a positive effect in my character. I am also thankful to Dr. Sinnott and Dr. Chernatynskiy. Dr. Sinnott has provided many valuable comments to my research during all the research meetings and Dr. Chernatynskiy is always there t o help me and everyone in our group. I am always amazed by his various ideas on a subject, from macro to micro. I would also like to extend my thanks to my committee members: Dr. Sinnott, Dr. Nino, Dr. Man uel, Dr. Tulenko and Dr. Taylor . I would also thank Dr. Jones and Dr. Webb, who have been very gracious to attend my defense. Since I first came to University of Florida, it has been my pleasure to work in the Computational Materials Focus Group, where I was always surrounded by many great people. I would like to thank all the previous and current members in the group. Finally, I would like to t hank my parents, who raised me up with devotion and love. It is their support and love that lead me this far. And my wife, who has make the last two and half years a very enjoyable time and has been very supportive during my dissertation writing.


5 TABLE OF CONTENTS page ACKNOWLEDGMENTS ................................ ................................ ................................ .. 4 LIST OF TABLES ................................ ................................ ................................ ............ 7 LIST OF FIGURES ................................ ................................ ................................ .......... 8 ABSTRACT ................................ ................................ ................................ ................... 11 CHAPTER 1 GENERAL INTRODUCTION ................................ ................................ .................. 13 Introduction to Heat Transfer in Crystalline Materials ................................ ............. 13 Materials of Interest: Silicon and Silica ................................ ................................ ... 15 Materials Defect Structure ................................ ................................ ...................... 16 Sc ope of This Study ................................ ................................ ................................ 18 2 THERMAL TRANSPORT IN DIELECTRIC MATERIALS ................................ ........ 22 Introduction to Thermal Conduction ................................ ................................ ........ 22 Phonon mediated Thermal Transport ................................ ................................ ..... 22 3 METHODOLOGY ................................ ................................ ................................ ... 30 Molecular Dynamics Simulation of Thermal Transport ................................ ........... 30 Anharmonic Correction ................................ ................................ ........................... 33 Non Equilibrium Molecular Dynamics ................................ ................................ ..... 34 Phonon Wave Packet Dynamics ................................ ................................ ............. 35 Calculating boundary conductance ................................ ................................ ......... 38 Computational Tool ................................ ................................ ................................ . 39 4 EFFECTS OF EDGE DISLOCATIONS ON THERMAL TRANSPORT IN UO 2 ....... 44 Introduction ................................ ................................ ................................ ............. 44 Simulation Setup ................................ ................................ ................................ ..... 48 Dependence of Thermal Conductivity on Temperature and Dislocation Density .... 49 Comparison of Various Dislocatio n Systems ................................ .......................... 51 Theoretical Model and Fuel Performance Code ................................ ..................... 51 Conclusions ................................ ................................ ................................ ............ 53 5 KAPITZA RESIST ANCE OF Si/SiO 2 INTERFACE ................................ ................. 63 Introduction ................................ ................................ ................................ ............. 63 Simulation Setup ................................ ................................ ................................ ..... 65


6 Phonon Scatterings at Interface ................................ ................................ .............. 67 Calculating Kapitza Resistance ................................ ................................ .............. 70 Conclusion ................................ ................................ ................................ .............. 72 6 THERMAL TRANSPORT AT A TWIST GRAIN BOUNDARY IN UO 2 ..................... 81 Introduction ................................ ................................ ................................ ............. 81 Kapitza Resistance by NEMD Method ................................ ................................ .... 83 Phonon Transmission by PWD Method ................................ ................................ .. 84 Conclusion ................................ ................................ ................................ .............. 88 7 LATTICE DYNMICS MODEL FOR INTERFACIAL THERMAL TRANSPORT ........ 96 Introduction ................................ ................................ ................................ ............. 96 One Dimensional Monatomic Chain ................................ ................................ ....... 98 One Dimensional Diatomic Chain ................................ ................................ ......... 100 Conclusions ................................ ................................ ................................ .......... 101 8 CONCLUSION AND OUTLOOK ................................ ................................ ........... 107 LI ST OF REFERENCES ................................ ................................ ............................. 111 BIOGRAPHICAL SKETCH ................................ ................................ .......................... 119


7 LIST OF TABLES Table page 1 1 Distinct sites in UO 2 crystal structure. ................................ ................................ . 19 3 1 Parameters for Busker potential. ................................ ................................ ........ 19 3 2 Parameters for extended SW potential. ................................ .............................. 41 4 1 Fitting results for different dislocation density models ................................ ........ 56 4 2 ................................ ................................ . 56


8 LIST OF FIGURES Figure page 1 1 A unit cell of fluorite structured UO 2 .. ................................ ................................ . 20 1 2 A unit cell of diamond structured crystalline Si. ................................ .................. 20 1 3 Example crystalline SiO 2 structures. ................................ ................................ ... 21 1 4 Illustration of edge and screw dislocations ................................ ......................... 21 1 5 Illustration of Umk lapp process and normal process . ................................ ......... 29 3 1 Phonon dispersion relation of silicon along [001]. ................................ ............... 41 3 2 Schematic of simulation box for calculating thermal conductivity using non equilibrium molecular dynamics. ................................ ................................ ......... 42 3 3 Illustration of a scattering event.. ................................ ................................ ........ 42 3 4 First Brillouin Zone sampling points of a diamond lattice in k z =0 plane. ............. 43 4 1 Three cases for dislocation core structures. ................................ ....................... 55 4 2 Increase in d islocation density with burnup and fits from the MFPR code .......... 56 4 3 Anharmonic correction of MD results. ................................ ................................ . 56 4 4 Change of thermal conductivity between d=0 and d= 1.68X10 16 m 2 . .................. 57 4 5 Temperature dependence of thermal conductivity of different dislocation density structures. ................................ ................................ .............................. 58 4 6 Dislocation dependence of resistivity at various temperatures. .......................... 59 4 7 Comparison of the thermal conductivity of {110}<110> and {100}<110> ............ 60 4 8 Comparison of dependence of thermal conductivity on dislocation density on with the Klemens Callaw ay model ................................ ................................ ...... 61 4 9 Comparison of temperature dependence of thermal conductivity of Klemens Callaway model and fitting ................................ ................................ .................. 62 5 1 Simulation setup for simulations . ................................ ................................ ........ 74 5 2 Transmission coefficients for LA and LO branches as a function of incident phonon frequency. ................................ ................................ .............................. 74


9 5 3 Transmission coefficients for all modes along k 00 direction for l sio2 =0.8 nm structure ................................ ................................ ................................ ............. 75 5 4 LA incident phonon transmission for l SiO2 =0.8 nm a nd l SiO2 =10.9 nm ................. 75 5 5 Relative transmitted energy as a function of k x and k y for four different incident phonon frequencies. ................................ ................................ .............. 76 5 6 Wavevector deviation as a function of incident phonon frequency for transmission and reflection phonons. ................................ ................................ . 76 5 7 Frequency distribution of wave packets for l SiO2 =0.8 nm structure. .................... 77 5 8 Contribution to the thermal conductivity as a function of k x and k y , with summation of all k z contribution. ................................ ................................ ......... 78 5 9 Extracting the Kapitza resistance of the Si/SiO 2 interface at 300 K. ................... 79 5 10 Contribution to Kapitza conductance by mode for l SiO2 =0.8 nm. ......................... 80 6 1 Schematic of a twist grain boundary. ................................ ................................ .. 89 6 2 Schematic illustration of determining t he rotation angle for CSL of 17 twist grain boundary ................................ ................................ ................................ .... 89 6 3 CSL distributi on of uranium dioxide sample. ................................ ....................... 90 6 4 Simulation configuration and temperature profile for calculating Kapitza resistance.. ................................ ................................ ................................ ......... 90 6 5 Temperature boundaries in UO 2 . ................................ ................................ ............................. 91 6 6 Avoided crossing for phonon dispersion relation in UO 2 .. ................................ ... 91 6 7 boundary as a function of phonon frequency.. ................................ .................... 92 6 8 Time evolution of atomic displacement with LA incident phonon wave packet ................................ ................................ ............... 93 6 9 Time evolution of atomic displacement with LO incident phonon wave packet ................................ ................................ ................ 94 6 10 Contribution of the un scattered mode to the energy transmission coefficient.. . 95 6 11 Phonon energy contributi on from inelastic scattering. ................................ ........ 95 7 1 Schematic of the model setup. ................................ ................................ ......... 103


10 7 2 Flow chart for solving LD model. ................................ ................................ ...... 103 7 3 One dimensional monoatomic chain model. ................................ ..................... 104 7 4 Energy transmission as a function of f 3 for the monatomic chain. .................... 104 7 5 Energy transmission for various extended interfaces. ................................ ...... 105 7 6 Illustration of diatomic model.. ................................ ................................ .......... 105 7 7 Phonon dispersion relation for the diatomic model. ................................ .......... 106 7 8 Energy transmission of diatomic model as a function of i ncident phonon frequency . ................................ ................................ ................................ ......... 106


11 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy THERMAL TRANSPORT AT DISLOCATION S AND INTERFACES BY ATOMISTIC SIMULATION By Bowen Deng Aug ust 2014 Chair: Simon R. Phillpot Major: Materials Science and Engineering The object of this thesis is to use atomic simulation technique s to investigate thermal transport in selected crystalline materials, especially to understand how defect structures affect thermal conduction. This task is accomplished primarily using the tools of molecular dynamics (MD) simulations and phonon wave packet dynamics (PWD). Non equilibrium m olecu lar dynamics is used to characterize the effects of dislocations on the thermal transport properties of UO 2 . The temperature dependence of the effect is obtained. The results are also compared with the classic Klemen s Callaway analysis. The effect of dislocation density is also quantified. The simulation results are fit to the pertinent part of the empirical formula for the thermal conductivity in the FRAPCON fuel performance code, which gives the overall effects of temperature and dislocation effects on thermal conductivity. The fitted results can be well described within this formalism, indicating that the results of molecular dynamics simulations can be used as a reliable source of parameters for models at longer l ength scales. A phonon wave packet dynamics method is used to in addition to the MD method to character ize the Kapitza resistance of twist grain boundaries in UO 2 and the Si/SiO 2


12 interface in a Si/SiO 2 /Si heterostructure. The wave packet dynamics method pr ovides additional detailed insights as to how phonon s scatter at a grain boundary and at an interface. The energy transmission coefficients for individual phonon mode is quantified. Also, an elementary chain model is used to simulate and illustrate phonon transport at a grain boundary. Connections to the results of phonon wave packet dynamics method are discussed .


13 CHAPTER 1 GENERAL INTRODUCTION Introduction to Heat Transfer in Crystalline Materials Heat transfer plays a critical role in our daily life, as well as in research in materials science and engineering. There are numerous applications involving heat transfer, which can be broadly categorized into two categories: insulation and dissipation. I n order to better control or utilize heat, it is crucial to understand the factors that c an affect heat transfer . In solid crystalline materials, a fundamental factor that regulate s heat transfer is the presence of defects. In order to better manage heat t ransfer, a good understanding of how individual defects affect the heat flow is needed. Crystallographic defects, including point defects, line defects, planar defects and bulk defects, act as obstacles to the heat carriers. Such structural defects general ly have at least one nanoscale dimension , or even atomic scale dimension , such as a lattice vacancy. Thus, a better understanding of the mechanisms and effects of the inhibiting processes at the atomic scale is a key building block to controlling heat tran sfer. Numerous experimental efforts have been made to understand the se mechanisms. However, the task remains very challenging due to the complex nature of real world materials structure s , especially in the presence of abundant defects. The heat carriers c an interact individually with each defect and any combinations of defects, thus making it difficult to separate each mechanism and effect. Also, certain processes such as the operation of a nuclear fuel pellet inside a nuclear reactor take place under extr eme environmental conditions, creating even more hurdles for research.


14 With the constant development of computational hardware [ 1 ] and algorithms [ 2 ] , computer simulation has become an indispensable contributor to the study of heat transfer problems in that it is able to focus on the interaction of only one defect with the heat carrier, observe phenomen a in sub nano time and space scale, and completely control environmental variables. Materials of Interest: Uranium Dioxide Uranium dioxide (UO 2 ) is ubiquitous as a nuclear fuel, and is used in all US commercial nuclear power reactors. [ 3 ] Among its advantages compared to metallic fuels are its high melting point and the fact that the fluorite structure is open enough to accommodate fission products without unacceptable levels of swelling [ 4 ] . However, one of the drawbacks of UO 2 fuel is its low thermal conductivity (~4 Wm 1 K 1 at 800 K [ 5 ] ) , which results in a large amount of thermal energy being stored in the pellet. Moreover, as the microstructure evolves during burn up, the thermal conductivity can further decrease by the increased scattering of phonons, the dominant heat carriers, from various microstructures elements. This leads not only to an overall decrease in the already low thermal conductivity, but also can lead to the development of relatively hotter and colder spots in the fuel, which can be expected to lead to further local variations in microstructure and further temperature variations (800 K over 0.5 cm [ 6 ] ) . UO 2 exhibits the fluorite type structure which belongs to the face cubic crystal system ( ). Its lattice constant is 5.4704 Ã… at 300 K [ 7 ] ; the unit cell is illustrated in Figure 1 1 . The uranium atoms are arranged FCC positions, while the oxygen atoms exhibit a simple cubic lattice; their respective sites is summarized in Table 1 1 .


15 Materials of Interest: Silicon and Silica Silicon is a well known and very widely used semiconducting material ; silicon di oxide , silica, is a very good dielectric material which is usually u sed as a substrate or protective layer in combination of silicon. For example, i n the silicon on insulator (SOI) technology in microelectronics, the SOI wafer structure can be made up of a thin layer of SiO 2 , functioning as a dielectrics , which prevents cu rrent leakage and is essential for the device to work. [ 8 ] Silicon has a very good thermal conductivity of 149 Wm 1 K 1 at 300 K [ 9 ] , which is comparable to some metals. For example , at 300K, the thermal conductivity of pure Fe is 80 Wm 1 K 1 ; W is 178 Wm 1 K 1 ; Al is 250 Wm 1 K 1 ; and Cu is 390 Wm 1 K 1 . [ 9 ] However, the overall cooling of the electronics by thermal conduction can be limited by both the introduction of the low thermal conduction oxide layer and the Si/SiO 2 interfaces themselves. Pure crystalline silicon has a diamond structure ( , space group of No.227, ). Its lattice constant is 5.4307 Ã… at 300 K [ 10 ] , with each silicon atom bonded to another four silicon atoms loc ated at the vertexes of a tetrahedron, as shown in Figure 1 2 . Silica, formally known as silicon dioxide, has various crystalline structures due to the flexibility of the O Si O bond , as illustrated in Figure 1 3 . Generally, the silicon atom is bonded to 4 nearby oxygen atoms, form ing a tetrahedral coordinat ed structure . A morphous silica also has the similar bonding since the big difference between crystalline and amorphous is only the absence of the long ranged order , with the coordination number of the Si atom being 4 and that of the O atom being 2 , with the


16 modification that, in reality, there are always over coordinated and sub coordinated atoms. The structural form used i cristobalite, because it can form a minimal lattice mismatch with silicon under finite dimension. It also has the same space group as crystal line silicon but with lattice constant of 7.16 Ã… at 300 K [ 11 ] , as shown in Figure 1 3 ( a ). Materials Defect Structure T here are always defects or imperfections in real crystalline materials. The defects, due to their inconsistency with the ideal lattice, break the local crystalline symmetry. One scheme to classify the defects is according to the dimensionality of the defect . Point defects are zero dimensional defects, hav ing no extension in space. There are two main kinds of intrinsic point defect s: vacancies and interstitials. As indicated by the ir names, vacancy defect refers to a vacant lattice site that would be otherwise occupied in an ideal crystal; an interstitial defect denotes an occupied position that is normally vacant in an ideal crystal. In an ionically bonded crystal (ionic crystal), ch arges are generally associated with each atom. In order to compensate the net charge created by an individual point defect, other point defects are normally present in order to maintain charge neutrality . There are two such charge neutral defect complexes : Schottky defect and Frenkel defect. A Schottky defect is a set of vacancies of zero net charge. However, this may or may not preserve the stoichiometry of structure. For example in SrTiO3, the Sr 2+ and O 2 pseudo Schottkys defect is charge neutral but do not preserve the stoichiometry of the system. A Frenkel defect, or Frenkel pair, is an interstitial defect and the vacancy it leaves behind.


17 However, there can also be cases where the ion defects exist in a non compensating manner, for instance, an off stoichiometr ic ion can co exist with an electron , such as in CeO 2 x . In metals and sometimes in ionics with two types of cations or anions, anti site defect s can be formed by two cations or two anions of different species interchang ing site s . E xtrinsic point defects involve impurity atoms as substitutional atoms or as interstitials . A dislocation is a one dimensional defect and it runs into a crystal in straight or curved paths . T he crystallographic alignment of a bulk structure terminates at the dislocation line. Dislocation defects can be classified into two primary categories: edge dislocations and screw dislocations. T he sketches in Figure 1 4 show the edge dislocation and screw dislocation in a simple cubic lattice. In both cases, the DC line indicates the dislocation line, which runs through the crystal. An e dge dislocation essentially contains one or more ms. As illustrated in the figure, only the half structure above the DC line is displaced along the direction of x axis , which is perpendicular to the DC line ; this mark s of the boundary between displaced and undisplaced crystal is the dislocation line. Sim ilarly, for screw dislocations, the dislocation line is also denoted as DC; however, the slip direction is now along the direction of z axis which is parallel to the dislocation line. These two categories of dislocations can exist individually, but in real materials they more coexist. For example a dislocation loop , which is a closed loop formed by dislocation lines , may be formed by one type of dislocation, or two types of dislocations (both edge, both screw or one of each) with different Burgers vectors [ 12 ] .


18 Planar defects are generally the exterior surface of a crystal or the internal interface separating two crystals or phases, or the atomic mismatch between two grains such as grain boundaries. A tilt grain boundary consists of a mutual rotation of two adjacent grains abo ut an axis parallel to the grain boundary. A twist grain boundary, on the other hand, has a rotation axis perpendicular to the grain boundary. More detailed description of the grain boundary structures will be provided in Chapter 6. Void and interstitial c lusters, as well as a second phase can be regarded as bulk defects. These clusters can be considered as a collective of vacancies or interstitials, where their size can vary from sub micron to 100 microns. Sc ope of This Study The objective of this work is determining how dislocations and grain boundaries affect the thermal conductivity of UO 2 and the interfacial thermal resistance at the Si/SiO 2 interface . The simulation methods are molecular dynamics and phon on wave packet dynamics. A lattice dynamics model is also developed to help elucidate the effects of a grain boundary on thermal transport. In particular, t he effect of the decrease in the thermal conductivity due to the presence of dislocation is quantified, and its dependence on temperature and dislocation density is investigated. The interfacial thermal conductance of the Si/ SiO 2 interface and of UO 2 twist grain boundaries are also determined and the detailed mechanisms of phonon scattering are identified.


19 Table 1 1 . Di stinct sites in UO 2 crystal structure. Atoms Multiplicity Wyckoff Site symmetry Coordinates O 8 c 0.25,0.25,0.25 U 4 a 0, 0, 0


20 Figure 1 1 . A unit cell of fluorite structured UO 2 . As illustrated, red spheres refer to oxygen atoms and green spheres refer to uranium atoms. Figure 1 2 . A unit cell of diamond structur ed crystalline Si.


21 Figure 1 3 . Example crystalline SiO 2 structure s . A ) cristobalite; B ) quartz; C) stishovite; D ) tridymite. (All generated from CrystalMaker [ 13 ] ) Figure 1 4 . Illustration of edge and screw dislocations. (From Ref. [ 14 ] )


22 CHAPTER 2 THERMAL TRANSPORT IN DIELECTRIC MATERIALS Introduction to Thermal Conduction Thermal conduction is one of the three major mechanism of heat transfer (the other two being therm al convection and thermal radiation). Conduction is the transfer of heat, without the occurrence of macroscopic mass movements, through the exchange of kinetic energy with direct contact . F or insulating solid materials, heat cond uction is by lattice vibrat ions. To measure the ability to conduct heat, the concept of thermal conductivity is introduced . When the temperature distribution in a material has reached steady state, the thermal conductivity can be determined from [ 15 ] : ( 2 1 ) where i s the heat current, is the thermal conductivity tensor, and T is the temperature. The S I units of thermal conductivity are W atts per meter Kelvin, Wm 1 K 1 . Phonon mediated Thermal Transport In electronic insulators, where there are few if any free electrons, the main heat carriers are lattice vibrations, or phonons which are the quantum equ ivalents of classical lattice vibrations. Much of the following discussion is based on Ref [ 16 19 ] . Latti ce Vibrations and Phonons At 0 K temperature, the atoms of a solid crystalline material reside at the equilibrium position, where the solid has a minimum energy. However at a finite temperature atoms in solids generally vibrate around their equilibrium position s although the displacements are rather small.


23 In 1907, Einstein proposed a theory regarding the heat capacity of solids, where he suggested that each atom in a solid vibrates independently about its equilibrium position with the same frequency . As a result, a solid with N atoms consists 3N harmonic oscillators. And each atom vibrates as if all other atoms are keeping still. However, this is a tool simple treatment of a solid. Due to the interactions among atoms in the solid, the solid is no longer treated as atoms that vibrate individually, but, in a collective manner, as lattice waves. Each of these travelling waves corresponds to a normal mode of vibration, and the thermal energies of a material are distributed amongst all these modes . A wave like solution that describes the atomic displacements, to the dynamic equation of the lattice is describe as Eq. 2 2 [ 16 ] . ( 2 2 ) Here, is referred to as normal coordinates or normal modes. Therefore, each atomic displacement can be represented by a linear combination of these normal coordinates; a phonon is a quantum of energy of a normal mode. Each normal is, the lattice vibration is expressed as a collection of harmonic oscillators. [ 16 ] This is an analo gy to the fact any signal, e.g. a rectangular function, can be expressed as a summation of cosine functions, with the help of Fourier transformation [ 20 ] . And in fact, the normal coordinates transformation is carried out by discrete Fourier transformation [ 21 ] . One of the important features of lattice waves is the dispersion relation, which describes the relationship between the frequency and wavevector of a wave. The group


24 velocity , , which describes the energy transportation rate along the wave, is defined as, ( 2 3 ) where, is the wave frequency and k is the wavevector. Modes can be classified according to orientation relationship between the propagation directions and the direction of atomic motion. In longitudinal modes, the direction of the atom movement (polarization direction) is parallel to the wave propagation dire ction , while in transverse modes atomic motion is perpendicular to the wave propagation direction. For a three dimensional lattice, where there is N atoms in a primitive cell, there are N longitudinal modes and 2N transverse modes. Along some high symmetry directions, some transverse modes are degenerate. [ 17 ] For cases where there is more than one atom in a primitive cell, depending on the relative motion of the atoms, the vibrational modes can be further distinguished as acoustic mode, when atoms move in phase with each other; and optical m odes, when atoms move out of phase. Lattice Thermal Transport One approach to understand the thermal transport through phonons is the containing a certain number of pho nons, which are considered as particles similar to gas molecules. The description of the thermal conductivity is then mapped onto the kinetic theory of gases. [ 17 ] The collisions of these phonons with crystal impurities and surfaces, as well as the collisions among phonons create thermal resistance inside the material resulting a finite thermal conductivity. One difference between a physical as


25 and a phonon gas is th at the number of phonons is not conserved while the number of particles in a gas is . To understand the thermal transport in the phonon picture, the concept of phonon phonon scattering needs to be introduced. Phonons are quasi particles, and any physical mass or physical momentum, since at equilibrium there is k , known as the quasi momentum, phonon momentum, or crystal momentum is assigned to phonons. Under the harmonic approximation, p honons are non interacting. For a perfect, harmonic, infinite lattice the thermal conductivity would be infinite since there would be no scattering processes to impede the movement of the phonons. However, in an anharmonic perfect infinite crystal lattice, the thermal conductivity is not infinite. The anharmonic terms in the potential part of the Hamiltonian, that is those terms higher than second order in spatial derivatives, describe the multi phonon processes. Among these anharmonic terms, the cubic and quartic terms are mo st important. There are two distinct kinds of phonon phonon scatterings: normal processes and Umklapp processes. Figure 2 1 illustrates a three p honon process, in which two phonons annihilated during the scattering and a third phonon is created. For the normal process, where the newly created phonon wavevector is inside the first Brillouin zone, no crystal momentum is created; for the Umklapp proce ss, as the two incident phonon wavevectors are close to the first Brillouin zone boundary, the new phonon wavevector is outside the first Brillouin zone, which is not allowed due to the confinement of the zone. Consequently, the extra momentum is transferr ed to the crystal. These two processes can be summarized in the following equation for a three phonon scattering:


26 ( 2 4 ) Where there is G =0 in normal process and G process . Now consider the thermal transport in crystals. Since the transfer of heat is transport of phonon energy, under a constant temperature gradient there is a net momentum (nonzero) of phonons in the direction of the temperature gradient. The normal p rocesses do not create resistance to the phonon transport since they do not change the net momentum; the Umklapp process result in a loss to the heat current since the total momentum is not conserved. This loss of heat current results in thermal resistance . Normal modes contribute indirectly to the thermal resistance as they provide mechanisms for the redistribution of the phonon distribution back to the equilibrium distribution. To link the above microscopic concepts to the macroscopic quantity of thermal conductivity, a phenomenological model from elementary kinetic theory of transport can be constructed by introducing the relaxation time and mean free path . The relaxation time is the average time interval between two phonon scatterings, the mean free path is the average distance travelled by a phonon between two scatterings. These two variables, in a very simple way, represent the detailed phonon scattering details. In this simple dynamic model, phonons are considered as rigid particles like gas molecu les. Thus, suppose in a simple one dimensional scenario, the thermal conductivity is simply given by [ 17 ] , ( 2 5 )


27 where is the specific heat of the material and is the average velocity of the phonon. Although this formula is rather simple, it often fits experiment very well and can be applied to both insulators and metals. [ 15 ] Interfacial Thermal Resistance Inte rfacial thermal resistance, or Kapitza resistance is the thermal boundary re s istance presented at an interface between two blocks of materials, where a temperature drop is observed. It has the unit of W 1 m 2 K , and is the inverse of interfacial thermal conductance . Assuming a heat flow perpendicular to the interface, then the thermal conductance can be expressed as [ 22 ] , ( 2 6 ) where G is the interfacial thermal conductance, J T is the temperature drop at the interface. The resistance has two basic origins: the structural disorder associated with the interface and the different properties of the heat carriers o n the two sides of the interface. One or other of the mechanisms is operative at all interfaces; both are present at some interfaces. In the case of electrical insulators, including most ionically bonded systems and many semiconductors, the dominant heat carriers are phonons. As a result, the Kapitza resistance is governed by details of the interfacial structure and availability of vibrational states on the two sides of the interface. There are two general theory account for the phonon mediated interfacial thermal transport: acoustic mismatch model (AMM) and diffuse mismatch model (DMM) . The AMM model [ 22 24 ] suggests that the phonon scatterings at the interface is a result of the mismatch in the acoustic impedance (analogous to the refractive index in optics) . The DMM model, on the other hand,


28 assumes that all the incident phonons are scattered randomly at the interface. The amount of energy of the incident phonon emitted to the two sides of the interface is proportional to the available states in the r espective material. [ 24 ]


29 Figure 2 1 Illustr ation of Umklapp process and normal process . From Ref [ 12 ] .


30 CHAPTER 3 METHODOLOGY Molecular Dynamics Simulation of Thermal Transport A tomistic simulation, as the name suggests, is used to study materials behavior s with an atomic scale resolution. All atomistic simulation s require a description of the interactions among the atoms. D epending on the description, the se can be classified as quantum or classical. In quantum descriptions , the electronic structure of the interacting atoms is considered explicitly , and treated by quantum mechanics . In the classical descriptions , atoms are generally treated as structureless, with the description of the interactions taking into account the electronic structure implicitly. The atomic motion is then determine of motion . In molecular dynamics (MD) system in time, with each atom mov ing under the forces applied by surrounding atoms . The interatomic potential energy function , or just potenti al, is a n empirical mathematical description of the interactions between atoms in real materials based on their relative locations and other factors, such as local coordination . For the case of UO 2 , t here are twenty or more such interatomic potentials in the literature; however, there i s no single potential that gives a good representation of all physical properties. [ 25 , 26 ] For this work , the Busker potential [ 27 ] is chosen because it prov ides a reasonably good representation of many physical properties of UO 2 and it gives a stable dislocation structure at high temperatures. Moreover, it shows a thermal conductivity that is reasonably consistent with experimental data after an anharmonic co rrection [ 7 ] . The Busker potential is a Bucking ham type potential where the interaction between two atoms are depicted by Eq. 3 1 :


31 ( 3 1 ) where , and are parameters for diff erent pairs of atomic species, is the interatomic distance between atom i and j. The values of these parameters are listed in Table 3 1 . T he less computational expensive charge neutralized real spa ce direct summation method [ 28 ] is us ed in preference to the Ewald summation to calculate the Coulomb sums. For the case of Si/SiO 2 interface, the extended Stillinger Weber potential developed by Watanabe et al. [ 11 ] is used to describ e the interatomic interactions. In body interaction part V 2 and a three body interaction part V 3 : ( 3 2 ) where i, j, k are the label for each atom in the system and eV is the energy unit that is the same as the original Stillinger Weber (SW) potential. The two body part is described by: ( 3 3 ) where , , , and are parameterized for different interaction pairs: S i Si, Si O and O O, and r ij denotes the distance between atom i and j, assuming the length unit Ã… in SW potential. The two body term differs from the SW potential by introdu cing a bond softening function , described by :


32 ( 3 4 ) ( 3 5 ) ( 3 6 ) ( 3 7 ) where m 1 , m 2 , m 3 , m 4 , m 5 , R and D are parameterized to reproduce th e cohesive energy of Si O bond. Clearly, the bond softening function is effective only when the pair interaction is considered between Si and O. The three body interaction is in the form of Eq. 3 8 : ( 3 8 ) where the function h loops over each atom of a triplet as center (e.g. in the triplet jik, atom i is considered as the center) and has a similar form as in SW potential: ( 3 9 ) where , and are parameters for different triplets as denoted by the subscripts and is the bond angle. All the parameters for the SW potential are summarized in Table 3 2 .


33 This potential has been proven capable of describing the Si/SiO 2 interface as well as its formation. [ 29 , 30 ] In the extended SW potential, the interaction function between Si atoms is exactly the same as the SW potential for pure Si [ 10 ] . As shown by the silicon phonon dispersion curve in Figure 3 1 , determined by lattice dynamics calculations, the SW potential gives a reasonably good representation of the LA mode. However, it is not able to reproduce the flatness of TA mode near the band edge due to the lack of long range interatomic interactions in the SW description . However, as will be shown later, the energy transmission coefficient for higher frequencies (>6 THz) phonons is almost fre quency independent; thus this aspect of the dispersion curve will not significantly affect the results. Further, the poor descriptions of the two optical modes (LO and TO) is not crucial to the work because they are thought not to be significant carriers of heat in Si [ 31 ] , a conclusio n later confirm ed in our simulation s . Anharmonic Correction Due to the fact that the classical potential used does not reflect all of the physics of the real system, there is considerable difference between the calc ulated thermal conductivity and the experimental value . Here, follow ing an anharmonicity analysis [ 32 ] , the results are systematically corrected . By combining the Gr üneisen relation of thermal expansion coefficient conductivity can be expressed in terms of a few simple materials constants , including the heat capacity , bulk modulus, Debye temperature and thermal expansion coefficient. The Gr üneisen rel ation states that the thermal expansion coefficient is a result of the anharmonicity of interatomic interactions. Also, as discussed above, the anharmonicity gives rise to the finite thermal conductivity of material. [ 7 ] Thus the thermal conductivity


34 and thermal expansion are linked. calculated value is then given by, ( 3 10 ) where is the thermal expansion coefficient, is the bulk modulus is the Debye temperature and is the constant volume specific heat. All of these parameters are calculated using General Utility Lattice Program (GULP) [ 33 , 34 ] , with experimental values taken from the IAEA Database [ 35 ] . All of the data shown below include this correction. Non Equilibrium Molecular Dynamics The direct method , first introduced by Jund and Jullien [ 36 ] , o f the non equilibrium MD (NEMD) is used calculate the thermal conductivity of a material. This method introduce s a steady heat current in the system, and the temperature profile is direct ly determined inside the simulating system. A typical simulation box is a p arallelepiped with dimensions ( L x , L y , L z ) , with L z much longer than the other dimensions as shown in Figure 3 2 . For the purposes of analysis , the structure is divided into several slices of equal thickness along L z , with each slice contains a few hundreds atoms. After the system is equilibrated at the desired temperature, a small amount of energy is inserted into the first or first few slice(s) (red region in Figure 3 2 ) by changing the kinetic energy of the atoms; the same amount of energy is taking out of the system from the slices located in the center of L z (indicated by blue region in Figure 3 2 ). In such a way, a steady energy flow is generated. The average temperature for each slice is captured by averaging the kinetic


35 energy of t he atoms [ 2 ] . After a large number of simulation steps (depending on the structure size, interatomic potential) , the temperature profile inside the system can be conductivity can be determined directly . Phonon Wave Packet Dynamics The phonon wave packet dynamics (PWD) approa ch is used to study the Kapitza resistance of interfaces. To generate the wave packet, the wave vector in the first Brillouin zone of the primi tive cell of the material of interest, also known as the k vector, need to be properly defined. In reciprocal space, the wave vector k is defined as (k x , k y , k z ), where the direction along the wave packet propagation is defined as k z ; while k x and k y are p erpendicular to k z but orthogonal to each other. The core idea of PWD is to form a wave packet of phonons from a single branch with the frequency being a narrow Gaussian distribution. [ 37 ] The following describe s the basics of generation and analysis of the phonon wave packet. Most of the descriptions are based on Ref. [ 24 , 37 , 38 ] and the physics background of the equations can be found in many introductory solid state physics books on phonons, for example Ref s . [ 16 , 19 ] . The wave packet is generate d as described in Eq. 3 11 , so that it is localized in both real space and reciprocal space. The initial atomic displacement is generated according by, ( 3 11 ) A subsequent inverse discrete Fourier transformation is used to obtain the normal coordinates , ,


36 ( 3 12 ) where and denote the displacement vector of the th atom in th primitive cell and the coordinates of the th primitive cell. is the total number of primitive cell and m is the mass of atom . is the amplitude of the phonon with wave vector and branch , and and are the correspond ing eigenvector and its complex conjugate for atom i respectively. T he initial displacement ( Eq. 3 11 ) can then be rewritten as : ( 3 13 ) w here the time dependence term is added so that the initial velocity can be calculated by d ifferentiation . T hen t he phonon velocity, , is found by merely the taking time derivative of the previous equation, ( 3 14 ) where is the phonon frequency. By setting t=0, Eqs. 3 14 and 3 17 are used to calculate the initial displacement initial coordinates, and initial velocity for each atom in a system, which serve as initial conditions for an MD simulation. T he incident phonon wave packet is generated close t o the defect of interests (e.g. i nterface, grain boundary or point defects) and is then launched towards it as a summation of propagating phonons. When the wave packet reaches the defect region , sca ttering events take place ; as i llustrated in Figure 3 3 , part of the wave packet is trans mitted while some is reflected.


37 By determining the energy of each atom, the energy ca rried by atoms at either side of the interface is determined, from which the energy transmission and reflection coefficients through the interface are determined. T he transmission coefficient, T , is the fraction of phonon energy transmitted through the in terface, and similarly for the R . T he energy being transmitted or reflected is determined by adding up the total energy of each atom in the trans mission side or reflection side. As the example in Figure 3 3 shows , the transmission part corresponds to the part to the right of the dashed line, while the reflection part corresponds to the part to the left of the dashed line. However, a problem may arise for this method, when, there is more than one atomic species in a system, especially in some ionic systems where the cation and anion have rather different energy. For instance, in a perfect UO 2 single crystalline system, after relaxation, the potent ial energy of a uranium atom is about 70.4 4 eV and that of an oxygen atom is about 17.02 eV. Thus, when performing the summation, the total number of atom, as well as their ratio on each side needs to be appropriately determined in order to avoid any err ors. Another way to determine the energies is to perform a n inverse Fourier transformation similar to Eq. 3 12 , since MD is able to capture a snapshot of a system at a given time step , which includes the position and velocity of each atom. This is done as the following equation: ( 3 15 ) where the t in the bracket denotes the variable at a given time t. Then the energy belongs to a certain mode at time t is determined by: ( 3 16 )


38 In doing so, the amount of energy being transmitted, as well as the vibrational mode to which it belongs to can both be determined. Calculating boundary conductance In order to calculate the Kapitza conductance of an interface from the results of phonon wave packet dynamics , a sampling scheme of the first Brillouin Zone (BZ) is required . In a supercell, where its dimensions are integer multiple (N x , N y , N z ) of a cubic unit cell with a lattice constant of L, the meaningful k points in the reciprocal space are , where i, j and k are integers smaller than N x , N y and N z respectively. In particular, the k x k y planes are sampled with an equal spaced mesh scheme with the . The k z direction is sampled with a finer mesh since the structure is much longer along the z direction . Moreover, the underlying symmetry of the structure is recognized here in order to reduce the computational time , as shown by Figure 3 4 , where only 1/8 of the first Brillouin zone is sampled. The total conductance , which is the inverse of resistance, at an interface can be calculated using the equation for an interface conduc tance by integration over the f irst BZ, [ 38 ] ( 3 17 ) where is the phonon transmission coefficient determined from the phonon wave packet dynamics simulations for phonon branch is th e reduced Planck constant and is the Bose Einstein distribution function at temperature T. The integration over the first BZ is converted to a sum over k points, with


39 the weighting factor, W( k ) assigned to each point determined by the associated k space volume, determined by a Voronoi construction over all the k points: ( 3 18 ) Computational Tool Several computational software tools are used for the study. For the part of calculation the effects of dislocations on UO 2 thermal transport and phonon wave packet dynamics on Si/SiO 2 interface , our in house molecular dynamics package i s used. Modifications were made in order to include the extended S tillinger W eber potential. For the phonon wave packet dynamics on UO 2 grain boundaries, the GULP [ 33 , 34 ] program is to obtain the phonon data , while the MD simulations are performed with LAMMPS [ 39 ] . A series of program s ha ve been written as interface between GULP and LAMMPS to automatically handle the data flow, as well as for generating initial displacement and normal mode analysis. The MATLAB ® [ 40 ] software is used to solve the lattice dynamics model. In addition , AtomEye [ 41 ] and CrystalMaker [ 13 ] are used for visualization purpose for the atomic simulations.


40 Table 3 1 Parameters for Busker potential. [eV] [ Ã… ] [eV Ã… 6 ] 0 0 0 1761.775 0.35643 0 9547.96 0.2192 32 Table 3 2 Parameters of extended SW potential from Ref. [ 11 ] 7.049556277 0.6022245584 4 0 115.364065913 0.9094442793 2.58759 2.39370 12.292427744 0 0 2.24432 1.8 1.4 1.25 m 1 0.097 m 2 1.6 m 3 0.3654 m 4 0.1344 m 5 6.4176 R 1.3 D 0.1 16.404 1.0473 1.8 1/3 2.9572 0.71773 1.4 0.6155 3.1892 0.3220 1.65 1/3 10.667 1.93973 1.9 1/3 0.25 1.4


41 Figure 3 1 Phonon dispersion relation of silicon along [001]. Solid lines are calculated from Stillinger Weber potential from lattice dynamics using GULP [ 33 ] ; dashed lines are reproduced from experiment [ 42 ] . 0 2 4 6 8 10 12 14 16 18 20 0 0.2 0.4 0.6 0.8 1 Phonon frequency [THz] kz [2 /a 0 ] TA LA LO TO


42 Figure 3 2 Schematic of simulation box for calculating thermal conductivity using non equilibrium molecular dynamics . Figure 3 3 Illustration of a scattering event. The vertical axis is the atomic displacement along z axis; the horizontal axis indicates the z coordination of the structure.


43 Figure 3 4 First Brillouin Zone sampling points of a diamond lattice in k z =0 plane. The area delineated by the blue line s is the irreducible region.


44 CHAPTER 4 EFFECTS OF EDGE DISLOCATIONS ON THERMAL TRANSPORT IN UO 2 1 Introduction Dislocation Structure in UO 2 Both edge and screw type dislocations can be present in UO 2 . As noted in Chapter 1, edge dislocation can typically be described as the insertion of an extra half plane of atoms. However, in the some ionic crystals, inlcuding UO 2 , two extra half planes of atoms are required in order to maintain charge neutrality. In the fluorite structure, the primary slip systems are {100}<100> and {110}<110> [ 44 ] . Unlike the case in metals, where the Burgers vector is determined by minimizing the elastic energy, the criterion of choosing the <110> Burgers vector in UO 2 is essentially minimizing the electrostatic energy, since it is much larger than the elastic energy. [ 45 ] Another feature of dislocations in UO 2 , as well as some ot her ionic crystals, is the effective charge around the dislocation cores. By counting the bonding around the core ions [ 44 ] , there are excess positive or negative charge on the {100}<110> configuration, whil e the {110}<110> configuration is un charged. These three cases are illustrated in Figure 4 1 . The screw dislocation s on all slip planes are uncharged, as they only involve displacement s along the dislocation line. T he most easily activated slip system depends on both temperature and stoi chiometry. At low temperatures (less than 1000 K), it is rep orted that {100}<110> is the dominant slip system , while the {110}<110> dislocation becomes a ctive as the 1 This chapter is based on 43. Deng, B., A. Chernatynskiy, P. Shukla, S.B. Sinnott, et al., Effects of edge dislocations on thermal transport in UO2. Journal of Nuclear Materials, 2013. 434 (1 3): p. 203 209.


45 temperature increases into the regime present in fuel during service . [ 46 ] Thus, the { 110 }< 110 > edge dislocation system is the main focus of this work . Dislocation Evolution The evolution of dislocation structures during burn up is well documented . In particular, Nogita and Une [ 47 ] examined dislocations in a UO 2 specimen in the peripheral region of fuel pellets, which had been irradiated for a wi de range of burn ups, from 6 to 83 GWd/t in an light water reactor ( LWR ) . The TEM image [ 48 ] showed tangled dislocations and interstitial type dislocation loop s at a burn up of 44 GWd/t (gigawatt days/ton of heavy metal) . Moreover, the density of dislocation increases with burn up, as indicated in Figure 4 2 . A t 83 GWd/t . The dislocation density , N ( in units of m/m 3 ) , was characterized as having an exponential increase over the range of 6 44 GWd/t burn up (Bu), characterized by the empirical relationship: [ 47 ] ( 4 1 ) Although the measured dislocation density at the burn up of 83 GWd/t was almost the same as that at 44 GWd/t, this value was thought to be an under estimate , since the method [ 49 ] used for measur ing dislocation density could not separate highly tangled dislocations. The extrapolated dislocation density at the 83 GWd/t burn up is 4.2X10 15 m 2 according to Eq. 4 1 . As also shown in Figure 4 2 , the dislocation density / burn up relationship has been successively modeled by the mechanistic MFPR (Module for Fission Product Release) code [ 50 ] . The FRAPCON fuel performance code [ 51 ] includes models for the thermal conductivity of the fuel rod. The functio nal form that considers the temperature and burnup for the thermal conductivity of UO 2 fuel pellets is given by: [ 52 ]


46 ( 4 2 ) The functional forms of , and are prescribed. The values of the parameters A , B , E , F and the parameters within , and are determined from experimental data. Considering Eq. 4 2 as containing two terms, the first term accounts for the phonon mediated heat transport. The A+BT contributio n accounts for the conductivity at the beginning of life (BOL), including the phonon phonon scattering (The BT term) and the scattering of phonons from microstructure presen t at BOL (the A term). The term accounts for the fission products in the crystal matri x; the term accounts for the irradia tion induced defects, while accounts for the temperature dependence of annealing of irradiation defects. At the BOL the only contribu tion that is present is the term. The second term describes the phonon electron contribution , which is only important at high temperatures (above Debye temperature) ; [ 6 ] this contribution cannot be characterized in MD simulations, and will not be considered further here. Theoretical Explan a tion Effects on Thermal Transport Klemens [ 53 ] provided a theore tical analysis of the effects of various defects on thermal conductivity arising from phonon defect scattering . In this approach , at high temperature (above the Debye temperature) it is appropriate to add the relaxation rate s of individual processes in an inverse manner. In particular, for phonon phonon and phonon dislocation scattering, the total relaxation rate can be written as


47 ( 4 3 ) where the subscript U and D indicates the phonon phonon scattering and phonon disloc ation scattering, respective [ 53 ] The contribution of the dislocation s can be further separated into scattering from the dislocation core and scattering from the dislocation strain field. In Klemens formalism, t he relaxation rates for the scattering due to Umklapp process, dislocation core, the strain field associated with screw dislocation s, and the strain field associated with edge dislocation s are given by [ 53 , 54 ] ( 4 4 ) ( 4 5 ) ( 4 6 ) ( 4 7 ) where is a factor indicating the orientation of temperature gradient with respect to the dislocation line (1 for perpendicular; 0 for parallel; 0.55 for random orientation ), is the Gr ü neisen anharmonicity parameter, is the angular fr e quency, is the Debye frequency, is the shear modulus, is the volume per atom, is the dislocation density, , , are the average, transverse and longitudinal velocities respectively, b is the magnitude of is the Poisson ratio . T he thermal


48 conductivity can then be dete rmined by incorporating Eqs. 4 4 ~ 4 7 formulation for the lattice conduct ivity [ 55 , 56 ] ( 4 8 ) Here is the Planck constant and is the D ebye te mperature. The resulting contributions to the thermal conductivity have distinctive temperature dependences. In particular, the Umklapp process es yield T 3 dependence at low temperatures, reflecting the well known temperature dependence of the spec ific heat in the quantum regime [ 53 ] . At high temperatures, this yields T 1 dependence , reflecting the dominant effe ct of phonon phonon scattering; this is captured in the BT term in Eq. 4 2 . All in reactor conditions correspond to this high temperature regime. The analys e s of both the core and strain fields of dislocations lead to contribution s to the phonon scattering rate that are independent of temperature. [ 53 ] These contributions can be related back to the A term in Eq. 4 2 . Simulation Setup The setup for the thermal transport simulations is a conventional 3 D periodic rectangular cylinder with dimensions L x L y L z and heat inflow and outflow similar in Figure 3 2 . T he dislocations are the only defects in the simulation cell and their cores are located half way between heat source and sink. The use of a periodic system imposes constraints on the struct ure that need to be considered. In particular, the dislocations are created in dipole pairs such that the total Burgers vector around the e ntire simulation cell is zero. One edge of the simulation cell is parall el to the dislocation line and the size of th e cell in this direction (L y ) does not affect parameters of the


49 simulations, such as dislocation density. Consequently, L y is set to be about 2.2 nm, which is the minimum required by the potential cutoff. This relatively small size essentially forces the d islocations to remain straight. The size of the simulation box in the two other directions (L x an L z ) can be used to control the dislocation density. However, there is a significant size effect on thermal conductivity associated with the length of the c ell L z along the temperature gradient direction. [ 57 ] Consequently, L z is set to about 155 nm and is kept the same for all structures. There is no such effect associated with the cross sectional area of the simulation box [ 57 ] ; This has also been verified in the system presented . For example, at 1600 K, by doubling L x and considering a dislocation free crystal , the thermal conductivity changes from 5.11 to 5.15 Wm 1 K 1 (before the anharmonic correction, see below). Therefore, the size of the system along the x axis is used to set the dislocation density. Ideally, this density would be set to correspond to the sp ecific burn up. However, this would make the simulation cell larger than is computationally practical. As a result the lowest dislocation densities considered are twice as large as the largest densities observed experimentally at high burnup. Dependence o f Thermal Conductivity on Temperature and Dislocation Density To determine the temperature dependence of the contribution of the dislocation to the thermal conductivity, the thermal conductivities of both a model containing dislocations and a defect free s ingle crystal are calculated for comparison . In Figure 4 3 , the temperature dependence of the thermal conductivities of the two systems as shown from 800 to 1800 K. It can be clearly seen that the presence of the dislocations results in a decrease in the thermal conductivity. To see this decrease more clearly, the thermal conductivity differences are shown in Figure 4 4 (a) , from which the magnitude


50 of the reduction decreases with increasing temperature can be identified . To quantitatively extract the impact of dislocations on thermal conductivity, the elementary kinetic theory for are applied , Eq. 4 3 . Since the thermal conductivity is proportional to the relaxation time , and the relaxation times add inverse ly as shown in Eq. 4 3 , the temperature dependence can be determined from the difference in the inverse thermal conductivities (the thermal resistivities) of the system s with and without dislocations. This decrease in inverse conductivity is given in Figure 4 4 (b), from which it can clearly be seen that the effect of the presence of the dislocations is independent of temperature. The error bars are larger at high temperature, due to the larger temperature fluctuations at high temperatures. This temperature independence is consi predictions [ 53 ] . To quantitatively analyze the effects of dislocation density, t hree structures with dislocation densities of 0.83×10 16 , 1.10×10 16 and 1.68×10 16 m 2 were built; their thermal conductivities at different temperatures are shown in Figure 4 5 . T he lowest of these densities corresponds to a dislocation density about twice that of a pellet at high burn up. As shown in Figure 4 6 , the resistivity shows a good lin ear relationship with the density of dislocation , with similar trends at all temperatures. Recalling the previous discussion on the reduction of resistivity ( Figure 4 4 resistivity consists of a temperature dependent intrinsic contribution and a temperature independent dislocation contri bution (Eq. 4 3 ) . The fact that all of the line fits are parallel to each other ( again show s that the effect of dislocations on thermal conductivity is temperature ind ependent , even at these high dislocation densities. B y av eraging the fits to the slopes of the fitting, it is


51 possible to obtain the thermal resistivity (R D ) per unit density (d) of dislocation to be R D =3.9 X10 18 W 1 m 3 K. This term in formed in analogous to the Kapitza resistance , where the Kapitza resistance [ 58 ] at flow per unit area . Here, to formally define this contribution, the term effective thermal conductivity is employed to descr ibe the effect of dislocations, ( 4 9 ) where is the thermal conductivity of the same structu re without dislocation Comparison of Various Dislocation S ystem s In addition to the {110}<110> dislocation, the UO 2 thermal conductivity with a set of {100}<110> edge dislocation s is determined . Figure 4 7 compares of the thermal conductivity of this structure with the equivalent structure of the {110}<110> dislocation. With similar dislocatio n densities (1.69X10 16 and 1.68X10 16 m 2 ), they show a comparable reduction in thermal conductivity. This is consistent with Klemens theory in that the major contribution of a dislocation to decreasing the thermal conductivity comes from the phonon scatter ing on the strain field rather than the scattering on dislocation core, which is similar to the scattering on point impurities due to the mass difference . [ 54 ] Since both dislocation structures have two half lattice plane s missing (or inserting) perpendicular to their slip direction (110) and both structures have similar system size, their strai n fields should also be similar. Theoretical Model and Fuel Performance Code The functio n of fuel performance codes is t o predict the evolution of a fuel pellet during burn up. Therefore, in this section, the results from the MD simulations of dislocation effects are used to parameterize the thermal conductivity temperature


52 formula used in FRAPCON code. In so doing, the formula descr ibing the temperature dependence of the thermal conductivity is parameterized . Of course, burn up has multiple effects on the microstructure, each of which in turn has an effect on the thermal conductivity. The objective here is merely to demonstrate that data obtained from simulation can be passed up to a fuel performance code, not to provide realistic parameters to the code itself. Further, this fitted model is compared to the numerical solution of Klemens 4 8 ) . The incorporation of MD results for single crystal thermal conductivity and thermal expansion into the FRAPCON code has previously been examined by Vega et al. [ 59 ] The MD simulation results of temperature dependen t thermal conductivity of various sets of dislocation density models are fit to ( 4 10 , which is the simplified formula of Eq. 4 2 . Specifically, the terms whose effects are not taken into consideration of the dislocation model are omitted , including the fission products effects (f(Bu)), defects annealing (h(T)) and electronic contrib ution. By relating the dislocation with burnup using the empirical relation in Eq. 4 1 and in the spirit of the thermal conductivit y functional in FRAPCON code, there i s , ( 4 10 ) As discussed above, the (A+BT) is the contribution of lattice thermal conductivity of pristine fuel [ 59 ] , whose structure contains many microstructural elements whose effects are included in A. Consequently, parameters A and B are fitted from the thermal conductivity data of a perfect structure, and with the f itted A and B, the function (Eq. 4 9 ) can be fur ther fitted with the thermal conductivity of all dislocation structures . The


53 term takes into account the effects of additional dislocations on thermal conductivity, which are produced during burn up. Thus, the thermal conductivity is obtained , which is linked with burnup. The fitting parameters are listed in Table 4 1 and the fitted results are shown in Figure 4 8 and F igure 4 9 . By doi ng this, an overall trend of temperature and dislocation effects on thermal conductivity is determined . Moreover, by inserting parameters into Klemen s Callaway model, there is the numerical solution of this model. The parameters us ed ( Table 4 2 ) are again calculated using GULP; the comparison of the two models is shown in Figure 4 8 and F igure 4 9 . As can be see n , both models have shown the same trend, however, the effects go t from M D fitting model is roughly 10% smaller. Conclusion s The general tr end of the dislocation density dependence, as well as the temperature independence of the effects is in agreement with Klemen s prediction ; but with an overall smaller effect. Overall, there is a 10 20 % decrease in thermal conductivity due to the presence of dislocations for the high densities u sed in the study. However, it is note d that the smallest dislocation density in the simulation is roughly twice larger than the saturate d dislocation concentration in nuclear reactors. Although simulations at lower densities would be desirable, the linearity of the resistivity over the taking this into a ccount, it seems that dislocations alone do not have a significant impact on thermal conductivity; in fact, according to the fitting model, there should be less than a 5% reduction in thermal conductivity at the saturated dislocation density. However,


54 poss ible effects of point defect segregation to dislocations on thermal transport have not been considered here. This work also shows that by incorporating the MD results into the thermal conductivity functional used in FRAPCON fuel performance code, a general formula for estimating the effects of dislocations on thermal conductivity in a larger scale can be developed .


55 Table 4 1 Fitting results for different dislocation density models compared to FRAPCON Parameters Values A (m K W 1 ) 0. 08076 B (m W 1 ) 1.839X10 4 C (m K W 1 ) 3.659X10 4 D 5.671 Goodness R 2 0. 9767 Table 4 2 Parameters used for Klemens model C 11 (GPa) 526 a (nm) 5.481 C 12 (GPa) 118 D (K) 500 C 44 (GPa) 118 1.55 Figure 4 1 Three cases for dislocation core structures. (The spacing is exaggerated to illustrate the atomic stacking; translucent atoms indicates lattice plane below.)


56 Figure 4 2 Increase in d islocation density with burnup (squares) and fits from the MFPR code. Reproduced from [ 50 ] . Figure 4 3 Anharmonic correction of MD results. Filled markers take into account the anharmonic correction.


57 Figure 4 4 Change of thermal conductivity between d=0 and d=1.68X10 16 m 2 .


58 Figure 4 5 Temperature dependence of thermal conductivity of different dislocation density structures. Lines indicate fits of the MD data to Eq. 4 9 .


59 Figure 4 6 Dislocation dependence of resistivity at various temperatures. Solid lines are the linear fits.


60 Figure 4 7 Comparison of the thermal conductivity of {110}<110> and {100}<110> dislocation structure (d=1.68X10 16 m 2 ).


61 Figure 4 8 Comparison of dependence of thermal conductivity on dislocation density on with the Klemens Callaway model (denoted Klemens) and by fitting the MD data to the FRPACON model (denoted Fitting).


62 F igure 4 9 Comparison of t emperature dependence o f thermal conductivity of Klemens Callaway model and fitting to the FRAPCON model for a dislocation densi ty of d= 1.68x10 16 m 2 .


63 CHAPTER 5 KAPITZA RESISTANCE OF Si/SiO 2 INTERFACE 2 Silicon is a widely used semi conducting material, and silicon oxide is a very good dielectric material which is usually used as substrates or protective layer in combination with silicon. For example, in the silicon on insulator (SOI) technology in microelectronics, the SOI wafer structure is can be made up of a thin layer of SiO 2 , used as dielectrics, sep arating single crystalline Si from current leakage as well as improving performance. [ 8 ] These Si/SiO 2 interfaces c an act as obstacles for heat flow in microelectronics, where a huge amount of heat is generated in very limited spaces. very important to understand the thermal resistance at the interface. Introduction There are s everal experimental studies o f the Kapitza resistance of Si/SiO 2 interface. Lee and Cahill [ 61 ] method. They measured the thermal conductivity at various SiO 2 layer thickness and determined the Kapitza resistance to be 2 X 10 8 m 2 K W 1 by attributing the drop of SiO 2 thermal conductivity entirely to the interface. Kato and Hatta [ 62 ] also measured the Kapitza resistance betwee n a silicon substrate and a thermally oxidized SiO 2 film to be 4.5 X 10 9 m 2 K W 1 but with a large scattering of data (9 X 10 9 m 2 K W 1 ). Hurley [ 63 ] have measured the Kapitza resistance at a Si/SiO 2 interface experimentally. Their defect region (4.5 nm). The defect region was characterized as a silicon oxide, which they believe to 1 This chapter is based on Ref. 60. Deng, B., M. Khafizov, D.H. Hurley, and S.R. Phillpot, Kapitza resistance of Si/SiO2 interface. Journal of Applied Physics, 2014. 115 (8): p. 084910.


64 be SiO 2 since the specimen was annealed at 1000 C . T hey determined the Kapitza resistance to be 2.3 X 10 9 m 2 K W 1 . There have been previous molecular dynamics (MD) simulations of the Kapitza resistance of Si/SiO 2 interfaces. Mahajan et al. [ 64 ] estimated the Kapitza resistance to be G K ~0.5X10 9 Km 2 W 1 using an extended Stillinger Weber potential. Lampin et al. [ 65 ] calculated the Si/SiO 2 boundary resistance to be 0.4X10 9 Km 2 W 1 at 500 K with the Tersoff potential [ 66 ] . The y also found that this boundary resistance is large enough to change the thermal properties in the case of ultra thin buried oxide layers. Chen et al. [ 67 ] focused on how the strength of the coupling across the interface affects the Si/SiO 2 interface resistance using non equilibrium molecular dynamics (NEMD). They found that in the weak interfacial coupling limit, the boundary resistance is sensitive to the details of the interfacial structure; in the strong coupling limit the boundary resis tance is not sensitive to the details of the interface structure. In this strong coupling limit, the Si/SiO 2 boundary resistance was found to be 0.9X10 9 Km 2 W 1 . Thus while these th ree th eoretical analysis give consistent values of G K ~0.4 0.9 X10 9 Km 2 W 1 , they do not provide any insights as to which phonon branches and wavelengths are involved. To provide such insights, i n this study the scattering s of phonons at the Si/SiO 2 interface are characterized using the phonon wave packet dynamics (PWD) techniqu e. This approach has been extensively used in the study of phonon scatterings in various silicon microstructures . [ 24 , 37 , 68 ] While most of the previous studies focused on the role of intrinsic and extrinsic defects in silicon, [ 24 , 37 , 68 ] here the concentration is on a detailed description of the phonon scattering at the Si/SiO 2 interface.


65 Simulation Setup To investigate the thermal transport at the Si/SiO 2 interface, the simulation structure is set up in a manner analogous to the exper iment structure of Hurley [ 63 ] , where a thin layer of SiO 2 film is sandwiched between two blocks of Si crystal. Here, the two Si crystals have the same crystallographic orientation. This differs from the exper imental situation of Hurley et al. in which the two Si crystals form a twist grain boundary. However, due to the presence of the SiO 2 layer, the presence or absence of a misorientation between the two Si crystals does not significantly affect phonon dynami cs ; moreover, the structure allows smaller cross sections to be used, which has a significantly lower computational load. The simulation setups for the wave packet dynamics and NEMD simulations are sketched in Figure 5 1 . Note that the silicon crystals are extremely long in the structure used for wave pac ket dynamics so that the transmitted and reflected wave packets can be fully analyzed . To build the sandwich structure, the blocks of crystalline Si are fully quenched to very low forces and stresses; the lattice constant in the cross sectional directions during annealing is fixed to its bulk value. Because periodic boundary conditions are a pplied, t he SiO 2 is prepared initially cristobalite and strained (about 1%) in order to fit the Si lattices in the (001) plane. The entire system is then annealed at 2000K to enable bonding at the interfaces. The entire structure is further quenched s o all of the forces are very small, less than 10 7 eV/Ã… per atom. This eliminates any excess structural energy that could affect the phonon scattering simulations. The structur e used for PWD simulation s has 4X4 silicon unit cells in the (001) cross section al plane, and is about 1600 nm long. Various thicknesses of the SiO 2 layer, from 0.8 to 24 nm are used. The same approach is


66 applied in preparing the structure used for the NEMD simulations, only with much shorter Si blocks (about 50 nm). Energy Transmissi on The study start s with the sandwich structure which has very thin layer of SiO 2 (0.8nm in thickness). To develop an understanding of phonon mediated interfacial resistance u sing the PWD method, simulations with various incident phonon wave packets for al l six branches , one longitudinal acoustic (LA) mode, two transverse acoustic (TA) modes, one longitudinal optical (LO) mode and two transverse optical (LO) modes, are performed. Figure 5 2 shows how the energy transmission coefficients of LA and LO phonons for all available k xy components change as the phonon frequency increases. It is note d from Figure 5 2 that there is no qualitative change in the transmission coefficient of the LO modes compared to the LA modes; that is, the transmission coefficient is largely just a function of the frequency regardless of the phonon symmetry . Moreover, it can also be seen that the transmission coefficient is not strongly sensitive to the phonon branch. This conclusion is reinforced in Figure 5 3 , which shows the energy transmission coefficients of only the k 00 phonons for all four non degenerate phonon branches. Klemens [ 53 ] has determined the frequency dependence of phonon scattering probability f or various defects structures using perturbation theory. Within this approach the probability is a power function of the phonon frequency, with the values of the exponent determined by the defect. Thus the reflection coefficient to frequency is fit to a po wer equation, as following,


67 ( 5 1 ) coefficients to fit. As shown in Figure 5 3 , there is a reasonable good fit with a=0.225 [ 53 ] predicts b=0 for the ideal grain boundary and b=4 for point defects. T hus b= 0.515 indicates a stronger scattering than an ideal grain boundary, which is attributed to the disorder in the SiO 2 layer. The SiO 2 layer thickness is also varied to determine the energy transmission coefficients. The energy transmission is further decreas ed with increased thickness of the SiO 2 layer : as shown in Figure 5 4 , there is a uniform decrease in energy transmission as the SiO 2 layer becomes thi cker (l SiO2 =0.8 nm to l SiO2 =10.9 nm). The effect of thickness in the next section. Phonon Scatterings at Interface The phonon packet launched at the SiO 2 region consists of vibrational modes of the Si lattice. Thus scattering can only take place at the Si/SiO 2 interface and in the SiO 2 block which has di fferent vibrational properties. Consequently, when the wave packet frequency is low , i.e., its wavelengt h is larger than the space dimensions of most defects, the phonon can easily travel through the SiO 2 region with very little reflection. A simple calculation shows that a wavelength equal to the SiO 2 layer thickness of 0.8nm corresponds to a frequency of ~ 10 THz. (the first point of k 00 line in Figure 5 2 ), the transmission coefficient approaches unity; for thes e long wavelengths, the relatively narrow SiO 2 regions do not present a significant obstacle to the phonon waves. However, as the incident phonon frequency increases, its wavelength becomes more and more comparable with the dimensions of the


68 defected regio n, yielding more intense phonon scatterings. This can be seen in Figure 5 5 , which shows how the transmitted phonon energy of initially pure k 00 phonon s scatters into various k xy points . As can ssee , as the incident frequency increases, the transmitted energy distribution in k xy is more widely spread. For instance, in the f=8 THz case, the dominant energy share in k 00 is no longer to be seen. To better quantify this dispersion statistically, the term weighted wavevector deviation (d k ) and is defined in Eq. 5 2 , ( 5 2 ) where k 0 is the peak of the incident phonon wavevector, k n and E n are the transmitted/reflected phonon wavevector and associate d energy respectively, and the summation is over all the available phonon wavevector for the transmission or reflection phonons. d k essentially measures how much the average phonon momentum deviate s from its initial value . The larger the deviation, the less character of the initial phonon features is preserved. d k is plot as a function of incident phonon frequency in Figure 5 6 . The figure shows the LA tr ansmission and reflection data for l SiO2 =0.8 nm structure . Clearly , scattering increases with increasing frequency causing that scattered phonons to lose forward momentum . The rest of the momentum is scattered to the non normal directions as shown in Figure 5 5 . There is also less energy transmitted as the phonons lose forward momentum, and more energy is reflected back. As a consequence there is a decreas e in energy transmission with increasing frequency.


69 In addition to scattering into different directions, the phonon mode is also able to k conversion after the scattering at the interface. Figure 5 7 shows the phonon wave packet energy partition into various phonon modes as a function of incident LA phonon energies. The pink dotted lines indicate the initially pure LA phonon wave packet energy di stribution for various incident phonon frequencies. The blue lines and red lines are the energy distributions of LA and TA modes, with dashed lines indicating reflection and solid line s indicating transmission. It can be directly see n that the energy shift (peak shift between blue and red) and mode conversion (phonon modes other than blue) as the incident phonon frequency increases; the LA TA conversion cannot take place at high frequency as t here are no TA modes available as can be seen in the dispersion c urve in Figure 3 1 . It is interesting that the LA modes still maintain their original frequency, but the frequency of the converted TA modes gradually shifts away fro m the incident phonon frequency ind icating more intense scattering events . Such phonon mode conversion is known to take place when phonons interact with interfaces, [ 69 , 70 ] a result of anharmonic phonon scattering at the interface. [ 71 ] This indicates that as phonon frequency increases (in the region that TA modes are still available), phonon mode conversion also increases, t hereby increasing the energy transmission coefficient. Thus, this provides evidence that anharmonic scattering s open up additional channels for phonon transmission thus reducing interfacial thermal resistance at high frequencies [ 72 , 73 ] .


70 Calculating Kapitza Resistance The energy transmission coefficients for individual modes alone are not sufficient to enable the K apitza resistance to be determined. The Kapitza resistance is therefore calculated using Eq. 5 3 , which has been introduced in Chapter 3. ( 5 3 ) T he factor W(k) in Eq. 5 3 is the volume weighting factor of the First Brillouin Zone . Since the mesh grid used in the k xy plane is equally spaced, the factor can be rewritten as . By first summing up the W z only, the relative contribution to the conductance at each (k x , k y ) is obtained , as shown in Figure 5 8 . Clearly, as k x or k y increas es, under the same frequency, their relative conductance drops rapidly because the k z component of the wave vector decreases. Although Figure 5 2 shows that all non perpendicular incident phonons have similar energy transmission coefficients, their contributions to the conductance can be different due to their different phonon group velocities and available k z dimension in the First Brillouin Zone. As d iscussed above, the system presented consists of two Si/SiO 2 interfaces and cristobalite SiO 2 T ) of the system can be written as [ 65 ] ( 5 4 ) where 1/ T is the total thermal resistance from simulation, 1/ K is the Kapitza resistance of Si/SiO 2 interface, the factor of 2 indicates the two interfaces, 1/ k SiO2 is the thermal resistivity of SiO 2 and l SiO2 is the length of the silica layer. T he pure Kapitza resistance of Si /SiO 2 is then determined in two steps. First, by varying the silica layer thickness l SiO2


71 and the total resistance is calculated at each l SiO2 value . Subsequently, a linear e xtrapolation is performed according to Eq. 5 4 , at l Sio2 =0 the Kapitza resistance of the individual interface is obtained. T he calculation of the first step is carried out using both NEMD and PWD methods. For t he PWD method, the full First Brillouin Zone sampling is only performed for the l SiO2 =0.8 nm and l SiO2 =2.2 nm structure. However, t he relative contribution to the total resistance from k 00 of both structure s are calculated and confirmed to be almost the same . Therefore, for the other thicknesses of SiO 2 layer, only the k 00 direction is sampled and the results are scaled with the same contribution obtained from l sio2 =0.8 nm k 00 data. This is justified by the dominant contribution a k 00 to the total conductance as shown in Figure 5 8 . The NEMD simulation is performed at 300 K, and by utilizing the Bose Einstein phonon distribution function ( Eq. 3 17 ), the Kapitza resistance at 300 K can be calculated using the PWD method. Figure 5 9 shows the thermal conductivity as a function of thickness of the SiO 2 layer as determined from the NEMD and PWD methods. The Si/SiO 2 Kapitza resistance, determined from the intercepts, is 1.4 8(±0.46) X10 9 m 2 K/W from NEMD and is 1.37(±0.42) X10 9 m 2 K/W from PWD. The good agreement of the results from PWD and NEMD methods suggests that they have captured the same physics of the interfacial thermal phonon scattering, discussed in the previous s ection. The relatively large error bars here are due to the linear regression since only a small number of data points are taken into account . The power law fit in the previous section is used to calculate th e Kapitza resistance. The transmission coefficie nt at any wavevector in the first Brillouin zone is determined by its corresponding frequency, following Eq. 5 1 . Therefore, by integrating all modes contributions, it gives the Kapitza resistance as 1.83(±0.09) X10 9


72 m 2 KW 1 . Although the value is not exactly the same from the atomistic method, it does provide a more convenient way to estimate the Kapitza resistance. Comparing with the multi magnitude measurements from experiments [ 65 ] , the results here are of the same magnitude. In particular, the atomistic model used in the work has a one to one correspondence with the experiment settings, thus the comparison is more reliable. In addition to the overall value, the PWD also provides a mode by mode contribution to the thermal conductivity. Using Fourier analysis, the mode wise contribution to Kapitza conductance is shown in Figure 5 10 . The acoustic modes contr ibute 88% to the total Kapitza conductance , while the optical modes contribute only 12%. This agrees with the conventional understanding that the acoustic phonons are the main heat carriers in silicon. [ 74 ] Conclusion In this work, th e Kapitza resistance at a Si/SiO 2 /Si sandwich structure is investigated by PWD and NEMD approaches. These methods provide similar values for the Si/SiO 2 Kapitza resistance; this value is also consistent with the experimental results and previous simulation s. It has demonstrated that the PWD approach provides far more detailed information about thermal transport through the interface. In particular, for the system considered here, the acoustic phonons are found to be the main contributor to the conductance t hrough the interface. The phonon anharmonic scattering provides additional channels for phonon transport. The phonon energy transmission at the interface is also analyzed using PWD. Frequency dependence of the reflection coefficient was found to be stronge 0.5 ideal planar defect (~ 0 ). This is associated with the disordered nature of the studied


73 interface. The fact that transmission coefficient is fitted very well by a power law can be used in the mesoscale m odels of the phonon transport [ 75 ] . However, it is noted that, as a classical simulation, the potential used will have an effect on the final values.


74 A B Figure 5 1 Simulation setup for A. Wave packet dynamics simulation; B. Non equilibr ium molecular dynamics simulation. The numbers in the figure indicate the length of each block and are in the unit of nm. Figure 5 2 Transmission coefficients for LA and LO branches as a function of incident phonon frequency. Points on the left of the dashed line belong to the LA branch while those on the right belong to LO branch. Lines are guides to an eye.


7 5 Figure 5 3 Transmission coefficients for all modes along k 00 direction for l sio2 =0.8 nm structure. The dashed line is the fit discussed in the text. Figure 5 4 LA incident phonon transmission for l SiO2 =0.8 nm (squares) and l SiO2 =10.9 nm (circles). 0.0 0.2 0.4 0.6 0.8 1.0 0 5 10 15 20 Transmission coefficient Frequency [THz] LA TA LO TO fit 0.0 0.2 0.4 0.6 0.8 1.0 0 2 4 6 8 10 12 14 Energy transmission Frequency [THz] LA(SiO2=10.9nm) LA(SiO2=0.8nm)


76 Figure 5 5 Relative transmitted energy as a function of k x and k y for four different incident phonon frequencies. All shown are for k 00 LA incident phonon on l SiO2 =0.8 nm. Note that the scal es are different in each figure and all kx and ky 0 Figure 5 6 Wavevector deviation as a function of incident phonon frequency for transmission and reflection phonons . 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 2 4 6 8 10 12 14 d k [2 /a 0 ] Frequency [THz] Transmission Reflection


77 Figure 5 7 Frequency distribution of wave packets for l SiO2 =0.8 nm structure. All incident phonons are in pure k 00 LA modes. All energies shown are relative and scaled to the same incident energy.


78 Figure 5 8 Contribution to the thermal conductivity as a function of k x and k y , with summation of all k z contribution. Data are fit by interpolation and line integration at each (k x , k y ). The region inside the white triangle indicates the region of the calculation, while the rest is plotted by exploiting the symmetry of the first Brillouin Zone


79 Figure 5 9 Extracting the Kapitza resistance of the Si/SiO 2 interface at 300 K. The extrapolated values at l SiO2 =0 are 2.96X10 9 and 2.74X10 9 m 2 K/W for NEMD and PWD respectively, which are twice the Kapitza resistanc e obtained from Eq. 3 14 . 0 2 4 6 8 10 0 5 10 15 20 25 30 Total Thermal Resistance [X10 9 W 1 Km 2 ] SiO 2 thickness l SiO2 [nm] NEMD PWD


80 Figure 5 10 Contribution to Kapitza conductance by mode for l SiO2 =0.8 nm . Both TA and TO include the contribution of two transverse mode s. The mode wise contribution is almost the same for l SiO2 =2.2 nm structure.


81 CHAPTER 6 THERMAL TRANSPORT AT A TWIST GRAIN BOUNDARY IN UO 2 In this chapter, the phonon wave packet dynamics method introduced in Chapter 3 is used to study the phonon transmission at a UO 2 twist grain boundary. As mentioned earlier, grain boundaries, as planar defects in crystal structures, act as obstacles to thermal transport as they break the periodicity of the crystal in the direction perpendicular to the boundaries. Introduction O n the experimental side, understanding the Kapitza resistance of a grain boundary remains challenging, because various microstructur al elements ; thus more comprehensive understanding of phonon transmi ssion through grain boundary is needed [ 2 , 76 ] . Consequently, most experiments calculate the a representative value for the Kapitza resistance of a n ensemble of grain boundar ies by averaging the effects from the thermal conductivity of a polycrystalline material [ 76 ] . However, by utilizing atomistic simulation methods, the effect of individual gra in boundar ies on ph onon transport can be dissected. A twist grain boundary, as its name suggests, consists a rotation of two grains about an axis that is perpendicular to the grain boundary. An idealization of the formation of a grain boundary is illustrat ed in Figure 6 1 . First, the crystal is cut into two grains, and then one grain is rotated by an angle with respect to the other ; the two grains are then stuck back t ogether. In reality, the rotation angle can be any value. However, in order to apply periodic boundary condition and have a computationally accessible system size, GB simulations are generally limited to boundaries which have a finite number of lattice sit es from the two grains that coincide, resulting in a finite unit


82 cross sectional area. Such structures are known as coincident site lattice (CSL) GBs [ 77 ] . The CSL concept is solely based on the geometry of the lattice, with the rotation angles limited to those that are able to bring two sites into coincidence, no matter how the atomic positions are configured. The rotation angle can be determined based on the lattice geometry. As shown in Figure 6 2 , the grey grid refers to a layer of a cubic lattice, with each small square being a unit cell. The dashed blue and red sq uare represent the two grains, and after the rotation indicates by the arrows, the periodic cell of the grain boundary is obtained as shown in the figure on the right. The four vertices of the new unit cell are the coincident sites, and by taking the perio dic condition into account, the new cell has one equation: ( 6 1 ) values shown in Figure 6 2 . For this particular example, the new cross twist grain bo to the original cell and it can be determined by the square sum of x and y. If x 2 +y 2 is an x 2 +y 2 (x 2 +y 2 )/2, for 2 coincident sites can be found in the new cell. ° is chosen ° . And both two struc ture has similar size in the cross sectional plane, being 24.46 Å X 24.46 Å Å X 22.54 Å . Such small


83 2 [ 78 ] , as shown in Figure 6 3 ; moreover, the small cross sectional area makes the computations quite rapid. Also, due to the four fold rotation axis along (001) for UO 2 (space group 225, see Chapter 2), the rotation angles are limited below 45 ° beyond 45 ° towards one direction is equivalent to rotating 90 Kapitza Resistance by NEMD Method The direct method [ 36 ] is used to calculate the Kapitza resistance due to the twist grain boundary. The simulation settings are similar to th ose used in Chapter 4, but with grain boundaries residing in place of dislocations, as indicated by the upper figure in Figure 6 4 . The two grain boundaries (required by periodic boundary condition along L z ) are located at 1/4 L z and 3/4 L z , and the two thin layer at the edge and center of L z act as the heat source and sink, thus creating the thermal current desired. The whole simulation box is divided into 64 equally thic k sampling regions, each region contains more than 200 atoms. The temperature profile is an average of the temperature in each region over 100 ns . temperature s of 400, 800, 120 0 and 1600 K, as shown in Figure 6 5 . Over the entire This may due to the fact a larger lattice mismatch. [55] Also, a temperature dependence of the Kapitza resistance is observed, with the resistance dropping as temperature increases. This trend is opposite that of the bulk thermal c onductivity of UO 2 in Chapter 4 ( Figure 4 3 ). This suggests that, anharmonicity , to at least some degree, facilitate s thermal transport at the grain boundary. Moreover, such a trend has also been observed in several other systems


84 [ 7 9 ] . However, it should be noted here that the decreasing trend is not uniform, and there is some error bar overlaps in some temperature ranges. Phonon Transmission by PWD Method To understand the phonon scatterings at the grain boundaries in more detail , the phonon wave packet dynamics method is employed to study how phonons with a single frequency and symmetry interact with a grain boundary. To prepare the structure for PWD simulations, the two grains of the GB structure are generated and rotated accord ing to the scheme described in the previous section ( Figure 6 2 ). The twist grain boundary is on the (001) plane, and each grain is 400 UO 2 unit cells in length along the direction perpendicular to the grain boundary. Such a large domain in this particular direction is to ensure that the wave packet can and propagate over an extended distance. The initial GB structure is first heated up to 2000 K and slowly cooled to 0K . A subsequent quench is performed so that the force on each atom is very small, less than 10 8 eV/Ã… per atom . This eliminates any excess structural energy that could affect the phonon scattering simulation. The interatomic potential used for this study i s Busker potential; as mentioned in Chapter 4, this potential has previously been used in studies of the effects of point defects [ 80 ] , dislocations [ 60 ] and grain size [ 7 ] on thermal conductivity. To perform the PWD simulation, the phonon dispersion relation must be obtained first. The phonon dispe rsion is calculated using GULP [ 33 ] along (001) direction. Although there is not an exact match with experiment [ 81 ] , the Busker potential does give a reasonable description of the acoustic phonon modes .


85 Preparing Wave Packet As introduced in Chapter 3, to generate the wave packet and perform the normal mode analysis, the available wavevectors need to be identified. This is trivial work for a single crystals whose dimensions are simply integer multiples of its unit cell. However, in a grain boundary system like in Figure 6 2 , where a rotation around the z axis has been performed, the wave vector should also be rotated accordingly [ 82 ] . The reciprocal lattice vector on the cross sectional plane for one grain boundary supercell is and , where L is the size of the supercell in x or y. Then the two vectors were multiplied by integers and a + ed. After the procedure, those combinations of which are located inside the first Brillouin zone of a UO 2 primitive cell are the wavevectors needed. The wave packets are constructed according to the procedure described in detail in Chapter 3. Unlike genera ting a wave packet in a system where there is only one atom, in UO 2 primitive cell, there are three atoms: one U and two O. Thus, there are a total of 9 eigenfrequency states, resulting in 81 components in the eigenvectors. Thus, when calculating the phono n dispersion curve along a high symmetry direction in reciprocal space, branches can cross each other; along a low symmetry direction, the irreducible representations of the symmetry group of the Hamiltonian. [ 83 ] Thus, it is crucial to differentiate the 9 different energy states , since GULP outputs energies in ascending order. To achieve this, the dispersion curve of each mode is distinguished by maintaining the smoothness of the dispersion curve. The second order derivatives of the eigenfrequency with respect to k z are used to detect any mo de crossing or avoided


86 crossing by maintaining the continuity of the curve. An example of th is analysis results is shown in Figure 6 6 . Here, the wave packet is initi ally located at the center of one grain, and is launch ed toward the grain boundary. When the wave packet reaches the grain boundary region, scattering events take place: part of the wave packet is transmitted through the boundary while some is reflected ba ck. By determining the energies of initial energy in the wave packet that transmitted, is obtained as a function of the input phonon frequency. Phonon Scatterings at G B To study the phonon scattering at the grain boundary, phonon wave packets 17 grain boundaries. Their transmission coefficients are shown in Figure 6 7 . It is clear that on incident phonon frequency. In particular, at low frequencies, corresponding to l ong wavelengths, the transmission is nearly 100%. This high transmission coefficient arises because the thickness of the grain boundary along the phonon propagation direction is much smaller than the wavelength. Such high transmission coefficients for low frequency LA modes have been seen for a number of other systems. [ 84 87 ] However, as the phonon frequency increases, the transmission coefficient gradually drops to very low value. To better illustrate the different scattering scenarios of low and high frequency phonons, Figure 6 7 illustrates a series of snapshots of two cases of different phono n


87 propagating through the grain boundary . The low frequency case is the 0.63 THz LA wave packet, corresponding to an average wavelength about 109 Ã…, which is much greater than the dimension of the grain boundary . As Figure 6 8 (a) shows the phonon wave packet passes through then GB (denoted by the vertical line in the center of each panel) essentially unscattered. This is consistent with Figure 6 7 , which shows that transmission coefficient for this phonon wave packet is close to unity. The high frequency cases is a 4.19 THz LA wave packet, with average wavelength of about 15.62 Ã… . Here there is strong scattering at the GB, with part of the wave packet being reflected back. Moreover, there is significant mode coupling , with the generation of TA phonons, as shown in green. However, for the case of the LO mode, when the frequency is ab ove 8 THz, the majority of the energy transmission comes from modes other than the initial wavepacket, indicating much stronger scattering. A similar snapshots analysis is shown for the case of LO mode in Figure 6 9 . To develop a more detailed understanding of the contributions from different vibrational modes to the energy transmission, normal mode analysis (see Chapter 3) is applied to the transmitted and reflected wave packet respectively. Figure 6 10 shows the fraction of energy transmitted through the grain boundary that still features the same vibrational proper ties as the incident wave packet. For the LA mode, when the frequency is less than 3 THz, the wave packet is hardly scattere d by the grain boundary at all. When the incident LA frequency is between 3 THz and 6.7 THz, a small amount of the energy is scatter ed into TA modes, as was clearly identified in Figure 6 8 (b). From 6.7 THz up to the maximum frequency of LA mode, the phonon is scattered into TO modes since TA mode s are no longer available . Moreover, mode conversion only


88 contributes a small fraction o f the transmitted energy for all LA modes . A s shown in Figure 6 11 , the amount o f energy transmitted and reflected due to the inelastic scatterings is almost the same, suggesting a typical DMM [ 2 ] behavior. From this, the snapshots of Figure 6 8 (b) c an be better understood: a large amount of the wave packet has been elastically scattered, where the original phonon vibrational mode is preserved; and part of the energy excited the vibr ational modes in the GB, and then reemission equally back to the left and right grains. Conclusion This work presented has shown that, by using NEMD simulation, the Kapitza conductance s of ies can be directly determined. Of th e two grain boundaries studied, the one with larger rotation angle and grain boundary energy has larger thermal resistance. Also, the Kapitza resistance tend s to drop with increasing temperature, suggesting anharmonicity helps with thermal transport. This is the first time that the PWD method has been applied to a system with more than one atomic species; moreover, it is an ionic system. The PWD analysis shows the detailed scattering picture at the grain boundary of well defined phonon wave packets. The pho non energy transmission coefficient is close to unity for LA mode at small incident frequencies; while it drops to a very small value with increasing frequencies. Also, the mode conversion, due to inelastic scatterings, only has limited contribution to the LA mode transport, but has significa nt contribution to LO mode transport.


89 Figure 6 1 Schematic of a twist grain boundary. Figure 6 2 Schematic illustra tion of determining the rotation angle for CSL of 17 twist grain boundary (cubic lattice).


90 Figure 6 3 CSL distribution of uranium dioxide sample. From reference [ 78 ] . Figure 6 4 Simulation configuration and temperature profile for calculating Kapitza resistance. The two GBs are located at L z =87 Ã… and L z =262 Ã… respectively.


91 Figure 6 5 boundaries in UO 2 . Figure 6 6 Avoided crossing for phonon dispersion relation in UO 2 . Here k x =0.1 and k y =0.2.


92 Figure 6 7 Phonon energy transmission coefficient (LA and LO) fo boundary as a function of phonon frequency. Blue color squares denote the LA mode and orange triangles denote the LO mode. The inset shows the same plot of transmission coefficient but as a function of wavevector.


93 Figure 6 8 Time evolution of atomic displacement with LA incident phonon wave packet indicate x and y displacement. The GB is i ndicated by the dashed li ne. A) The case of k z =0.15*2 /a 0 ; B) the case of k z =0.35*2 /a 0 .


94 Figure 6 9 Time evolution of atomic displacement with LO incident phonon wave packet indicate x and y displacement. The GB is indicated by the dashed line. (a) Shows the case of 0.63 THz fre quency; (b) shows the case of 4.19 THz frequency.


95 Figure 6 10 Contribution of the un scattered mode to the energy transmission coefficient. The dashed line is an indication of the total transmission coef ficient shown in Figure 6 7 . The inset shows the same plot in wavevector space. Figure 6 11 Phonon energy contribution from inelastic scattering. 0.00 0.05 0.10 0.15 0.20 0 2 4 6 8 10 Percentage Frequency [THz] Reflection Transmission


96 CHAPTER 7 LATTICE DYNMICS MODEL FOR INTERFACIAL THERMAL TRANSPORT In this chapter, an atomistic lattice dynamics (LD) approach is used to study behavior of interfacial thermal transport at a phenomenological le vel. First, a one dimensional monoatomic chain model is used illustrate the mechanisms of the method and how some key variables in the model can affect the energy transmission. Then the one dimensional diatomic chain model is used to illustrate the transmi ssion of longitudinal acoustic and longitudinal optical vibration modes at an interface. Introduction A general outline of the model is shown in Figure 7 1 . Two leads of material are join ed together with an interface (or any linking material) presented in the center, and as indicated by the arrow, the propagation direction for the energy is assumed to be from right to left. The model assumes that the left and right lead infinite [ 88 ] , such that the motion of the atoms belonging to the leads can be described using the wave representation of atomic displacement. That is, the displacement of the atoms in t he right lead is the displacement from the incident vibrational mode, plus that from the reflected vibrational modes; and the displacement of the atoms in the left leads is just the sum from the transmitted vibrational modes, as expressed in Eq. 7 1 and 7 2 [ 89 91 ] , ( 7 1 ) ( 7 2 ) where on the left term is the total displacement of the i th atom in the l th unit cell on the left or right. 0 ,k 0 ),


97 s the incident vibration mode, the vibration mode resulting from R e T r are the coefficient s of amplitude reflection and coefficient of amplitude transmission. The underlying crystal symmetry is lost for atoms in the interface region; thus, the wave representation is no longer a suitable descriptor and their motions are better described by a series of dynamical equations of motion. An example of a one dimensional form of the lattice dynamics equation with the assumption of nearest neighbor interactions only will be illustrated in the next section . Many lattice dynamics approache s have been used to study the energy transmission at an interface. [ 92 95 ] Wang [ 89 ] have generalized these special cases and proposed the scattering boundary method within the lattice dynamics approach. Here, a similar method is applied with the help of MATLAB ® symbolic function capabilites, and the process is summarized in Figure 7 2 . The main focus is solving for the coefficients of amplitude transmission/reflection. The total energy reflection coefficient and transmission coefficient of an 0 ,k 0 ) can th en be calculated [ 89 ] as , ( 7 3 ) ( 7 4 ) w here v g denotes the group velocity of a particular phonon. With these two coefficients, the interfacial thermal conductance can the n be calculated as (see Eq. 3 18):


98 ( 7 5 ) One Dimensional Monatomic Chain There have been several studies on the one dimensional model previously [ 92 , 96 , 97 ] . As used here , the model aims to simulate the situation at a simple interface and validate the procedure. A schematic of the one dimensional monoatomic chain is shown in Figure 7 3 (a). The atoms are numbered with integer N ascending from left to right where the left lead (N< 1) and right lead (N>0) contains a large number of atoms, and the interface region has two atoms (N= 1, 0). Each atom is connected to its nearest neighbor locat ed on the left and right by a linear spring. Here, all the atoms are assumed to have the same atomic mass, that m 1 =m 2 , a nd the spring constant in both lead s are assume d to be the same, f 1 =f 2 . The only different spring constant is f 3 , which represents the irregular bonding at the interface. From the simple nature of the model, resulting from these assumptions, the dynamics of the system (atomic displacement of the whole system) can be obtained and the amplitude transmission and reflection coefficient s can b e determined by following the flow chart in previous section, as described in the following. (The physics of these can be found in many introductory books on solid state physics, such as [ 16 ] ) In a semi infinite harmonic one dimensional monoatomic system, such as the left and right leads, the wave representation of an atomic displacement is given by, ( 7 6 )


99 is the phonon frequency. The lead atoms next to the interfaces are N= 2 and N=1, and their displacement can be written according to Eqs. 7 1 , 7 2 and 7 6 , ( 7 7 ) ( 7 8 ) . 7 1 ~ 7 4 are not present, since there is only one vibration mode (LA), and the wavevector can only be along left or right. In a harmonic system, the dynamic equation of the system can always be given as, ( 7 9 ) where f L and f R denotes the spring constant of the left and right spring of atom N. Eq. 7 9 can be rewritten as, ( 7 10 ) ( 7 11 ) Clearly, u N 1 /u N+1 can be solved by knowing u N and u N+1 /u N 1 . Here, the operation is defined as S L and S R as expressed in the equations , which are carried out by MATLAB ® . Thus, after a few iterations, u 1 and u 2 can be obtained as, ( 7 12 ) ( 7 13 )

PAGE 100

100 Then the last step is inserting Eqs. 7 7 and 7 8 to the results of 7 10 and 7 11 . The amplitude transmission coefficient can be solved, as well as the total energy transmission coefficient. As a test case, the mass of the atoms are set to m=1, and the spring constant inside the leads are set to k=1. The total energy transmission is plot in Figure 7 4 as a function of wavevector for different k 3 . The dashed line shows a case of k 3 =1 as a validation, which suggests a homogeneous system, the transmission coefficient is unity for all wavevectors. For all cases when k 3 , the transmission coefficient is unity at small wavevectors, which agrees with the long wavelength limit [ 16 ] , that the phonons are little disturbed by the interface. As the discrepancy between the bonding at the interface and the bonding at the interface increases, the total energy transmissi on decreases rapidly. In some cases, the interface may actually be an extended region, as the atomic bindings next to the interface are also altered. As sketched in Figure 7 3 (b), the extended interface region is simulated by assigning more than one irregular spring. Figure 7 5 summarizes the cases with 1~3 irregular springs. One feature arises from the multiple irregular spring configuration is the oscillation in the energy transmission coefficient, which is the result of multiple scatterings [ 98 , 99 ] . One Dimensional Diatomic Chain In order to study the transmission behavior of optic al modes, the model is extended to a one dimensional diatomic chain. As shown in Figure 7 6 (a), each unit cell now contains two atoms m 1 and m 2 (assuming m 1 2 ), and assuming m 1 and m 2 are equally spaced and only the nearest neighbors interactions are present. In order to be somewhat analogous to a UO 2 structure as indicated by Figure 7 6 (b), the atomic

PAGE 101

101 masses are assumed to be m 1 =238.03 and m 2 =32.00 (representing both oxygen atoms). The spring constant in the lead is properly chosen, so that the (001) dispersion relation of LA mode in the model matches that in the UO 2 with Busker interatomic potential (Chapter 3). As shown in Figure 7 7 , the two LA modes matches very well, but the LO mode of the model cannot be fit UO 2 dispersion well at the sa me time, as there are three atoms in a UO 2 primitive cell. The energy transmission from the model is fitted to the LA transmission from PWD data. Here, the assumption that the GB only disturbs the bonding locally is made. So the only fitting parameter here is the spring constant in the center (f 3 as in Figure 7 6 (a)), which represents the irregular bonding condition at the UO 2 grain boundary. F itting is performed using a least squares process [ 100 ] . As illustrated in Figure 7 8 , the model generates the correct general trend of the energy transmission of the two modes. In particular, the match of the LA mode is rather good, but the LO mode is not. This is because that for the LA mode, all atoms move in phase, thus the model here is more suitable for the LA plane wave; for the LO mode, since there are three atoms in a UO 2 primitive cell, the one dimensional representation is not sophisticated enough . A simple parametric study is carried out, where the fitting parameter, f 3 , is increased/decreased by 20% and 40% respectively. As a result, as shown in the inset of Figure 7 8 , the LA transmission curve shifts outward when f 3 is more close to the spring constant in the bulk, indicating less lattice mismatch; the curve shifts inward when f 3 is more deviate from the bulk spring constant, indicating more lattice mismatch . Conclusions This work illustrated the energy transmission at a GB using a simple chain model with lattice dynamics. The match with the PWD data for the longitudinal acoustic mode

PAGE 102

102 illustrated its potential for use in the study more complex systems. The mo del is able to provide a fast way to estimate the energy transmission coefficients, and therefore could be used in higher scale simulations. In order to include transverse modes (TA and TO), a higher dimensionality in the lattice is required. In this way, mode conversion could be included, which was seen to take place in the PWD simulations in Chapters 5 and 6. Moreover , a three mass model may help improve the optical mode behavior in this model. In addition, a framework of the coding was written in MATLAB ® , and with a structure that facilitates the introduction of more complex physics.

PAGE 103

103 Figure 7 1 Schematic of the model set up . Figure 7 2 Flow chart for solving LD model.

PAGE 104

104 A B Figure 7 3 One dimensional monoatomic chain model. A) Single spring model. B) Extended spring model. Figure 7 4 Energy transmission as a function of f 3 for the monatomic chain .

PAGE 105

105 Figure 7 5 Energy transmission for various extended interface s. A (b) Figure 7 6 Illustration of diatomic model. A) One dime nsional diatomic chain model. B ) M apping to th e UO 2 structure.

PAGE 106

106 Figure 7 7 Phonon dispersion relation for the diatomic model. Blue dashed line shows the LA dispersion from UO 2 . Red and purple dashed lines shows the two LO dispersion from UO 2 . Figure 7 8 Energy transmission of diatomic model as a function of in cident phonon frequency (lines); the data points denote the PWD results . The inset shows the transmission change as a function of f 3 variation.

PAGE 107

107 CHAPTER 8 CONCLUSION AND OUTLOOK In this work, investigations o f how defect structures affect the thermal transport of UO 2 and Si have been carried out using atomistic simulation techniques. The effects of UO 2 thermal transport drop due to the presence of edge dislocations was quantified and illustrated in Chapter 4. The temperature and dislocation density dependence of the reduction was also established. In addition, attempts were made to extract parameters of empirical model for UO 2 fuel thermal transport from the MD simulation results. In Chapter 5, the Kapitza resistance of a Si/SiO 2 interface from a Si/SiO2/Si hetero system was calculated. The phonon interface scatterings from individual phonon mode s , sampled from the first Brillouin zone, were investigated and the contribution to the interfacial resistance from each phonon mode was also determined . The Kapitza 2 twist grain boundaries were determined in Chapter 6 . In addition, the phonon scatterings at the GB were also analyzed, where different behaviors of LA and LO phono ns are recognized . A MATLAB ® code for linear lattice dynamics chain model has been developed in Chapter 7. Th is one dimensional model is capable of studying longitudinal phonon energy transmission at simple interfaces and has potential for the application in more complex systems. The MD, PWD and LD methods were all applied to study the interfacial thermal transport in some way or another in this work. The molecular dynamics method is capable of calculating thermal transport in materials where phonons are th e main heat carrier. It can directly predict the effects on thermal conductivity reduction due to various defects at different temperatures in relative large system scales. However, post -

PAGE 108

108 simulation corrections for the results are sometimes needed at elevat ed temperatures in order to match experimental data, due to the anharmonic effects and size effects of the MD simulation. The phonon wave packet dynamics is an attractive method to study the phonon transport at various defect structures. This method provi des more than just a number, as the direct MD method does. There are a few issues relating the method that should be noted. The whole simulation is performed at close to zero temperature, for the system is deeply quenched to minimize any thermally activat ed phonons, and the energy contained in the wave packet is extremely small. Conversely, if the above requirement is not satisfied, for instance, a system running at finite temperature, the scatterings between the wave packet pulse and other thermal lattice vibrations would be so strong, that it would be nearly impossible to distinguish the wave packet signal. Therefore, when calculating the Kapitza resistance or thermal conductivity from the PWD results, all the temperature effect come into play with the in troduction of Bose Einstein statistics; that is, the phonon phonon interactions are neglected. The development of a method that can look at individual phonon interface scattering events in the presence of a thermal distribution of phonons remains a challen ge for the future Another con straint on the method is computational cost. Firstly, to provide a complete phonon transport picture, a sampling of the first Brillouin Zone is needed. Some phonon modes have very small group velocities, requiring rather long s imulation time. Secondly, a rather large dimension along the direction of wave packet propagation is required. On the one hand, this requirement assures that the wave packet can be

PAGE 109

109 localized On the other hand, when generating the initial wave packet or ana lyzing the final system displacement, the phonon data of the bulk crystal is incorporated, all wave packet signal should keep an adequate distance from the defects, where structural disorder is likely to exist. This long dimension requirement create comput ational obstacles when studying a structure where a moderate cross sectional area is needed, for example, dislocations. The lattice dynamics model, which provides a rapid way to investigate and modify some key variables in the real system regarding interfa cial thermal transport. The model is expected to become much more complicated when studying a more real istic model. In fact, the mechanism of this model, in some ways, is very similar to the phonon wave packet dynamics, if, detailed atomistic bonding as t he PWD is taking into account in the model and same interatomic potential is used (if linearization of the potential is possible) . This work has also provided some directions of future work. Firstly, as various effects of defects on thermal conductivity have been identified individually by some previous work, as well as this one, it would be desirable to analyze complex effects arising from multiple scattering processes. Secondly, it would be interesting to introduce the phonon phonon interactions at the interface when using the PWD method. Thirdly, the model could be extended to more variables, for example, including longer range interactions, anharmonic effects. Fourthly, it would be useful to compare the results of MD, PWD and LD model for a same syste m setting (e.g., same structure and potential), and extract the underlying physics from the differences. Finally, it would be

PAGE 110

110 valuable to systematically investigate possible means to pass the results from the atomic level simulations to larger scale models , such as mesoscale models.

PAGE 111

111 LIST OF REFERENCES 1. Schaller, R.R., Moore's law: past, present and future. Spectrum, IEEE, 1997. 34 (6): p. 52 59. 2. Cahill, D.G., W.K. Ford, K.E. Goodson, G.D. Mahan, et al., Nanoscale thermal transport. Journal of Applied Physics, 2003. 93 (2): p. 793 818. 3. Tulenko, J.S., Nuclear Reactor Materials and Fuels. Nuclear Energy, 2013: p. 203 214. 4. Richerson, D., Modern Ceramic Engineering: Properties, Processing, and Use in Design, Third Edition . 2005: Tayl or & Francis. 5. Fink, J.K., Thermophysical properties of uranium dioxide. Journal of Nuclear Materials, 2000. 279 (1): p. 1 18. 6. Phillpot, S.R., A. El Azab, A. Chernatynskiy, and J.S. Tulenko, Thermal Conductivity of UO 2 Fuel: Predicting Fuel Performance from Simulation. Jom, 2011. 63 (8): p. 77 83. 7. Watanabe, T., S.B. Sinnott, J.S. Tulenko, R.W. Grimes, et al., Thermal transport properties of uranium dioxide by molecular dynamics simulations. Journal of Nuclear Materials, 2008. 375 (3): p. 388 396. 8. Celler, G.K. and S. Cristoloveanu, Frontiers of silicon on insulator. Journal of Applied Physics, 2003. 93 (9): p. 4955 4978. 9. Kasap, S.O., Principles of electronic materials and devices . 2006: McGraw Hill. 10. Stillinger, F.H. and T.A. Weber, Computer si mulation of local order in condensed phases of silicon. Phys Rev B Condens Matter, 1985. 31 (8): p. 5262 5271. 11. Watanabe, T., H. Fujiwara, H. Noguchi, T. Hoshino, et al., Novel interatomic potential energy function for Si, O mixed systems. Japanese Journ al of Applied Physics Part 2 Letters, 1999. 38 (4A): p. L366 L369. 12. V., B.V. and C. Wei, Computer Simulation of Dislocations . Oxford Series on Materials Modeling. 2006: Oxford University Press. 13. CrystalMaker . 2014, CrystalMaker Software Ltd.,: Oxfordshire, UK. 14. Allen, S.M. and E.L. Thomas, The structure of materials . MIT series in materials science and engineering. 1999, New York: John Wiley. 15. Phillpot, S.R. and A.J.H. McGaughey, Introduction to thermal transport. Materials Today, 2005. 8 ( 6): p. 18 20.

PAGE 112

112 16. Kantorovich, L., Quantum Theory of the Solid State: An Introduction . 2004: Kluwer Academic Publishers. 17. Ziman, J.M., Electrons and Phonons: The Theory of Transport Phenomena in Solids . 1960: Clarendon Press. 18. Callaway, J., Quantum t heory of the solid state . 1974: Academic Press. 19. Ashcroft, N.W. and N.D. Mermin, Solid state physics . 1976: Saunders College. 20. Weisstein, E.W. Fourier Transform . Available from: http: // . 21. Weisstein, E.W. Discrete Fourier Transform . Available from: . 22. Cahill , D.G., P.V. Braun, G. Chen, D.R. Clarke, et al., Nanoscale thermal transport. II. 2003 2012. Applied Physics Reviews, 2014. 1 (1): p. 011305 23. Swartz, E.T. and R.O. Pohl, Thermal boundary resistance. Reviews of Modern Physics, 1989. 61 (3): p. 605 668. 24 . Schelling, P.K., S.R. Phillpot, and P. Keblinski, Kapitza conductance and phonon scattering at grain boundaries by simulation. Journal of Applied Physics, 2004. 95 (11): p. 6082 6091. 25. Govers, K., S. Lemehov, M. Hou, and M. Verwerft, Comparison of Inte ratomic Potentials for UO 2 . Part I: Static Calculations. Journal of Nuclear Materials, 2007. 366 (1 2): p. 161 177. 26. Govers, K., S. Lemehov, M. Hou, and M. Verwerft, Comparison of Interatomic Potentials for UO 2 Part II: Molecular Dynamics Simulations. Jo urnal of Nuclear Materials, 2008. 376 (1): p. 66 77. 27. G., B., Solution and Migration of Impurity Ions in UO 2 , in Department of Materials . 2002, Imperical College of Science, Technology and Medicine. 28. Wolf, D., P. Keblinski, S.R. Phillpot, and J. Eggeb recht, Exact Method for the Simulation of Coulombic Systems by Spherically Truncated, Pairwise r 1 Summation. Journal of Chemical Physics, 1999. 110 (17): p. 8254 8282. 29. Watanabe, T. and I. Ohdomari, Modeling of SiO2/Si(100) interface structure by using extended Stillinger Weber potential. Thin Solid Films, 1999. 343 : p. 370 373. 30. Watanabe, T., K. Tatsumura, and I. Ohdomari, SiO(2)/Si interface structure and its formation studied by large scale molecular dynamics simulation. Applied Surface Science, 2 004. 237 (1 4): p. 125 133.

PAGE 113

113 31. Broido, D.A., A. Ward, and N. Mingo, Lattice thermal conductivity of silicon from empirical interatomic potentials. Physical Review B, 2005. 72 (1). 32. Chernatynskiy, A., C. Flint, S. Sinnott, and S. Phillpot, Critical assess ment of UO<sub>2</sub> classical potentials for thermal conductivity calculations. Journal of Materials Science: p. 1 10. 33. Gale, J.D., GULP: A Computer Program for the Symmetry Adapted Simulation of Solids. Journal of the Chemical Society Fa raday Transactions, 1997. 93 (4): p. 629 637. 34. Gale, J.D. and A.L. Rohl, The General Utility Lattice Program (GULP). Molecular Simulation, 2003. 29 (5): p. 291 341. 35. IAEA, Thermophysical Properties Database of Materials for Light Water Reactor and Heav y Water Reactors . www 36. Jund, P. and R. Jullien, Molecular dynamics calculation of the thermal conductivity of vitreous silica. Physical Review B, 1999. 59 (21): p. 13707 13711. 37. Schelling, P.K., S.R. Phillpot, and P. Keblinski, Phonon wa ve packet dynamics at semiconductor interfaces by molecular dynamics simulation. Applied Physics Letters, 2002. 80 (14): p. 2484 2486. 38. Kimmer, C., S. Aubry, A. Skye, and P.K. Schelling, Scattering of phonons from a high energy grain boundary in silicon: Dependence on angle of incidence. Physical Review B, 2007. 75 (14). 39. Plimpton, S., Fast Parallel Algorithms for Short Range Molecular Dynamics. Journal of Computational Physics, 1995. 117 (1): p. 1 19. 40. MATLAB . 2014, The MathWorks Inc.,: Natick, Massa chusetts. 41. Ju, L., AtomEye: an efficient atomistic configuration viewer. Modelling and Simulation in Materials Science and Engineering, 2003. 11 (2): p. 173. 42. Dolling, G. Inelastic scattering neutrons in solids and liquids . in Proc. Symp. 2nd Chalk River . 1962(1963). Can. 43. Deng, B., A. Chernatynskiy, P. Shukla, S.B. Sinnott, et al., Effects of edge dislocations on thermal transport in UO2. Journal of Nuclear Materials, 2013. 434 (1 3): p. 203 209. 44. Evans, A.G. and P.L. Pratt, Dislocations in The Fluorite Structure. Philosophical Magazine, 1970. 21 (171): p. 645 645. 45. Sawbridge, P.T. and E.C. Sykes, Electrostatic charges on dislocations in uranium dioxide. Journal of Nuclear Materials, 1970. 35 (1): p. 122 125.

PAGE 114

114 46. Ashbee, K.H.G., Stacking faults in uranium dioxide. Proceedings of the Royal Society of London Series a Mathematical and Physical Sciences, 1964. 280 (1380): p. 37 46. 47. Nogita, K. and K. Une, Radiation induced microstructural change in high burnup UO 2 fuel pellets. Nuclear Instruments & Methods in Physics Research Section B Beam Interactions with Materials and Atoms, 1994. 91 (1 4): p. 301 306. 48. Une, K., K. Nogita, S. Kashibe, and M. Imamura, Microstructural change and its influence on fission gas release in high burnup UO2 fuel , in Nuclear Materials for Fission Reactors , H. Matzke and G. Schumacher, Editors. 1992, Elsevier: Oxford. p. 65 72. 49. Ham, R.K., The Determination of Dislocation Densities in Thin Films. Philosophical Magazine, 1961. 6 (69): p. 1183 1184. 50. Veshchunov, M.S. and V.E. Shestak, Model for evolution of crystal defects in UO 2 under irradiation up to high burn ups. Journal of Nuclear Materials, 2009. 384 (1): p. 12 18. 51. Geelhood, K.J. and W.G. Luscher, FRAPCON 3.5: A Computer Code for the Calculation of Steady State, Thermal Mechanical Behavior of Oxide Fuel Rods for High Burnup , P.N.N. Laboratory, Editor. 2014: NUREG/CR 7022. 52. Geelhood, K.J., W.G. Luscher, and C.E. Beyer, FRAPCON 3.4: A Computer Code for the Calculation of Steady State Thermal Mechani cal Behavior of Oxide Fuel Rods for High Burnup . 53. Klemens, P.G., THERMAL CONDUCTIVITY AND LATTICE VIBRATIONAL MODES. Solid State Physics Advances in Research and Applications, 1958. 7 : p. 1 98. 54. Klemens, P.G., The scattering of low frequency lattice waves by static imperfections. Proceedings of the Physical Society of London Section A, 1955. 68 (12): p. 1113 1128. 55. Callaway, J., Model for lattice thermal conductivity at low temperatures. Physical Review, 1959. 113 (4): p. 1046 1051. 56. Zou, J., D. K otchetkov, A.A. Balandin, D.I. Florescu, et al., Thermal conductivity of GaN films: Effects of impurities and dislocations. Journal of Applied Physics, 2002. 92 (5): p. 2534 2539. 57. Schelling, P.K., S.R. Phillpot, and P. Keblinski, Comparison of Atomic Le vel Simulation Methods for Computing Thermal Conductivity. Physical Review B, 2002. 65 (14): p. 144306.

PAGE 115

115 58. Hu, M., P. Keblinski, J. S. Wang, and N. Raravikar, Interfacial thermal conductance between silicon and a vertical carbon nanotube. Journal of Applie d Physics, 2008. 104 (8): p. 083503 . 59. Vega, D.A., Atomistically informed fuel performance codes: a proof of principle using molecular dynamics and FRAPCON simulation , in Nuclear and Radiological Engineering . 2008, University of Florida: Gainesville, Fla. 60. Deng, B., M. Khafizov, D.H. Hurley, and S.R. Phillpot, Kapitza resistance of Si/SiO2 interface. Journal of Applied Physics, 2014. 115 (8): p. 084910. 61. Lee, S. M. and D.G. Cahill, Heat transport in thin dielectric films. Journal of Applied Physics, 1 997. 81 (6): p. 2590 2595. 62. Kato, R. and I. Hatta, Thermal Conductivity and Interfacial Thermal Resistance: Measurements of Thermally Oxidized SiO2 Films on a Silicon Wafer Using a Thermo Reflectance Technique. International Journal of Thermophysics, 200 8. 29 (6): p. 2062 2071. 63. Hurley, D.H., M. Khafizov, and S.L. Shinde, Measurement of the Kapitza resistance across a bicrystal interface. Journal of Applied Physics, 2011. 109 (8). 64. Mahajan, S.S., G. Subbarayan, and B.G. Sammakia. Estimating Kapitza re sistance between Si SiO 2 interface using molecular dynamics simulations . in Thermal and Thermomechanical Phenomena in Electronic Systems, 2008. ITHERM 2008. 11th Intersociety Conference on . 2008. 65. Lampin, E., Q.H. Nguyen, P.A. Francioso, and F. Cleri, T hermal boundary resistance at silicon silica interfaces by molecular dynamics simulations. Applied Physics Letters, 2012. 100 (13): p. 131906. 66. Munetoh, S., T. Motooka, K. Moriguchi, and A. Shintani, Interatomic potential for Si O systems using Tersoff p arameterization. Computational Materials Science, 2007. 39 (2): p. 334 339. 67. Chen, J., G. Zhang, and B.W. Li, Thermal contact resistance across nanoscale silicon dioxide and silicon interface. Journal of Applied Physics, 2012. 112 (6): p. 064319 064319. 6 8. Yao, M., T. Watanabe, P.K. Schelling, P. Keblinski, et al., Phonon defect scattering in doped silicon by molecular dynamics simulation. Journal of Applied Physics, 2008. 104 (2). 69. Li, X.B. and R.G. Yang, Effect of lattice mismatch on phonon transmissi on and interface thermal conductance across dissimilar material interfaces. Physical Review B, 2012. 86 (5).

PAGE 116

116 70. Chen, G., T. Zeng, T. Borca Tasciuc, and D. Song, Phonon engineering in nanostructures for solid state energy conversion. Materials Science and Engineering a Structural Materials Properties Microstructure and Processing, 2000. 292 (2): p. 155 161. 71. Wang, J. and J.S. Wang, Single mode phonon transmission in symmetry broken carbon nanotubes: Role of phonon symmetries. Journal of Applied Physics, 2009. 105 (6): p. 063509. 72. Ge, Y.Z. and S.G. Chen, Effect of Anharmonicity on Kapitza Resistance. Solid State Communications, 1991. 77 (4): p. 313 315. 73. Hu, M., P. Keblinski, and P.K. Schelling, Kapitza conductance of silic on amorphous polyethylene interfaces by molecular dynamics simulations. Physical Review B, 2009. 79 (10): p. 7. 74. Mazumder, S. and A. Majumdar, Monte Carlo study of phonon transport in solid thin films including dispersion and polarization. Journal of Hea t Transfer Transactions of the Asme, 2001. 123 (4): p. 749 759. 75. McGaughey, A.J.H. and A. Jain, Nanostructure thermal conductivity prediction by Monte Carlo sampling of phonon free paths. Applied Physics Letters, 2012. 100 (6): p. 061911. 76. Cahill, D.G. , P.V. Braun, G. Chen, D.R. Clarke, et al., Nanoscale thermal transport . II . 2003 2012. 2014. 011305 . 77. Grimmer, H., W. Bollmann, and D.H. Warrington, Coincidence site lattices and complete pattern shift in cubic crystals. Acta Crystallographica Sect ion A, 1974. 30 (2): p. 197 207. 78. Nerikar, P.V., K. Rudman, T.G. Desai, D. Byler, et al., Grain boundaries in uranium dioxide: Scanning electron microscopy experiments and atomistic simulations. Journal of the American Ceramic Society, 2011. 94 (6): p. 18 93 1900. 79. Ju, S. H. and X. G. Liang, Investigation on interfacial thermal resistance and phonon scattering at twist boundary of silicon. Journal of Applied Physics, 2013. 113 (5): p. 053513 053513. 80. Watanabe, T., S.G. Srivilliputhur, P.K. Schelling, J .S. Tulenko, et al., Thermal Transport in Off Stoichiometric Uranium Dioxide by Atomic Level Simulation. Journal of the American Ceramic Society, 2009. 92 (4): p. 850 856. 81. Pang, J.W.L., A. Chernatynskiy, B.C. Larson, W.J.L. Buyers, et al., Phonon densit y of states and anharmonicity of UO 2 . Physical Review B, 2014. 89 (11): p. 115132.

PAGE 117

117 82. Hicks, J., M. Sprinkle, K. Shepperd, F. Wang, et al., Symmetry breaking in commensurate graphene rotational stacking: Comparison of theory and experiment. Physical Review B, 2011. 83 (20): p. 205403. 83. Chernatynskiy, A., Private communications . 2014. 84. Schelling, P.K., Kapitza conductance and phonon scattering at grain boundaries by simulation. Journal of Applied Physics, 2004. 95 (11): p. 6082 6082. 85. Becker, B., P.K. Schelling, and S.R. Phillpot, Interfacial phonon scattering in semiconductor nanowires by molecular dynamics simulation. Journal of Applied Physics, 2006. 99 (12): p. 123715 123715. 86. Roberts, N.A., Phonon Wave Packet Simulations of Ar / Kr Interfaces. 2 006. 87. Yao, M., T. Watanabe, P.K. Schelling, P. Keblinski, et al., Phonon defect scattering in doped silicon by molecular dynamics simulation. Journal of Applied Physics, 2008. 104 (2): p. 024905 024905. 88. Feuchtwang, T.E., Dynamics of a Semi Infinite C rystal Lattice in a Quasiharmonic Approximation. II. The Normal Mode Analysis of a Semi Infinite Lattice. Physical Review, 1967. 155 (3): p. 731 744. 89. Wang, J. and J. S. Wang, Mode dependent energy transmission across nanotube junctions calculated with a lattice dynamics approach. Physical Review B, 2006. 74 (5): p. 054303 054303. 90. Wang, J. and J. S. Wang, Characteristics of phonon transmission across epitaxial interfaces: a lattice dynamic study. Journal of Physics: Condensed Matter, 2007. 19 (23): p. 2 36211 236211. 91. Wang, J.S., J. Wang, and J.T. Lü, Quantum thermal transport in nanostructures. The European Physical Journal B, 2008. 62 (4): p. 381 404. 92. Jex, H., The transmission and reflection of acoustic and optic phonons from a solid solid interface treated in a linear chain model. Zeitschrift fr Physik B Condensed Matter, 1986. 63 (1): p. 91 95. 93. Paranjape, B.V., N. Arimitsu, and E.S. Krebes, Reflection and transmission of ultrasound from a planar interface. Journal of Applied Physics, 19 87. 61 (3): p. 888 888. 94. Young, D.A. and H.J. Maris, Lattice dynamical calculation of the Kapitza resistance between fcc lattices. Physical review B, 1989. 40 (6): p. 3685 3693. 95. Zhao, H. and J.B. Freund, Lattice dynamical calculation of phonon scatter ing at ideal Si Ge interfaces. Journal of Applied Physics, 2005. 97 (2): p. 024903 024903.

PAGE 118

118 96. Steinbrüchel, C., The scattering of phonons of arbitrary wavelength at a solid solid interface: Model calculation and applications. Zeitschrift für Physik B Conde nsed Matter, 1976. 299 : p. 293 299. 97. Lumpkin, M.E., W.M. Saslow, and W.M. Visscher, One dimensional Kapitza conductance: Comparison of the phonon mismatch theory with computer experiments. Phys. Rev., B;(United States), 1978. 17 : p. 4295 4302. 98. Zhang , L., P. Keblinski, J. S. Wang, and B. Li, Interfacial thermal transport in atomic junctions. Physical Review B, 2011. 83 (6): p. 064303 064303. 99. Han, H., L.G. Potyomina, A.a. Darinskii, S. Volz, et al., Phonon interference and thermal conductance reduct ion in atomic scale metamaterials. Physical Review B, 2014. 89 (18): p. 180301 180301. 100. Weisstein, E.W. Least Squares Fitting . Available from: .

PAGE 119

119 BIOGRAPHICAL SKETCH Bowen Deng was born and raised degree in polymer materials science and engineering from Beijing University of Chemical Technology in Beijing, China in 2009. In the August the same year, he began his graduate study at the Department of Materials Science and Engineering at University of Florida under the advising of Dr. Simon Phillpot. He is expecting to get his Ph.D. degree in August , 2014.