CC BY
81
Расчет термодинамических свойств наноструктур методом молекулярной динамики
И.Ф. Головнев, Е.И. Головнева, В.М. Фомин
Институт теоретической и прикладной механики им. С.А. Христиановича СО РАН, Новосибирск, 630090, Россия
В работе проведено исследование термодинамических характеристик металлического нанокластера методом молекулярной динамики для изобарического и изохорического процессов. Показано, что свойства наноструктуры и свойства макротел различаются влиянием поверхностных атомов в наноструктурах. Предложенная методика позволяет рассчитывать термомеханические свойства и для макротел. Для этого необходимо увеличить размеры наносистемы до масштаба, при котором можно пренебречь влиянием поверхностных атомов.
Molecular dynamics calculation of thermodynamic properties of nanostructures
I.F. Golovnev, E.I. Golovneva and V.M. Fomin
S.A. Khristianovich Institute of Theoretical and Applied Mechanics SB RAS, Novosibirsk, 630090, Russia
The molecular dynamics method is applied to study the thermodynamic characteristics of a metal nanocluster for the isobaric and isochoric processes. It is shown that the nanostructure properties differ from the macrobody properties due to the influence of surface atoms in nanostructures. The proposed technique allows the thermomechanical properties to be also calculated for macrobodies. To do this, the nanosystem size should be enlarged to the scale when the influence of surface atoms could be neglected.
1. Введение
В связи с бурным развитием нанотехнологий в настоящее время возникла потребность в определении физических и механических свойств наноструктур и, в частности, уравнения состояния в широком диапазоне температур и давлений. Проблема расчета уравнения состояния является важнейшей и для макротел, но к настоящему времени даже эту задачу удается решить только полуэмпирическими методами (см., например, обзор [1]).
В случае наноструктур проблема еще более усложняется. В частности, невозможно использовать приближение трансляционной симметрии, на котором построена вся теория твердого тела.
Наноструктуры отличаются от обычных материалов тем, что количество поверхностных атомов, их энергия и энергия связи с объемными атомами сопоставима с
аналогичными характеристиками объемных атомов. В то же время, свойства поверхностных и объемных атомов сильно различаются. Наиболее прямой путь — это построение термодинамики систем, состоящих из двух подсистем: объемной и поверхностной. Однако традиционный подход в случае наноструктур не применим по следующей причине. Термодинамика содержит в себе постулат аддитивности: энергия термодинамической системы есть сумма энергий ее отдельных частей.
Полную энергию наноструктуры Ut можно представить в виде суммы:
Ut = Usa + Uv + Ub, где UV — энергия объемных атомов; Usa — энергия поверхностных атомов; Ub — энергия связи этих подсистем. Для применения традиционного подхода, т.е. использования постулата аддитивности, необходимо, чтобы энергия связи была много меньше энергий UV
© Головнев И.Ф., Головнева Е.И., Фомин В.М., 2007
и Ц,а. Однако для наносистем это условие не выполняется.
Все это обусловило проведение исследований термодинамических свойств наноструктур с детализированным выделением подсистем объемных и поверхностных атомов.
2. Применение молекулярно-динамического подхода к построению уравнения состояния наноструктур
В методе молекулярной динамики процессы исследуются на микроуровне, т.е. рассчитываются траектории атомов в фазовом пространстве. На основе этой наиболее полной информации далее рассчитываются все необходимые характеристики.
В рамках статистической физики получено выражение для свободной и внутренней энергии твердого тела в гармоническом приближении [2]:
3^ к(ії
F = и„(V) + кТ X 1п —і
і=1 кТ
(1)
Е = и0(7) + 3МТ, (2)
где и0 (V) — потенциальная энергия атомов твердого тела при Т = 0.
Зная свободную энергию, давление можно выразить в виде:
\ =_ЭЦс _кт ^ ± Эю,
ЭV 1Т дV “1 ®г- дV'
Р = -
(3)
Здесь -Эи°
— «холодное» давление;
ш 1 йт
-кТ 5 —
І=1 ю,- ду
дУ
«тепловое» давление.
Так как расчет зависимости ш( V) даже в простейших случаях идеальных кристаллов представляет чрезвычайно сложную задачу, сумму в выражении для давления представляют в виде произведения числа частот на не-
1
которое среднее значение '
ыдУ
1 Эш
Р = ^ _ кТ= Рс - кТ3№
ду £Ю,- ду с \адУ
(4)
или в безразмерном виде: 3ЫкТ IУ Эю \юЭУ
Р = Р --
у
Определив величину у как у = -
1 А1ТГ (5)
\rndV /
уравнение состояния записывается в форме Ми-Грю-найзена [3, 4]:
3МТ
Р - Рс =У:
у
(6)
Часто давление связывают с внутренней энергией твердого тела:
Р - Р =
У( Е - Ес) У '
Полученные выражения послужили в настоящей работе основой для общего подхода к расчету калорического Е = Е(V, Т) и термического Р = P(V, Т) уравнений для наноструктур с помощью метода молекулярной динамики.
В работе в качестве исследуемого объекта был выбран медный нанокластер сферической формы радиусом 2 нм. Выбор такой геометрической структуры исследуемого физического объекта обусловлен тем, что в наноструктурах, имеющих форму прямоугольного параллелепипеда, поверхностные атомы, расположенные на гранях, ребрах и в вершинах, имеют большой разброс по когезионной энергии, т.е. очень неравноправны, в отличие от кластеров сферической формы.
Основная суть методологии учета влияния поверхностных атомов основана на учете координационного числа каждого атома. В ходе численного моделирования ведется учет количества ближайших соседей у каждого атома и проверяется условие равенства этого количества координационному числу идеальной ГЦК-структуры. Тем самым выделяются поверхностные и объемные атомы, для которых проводится детализированное изучение термомеханических характеристик. Кроме того, параллельно рассчитывались такие же характеристики для всей системы в целом.
Для учета взаимодействия атомов внутри сферы использовался потенциал Вотера для меди [5].
Для подготовки начальных данных вся система охлаждалась с помощью метода искусственной вязкости, в результате чего ансамбль атомов в конце расчета находился в минимуме потенциальной энергии. Конечные массивы координат и импульсов, полученные в процессе охлаждения, использовались далее в качестве начальных данных при моделировании квазистатичес-кого сжатия и растяжения сферы.
3. Расчет «холодных» характеристик медной наносферы
На первом этапе необходимо провести расчеты холодных характеристик при температуре Т = 0 К. Давление У0 и энергия наноструктуры для нулевого внешнего давления выбираются в качестве начала отсчета. Далее, используя метод искусственной вязкости, проводятся расчеты «холодных» энергии и давления структуры для заданного набора внешних давлений Рсг-. Этот метод позволяет в конце расчетного интервала иметь статически равновесные состояния при Т = 0 К, т.е. получить набор значений У1, рС, (V,), исг- (V,) в необходимом диапазоне давлений или изменений объема х, = У1/У0.
В качестве внешнего параметра задавалось давление Р0 на поверхности. Расчет сил /а, действующих на поверхностные атомы и направленных вдоль радиуса
Ди„, 10-17 Дж
0 100000 200000 1ЧТ
Рис. 1. Относительное изменение полной энергии объемной части системы как функция числа шагов по времени для трех разных значений внешних сжимающих давлений
-0.1 0 0.1 х,
Рис. 3. Зависимость внешнего давления от относительного изменения полного объема для интервала давлений от -20 до 12 ГПа
от центра (растяжение сферы) или к центру сферы (сжатие сферы), проводили по формуле:
/ = (ЗД/Ж8, (8)
где — площадь поверхности; N. — полное число поверхностных атомов.
Для гашения ударно-волновых процессов в сфере применяли метод линейного нарастания сил /а от нуля до заданного значения и искусственную вязкость. В результате в конце расчета параметры наносферы имели статически равновесное состояние при температуре Т = = 0 К (рис. 1). Далее конечные координаты системы использовали для расчета необходимых характеристик.
В качестве независимой переменной для удобства выбрано относительное изменение объема
*«= (У - Уг)/Ут- (9)
Для сжатия сферы эта величина отрицательна, для растяжения — положительна.
Расчеты растяжения сферы проведены в интервале от 0.1 до 12 ГПа, сжатия — от 0.1 до 20 ГПа.
При этом проводился детализированный анализ состояния сферы как единого целого, объемной и поверхностной частей. Индекс а в этом случае принимает обозначения ^ V, s соответственно.
Физический анализ энергетических характеристик удобно проводить для изменений потенциальных энер-
гий по отношению к начальной потенциальной энергии, в качестве которой выступает энергия охлажденной сферы
Аиа= иа- иай. (10)
На рис. 2 приведена такая зависимость для полной энергии системы от относительного изменения объема. Зависимость близка к квадратичной, хотя и несколько отличается от нее. Это обусловлено тем, что для больших значений внешнего давления отклонения атомов от положения равновесия уже велики и сказывается нелинейная часть потенциала межатомного взаимодействия. Аналогичные характеристики рассчитывались для всех подсистем, описанных во введении.
Таким образом, в результате проведенных расчетов получена зависимость «холодной» внутренней энергии от объема.
Второй важнейшей «холодной» характеристикой является зависимость давления от объема. На рис. 3 приведена зависимость внешнего давления от относительного изменения полного объема для интервала давлений от -20 до 12 ГПа. Нелинейность, как и в случае энергий, обусловлена нелинейностью потенциала межатомного взаимодействия. Для удобства на рис. 4 такой же график приведен для внешних давлений, изменяющихся от -1 до 1 ГПа. Зависимость линейна и удобна для проведения оценок.
Ли,, 10-17 Дж
Рис. 2. Изменение потенциальной энергии системы по отношению к
начальной потенциальной энергии как функция относительного изме- Рис. 4. Зависимость внешнего давления от относительного изменения
нения объема полного объема для интервала давлений от -1 до 1 ГПа
Полученные численно характеристики полностью определяют «холодную» составляющую уравнения состояния наноструктуры.
4. Расчет «тепловых» характеристик нанокластера
Второй этап заключается в получении тепловых составляющих давления и энергии. В соответствии с (2), (3), если провести разогрев системы при V = const, т.е. для изохорического процесса, то можно получить численно функции
PT = Y(T, V = const)
3NkT V ’
ET (T, V- = const)
без каких-либо упрощающих предположений об ангармонических вкладах потенциала и влияния объема на константу Грюнайзена и тем самым построить калорическое и термическое уравнения для наноструктуры из первых принципов. Моделирование изохорического процесса проводили следующим способом. В качестве начальных данных брали координаты и импульсы атомов в конечных статически равновесных состояниях для фиксированного внешнего давления Рсі. Далее поверхностные атомы помещали во внешний потенциал, который препятствовал увеличению расстояния от атома до центра сферы. Этот прием обеспечивал постоянство объема всей наносферы.
Численный разогрев наноструктуры проводили в интервале температур от 0 до 1000 К. Для этого использовали метод стохастических сил [6]. Для построения зависимости всех характеристик от температуры через каждые 25 К стохастическая сила отключалась и система релаксировала к термодинамически равновесному состоянию в течение 10_12 с. На этом же временном интервале проводили усреднение всех характеристик по тепловым флуктуациям. В результате таких расчетов найдены все необходимые термодинамические характеристики системы, в том числе тепловые доли давления и энергии на численной сетке (Тк, У,). На этом этапе температуру рассчитывали следующим образом:
3 N Дп2
- MTkm = Ekin 2 i=i 2m
(15)
где Ар, = р, - т,V, — хаотическая составляющая импульса г-го атома; р, — импульс г-го атома в лабораторной системе координат; т1 — его масса; V, — скорость центра масс системы. Так как разогрев проводился в интервале от 25 до 1000 К, далее используется индекс расчетной температурной сетки Т, по которому легко определяется соответствующая хаотическая температура системы Тк1п = Т ■ 25 К.
Данный процесс является изохорическим для всей системы (рис. 5, 1 ГПа). В то же время, видно, что объем поверхностной системы увеличивается, а объем внут-
Рис. 5. Изменение объема разных подсистем в зависимости от индекса температурной сетки (вся наносфера — х^ объемная подсистема — ху, поверхностная подсистема — Xs)
ренней части системы уменьшается, и для подсистем процесс неизохорический. Таким образом, этим способом можно моделировать изохорические процессы для всей системы в целом.
4.1. Анализ внутренней энергии
Наиболее фундаментальной величиной в термодинамике является полная внутренняя энергия системы Е1п, имеющая механический аналог и включающая в себя как потенциальную энергию взаимодействия атомов (в том числе и «холодную»), так и энергию хаотического движения атомов Ек1п.
В качестве примера на рис. 6 приведена зависимость приращений энергий Ди = и(Т, V) - Е0 и Ект всей системы от индекса температурной сетки для случая начального сжатия сферы. Их значения близки, что говорит о наличии равновесного термодинамического состояния. Особенно важно отметить, что полученная зависимость линейна (теплоемкость — константа). При
Рис. 6. Зависимость приращения потенциальной энергии Аи и энергии теплового движения атомов всей системы £йп от индекса температурной сетки для случая начального сжатия сферы (Ди — сплошная линия; Еш — штриховая линия)
Рис. 7. Зависимость полной внутренней энергии от индекса температурной сетки для внешних сжимающих давлений: 0.5 (штрихпунк-тирная линия), 1 (сплошная линия), 5 ГПа (штриховая линия)
Рис. 9. Зависимость температуры Ту, определенной через выражение (16), для полной системы (сплошная линия) и подсистем объемных (штрихпунктирная линия) и поверхностных атомов (штриховая линия) от индекса температурной сетки
малых температурах значения AU и Ekin совпадают, а при увеличении температуры наблюдается расхождение зависимостей, что связано с увеличением влияния ан-гармонизма (поверхностной части).
Зависимость полной внутренней энергии от индекса температурной сетки приведена на рис. 7 для внешних сжимающих давлений 0.5, 1, 5 ГПа. Зависимость внутренней энергии от температуры также линейна и графики совпадают для всего интервала давлений.
Как известно, для не очень высоких температур, пока не сказывается ангармоническая часть потенциала, кристалл рассматривается как система ^трехмерных гармонических осцилляторов и температура может определяться по формуле:
3NkTv = Ein. (16)
Для анализа выбора выражений для определения температуры системы приведена зависимость Tkin (рис. 8) и Tv (рис. 9) для полной системы и подсистем объемных и поверхностных атомов. Как видно, температуры подсистем, определенные по выражению (15),
совпадают с большой точностью. В то же время значения температуры Ту, рассчитанной по выражению (16), для поверхности и объема уже начинают различаться для температур порядка 250 К. Таким образом, несмотря на установление термодинамического равновесия во всей системе температуры Ту подсистем отличаются, что противоречит второму началу термодинамики. Поэтому везде ниже под температурой будем понимать величину, найденную по выражению (15).
Дополнительные ограничения, которые необходимо учитывать при анализе термодинамических явлений в наносистемах, связаны с использованием классической механики. При температурах, меньших температуры Дебая, необходимо учитывать квантовые эффекты. Следовательно, полученные в работе выводы верны для температур Т > 0О. Для меди 0О = 315 К.
Полученные результаты позволяют построить уравнение состояния в форме Ми-Грюнайзена (7). На рис. 10 приведено семейство кривых, дающих зависимость теплового давления от внутренней энергии для
Рис. 8. Зависимость температуры Гкіп, рассчитанная через £кіп, для Рис. 10. Зависимость теплового давления Р от внутренней энергии
атомов объема, атомов поверхности и для всей системы от индекса Еіп^ для различных значений начального «холодного» давления
температурной сетки (или от объема наноструктуры)
Рис. 11. Изменение теплового давления АР = Р - Рс с температурой для различных начальных «холодных» давлений
Рис. 12. Зависимость коэффициента Грюнайзена у от температуры Т для различных начальных «холодных» давлений (или объемов)
различных значений начального «холодного» давления или от объема наноструктуры (случай сжатия). Видно, что с увеличением внутренней энергии, или, что то же самое, с увеличением температуры, коэффициент Грю-найзена уменьшается.
Аналогичные исследования можно провести и для уравнения состояния (6). Так как объем полной системы при разогреве сохраняется, то изменение давления обусловлено только повышением температуры. На рис. 11 приведено изменение теплового давления АР = Р - Рс с температурой для различных начальных «холодных» давлений. В интервале от 0.1 до 1 ГПа графики практически совпадают, а при более высоких давлениях видна параметрическая зависимость от начального «холодного» давления или от объема сферы.
Используя (6), коэффициент Грюнайзена определяли в соответствии с выражением
у = ДР/ (3ШТ/¥).
На рис. 12 приведено семейство кривых у от Т для различных начальных «холодных» давлений (или объемов). Видно, что у зависит как от давления, так и от температуры. Для Т> 300 К (для которых справедливы настоящие расчеты) коэффициент Грюнайзена изменяется от 2.0 до 2.18 ГПа в интервале начальных давлений от 0.1 до 10 ГПа. Экспериментальное значение коэффициента Грюнайзена для макроскопических образцов меди — 1.84.
5. Выводы
Таким образом, в настоящей работе предложена методика расчета термодинамических свойств наноструктур методом молекулярной динамики, основанная на
первых принципах, которая апробирована для случая сферических нанокластеров меди.
Исследование показало, что результаты для объемной части системы совпадают с результатами для макрообъектов. Отличие свойств наноструктуры от свойств макротел вызвано значительным влиянием поверхностных атомов в наноструктурах.
Предложенная методика позволяет рассчитывать термомеханические свойства и для макротел. Для этого необходимо увеличить размеры наносистемы до масштаба, при котором можно пренебречь влиянием поверхностных атомов.
Работа выполнена при поддержке научной школы НШ-8732.2006.1, гранта РФФИ № 05-01-00211а, Фонда содействия отечественной науке, гранта молодежных проектов Лаврентьевского конкурса СО РАН, в рамках интеграционного проекта СО РАН № 28.
Литература
1. Фомин В.М., Гулидов А.И., Сапожников Г.А. и др. Высокоскоростное взаимодействие тел. - Новосибирск: Изд-во СО РАН, 1999. -600 с.
2. Ландау Л.Д., Лифшиц Е.М. Теоретическая физика. T. V. Статистическая физика. - М.: Наука, 1964. - 567 с.
3. Mie G. Zur kinetischen Theorie der einatomigen Körper // Ann. Phys. -
1903. - V. 316. - No. 8. - P. 657-697.
4. Grüneisen E. Zustand des festen Körpers // Handbuch der Physik. -Bd. 10. Thermische Eigenschaften der Stoffe. - Berlin: Springer Verlag, 1926.- P. 1-59.
5. Voter A.F. Embedded Atom Method Potentials for Seven FCC Metals: Ni, Pd, Pt, Cu, Ag, Au, and Al // Los Alamos Unclassified Technical Report LA-UR-93-3901.
6. Болеста А.В., Головнев И.Ф., Фомин В.М. Исследование процесса
столкновения сферического кластера меди с жесткой стенкой методом молекулярной динамики // Физ. мезомех. - 2001. - Т. 4. -№ 1. - С. 5-10.
Поступила в редакцию G2.GS.2GG7 г.
