CC BY
49
УДК 539.3, 538.9
МОЛЕКУЛЯРНО-ДИНАМИЧЕСКИЙ РАСЧЕТ КОЭФФИЦИЕНТА ТЕПЛОВОГО РАСШИРЕНИЯ НАНОКЛАСТЕРОВ МЕДИ
Е. И. Головнева1, И. Ф. Головнев1, В. М. Фомин1
1 Учреждение Российской академии наук Институт теоретической и прикладной механики им. С. А. Христиановича Сибирского отделения РАН
elena@itam.nsc.ru, golovnev@itam.nsc.ru, fomin@itam.nsc.ru PACS 31.15.xv, 65.80.-g, 82.60.Qr
В работе проведено молекулярно-динамическое исследование коэффициента теплового расширения (ктр) на-нокластеров меди сферической и кубической формы в широком диапазоне размеров. Для этого был проведен нагрев нанокластера методом стохастических сил. При этом через определенное число шагов по времени проводился расчет ряда характеристик, в том числе найдены линейные размеры исследуемой наноструктуры, ее объем и температура. Далее были найдены значения ктр, построены его зависимости от размера и формы наностуктуры. Показано, что коэффициент теплового расширения зависит от размера наноструктуры и с увеличением размера приближается к экспериментальному значению для макрообразцов.
Ключевые слова: молекулярно-динамическое моделирование, медный нанокластер, коэффициент теплового расширения.
1. Введение
Коэффициент теплового расширения является одной из важнейших термодинамических характеристик материалов и до недавнего времени его находили исключительно из экспериментальных исследований. С появлением высокопроизводительных (мощных) электронно-вычислительных систем получил развитие метод молекулярной динамики, который, при наличии потенциалов межатомного взаимодействия для вещества, позволяет численно рассчитывать и эту характеристику материалов. До настоящего времени коэффициент теплового расширения находился методом молекулярной динамики для систем с использованием периодических граничных условий (см., например, [1-3]). А это означает, что моделировалась макросистема и найденный коэффициент не описывает свойства наноструктур. В то же время, авторы работы [4] на основе анализа экспериментальных данных сделали предположение об определяющей роли различия коэффициентов теплового расширения твердой и жидкой фаз металла при кристаллизации. Ими был предложен следующий механизм кристаллизации. Кристаллические кластеры с размерами порядка нанометров, являющиеся зародышами кристаллизации [5], имеют коэффициенты теплового расширения примерно на порядок больше аналогичной характеристики жидкой фазы. При понижении температуры расплава кристаллические кластеры уменьшают свои размеры больше, чем окружающая жидкая фаза. Это обусловливает возникновение градиента напряжений в объеме, изменяющего процесс диффузии атомов из жидкой фазы к кристаллическому кластеру, что в свою очередь значительно влияет на сам процесс кристаллизации.
Поставленная проблема определила необходимость исследования зависимости коэффициента теплового расширения металлических нанокластеров от размеров и формы. В
работе проведены расчеты для сферических и кубических нанокластеров меди с размерами до 9 нм.
2. Физическая система
Для исследования влияния размера и формы наноструктур на термодинамические свойства проведен молекулярно-динамический расчет коэффициента теплового расширения кластеров меди сферической и кубической формы, имеющих ГЦК-структуру. Выбрана
о
ориентация кристаллов [100]. Размеры радиусов наносфер варьируются от 20 до 90 A, а
о
размеры ребер кубов —от 40 до 120 A. Как известно, свойства металлов можно описать только с помощью многочастичных потенциалов межатомного взаимодействия. В связи с этим Доу и Баскесом [6,7] был предложен многочастичный потенциал взаимодействия атомов в металлах с гранецентрированной кубической решеткой, который рассчитывается в приближении так называемого «метода внедренного атома» (МВА). В данном случае выражение для потенциальной энергии имеет следующую форму:
i i j=i
где pi — электронная плотность в месте расположения атома с номером г, вызванная присутствием остальных атомов в системе. Для ее расчета использовались усредненные по углу электронные плотности свободных атомов (г), вычисленные при помощи теории Хартри—Фока. Таким образом электронная плотность будет иметь вид:
р^ = & (ra)-
j = i
Парные потенциалы описывают отталкивание атомов:
(г)(г)
ФаЬ (Г) = -,
Г
а эффективный заряд атома Za (г) монотонно уменьшается с увеличением расстояния за счет экранировки. Зависимость F (р) вытекает из условия удовлетворения системы универсальному уравнению состояния, предложенному в [8], согласно которому связь между когезионной энергией и объемом ячейки Вигнера-Зейца для широкого класса веществ описывается универсальной зависимостью с привлечением всего двух параметров.
Развитием идей Доу и Баскеса явилась разработка значительного количества потенциальных функций типа МВА для целого ряда ГЦК металлов и их сплавов, что дало возможность исследовать их объемные, поверхностные свойства, такие как поверхностная релаксация решетки, свойства жидкого состояния [9]. Авторы [10] придали МВА потенциалу твердый теоретический базис.
Однако, для более точного описания энергетики дефектов упаковки в кристалле, потенциальные функции должны учитывать взаимодействие атомов в первых трех координационных сферах. Это связано, например, с тем, что атомы дефекта упаковки в ГЦК кристалле локально находятся в структуре гексагональной плотноупакованной (ГПУ) решетки. Отличие же функций радиального распределения ГЦК и ГПУ кристаллов начинается только с третьей координационной сферы. Из этих соображений в работе использовался потенциал взаимодействия Воутера [11], который учитывает взаимодействие атомов в первых трех координационных сферах. Функции f (г), ф(г) и F(р), рассчитанные для кристаллов Ni, Pd, Pt, Cu, Ag, Al, даны авторами в табличном виде. Для их аналитического представления в настоящей работе использовались локальные сплайны (В-сплайны) третьей степени. При таком представлении итоговый потенциал имеет непрерывную вторую производную, что,
безусловно, является необходимым при использовании численной схемы второго порядка точности по временному шагу.
2.1. Охлаждение системы
Все расчеты проведены методом молекулярной динамики с использованием скоростной модификации Верле [12]. Шаг по времени при этом был равен т = 10_16с, а полный процесс охлаждения длился 250000 шагов по времени. Далее, находится минимум потенциальной энергии системы с помощью метода искусственной вязкости [13]. Таким образом, находилось состояние системы при температуре, близкой к абсолютному нулю, что соответствует статистически равновесному состоянию кристалла.
Суть метода искусственной вязкости состоит в следующем. На первом этапе моделирования выбираются форма и размер наноструктуры, которые определяются поставленной задачей. Однако, расположение атомов или молекул в такой системе соответствуют идеальной кристаллической структуре для которой и находился потенциал межатомного взаимодействия. В связи с этим, вся система не находится в минимуме потенциальной энергии из-за наличия развитой поверхности. Возникает задача нахождения расположения атомов в минимуме потенциальной энергии. Традиционные численные методы здесь не применимы, т.к. мы имеем системы включающие десятки и сотни тысяч атомов или, то же самое, аргументов в потенциальной функции системы и (г! ,г2 ,■■■■). В связи с этим используется на практике динамический поиск минимума потенциальной энергии. Он заключается в том, что системе предоставляется возможность эволюционировать в соответствии с динамическими законами механики, но при этом включается диссипативный процесс, который приводит к уменьшению первоначального избытка потенциальной энергии системы. В этом случае на атомы все время действуют диссипативные силы ^ = , где & - импульсы атомов, а и — коэффициент искусственной вязкости. Вариация этой величины позволяет ускорить или замедлить процесс эволюции системы в состояние с минимумом потенциальной энергии. Температура определялась через хаотическую составляющую кинетической энергии атомов системы.
Зависимость температуры системы в процессе охлаждения показана на рис. 1, а на рис. 2 — зависимость изменения потенциальной энергии системы. Полученные координаты и импульсы атомов использовались далее в качестве начальных данных.
2.2. Разогрев системы
После охлаждения системы, полученные координаты и импульсы атомов использовались в качестве начальных данных для разогрева наноструктуры с помощью метода стохастических сил [14]. В рамках этого метода решается динамическая задача движения атомов под действием стохастических сил, определяемых с помощью датчика случайных чисел. Параметрами действующей случайной силы являются: средняя частота ударов ш вдоль каждой декартовой оси и амплитуда /о, обеспечивающая прирост компоненты импульса атома при одном ударе Ар = /0 т.
Так как координаты атомов при ударе не меняются, скорость разогрева в данном случае определяется приростом кинетической энергии:
АЕ _ . (р? + др)2 (р°Л _ ш
^ = £ ^^^--Щ- = +(АЛ
а=1 г=1 \ / а=1 г=1
с последующим переносом энергии из кинетической в потенциальную составляющую. Первое слагаемое в скобках в правой части этого выражения при суммировании обращается
РИС. 1. Зависимость температуры от числа шагов по времени
РИС. 2. Зависимость изменения потенциальной энергии системы от числа шагов по времени
в нуль, так как случайные силы и импульсы не коррелированны, и в итоге имеем:
АЕ 3иМ 2
Таким образом, скорость разогрева не зависит от энергии кластера, в этом и заключается главное преимущество данного способа разогрева по сравнению с широко используемым масштабированием импульсов. Параметры случайной силы подбирались в соответствии с требованием равновесного состояния кластера в каждый момент процесса разогрева.
При этом через каждые 25 К разогрев прекращается, и система релаксирует в течение 10"12 с к термодинамическому равновесному состоянию, которое определяется по А-критерию Колмогорова [15]. В определенные моменты времени для всей совокупности атомов находился параметр
А= у/Ы тах \Р* (р)-^(р) |
где Р* (р) — статистическая функция распределения значений импульсов, а Р (р) — соответствующая теоретическая функция распределения. С помощью параметра А определялась вероятность реализации гипотезы о нормальном распределении
Рх = 1 - ^ (-1)к ехр(—2к2А2).
к=—<х
Численные исследования показали, что для значений Р\ ^ 0.8 распределение импульсов описывается распределением Максвелла.
После того, как система отрелаксировала к равновесному состоянию, все необходимые характеристики усреднялись по тепловым флуктуациям на интервале 10—12 с. Таким образом, была найдена зависимость радиусов, длин ребер и объемов кластеров от температуры.
2.3. Расчет размеров системы
Для нахождения коэффициента теплового расширения необходимо было определить линейные размеры кластеров. Он проводился по следующей схеме. Центр масс как кубического, так и сферического нанокластеров помещался в начало координат.
В нанокластерах кубической формы грани располагались перпендикулярно осям координат. Поэтому средняя координата где ^ = х, у, г, каждой грани куба вычислялась как среднее значение координат атомов ^ соответствующей грани. Линейный размер кластера вдоль той или иной оси координат находился как разница между средними координатами внешних атомных плоскостей, перпендикулярных выбранной оси.
Нанокластеры сферической формы в действительности не являются идеальной сферой, а представляют собой многогранник, каждая грань которого — внешняя атомная плоскость. Поэтому в качестве радиуса сферы использовалось расстояние от центра масс сферы до самого удаленного атома на внешней грани.
3. Результаты расчетов
Полученная в результате расчетов зависимость относительных изменений радиусов АД = (К — К0)/К0 от температуры для разных начальных размеров кластеров приведена на рис. 3. Необходимо отметить, что в случае сферических кластеров наблюдается силь-
о
ная зависимость тангенса угла наклона от начального размера структур вплоть до 80 А. При этом угол наклона уменьшается с увеличением радиуса. Так, для сфер с радиусами
о
80 и 90 А зависимости относительного изменения радиуса от температуры практически совпадают. Таким образом, можно сделать вывод, что при дальнейшем увеличении радиусов, эта величина выйдет на асимптоту, соответствующую аналогичной характеристике для макротел.
На следующем этапе расчетов была проведена линейная интерполяция по температуре для каждого фиксированного начального значения линейного размера кластеров. Для нахождения линейного коэффициента теплового расширения было использовано известное выражение А К = (К — Ко )/К0 = &В.АТ. С помощью метода наименьших квадратов были построены соответствующие интерполяционные полиномы, и была найдена зависимость линейного коэффициента теплового расширения от начального размера наноструктур.
Коэффициенты теплового расширения для нанокластеров кубической формы находились аналогично. На рис. 4 приведена зависимость относительного изменения длины ребра куба от температуры для разных размеров нанокластеров. Для иллюстрации на эту зависимость нанесена интерполяционная прямая. По рисунку видно, что она хорошо передает зависимость относительного изменения длины ребра куба для всего исследуемого интервала размеров.
Окончательные результаты расчетов КТР для сферических и кубических кластеров представлены на рис. 5. Для удобства на этом же рисунке приведен интервал (1.65 ^ 1.75 • 10—5К—1) экспериментальных значений коэффициента теплового расширения меди, найденных для макроскопических образцов.
Как видно, размеры и форма кластера оказывают сильное влияние на процесс теплового расширения. Для сферических кластеров наблюдается уменьшение коэффициента теплового расширения с увеличением размеров. При этом наблюдается приближение расчетного значения линейного коэффициента теплового расширения к экспериментальным значениям.
Аналогичным образом проведены молекулярно-динамические расчеты объемного коэффициента теплового расширения для нанокластеров сферической и кубической формы.
Рис. 3. Зависимость относительного изменения радиуса нанокластера от температуры для различных начальных радиусов. Сплошная
о
линия —R=20 А, штри-
о
ховая линия — R=40 А, штрих-пунктирная
о
линия — R=80 А, пунк-
о
тирная линия — R=90 А
РИС. 4. Зависимость относительного изменения длины ребра куба от температуры для кубов разных размеров. Сплошная
о
линия — L=40 А, штрих-пунктирная линия —
о
L=80 А, штриховая
о
линия —L=100 А, пунк-
о
тирная линия — L=120 А
РИС. 5. Зависимость линейного коэффициента теплового расширения от размера кластера (Ь — длина ребра куба, И — диаметр сферы). Сплошная линия—линейный коэффициент теплового расширения кубического кластера от длины ребра куба. Штриховая линия — линейный коэффициент теплового расширения сферического кластера от его диаметра. Штрих-пунктирная линия — интервал экспериментальных значений медных макро-образцов
Из физики твердого тела известно, что линейный коэффициент теплового расширения примерно в три раза меньше объемного коэффициента теплового расширения. Полученные результаты молекулярно-динамических расчетов отражают подобное соотношение для всех исследуемых в данной работе кластеров.
4. Заключение
На основе полученных результатов сделан ряд выводов:
• Как линейный, так и объемный коэффициенты теплового расширения металлических нанокластеров сферической формы зависят от размера исследуемой системы.
• С увеличением размера сферического нанокластера наблюдается приближение значений коэффициента теплового расширения к экспериментальным значениям, полученным для макротел (в данном случае - медь).
• Для нанокластера кубической формы не наблюдается зависимость от размера в пределах погрешности расчетов. Полученные значения как линейного, так и объемного коэффициента теплового расширения для кубов, близки к экспериментальному значению, полученному для макротел.
• Значение ктр для наноразмерных объектов зависит от формы этих объектов.
• По предложенной методике могут быть рассчитаны коэффициенты теплового расширения для любых металлов, для которых имеется хорошо апробированный потенциал межатомного взаимодействия.
Работа выполнена при поддержке интеграционного проекта ОЭММПУ РАН № 13 «Трибологические и прочностные свойства структурированных материалов и поверхностных слоев».
Литература
[1] Sekkal W., Bouhafs B., Aourag H., Certier M. Molecular-dynamics simulation of structural and thermodynamic properties of boron nitrid // J. Phys. Condens. Matter. — 1998. — V. 10. — P. 4975-4984.
[2] Moon Won Ha, Hwang Ho Jung. Structural and thermodynamical properties of GaN: a molecular dynamics simulation // Physics Letters A. — 2003. — V. 315. — P. 319-324.
[3] Berroukche A., Soudini B., Amara K. Molecular dynamics simulation study of structural, elastic and thermodynamic properties of tin below 286 K // Int. J. Nanoelectronics and Materials. — 2008. — V. 1. — P. 41-51.
[4] Губернаторов В.В., Сычева Т.С., Романов Е.П., Владимиров Л.Р. Роль теплового расширения фаз в процессе кристаллизации и рекристаллизации металлов // Докл. акад. наук. — 2007. — T. 413, № 1. — C. 41-44.
[5] Гельд П.В., Сидоренко Ф.А. Силициды переходных металлов четвертого периода. — М.: Металлургия, 1971. — 582 с.
[6] Foiles S.M. Calculation of the surface segregation of Ni-Cu alloys with the use of the embedded-atom method // Phys. Rev. B.. — 1985. — V. 32. — P. 7685-7693.
[7] Foiles S.M., Baskes M.I., Daw M.S. Embedded-atom method functions for the fcc metals Cu, Ag, Au, Ni, Pd, Pt, and their alloys // Phys. Rev. B. — 1986. — V. 33. — P. 7983-7991.
[8] Rose J.H., Smith J.R., Guinea F., Ferrante J. Universal features of the equation of state of metals // Phys. Rev. B. — 1984. — V. 29. — P. 2963-2969.
[9] Garrison B.J., Srivastava D. Potential energy surfaces for chemical reactions at solid surfaces // Annu. Rev. Phys. Chem. — 1995. — V. 46. — P. 373-394.
[10] Jacobsen K.W., Norskov J.K., Puska M.J. Interatomic interactions in the effective-medium theory // Phys. Rev. B. — 1987. — V. 35. — P. 7423-7442.
[11] Voter A.F. Embedded Atom Method Potentials for Seven FCC Metals: Ni, Pd, Pt, Cu, Ag, and Al. Los Alamos Unclassified Technical Report // LA-UR-93-3901, 1993. — 9 p.
[12] Allen M.P., Tildesley D.J. Computer simulation of liquids. — Oxford, N.Y.: Clarendon Press, 1987. — 385 p.
[13] Golovnev I.F., Golovneva E.I., Fomin V.M. Simulation of quasi-static processes in the crystals by molecular dynamics method // Physical mesomehanics. — 2003. — V. 6, № 5-6. — P. 41-45.
[14] Bolesta A.V., Golovnev I.F., Fomin V.M. Contact melting of nickel cluster at collision with rigid wall // Physical mesomehanics. — 2001. — V. 4, № 1. — P. 5-10.
[15] Колмогоров А.Н. Теория вероятностей и математическая статистика. — М.: Наука, 1986. — 534 с.
