Molecular Dynamics in Nanoscale Heat Transfer

Get Complete Project Material File(s) Now! »

Molecular Dynamics in Nanoscale Heat Transfer

MD simulation is a real space analysis, dealing with atomic/molecular structures of materials, currently considered as a very powerful tool for investigating the details of nanoscale phenomena. Despite MD being a generic tool for understanding numerous physical aspects of matter, our discussion in this thesis is particularly limited to the heat conduction process. Depending on the strategy to evaluate interatomic forces, MD methods can be differentiated as ab-initio and classical approaches where quantum theory-based sub-atomic particle level evaluations and experimental data-based empirical interatomic potentials are used, respectively [68]. The use of the ab-initio MD method is restricted due to the higher computational demands, therefore, currently it is limited for systems with few atoms [69].  Nevertheless, classical MD simulations can be used for wide range of applications and it is very convenient and practical with respect to the current computational power. However, use of classical MD is limited in cases where quantum effects are critical such as in low temperature (relative to material Debye temperature) [70].Therefore, some works [70]–[72] introduced quantum corrections for classical MD simulations resulting if the system temperature is not near or above the Debye temperature. The fundamental inputs needed for MD simulations (from this point onward MD simulation refers to its classical form) are the atomic structure of the matter and interatomic potentials which satisfactorily characterize the potential fields among atoms. Newton‟s second law is applied for translational motion of each and every atom in the structure as given below. where m_i, r_i,f_i and represent the mass, position vector, force vector acting on i^th atom and resultant potential field of the whole system, respectively.  If the system consisting of P atoms and Cartesian components in 3D space is considered, 3P equations have to be solved to get the particle‟s position and momentum trajectories. Several algorithms are introduced for integrating Newton‟s equation to improve the computational efficiency while minimising numerical errors. Computationally efficient algoritFhms suggest computationally faster and lower in memory storage. From the physics point of view, selected algorithms must satisfy the time reversible property and thereby energy conservation which Newton‟s equation possesses. Even though literally none of the available algorithms satisfy the above requirements, it is widely accepted Verlet algorithms are a very good candidate. The Velocity Verlet algorithm is well-known for its easy implementation among other Verlet algorithms. The evaluation of proper interatomic potential function is the most challenging task of MD implementations since its success substantially affects the accuracy of results. The accurate prediction of potential functions is currently a combined effort of real experiments and density functional analysis [74]. It is a rather acceptable practice to decompose the combined potential energy in a P-atom system as a summation of one body, two body, three body and so on [75]. Thus, the total potential energy  is in the following form. Here one body potential represents the effects from walls and external electrical fields on the system. To be useful, it is a necessary condition that converges to zero immediately P increases. Lennard-Jones (LJ) potential [76] is the most widely used pair potential function. Among the different variants of LJ potential, LJ 12-6 form is well recognized as a suitable potential to model the properties of noble elements such as Argon [77]. In equation 2-3,  and σ have the dimensions of energy and length, respectively. The repulsive forces due to the overlapping of inner shell electrons and ions are modelled by the first term on the right hand side of equation 2-3 while second term models the attractive forces due to electrostatic effects. A large number of studies used this potential to demonstrate the heat conduction of matter [49], [61], [78]–[81]. Although, LJ potential is successful in modelling the behaviour of noble elements and some crystals, the stability of diamond structure of tetrahedral semiconductors such as Si and Ge cannot be described in a proper way. Thus, Stillinger-Weber (SW) potential [75] which comprises both two-body and three-body terms, was introduced for diamond lattices such as Si, C, Ge and GaAs. This potential was proved to be capable of replicating mechanical and thermal properties of Si. For example, the heat conduction of Si thin films [79] and nanowires [82], [83] was investigated using SW potential. Tersoff [84], [85] proposed a many-body potential called the Tersoff potential for materials such as Si, Ge. The Tersoff potential has been widely used to predict trends in the heat conduction process of materials, for example, crystalline/amorphous core/shell Si nanowires [86], Si/Ge superlattice nanowires [87], polycrystalline Si nanowires [88], carbon nanotubes [89], etc.

Direct Method for Evaluating Thermal Conductivity

The direct method is based on non-equilibrium molecular dynamics (NEMD) and is closely analogous with the experimental evaluation of TC. In this method, the structure in which TC must be calculated is sandwiched by heat source and sink similar to the experimental set up. Subsequently two fixed regions are attached to the two extremes of the simulation cell including the structure and two heat reservoirsin order to keep atoms attached in reservoirs.Basically there are two different implementation procedures in the direct method: in one precedure the heat flux is imposed while in the other the temperature gradient is imposed. In the former, a prescribed heat flux is added to the heat source while same amount is removed from the sink in the form of kinetic energy and the system is allowed to develop a steady state temperature gradient. In the latter method, the temperatures of the reservoirs are prescribed and the system is allowed to create a steady state heat flux. However, Lukes et al. [90] showed that the latter method takes a long simulation time for the system to approach a converged heat flux, therefore the former method is employed in this study.  Once the steady state temperature gradient is developed in the system, the TC (K) is evaluated using Fourier‟s law as follows, where q is the prescribed heat flux. Also it is important to note that the prescribed q must be big enough to create a large enough temperature difference between the heat reservoirs compared to the statistical temperature fluctuation of the system and small enough to create a linear temperature profile along the structure [91]. Even though the direct method uses finite sized structures to evaluate the TC, techniques were developed to infer the TC of infinitely large structures [92]. A large number of studies can be found in the literature which applied the direct method to evaluate the TC of structures including LJ materials [78], Si,diamond [92] and  Si, Si/Ge nanowires [87], [93] etc.

Conductivity

The Green-Kubo (GK) method is taken as a key approach to evaluate transport phenomena such as electric current, momentum, heat, etc. [94]. The underlying principle in this approach is to relate the currents due to an externally applied field in a system governed by Hamiltonian dynamics to equilibrium correlation functions of the current. The molecular chaos assumption employed in the Boltzmann transport equation (see section 2.6 ) is not applied here, therefore, the GK formulation can be applied for both dilute and dense systems. When transportation takes place in a system such as heat conduction, the change of the Hamiltonian is taken as a perturbation to the initial equilibrium Hamiltonian.  Thus, the external disturbance which perturbs the initial equilibrium Hamiltonian leads to a small disturbance of the particle density function mentioned above. Linear response theory facilitates establishing the linear relation via the transport coefficient between the external force field and the flux arising due to the disturbance to the particle distribution function. The works of Green [95] and Kubo [96], [97] led to an exact expression based on the fluctuationdissipation theorem for the linear transport coefficient in terms of time-dependent equilibrium fluctuations of flux.  In the case of the heat transport process, the external force field becomes the temperature gradient even though it is internal, while the flux becomes the heat flux vector. Therefore, using the GK formulation the thermal conductivity coefficient is expressed in terms of autocorrelation of the heat flux vector which fluctuates around zero at equilibrium. Therefore, the GK method is based on equilibrium molecular dynamics and the parameters described in the following equations are evaluated using MD simulations. The heat flux vector is defined as, where V, r_i and E_i stand for volume, i^th atom‟s position vector and total energy, respectively. The total energy E_i is the summation of both kinetic and potential energies as given in the following equation wherem_i,v_i and  are mass, velocity of i^th atom and pair-wise interaction of i^th atom with j^th atom when they are apart from the distance r_ij, respectively. The lattice thermal conductivity tensor is given by where K_αβ is the rank-two thermal conductivity tensor. A large body of works can be found which used the GK formulation to evaluate the thermal conductivity of different materials, for example, LJ argon [98], [99], diamond [100] , crystalline Si [92], [101], amorphous Si [72] and Si nanowire [102]. The GK formulation can yield all the components of  the thermal conductivity tensor following a single simulation in contrast to the direct method based on NEMD. However, one of the drawbacks of this method is the slow decay of the heat flux autocorrelation function for higher thermal conductivity material compared to lower ones [102].

Lattice Dynamics

In this section some important concepts in lattice dynamics such as normal mode coordinates, dynamical matrix, Hamiltonian etc. are presented in the limit of harmonic approximation where the potential energy function is approximated to a quadratic form, therefore, the restoring force is proportional to the square of the displacement from the equilibrium position.

Modal Analysis of the Thermal Conductivity of Nanowires

Introduction

Investigating the phonon-driven heat transport mechanism of quasi one dimensional (1D) structures such as nanowires, nanotubes and nanoribbons is important because of their potential uses in thermoelectric (TE) applications [40]. The key differences of thermal transport mechanism in low dimensional structures over their bulk counterparts emerges due to phonon confinement and extra phonon-surface scatterings [8], [55], [112]. The origin of the former effect is connected with the relative dimensions of characteristic length of the structure and phonon wavelength [112], [113]. This confinement effect of quasi 1D structures modifies the phonon spectrum, thereby increasing the dominance of long wavelength phonons and reducing phonon group velocities [87], [114]. These distinctive characteristics of phonons give rise to anomalous thermal behaviours over bulk materials. For instance, a universal quantum conductance, quantized limiting value for the conductance, was experimentally demonstrated [115] for a 1D ballistic channel at low temperature as ?²?_?²?/ Furthermore, the problem of “long wavelengthphonons” and convergence issues of quasi 1D thermal conductivity (TC) with the total length have long been the subject of debate. As discussed in Appendix A, the combined effects of 1D Brillouin zone (BZ), linear dispersion and first-order 3-phonon scatterings inevitably lead to the divergence of TC of 1D structures. Several studies [56], [58], [89], [102], [116]–[121] addressed the convergence issue of quasi 1D structures with different material assemblies and techniques which sometimes led to contradictory conclusions. These contradictory conclusions would be consequences of the material structures analysed and the limitations of the techniques employed in the studies. The purpose of this chapter partly is to shed some light on this convergence issue using the techniques based on the Boltzmann transport equation (BTE) in the limit of single mode relaxation time approximation (SMRTA). Among the different varieties of quasi 1D structures, our study is limited for nanowires. In this study, the TC of nanowires is investigated employing normal mode decomposition (NMD) simulations [81], [107], [122] which are based on BTE-SMRTA and subsequently, equilibrium molecular dynamics. The Green-Kubo (GK) method is used for the verification of the results given by NMD. For the investigation, Lennard-Jones (LJ) materials are used to examine thermal transport of solids with a significantly more reduced computational demand than real materials such as Si and Ge [61], [81], [107], [122]. LJ potential is considered as a suitable inter-atomic potential to evaluate the physical and thermal properties of Argon including solid phase. Tretiakov and Scandolo [99] demonstrated that TC evaluated for solid Argon by MD simulations using LJ potential is in good agreement with experimental values. By varying energy depth (ε) and zero-potential pair separation (σ), materials with different properties (soft or hard) can be derived. LJ potential has been widely used to investigate thermal performance of materials due to less computational demand than real materials. LJ potential is popular among researchers for comprehensive analysis of material properties including thermal performance. Moreover the investigation of the non-monotonic trends in TC of nanowires with crosssectional width is another case of this chapter which was addressed before [87], [102], [121] with seemingly different explanations. For example, Karamitaheri et al. [121] attributed the increasing TC trend as cross-sectional width is decreased to the emergence of low frequency modes while Zhou et al. [102] attributed that phenomenon to the dominance of the normal process of acoustic phonons which leads to phonon hydrodynamic flow. We revisit the problem using the system derived from LJ materials mentioned above and suggest possible causes for that non-monotonic trend in the limits of the computational procedure employed here. Seeking further reduction of TC of nanowires via doping or alloying the entire structure with heavier mass atoms has been tested before [123], [124] experimentally and computationally. Even though significant reduction of TC is achieved, it might adversely affect the electrical conductivity of the structure, and thereby the TE figure of merit [125]. Consequently, alloying or doping was confined only to the shell region [125]–[127], thus, the electrons are free to move through the core region while phonons are disrupted with the shell modifications. Despite the significant reduction of TC of shell alloyed/doping nanowires being attributed to the mode localisation [125], [128] depicted by the mode participation ratio (PR), comprehensive analysis is still lacking.Thus, in this chapter, we develop computational techniques which facilitate analyzing and decomposing the contribution of different vibrational modes to the TC of nanowires in which shell is alloyed with heavier mass atoms. This chapter is organized as follows. Firstly, the theoretical and computational formalisms of the NMD method are presented including lattice dynamics calculations at the harmonic limit. Subsequently, the NMD implementation is demonstrated and compared with GK simulations using pristine and superlattice nanowires with a discussion pertaining to the convergence of TC of the structures considered here. Subsequently, the non-monotonic behaviour of TC of nanowires is addressed and an attempt is made to explain the trend following the insight given by the NMD method. Then we investigate the effect of shell mass disorder on the TC of nanowires by developing techniques to decompose the influence of different vibrational modes such as propagons and diffusions.

1. Introduction
1.1 Nanostructures
1.2 Nanostructures in Thermoelectric Applications
1.3 Heat Transfer in Nanostructures
1.4 Motivation
1.5 Outline of The Thesis
1.6 List of Articles Prepared
2. Theoretical Background
2.1 Outline
2.2 Molecular Dynamics in Nanoscale Heat Transfer
2.3 Direct Method for Evaluating Thermal Conductivity
2.4 Green-Kubo Method for Evaluating Thermal Conductivity
2.5 Lattice Dynamics
2.6 Boltzmann Transport Equation
2.7 Normal Mode Decomposition
2.8 Allen-Feldman Theory
3. Modal Analysis of the Thermal Conductivity of Nanowires
3.1 Introduction
3.2 Theoretical and Computational Formalism
3.3 Results and Discussion
4. Thermal Performance of Si/Ge Random Layer Nanowires
4.1 Introduction
4.2 Simulation Procedures
4.3 Results and Discussion
5. Thermal Performance of  Nanotwinned Random Layer Structures
5.1 Introduction
5.2 Computational Procedure
5.3 Results and Discussion
6. Conclusions, Limitations and Future Work
6.1 Conclusions
6.2 Limitations
6.3 Future Work
GET THE COMPLETE PROJECT
Phonon-Governed Heat Conduction in Nanostructures