Молекулярно-динамическое исследование термической неустойчивости в наносистемах Текст научной статьи по специальности «Физика»

УДК 539, 536, 51-72

Молекулярно-динамическое исследование термической неустойчивости в наносистемах

И.Ф. Головнев, Е.И. Головнева, В.М. Фомин

Институт теоретической и прикладной механики им. С.А. Христиановича СО РАН, Новосибирск, 630090, Россия

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

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

Molecular dynamics study of thermal instability in nanosystems

I.F. Golovnev, E.I. Golovneva, and V.M. Fomin

Khristianovich Institute of Theoretical and Applied Mechanics SB RAS, Novosibirsk, 630090, Russia

Molecular dynamics simulation has been used to study the behavior of nanostructures at elevated temperatures. The procedure of heating a system up to a given temperature is detailed. A temperature region has been revealed in which the system is thermally unstable, its crystal structure is disrupted, but the system remains in the solid phase, no melting occurs, and the system temperature rises spontaneously.

Keywords: thermal instability, metals, nanostructure, molecular dynamics simulation, temperature effect

1. Введение

Исследованию влияния температуры на механические характеристики материалов с помощью метода молекулярной динамики уделяется достаточно большое внимание. Например, в работе [1] исследуется влияние температуры на предел сдвиговой прочности А1 и Си. В работах [2, 3] исследуется влияние температуры на разрушение углеродных нанотрубок и графена. Авторы статьи [4] исследовали влияние тепловых флуктуаций на критическое напряжение и разрушение микроструктуры при постоянной внешней нагрузке твердых тел с парным межатомным потенциалом Леннарда-Джонса, который не описывает металлические системы. Работа [5] посвящена исследованию влияния температуры на механические свойства биметаллической нанопроволо-ки бесконечной длины (периодические граничные условия вдоль оси деформации). При этом моделируется квазистатическая деформация, что не позволяет исследовать динамическое разрушение. В работах [6, 7] ис-

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

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

© Головнев И.Ф., Головнева Е.И., Фомин В.М., 2017

ным размером при криогенных температурах; 2) разогрев структуры до заданной температуры и релаксация системы для формирования термодинамического равновесного состояния.

2. Генерация начального состояния системы

Поведение наноструктур при нагревании исследовалось на примере нанопроволоки. Для моделирования нанопроволоки строился идеальный кристалл меди в форме прямоугольного параллелепипеда с числом кристаллических ячеек пх = 50, пу = п2 = 5 вдоль соответствующих осей. Выбрана ориентация кристалла (1,0, 0). Первоначально атомы находились в узлах идеальной кристаллической ГЦК-решетки. Взаимодействие атомов описывалось потенциалом, полученным по методу внедренного атома. В данной работе выбран потенциал Воутера [8].

В работе использовалась широко известная скоростная модификация Верле второго порядка точности с шагом по времени 10-16 с [9]. Отметим, что в случае изолированной системы погрешность по энергии не превышает 10-5 % для интервала времени до 50 пс.

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

3. Разогрев структуры до заданной температуры

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

В отличие от работ [13-18], в которых имитируется действие случайных сил на атомы, в развиваемом подходе проводится имитация поглощения атомами квантов энергии со случайными импульсами (аналог поглощения теплового излучения). Атомы системы через временной интервал т(Ь поглощают дополнительный импульс Арг-. Временной интервал можно выразить через шаг по времени в схеме и число шагов, через которые определяются случайные силы :

Ть = Аь(1)

С помощью датчика случайных чисел были найдены компоненты приращения импульса атомов по схеме

Арц =Аро(0.5 - f), 1 = х, у, г, (2)

где / — случайные числа в интервале (0, 1). Поскольку число атомов наноструктуры N конечно, то для того чтобы полный импульс системы не изменялся, в про-

грамме использовался прием, позволяющий перейти к приращениям компонент импульса Ар ^, сумма которых по всей системе равна нулю:

N

, (3)

г=1

Щ =ЬРъ - Р\ • (4)

Кроме того, чтобы система не приобретала дополнительный момент импульса МсЬ, после пересчета всех импульсов со случайными добавками рассчиты-

вался полный момент импульса

и N

мсЬ =Е[Г-Рг- ] (5)

1=1

и находились (с учетом что система — кластер и начало координат совпадает с центром масс) компоненты угловой скорости , которые появились за счет конечного и не всегда очень большого числа атомов системы:

М f =1 .

(6)

После этого пересчитывались скорости всех атомов по формуле

V= V- [Пг,- ] (7)

и определялись конечные импульсы атомов.

Приращение полной кинетической энергии системы за временной интервал т(Ь

AE¿ = £ i=i

(Pi +Ap ¿ Г 2m

N p-Pi =

i=i 2m

= £ рг- Api Ap¿ i=i m i=i 2m

(8)

Поскольку импульсы атомов и случайные приращения импульсов не коррелированы, то первое слагаемое в правой части (8) равно нулю и приращение кинетической энергии имеет вид N A62

AEk =1^- (9)

i=i 2m

Если учесть, что приращение пропорционально амплитуде Ap0 (т.е. Api = Ap0£i, где £i — случайный вектор с длиной меньшей единицы), то производная кинетической энергии по времени может быть записана в виде

Ae¿ =^PL £ JL.

At NTh т i=12m' Таким образом, скорость разогрева определяется двумя параметрами: амплитудой случайного импульса Ap0 и числом шагов по времени, через которые происходит разогрев системы.

Далее будет более удобно пользоваться системой единиц измерения физических величин, используемых в расчетной программе: единица длины L — 10-8 см, единица времени T — 10-13 с, единица массы m — 10-24 г. Следовательно, единица скорости

v =

L-10-8 см

T-10-13 с

™ _ . „5 СМ _

-= V -10 — = v км/с,

Рис. 1. Зависимость температуры системы от времени

единица силы F =

mL т - 10-27% - 10-10

Т2 • 10-26

кг • м/ с2 = F-10-11 Н.

Отсюда можно получить единицу импульса р = р -10-24 Н-с.

В расчетах использовался шаг по времени 10-16 с. Следовательно, т = 0.001.

Ниже понадобится единица энергии

Е =

т -

т -10-24 %2 -10-16 г - см2 Т2 -10-26 с2

= Е -10-14 эрг = 10-21 Дж.

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

В первом блоке происходит разогрев системы на величину АТС по описанному выше методу. Численные

Рис. 2. Зависимость изменения потенциальной энергии системы от температуры

Е, 10"21 Дж 120000 8000040000 0

Т, К

Рис. 3. Зависимость изменения внутренней энергии системы: потенциальной энергии (1) и хаотической составляющей кинетической энергии атомов (2) от температуры

исследования показали, что наиболее оптимальным является выбор АТС от 5 до 25 К.

Во втором блоке система эволюционирует свободно и релаксирует к термодинамическому равновесию в течение времени Ыге1 т.

В третьем блоке система также свободна и проводится усреднение по тепловым флуктуациям в течение времени Ыау т. После вывода информации расчет возвращается к первому блоку.

В проведенной работе система «разогревалась» со следующими значениями параметров: амплитуда случайного импульса Ар0 = 4, число шагов по времени между стохастическим прогревом = 20, время частичной релаксации (число шагов по времени) Ыге1 = = 500, время усреднения Ыау = 500.

4. Разогрев нанопроволоки

Для приведенных параметров схемы на рис. 1 приведена зависимость температуры системы от времени, а на рис. 2 — зависимость изменения потенциальной энергии от температуры. Потенциальная энергия отсчи-

Рис. 4. Зависимость температуры системы в процессе релаксации от времени

Е, 10"21 Дж

80000 70000 60000 50000

0

20 40

60

80 Ь пс

Рис. 5. Зависимость изменения внутренней энергии системы: потенциальной энергии (1) и хаотической составляющей кинетической энергии атомов (2) от времени. Температура разогрева 500 К

Рис. 6. Вид наноструктуры в плоскости XX после релаксации. Температура после разогрева 500 К

Рис. 7. Вид наноструктуры в плоскости XX после охлаждения

тывается от начального состояния системы при температуре около 0 К. В районе 1300 К наблюдается перегиб, соответствующий фазовому переходу «твердое тело - жидкость». Очевидно, что выше 1000 К исследования термомеханических свойств твердой фазы проводить не имеет смысла.

Анализ наличия термодинамического равновесного состояния наиболее удобно провести путем сравнения компонент внутренней энергии структуры — изменения потенциальной энергии системы А и и хаотической составляющей кинетической энергии атомов Е^п. В случае когда гармоническое приближение для колебательных степеней свободы выполняется (что верно до 1000 К с достаточной точностью), эти составляющие равны. На рис. 3 приведено такое сравнение.

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

5. Релаксация системы для формирования термодинамического равновесного состояния

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

На рис. 4 приведена зависимость температуры системы, а на рис. 5 — изменения потенциальной энергии системы Аи и хаотической составляющей кинетической энергии атомов Екп от времени в процессе релаксации для начальной температуры после разогрева, равной 500 К. Изменение потенциальной энергии, как и рань-

ше, отсчитывается от начального состояния системы при температуре около 0 К.

Как видно, в пределах температурных флуктуаций эти параметры системы остаются постоянными. Следовательно, система уже в процессе разогрева достигла стационарного, термодинамически равновесного состояния. Выясним, в чем причина отличия между компонентами внутренней энергии. На рис. 6 приведен вид наноструктуры после релаксации в плоскостиXX. Видно сильное отклонение от идеальной кристаллической структуры при криогенной температуре (рис. 7).

В связи с этим необходимо посмотреть на пространственное распределение атомов в нанопроволоке после релаксации. На рис. 8 приведен вид парной функции распределения для системы при температуре 500 К. Там же для сравнения эта функция приведена для этой же системы при криогенной температуре. Как видно, нагрев только до 500 К приводит к значительной амор-физации системы. Следовательно, разница между компонентами внутренней энергии обусловлена значительным влиянием дефектов кристаллической структуры. Несмотря на это, структура находится в стационарном состоянии.

Проведение исследования релаксации при более высоких температурах разогрева привело к неожиданному результату. В качестве иллюстрации ниже приведен пример для температуры разогрева 700 К. На рис. 9 приведена зависимость изменения полной и потенциальной энергии системы от времени в процессе релаксации. Для большей наглядности изменения отсчи-

0.2 0.3 0.4 0.5 0.6 г, нм

Рис. 8. Парная функция системы после релаксации — черная линия. Температура 500 К. Серая линия — идеальный ГЦК-кристалл

Аи, 10"21 Дж

-20000

0

-10000- Лк л.

0

20

40 60 80 г, пс

Рис. 9. Зависимость изменения полной (1) и потенциальной энергии (2) в процессе релаксации от времени. Начальная температура после разогрева 700 К

тываются от начального момента релаксации. Видно, что через 5-10 пс потенциальная энергия системы начинает понижаться, в то время как полная энергия с большой точностью остается постоянной. Процесс перестройки длится примерно до 40 пс, а затем потенциальная энергия, с точностью до тепловых флуктуаций, остается постоянной, т.е. система находится в стационарном состоянии с более низкой потенциальной энергией. При этом температура системы в процессе перестройки растет (рис. 10) более чем на 100 К.

На рис. 11, а приведен вид структуры в плоскости XX сразу после разогрева до 700 К, а на рис. 11, б — после релаксации. Как видно, система перешла из состояния, в котором еще сохранялись атомные плоскости, в состояние, похожее на жидкий кристалл с нитевидными кластерами, ориентированными вдоль оси X (рис. 12).

Этот переход можно объяснить наличием термического неустойчивого состояния системы, которое наступает при определенной температуре. Для обнаружения этой критической температуры были проведены детализированные расчеты релаксационного процесса через 25 К, начиная с 500 К. На рис. 13 представлены результаты изменения потенциальной энергии и температуры системы от времени.

Как видно, критическое значение температуры, при которой возникает неустойчивость в нанопроволоке с поперечным размером 5x5 кристаллических ячеек, — 525 К. Подробное исследование процесса релаксации после разогрева показало, что температура, аналогично вышеприведенной иллюстрации, растет на определенную величину АТСГ в течение временного интервала от 10 до 40 пс вплоть до температур разогрева, превышающих 2000 К. Величина АТСГ с ростом температуры

-4 0 4 Х,ям

Рис. 11. Вид наноструктуры в плоскости XX после разогрева до температуры 700 К (а) и после релаксации (б)

разогрева уменьшается до нуля. Принципиальное отличие интервала температур до фазового перехода «твердое тело - жидкость» состоит в том, что в процессе релаксации нанопроволока сохраняет геометрическую форму, хотя бы отдаленно напоминающую прямоугольный параллелепипед, а выше температуры фазового перехода образец, изначально имеющий форму прямоугольного параллелепипеда, приобретает сферическую форму. Для температур превышающих 2000 К структура приобретает сферическую форму уже в процессе разогрева. Эти свойства принципиально разделяют явление термической неустойчивости от явления фазового перехода «твердое тело - жидкость».

Рис. 10. Зависимость температуры от времени в процессе релаксации. Температура разогрева 700 К

Рис. 12. Вид наноструктуры в плоскости УХ после релаксации. Температура разогрева 700 К

AU, 10-21 Дж 8000

0

T, K

600-

500-

1 | а

"IF'

0 20 40 60 80 t, пс

- 1

2

' ■ ■ г " 1 1л '

0

20

40

60

t, пс

Рис. 13. Зависимость изменения потенциальной энергии (а) и температуры структуры (б) от времени в процессе релаксации. Температура разогрева: 500 (1), 525 (2), 550 К (5)

Таким образом, благодаря термической неустойчивости возникла «запретная» зона в температурном диапазоне примерно от 525 до 650 К. При нагревании нанопроволоки до заданной температуры в этом интервале происходит спонтанное дополнительное увеличение температуры системы примерно на 100 К. Это приводит к тому, что в интервале температур от 525 до 650 К при молекулярно-динамическом моделировании невозможно задать определенную температуру для на-норазмерной структуры.

6. Заключение

В результате было показано, что уже при температурах порядка 500 К идеальная кристаллическая структура в нанопроволоке с поперечным размером 5x5 кристаллических ячеек переходит в аморфную фазу.

При температурах около 500 К отличие потенциальной и кинетической составляющих внутренней энергии превышает 5 % , а при достижении температуры 1000 К отличие достигает 25 %, что связано с сильным ангар-монизмом колебаний.

Обнаружено, что при температурах выше 525 К наблюдается термическая неустойчивость — система спонтанно переходит в состояние с более глубоким потенциальным минимумом, что сопровождается ростом температуры примерно на 100 К. При этом конечная стационарная структура отдаленно напоминает жидкий кристалл.

В результате термической неустойчивости возникает «запретная» зона в температурном диапазоне примерно

от 525 до 650 K. В нанопроволоке происходит спонтанное увеличение температуры примерно на 100 K.

Работа выполнена при финансовой поддержке РФФИ (грант № 14-01-00465).

Литература

1. Iskandarov A.M., Dmitriev S.V., Umeno Y. Temperature effect on ideal

shear strength of Al and Cu // Phys. Rev. B. - 2011. - V. 84. -P. 224118.

2. Tang Ch., Guo W., Chen Ch. Structural and mechanical properties of partially unzipped carbon nanotubes // Phys. Rev. B. - 2011. - V. 83.-P. 075410.

3. Zhao H., Aluru N.R. Temperature and strain-rate dependent fracture strength of graphene // J. Appl. Phys. - 2010. - V. 108. - P. 064321.

4. Yamamoto A., Kun F., Yukawa S. Microstructure of damage in thermally activated fracture of Lennard-Jones systems // Phys. Rev. E. -2011. - V. 83. - P. 066108.

5. Subramanian K.R., Sankaranarayanan S., Venkat R. Bhethanabotla, Babu Joseph. Molecular dynamics simulation of temperature and strain rate effects on elastic properties of bimetallic Pd-Pt nanowires // Phys. Rev. B. - 2007. - V. 76. - P. 134117.

6. Liu Y., Wang F., Zhao J., Jiang L., Kiguchi M., Murakoshi K. Theoretical investigation on the influence of temperature and crystallo-graphic orientation on the breaking behavior of copper nanowire // Phys. Chem. Chem. Phys. - 2009. - V. 11 - P. 6514-6519.

7. Wang F., Sun W., Gao Y., Zhao J., Sun Ch. Investigation on the most probable breaking behaviors of cupper nanowires with the dependence of temperature // Comp. Mat. Sci. - 2013. - V. 67. - P. 182-187.

8. 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. - 1993.

9. Allen M.P., Tildesley D.J. Computer Simulation of Liquids. - New York: Oxford University Press, 1987. - 385 p.

10. Головнева Е.И., Головнев И.Ф., Фомин В.М. Моделирование квазистатических процессов в кристаллах методом молекулярной динамики // Физ. мезомех. - 2003. - Т. 6. - № 6. - C. 5-10.

11. Головнева Е.И., Головнев И.Ф., Фомин В.М. Молекулярно-динами-ческий расчет коэффициента теплового расширения для наноклас-теров меди // Наносистемы: физика, химия, математика. - 2011. -Т. 2. - № 3. - C. 71-78.

12. Болеста А.В., Головнев И.Ф., Фомин В.М. Плавление на контакте при соударении кластера никеля с жесткой стенкой // Физ. мезомех. - 2001. - Т. 4. - № 1. - С. 5-10.

13. Schneider T., Stoll E. Molecular-dynamics study of a three-dimensional one-component model for distortive phase transitions // Phys. Rev. B. - 1997. - V. 17. - No. 3. - P. 1302-1322.

14. Berendsen H.J.C., Postma J.P.M., van Gunsteren W.F., DiNola A., Haak J.R. Molecular dynamics with coupling to an external bath // J. Chem. Phys. - 1984. - V. 81. - No. 8. - P. 3684-3690.

15. Andersen H.C. Molecular dynamics simulations at constant pressure and/or temperature // J. Chem. Phys. - 1980. - V. 72. - P. 2384-2393.

16. Nose S. A unified formulation of the constant temperature molecular dynamics methods // J. Chem. Phys. - 1984. - V. 81. - No. 1. - P. 511.

17. Hoover W.G. Canonical dynamics: Equilibrium phase-space distributions // Phys. Rev. A. - 1985. - V. 31. - No. 3. - P. 1695-1697.

18. Martyna G.J., Klein M.L., Tuckerman M. Nose-Hoover chains: The canonical ensemble via continuous dynamics // J. Chem. Phys. -1992. - V. 97. - No. 4. - P. 2635-2643.

Поступила в редакцию

__01.09.2016 г.

Сведения об авторах

Головнев Игорь Федорович, к.ф.-м.н., снс, снс ИТПМ СО РАН, golovnev@itam.nsc.ru Головнева Елена Игоревна, к.ф.-м.н., снс ИТПМ СО РАН, elena@itam.nsc.ru

Фомин Василий Михайлович, д.ф.-м.н., акад. РАН, научн. рук., гнс, зав. лаб. ИТПМ СО РАН, fomin@itam.nsc.ru