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

УДК 539.3, 51-72

Молекулярно-динамическое исследование роли поверхности в процессе разрушения наноструктур

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

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

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

Ключевые слова: наноструктура, поверхность, разрушение, растягивающее напряжение, молекулярно-динамическое моделирование

Molecular dynamics study into the role of the surface in fracture

of nanostructures

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

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

The paper reports on a molecular dynamics study of the fracture of ideal nanostructures and processes occurring in surface and bulk atomic systems. It is shown from the first principles that fracture begins in the surface layer and then a crack is generated in the bulk.

Keywords: nanostructure, surface, fracture, tensile stress, molecular dynamics simulation

1. Введение

Впервые, на основании экспериментальных исследований, академик А.Ф. Иоффе [1] сделал вывод, что разрушение твердых тел под влиянием внешних нагрузок начинается с разрастания микротрещин и дефектов на поверхности, которые затем приводят к образованию трещины, распространяющейся в объем. Более подробный анализ этого явления он привел в монографии [2].

Современные экспериментальные и теоретические исследования влияния поверхности на характеристики материалов [3, 4] углубили представление о роли поверхности и границ раздела и подтвердили выводы А.Ф. Иоффе.

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

ки, опирающегося на экспериментальные данные о строении вещества и потенциалы, описывающие взаимодействия атомов в нем.

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

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

Система поверхностных и объемных атомов разбивается на мезообъемы, в которых рассчитываются такие

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

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

2. Физико-математическая модель

Методика и комплекс программ по молекулярно-динамическому исследованию процессов распространения возмущений в системе поверхностных и объемных атомов проходили отладку на примере медного кластера, имеющего форму прямоугольного параллелепипеда, размеры которого составляли 50x5x5 кристаллических ячеек вдоль осей координатX, Y, Z соответственно. Выбрана ориентация кристалла (1, 0, 0) для большей прозрачности получаемых результатов. Для описания межатомного взаимодействия использован широко апробированный многочастичный потенциал, полученный в рамках метода внедренного атома [6] для ГЦК-металлов.

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

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

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

fax =-СТ0 Slnf (1)

где S — площадь левой грани кристалла; nf — число атомов на этой грани.

Подвижный зажим моделировали следующим образом. Атомы крайней левой атомной плоскости полагали свободными, и на них действовали внешние силы (1), направленные в отрицательном направлении оси X

(растягивающее напряжение). Поскольку площадь грани в процессе растяжения может меняться, ее пересчитывали на каждом временном шаге, как и внешние силы, действующие на атомы.

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

V(r) = x -x0)2,

(2)

/ 0 0 0ч . « «

где (х, , у, , zi) — координаты г-го атома крайней правой грани после охлаждения.

Расчеты проведены для интервала внешних возмущений от 5 до 30 ГПа. Продолжительность импульса Гс не ограничивали, и напряжение действовало вплоть до разрушения структуры.

3. Мезоаиализ распространения волны в кристалле

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

Хм + X, X, + Xi_

AL, =

i-1

(3)

2 2' У ' где х, — среднее значение координат х атомов г-й плоскости. Знание массива длин мезообъемов в начальный и произвольный момент времени позволяет определить распределение в пространстве локальных деформаций

Ех (4)

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

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

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

4. Мезоаиализ распространения волн в подсистемах объемных и поверхностных атомов

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

нм2 2.8

2.0

1.20.4

0.4 0.0

0 0.4 0.8 1.2 1.6 Y, нм

Рис. 1. Расположение атомов в плоскости, перпендикулярной оси X

но 12, как у идеальной бесконечной ГЦК-решетки. Если число ближайших соседей меньше 12, то по построению это поверхностные атомы. В качестве примера на рис. 1 приведен пример одной из атомных плоскостей. Черные кружки соответствуют атомам с координационным числом 12, а серые — с меньшим 12. Как видно, рисунок подтверждает правильность выделения поверхности.

Второй проблемой является выбор границы системы объемных и поверхностных атомов. Для выбора границы между объемными и поверхностными атомами использовался принцип геометрической равноправности объемов этих атомов. Границу в плоскости YZ проводили ровно посредине между атомами этих видов. Для определения границы поверхности использовали условие совпадения расчетной и экспериментальной плотности вещества (в данном случае меди). Оказалось, что внешнюю границу надо отнести на 0.025 нм от координат поверхностных атомов (точнее ядер этих атомов). Поэтому алгоритм расчета площади поперечных сечений был следующим. Находили средние значения координат у и г крайних объемных и поверхностных атомов. Координаты разделяющих поверхностей находили как полусумму этих координат. Это позволило определить площадь поперечного сечения системы объемных атомов. Полную площадь сечения структуры находили на основании знания координат поверхности (средние значения координат поверхностных атомов в плоскости и добавленный отступ). Площадь сечения поверхностных атомов определяли как разность полного сечения и сечения объемных атомов. На рис. 2 в качестве примера приведено распределение размеров площадей сечения в начальный момент времени.

Как видно, их величина отличается примерно в 3 раза. По этой причине на первом этапе представляет интерес проведение анализа сил, действующих в подсистемах атомов, т.к. они не опираются на определение

Рис. 2. Зависимость площади сечения атомных плоскостей от координаты атомной плоскости для поверхностных (1) и объемных атомов (2)

площадей. Силы рассчитывали следующим образом. Определяли координату хр1 между двумя последовательными атомными плоскостями [р и [р +1 и на основании известного потенциала [6] рассчитывали силы, действующие на атомы с координатой меньшей хр1 со стороны атомов с координатой большей хр1. Далее для каждого набора плоскостей силы суммировали отдельно для объемной и поверхностной системы атомов. В качестве примера на рис. 3 приведено распределение этих сил вдоль наноструктуры.

Первое, что необходимо отметить, что х-проекции сил, действующих на системы поверхностных и объемных атомов, лежащих в одной атомной плоскости, точно совпадают по модулю и противоположны по направлению. Таким образом, равенство нулю суммы х-проекций сил обеспечивают разные подсистемы атомов. Второе — несимметричность на концах структуры связана с несимметричностью расчета сил («действие справа налево») в сочетании с использованием радиуса сферы взаимодействия потенциала Воутера [6], которое усиливает несимметричность. С учетом отмеченного выше несовпадения площадей подсистем на рис. 4 приведено распределение напряжений ахх в поверхностном слое и в объеме. Эту характеристику определяли по модели Коши как отношение силы к площади поверхности.

Как видно, компоненты имеют разные знаки и по модулю отличаются примерно в 3 раза. Детализированный расчет характеристик, описанных выше, проводили в течение всего интервала исследований.

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

рх, 10-8 н

1 0 --1-2

0 4 8 12 16 X, нм Рис. 3. Распределение сил в пространстве для поверхностных (1) и объемных атомов (2)

2

1

12 16 X, нм

20 10 0 --10

t = 0

2

Г

0 4 8 12 16 X, нм

Рис. 4. Распределение компоненты тензора напряжений а вдоль оси X для поверхностных (1) и объемных атомов (2)

нии структуры, атомы перемещаются в междоузлия, и идеальные атомные плоскости разрушаются.

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

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

5. Макро- и мезоанализ разрушения наноструктуры

Под макроанализом здесь понимается определение характеристик, присущих кристаллу как целому. На пер-

Рис. 5. Вид структуры в плоскости XX. Внешнее напряжение а 0 = 16 ГПа, момент времени t = 3 пс

Рис. 6. Зависимость внешней силы со стороны неподвижного зажима (а) и изменение полной потенциальной энергии (б) от времени

вом этапе представляет интерес общий вид структуры после наступления разрушения. Ниже приведены результаты на момент времени действия импульса силы t = 3 пс, что меньше времени прихода волны на закрепленный зажим. Это означает, что волна, отраженная от неподвижного зажима, не могла вернуться на подвижный зажим. Следовательно, на первом этапе рассматривается наиболее простой случай разрушения структуры. Внешний вид разрушившегося кристалла на этот момент времени наблюдается при наименьшем значении напряжения а0 = 16 ГПа (рис. 5). Далее будет проведен анализ, не является ли это значение критическим для описанного выше внешнего возмущения.

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

Таблица 1

Номер точки

1

Момент времени, пс

0.01

0.51

0.80

1.14

2.25

2.69

2.98

ахх> ГПа

В табл. 1 приведены моменты времени, соответствующие номерам точек.

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

На интервале 4-5 сила начинает возрастать быстрее, чем на упругом интервале 1-2, а на интервале 56-7 выходит на плато, причем в точке 6 по модулю она достигает минимального значения. На этом же интервале потенциальная энергия растет слабо, что говорит о пластических необратимых деформациях, а после точки 7 начинает уменьшаться. Это говорит об интенсивном разрушении структуры. При этом важнейшим оказывается вопрос о применимости приближения атомных плоскостей для мезоанализа.

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

мере, для первых плоскостей. Как показало подробное исследование — это первые пять атомных плоскостей.

После того как проведен анализ используемых приближений для расчета мезохарактеристик, можно приступать к исследованию процессов в системах поверхностных и объемных атомов. На рис. 8 приведено распределение изменений локальных напряжений Аахх = = ахх - а0^ по отношению к начальному моменту времени на поверхности и в объеме.

Для начального интервала времени (рис. 8, а, б) пик возмущения в структуре, вызванного внешним приложенным напряжением, распространяется быстрее в системе поверхностных атомов. Затем возмущение (волна), распространяющееся в объемной системе, догоняет поверхностную волну. Однако приращение напряжения в поверхностных атомах больше примерно на 2 ГПа. Необходимо отметить, что со временем ширина фронта возрастает от 1 до 7 нм. Сама величина приращения напряжения примерно равна внешнему приложенному напряжению (16 ГПа), а осцилляции соответствуют колебаниям атомных плоскостей. В то же время представляет интерес полная величина напряжения в подсис-

15 X нм

Рис. 7. Распределение дисперсии координат атомов в атомных плоскостях по структуре в различные моменты времени. Внешнее растягивающее напряжение а 0 = 16 ГПа

15 X, нм

Рис. 8. Распределение локального напряжения в объеме (черная линия) и на поверхности (серая линия) по структуре в различные моменты времени. Внешнее растягивающее напряжение а 0 = 16 ГПа

Аахх, ГПа

0 5 10 15 X, нм

Рис. 9. Распределение локального напряжения в объеме (черная линия) и на поверхности (серая линия) по структуре в момент времени 1.14 пс. Внешнее растягивающее напряжение 16 ГПа

темах атомов. На рис. 9 приведена иллюстрация этих характеристик в критический момент времени.

Как видно, напряжение в поверхностном слое примерно в 3 раза выше, чем в объемном, и превышает критическое значение напряжения, полученное выше (а 0 = 16 ГПа). Для дальнейшего представляет интерес распределение температур в наноструктуре вдоль оси X.

Как указывалось выше, мезоанализ в области разрушения проводится в неподвижных в пространстве мезо-ячейках. В качестве примера на рис. 10 приведено распределение температур в пространстве в разные моменты времени. Температуры рассчитывали усреднением по общей системе поверхностных и объемных атомов. Для лучшего разрешения области разрушения результаты приведены для мезоячеек с координатами до 4 нм.

В конце интервала линейного растяжения структуры (точка 2, I = 0.51 пс) температура в первой мезоячейке превышает 300 К. На интервале развития пластической деформации (2- 4, t = 1.14 пс) наблюдается резкое увеличение температуры до 1100 К. Сам процесс разрушения (удаление фрагмента на рис. 10, в) сопровождается повышением температуры до 3000 К. Далее за счет теплообмена температура в области разрушения падает до 1600 К, с одновременным ее ростом в соседних мезообъемах.

Анализ развития разрушения проводился на основе расчета количества атомов объемной и поверхностной

4000-

2000-

800 0

г = 2.34 пс :_^и_

...... > к у г = 2.69 пс 1 1

-4

-2

0

X, нм

Рис. 10. Распределение температур в объеме (жирная линия) и в поверхности (тонкая линия) в мезоячейках структуры в критические моменты времени

систем в неподвижных мезоячейках. Для того чтобы сравнить эту характеристику для подсистем, использовали их отношение к начальному количеству атомов. Подробный анализ этой характеристики через 0.01 пс показал, что вначале трещина появилась в системе поверхностных атомов (рис. 11, а). Для лучшего разрешения области разрушения результаты приводятся в пространственном интервале от -3 до +3 нм.

В момент времени 2 пс трещина образовалась в объемной части. В этой же ячейке поверхность разру-

X, нм

X, нм

Рис. 11. Распределение отношения числа атомов в поверхностной (серая линия) и объемной (черная линия) подсистемах к их начальному значению в этих подсистемах по мезоячейкам в момент времени 1.45 (а), 2.00 (б), 2.02 (в), 2.34 пс (г)

шена лишь частично (рис. 11, б). Произошло так называемое разрушение под поверхностью.

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

6. Выводы

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

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

Анализ этих характеристик позволил выделить как упругий режим растяжения наноструктуры, продолжающийся до достижения значений напряжения 16 ГПа и локальной деформации 0.25, так и три необратимых (пластических, но без наличия дислокаций) режима деформации. На первом интервале сила начинает прирастать медленнее со временем, чем на упругом интервале, а потенциальная энергия начинает линейно возрастать с ростом времени нагружения. На втором интервале прирост силы можно приблизительно считать линейным, но угол наклона много больше, чем на упругом и на первом пластическом интервале. Этот интервал характеризуется линейным, но более медленным приростом потенциальной энергии системы (образование дефектов). Третий интервал характеризуется «плато» (или почти постоянным значением силы, несмотря на рост деформации структуры) и очень слабым нарастанием потенциальной энергии. Этот интервал характеризуется развитием процесса разрушения.

Критическое значение усредненной по системам поверхностных и объемных атомов компоненты прироста напряжения соответствует внешнему приложенному напряжению 16 ГПа. Изменение значения компоненты прироста напряжения в системах поверхностных и объемных атомов близко к внешнему приложенному напряжению 16 ГПа, однако для поверхности эта величина больше примерно на 2 ГПа. Абсолютное значение напряжения на поверхности превышает критическое значение внешнего напряжения примерно в 2 раза и примерно в 3 раза выше значения напряжения в объеме. Температура в мезообъемах, в которых происходит разрушение, превышает температуру плавления.

Первая трещина появляется в поверхностном слое. Это равнозначно факту, что разрушение начинается с поверхности.

Работа выполнена при поддержке программы фундаментальных исследований ОЭММПУ РАН № 12.

Литература

1. Иоффе А.Ф. Физика кристаллов. - М.-Л.: ОГИЗ, 1929. - 250 с. Ioffe A.F. The Physics of Crystals. - New York: McGraw-Hill Book Company, Inc., 1928. - 198 p.

2. Иоффе А.Ф. Отчет о работе Физико-технического института // УФН. - 1936. - Т. XVI. - № 7. - С. 847-871.

Ioffe A.F. A progress report from Physical-Technical Institute // UFN. -1936. - V. XVI. - No. 7. - P. 847-871.

3. Поверхностные слои и внутренние границы раздела в гетеро-генныж материалах / Под ред. В.Е. Панина. - Новосибирск: Изд-во СО РАН, 2006. - 520 с.

Surface Layers and Internal Interfaces in Heterogeneous Materials / Ed. by V.E. Panin. - Novosibirsk: Izd-vo SO RAN, 2006. - 520 p.

4. Панин В.Е., Сергеев В.П., Панин А.В. Наноструктуирование поверх-

ностных слоев конструкционных материалов и нанесение нано-структурных покрытий. - Томск: Изд-во ТПУ, 2009. - 285 с. Panin V.E., Sergeev V.P., Panin A.V. Nanostructuring of Surface Layers in Structural Materials and Deposition of Nanostructured Coatings. - Tomsk: Izd-vo TPU, 2009. - 285 p.

5. Golovnev I.F., Golovneva E.I., Fomin V.M. The influence of the surface on the fracture process of nanostructures under dynamic loads // Comput. Mater. Sci. - 2015. - V 97. - P. 109-115.

6. 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.

7. Allen M.P., Tildesley D.J. Computer Simulation of Liquids. - New York: Clarendоn Press, 1987. - 385 p.

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

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

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

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

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