CC BY
0
EDN: VQWXRW УДК 534.16
MD Investigations of of Heat Flow throw Interfaces in 1D Systems
Alexander S. Fedorov* Maxim A. Visotin
Kirensky Institute of Physics Federal Research Center KSC SB RAS Krasnoyarsk, Russian Federation Siberian Federal University Krasnoyarsk, Russian Federation
Oleg A. Sosedkin Egor V. Eremkin^
Siberian Federal University Krasnoyarsk, Russian Federation
Received 10.11.2022, received in revised form 03.02.2023, accepted 20.04.2023
Abstract. Molecular dynamic calculations (MD) of heterogeneous 1D periodical systems are presented. It is proposed the new technique of direct calculations of thermal conductivity, where there is only one thermostat in one piece of unit cell as well as another piece where artificial friction forces act on atoms. With the help of this scheme, calculations of 1D heterogeneous systems having regions with atoms of different atomic masses are presented. It is shown that the difference in atomic masses in adjacent regions of the systems leads to a significant temperature jump at interfaces between these regions. This temperature jump exists independently of the mass ratio on both sides of the interface.The reasons for these jumps are discussed. It is also shown that, by changing the alternation of regions with different masses of atoms, it is possible to reduce the total thermal conductivity of the system by several times. On the base of these results, we can hope that for three-dimensional structures also, the thermal conductivity can be significantly reduced.
Keywords: molecular dynamic, thermal conductivity, interface, temperature jump.
Citation: A.S. Fedorov, M.A. Visotin, O.A. Sosedkin, E.V. Eremkin, MD Investigations of of Heat Flow throw Interfaces in 1D Systems, J. Sib. Fed. Univ. Math. Phys., 2023, 16(3), 385-396. EDN: VQWXRW.
Introduction
In recent years, interest in studying the thermal conductivity of various materials, especially nanomaterials, has sharply increased. This is explained by the various possible applications of these properties in practice.
In particular, thermoelectric converters, based on thermoelectric rectification effects, are widely used to generate electrical currents using thermal sources [1,2]. For good efficiency of the thermoelectric converters applying, it is necessary to create thermoelectric materials having a great dimensionless thermoelectric figure of merit ZT = S2aT/x, where T is temperature, S is the Seebeck coefficient (also known as thermopower), a is the electrical conductivity, and x is the thermal conductivity [3]. As it can be seen, to increase ZT value, one can decrease the thermal
*alex99<8iph.krasn.ru https://orcid.org/0000-0002-7911-3301 teremkin.e.v@yandex.ru https://orcid.org/0000-0003-2118-3435 © Siberian Federal University. All rights reserved
conductivity x- At the moment, this method of the thermal conductivity decreasing is considered the method of "phonon engineering" as one of the most promising- Controlling this value of nanostructured or reduced dimensional materials could enable new generations of thermoelectric materials and thermal insulators. Recently it was shown [4,5] that thermal conductivity of multilayer thin films of metals and oxides can be reduced at several times in comparison with thermal conductivity of the homogeneous amorphous oxide. In [6] authors experimentally and theoretically investigated the thermal conductivity of disordered thin films of the layered crystal WSe2 and established it can be as small as 0-05 W/(mK, that is minima value among ever observed fully dense solids. Potentially very promising elements are "thermal rectifiers" [7], the thermal analogue of the electrical diode, which let the heat flow in one direction easier than in the opposite. Thermal rectification is impossible in the framework of the linear theory of thermal conductivity, based on consideration of the propagation of quants of thermal excitation (phonons) through the system. But it is possible in the presence of non linearity (anharmonic-ity) and therefore symmetry breaking. These an harmonic effects have been investigated in some theoretical and experimental works [7-13]. Also heat conductivity was discussed as base for different thermal components such as thermal memory [14], thermal transistors [8,15,16] and even thermal logic gates [14]. These elements have potential applications in many technologies including thermal management systems for nanoelectronics, energy-harvesting modules and thermal circuits that use heat instead of electricity for different operations [17-19].
The main objective of the work is to study general regularity of interface influence on the thermal conductivity and temperature distribution in 1D dielectric systems contained adjacent sectors having different properties. In dielectric systems, where the valence electrons are tightly bound to the atomic nuclei, lattice vibrations are dominate the thermal transport. At that the thermal transport analysis is typically done in the momentum space by introducing of phonon excitations which corresponds to a set of independent harmonic oscillators and characterized by a wave vector q and a zone number n. It is usually done by calculating the force matrix using quantum chemical calculations and the frozen phonon method. In this method the force constant matrix is calculated by displacing each atom in a unit cell and calculating the resulting forces on every other atom.
Unfortunately, the only harmonic vibrations are described by this technique, so this theory is correct only at zero temperature. As temperature increases, an-harmonic effects of phonon interaction, which are difficult to model theoretically, become important.
Another problem is the description of heat transport in complex structures, or in structures containing any defects (impurities, interfaces between different layers, etc.). In such systems, thermal transport even at low temperatures is also difficult to describe by finding individual and orthogonal to each other phonons due to the complexity of their calculation.
The molecular dynamic (MD) simulations are a suitable tool for the analysis of thermal flow analysis at finite temperature in complex systems. At that these simulations allow for the natural inclusion of an-harmonic effects of phonon-phonon scattering. By ignoring electrons in non ab initio (empirical) MD scheme and description of atomic dynamics only the computational demands of MD are greatly reduced. Due to the simplicity of the method, MD simulations run in parallel and can handle systems with millions or even billions of atoms [20] on on the most powerful supercomputers.
For many practical applications, it is necessary optimize the heat flux in one direction only. For that it is possible to produce a superlattices consisting of different layers perpendicular to this direction.
And although all such superlattices are three-dimensional, it can be assumed that for modeling of heat flow it is possible to substitute 3D structures to 1D ones, where the thermal properties of each homogeneous superlattice layer can be described by the corresponding properties of 1D sector. For example, in [21] it was shown the eigen vectors of phonons vibrations in monatomic
3D case can be reduced to the 1D case of one-component wave vector of vibrations k = [0,0, k]
Mujxi = D(k)xi; D(k)=J2 D(R)e-ikR, (1)
R
where x is eigenvector of atoms vibrations for made i, J2 is eigenvalue (frequency) for made i, D(k) is dynamic matrix for wave vector k, D(R) is force matrix.
1. Systems description
In this work it is investigated the thermal conductivity of periodic one-dimensional systems with help of empirical molecular dynamic (MD) simulations. More specifically, 1D systems including sectors having different properties were investigated. Without reducing the generality, only the masses of atoms belonging to different sectors are changed without changing the elastic constants of the interacting neighboring atoms.
In more detail, the systems consist of periodically repeated chains of 320, 640, and 1280 argon atoms were investigated, where boundary atoms interacted with the atoms of adjacent chains. Atom interaction was described by the Lennard-Jones (LJ) potential U(r) = 4e((a/r)6-(a/r)12), where the potential parameters e = 103.2 eV and a = 3.405 A_parameters were taken from [22].
In all calculations each periodically repeated chain of length N atoms consisted of the sector with an atomic mass m1 of length N/4, the sector with an atomic mass of m2 length N/2 and the sector with atomic mass m1 of length N/2, see Fig. 1.
Fig. 1. Geometry of periodic systems consisting of sectors with different atomic mass
MD equations of atoms dynamics were integrated using a time step AT = 0.5 fsec and a number of time steps Ntime « 2 * 105. Changes in the atoms coordinates r and velocities v were carried out using Verlet algorithm in speed form [23]:
r(t + At)
, At
v(t + T)
a(t)
r(t)+ v(t)At + 1 a(t)At2,
v(t) + v(t)At + 1 a(t)At, VU (r(t))
2. New technique of thermo-conductivity direct calculations
Here we propose a new technique of thermo-conductivity calculation. Instead of the standard direct calculation method, where the temperature difference is maintained in two pieces of the system by two thermostats having temperatures T and T2, there is only one thermostat in this scheme. Instead of another thermostat in the second piece, it introduces the forces of viscous friction Ffr,i acting on all atoms i in this region. We postulate these forces are proportional to the velocity of the atoms: Ffr,i = aVi. These forces leads to energy dissipation and decreasing the average kinetic energy and the temperature T2 of atoms in the second piece.
m
To check this technique, we made MD simulations of the periodical chain of argon atoms, where the unit cell have included 320 atoms. In the first half of the unit cell there was a piece included 8 consecutive atoms where the thermostat maintained the T temperature. In the another half of the unit cell the friction forces described by the coefficient of viscous friction equal to a = 0,05 m/t acted on all atoms. For comparison we have made MD simulations of the system by standard direct calculation method using two thermostats having temperatures T and T2, see Fig. 2. As shown in the figure, the temperature profile of the proposed modified scheme using one thermostat quantitatively coincides with the profile of the classical direct method where the cold and hot pieces are controlled by two thermostats.
Fig. 2. Comparison of the temperature profiles calculated by standard direct method using two thermostats (red line) and modified scheme using one thermostat (black dotted line)
3. Investigations of heterogeneous periodical system thermal conductivity
To determine the optimal configuration of the system, a study is made of the dependence of the thermal conductivity coefficient on the masses m1 and m2. The results are presented in the form of a surface (Fig. 3).
From the graph (Fig. 3) it is seen that in order to obtain the minimum thermal conductivity, the masses should be as different as possible.
4. Temperature profile
In view of the fact that the temperature of the structure is different for different mass layers, the temperature surfaces were considered: the dependence of temperature on the atom number of the cell and the mass ratio of the layers m\/m2.
It can be seen from Fig. 4 that upon a transition from one mass layer to another, a jump in temperature is observed for all mass ratios m\/m2, except for the case m\/m2 = 1. The temperature in atomic layers of the same mass varies linearly with the exception of the thermostat region and the region of absorption of atomic energy by resistance forces. The larger the mass ratio m\/m2 differs from 1, the more temperature jumps.
ml
Fig. 3. dependence of thermal conductivity on atomic masses in periodic segments of the system
Fig. 4. Dependence of temperature on the number of atoms N on the ratio of atomic masses in periodic layers
5. Correlations of atomic displacements in the structure
To represent the motion of atoms relative to each other in a two-layer structure with the most different masses of layers (mi = 2, m2 = 50), correlation coefficients were calculated
n (■ = 0(rit - fi)(r3t - ) (o\
Corr(i,j)= , (2)
yY,u=o (rit - ri)2(rjt - rj)2
where rit, rjt are the coordinates of the atoms i and j at time t; ri =Y1, rit/Nt, fj =Y1 r'jt/Nt is
t t
the average value of atomic coordinates.
j
Fig. 5. Dependence of the correlation coefficients between atoms i and j, where i varies from 0 to 2n, j varies from i to i + 2n
In view of the symmetry of the unit cell, the correlation coefficient was calculated between each atom of the first half of the cell and 2n subsequent atoms. Thus, we exclude repetitions. Fig. 3 shows that the correlation of adjacent atoms is positive and decreases as they move away from each other. Also, the values of the correlation coefficient at the boundary of the barrier have a maximum value, and zero is reached at the largest possible distance between the atoms relative to the other fixed positions of atom i. This suggests that the atoms located near the barrier correlate with a large number of atoms around themselves compared with atoms located in the middle of a layer of the same mass.
6. Autocorrelation of atomic displacements in the structure
An autocorrelation function is used to study the behavior of phonons at the mass boundary
/w
r(t)r*(t - T)dT, (3)
-w
where r(t) is a function of the coordinate of time t, r(t-T) floor is a function shifted by the value of T The autocorrelation shows the behavior of phonons in the system. We represent the result as a 3D surface.
The graph (Fig. 6) shows pronounced fluctuations. Long-wave oscillations are visible on the entire surface, short-wave ones are observed only in a layer of mass m1. It can be assumed that they do not pass. We verify this by cutting off the contribution of long-wave oscillations. The result is represented as a 3D surface.
From the graph in the Fig. 7 it is seen how the phonons pass from the layer m1 to the layer m2. Long-wave phonons pass with increasing amplitude, but preserving the period. Short-wave phonons do not pass into the zone with mass m2.
To test the effect of the size of the system on the behavior of atoms, a series of tests was carried out. Each iteration, the length of the system became larger. The results in the form of autocorrelation functions are displayed on the graphs (Fig. 8).
Fig. 6. Dependence of the autocorrelation function on the number of the atom i, which is close at the mass barrier at T = 300 K
3000
2000
time (50*fsec) 600 60 65 70
Fig. 7. Autocorrelation for atoms at the mass barrier without taking into account the contribution of long-wave phonons at T = 300 K
Analyzing them, one can see that atomic vibrations are observed in the region of the mass barrier. As the size of the system increases, the frequency becomes lower, but the amplitude increases. Based on this, it can be assumed that, at small cell sizes, long-wave oscillations are partially reflected on the barrier, and as a result, the reflected and reflected oscillations are superimposed. For large cell sizes, these oscillations lose more energy to the barrier, are reflected weakened, or are not reflected at all.
c)
Fig. 8. Autocorrelation surface function for systems of different lengths N at T = 300 K: a) N = 640, b) N = 1280, c) N = 2560
7. Green-Kubo method
The Green-Kubo method [24,25] is an equilibrium method for calculating the thermal conductivity coefficient k. System always in linear response mode, because the simulations are done in equilibrium and there is no driving force. We use the Kubo method to calculate k of our one-dimensional system and compare it with our new direct method scheme.
The Green-Kubo expression is defined through the equilibrium current-current autocorrelation function
1 r TM
k^u(TM) = qkbT2 j <mt)jv(°) >dт, (4)
where Q is the system volume, kB is the Boltzmann constant, T is the system temperature, and the angular brackets is the heat current autocorrelation function (HCACF). Since the simulation is performed for discrete steps of MD length dT the integration of Eq. (4) goes into the summation.
Including time averaging, what we actually compute is
M N-m
(tm) = q^zY^ (N - m) 1 Y J^(n + m) Ju(n) (5)
B m=1 n=1
where tm is given by MAt and JM(n + m) is ^th component of the heat current at MD timestep n + m. The number of integration steps M should be less than the steps of the MD simulation N to obtain a good averaging value. For example, in the presented work, the number of MD steps of the simulation N = 2 x 106 (1 ns for a 0.5 fs MD step), while the number of integration steps M = 2 x 105 (tm = 100 ps).
The heat current can be found as
J = dt ? n(t)£i (t), (6)
i
where ri (t) is the time-dependent coordinate of atom i and ei(t) is the local energy. For a pair potential, where the total potential energy is written in terms of the pairwise interactions, the local energy can be fou
£i = ^ miv1 + 2 Y u(rij). (7)
It is easy to show that for a given definition of local energy, the thermal current is defined as
J(t) = X)Vi£i +12 rij(FijVj), (8)
i i,j i=j
where Fij is the force on atom i due to its neighbor j from the pair potential.
The Green-Kubo method has many advantages over the direct method described in [26]. But for the case of one-dimensional superlattices, the direct method is more suitable, since in the GK method it is impossible to predict the boundary thermal resistance. As seen in Fig. 9, the GK method requires a lot of time steps for the HCACF to converge to zero, and the time-dependent thermal conductivity "plateau". And also, when the ratio of the masses of different layers increases, the difference between the thermal conductivity of the two different sets of masses decreases, because the HCACF makes fluctuations around zero, comparable to the difference between these coefficients.
t, pc
Fig. 9. Thermal conductivity at T = 300 K for different mass layers found by integrating the current-current correlation function using Eq. (5) as a function of the upper integration limit tm
Conclusions
Within the framework of the method of molecular dynamics calculations, a new modification of direct calculations of thermal conductivity was implemented. Using this scheme, calculations are presented for one-dimensional heterogeneous systems having regions with atoms of different atomic masses. Calculations show that the difference in atomic masses in neighboring regions of systems leads to a significant jump in temperature. This jump in temperature at the interfaces can be observed regardless of the mass ratio on either side of the interface. It is also shown that, depending on the order of layers with atoms of different masses, the thermal conductivity of systems can differ several times. On the basis of these results, it can be hoped that thermal conductivity can be significantly reduced for three-dimensional structures as well.
References
[1] D.M.-T.Kuo, Y.-c.Chang, Thermoelectric and thermal rectification properties of quantum dot junctions, Physical Review B, 81(2010), no. 20, 205321.
DOI: 10.1103/PhysRevB.81.205321
[2] Z.H.Zhang, Y.S.Gui, L.Fu, X.L.Fan, J.W.Cao, D.S.Xue, P.P.Freitas, D.Houssameddine, S.Hemour, K.Wu, C.-M.Hu, Seebeck Rectification Enabled by Intrinsic Thermoelectrical Coupling in Magnetic Tunneling Junctions, Physical Review Letters, 109(2012), no. 3, 37206. DOI: 10.1103/PhysRevLett.109.037206
[3] A.F.Ioffe, Semiconductor thermoelements and thermoelectric cooling, Infosearch Limited, London, 1957.
[4] R.M.Costescu, Ultra-Low Thermal Conductivity in W/Al2O3 Nanolaminates, Science, 303(2004), no. 5660, 989-990. DOI: 10.1126/science.1093711
[5] Y.S.Ju, M.T.Hung, M.J.Carey, M.C.Cyrille, J.R.Childress, Nanoscale heat conduction across tunnel junctions, Applied Physics Letters, 86(2005), no. 20, 1-3.
DOI: 10.1063/1.1931827
[6] C.Chiritescu, D.G.Cahill, N.Nguyen, D.Johnson, A.Bodapati, P.Keblinski, P.Zschack, Ul-tralow thermal conductivity in disordered, layered WSe2 crystals., Science (New York, N.Y.), 315(2007), no. 5810, 351-353. DOI: 10.1126/science.1136494
[7] M.Terraneo, M.Peyrard, G.Casati, Controlling the energy flow in nonlinear lattices: A model for a thermal rectifier, Physical Review Letters, 88(2002), no. 9, 943021-943024.
DOI: 10.1103/PhysRevLett.88.094302
[8] B.Li, L.Wang, G.Casati, Thermal diode: Rectification of heat flux, Physical Review Letters, 93(2004), no. 18. DOI: 10.1103/PhysRevLett.93.184301
[9] D.Segal, A.Nitzan, Spin-boson thermal rectifier,Physical Review Letters,94(2005), no. 3. DOI: 10.1103/PhysRevLett.94.034301
[10] C.W.Chang, D.Okawa, A.Majumdar, A.Zettl, Solid-state thermal rectifier, Science (New York, N.Y.), 314(2006), no. 5802, 1121-1124. DOI: 10.1126/science.1132898
[11] D.Segal, Single mode heat rectifier: Controlling energy flow between electronic conductors, Physical Review Letters, 100(2008), no. 10. DOI: 10.1103/PhysRevLett.100.105901
[12] A.L.Cottrill, M.S.Strano, Analysis of Thermal Diodes Enabled by Junctions of Phase Change Materials, Advanced Energy Materials, 5(2015), no. 23. DOI: 10.1002/aenm.201500921
[13] S.Wang, A.L.Cottrill, Y.Kunai, A.R.Toland, P.Liu, W.-J.Wang, M.S.Strano, Microscale solid-state thermal diodes enabling ambient temperature thermal circuits for energy applications, Physical chemistry chemical physics: PCCP, 19(2017), no. 20, 13172-13181. DOI: 10.1039/c7cp02445b
14] L.Wang, B.Li, Thermal logic gates: Computation with phonons, Physical Review Letters, 99(2007), no. 1. DOI: 10.1103/PhysRevLett.99.177208
15] P.Ben-Abdallah, S.A.Biehs, Near-field thermal transistor, Physical Review Letter, 112(2014), no. 4. DOI: 10.1103/PhysRevLett.112.044301
16] K.Joulain, J.Drevillon, Y.Ezzahri, J.Ordonez-Miranda, Quantum Thermal Transistor, Physical Review Letters, 116(2016), no. 20. DOI: 10.1103/PhysRevLett.116.200601
17] N.Li, J.Ren, L.Wang, G.Zhang, P.Hänggi, B.Li, Colloquium: Phononics: Manipulating heat flow with electronic analogs and beyond, Reviews of Modern Physics, 84(2012), no. 3, 1045-1066. DOI: 10.1103/RevModPhys.84.1045
18] S.Narayana, Y.Sato, Heat flux manipulation with engineered thermal materials, Physical Review Letters, 108(2012), no. 21. DOI: 10.1103/PhysRevLett.108.214303
19] M Maldovan, Sound and heat revolutions in phononics, Natur, 503(2013), no. 7475, 209-217. DOI: 10.1038/nature12608
20] K.Kadau, T.C.Germann, P.S.Lomdahl, Large-scale molecular-dynamics simulation of 19 billion particles, International Journal of Modern Physics, 15(2004), no. 1, 193-201. DOI: 10.1142/S0129183104005590
21] W.A.Harrison, Solid State Theory, Dover Books on Physics, Dover Publications, 2012.
22] Y.Chen, J.R.Lukes, D.Li, J.Yang, Y.Wu, Thermal expansion and impurity effects on lattice thermal conductivity of solid argon, The Journal of Chemical Physics, 120(2004), no. 8, 3841-3846. DOI: 10.1063/1.1643725
23] L.Verlet, Computer "experiments" on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules, Physical Review, 159(1967), no. 1, 98-103. DOI: 10.1103/PhysRev.159.98
24] M.S.Green, Markoff random processes and the statistical mechanics of time-dependent phenomena. II. Irreversible processes in fluids, The Journal of Chemical Physics, 22(1954), no. 3, 398-413. DOI: 10.1063/1.1740082
25] R.Kubo, Statistical'Mechanical Theory of Irreversible Processes. I. General Theory and Simple Applications to Magnetic and Conduction Problems, Journal of the Physical Society of Japan, 12(1957), no. 6, 570-586. DOI: 10.1143/JPSJ.12.570
26] A.J.McGaughey, M.Kaviany, Thermal conductivity decomposition and analysis using molecular dynamics simulations. Part I. Lennard-Jones argon, International Journal of Heat and Mass Transfer, 47(2004), no. 8-9, 1783-1798. DOI: 10.1016/j.ijheatmasstransfer.2003.11.002
МД исследования границ раздела теплового потока в оДномерных системах
Александр С. Федоров Максим А. Высотин
Институт физики им. Л. В. Киренского СО РАН Федеральный исследовательский центр КНЦ СО РАН Красноярск, Российская Федерация Сибирский федеральный университет Красноярск, Российская Федерация
Олег А. Соседкин Егор В. Еремкин
Сибирский федеральный университет Красноярск, Российская Федерация
Аннотация. Представлены молекулярно-динамические расчеты (МД) гетерогенных одномерных периодических систем. Предлагается новая методика прямых расчетов теплопроводности, при которой в одном элементе элементарной ячейки находится только один термостат, а в другом элементе действуют силы искусственного трения на атомы. С помощью этой схемы представлены расчеты одномерных гетерогенных систем, имеющих области с атомами разной атомной массы. Показано, что различие атомных масс в соседних областях систем приводит к значительному скачку температуры на границах раздела между этими областями. Этот скачок температуры существует независимо от отношения масс по обе стороны от границы раздела. Обсуждаются причины этих скачков. Также показано, что, изменяя чередование областей с разной массой атомов, можно в несколько раз уменьшить общую теплопроводность системы. На основании этих результатов можно надеяться, что и для трехмерных структур теплопроводность может быть значительно снижена.
Ключевые слова: молекулярная динамика, теплопроводность, скачок температуры интерфейса.
