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

УДК 539.3, 51-72

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

И.Ф. Головнев, Е.И. Головнева, Л.А. Мержиевский1, В.М. Фомин, В.Е. Панин2

Институт теоретической и прикладной механики им. С.А. Христиановича СО РАН, Новосибирск, 630090, Россия 1 Институт гидродинамики им. М.А. Лаврентьева СО РАН, Новосибирск, 630090, Россия 2 Институт физики прочности и материаловедения СО РАН, Томск, 634055, Россия

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

Molecular dynamics study of cluster structure and rotational wave properties

in solid-state nanostructures

I.F. Golovnev, E.I. Golovneva, L.A. Merzhievsky1, V.M. Fomin, and V.E. Panin2

Khristianovich Institute of Theoretical and Applied Mechanics SB RAS, Novosibirsk, 630090, Russia 1 Lavrentiev Institute of Hydrodynamics SB RAS, Novosibirsk, 630090, Russia 2 Institute of Strength Physics and Materials Science SB RAS, Tomsk, 634055, Russia

Molecular dynamics simulation was performed to study the formation of cluster structure, interfaces, and surfaces differing in curvature radius in a perfect nanocrystal when passed through by a nonlinear wave. It is shown that this process is a type of nanostructure self-organization in response to external energy flow with subsequent development of a strong rotational field.

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

Аналитическая теория нелинейных волн структурных трансформаций развита в работах [4-6]. Согласно [7], управляющим параметром нелинейных волн локальных структурных трансформаций является локальная кривизна в сильнонеравновесной структуре. Другими словами, распространение нелинейных волн локализованного пластического течения недислокационной природы является трансляционно-ротационным процессом. Экспериментальная проверка [8] этого положения при одноосном растяжении плоских образцов металлических материалов с наноструктурированными поверхностными слоями подтвердила предсказания нелинейной волновой теории.

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

© Головнев И.Ф., Головнева Е.И., Мержиевский Л.А., Фомин В.М., Панин В.Е., 2014

трансформаций в монокристалле меди под действием внешних приложенных напряжений.

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

Исследование процессов зарождения и развития волновых процессов в наноструктурах проводится на примере медного кластера, имеющего форму прямоугольного параллелепипеда, размеры которого составляют 50x5x5 кристаллических ячеек вдоль осей координат X, У, 2 соответственно. Выбрана ориентация кристалла (1,0,0) для большей ясности получаемых результатов. Для описания межатомного взаимодействия использован широко апробированный потенциал Воутера, полученный в рамках метода внедренного атома [9] для металлов.

Первоначально атомы размещаются в узлах идеальной кристаллической ГЦК-решетки. Затем с помощью метода искусственной вязкости определяются координаты и импульсы системы в состоянии минимума потенциальной энергии. Далее найденные значения координат и импульсов используются в качестве начальных данных.

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

fax =а0 Slnf > С1)

где S — площадь левой грани кристалла; nf — число атомов на этой грани. Так как площадь грани в процессе сжатия может меняться, она пересчитывается на каждом временном шаге, как и внешние силы, действующие на атомы при условии, что а0 = const. Продолжительность импульса T также является внешним контролируемым параметром, как и величина внешнего напряжения а0. Упругая волна достигает правой свободной границы за 6 пс. Для того чтобы явление отражения от свободной границы не усложняло картину процессов, на данном этапе работы исследование всех явлений ограничено временным интервалом T в 6 пс.

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

_ Xi+1 + Xi - Xi+Xi-i 2 2

Рис. 1. Схема формирования мезообъемов

где — среднее значение координат х атомов г&-й плоскости. Знание массива длин мезообъемов в начальный и произвольный момент времени позволяет определять распределение в пространстве локальных относительных удлинений 8 х в любой момент времени.

Поперечные размеры мезообъема рассчитываются, исходя из средних значений координат у и z крайних атомов, лежащих в г&-й плоскости.

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

Главное отличие этого подхода состоит в том, что за атомом не сохраняется номер атомной плоскости. Все пространство разбивается неподвижными плоскостями, перпендикулярными оси X, вдоль которой движется волновой фронт. Расстояние между плоскостями — а 2, где а = 0.3615 нм — размер кристаллической ячейки меди. При таком выборе размеров в невозмущенном состоянии (в начальный момент времени) в мезоячейке находилась одна атомная плоскость (рис. 1).

Длина кристалла составляет 18.075 нм (50 кристаллических ячеек вдоль оси X). В программе используются 1010 ячеек (182.5575 нм), что позволяет учесть смещение кристалла как целого в пространстве.

Рис. 2. Схема элементарного мезообъема для расчета момента импульса подсистемы атомов

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

Для поставленной задачи исследования ротационных свойств волны проводилось дополнительное разбиение описанных выше мезообъемов на мезообъемы в форме прямоугольных параллелепипедов (рис. 2).

Для исследования вращения вдоль оси X все пространство разбивалось на мезоячейки того же объема по схеме, приведенной на рис. 3, а. Мезообъемы, изображенные на рис. 1, дополнительно разбивались на 15 прямоугольных параллелепипедов с образующими, параллельными осям У (рис. 3, б) и 2 (рис. 3, в).

Это позволило рассчитать распределение компонент моментов импульса относительно осей, проходящих через центры мезообъемов с координатами (хс, ). Следует отметить, что эти оси параллельны соответствующим осям координат. Расчет ведется по обычным формулам:

МХ = Е ((У - Ус )Ръ! - (Ъ - ) Ру! X

МСУ = Е ((Ъ - )Рх! - (Х- - Хс )Ръ! X МС2 = Е ((Х- - хс )Ру! - (У! - Ус )Рх!).

Из этих выражений при необходимости можно вывести формулы для угловых скоростей вихрей в соответствующих мезообъемах.

В работе использовалась широко известная скоростная модификация Верле второго порядка точности с шагом по времени 10-16 с. После отключения внешнего

Рис. 3. Построение мезообъемов для исследования ротационного движения вдоль осей X (а), У (б), 2 (в)

напряжения в момент Та система становится замкнутой и изолированной. С этого момента проводился расчет ошибки интегрирования:

Et - ЕТа

егг = •

•100%.

Здесь Ег и ЕТа — энергия системы в момент времени t и Та соответственно. В работе [1] было показано, что с увеличением внешнего напряжения от 0.1 до 4.0 ГПа ошибка возрастает на 9 порядков от 10-10 до 10-1 . При этом необходимо отметить, что в расчетах использовался список Верле с радиусом, превышающим радиус обрезания потенциала Воутера на 0.1 нм, а шаг обновления составлял 10-14 с (или сто шагов по времени). Однако в результате более поздних исследований было обнаружено, что с увеличением внешнего напряжения наблюдаются высокочастотные осцилляции в объеме с периодом, меньшим 10-14 с. В связи с этим было проведено уточнение численной схемы с помощью значительного увеличения радиуса Верле, вплоть до превышения радиуса обрезания потенциала на 0.3 нм. Все дальнейшие результаты расчетов приведены для этого случая. В качестве примера на рис. 4 представлено поведение ошибки на интервале до 50 пс для трех значений внешнего напряжения. Погрешность по энергии в течение интервала в 50 пс (500000 шагов по времени) не превышала 10-5 %.

Рис. 4. Зависимость логарифма ошибки по энергии от времени, внешнее напряжение 0.1 (1), 2.0 (2), 4.0 ГПа (3)

Иллюстрация первых выводов проведена для внешнего воздействия а0 = 4 ГПа в момент времени 6 пс. Вид расположения атомов кристалла в плоскости Х2 приведен на рис. 5.

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

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

После этого был проведен анализ моментов импульсов в описанных выше мезообъемах. Прежде всего необходимо было выяснить, как соотносятся значения

Рис. 6. Расположение атомов структуры в плоскости XZ: t = = 6 пс, а =3.6 ГПа. Черные кружки — структура в состоянии динамического движения, белые кружки — после охлаждения

моментов М с, М£ в перпендикулярных плоскостях. На рис. 7 приведено распределение моментов вдоль оси X. Видно, что значения моментов М£, М£ в мезообъе-мах X совпадают. Кроме того, необходимо отметить, что передний фронт волны с дефектами совпадает с фронтом значительного ротационного движения. В связи с этим ниже будет употребляться термин «ротационная волна».

Иллюстрация распределений этих моментов при фиксированном значении координаты х приведена для двух мезообъемов X с координатами 4.07 и 9.50 нм. Они соответствуют центральной части и фронту ротационной волны (рис. 8).

Как видно, распределения моментов в плоскостях, перпендикулярных оси X, совпадают. Дополнительные исследования показали, что суммы по каждому мезо-объему X равны нулю:

Е м\\ = о, Е мс2 = о. (5)

Рис. 5. Расположение атомов структуры в плоскости XZ: г = = 6 пс, а = 4 ГПа, Та= 6 пс

Рис. 7. Распределение моментов вращения по мезоячейкам кристалла вдоль оси X. Белые кружки — моменты М2, черные — моменты М ¥, г = 6 пс, а = 4 ГПа, Та= 6 пс (момент отключения внешнего напряжения)

-100 -150

-0.5 0.0 0.5 1.0 1.5 2.0

Рис. 8. Распределение моментов вращения по мезоячейкам кристалла вдоль осей У, 2 при х = 4.07 (а), 9.50 нм (б). Белые кружки — моменты М2, черные — моменты MY, г = 6 пс, а = 4 ГПа, Та = 6 пс

В связи с этим далее будет приводиться иллюстрация только по моментам Мг. Распределение моментов Мх приведено на рис. 9. Необходимо отметить, что они значительно меньше по величине моментов

мс, мсх .

Представляет интерес эволюция ротационной волны в интервале времени до 6 пс, когда поток внешней энергии отличен от нуля. Ниже на рис. 10-13 приведена

Рис. 9. Распределение моментов вращения Мх по мезоячейкам кристалла вдоль оси X, г = 6 пс, а = 4 ГПа, Та = 6 пс

Рис. 10. Расположение атомов структуры в плоскости XZ, г = = 2 пс, а = 4 ГПа, Та= 6 пс

Рис. 11. Распределение моментов вращения по мезоячейкам кристалла вдоль оси X, г = 2 пс, а = 4 ГПа, Та= 6 пс

Рис. 12. Расположение атомов структуры в плоскости XZ, г = = 5 пс, а = 4 ГПа, Та= 6 пс

Рис. 13. Распределение моментов вращения по мезоячейкам кристалла вдоль оси X, г = 5 пс, а = 4 ГПа, Та = 6 пс

§ о.з- и 2 02&

^0.10.0-: С

-1- 1 4 -1-1- 8 12 X, нм 16

Рис. 14. Состояние кристалла в момент времени t = 2 пс, а = = 4 ГПа, Га = 6 пс. Распределение скорости ¥сх центров масс (а) и температуры мезоячеек по оси X (б)

Рис. 15. Состояние кристалла в момент времени t = 5 пс, а = = 4 ГПа, Га = 6 пс. Распределение скорости ¥сх центров масс (а) и температуры мезоячеек по оси X(б)

динамическая картина развития ротационной волны в кристалле в моменты времени 2 и 5 пс. На рис. 10, 12 приведена структура расположения атомов в пространстве в плоскостиXZ, а рис. 11, 13 — распределение моментов импульсов по Х-мезообъемам вдоль оси X.

Для определения скорости распространения ротационной волны наиболее удобными характеристиками являются распределение температуры и скорости центров масс атомов в мезоячейках вдоль оси X по кристаллу. На рис. 14, 15 представлены эти распределения в моменты времени 2 и 5 пс в интервале до 6 пс. Зная координату определенной точки волны в определенный момент времени, можно определить скорость распространения волны под действием внешнего напряжения.

Средняя скорость волны (Ут), найденная по температурному пику в интервале от 1 до 6 пс, составляет 1.3 км/с. Проведенный анализ позволяет кроме определения скорости исследовать релаксационный период формирования ротационной волны. До 2 пс на фронте формируется плато шириной 4нм с линейным возрастанием скорости вещества. Далее на расстоянии примерно в 1 нм скорость увеличивается примерно в 4 раза. Одновременно над плато формируется пик температуры до 50 К. Далее до 4 пс плато массовой скорости переходит в S-образную структуру с отрицательной скоростью. Одновременно формируется температурная волна с острым пиком свыше 300 К.

Общепринятая теория деформационных дефектов в твердых телах разработана для термодинамически равновесной трансляционно-инвариантной структуры трехмерного кристалла. Наноструктурные твердые тела являются термодинамически неравновесными системами и к ним традиционная теория деформационных дефектов неприменима. Образование кластерной структуры в нанокристаллах понижает внутреннюю энергию таких сильнонеравновесных систем и обусловливает развитие локализованного пластического течения в виде вихревых нелинейных волн локальных кластерных трансформаций [7].

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

На рис. 16 приведена кривая зависимости термодинамического потенциала Гиббса от функции распределения кластеров различных атомных конфигураций в области наноструктурных состояний [2]. Представленный на рис. 1 исходный нанокристалл соответствует сильнонеравновесному состоянию с функцией распреРис. 16. Зависимость термодинамического потенциала Гиббса от функции распределения п = f (V), где V — молярный объем, в области наноструктурных состояний

деления п0 = f (V) на рис. 16. При наложении внешнего поля распад неравновесного состояния п0 = f (V) может происходить в рамках широкого спектра метаста-бильных состояний. Другими словами, сильнонеравновесная наноструктура может переходить в равновесное состояние (с учетом внешнего поля) по многочисленным каналам. Это и предсказывает молекулярная динамика за фронтом распространения нелинейной волны: возникновение кластерной структуры, многочисленных границ раздела и поверхностей с различными радиусами кривизны.

Аналитическое рассмотрение механизмов структурного распада сильнонеравновесного кристалла во внешнем поле проведено в [3-6]. В соответствии с [6], волновые уравнения для безразмерных величин потока дефектов J и плотности дефектов а имеют вид:

и и а

д 2 ^ и и а

с 2 дг1 дх2

д 1п иа (х, г)

д2 аа

с2 дг2

Е дх,

^ С-—рв с

д2 ос!а

дх2 д2

/ с!v — с

дху дхх ав дхх

при условии совместности источников: = 0,

дН, „ дG

дг Ыт дх/

где G — правая часть выражения (6); Н — правая часть выражения (7); и(х, г) — неупругие смещения в концентрационной волне.

Правая часть уравнения (6) характеризует источники потока дефектов в концентрационной волне. Они определяются скоростью квазиупругой деформации

| (Е, Е - Ев Сав) Е

В скобках представлена разность внутренних напряжений сжатия и сдвига на фронте концентрационной волны. Релаксационные процессы перестроения дефектов типа нанокластеров различных атомных конфигураций представлены в (7) членом 1/ Е Р^С)^ •

Правая часть уравнения (7) характеризует источник плотности деформационных дефектов в концентрационной волне. Им является завихренность

К = Эх Е - Р )С^ Е сдвиговой деформации, вызванной релаксацией сдвиговых напряжений в локальных зонах гидростатического растяжения. Характер волновых потоков деформационных дефектов определяется правыми частями уравнений (6), (7). Величина завихренности к определяет ротационную составляющую нелинейного волнового потока и формирование кластерной структуры среды при дисперсии нелинейных волн.

Сопоставление аналитической теории [4, 6] и результатов молекулярно-динамического моделирования поведения нанокристалла меди в условиях приложения внешнего напряжения показывает полное их соответствие. В то же время, особого внимания заслуживает результат моделирования, предсказывающий возникновение кластерной структуры в сильнонеравновесных нанокристаллах в полях внешних воздействий. Этот очень важный эффект в литературе обычно не учитывается.

Молекулярно-динамические исследования позволили выявить ряд особенностей поведения наноструктуры меди при воздействии внешнего сжимающего напряжения.

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

Показано, что передний фронт ротационной волны является комплексом сформированного множества границ раздела с различными радиусами кривизны.

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

Рассчитана скорость движения фронта ротационной волны, которая составляет примерно 1300 м/с.

Показано, что фронт ротационной волны составляет примерно 7 нм и имеет сложную структуру: в передней части фронта массовая скорость отрицательна (вещество движется против движения волны, что связано с формированием кластерной структуры), а затем скорость становится положительной (400 м/с) и по амплитуде превышает скорость на переднем фронте примерно в 2 раза (-200 м/с). Над фронтом с отрицательной скоростью формируется температурный пик, в котором температура выше 300 K, что примерно в 3 раза выше, чем за фронтом (около 100 K).

Работа выполнена при финансовой поддержке междисциплинарного Интеграционного проекта СО РАН №78.

Литература

Golovnev I.F., Golovneva E.I., Merzhievsky L.A., Fomin V.M. Defect generation as a phenomenon of structure self-organization under external loads // Phys. Mesomech. - 2013. - V. 16. - No. 4. - P. 294302.

Infeld E., Rowlands G. Nonlinear Waves, Solitons and Chaos. - Moscow: Fizmatlit, 2006. - 478 p.

Egorushkin V.E., Panin V.E., Savushkin E.V., Khon Yu.A. Highly excited states in crystals // Russ. Phys. J. - 1987. - V. 30. - No. 1. - P. 5.

волновые структурные носители пластической деформации и разрушения // Физ. мезомех. - 2013. - Т. 16. - № 3. - C. 7-26. Panin V.E., Egorushkin V.E. Curvature solitons as generalized wave structural carriers of plastic deformation and fracture // Phys. Meso-mech. - 2013. - V. 16. - No. 4. - P. 267-286.

Panin V.E., Egorushkin V.E., Panin A.V. Nonlinear Wave Processes in a Deformable Solid as a Multiscale Hierarchically Organized System // Usp. Fiz. Nauk. - 2012. - V. 182. - No. 12. - P. 1351-1357.

Поступила в редакцию 29.11.2013 г., после переработки 26.05.2014 г.

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

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

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

Мержиевский Лев Алексеевич, д.ф.-м.н., проф., гнс ИГиЛ СО РАН, merzh@hydro.nsc.ru

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

Панин Виктор Евгеньевич, д.ф.-м.н., акад. РАН, советник РАН, зав. лаб. ИФПМ СО РАН, paninve@ispms.tsc.ru