UFDC Home  myUFDC Home  Help 



Full Text  
PAGE 1 1 MOLECULAR DYNAMICS STUDY OF THERMAL CONDUCTIVITY OF BISMUTH TELLURIDE By SHIKAI WANG A THESIS PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE UNIVERSITY OF FLORIDA 2013 PAGE 2 2 2013 Shikai Wang PAGE 3 3 To my Mom, for her kindness. To my Dad, for his tolerance. To my sister, for her loveliness. PAGE 4 4 ACKNOWLEDGMENTS I would like to thank my parents and my little sister at this very beginning. I cannot enjoy my life without their both material and spiritual support all along. They show me one truth that no matter what trouble I have, families are always my solid suppor ter and biggest fans. With no doubt, I must than k Dr. Youping Chen who prov ides me this opportunity to experience the happiness and sorrow of doing research and who will also fully support me to travel to IMECE conference in November. I also thank Dr. Curt is Taylor for his interest in molecular dynamics simulation and willingness be my committee member and review my thesis work. I t hank my lab colleges: Liming Xiong, Ning Zhang, Shengfeng Yang, Xiang Chen, Rui C he, Chen Zhang, and Zexi Zheng for their techn ic supports and advices. Special thank s is given to my ex girlfriend for her sustained encouragement during the time we spen t together. And I truly hope she will find her Mr. Right someday. At last, I declare that I am an agnostic. But if there is God, for creating this wonderful world, I sincerely thank Him. PAGE 5 5 TABLE OF CONTENTS page ACKNOWLEDGMENTS ................................ ................................ ................................ .. 4 LIST OF TABLES ................................ ................................ ................................ ............ 6 LIST OF FIGURES ................................ ................................ ................................ .......... 7 LIST OF ABBREVIATIONS ................................ ................................ ............................. 8 ABSTRACT ................................ ................................ ................................ ..................... 9 CHAPTER 1 INTRODUCTION ................................ ................................ ................................ .... 11 Bismuth Telluride ................................ ................................ ................................ .... 11 Molecular Dynamics ................................ ................................ ................................ 12 Software ................................ ................................ ................................ .................. 13 2 MODELING WITH TWO POTENTIALS ................................ ................................ .. 16 Simulation Parameters ................................ ................................ ............................ 16 Orthogonal Simulation Box Model ................................ ................................ .......... 17 Two Avai lable Potentials ................................ ................................ ......................... 18 3 MD SIMULATION OF SINGLE CRYSTAL MODEL ................................ ................ 23 Relaxation ................................ ................................ ................................ ............... 23 Direct Method ................................ ................................ ................................ ......... 23 Background ................................ ................................ ................................ ...... 23 Simulation of Nano Scale Specimen ................................ ................................ 25 Finite Size Effect ................................ ................................ .............................. 26 Result ................................ ................................ ................................ ...................... 28 4 MD SIMULATION OF BI CRYSTAL MODEL ................................ .......................... 32 Simulation Procedure ................................ ................................ .............................. 32 Result ................................ ................................ ................................ ...................... 3 2 5 CONCLUSION AND DISCU SSION ................................ ................................ ........ 36 LIST OF REFERENCES ................................ ................................ ............................... 38 BIOGRAPHICAL SKETCH ................................ ................................ ............................ 39 PAGE 6 6 LIST OF TABLES Table page 2 1 ................................ ................. 22 3 1 Table of thermal conductivities of single crystal bismuth. ................................ ... 31 PAGE 7 7 LIST OF FIGURES Figure page 1 1 Figures show the relationship between a hexagonal lattice and a rhombohedral lattice of Bi 2 Te 3 ................................ ................................ ........... 14 1 2 Figures illustrates the relationship between Bi 2 Te 3 self made orthorhombic Lattice. ................................ ................................ ......... 15 2 1 Diagrams above give the temperature di stributions under different heat currents applied. ................................ ................................ ................................ 20 2 2 Figures show an 18nm atomic model applied in the simulation. ......................... 21 2 3 The color atoms are a part of the bi crystal model. ................................ ............. 21 2 4 Figures show the brief procedure to search angles. ................................ ........... 22 3 1 Diagram shows the stabilizing temperature. ................................ ....................... 29 3 2 This diagram gives that results of adding temperature and adding heat flux are very close. ................................ ................................ ................................ .... 29 3 3 A sketch of direct method under periodic boundary condition. ........................... 30 3 4 This diagram shows the comparison of two averages in a 24nm model. ............ 30 3 5 Diagram gives the extrapolation of bulk thermal conductivity. ............................ 31 4 1 Figure above presents the atom arrangement of a PBC bi crystal model 8 after relaxation. ................................ ................ 34 4 2 Similar to Figure 3 3, this figure illustrates the direct method model with grain boundary, where the black be lts are the location of twin boundaries. ................. 34 4 3 From the comparison of temperature distributions of two 24nm models, the temperature gradients are different near the twin boundary. .............................. 34 4 4 Similar to the 24nm model comparison in Figure 4 3, slightly d ifferences are generated by the twin boundary. ................................ ................................ ........ 35 PAGE 8 8 LIST OF ABBREVIATIONS EMD Equilibrium Molecular Dynamics MD Molecular Dynamics NEMD Non equilibrium Molecular Dynamics PBC Periodic Boundary Condition PMFP Phonon Mean Free Path PAGE 9 9 Abstract of Thesis Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Master of Science MOLECULAR DYNAMICS STUDY OF THERMAL CONDUCTIVITY OF BISMUTH TELLURIDE By Shikai Wang May 2013 Chair: Youping Chen Major: Mechanical Engineering Bismuth telluride (Bi 2 Te 3 ) is one of the best thermoelectric materials at room temperature. Since the crystal structure of Bi 2 Te 3 is very complicated and hard to be modeled, a new method based on a self made orthorhombic lattice is applied to build atomic specimen in orthogonal bo xes easily. Generally in a molecular dynamics (MD) simulation, the thermal conductivity of Bi 2 Te 3 is measured by using Green Kubo method which is an equilibrium method and able to calculate the whole thermal conductivity tensor in just one simulation. How ever, more detailed information could be obtained by applying direct method instead which is a non equilibrium MD method, such as temperature distribution in a specimen or even near a grain boundary in a heat transfer which could give us further understand ings of the thermal properties of Bi 2 Te 3 Due to the large temperature gradient caused by the low thermal conductivity, the direct method based MD simulation of Bi 2 Te 3 could be challenging. In this paper, with the interatomic potentials studied by Huang et al 2 and Qiu et al 8 appropriate heat fluxes are applied to the tests. And an analysis of finite size effect shows the result which agrees with both experiment data and the result of Green Kubo method quite well. PAGE 10 10 Finally the temperature distribution and Kapitza conductance of bi crystal specimens with twin boundaries is studied which gives the conclusion that twin boundaries may not be able to generate very much heat resistance in Bi 2 Te 3 PAGE 11 11 CHAPTER 1 INTRODUCTION Bismuth Telluride Nowadays, the unstoppable increase of energy consumption has be en becoming a severe problem as fossil fuel, one of the main energy source, is less and less. For decades, people are seeking more alternative energy resources an d thermoelectric material could be a good choice to collect energy which has been widely applied to solid state cooling devices to convert electricity to heat and thermoelectric devices inversely to convert heat into electrical energy. Usually a good therm oelectric material should have not only high electrical conductivity but also low thermal conductivity. Bismuth Telluride ( Bi 2 Te 3 ), considered as one of the best thermoelectric materials has a very low thermal conductivity A nd its alloys give high thermo electric figure of merit, defined as where is the Seebeck coefficient, is the electrical conductivity, is the thermal conductivity, and is the absolute temperature 12 A material that can afford competitive performance must have a value higher than 1 and a high merit of 1.4 at 100 was reported in a p type nanoc rystalline BiSbTe bulk alloy 1 Besides, at near room temperature (0 to 250 ), the best commercial thermoelectric materials for applications are still the Bi 2 Te 3 b ased materials 12 As stated in other publications 2,3 the thermal conductivity of Bi 2 Te 3 is around 1~2W/ mK based on the results of experiments and simulations. The lattice of Bi 2 Te 3 is mostly visualized as a layer structure 5 ( Figure 1 1) which has a rhombohedral crystal structure and includes five atoms along the trigona l axis in the sequence of Bi Te (1) Te (2) Te (1) Bi. The rhombohedral unit cell parameters are and 2 Conventionally, Bi 2 Te 3 also can be indexed in terms PAGE 12 12 of a hexagonal unit cell with lattice parameters and 4 which contains three blocks and each block is a five layer packet 2 with the atom sequence of Te (1) Bi Te (2) Bi Te (1) Chemical bonds vary between layers. Te (2) and Bi atoms in adjacent layers are connected by covalent bonds. And Te (1) and Bi atoms in adjacent layers are linked by both covalent and ionic bonds. These two bonds are quite strong comparing to the bonding between Te (1) and Te (1) atoms in adjacent lay ers, which is just bonded by Van der Waals force 6 For modeling convenience, a self made unit ( Figure 1 2 ) with orthogonal repeatable boundaries is applied. This base centered orthorhombic crystal lattice also contains three blocks and each block has a fiv e layer atom structure. Totally there are 30 atoms in this lattice structure and 2 atoms in each layer. Molecular Dynamics Regarding an atom as the basic rigid body, m olecular dynamics (MD) is an approach that calculating the movements of a specific group of atoms to study material behavior by computer simulations Generally, only molecular models and interatomic potentials are needed to achieve an MD simulation After years of development, especially the enhancing performance of computers, MD simulation ha s been widely used in many areas, including material science. This paper mainly focus on the calculation of thermal conductivity by MD simulations. Generally, t he re are two major approaches to compute thermal conductivity : one is nonequilibrium MD (NEMD) method and the other one is equilibrium MD (EMD) method such as the Green Kubo method respectively 7 which are the two most popular methods for computing thermal conductivity 7 The direct method is applied by generate a PAGE 13 13 temperature gradient across the model which is very similar to a n experimental case 7 On the contrary, the Green Kubo method calculate s thermal conductivity through the heat fluctuation inside a simulation specimen which does not have a temperature gradient 7 In this paper all the conductivity computations are via the direct method since the Green Kubo method has been tested on Bi 2 Te 3 w hich gave rather good results 2,8 and the direct method application is still rare. The large temperature gradient as a necessity in the direct method may cause problems However, as one benefit of this method by imposing temperature gradient, detailed information could be demonstrated directly, such as temperature distribution along the heat flux direction and near grain boundaries. Based on the potentials developed by Huang et al 2 and Qiu et al 8 the simulations will be applied to a group of model s of different lengths and then calculate the thermal conductivity of infinite length (bulk thermal conductivity) by finite size effect analysis. Software Lammps and VMD are two useful softwares when people deal with MD simulations. Lammps is a powerful MD simulator which is distributed by Sandia National Laboratories as an open source code. Generally, Lammps is able to solve the problem w hen the potential, boundary conditions and loads of a n atomic model are set A lthough large model which contains millions of atoms could be time consuming to reach convergence As one essential advantage, Lammps is able to run eithe r on single processors or in parallel which makes massive simulation s be available. PAGE 14 14 VMD (Visual Molecular Dynamics) is a visualization program developed by the Theoretical and Computational Biophysics Group in University of Illinois at Urbana Champaign. The purpose of using VMD is that this software is able to read output files i.e. xyz files, for displaying, animating, and analyzing. In this paper, most of the atomic schematics are output from VMD. Figure 1 1 Figures show the relationship between a hexagonal lattice and a rhombohedral lattice of Bi 2 Te 3 A) A conventional crystal structure in a hexagonal lattice which contains three blocks and each block has five layers. Lattice con stants are marked in the figure. B) Hexagonal lattice is shown with a coordinate system and a corresponding rhombohedral lattice. C) The rhombohedral lattice of Bi 2 Te 3 is shown with lattice constants. PAGE 15 15 Figure 1 2. Figures illustrates the relationship between Bi 2 Te 3 self made orthorhombic Lattice. A) A hexagonal lattice of bismuth telluride. B) A coordinate system and self made lattice are shown in a hexagonal lattice. C) The structure of the self made lattice which has a structure of b ase centered o rthorhombic crystal. PAGE 16 16 CHAPTER 2 MODELING WITH TWO POTENTIALS Simulation Parameters While modeling specimens, the first step is to determine the geometry parameters, such as length and cross section. Different lengths and cross sections may lead to different results which will be discussed The simulations in this paper are done with several specimens of different sizes. The thermal conductivity is calcula ted by the slope fit ted from the temperature distribution data collected from slices along the length directio n of the specimen And each slice gi ve s an average temperature of all the atoms within it. Although as stated in Ref. 7, the cross section did not show strong effect on simulation result. But in my simulations, if a too small cross section is applied, temperature gradient could be obviously affected. Specimen with small cross section could not provide enough atoms in each sampling slices for da ta co llection. On the contrary big one could consume much more computation resources. As tested for many times, the cross section of 5.98nm 2 should be a good choice. Except cross section, model length will influent the result as well. Short specimen gives les s sampling slices which may not have an accurate result but makes large temperature gradient available On the other hand, long specimen gives enough sampling slices but temperature gradient would remain very small since small temperature gradient represe nts small heat flux which could be easily affected by system error. Besides, in order to employ finite size effect, the total length should not be too much bigger than the phonon mean free path of Bi 2 Te 3 which is only around 2 nm 9 PAGE 17 17 at 300K Under these considerations, a range of 18nm to 60nm was appropriately selected. Prior work 7 has indica ted that loading different heat flux could also affect the MD simulation result. Low heat flux could be affected by comparable system noise and lose lot s of accuracy. However, high heat flux model could also be inaccurate that generate negative energy value which results in simulation fail. Shown in Figure 2 1 larger hea t flux result s in smoother curve. But in the last diagram, the temperature in heat source is beyond the potential limit 500K. The reason that heat current instead of heat flux is used in the diagrams is that all the models have the same cross section and t hus heat current comparison is identical to heat flux comparison. Bi 2 Te 3 based materials dominate at relatively near room temperature (0 to 250 ) 1 and in addition, the Debye temperature of Bi 2 Te 3 is 162K 10 The simulation temperature should be higher than the Debye temperature and thus a temperature of 400K (127 ) was chosen during all the simulations. Orthogonal Simulation Box Model In the prior works about Bi 2 Te 3 M D simulation, orthogonal specimens were hardly been reported mainly because of the difficulties in Bi 2 Te 3 modeling which brought by its complicated crystal structure. In many cases, an orthogonal simulation box could be much simpler to model and control. Based on the self made b ase centered orth orhombic crystal s tructure ( Figure 1 2 ) specimens of different lengths along z direction from 18nm to 60nm and the same cross section of 5.98 nm 2 were built with repeatable orthogonal boundaries. ( Figure 2 2 ) PAGE 18 18 Besides, bi crystal models with twin boundaries were built with the orthogonal unit cell as well. The only difference between the bi crystal and single crystal models was section remain s the same and lengths along z direction were selected to be 24nm and 36nm. In Figure 2 3 the overlapping of scanned image and atomic model is quite precise and hence the bi crystal model is ready to use Potentials Since Bi 2 Te 3 is a relatively new material for MD study, there are only two available potentials so far which are from the previous works of Huang et al 2 and Qiu et al 8 Particularly the latter is based on Morse potential and Coulomb force together and the former potential combin es Morse potential, cosine square angular potential, and Coulomb force. The angular potential also makes the potential developed by Huang et al 2 more difficult to be employed to atomic models, especially in simulation boxes of periodic boundary condition. It is important that the cutoff of each pair was chosen to be a little larger than the specific atom distance between and the cutoff radius of the electrostatic terms is ad justed from 12 to 10 as recommended by Huang in our private contact. Table 2 1 gives the detail information of the cutoffs of each pair listed in his potential. almost r emain stable. To deal with the adjacent/same layers pair, the atoms in each block in a sequence of Te (1) Bi Te (2) Bi Te (1) were renamed as Te (1 )( 1) Bi (1) Te (2) Bi (2) Te (1 )( 2) and hence there were five types of atoms Therefore, to distinguish adjacent / same layers pair became very easy. PAGE 19 19 When applying angular potentials, some more effort are neede d. In Lammps, to make the angles readable to the program, all the angles must be listed by the three indexes of the atoms in that specific angle. In order to lis t these indexed out, the method that I applied is to find the center atom one by one, then search the nearest atoms that belong to the angles for each center atom and finally all the angles could be located in such approach The larger the model is, the s ignificant longer searching time will be needed. (Fig ure 2 4 ) In most cases, pe riodic boundary condition (PBC) should be employed in an MD thermal simulation. To apply angular potential in a PBC case, First I wrap the whole specimen with atoms which are l ocated in the nei ghbor boxes Then mark these atoms with the index information of their mapping atoms in the current simulation box. At last, use the method stated in last paragraph and hence all the angles would be listed in a data file which could be rea dy to associate with the original model which does not have wrapping atoms. The total angle amount in a PBC mode l could be expressed as where is the amount of angles and is the amount of total atoms. PAGE 20 20 Figure 2 1 Diagrams above give the temperature distributions under different heat currents applied. All of these models share the same length and cross section. In each diagram, the highest partial indicates the location of heat source and the lowest one indicates the location of heat sink. 380 390 400 410 420 Temperature (K) A 0.11ev/ps 350 375 400 425 450 Temperature (K) B 0.21ev/ps 300 350 400 450 500 Temperature (K) C 0.53ev/ps 200 300 400 500 600 0 4 8 12 16 20 24 Temperature (K) Length (nm) D 1.06ev/ps PAGE 21 21 Figure 2 2 Figures show a n 18nm atomic model applied in the simulation. A) X direction view of an 18nm single crystal model. B) Y direction view of the same model. Figure 2 3 The color atoms are a part of the bi crystal model. The blur background picture is HAADF STEM image of the basal twin projected along direction 4 A B PAGE 22 22 Figure 2 4 Figures show the brief procedure to search angles. A) The green atom is the center atom, the bonded red atoms are to be found B) One possible angle type (adjacent layer angles). C) The other possible angle type (three adjacent layer angles). Table 2 1. Table of a dditional cutoff s for potential 2 Interaction Cutoff ( ) Te (1) Bi (adjacent layers) 4.5 Te (2) Bi (adjacent layers) 4.5 Te (1) Te (1) (adjacent layers) 4.5 Bi Bi (same layer) 5.5 Coulomb force (All atoms) 10 A B C PAGE 23 23 CHAPTER 3 MD SIMULATION OF SINGLE CRYSTAL MODEL Relaxation The first step of an MD simulation after modeling is to relax the system which is a processing that minimize s the total energy of the whole system and let the system get to an equilibrium condition. The meth od of applying relaxation process is quite simple which is to run the simulation without any load for a period such as 1ns with NVE ensemble Generally, the model could not perfectly consist with the potential A s shown in Figure 3 1, the system temperat ure fluctuate very severely at beginning and become sta b le after a while which gives a typical relaxation process. Usually the time for relaxation varies in different situations. In this paper, the relaxation time for all simulation cases is 0.25ns and the time step is set to 0.25fs, hence totally one million time steps will be undertaken during the relaxation Then the temperature of the specimen will be increasing to 400K gradually with NVT ensemble After the system temperature reach 400K, keep the simulation running with NVE ensemble for another 0.25ns to eliminate probable unstable factors as many as possible. Direct Method Background As introduced in chapter 1, direct method is an NEMD method to calculate thermal conductivity and it is very similar to an experimental measurement While employing direct method either temperature or heat flux could be loaded on the system to generate an appropriate temperature gradient. When the two systems both reach stable status, t he same re sult should the simulations give based on these two loading PAGE 24 24 approaches ( Figure 3 2). In my thesis work, all simulations are employing the energy loading approach, as it is easier to compute the heat flux in the specimen and the reference temperature could be easily determined as 400K which is initially fixed As shown in Figure 3 3, same amount of energy is added to the heat source and subtracted from the heat sink every time step and hence the whole system remains energy conserved. is the heat flux which is measured in [W/m 2 ] and defined as: ( 3 1) Where which is measured in [ev] is the cross section of the specimen, is the time step used in the simulation, and hence is the added/subtracted energy eve ry time step as a convenience to edit input file The coefficient in this equation means the unit conversion from [ev] to [J]. T he arrows in F igure 3 3 show the direction of the heat flux and t ( 3 2) Where is thermal conductivity and is temperature gradient. The negative sign in this formula means the directions of heat flux and temperature gradient are opposite. Thus, the thermal conductivity could be expressed in this way: ( 3 3) Where the absolute value of temperature gradient is used to ensure the thermal conductivity remain a positive value It is notable that while a large temperature gradient applied to the system, strong nonlinear behavior will appear and which could make this method invalid. Even with PAGE 25 25 relative small temperature gradient, there still will be nonlinear behaviors at boundaries of source and sink caused by boundary scattering 7 In some materials, i.e. silicon, this boundary non linearity is very significant. Simulation of Nano Scale Specimen As mentioned in the first chapter, a t near room temperature (0 to 250 ), it is hard to find a competitor to Bi 2 Te 3 b ased material for its outstanding thermoelectric property. Thus, the reference temperature of 400K was selected which is also higher than the Debye temperature 162K 10 During the simulation, the kinetic energy of all atoms are collected slice by slice and the thickness of each slice is 5nm which contains 90 atoms in total T he temperature of each slice could be computed with the relationship between kinetic energy and temperature: ( 3 4) Where is the kinetic energy of all the atoms in each slice constant, and is the needed temperature. In thermal simulation with direct method, long time average is usually needed to reduce statistically fluctuation 11 after stable status achieved First the relaxed system runs for 2 million time step s (0.5ns) to reach a stable condition, then two data collections are applied during the following 4 million time step s, each data collection runs for 2 million time step s respectively. In F igure 3 4 the orange dots and the blue dots show the temperature averaged from the two collections. The well overlapped pattern indicate that the system has been reach a stable status after the 0.5ns running before the collection process starts PAGE 26 26 The temperature gradien t needs to be calculated before the th ermal conductivity computation. Very obvious nonlinear curves present near the heat source and heat sink boundaries which are the result of strong phonon boundary scattering near that regions ( Figure 3 2 and 3 4) but respectively at intermediate part of the specimen and the part near the periodic boundary along the length direction. The temperature gradient is given by the two fitted slopes from these two regions. Apply with the formulas which introduced in the last section, the thermal conductivities of all the single crystal models are calculated and listed in Table 3 1. In the next section, finite size effect will be introduced in detail, and which could give an approach to approximate the thermal conductivity of infinite specimen which could be regarded as the Bi 2 Te 3 Finite Size Effect T he thermal conductivity could be affected by the size of the simulation box along the heat flux direction when the specimen length is not significantly longer than the materia l s phonon mean free path this phenomenon is named finite size effect also know n as the Casimir limit 12 Generally smaller specimen dimension gives lower thermal conductivity In ideal single crystal, t h e limited effective thermal conductivity could be expressed in such way 13 : ( 3 5) Where is effective thermal conductivity and its inverse is the effective thermal resistance per unit length, is the material bulk thermal resistance per unit length which PAGE 27 27 is a material dependent property and is the boundary thermal resistance, including both heat source and sink boundary resistance, and is the simulation box length. Thus, this expression indicate s that the total e ffective thermal resistance equals to the sum of material thermal res istance and boundary thermal resistance and that obviously the inverse of the effective thermal conductivity has a linear relation ship to the inverse of the specimen When the length of the specimen gets infinity, could clearly be identical with In further detail, k inetic theory 7 give an equation to determine thermal conductivity: ( 3 6) Where is the material specific heat, is the sound velocity in the material, and is the phonon mean free path. Insert this expression into Equation 3 5 gives: ( 3 7) Where is the effective phonon mean free path, is the bulk phonon mean free path. In Ref. 7 it has been indicated that And the specific heat could be expressed as: ( 3 8) Where is the number density of atoms in the system and in Bi 2 Te 3 nm 3 Thus, substituting Equation 3 6 and Equation 3 8 into Equation 3 5 gives: PAGE 28 28 ( 3 9) Result Thermal conductivity computed from different length specimens are listed in Table 3 1. Obviously, longer specimen gives larger thermal conductivity as finite size effect said 2 gives larger results. Based on the finite size effect, the bulk thermal conductivity could be extrapolated from the finite si ze results of different lengths In Figure 3 5, when the inverse of is zero, that is, gets to infinity, the inverse of the value on the vertical axis (the intercept) would be the bulk thermal conductivity which is 1.566 W/ mK for Huang s potential 2 a nd 0.609W/mk for Qiu s potential 8 Besides, t he result of Green Kubo method is around 1W/ mK and result of experiment is about 2W/ mK 2,3 Compared with these results, this value is quite close to those In other words, the accurac ies of the two potential s are acceptable In the next Chapter, a set of bi crystal specimen is employed with the same simulation procedure and the influence of twin boundary to the temperature distribution is discussed. PAGE 29 29 Figure 3 1 Diagram shows the stabilizing temperature. After 1000,000 steps, the system temperature has been nearly constant around 27K. Figure 3 2 This diagram gives that results of adding temperature and adding heat flux is 5.7K which could be mo dified by adjusting the load of temperature control method. 0 10 20 30 40 50 60 0 200000 400000 600000 800000 1000000 Temperature (K) Time Step 255.7 305.7 355.7 405.7 455.7 505.7 555.7 250 300 350 400 450 500 550 0 4 8 12 16 20 24 Temperature of heat flux control model (K) Temperature of temperature control model (K) Length (nm) Temperature control Heat flux control PAGE 30 30 Figure 3 3 A sketch of direct method under periodic boundary condition. Where is heat flux energy of is added and subtracted from the heat source and sink respectively. Thus a temperature gradient could be established. Figure 3 4 This diagram shows the comparison of two averages in a 24nm model. Each average is computed through 2 million time st eps. 250 300 350 400 450 500 550 0 4 8 12 16 20 24 Temperature (K) Length (nm) 1st average 2nd average PAGE 31 31 Figure 3 5 Diagram gives the extrapolation of bulk thermal conductivity. The blue dots are the MD result of finite size thermal conductivity and the red dot is the extrapolated point. Note that the horizontal and vertical axes are labeled in the inverse of length and thermal conductivity. These data are based on the potential developed by Huang et al 2 Table 3 1. Table of thermal conductivit ies of single crystal bismuth Length (nm) Thermal conductivity (W/mK) 2 18 0.550 24 0.672 30 0.800 36 0.789 60 1.003 + 1.566 8 18 0.206 24 0.244 30 0.273 36 0.314 + 0.609 0 0.2 0.4 0.6 0.8 1 1.2 0 0.01 0.02 0.03 0.04 0.05 0.06 1/thermal conductivity (mK/W) 1/length (nm 1 ) PAGE 32 32 CHAPTER 4 MD SIMULATION OF BI CRYSTAL MODEL Simulation Procedure The relaxation process applied on bi crystal model is the same as the one did on single crystal model. In Figure 4 1, atoms near twin boundary become messy after relaxation in the model associated with Q 8 However, this twin boundary structure based on 2 keeps a stable status during the simulation. Hence, the bi crystal 2 only. Still, the direct meth od is employed here (Figure 4 2 ) The goal of this simulation is to observe some possible effects on temperature distribution generated by the twin boundary. Theoretically, phonon will be scattered when going through a grain boundary Locally at this inter face, there will be a n additional thermal resistance and hence thermal conductance as its inverse, named Kapitza conductance which have this following expression 14 : (4 1 ) Where is the heat flux, ( not ) is the temperature discontinuity at the interface, and is the Kapitza conductance whose inverse is Kapitza resistance Result The measurement applied on bi crystal model is the same as the one did on single crystal model except for the finite size effect analysis With 120 and 180 slices along the length direction, the results of collected temperature are shown in Figure 4 3 and 4 4 PAGE 33 33 Very slight difference is found between the curves of single and bi crystal models. However, the discontinuous slope of temperature could be observed. The intercepts at twin boundary of fitted line s are calculated which are derived from the temperature data on both left and right side of grain boundary. The difference of the two intercepts is regarded as the temperature discontinuity at the interface. The heat flux inside the specimen is 14.3GW/m 2 in both two twin boundary models. Therefore, the Kapitza conduct ance could be computed from Equation 4 1 which gives 1.94GW/m 2 K and 1.19GW/m 2 K respectively in 24nm and 36nm models. Recall the Equation 3 5, similar equation could be achieved for twin boundary crystal : (4 2 ) Where is the resultant thermal conductivity of specimen with twin boundaries, the summation of first and second terms on the right hand side is identical to the inverse of in the Equation 3 5, and the numerator 2 in the third term accounts for the fact that there are two grain boundaries in the specimen of length Therefore, the idea of this equation is that the overall thermal resistance is equal to the summation o f material thermal resistance, finite size boundary resistance and grain boundary thermal resistance. PAGE 34 34 Figure 4 1 Figure above presents the atom arrangement of a PBC bi crystal model 8 after relaxation. The twin bounda ry (in the middle of the figure) structure is not stable. Figure 4 2 Similar to Figure 3 3, this figure illustrates the direct method model with grain boundary, where the black belts are the location of twin boundaries. Figure 4 3 From the comparison of temperature distributions of two 24nm models, the temperature gradients are different near the twin boundary. The boundary also contributes to higher heat source temperature and lower heat sink temperature. Thus generally, the twin bounda ry do affects the thermal transport. 250 300 350 400 450 500 550 0 6 12 18 24 Temperature (K) Length (nm) bicrystal singlecrystal PAGE 35 35 Figure 4 4 Similar to the 24nm model comparison in Figure 4 3, slightly differences are generated by the twin boundary. 250 300 350 400 450 500 550 0 9 18 27 36 Temperature (K) Length (nm) bicrystal singlecrystal PAGE 36 36 CHAPTER 5 CONCLUSION AND DISCUSSION So far, potential from Ref. 2 and Ref. 8 have been tested in direct method simulations. Compare with the Bi 2 Te 3 public ations, i t is no doubt that both of the potentia l s give reasonable result s In Ref. 2, it is said the cutoff of Coulomb force is 12 angstrom. But actually this cutoff could result in a disorder model that all the blocks are dislocated along each other. By private contact, the cutoff was suggested as 10 angstrom and which stabilize the model. The co mplexity of this potential 2 may bring issue that the running time is longer than the potential in Ref. 8 The compute speed could be increased a little by pre building the neighbor l ist, which was suggested by Huang as well. With these modification, this potential could be promising In this paper, a low cross plane thermal conductivity of bismuth telluride has been achieved. But some error could still exist which may due to the small amount of sampling slices or relatively low heat flux. Also the error may come from the potential itself. Generally, comparing with the results of experiment and Green Kubo simulation, the result in this paper is acceptable. In Ref. 16, a relatively high thermal conductivity of single crystal model associated with potential in Ref. 8 was given. It is i nteresting to notice a very large cutoff of pair interactions is chosen as 30 and they said this value is from the Ref. 8 of this paper. However I cannot find t he words about this large cutoff in Ref. 8. Their result seems questionable. This MD simulatio n gives the result that tw in boundary may not be able to strongly influence the heat conduction of bismuth telluride. It is interesting to see that PAGE 37 37 the Kapitza conductance value is very close to the one of some silicon grain boundary stated in Ref. 14 wher e very significant temperature drops appear at the grain boundaries. In addition, from Equation 4 2 we have ( 5 1 ) Where is the single crystal thermal conductivity which is numerically equals to the in Equation 4 5. From this equation, it is understandable that thermal conductivity will undoubtedly de crease when grain boundary exist since the ratio on the left hand side will always larger than one. And how much thermal effect could grain boundaries have depends on the ration of material thermal conductivity and boundary thermal resistance which is exactly the second term on the right hand side. In silicon, the grain boundaries mentioned in Ref. 14 could have a very large ratio In bismuth telluride, however, this ratio remains much smaller than one and hence there is An alternate expression of this equation has been stated in Ref. 15 When applying direct method, heat flux could be approximately calculated b efore ec ted temperature gradient and experimental thermal conductivity. However if the potential is not good enough, this method could bring trouble such as an unexpected temperature gradient. Also, the strong nonlinearity regions at heat source and sink could bring difficulties to finding proper temperature gradient. But pre determined heat flux could always be a choice. PAGE 38 38 LIST OF REFERENCES 1. B. Poudel, Q. Hao, Y. Ma, Y. C. Lan, A. Minnich, B. Yu, X. Yan D. Z. Wang, A. Muto, D. Vashaee, X. Y. Chen, J. M. Liu, M. S. Dresselhaus, G. Chen, and Z. F. Ren, Science 320, 634 (2008) 2. B. L. Huang and M. Ka viany, Phys. Rev. B 77, 125209 (2008) 3. P.A. Sharma, A.L.L. Sharma, D.L. Medlin, A.M. Morales, N. Yang, M. Barne y, J. He, F. Drymiotis, J. Turner, and T.M. Tritt, Phys. Rev. B 83, 235209 (2011). 4. D. L. Medlin, Q. M. Ramasse, C. D. Spataru, and N. Y. C. Yang, J Appl. Phys. 108, 043517 (2010). 5. J. Drabble and C. Goodman, J. Phys. Chem. Solids 5, 142 (1958). 6. Y. Tong, F. J. Yi, L. S. Lui, P. C. Zhai, and Q. J. Zhang, Comput. Mater. Sci. 48, 343 (2010). 7. P.K. Schelling, S.R. Phillpot, P. Keblinski, Phys. Rev. B 65, 144306 (2002). 8. B. Qiu and X. Ruan, Phys. Rev. B 80, 165203 (2009). 9. N. Peranio, O. Eibl, and J. Nurnus, J. Appl. Phys. 100, 114306 (2006). 10. G. E. Shoemake, J. A. Rayne, and R. W. Ure, Phys. Rev. 185, 1046 (1969). 11. G. Chen, Nanoscale Energy Transport and Conversion: A Parallel Treatment of Electrons, Molecules, Phonons, and Photons (Oxford University Press, New York, N ew York, 2005). 12. H. B. G. Casimir, Physica 5, 495 (1938). 13. J. Michalski, Phys. Rev. B 45, 7054 (1992). 14. P. K. Schelling, S. R. Phillpot, and P. Keblinski, J. Appl. Phys. 95, 6082 (2004). 15. S. Aubry, C. J. Kimmer, A. Skye, and P. K. Schelling, Phys. Rev. B 78, 0 64112 (2008). 16. K Termentzidis, O Pokropyvnyy M Woda S Xiong Y Chumakov P Cortona and S Volz J. Appl. Phys. 113, 013506 (2013) PAGE 39 39 BIOGRAPHICAL SKETCH Shikai Wang was born in Weinan, Shaanxi, China in 1989. From 2 007 to 2011, he studied in Hua zhong University of Science and Technology Wuhan, China. And in 2001, He graduated with a degree of Bachelor of Science in Energy and Power Engineering. From 2011 to 2013, Shikai was enrolled in University of Florida, FL, US to pursue a master degree in the Department of Mechanical and Aerospace Engineering. Since 2012, Shikai has lab doing research on molecular dynamics simulation. 