CC BY
64
УДК 536.2.02
РАСЧЕТ МАКРОХАРАКТЕРИСТИК НАНОСИСТЕМ.
ЧАСТЬ 1. КОЭФФИЦИЕНТ ТЕПЛОПРОВОДНОСТИ ОДНОРОДНЫХ НАНОСИСТЕМ
ВАХРУШЕВ А. В., СЕВЕРЮХИН А. В., СЕВЕРЮХИНА О. Ю.
Институт механики Уральского отделения РАН, 426067, г. Ижевск, ул. Т. Барамзиной, 34
АННОТАЦИЯ. В данной работе приводятся физические основы, а также численные методики расчета коэффициента теплопроводности однородных наносистем. Рассмотрены уравнения, описывающие многочастичные потенциалы МЕАМ, ЕБ1Р и др. Определены температурные зависимости коэффициента теплопроводности для различных видов материалов. Выполнены расчеты теплофизических характеристик однородных наносистем на основе кремния и золота. Представлены кривые температурной зависимости коэффициента теплопроводности для систем различной размерности. Проведено сравнение данных, полученных с использованием потенциалов МЕАМ и ЕБ1Р с экспериментальными данными. Выявлено, что характер кривых и значения, полученные при моделировании, хорошо согласуются с экспериментальными данными.
КЛЮЧЕВЫЕ СЛОВА: наноматериалы, теплопроводность, моделирование, молекулярная динамика.
ВВЕДЕНИЕ
В конце 20-го столетия сформировалась новая междисциплинарная область науки и техники, получившая название "Нанотехнологии", связанная с изучением и управлением групп атомов и молекул, характерные размеры которых (или размер по какому-либо одному или двум направлениям) не превышают 100 нм (1 нм = 10-9 м). При этом, как показывают экспериментальные исследования, механические, электромагнитные, химические и другие свойства наноэлементов в несколько раз отличаются от свойств материалов, в которых атомы и молекулы организованы обычным, давно известным из классической физики, образом. Это требует глубокого понимания протекающих в наномасштабе физико-химических процессов [1], однако малый масштаб процессов затрудняет их исследование экспериментальными методами и приводит к необходимости использования математического моделирования.
Математическое моделирование является мощным инструментом при проектировании различного типа наносистем и анализе протекающих в них процессов [2 - 6]. Особенно важно использовать методы математического моделирования в новых, "пионерских" областях науки и техники, в которых ещё не накоплен опыт эксплуатации и экспериментальные данные. Главными задачами математического моделирования в наносистемах являются: формирование наноэлементов (форма, структура, свойства), взаимодействие отдельных элементов наносистемы (формирование и нагружение), определение структуры изолированной наносистемы в квазистатическом и динамическом состояниях, расчет параметров наносистемы при взаимодействии с окружающей средой и расчет макропараметров наносистемы. Отметим, что последняя в списке задач, а именно расчет макропараметров наносистем в зависимости от их размера, структуры и химического состава, является наиболее актуальной в настоящий момент. Это обусловлено совершающимся сейчас быстрым переходом от исследования фундаментальных свойств наносистем к решению задач их получения и применения, в которых макропараметры наносистем имеют определяющее значение.
Данной работой мы начинаем последовательное изложение теоретических основ, методов моделирования и результатов расчетов макрохарактеристик наносистем, основанных на работах по моделированию процессов формирования и структуры различных наносистем. Будут рассмотрены последовательно теплофизические параметры, упругие
и прочностные характеристики, реологические параметры, параметры адсорбции и другие макрохарактеристики наносистем.
В представленной статье, в развитие предыдущих работ авторов [7 - 8], рассмотрена задача расчета коэффициента теплопроводности однородных наносистем.
ТЕОРЕТИЧЕСКИЕ ОСНОВЫ
Вопрос переноса тепла в твердых телах широко исследуется еще с XVIII века. Так в 1807 году французский физик и математик Жан Батист Жозеф Фурье представил свою работу «О распространении тепла в твёрдом теле» [9]. Согласно закону Фурье, тепловой поток может быть определен исходя из формулы:
q = -kVT. (1)
В уравнении (1) символом q обозначен вектор теплового потока, характеризующий поток тепла, проходящий через единицу площади в единицу времени, k - коэффициент
теплопроводности, а под VT понимается градиент температуры. Несмотря на то, что с его открытия прошло более 200 лет, закон Фурье и по сей день обеспечивает основу для моделирования процессов переноса тепла в твердых телах.
Физика теплообмена претерпела значительную эволюцию в течение прошлых двух столетий. В начале XIX века, считалось, что объект, называемый «калория» переносит тепло от горячего тела к холодному [10]. Классическая термодинамика была разработана главным образом на основе этого понимания. Однако в 1914 году голландский физик Петр Йозеф Вильгельм Дебай предложил модель, в соответствии с которой передача тепла происходит за счет распространения колебательной энергии. Это была первая корреляция теплообмена с кинематикой фундаментальных частиц. В модели Дебая колебания решетки рассматриваются как газ фононов - квазичастиц, теоретически обоснованных советским физиком Игорем Евгеньевичем Таммом в 1932 году. Такой вид процесса передачи тепла принято называть решеточной или фононной теплопроводностью.
Перенос тепла может быть осуществлён посредством двух механизмов: посредством свободных электронов и за счёт тепловыми колебаниями атомов (фононная теплопроводность). За последнее столетие, строгий эксперимент показал, что тепло передается посредством стохастического движения элементарных частиц, таких как электроны, атомы и молекулы.
В твердых телах, как отмечалось ранее, тепловые колебания атомов обуславливают данный вид теплопроводности. С увеличением температуры, как известно, возрастает и внутренняя энергия кристаллической решетки материала и, как следствие, растет и амплитуда колебаний атомов. Поскольку связь между атомами решетки достаточно велика, то в случае возникновения теплового возбуждения в одном месте решетки, оно будет передаваться в виде некой упругой волны между ними. Упругая волна, достигнув поверхности тела, отражается от неё. Затем происходит образование стоячей волны, как следствие наложения двух волн (прямой и отраженной). Такая стоячая волна называется нормальным колебанием и имеет некоторую частоту (. Если рассматривать систему, содержащую N атомов, то в такой системе возможно возбудить 3N нормальных колебаний. Соответственно частоты образующихся колебаний можно обозначить (, где ie1...3N.
Необходимо учитывать, что на поверхности должны находиться узловые точки, в которых равняется нулю амплитуда колебаний стоячей волны. Причем расстояние между узлами должно быть таковым, чтоб в нем помещалось целое число длин полуволн (см. рис. 1).
Тогда максимальная частота колебаний будет рассчитываться по формуле:
w = 2pv // @pv Id,
max sei mm зв I >
где V3e - скорость звука, а расстояние между атомами - d .
Минимальная частота колебаний будет определяться исходя из линейного размера кристалла - L :
О = 2pu ll =pu IL.
mm зв I max зв I
L^b_;_J
a) 6)
Рис. 1. Модель твердого тела одномерная (а) — цепочка атомов длины L , б) — стоячие волны с длинной волны l:
1 - l = Imax = 2L , 2 - l = L, 3 - l = 2L/3, 4 - l = lmm = 2d
Тогда максимальная частота колебаний будет рассчитываться по формуле:
О = 2pu ll Id,
max зв I mm зв I >
где изв - скорость звука, а расстояние между атомами - d .
Минимальная частота колебаний будет определяться исходя из линейного размера кристалла - L :
О = 2pu ll =pu IL.
mm зв I max зв I
Исходя из квантования энергии нормальных колебаний и пропорциональности смещений атомов силам межатомного взаимодействия (гармоническое приближение), и как следствие отсутствие взаимодействия нормальных колебаний между собой, их энергия может быть записана в виде:
Ег = h( ^n + , n = 0,1,2,...,
где h = -h--приведенная постоянная Планка, ( - частота г -го нормального колебания.
Тогда под фононом будет пониматься минимальный квант энергии нормально колебания атомов решетки.
Теплопроводность диэлектриков (решеточная теплопроводность)
Объяснить температурную зависимость и конечное значение теплопроводности решетки можно, только отказавшись от представления фононного газа идеальным.
Если рассматривать идеальный фононный газ, то даже при малейшем изменении температуры должен возникать поток фононов, переносящих тепловую энергию. Рассмотрим уравнение закона Фурье для одномерного случая - распространения тепла вдоль некоторой оси. Он может быть записан следующим образом:
q = ~kdTi, (2)
dx
где dT представляет собой градиент температуры вдоль выбранной оси. dx
Из формулы (2) видно, что для того, чтобы величина теплового потока при бесконечно малом изменении температуры (dT ® 0) была конечна, необходимо, чтобы значение коэффициента теплопроводности было бесконечно большим. Однако, коэффициент теплопроводности в реальных телах имеет конечное значение. Причиной этого служит взаимодействие между фононами (рассеяние фононов на фононах), в результате которого
часть энергии фононов передается решетке, что обеспечивает установление теплового равновесия и ограничивает величину теплопроводности решетки.
Если применить молекулярно-кинетическую теорию газов к фононному газу, то коэффициент теплопроводности для фононного газа будет записываться в виде:
к » 3 РСЛ1) и,
где р - плотность газа, Су - его удельная теплоемкость, т.е. теплоемкость единицы объема фононного газа, (/) - средняя длина свободного пробега фононов, Ц/ - средняя скорость
фононов, которая будет равна скорости звука в твердом теле.
В соответствии с законом Дюлонга-Пти теплоемкость твердых тел не зависит от температуры, однако изменение концентрации фононов будет пропорционально температуре
в диапазоне температур Т » 0О (где 0О = - характеристическая температура Дебая,
кв
кв - константа Больцмана). Таким образом, коэффициент теплопроводности при температуре Т » 0О будет зависеть только от длины свободного пробега фононов и, как
следствие, будет обратно пропорционален температуре и концентрации фононов, т.е. к ~ — .
Это согласуется с экспериментальными данными. На рис. 2 приведен график зависимости коэффициента теплопроводности монокристалла искусственного сапфира (Л120з) от температуры.
Рис. 2. График зависимости коэффициента теплопроводности монокристалла искусственного сапфира (А1203) от температуры
В диапазоне же температур Т « 0О/20 длина свободного пробега фононов не будет
зависеть от температуры, а будет обуславливаться непосредственно размерами рассматриваемого тела. Таким образом, температурная зависимость коэффициента теплопроводности сведется к к ~ Т3, что качественно согласуется с экспериментальными данными.
Теплопроводность металлов
В отличие от диэлектриков в металлах работают оба механизма переноса тепла: как засечёт электронов, так и посредством фононного газа. Таким образом, коэффициент теплопроводности для данного класса веществ может быть записан в виде:
к = ке + крк.
В чистых металлах процесс переноса тепла осуществляется в основном посредством свободных электронов, концентрация которых на единицу объема в металлах весьма велика. Для определения качественной зависимости электронной составляющей коэффициента
теплопроводности от температуры можно представить всю совокупность электронов металла в виде электронного газа, по аналогии с фононным газом описанным выше. Тогда воспользовавшись молекулярно-кинетической теорией газов получим:
ке = 3 СеУР/е
где под Се следует понимать удельную теплоемкость электронного газа, равную ж2Тпк
Се =-в ,где кв - константа Больцмана, п - концентрация свободных электронов;
2ЕР
под ьр - скорость движения электронного газа, соответствующую энергии Ферми Ер, и определяемую как у 2Ер/т ; под /е - длину свободного пробега электронов. Таким образом, получим:
ж2пк„Т/„
к=
"В^ У
ЗтЬр
При температуре Т» 0О длина свободного пробега электронов будет определяться исходя из рассеяния электронов на фононах. Причем концентрация фононов будет прямо пропорциональна температуре. Тогда длина свободного пробега электронов будет обратно пропорциональна температуре (/е ~ 1/Т), и коэффициент теплопроводности не будет
зависеть от температуры. Это положение хорошо согласуется с данными экспериментов представленными на рис. 3.
Что касается низкотемпературного диапазона, там наблюдается иная картина: в соответствии с экспериментальными данными (/е ~ 1/Т3). Исходя из этого, получим
ке ~ 1/Т2. При температурах близких к абсолютному нулю (Т ® 0) концентрация фононов будет мала, вследствие чего, коэффициент теплопроводности будет прямо пропорционален температуре (ке ~ Т).
Следует отметить, что в чистых металлах при нормальных условиях теплопроводность электронного газа много больше решеточной теплопроводности.
МЕТОДИКА МОДЕЛИРОВАНИЯ
В качестве метода моделирования использовался аппарат молекулярной динамики. Метод МД получил широкое распространение при моделировании поведения наносистем благодаря простоте реализации и высокой точности моделирования.
В молекулярно-динамических расчетах величину коэффициента теплопроводности можно вычислить различными способами [11].
В данной работе используется формализм Грина-Кубо (Огееи-КиЬо), который связывает автокорреляционную функцию теплового потока с коэффициентом теплопроводности. Тепловой поток может быть рассчитан из колебаний потенциальной и кинетической энергии атома и тензора напряжения атома в стационарном уравновешенном моделировании. Это отличие от двух предыдущих неравновесных методов молекулярной динамики (КЕМО), где энергия течет непрерывно между горячими и холодными областями в моделируемом образце.
Коэффициент теплопроводности в модели Грина-Кубо рассчитывается по следующей формуле [12, 13]:
к = Ь Й ТГ^ ^ 3 (')3 (0)> ^ (3)
где к - коэффициент теплопроводности ё -мерной системы с линейным размером Ь, Т - температура, кв - постоянная Больцмана, 3 - компонента потока тепла.
Постоянная времени Т - это минимальное время, необходимое для того, чтобы автокорреляционная функция теплового потока затухала до нуля. Этот метод более подходит для анизотропных систем. Основным недостатком этого метода является то, что для затухания до нуля автокорреляционной функции теплового потока требуется много времени, а получаемые значения теплопроводности зависят от размера системы.
Автокорреляционные функции справа в формуле (3) оцениваются в равновесии, без градиента температуры. Автокорреляция — статистическая взаимосвязь между случайными величинами из одного ряда, но взятыми со сдвигом, например, для случайного процесса — со сдвигом по времени. Автокорреляционная функция может быть определена как:
ВД = \ / ($)/($ -Т)Л.
Общий тепловой поток в системе вычисляется как
3 (0 = | ](х, г)ск,
где ] (х, () - плотность теплового потока.
Порядок пределов в формуле (3) имеет большое значение. При правильных порядках пределов можно вычислить корреляционные функции с произвольными граничными условиями и применять формулу (3). Существуют различные формы записи уравнения (1) [12 - 18].
Существуют ситуации, когда формула (3) не применима. Во-первых, для маленьких систем, которые изучаются в мезоскопической физике, термодинамический предел не имеет смысла. Во-вторых, во многих низко-размерных системах перенос тепла аномальный и теплопроводность существенно отклоняется от экспериментальных значений [19 - 20]. В таких случаях невозможно взять пределы в формуле (3). Поэтому тепловую проводимость рассматривают как функцию от длины Ь. В литературе, посвященной этой тематике [19 - 20] обычно в формуле (3) меняют верхний предел интегрирования 1с на Ь. Другой способ использования формулы Огееи-КиЬо для конечных систем заключается в том, чтобы внедрять бесконечные резервуары, как это сделано в работах [21 - 22].
В данной работе для расчета теплопроводности кремниевых наноструктур был выбран метод Огееи-КиЬо. Для вычисления коэффициента теплопроводности использовалась следующая запись формулы (3) [23]:
1 ¥ 1 1- К 3, (0) 3х (* ))*=-1
к = щТ1 Ц 3 (0) 3х (')^ =:МТ* Ц1 (0)'1 (')^, (4)
где V - объем системы.
Уравнения, описывающие рассматриваемую атомарную структуру, представляют собой систему дифференциальных уравнений, которая определяет движение всех атомов системы [24]:
— 2г. г
щ-^ = К (г,Г(0),* = 1,2,...,N, (5)
^ = 0, Г (^)= Го, = и (') = и0, * = 1,2,..,N, (6)
где тг - масса г -го атома, N - число атомов в системе, г0 - начальный радиус-вектор г -го атома, Г - текущий радиус-вектор г -го атома, К (г, Г (г)) - суммарная сила, действующая на г -й атом. Выражение (6) задает начальные условия для рассматриваемой системы, где Ц и Ц0 - текущая и начальная скорость г -го атома соответственно.
Силы К (г, г (г)) в уравнении (5), определяются из соотношения:
п. , ч. ди (г (г))
К (г,Г(г)) =--^^,г = 1,2,...,N , (7)
,\> \ П д г4 (г) ' ' ' ' ''
где Г (г) = {г15 Г2,..., ^ } , и(г (г)) - некоторая потенциальная функция, описывающая
взаимодействие всех атомов системы.
Для решения рассматриваемой задачи потенциал взаимодействия системы может быть записан в виде:
и (Г (г )) = Е,
где Е - потенциал взаимодействия.
Очень важную роль в аппарате молекулярной динамики является выбор потенциала взаимодействия. Особенно эта проблема стоит остро для металлов и полупроводников. Кремний относится к материалам с ковалентными типами связей. В зависимости от температуры он обладает различной структурой и свойствами. Так при нормальных условиях для него характерна алмазная структура. С ростом давления в нем могут формироваться новые кристаллографические структуры: простая кубическая, гранецентрированная кубическая. В данном случае наблюдается качественный переход от одной структуры к другой, происходит увеличение координационного числа. Наличие описанных особенностей делает задачу выбора и построения межатомного потенциала чрезвычайно сложной.
В пакете ЬЛММРБ для кремния доступно множество потенциалов: Ьеппагё-1опе8, МЕЛМ_1^, МЕЛМ_2^, Б1, Те^ой", ВОР, МЕЛМ_8РЬШЕ_1, МЕЛМ_8РЬШЕ_2, МЕЛМ_81, ЕБ1Р.
При учете только парного межатомного взаимодействия в математическом моделировании металлических и/или полупроводниковых систем возникает ряд проблем. В работе [25] было показано, что при использовании только парного потенциала взаимодействия в системах металл и/или полупроводник выполняется нефизическое соотношение для коэффициентов Коши (С12 = С44).
Парные потенциалы не могут обеспечить реалистичных значений физических характеристик материала [26]. Для корректного описания свойств твердых тел необходимо использовать многочастичные потенциалы. Известно, что ни один из существующих потенциалов не способен воспроизвести полный набор характеристик твердых веществ. Таким образом, выбор потенциала для математического моделирования - сложная комплексная задача. Большинство эмпирических потенциалов хорошо описывают объемные свойства материалов, но, тем не менее, некоторые с успехом используются для описания и поверхностных свойств.
При моделировании металлических и полупроводниковых систем наибольшее распространение получили следующие подходы, учитывающие многочастичное взаимодействие:
• потенциал Стиллинжера-Вебера [27];
• потенциал Абеля-Терсоффа [28];
• метод погруженного атома [25, 29];
• модифицированный метод погруженного атома [30].
В методе погруженного атома (ЕАМ) энергия связи системы атомов представлена в виде:
E = Z F
Z р j (R)
+ -2
1Z jj (R), (8)
где Z Fi
Z p j (R)
j
- функция погружения атома i, зависящая от суммарной электронной
плотности в области погружения i-го атома; jij- (R) - энергия парного взаимодействия.
Вывод данного уравнения с использованием теории функционала электронной плотности (DFT) можно найти в работе [29].
Каждый атом системы рассматривается как частица, погруженная в электронный газ, создаваемый остальными атомами моделируемой системы. Энергия, необходимая для погружения, зависит от электронной плотности в точке погружения. Введенная таким образом функция погружения позволяет определить обменную и корреляционную энергии электронного газа системы.
Смысл функции погружения может быть определен как энергия, необходимая для погружения одного атома в однородный электронный газ с плотностью p. Однако существуют и другие преобразования [31], позволяющие изменять функции (8) с тем условием, что результирующие энергия и межатомные силы не изменятся. В ЕАМ используются следующие приближения:
1) Функция электронной плотности одного атома является сферически симметричной функцией, зависящей только от расстояния между атомами. Данное приближение существенно ограничивает область применения ЕАМ и позволяет рассматривать системы, в которых направленностью ковалентной составляющей связи можно пренебречь.
2) Электронная плотность в области погружения атома i определяется как линейная
суперпозиция электронных плотностей остальных атомов системы Z Р j (R). Данное
j *i
приближение существенно упрощает вычисление электронной плотности.
3) Величина Z Р j (R) в металлических системах в области расположения атома i
j*i
меняется слабо по сравнению с электронной плотностью самого атома рг.. Таким образом,
Z Р j (R) в области расположения атома i заменяется константой p [29]. Энергия
j *i
электронного газа аппроксимируется функцией, зависящей только от величины среднего значения электронной плотности в области погружения, а не сложными функционалами, как в методе DFT.
В настоящее время выведены потенциалы EAM для большинства металлов и некоторых бинарных систем. Также рассчитаны потенциалы для тройных систем [32]. Однако такие «тройные» потенциалы качественно не воспроизводят физические свойства материалов.
На основе описанных методик был предложен полуэмпирический подход, объединяющий преимущества многочастичных потенциалов и метода погруженного атома. Теория модифицированного метода погруженного атома (MEAM - modified embedded-atom method) выведена с применением теории функционала электронной плотности (DFT) [33]. Метод DFT в настоящее время считается наиболее признанным подходом к описанию электронных свойств твердых тел. В методе EAM полная электронная плотность
представляется в виде линейной суперпозиции сферически усредненных функций. Данный недостаток устранен в модифицированном методе погруженного атома.
В методе МЕЛМ энергия связи системы записывается в следующем виде:
( ( р Л
Е = I
К
V V ^г У
р' + ¿Еф (Я,)
где Е - энергия атома г; Кг - функция погружения для атома г, погруженного в электронную плотность рг; 2{ - число ближайших соседей атома г в его референтной кристаллической структуре; ф, - парный потенциал между атомами г,,, находящимися на расстоянии Я,.
В МЕЛМ функция погружения К (р) определяется как
К (р) = АЕср 1пр , где А - регулируемый параметр; Ес - энергия связи.
Парный потенциал между атомами г,, определяется по формуле
ф, (Я ) = | \ Е (Я)- К
( 7^0
р0{Я)
V ^ У
где рг - электронная плотность.
Полная электронная плотность в точке погружения включает в себя угловые зависимости и записывается в виде:
р = р(0)о (г).
Существует много разновидностей для функции О (Г). Однако, наибольшее распространение получила форма записи в виде:
О (Г) = л/Т+Г. Функция Г вычисляется по формуле:
3 ,,( р(й) Л2
г = I г(И) I р
И=1
р
(0)
где И = 0 - 3, соответствуют я, р, -, / симметрии; г(И) - весовые множители; р(И) - величины, определяющие отклонение распределения электронной плотности от распределения в идеальном кристалле кубической сингонии р(0) :
(И = 0): р(0)= I ра(0)(г г)
Р(И = 1): (р(1))2 = I
I ра(1)( Г )■"
-(И = 2): (р(2))2 = I
а,Ь
г'г1
I р"р>( Г)Га Г
■11
1
I р* 2)( Г)
/(И = 3): (р<3')2 = I
а,р,у
I ра(3)( Г)
Га Г Г
Здесь ра(И) - радиальные функции, которые представляют уменьшение вклада расстояний Г, верхний индекс г указывает ближайшие атомы, а, Ь, у - индексы суммирования по каждому из трех возможных направлений. Наконец, индивидуальный вклад вычисляется по формуле:
ра( И)( г ) = рс
-р(И) IГ-1
я
2
а
2
2
Г
2
Г
е
Потенциал EDIP (Environment-Dependent Interatomic Potential) был впервые предложен в работе [34]. По сути, он представляет собой объединение формализмов Стиллинджера-Вебера и Терсоффа и включает двухчастичные и трехчастичные взаимодействия:
1 ( ^ E = 2ZIZV (,z) + Z Z V (r,rik,z)
2 i V j*i j*i k*i,k> j у
V (r, Z ) = A
— | - e r
bz2
exp
s
r - a
V3 (rj, Гк, zi) = g (rij) g (Гк)h (cos 0yk, zi),
z = Z f (Rim ),
f (r ) =
1, r < c
exp
a
1 - X~
0, r > a
,c < r < a,
(1 -e-e(z)(l+t(z))2 ) + (z)(i + T(z))
-mz
к (1,7 ) = 1
й (7 ) = йе X (7) = и1 + и2 (ы3в~"47 - в~2"47), где Е - полная энергия системы, приходящаяся на атом, 0г]к - угол между гранями у и гк, Гу - единичный вектор, направленный от атома I к атому у, а - радиус обрезания, где 7г - координационное число, К, (г, 7) - функция двухчастичного взаимодействия, У3 (Гу, гк, 7г) - функция трехчастичного взаимодействия, / ( г ), g (г ), к (I,7) - функция
обрезания, радиальная и угловая функции соответственно.
Как показано в работе [35], данный вид потенциала лучше описывает аморфный кремний, чем потенциал Стиллинджера-Вебера.
Потенциал ЕБ1Р имеет также ряд преимуществ: хорошо воспроизводит эластические характеристики кремния и предсказывает температуру плавления близкую к экспериментальным данным.
РЕЗУЛЬТАТЫ РАСЧЕТОВ
В ходе моделирования были рассмотрены системы различной размерности от 4*4*4 до 4x4x144 (в элементарных ячейках). Проводились расчеты для полупроводниковых систем, на основе монокристаллов кремния и металлических систем на основе чистого золота. В ходе моделирования были использованы такие виды потенциалов взаимодействия, как МЕАМ и ЕБГР. Несмотря на то, что теплофизические характеристики кремниевых систем является предметом интереса многих ученых, данные полученные различными группами исследователей порой противоречивы. Так в работе [36] коэффициент теплопроводности равен 20 Вт/(м-К) для системы 4*4*144 элементарные ячейки. В работе [37] коэффициент теплопроводности равен 235 Вт/(м-К). Квантовые расчеты дают значение 357 Вт/(м-К). В работе [38] коэффициент теплопроводности (9*9*302, потенциал Б') при температуре 1000 К равен 20 Вт/(м- К). В работе [39] коэффициент теплопроводности для системы 7*7*7 при 600 К получился равным 37 Вт/(м-К), 10*10*10 - 43 Вт/(м-К). В то же время, эксперимент показывает значение 64 Вт/(м-К).
Для проверки корректности работы, предложенной модели были построены графики автокорреляционной функции теплового потока, а значения функции в формуле Грина-Кубо
были проверены на сходимость. В результате данных операций был сделан вывод о корректности предложенной модели и возможности ее использования для расчета теплофизических характеристик наносистем.
На рис. 4 представлена температурная зависимость коэффициента теплопроводности кремния для систем разной размерности.
юоо
0 100 200 300 400 500 600 700 800 900 1000
Т,к
О Е01Р 4x4x4 -М-Е01Р 4x4x40 -£—Е01Р 4x4x144 -+- МЕАМ 4x4x4 -В—[13] -«-[14] Рис. 4. Температурная зависимость коэффициента теплопроводности кремния, полученная при
моделировании и из экспериментов
Как видно из рисунка характер кривых, полученных в ходе моделирования, соответствует экспериментальным данным. Наиболее четкое соответствие получается при температуре в области температуры Дебая и выше неё. Это позволяет говорить о возможности использования предложенной модели для расчета теплофизических характеристик наносистем, а также прогнозирования значений рассчитываемых параметров для макросистем.
Следует отметить, что при низких температурах размер системы влияет на значение коэффициента теплопроводности, что может быть обусловлено использованием метода Грина-Кубо, однако с ростом температуры влияние размерного фактора нивелируется.
Использование различных типов потенциалов также сказывается на значении коэффициента теплопроводности. Это подтверждается графиками на рис. 5.
500
О 100 200 300 400 500 600 700 800 900 1000
Т,к
О ЕС11Р 4x4x4 -Ч-МЕАМ 4x4x4
Рис. 5. Сравнение значений коэффициента теплопроводности системы 4x4x4, полученных с использованием потенциалов МЕАМ и ЕБ1Р
Исходя из представленных данных, можно сделать вывод, что при температурах выше температуры Дебая (645±5 К для кремния) различные потенциалы дают схожие значения коэффициента теплопроводности. Однако, в области низких температур выбор потенциала оказывает значительное влияние на получаемых значения, данной теплофизической характеристики.
По результатам моделирования процесса теплопроводности в металлических системах, в частности в кристаллах золота, были построены аналогичные зависимости. Кривая температурной зависимости золота представлена на рис. 6.
1
0,9 0.8 0:7 0.6
а 0,5 J 0,4 0,3 0,2 0,1
0 -t-,-,-,-,-,-,-,-,-,-,-,-,-,-,-,-,-,-,-,-
0 100 200 300 400 500 600 700 800 900 1000
Т,к
-•-4x4x4 -и-~1ЯЛ2
Рис. 6. Температурная зависимость коэффициента теплопроводности золота
Поведение кривой соответствует описанным выше закономерностям для металлов, что свидетельствует о хорошей согласованности с экспериментальными данными и теоретическими выкладками.
ВЫВОДЫ
В данной работе приведены методики расчета теплофизической характеристики вещества, в частности коэффициента теплопроводности, для различных веществ. Описаны теоретические температурные зависимости коэффициента теплопроводности металлов и диэлектриков. Представлена математическая модель процесса теплопроводности в наносистемах. Изложены методики расчета характеристик данного процесса. Также в работе показаны особенности и различия потенциалов взаимодействия, используемых в молекулярно-динамических расчетах.
В ходе моделирования были рассчитаны коэффициенты теплопроводности систем с использованием таких материалов как кремний и золото. Построенные температурные зависимости коэффициентов теплопроводности кремния показали хорошую согласованность с экспериментальными данным и работами других исследователей. Характер кривой зависимости коэффициента теплопроводности золота от температуры соответствует представленным теоретическим выкладкам для металлов, что в свою очередь свидетельствует о возможности применения представленных методик моделирования для прогнозирования теплофизических характеристик различных веществ.
СПИСОК ЛИТЕРАТУРЫ
1. Суздалев И. П. Нанотехнология: физико-химия нанокластеров, наноструктур и наноматериалов. М.: КомКнига, 2006. 592 с.
2. Steinhauser M. O. Computational Multiscale Modeling of Fluids and Solids. Theory and Application. BerlinHeidelberg: Springer-Verlag, 2008. 427 p.
3. Вахрушев А. В., Липанов А. М., Суетин М. В. Моделирование процессов аккумуляции водорода и углеводородов наноструктурами. М.-Ижевск: Институт компьютерных исследований; НИЦ «Регулярная и хаотическая динамика», 2008. 118 с.
4. Аликин В. Н., Вахрушев А. В., Голубчиков В. Б., Липанов А. М., Серебренников С. Ю. Разработка и исслдование аэрозольных нанотехнологий. Топлива. Заряды. Двигатели. Т. 3. М.: Машиностроение, 2010. 196 с.
5. Suyetin M. V., Vakhrushev A. V. Nanocapsule for safe and effective methane storage // Nanoscale Research Letters, 2009, vol. 4, no. 11, pp. 1267-1270.
6. Вахрушев А. В., Федотов А. Ю. Моделирование формирования композиционных наночастиц из газовой фазы // Международный научный журнал Альтернативная энергетика и экология. 2007. № 10. С. 22-26.
7. Вахрушев А. В., Северюхин А. В., Северюхина О. Ю., Федотов А. Ю. Исследование теплофизических свойств наноматериалов на основе кремния методом Green-Kubo с использованием потенциала EDIP // Химическая физика и мезоскопия. 2016. Т. 18, № 2. С. 187-198.
8. Северюхин А. В., Северюхина О. Ю., Вахрушев А. В., Федотов А. Ю. Исследование теплофизических свойств кремниевых наноматериалов методом Green-Kubo // Труды Института механики УрО РАН «Проблемы механики и материаловедения». Ижевск: ИМ УрО РАН, 2016. С. 210-223.
9. Fourier J. Théorie de la Propagation de la Chaleur Dans Les Solides. Manuscrit présenté à l'Institut de France,
1807.
10. Chen G. Nanoscale Energy Transport and Conversion: A Parallel Treatment of Electrons, Molecules, Phonons, and Photons. Oxford University Press, USA, 2005. 560 p.
11. URL: http://lammps.sandia.goV/doc/Section_howto.html#calculating-thermal-conductivity (дата обращения 6 апреля 2016).
12. Green M. S. Markoff Random Processes and the Statistical Mechanics of Time-Dependent Phenomena. II. Irreversible Processes in Fluids // Journal of Chemical Physics, 1954, vol. 22, pp. 398-413.
13. Kubo R., Yokota M. and Nakajima S. Statistical-Mechanical Theory of Irreversible Processes. II. Response to Thermal Disturbance // Journal of the Physical Society of Japan, 1957, vol. 12, pp. 1203-1211.
14. Mori H. Statistical-Mechanical Theory of Transport in Fluids // Physical Review, 1958, vol. 112, iss. 6, pp. 1829-1842.
15. Green M. S. Comment on a Paper of Mori on Time-Correlation Expressions for Transport Properties // Physical Review, 1960, vol. 119, iss. 3, pp. 829-830.
16. Kadanoff L.P. and Martin P.C. Hydrodynamic equations and correlation functions // Annals of Physics, 1963, vol. 24, pp. 419-469.
17. Luttinger J. M. Theory of Thermal Transport Coefficients // Physical Review, 1964, vol. 135, pp. A1505-
A1514.
18. Visscher W. M. Transport processes in solids and linear-response theory // Physical Review A, 1974, vol. 10, pp. 2461-2472.
19. Lepri S., Livi R. and Politi A. Thermal conduction in classical low-dimensional lattices // Physics Reports, 2003, vol. 377, pp. 1-80.
20. Kundu A., Dhar A. and Narayan O. The Green-Kubo formula for heat conduction in open systems // Journal of Statistical Mechanics: Theory and Experiment, 2009, iss. 3, pp. L03001.
21. Allen K. R. and Ford J. Lattice Thermal Conductivity for a One-Dimensional, Harmonic, Isotopically Disordered Crystal // Physical Review, 1968, vol. 176, pp. 1046-1055.
22. Fisher D. S. and Lee P. A. Relation between conductivity and transmission matrix // Physical Review B, 1981, vol. 23, pp. 6851-6854.
23. URL: http://lammps.sandia.gov/doc/compute_heat_flux.html (дата обращения 06.04. 2016).
24. Vakhrushev A. V., Severyukhina O. Yu., Severyukhin A. V., Vakhrushev A. A., Galkin N. G. Simulation of the processes of formation of quantum dots on the basis of the transition metals // Nanoscience and Technology: An International Journal, 2012, vol. 3, iss. 4. pp. 51-75.
25. Daw M. S., Baskes M. I. Embedded-atom method: derivation and application to impurities, surfaces, and other defects in metals // Physical Review B, 1984, vol. 29, no. 12, pp. 6443-6453.
26. URL: http://www.fisica.uniud.it/~ercolessi/ (дата обращения 06.04. 2016).
27. Stillinger F. H., Weber T. A. Computer simulation of local order in condensed phases of silicon // Physical Review B, 1985, vol. 31, pp. 5262-5271.
28. Tersoff J. New empirical approach for the structure and energy of covalent systems // Physical Review B, 1988, vol. 37, no. 12, pp. 6991-7000.
29. Daw M. S. Model of metallic cohesion: The embedded-atom method // Physical Review B, 1989, vol. 39, no. 11, pp. 7441-7452.
30. Baskes M. I. Modified embedded-atom potentials for cubic materials and impurities // Physical Review B, 1992, vol. 46, pp. 2727-2742.
31. Ruda M., Farkas D., Abriata J. Interatomic potentials for carbon interstitials in metals and intermetallics // Scripta Materialia, 2002, vol. 46, no. 5, pp. 349-355.
32. Tomar V., Zhou M. Classical molecular-dynamics potential for the mechanical strength of nanocrystalline composite fcc Al+ a -Fe2O3 // Physical Review B, 2006, vol. 73, no. 17, 174116.
33. Hohenberg P., Kohn W. Inhomogeneous electron gas // Physical Review, 1964, vol. 136, no.3B, pp. 864-871.
34. Justo J. F., Bazant M. Z. Interatomic potential for silicon defects and disordered phases // Physical Review B, 1998, vol. 58, no. 5, pp. 2539-2550.
35. Allred C. L., Yuan X., Bazant M. Z. and Hobbs L. W. Elastic constants of defected and amorphous silicon with the environment-dependent interatomic potential // Physical Review B, 2004, vol. 70, pp. 134113.
36. Schelling P. K., Phillpot S. R. and Keblinski P. Comparison of atomic-level simulation methods for computing thermal conductivity // Physical Review B, 2002, vol. 65, pp. 144306(12).
37. Lee Y., Lee S. and Hwang G. S. Effects of vacancy defects on thermal conductivity in crystalline silicon: A nonequilibrium molecular dynamics study // Physical Review B, 2011, vol. 83, pp. 125202(7).
38. Hu M., Zhang X., Giapis K. P. and Poulikakos D. Thermal conductivity reduction in core-shell nanowires // Physical Review B, 2011, vol. 84, pp. 085442(9).
39. Esfarjani K. and Chen G. Heat transport in silicon from first-principles calculations // Physical Review B, 2011, vol. 84, pp. 085204(11).
CALCULATION OF MACRO CHARACTERISTICS OF NANOSYSTEMS.
PART 1. COEFFICIENT OF THERMAL CONDUCTIVITY OF HOMOGENEOUS NANOSYSTEMS
Vakhrouchev A. V., Severyukhin A. V., Severyukhina O. Yu.
Institute of Mechanics, Ural Branch of the Russian Academy of Science, Izhevsk, Russia
SUMMARY. This paper provides a physical basis and numerical methods of calculation of thermal conductivity of homogeneous nanosystems. Equations describing multiparticle potentials MEAM, EDIP, etc. are considered. The features and differences of the interaction potentials used in molecular dynamics calculations are shown in the paper. The temperature dependences of the heat conductivity coefficient for various types of materials are determined. Calculations of thermophysical characteristics of homogeneous nanosystems based on silicon and gold are performed. In this work the formalism of Green-Kubo, which connects an autocorrelation function of a heat flux with a thermal conductivity, is used. The simulation was performed using the software package LAMMPS. The curves of the temperature dependence of the thermal conductivity coefficient for systems of various dimensions are presented. Comparison of the data, obtained with use of capacities of MEAM and EDIP, with experimental data is carried out. It is revealed that the shape of the curves and the values obtained in the simulation are in good agreement with experimental data. The nature of the curve of temperature dependence of coefficient of thermal conductivity of gold corresponds to the submitted theoretical calculations for metals. This indicates the possibility of using the presented modeling techniques for predicting the thermophysical characteristics of various substances.
KEYWORDS: nanomaterials, thermal conductivity, simulation, molecular dynamics.
REFERENCES
1. Suzdalev I. P. Nanotekhnologiya: fiziko-khimiya nanoklasterov, nanostruktur i nanomaterialov [Nanotech-nology: physical chemistry of nanoclusters, nanostructures and nanomaterials]. Moscow: KomKniga Publ., 2006. 592 p.
2. Steinhauser M. O. Computational Multiscale Modeling of Fluids and Solids. Theory and Application. BerlinHeidelberg: Springer-Verlag, 2008. 427 p. doi: 10.1007/978-3-662-53224-9
3. Vakhrushev A. V., Lipanov A. M., Suetin M. V. Modelirovanie protsessov akkumulyatsii vodoroda i uglevodorodov nanostrukturami [Modeling of processes of accumulation of hydrogen and hydrocarbons by nanostructures]. Moscow-Izhevsk: Institut komp'yuternykh issledovaniy; NITs Regulyarnaya i khaoticheskaya dinamika Publ., 2008. 118 p.
4. Alikin V. N., Vakhrushev A. V., Golubchikov V. B., Lipanov A. M., Serebrennikov S. Yu. Razrabotka i issledovanie aerozol'nykh nanotekhnologiy. Topliva. Zaryady. Dvigateli. Tom 3 [Development and research of aerosol nanotechnologies. Fuel. Charges. Engines. Vol. 3]. Moscow: Mashinostroenie Publ., 2010. 196 p.
5. Suyetin M. V., Vakhrushev A. V. Nanocapsule for safe and effective methane storage. Nanoscale Research Letters, 2009, vol. 4, no. 11, pp. 1267-1270. doi: 10.1007/s11671-009-9391-x
6. Vakhrushev A. V., Fedotov A. Yu. Modelirovanie formirovaniya kompozitsionnykh nanochastits iz gazovoy fazy [Modelling of composite nanoparticle formation from a gas phase]. Mezhdunarodnyy nauchnyy zhurnal Al'ternativnaya energetika i ekologiya [International Scientific Journal for Alternative Energy and Ecology], 2007, no. 10, pp. 22-26.
7. Vakhrushev A. V., Severyukhin A. V., Severyukhina O. Yu., Fedotov A. Yu. Issledovanie teplofizicheskikh svoystv nanomaterialov na osnove kremniya metodom Green-Kubo s ispol'zovaniem potentsiala EDIP [Research of thermal properties of nanomaterials on the silicon basis by the Green-Kubo method using the EDIP potential]. Khimicheskaya fizika i mezoskopiya [Chemical Physics and Mesoscopy], 2016, vol. 18, no. 2, pp. 187-198.
8. Severyukhin A. V., Severyukhina O. Yu., Vakhrushev A. V., Fedotov A. Yu. Issledovanie teplofizicheskikh svoystv kremnievykh nanomaterialov metodom Green-Kubo [Study of thermophysical properties of silicon nanomaterials by the Green-Kubo method V sbornike Trudy Instituta mekhaniki UrO RAN. Problemy mekhaniki i materialovedeniya [Problems of Mechanics and Materials Science, Proc. IM, UB RAS. Izhevsk: IM UrO RAN Publ., 2016, pp. 10-223.
9. Fourier J. Théorie de la Propagation de la Chaleur Dans Les Solides. Manuscrit présenté à l'Institut de France, 1807.
10. Chen G. Nanoscale Energy Transport and Conversion: A Parallel Treatment of Electrons, Molecules, Phonons, and Photons. Oxford University Press, USA, 2005. 560 p.
11. http://lammps.sandia.gov/doc/Section_howto.html#calculating-thermal-conductivity (accessed April 6, 2016).
12. Green M.S. Markoff Random Processes and the Statistical Mechanics of Time-Dependent Phenomena. II. Irreversible Processes in Fluids. Journal of Chemical Physics, 1954, vol. 22, pp. 398-413. doi: 10.1063/1.1740082
13. Kubo R., Yokota M. and Nakajima S. Statistical-Mechanical Theory of Irreversible Processes. II. Response to Thermal Disturbance. Journal ofthe Physical Society ofJapan, 1957, vol. 12, pp. 1203-1211.
14. Mori H. Statistical-Mechanical Theory of Transport in Fluids. Physical Review, 1958, vol. 112, iss. 6, pp. 1829-1842. doi: 10.1103/PhysRev.112.1829
15. Green M. S. Comment on a Paper of Mori on Time-Correlation Expressions for Transport Properties // Physical Review, 1960, vol. 119, iss. 3, pp. 829-830. do: 10.1103/PhysRev.119.827
16. Kadanoff L.P. and Martin P.C. Hydrodynamic equations and correlation functions. Annals of Physics, 1963, vol. 24, pp. 419-469. doi: 10.1016/0003-4916(63)90078-2
17. Luttinger J. M. Theory of Thermal Transport Coefficients. Physical Review, 1964, vol. 135, pp. A1505-A1514. doi: 10.1103/PhysRev.135.A1505
18. Visscher W. M. Transport processes in solids and linear-response theory. Physical Review A, 1974, vol. 10, pp. 2461-2472. doi: 10.1103/PhysRevA.10.2461
19. Lepri S., Livi R. and Politi A. Thermal conduction in classical low-dimensional lattices. Physics Reports, 2003, vol. 377, pp. 1-80. doi: 10.1016/S0370-1573(02)00558-6
20. Kundu A., Dhar A. and Narayan O. The Green-Kubo formula for heat conduction in open systems. Journal of Statistical Mechanics: Theory and Experiment, 2009, iss. 3, pp. L03001. doi: 10.1088/1742-5468/2009/03/L03001
21. Allen K. R. and Ford J. Lattice Thermal Conductivity for a One-Dimensional, Harmonic, Isotopically Disordered Crystal. Physical Review, 1968, vol. 176, pp. 1046-1055. doi: 10.1103/PhysRev.176.1046
22. Fisher D. S. and Lee P. A. Relation between conductivity and transmission matrix. Physical Review B, 1981, vol. 23, pp. 6851-6854. doi: 10.1103/PhysRevB.23.6851
23. http://lammps.sandia.gov/doc/compute_heat_flux.html (accessed April 6, 2016).
24. Vakhrushev A. V., Severyukhina O. Yu., Severyukhin A. V., Vakhrushev A. A., Galkin N. G. Simulation of the processes of formation of quantum dots on the basis of the transition metals. Nanoscience and Technology: An International Journal, 2012, vol. 3, iss. 4, pp. 51-75. doi: 10.1615/NanomechanicsSciTechnolIntJ.v3.i1.30
25. Daw M. S., Baskes M. I. Embedded-atom method: derivation and application to impurities, surfaces, and other defects in metals. Physical Review B, 1984, vol. 29, no. 12, pp. 6443-6453. doi: 10.1103/PhysRevB.29.6443
26. http://www.fisica.uniud.it/~ercolessi/ (accessed April 6, 2016).
27. Stillinger F. H., Weber T. A. Computer simulation of local order in condensed phases of silicon. Physical Review B, 1985, vol. 31, pp. 5262-5271. doi: 10.1103/PhysRevB.31.5262
28. Tersoff J. New empirical approach for the structure and energy of covalent systems. Physical Review B, 1988, vol. 37, no. 12, pp. 6991-7000. doi: 10.1103/PhysRevB.37.6991
29. Daw M. S. Model of metallic cohesion: The embedded-atom method. Physical Review B, 1989, vol. 39, no. 11, pp. 7441-7452. doi: 10.1103/PhysRevB.39.7441
30. Baskes M. I. Modified embedded-atom potentials for cubic materials and impurities. Physical Review B, 1992, vol. 46, pp. 2727-2742. doi.org/10.1103/PhysRevB.46.2727
31. Ruda M., Farkas D., Abriata J. Interatomic potentials for carbon interstitials in metals and intermetallics. Scripta Materialia, 2002, vol. 46, no. 5, pp. 349-355. doi: 10.1016/S1359-6462(01)01250-7
32. Tomar V., Zhou M. Classical molecular-dynamics potential for the mechanical strength of nanocrystalline composite fcc Al+a-Fe2O3. Physical Review B, 2006, vol. 73, no. 17, 174116. doi: 10.1103/PhysRevB.73.174116
33. Hohenberg P., Kohn W. Inhomogeneous electron gas. Physical Review, 1964, vol. 136, no. 3B, pp. 864-871. doi: 10.1103/PhysRev.136.B864
34. Justo J. F., Bazant M. Z. Interatomic potential for silicon defects and disordered phases. Physical Review B, 1998, vol. 58, no. 5, pp. 2539-2550. doi: 10.1103/PhysRevB.58.2539
35. Allred C. L., Yuan X., Bazant M. Z. and Hobbs L. W. Elastic constants of defected and amorphous silicon with the environment-dependent interatomic potential // Physical Review B, 2004, vol. 70, pp. 134113. doi: 10.1103/PhysRevB.70.134113
36. Schelling P. K., Phillpot S. R. and Keblinski P. Comparison of atomic-level simulation methods for computing thermal conductivity. Physical Review B, 2002, vol. 65, pp. 144306(12). doi: 10.1103/PhysRevB.65.144306
37. Lee Y., Lee S. and Hwang G. S. Effects of vacancy defects on thermal conductivity in crystalline silicon: A nonequilibrium molecular dynamics study. Physical Review B, 2011, vol. 83, pp. 125202(7). doi: 10.1103/PhysRevB.83.125202
38. Hu M., Zhang X., Giapis K. P. and Poulikakos D. Thermal conductivity reduction in core-shell nanowires. Physical Review B, 2011, vol. 84, pp. 085442(9). doi: 10.1103/PhysRevB.84.085442
39. Esfarjani K. and Chen G. Heat transport in silicon from first-principles calculations. Physical Review B, 2011, vol. 84, pp. 085204(11). doi: 10.1103/PhysRevB.84.085204
Вахрушев Александр Васильевич, доктор физико-математических наук, профессор, заведующий лабораторией механики наноструктур ИМ УрО РАН, тел. (3412) 21-45-83, e-mail: vakhrushev-a@yandex. ru
Северюхин Александр Валерьевич, кандидат физико-математических наук, ученый секретарь ИМ УрО РАН, тел. (3412) 50-88-10, e-mail: severfam@mail. ru
Северюхина Олеся Юрьевна, кандидат физико-математических наук, научный сотрудник ИМ УрО РАН, e-mail: lesienok@mail. ru
