Научная статья на тему 'МЕТОДЫ ВЫЧИСЛЕНИЯ РЕШЕТОЧНОЙ ТЕПЛОПРОВОДНОСТИ МЕТАЛЛОВ ПРИ ВЫСОКИХ И НИЗКИХ ТЕМПЕРАТУРАХ'

МЕТОДЫ ВЫЧИСЛЕНИЯ РЕШЕТОЧНОЙ ТЕПЛОПРОВОДНОСТИ МЕТАЛЛОВ ПРИ ВЫСОКИХ И НИЗКИХ ТЕМПЕРАТУРАХ Текст научной статьи по специальности «Физика»

CC BY
81
9
i Надоели баннеры? Вы всегда можете отключить рекламу.
Область наук
Ключевые слова
МЕТОД НЕРАВНОВЕСНОЙ МОЛЕКУЛЯРНОЙ ДИНАМИКИ / РАСПРЕДЕЛЕНИЕ ГАУССА / РЕШЕТОЧНАЯ ТЕПЛОПРОВОДНОСТЬ / ТЕМПЕРАТУРОПРОВОДНОСТЬ

Аннотация научной статьи по физике, автор научной работы — Саламатов Евгений Иванович, Долгушева Елена Борисовна

В статье рассмотрены вопросы теоретического описания решеточного вклада в теплопроводность металлов на примере алюминия. Проведен численный анализ решеточной теплопроводности, опирающийся на метод неравновесной молекулярной динамики. Получена температурная зависимость коэффициента решеточной теплопроводности алюминия, согласующаяся с результатами первопринципных расчетов при высоких температурах. При низких температурах получено значение коэффициента температуропроводности близкое по значению к экспериментальным результатам, знание которого позволяет вычислить фононную теплопроводность.

i Надоели баннеры? Вы всегда можете отключить рекламу.

Похожие темы научных работ по физике , автор научной работы — Саламатов Евгений Иванович, Долгушева Елена Борисовна

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
i Надоели баннеры? Вы всегда можете отключить рекламу.

METHODS FOR CALCULATING THE LATTICE THERMAL CONDUCTIVITY OF METALS AT HIGH AND LOW TEMPERATURES

The Molecular Dynamics (MD) method seems to be the most promising for determining the lattice contribution to the overall thermal conductivity of metals and metal alloys. In the paper, the method is used for studying the lattice thermal conductivity of aluminum at high and low temperatures with a proven potential. It is shown that at high temperatures, standard algorithms are more convenient for calculating the lattice thermal conductivity coefficient. In this case, the thermal conductivity coefficient is calculated using the Fourier equation, and the MD calculation is used to simulate a stationary non-equilibrium state with a linear temperature gradient at a length comparable to the size of the calculated cell. This approach is shown to give the values of the lattice thermal conductivity coefficient in good agreement with the results of the first-principles calculations. With a decrease in the size of the base crystallite, the thermal conductivity coefficient decreases due to the depletion of the low-frequency section of the phonon spectrum, the contribution of which to thermal conductivity becomes insignificant with increasing temperature. At high temperatures, the thermal conductivity coefficient does not depend on the size of the crystallite and coincides with that obtained from the first-principles calculations. To calculate the thermal conductivity at low temperatures, a new method is proposed, based on the homogeneous heat equation given on an infinite line. In this case, the task of the MD method is to obtain a stationary non-equilibrium temperature distribution in the system in the form of a Gaussian curve, which is the fundamental solution of the equation. The approximation of the system relaxation from the non-equilibrium to equilibrium state allows finding the thermal diffusivity coefficient related to the thermal conductivity coefficient. The test calculations carried out for a thin aluminum film at low temperatures with different initial conditions showed that the obtained thermal diffusivity coefficient does not depend on the parameters of the initial Gaussian distribution, which suggests the suitability of the proposed method for studying the lattice thermal conductivity.

Текст научной работы на тему «МЕТОДЫ ВЫЧИСЛЕНИЯ РЕШЕТОЧНОЙ ТЕПЛОПРОВОДНОСТИ МЕТАЛЛОВ ПРИ ВЫСОКИХ И НИЗКИХ ТЕМПЕРАТУРАХ»

https://doi.org/10.15350/17270529.2021.4.44

УДК 538.913+534-6

Методы вычисления решеточной теплопроводности металлов при высоких и низких температурах

Е. И. Саламатов, Е. Б. Долгушева

Удмуртский федеральный исследовательский центр УрО РАН, Россия, 426067, Ижевск, ул. Т. Барамзиной, д. 34

Аннотация. В статье рассмотрены вопросы теоретического описания решеточного вклада в теплопроводность металлов на примере алюминия. Проведен численный анализ решеточной теплопроводности, опирающийся на метод неравновесной молекулярной динамики. Получена температурная зависимость коэффициента решеточной теплопроводности алюминия, согласующаяся с результатами первопринципных расчетов при высоких температурах. При низких температурах получено значение коэффициента температуропроводности близкое по значению к экспериментальным результатам, знание которого позволяет вычислить фононную теплопроводность.

Ключевые слова: метод неравновесной молекулярной динамики, распределение Гаусса, решеточная теплопроводность, температуропроводность.

И Елена Долгушева, e-mail: elena@udman.ru

Methods for Calculating the Lattice Thermal Conductivity of Metals at High and Low Temperatures

Evgeny I. Salamatov, Elena B. Dolgusheva

Udmurt Federal Research Center UB RAS (34, T. Baramzina St., Izhevsk, 426067, Russian Federation)

Summary. The Molecular Dynamics (MD) method seems to be the most promising for determining the lattice contribution to the overall thermal conductivity of metals and metal alloys. In the paper, the method is used for studying the lattice thermal conductivity of aluminum at high and low temperatures with a proven potential. It is shown that at high temperatures, standard algorithms are more convenient for calculating the lattice thermal conductivity coefficient. In this case, the thermal conductivity coefficient is calculated using the Fourier equation, and the MD calculation is used to simulate a stationary non-equilibrium state with a linear temperature gradient at a length comparable to the size of the calculated cell. This approach is shown to give the values of the lattice thermal conductivity coefficient in good agreement with the results of the first-principles calculations. With a decrease in the size of the base crystallite, the thermal conductivity coefficient decreases due to the depletion of the low-frequency section of the phonon spectrum, the contribution of which to thermal conductivity becomes insignificant with increasing temperature. At high temperatures, the thermal conductivity coefficient does not depend on the size of the crystallite and coincides with that obtained from the first-principles calculations. To calculate the thermal conductivity at low temperatures, a new method is proposed, based on the homogeneous heat equation given on an infinite line. In this case, the task of the MD method is to obtain a stationary non-equilibrium temperature distribution in the system in the form of a Gaussian curve, which is the fundamental solution of the equation. The approximation of the system relaxation from the non-equilibrium to equilibrium state allows finding the thermal diffusivity coefficient related to the thermal conductivity coefficient. The test calculations carried out for a thin aluminum film at low temperatures with different initial conditions showed that the obtained thermal diffusivity coefficient does not depend on the parameters of the initial Gaussian distribution, which suggests the suitability of the proposed method for studying the lattice thermal conductivity.

Keywords: molecular dynamics method, Gaussian distribution, lattice thermal conductivity, thermal diffusivity.

El Elena Dolgusheva, e-mail: elena@udman. ru

ВВЕДЕНИЕ

Процессы переноса тепла решеткой в металлах и металлических соединениях до сих пор остаются мало изученными, поскольку в большинстве металлов электронная составляющая теплопроводности существенно преобладает над фононной. К настоящему времени нет точного количественного описания зависимости решеточной (фононной) теплопроводности от температуры в металлах и интерметаллидах. Оценка решеточной теплопроводности по электросопротивлению затруднительна из-за чувствительности эксперимента к дефектам решетки при низких температурах. Строгие первопринципные расчеты ограничены нулевой температурой [1], а если вводится температура, то через какие-либо параметры. Перспективным методом вычисления фононной теплопроводности является метод неравновесной молекулярной динамики (NEMD) [2, 3]. В настоящее время для молекулярных и атомных симуляций широко используются стандартные пакеты: XMD, LAMMPS и др. Существует несколько алгоритмов расчета коэффициента теплопроводности, вшитые в пакет LAMMPS [4], но все они корректно работают только при температурах от 100 К и выше. Поэтому требуется развивать методы расчета теплопроводности при низких температурах.

Для оценки теплопроводности металла при низких температурах в работе предлагается оригинальный метод получения коэффициента температуропроводности, зная который можно определить фононную теплопроводность.

Вычисления решеточной теплопроводности были проведены для кристаллитов алюминия. Для описания межатомного взаимодействия использовался хорошо проверенный и апробированный потенциал алюминия из работы [5]. Ранее с этим потенциалом нами были получены различные физические характеристики алюминия: параметры решетки, упругие модули, фононные спектры, тепловое расширение, теплоемкость [6, 7]. Сравнение всех этих величин с экспериментальными данными показало хорошее согласие.

СХЕМА РАСЧЕТОВ ТЕПЛОПРОВОДНОСТИ ПРИ Т > 100 K

Моделирование потока тепла проводилось на кристаллитах с сечением 12*12 (xy) элементарных ячеек (u.c.) ГЦК алюминия и длинами вдоль оси г 48 и 96 ячеек, что соответствует размерам, примерно 5*5*20 nm и 5*5*40 nm. Вдоль всех осей задавались периодические граничные условия. Тепловой поток формировался вдоль направления [001].

В выражении, используемом при расчетах решеточной теплопроводности (kp) [2] :

dl=-k S

dt p AZ

S - усредненное сечение, dQ/dt - поток энергии в единицу времени, ДТ/AZ - градиент температуры. По сути, это закон теплопроводности Фурье, который, вообще говоря, предполагает линейное изменение температуры с расстоянием, а это можно достигнуть только в сильно разряженной среде, например, в газах. Мы же хотим посмотреть решеточную (фононную) теплопроводность металла с помощью пакета LAMMPS. Для расчета коэффициента решеточной теплопроводности был выбран алгоритм "eHEX" [8], основанный на методе неравновесной молекулярной динамики. Этот алгоритм встроен в пакет LAMMPS и позволяет моделировать стационарный тепловой поток в течение временного интервала порядка 100 ns без потерь энергии.

Вначале проводилась оптимизация параметров решетки и минимизация энергии при определенной температуре и атмосферном давлении в условиях NPT-ансамбля. После проведения релаксации, как только устанавливалось равновесное состояние, фиксировался параметр решетки, и уже в условиях NVT-ансамбля в кристаллите создавались зоны "горячего" и "холодного" термостатов путем добавления и отвода одинаковой порции энергии (dQ) на каждом временном шаге. Шаг по времени равен 1 fs. Поскольку граничные условия заданы периодические, следовательно, "горячий" и "холодный" термостаты

dQ = -kp^ AT, (1)

задаются не на границе кристаллита, а на расстоянии одной четвертой части длины от обоих концов кристаллита и, соответственно, между ними расстояние равное половине длины кристаллита. Возникающий градиент температуры между этими термостатами вдоль направления [001] приводит к появлению теплового потока, который контролировался и поддерживался постоянным. Кристаллит разбивался на слои поперек направления распространения теплового потока. Толщина каждого слоя равнялась размеру элементарной ячейки, и в каждом из слоев, содержащем по 576 атомов, вычислялась мгновенная средняя температура по слою. Распределение температуры по слоям отслеживалось в течение определенного периода моделирования до реализации стационарного режима.

На рис. 1 показано пространственное распределение температуры по слоям кристаллита Al длиной 96 u.c. от "горячего" термостата к "холодному" для двух значений температуры 70 K (а) и 700 K (b). Каждая точка соответствует значению температуры, полученной из усредненной кинетической энергии атомов одного слоя. Усреднение температуры в каждом из слоев проводилось на временном интервале равном 1 ns, после установления стационарного режима, который в данных примерах установился через время равное 3 ns. Ошибка - это стандартное отклонение мгновенных значений температуры. Из рисунка видно, что при Т = 70 К температурный профиль очень далек от линейного распределения, для которого собственно можно было бы применять закон Фурье. Особенно сильно это наблюдается непосредственно вблизи расположения термостатов: "горячий" - 24 слой и "холодный" - 72. При температуре Т = 700 К, как видно из рис. 1, b, распределение температуры уже приближается к линейному, хотя небольшие отклонения около "термостатов" все равно есть. Поэтому значения коэффициента решеточной теплопроводности kp вычислялись только на линейном участке распределения температуры (kpj). Линейный участок, полученный аппроксимацией методом наименьших квадратов, на рисунках (а) и (b) показан прямой сплошной линией, он находится между 35 и 60 слоями. Этот интервал был определен при наименьшей из рассмотренных температур, и для всех остальных температур не изменялся.

О 10 20 30 40 50 60 70 80 90 о 10 20 30 40 50 60 70 80 90

Layer Number Layer Number

Рис. 1. Температурные профили в стационарном режиме для кристаллита Al (12x12x96 u.c.) при значениях температуры T = 70 K - (а) и T = 700 K - (b).

Ошибка показывает стандартное отклонение мгновенных значений температуры

Fig. 1. Temperature profiles in stationary mode for Al (12x12x96 u.c.) crystallite at temperature values T = 70 K - (a) and T = 700 K - (b). The error shows the standard deviation of the instantaneous temperature values

В формуле (1) значения всех величин определяются из результатов моделирования, кроме параметра dQ, который необходимо задать. Выбор параметра dQ ограничен следующими условиями: во-первых - этот параметр не может быть слишком большим, т.к. в этом случае кристаллит перегреется и может измениться его структура; во-вторых - при слишком малом значении dQ мгновенные значения "горячего" и "холодного" термостатов будут перекрываться; в-третьих - надо учитывать, что поскольку к термостатам подводится и отводится одно и то же количество энергии равное dQ, температуры "горячего" и "холодного" термостатов должны отличаться от заданной температуры на одно и то же

значение (не более, чем величина статистической ошибки). Таким образом, от выбора dQ зависит, выйдет ли система на стационарный режим без нарушения структуры, после чего можно проводить процедуру усреднения для получения температурного профиля. При изменении значения параметра dQ необходимо повторять всю процедуру заново, и чем меньше заданная температура системы, тем сложнее подобрать подходящее значение этого параметра.

На рис. 2 показано как выбор параметра dQ влияет на значения температурного профиля в алюминиевом образце длиной 96 элементарных ячеек. При температуре Т = 70 К было выполнено несколько расчетов с разными значениями параметра dQ(eV): 2.2, 2.0, 1.8, 1.6, 1.4. Чтобы разница между результатами была лучше видна, на рисунке представлена только "линейная" область профиля, по которой затем вычисляется значение коэффициента теплопроводности. Отметим, что во всех случаях процедура релаксации, выход на стационарный режим и усреднение температуры выполнено в одних и тех же временных интервалах. Из рисунка видно, что при уменьшении параметра dQ изменяется наклон прямых, следовательно, и разность температур между крайними точками линейного интервала. Отметим, что при Т = 70 К величина ошибки ±0.9 К (см рис. 1, а). Из рис. 2 видно, что при задании dQ = 2.2 нарушается третье условие - равенства отклонений от заданной температуры, что видно даже на рис. 2. А при dQ = 1.4 происходит перекрывание мгновенных значений температур термостатов. Таким образом, результаты с параметрами dQ = 2.2 и 1.4 не удовлетворяют перечисленным условиям. Среднее значение по оставшимся трем реализациям кр1 = 24.386 (Ш/(тК)).

68.5 -1-'-

35 40 45 50 55 60

Layer Number

Рис. 2. Линейные области температурного профиля для кристаллита Al (12x12x96 u.c.) при T = 70 K с различными значениями параметра dQ

Fig. 2. Linear regions of the temperature profile for the Al crystallite (12*12x96 u.c.) at T =70 K with different values of the parameter dQ

В таблице показаны результаты расчета kpi при допустимых значениях параметра dQ, т.е. при выходе на стационарный режим теплового потока. Из физики задачи следует, что значение коэффициента теплопроводности для фиксированной температуры не должно зависеть от отношения dQ/A Т.

Таблица - Теплопроводность в зависимости от выбора параметра dQ

Table - Thermal conductivity depending on the choice of parameter dQ

dQ, meV AT и K dQ/AT,, eV/K kp,, W/(m K)

2.2 3.263 0.674 23.123

2.0 2.761 0.724 24.850

1.8 2.562 0.702 24.097

1.6 2.267 0.705 24.211

1.4 2.050 0.682 23.409

Таким образом, учитывая все перечисленные особенности расчетов, были получены температурные зависимости коэффициента решеточной теплопроводности алюминия для кристаллитов двух размеров вдоль направления [001]. Результаты приведены на рис. 3. Здесь линии с темными кружками соответствуют расчетам для кристаллита длиной 96 u.c., линии с темными квадратиками - 48 u.c. Как видно из рисунка, при T > 600 K эти линии практически сливаются, а при T > 700 K значения £Р;(Т) для кристаллитов разной длины совпадают. Это связано с особенностями фононного спектра рассчитанного в [9, 10] - минимальная частота фононов в длинном кристаллите в два раза меньше, чем в коротком, т.е. низкочастотная область спектра короткого образца обеднена. При низких температурах, когда основной вклад в коэффициент теплопроводности дают длинноволновые фононы, эта разница относительно велика, а при высоких становится пренебрежительно малой.

Для сравнения на этом же графике пунктирными линиями представлены расчеты решеточной теплопроводности с учетом только фонон-фононного рассеяния для алюминия, выполненные из первых принципов на основе теории функционала плотности (DFT) с использованием разных функций для обменно-корреляционного потенциала и различных подгоночных параметров. Линия 1, вычисленная в приближении локальной плотности (LDA) и линия 2 - в приближении обобщенного градиента (GGA), взяты из работы [11]. Линия 3 показывает результаты вычислений температурной зависимости kp(T) алюминия с помощью теории возмущения функционала плотности (DFPT) и динамики решетки в сочетании с решением транспортного уравнения Больцмана (BTE), в котором скорости вычисляются в приближении времени релаксации [12]. Отметим, что наши результаты и результаты, полученные из первых принципов, совпадают только при высоких значениях температуры, когда рассеяние оказывает заметное влияние и длина свободного пробега фононов уменьшается.

о -1-1-1-»-1-1-

0 100 200 300 400 500 600 700 800 900 Temperature, К

Рис. 3. Температурная зависимость решеточной теплопроводности Al для кристаллитов сечением 12x12 u.c. и длиной вдоль направления [001]: 96 u.c. - линии с кружками;

48 u.c. - линии с треугольниками. Пунктирными линиями - расчеты решеточной теплопроводности Al из первых принципов: 1 и 2 из работы [11], 3 - из [12]

Fig. 3. Temperature dependence of the lattice thermal conductivity of Al for crystallites with a cross section of 12x12 u.c. and a length along the direction [001]: 96 u.c. - lines with circles; 48 u.c. - lines with triangles.

Dotted lines - calculations of the lattice thermal conductivity of Al from the first principles: 1 and 2 - [11], 3 - [12]

Известно, что при повышении температуры фонон-фононное рассеяние становится определяющим фактором в теплопроводности, в результате которого коэффициент теплопроводности уменьшается обратно пропорционально температуре. Наши расчеты показывают такую же тенденцию и полученные значения kpt(T) находятся в разумных пределах в сравнении с первопринципными, а при высоких температурах совпадают.

Таким образом, применяя метод неравновесной молекулярной динамики в предложенной схеме расчета, т.е. учитывая только линейную область распределения температуры, были получены значения коэффициента решеточной теплопроводности для А1 в зависимости от температуры. Тем не менее, надо представлять все недостатки этого метода и учитывать их при расчетах теплопроводности металлов, находящихся в твердом состоянии:

1) сложности связанные с выбором корректного потенциала;

2) зависимость результатов расчета от размеров базового кристаллита;

3) сложности с выбором параметра dQ, особенно при низких температурах;

4) невозможность получения результатов при температурах ниже 50 К.

МЕТОД РАСЧЕТА РЕШЕТОЧНОЙ ТЕПЛОПРОВОДНОСТИ ПРИ НИЗКИХ ТЕМПЕРАТУРАХ

В случае низких температур, наиболее интересный для низкоразмерных систем в связи с условием их эксплуатации будем исходить из однородного уравнения теплопроводности, заданного на бесконечной прямой [13]

^ = ВАТ (z, t) (2)

ot

с начальными условиями (задача Коши): T(z,0) = ф(г). Здесь z - координата -<ж< z <<ж, D - коэффициент температуропроводности, связанный с коэффициентом теплопроводности (к) по формуле к = CpD; где C - теплоемкость, р - плотность материала. Фундаментальное решение f(z,t) уравнения (2) с начальными условиями ф(г) = S(z-z0), представляет собой функцию Гаусса с дисперсией о2 = 2D t и математическим ожиданием z0:

f ( z ) = —^exp С 2ж

f (z - Zo)n

V

2c2

(3)

^ 2 В левой части рис. 4 продемонстрировано поведение Т(£) для значения а = 1,

а в правой - некоторые траектории релаксации температуры в центре распределения (3) для

той же дисперсии и различных Д которые описываются уравнением:

Щ z = ю) = 0.5Т0 (2пЫуш. (4)

z t

Рис. 4. Слева - качественное поведение T(z,ffl = 1) в произвольных единицах. Справа - траектории релаксации температуры в центре распределения T(t,z = z0) при увеличении D (сверху вниз)

Fig. 4. On the left is the qualitative behavior of T(z,o2=1) in arbitrary units. On the right are the temperature relaxation trajectories in the center of the distribution T(t,z=z0) with increasing D (from top to bottom)

Из рис. 4 следует, что для любой дисперсии а существует множество траекторий релаксации. Абсциссы точки пересечения горизонтальной прямой, проходящей от максимума распределения, с этими кривыми определяют время старта релаксации ¿0 при некотором Б. Значение коэффициента Б для конкретной системы можно определить, аппроксимируя выражением (4) полученную траекторию релаксации.

Таким образом, для нахождения коэффициента температуропроводности в рамках метода молекулярной динамики необходимо смоделировать в системе стационарное неравновесное распределение температуры в виде кривой Гаусса (3), которое полностью определяется дисперсией, т.е. произведением D•t0, где ¿0 - это время соответствующее началу релаксации системы.

После достижения стационарного состояния необходимо дать системе свободно релаксировать к равновесному состоянию. Заметим, что для выполнения условия квазисвободной релаксации при периодических граничных условиях (по у и г) требуется достаточно большой размер базового кристаллита или малое время счета, за которое возмущение не распространяется до соседних кристаллитов, транслируемых граничными условиями. В случае выполнения этих требований, процесс релаксации "контролируется" только коэффициентом температуропроводности Б, что позволяет найти его значение независимо от ¿0.

В данной работе, предлагаемый метод был использован для нахождения коэффициента температуропроводности тонкой алюминиевой пленки при низких температурах. Параметры пленки были такие же, как в работах [9, 10], посвященных динамической устойчивости пленки с увеличенными размерами расчетной ячейки вдоль направлений у и г до 60 атомных слоев (238.3 А).

Для получения нужного распределения температуры по Гауссу вдоль оси г в центре расчетной ячейки выделим достаточно узкий слой ёг, состоящий из некоторого числа атомных слоев, перпендикулярных оси г. В начальный момент времени каждый атом в этом слое имеет кинетическую энергию Тн = Т0 + ёТ, а вне его - Т0. Это значит, что исходное распределение является прямоугольным (рис. 5). После чего ослабим требование к распределению энергии, считая, что Тн и Т0 относятся не к каждому атому, а к среднему значению температуры в слое ёг и вне его, соответственно. Это приведет к перестройке температурного поля до установления неравновесного стационарного состояния. Причем оказывается, что распределение температуры описывается выражением (3), то есть функцией Гаусса с дисперсией а 0 = 2 численное значение которой находится подгонкой.

Изменение распределения Т(г) при релаксации демонстрирует рис. 6.

T-original — 1— 1 1

"calc • —

gauss fr

■ M «t** •• , \ . . . •• . - .

----•----V-v- " "

• i i J_ i i

О 10 20 30 40 50 60

z,Layer Number (In)

Рис. 5. Схема построения распределения Гаусса. Т0 = 8 К, dT = 6 K, dz = 2 атомных слоя

Fig. 5. Scheme for constructing a Gaussian distribution. Т0 = 8 К, dT = 6 K, dz = 2 atomic layers

z , In

Рис. 6. Изменение распределения Гаусса в процессе релаксации, f - время релаксации (f = t -10)

Fig. 6. Change in the Gaussian distribution during relaxation, t is the relaxation time (f= t - t0)

Для проведения анализа удобней следить за температурой точки максимума распределения T(t,z = z0). Рис. 7 показывает поведение этой величины во время всего моделирования - от стабилизации начального условия (t < t0) до окончания процесса релаксации. Выражение (4) можно записать в виде удобном для нахождения D:

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

T(t,z = zo) = 0.5То (2nD(f+to ))-1/2 , (5)

х 2 где t = t - t0, из которого, уже зная о 0 = 2Dt0, можно найти t0 и D.

Рис. 7. Временная зависимость температуры центра распределения за все время счета. Пунктирная кривая является аппроксимацией кривой релаксации выражением (4) с D = 1.0-10-4 (m2/s)

Fig. 7. Time dependence of the temperature of the distribution center for the entire calculation time. The dotted curve is the approximation of the relaxation curve by expression (4) with D=1.0-10-4 (m2/s)

Полученное значение коэффициента решеточной температуропроводности D = 1.0 10-4 (m2/s) невозможно сравнить с экспериментальными результатами, в виду их отсутствия для пленок при низких температурах, но по порядку величины он совпадает со значениями коэффициента температуропроводности для объемного образца алюминия [14].

На рис. 8 показаны релаксационные кривые, полученные при различных начальных распределениях температуры. Из рисунка видно, что все они описываются одной

^ ^ 4 2 „

теоретической кривой со значением D = 1.0-10" (m/s). Таким образом, физическое свойство материала не зависит от выбора начальных условий расчета и демонстрирует разумность предложенного подхода для оценки решеточного вклада в теплопроводность.

Рис. 8. Релаксационные кривые, полученные для начальных условий: T0 = 8 K, 1) dz = 2 (ln), dT = 6 К; 2) dz = 2(ln); dT = 8 K; 3) dz = 1(ln), dT = 8 K. Пунктирная кривая - та же, что и на рис. 7. (Результаты расчетов в случаях 1 и 3 перекрываются)

Fig. 8. Relaxation curves obtained for the initial conditions: T0=8 K, 1) dz=2 (ln), dT=6 К; 2) dz=2(ln); dT=8 K; 3) dz=1(ln), dT=8 K. The dotted curve is the same as in Fig. 7. (The calculation results in cases 1 and 3 overlap)

ЗАКЛЮЧЕНИЕ

В работе показано, что для определения решеточного вклада в общую теплопроводность металлов и металлических сплавов возможно применение метода МД. С помощью этого метода, используя апробированный потенциал, проведены исследования решеточной теплопроводности алюминия при высоких и низких температурах.

Показано, что для вычисления коэффициента решеточной теплопроводности при высоких температурах целесообразно использовать метод неравновесной молекулярной динамики. В этих случаях исходным уравнением для нахождения коэффициента теплопроводности является уравнение Фурье и задачей МД расчетов становится получение стационарного неравновесного состояния с линейным градиентом температуры на длине, сравнимой с размерами расчетной ячейки. Проведенные в работе исследования показали, что такой подход приводит к получению значений коэффициента решеточной теплопроводности, согласующихся с результатами первопринципных расчетов. С уменьшением размера базового кристаллита происходит уменьшение коэффициента теплопроводности, которое обусловлено обеднением низкочастотного участка фононного спектра, вклад которого в теплопроводность с ростом температуры становится незначительным. При высоких температурах коэффициент теплопроводности не зависит от размера кристаллита и совпадает с полученным из первопринципных расчетов.

Для вычисления теплопроводности при низких температурах в работе предложен оригинальный метод, исходящий из уравнения теплопроводности, заданного на бесконечной прямой. В этом случае задачей моделирования является получение в системе стационарного неравновесного распределения температуры в виде кривой Гаусса, которое является фундаментальным решением этого уравнения. Аппроксимация релаксации системы из неравновесного состояния к равновесному позволяет найти коэффициент температуропроводности, связанный с коэффициентом теплопроводности.

Тестовые расчеты, проведенные для тонкой пленки алюминия при низких температурах с различными начальными условиями, показали, что полученный коэффициент температуропроводности не зависит от параметров начального распределения Гаусса, что указывает на пригодность предлагаемого метода для изучения решеточной теплопроводности.

Работа выполнена с использованием кластера "Уран" СКЦ ИММ УрО РАН в рамках темы НИР УдмФИЦ УрО РАН "Теоретические исследования фазовых состояний, спектральных и кинетических свойств электронов и фононов в системах с пониженной размерностью" № 121030100005-1.

The work was carried out using the Uranium cluster of the SCC IMM UB RAS within the framework of the R&D theme of the UdmFRC UB RAS "Theoretical studies of phase states, spectral and kinetic properties of electrons and phonons in systems with reduced dimension" No. 121030100005-1.

СПИСОК ЛИТЕРАТУРЫ

1. Bhalla P., Kumar P., Das N., Singh N. Theory of the dynamical thermal conductivity of metals // Physical Review B, 2016, vol. 94, iss. 11, article: 115114, 12 p. https://doi.org/10.1103/PhysRevB.94.115114

2. MuUer-Plathe F. A simple nonequilibrium molecular dynamics method for calculating the thermal conductivity // Journal of Chemical Physics, 1997, vol. 106, no. 14, pp. 6082-6085. https://doi.org/10.1063/L473271

3. Evans D. J., Hoover W. G. Flows Far From Equilibrium Via Molecular Dynamics // Annual Review of Fluid Mechanics, 1986, vol. 18, pp. 243-264. https://doi.org/10.1146/annurev.fl.18.010186.001331

4. Plimpton S. Fast Parallel Algorithms for Short-Range Molecular Dynamics // Journal of Computational Physics, 1995, vol. 117, no. 1, pp. 1-19. https://doi.org/10.1006/JCPH.1995.1039

5. Zope R. R., Mishin Y. Interatomic potentials for atomistic simulations of the Ti-Al system // Physical Review B, 2003, vol. 68, article: 024102, 14 p. http://dx.doi.org/10.1103/PhysRevB.68.024102

6. Dolgusheva E. B., Trubitsin V. Y. Lattice Heat Capacity of Nanostructured Materials Based on Titanium/Zirconium and Aluminum // Physics of the Solid State, 2018, vol. 60, pp. 837-846. https://doi.org/10.1134/S1063783418050074

7. Dolgusheva E. B. Thermal properties of fcc titanium and aluminum thin films // Computational Materials Science, 2018, vol. 155, pp. 55-62. https://doi.org/10.1016/j.commatsci.2018.08.033

8. Wirnsberger P., Frenkel D., Dellago C. An enhanced version of the heat exchange algorithm with excellent energy conservation properties // Journal of Chemical Physics, 2015, vol. 143, article: 124104, 8 p. https://doi.org/10.1063/1.4931597

9. Саламатов Е. И., Долгушева Е. Б. Атомарный механизм разрушения тонкой пленки при высоких температурах. Метод молекулярной динамики // Химическая физика и мезоскопия. 2020. Т. 22, № 3. C. 281-288. https://doi.org/10.15350/17270529.2020.3.27

10. Salamatov E. I., Dolgusheva E. B. The role of long-wave bending vibrations in the destruction of ultrathin Al films // Physica Status Solidi B, 2021, vol. 258, iss. 3, article: 2000484, 7 p. https://doi.org/10.1002/pssb.202000484

11. Wang Y., Lu Z., Ruan X. First principles calculation of lattice thermal conductivity of metals considering phonon-phonon and phonon-electron scattering // Journal of Applied Physics, 2016, vol. 119, iss. 22, article: 225109, 10 p. https://doi.org/10.1063/L4953366

12. Jain A., McGaughey A. J. H. Thermal transport by phonons and electrons in aluminum, silver, and gold from first principles // Physical Review B, 2016, vol. 93, article: 081206(R), 5 p. http://dx.doi.org/10.1103/PhysRevB.93.081206

13. Тихонов А. Н., Самарский А. А. Уравнения математической физики. М.: Наука, 1977. 736 с.

14. Зиновьев В. Е. Теплофизические свойства металлов при высоких температурах. Справочник. М.: Металлургия, 1989. 384 с.

REFERENCES

1. Bhalla P., Kumar P., Das N., Singh N. Theory of the dynamical thermal conductivity of metals. Physical Review B, 2016, vol. 94, iss. 11, article: 115114, 12 p. https://doi.org/10.1103/PhysRevB.94.115114

2. Muller-Plathe F. A simple nonequilibrium molecular dynamics method for calculating the thermal conductivity. Journal of Chemical Physics, 1997, vol. 106, no. 14, pp. 6082-6085. https://doi.org/10.1063/L473271

3. Evans D. J., Hoover W. G. Flows Far From Equilibrium Via Molecular Dynamics. Annual Review of Fluid Mechanics, 1986, vol. 18, pp. 243-264. https://doi.org/10.1146/annurev.fl.18.010186.001331

4. Plimpton S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. Journal of Computational Physics, 1995, vol. 117, no. 1, pp. 1-19. https://doi.org/10.1006/JCPH.1995.1039

5. Zope R. R., Mishin Y. Interatomic potentials for atomistic simulations of the Ti-Al system. Physical Review B, 2003, vol. 68, article: 024102, 14 p. http://dx.doi.org/10.1103/PhysRevB.68.024102

6. Dolgusheva E. B., Trubitsin V. Y. Lattice Heat Capacity of Nanostructured Materials Based on Titanium/Zirconium and Aluminum. Physics of the Solid State, 2018, vol. 60, no. 5, pp. 837-846. https://doi.org/10.1134/S1063783418050074

7. Dolgusheva E. B. Thermal properties of fcc titanium and aluminum thin films. Computational Materials Science, 2018, vol. 155, pp. 55-62. https://doi.org/10.1016/j.commatsci.2018.08.033

8. Wirnsberger P., Frenkel D., Dellago C. An enhanced version of the heat exchange algorithm with excellent energy conservation properties. Journal of Chemical Physics, 2015, vol. 143, article: 124104, 8 p. https://doi.org/10.1063/1.4931597

9. Salamatov E. I., Dolgusheva E. B. Atomarnyy mekhanizm razrusheniya tonkoy plenki pri vysokikh temperaturakh. Metod molekulyarnoy dinamiki [Atomic Mechanism of Thin Film Destruction at High Temperatures. Molecular Dynamics Method]. Khimicheskaya fizika i mezoskopiya [Chemical Physics and Mesoscopy], 2020, vol. 22, no. 3, pp. 281-288. (In Russian). https://doi.org/10.15350/17270529.2020.3.27

10. Salamatov E. I., Dolgusheva E. B. The role of long-wave bending vibrations in the destruction of ultrathin Al films. Physica Status Solidi B, 2021, vol. 258, iss. 3, article: 2000484, 7 p. https://doi.org/10.1002/pssb.202000484

11. Wang Y., Lu Z., Ruan X. First principles calculation of lattice thermal conductivity of metals considering phonon-phonon and phonon-electron scattering. Journal of Applied Physics, 2016, vol. 119, iss. 22, article: 225109, 10 p. https://doi.org/10.1063/L4953366

12. Jain A., McGaughey A. J. H. Thermal transport by phonons and electrons in aluminum, silver, and gold from first principles. Physical Review B, 2016, vol. 93, article: 081206(R), 5 p. http://dx.doi.org/10.1103/PhysRevB.93.081206

13. Tikhonov A. N., Samarskiy A. A. Uravneniya matematicheskoiy fiziki [Equations of mathematical physics]. Moscow: Nauka Publ., 1977. 736 c.

14. Zinovev V. E. Teplofizicheskie svoistva metallov pri vysokikh temperaturakh [Thermophysical properties of metals at high temperatures]. Spravochnik. Moscow: Metallurgiya Publ., 1989. 384 c.

Поступила 22.11.2021; после доработки 06.12.2021; принята к опубликованию 10.12.2021 Received 22 November 2021; received in revised form 06 December 2021; accepted 10 December 2021

Саламатов Евгений Иванович, доктор физико- Evgeny I. Salamatov, Dr. Sci. (Phys.-Math.), Chief

математических наук, главный научный сотрудник Researcher, Udmurt Federal Research Center UB RAS,

УдмФИЦ УрО РАН, Ижевск, Российская Федерация Izhevsk, Russian Federation

Долгушева Елена Борисовна, кандидат физико- Elena B. Dolgusheva, Cand. Sci. (Phys.-Math.), Senior

математических наук, старший научный сотрудник Researcher, Udmurt Federal Research Center UB RAS,

УдмФИЦ УрО РАН, Ижевск, Российская Федерация, Izhevsk, Russian Federation, e-mail: elena@udman.ru e-mail: elena@udman.ru

i Надоели баннеры? Вы всегда можете отключить рекламу.