Особенности применения методов механики сплошных сред для описания наноструктур

Особенности применения методов механики сплошных сред для описания наноструктур

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

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

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

Peculiarities of application of continuum mechanics methods to the description of nanostructures

E.I. Golovneva, I.F. Golovnev, and VM. Fomin

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

Using a molecular dynamics method we study the possibility of applying continuum mechanics to the description of nanostructures. It is demonstrated that the nanostructure should be considered as a heterostructure consisting of two subsystems. One subsystem presents atoms of the surface, and the other atoms of the bulk. In this connection continuum mechanics methods can be applied for either the entire nanostructure considered as a single “elementary” volume or a heterostructure consisting of subsystems with different properties. In both cases, mechanical characteristics are to be calculated in the framework of a molecular dynamics method.

1. Введение

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

[1]. Следовательно, эти структуры имеют чрезвычайно

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

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

[2].

В качестве внешнего воздействия рассматривалось одноосное растяжение кристалла (подробно методика изложена в работах [3, 4]).

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

Для получения конечного статически равновесного состояния в расчетах использовалась искусственная вязкость, когда на каждый атом действовала сила трения Г = -^р, где р — импульс атома. Это позволило в рамках метода молекулярной динамики смоделировать конечное состояние, в котором сумма сил, действующих на атом, равна нулю [5].

2. Исследование взаимосвязи характеристик микро- и мезомасштабного уровней

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

Для расчета этих характеристик необходимо, прежде всего, знание площади крайних граней и площади поперечного сечения кристалла, определение длины кристалла Ь и введение параметра I — характерной длины области внутри кристалла, где рассчитывается величина ст хх. Первая проблема, которая при этом возникает, каким образом определить границу кристалла.

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

а а а

L0 = пга +-------+— = пга +—,

0 x 4 4 x 2

где пх — число кристаллических ячеек вдоль оси X.

Длина кристалла под воздействием внешней нагрузки меняется и в расчетах определяется следующим образом. Находились средние значения х — координат атомов левой и правой крайних граней:

Е X

(x) =-

где суммирование ведется по координатам атомов грани, а п — количество атомов крайних граней. Тогда длина кристалла

L = к)-< xi) + ^

индексы 1 и г относятся к левой и правой граням. Относительное удлинение кристалла как целого («элементарный объем»):

L - L0 Lo

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

5 =

Тогда полный объем кристалла:

к = ЬБ = х) г -(х) 1 + - * -(уг) 2 -(уг) 1 +

Х[(2г)2 -(1 + §] .

При исследовании локальных характеристик (внутри кристалла) площадь поперечного сечения между гранями с номерами а и а + 1, перпендикулярными оси X, рассчитывается аналогично (5). Длина параллелепипеда, включающего в себя эти грани, согласно требованиям соответствия собственному объему атомов и осреднения массы атомов по всему мезообъему, определяется следующим образом:

(х,Л +(хЛ „ (хЛ +(хЛ

/ _ \ Ч а+1 \ ч а+2 \ ч а-1 \ ч а т\

а,а+1 _ 2 2 ’

Тогда локальная деформация рассчитывается по формуле:

є, =-

І - І

*а,а+І *а,а+І

а,а+І

где /а а+! — начальная длина контрольной мезоячейки.

Сравним выражения для относительного удлинения всего кристалла єь и локального удлинения є /. Пусть (х^ — средняя координата а-й грани. Если кристалл не деформирован, то

2

В случае упругой одноосной деформации кристалла вдоль оси X, когда левая грань закреплена:

(x) а=( xi) l+Axa

где

Лха = C

есть смещение под воздействием внешней силы (С = = const и определяется модулем Юнга). Подставляя эти выфажения в формулу для локальной деформации, получаем:

2C

= —.

Кг =

L

n

f

а

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

2С-1 1

є ь =---1 1----

ь а , N

= є /1 1 —

1

N

где N — число атомных плоскостей в нанокристалле, перпендикулярных оси X. Если число кристаллических ячеек вдоль оси X равно пх, то N = 2пх +1.

Видно, что

2С

11Ш еь = е 1 =-,

N а

и для больших кристаллов эти характеристики совпадают. В частности, для исследуемого кристалла пх = 3

и е Ь = 6ег/7.

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

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

и = г — Г0 =Дгг (Г0).

Далее находятся соответствующие производные и компоненты тензора деформации по известной формуле континуальной механики

1

игк = 2

диі

дхк

ди

к + ди/ ди/

дХ: дХ: дх,

или для малых деформаций:

диі

дхк

ди

дХ:

к )

Особо следует отметить, что данное определение не требует введения границы кристалла. Таким образом, в рамках молекулярной динамики рассчитываются все геометрические характеристики кристалла на основе данных для координат атомов (дискретной системы). С другой стороны, движение атомов описывается динамическими законами системы материальных точек (закон Ньютона, уравнения Лагранжа или уравнения Гамильтона). Для сопоставления с результатами механики сплошных сред необходимо перейти к понятию напряжений, т.е. провести некоторое усреднение сил, действующих между точечными частицами, отнеся их к единице площади.

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

= оь£

где щ — количество атомов этой грани, а 5 определялась по формуле (5).

Расчет локальной компоненты тензора напряжений стхх для одноосного растяжения кристалла проводился двумя способами. Согласно определению это сила в направлении оси X, действующая на единицу площади, перпендикулярной оси X. Пусть /х — полная сила, действующая на атомы, находящиеся по одну сторону сечения площади 5 со стороны атомов, расположенных по другую сторону этого сечения. Тогда

ст = А

хх 5 '

Полная сила находилась по формуле

/х = ЕР*,

^ 3

7—г х ^ ^

где р — х-я компонента силы, действующей на ^ю частицу со стороны 3-й, частицы находятся по разные стороны контрольного сечения.

Второй способ — использование формулы, полученной в рамках кинетической теории [6]:

-’ар

(

Е

Га Гв 1 * +- Е ГИ

т

іф і

Здесь Кс — объем исследуемой системы; р“ — а-я

компонента импульса ьй частицы; т

масса атома;

Г

— Р-я компонента радиус-вектора между ^й и 3-й частицами. Суммирование ведется по всем частицам, находящимся в исследуемом объеме. В случае статического равновесия эта формула приобретает вид:

ав

Если в качестве объема Ус рассматривать объем, включающий только две кристаллические грани, то, используя предыдущие обозначения, получаем:

(

-Е Fxrx 2 /-і іі 1

2 іФ 1 )

1 / /

= —2 /х/ = —. 2К К

Здесь I — расстояние между атомными плоскостями. С другой стороны, объем, соответствующий паре атомных плоскостей, можно выразить как

Кс = 521,

следовательно,

ст = /х

2Б'

т.е. получаем в два раза меньший результат, чем при использовании прямого определения (19) для локального напряжения.

Если формула (22) применяется для объема кристалла в целом (как к элементарному объему), то

п

+

и

О 1СО ЭТИ

Рис. 1. Зависимость относительного удлинения всего кристалла е ь (сплошная линия) и относительного локального удлинения е/ (пунктирная линия) от числа шагов по времени N т (103). Величина шага т = 10-16 с. Внешнее напряжение ст ь = 05 ГПа

т.е. получаем то же соотношение с локальным определением компоненты тензора напряжений (22), что и для компонент тензора деформаций (13).

3. Результаты расчета на микроуровне

В расчетах атомы левой грани кристаллического куба

0 1ГО 14,

Рис. 2. Зависимость внешнего приложенного напряжения сть (сплошная линия) и локального напряжения ст/ (пунктирная линия) от числа шагов по времени Nz (103). Величина шага т = 10-16 с. Внешнее напряжение ст ь = 0.5 ГПа

Р

ас

0 Ш МО N.

Рис. 3. Зависимость внешней приложенной силы Р (сплошная линия) и локальной силы во внутренней мезоячейке (пунктирная линия) от числа шагов по времени Nz (103). Величина шага т = 10-16 с. Внешнее напряжение ст ь = 0.5 ГПа

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

Далее для различных значений внешнего напряжения ст Ь была найдены значения относительного удлинения кристалла в состоянии статического равновесия (рис. 5). Аналогичные зависимости были получены и для локальных характеристик е / и ст/. По установившимся значениям относительного удлинения и известным значениям напряжения были построены о-е-диа-граммы (рис. 6). Как видно, эти кривые не совпадают. Например, модули Юнга для этих диаграмм имеют значения ЕЬ = 116.51, Е/ = 96.07. Таким образом, здесь особенно ярко выделяются две концепции подхода механики деформируемого твердого тела к описанию поведения наноструктур, в которых дискретность и поверхностные эффекты имеют большое значение. Если рассматривать систему как «элементарный» объем, то

а ,,

1»

150-  - -

й 1<Ю ИЮ 14,

Рис. 4. Зависимость площади крайней правой грани (10-16 см2) (сплошная линия) и площади поперечного сечения внутренней мезо-ячейки (пунктирная линия) от числа шагов по времени Nx (103). Величина шага т = 10-16 с. Внешнее напряжение а ь = 0.5 ГПа

Рис. 5. Зависимость относительного удлинения всего кристалла еь от числа шагов по времени Nx

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

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

ди ди у

х. = 0.0049, у

ди.

дх

ду д,

= —0.00199,

и3з = 0 при i ф].

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

Рис. 6. а-£-диаграммы бездефектного нанокристалла меди: аь - єь-сплошная линия; а/ - Є/ — пунктирная линия

Рис. 7. Зависимость компонент векторов смешения (Е) атомов в кристалле для одноосного растяжения. сть = 0.5 ГПа

ложение о разбиении всего стержня на «элементарные» объемы. Относительная деформация вдоль оси X, найденная в соответствии с выражением е = их = 0.00495, с большой точностью совпадает с локальным значением е/ = 0.0043, найденным по формуле (8).

Из выражения ст хх = Еихх находим модуль Юнга Е = = 99.8 ГПа, что также близко к значению «локального» модуля Е/.

Таким образом, при переходе от дискретной модели на микроуровне к континуальной модели среды возможны две концепции. Первая — это рассмотрение наноструктуры как «элементарного» неделимого объема с соответствующими механическими параметрами, описывающими адекватно отклик этого объекта на внешнее механическое возмущение. Второй — применение основных положений механики деформируемого тела к расчету процессов внутри нанообъекгов. При этом необходимо использовать «локальные» механические параметры, найденные также в молекулярно-динамических расчетах. Результаты, полученные в рамках обоих подходов (см. (14), (26)), совпадают, когда число кристаллических ячеек в наноструктуре становится много больше единицы.

4. Влияние поверхности на свойства наноструктур

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

Рис. 8. Зависимость изменения потенциальной энергии единицы объема кристалла(108 Дж/м3) от величины е Ц2 (10-5) при одноосном растяжении

ческих характеристик наноструктур. Однако не менее важным фактором, влияющим на свойства наноструктур, являются поверхностные эффекты. В качестве примера проведем дополнительный энергетический анализ одноосного растяжения кристалла (см. п. 3). В дополнение к расчету геометрических и силовых характеристик определяли изменение потенциальной энергии Д и = и — и0, где и — потенциальная энергия кристалла [2] в состоянии статического равновесия для заданного внешнего напряжения аь, а и0 — та же энергия недеформированного кристалла. С другой стороны, независимо от этого определялись, как описано выше, относительные удлинения є ь. Для проверки известного выражения

— = Е ^

V 2

построена зависимость Аи/У от величины е2Ь/2 (рис. 8). Как видно эта зависимость линейна, а модуль Юнга Е = 114.3 ГПа, что с большой точностью совпадает с ЕЬ, найденным выше для кристалла как единого «элементарного» объема. Здесь особо необходимо отметить, что величина А и относится ко всему кристаллу в целом и содержит в себе как энергию поверхностных атомов, так и объемных, т.е. это усредненная по подсистемам характеристика.

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

В работе исследовались сферические кластеры меди с радиусами до 70 Е. На первом этапе находились координаты атомов, соответствующие идеальной ГЦК-структуре. Далее эти координаты использовались как начальные данные для охлаждения кристалла (поиска состояния, в котором полная энергия системы минимальна) [5]. Полученные в результате расчета конечные

Рис. 9. Зависимость отношения числа атомов на поверхности Д8Г к полному числу атомов N в сферическом нанокластере меди от радиуса сфер

координаты использовались для анализа влияния поверхностных эффектов.

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

На рис. 9 приведена важнейшая характеристика наноструктур — зависимость отношения числа атомов на поверхности к полному числу атомов от радиуса сферы. Как видно, даже для сфер с радиусом в 70 Е это отношение близко к 0.1, а не стремится к нулю, как для макрообъектов.

Для анализа энергии подсистем атомов, находящихся на поверхности и в объеме, необходимо использовать наиболее общее выражение:

и г = и 8 + иу + и5У.

Здесь [2]

и г =Е Р

2

іФ 1

есть полная энергия кристалла;

и5 =1 Р

iGs г, jGs

^ 3

есть энергия всех атомов поверхности, рассчитанная в отсутствии атомов объема;

и у =1Р

г'еу 2 г,з'еу

3

есть полная энергия атомов в объеме, найденная при отсутствии атомов поверхности; и — энергия связи подсистем (поверхность - объем).

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

Рис. 10. Зависимость средней энергии

риментальному значению для когезионной энергии и для сфер с радиусом 70 Е отличаются от последнего на 1.2 %.

В качестве дополнительной характеристики наноструктур может служить отношение энергии поверхности к полной энергии кластера и к энергии объемных атомов (рис. 11). Для сфер радиуса 10 Е энергия поверхности достигает половины полной энергии кристалла и даже превышает энергию объема. С увеличением радиуса до 70 Е доля энергии атомов поверхности уменьшается, но составляет 10 % от полной энергии. Это доказывает, что нельзя пренебрегать поверхностными эффектами в наноструктурах с размерами до 100 Е.

Поверхностную энергию, в соответствии с (28), определим как сумму энергий атомов поверхности и энергию связи с атомами объема:

или

и5г = и - иу.

На рис. 12 приведена зависимость отношения этой величины к площади поверхности сферы (удельная по-

Рис. 12. Зависимость отношения поверхностной энергии к площади поверхности

верхностная энергия) от радиуса сферы. Как видно, этот коэффициент достаточно универсален, т.к. почти не зависит от радиуса сфер и равен 10.8 Дж/м2.

Коэффициент поверхностного натяжения определялся в соответствии с наиболее общим выражением:

У =

и

д5

в котором изменение площади вызвано приложенным внешним напряжением. В связи с этим было проведен следующий численный эксперимент. Сфера радиуса 20 Е нагружалась внешним контролируемым давлением, которое изменялось в интервале от нуля до 2 ГПа. Зная площадь, количество поверхностных атомов и давление, находились внешние силы, приложенные к поверхностным атомам и направленные по радиусу сферы. В соответствии с описанной выше методикой находилось статически равновесное состояние, в котором определялись радиус сферы, ее площадь и объем, а также изменение полной энергии, энергии объемных атомов и поверхностной энергии (33). Это позволило построить зависимость изменения поверхностной энергии от изменения площади поверхности (рис. 13) и найти коэффициент поверхностного натяжения у На рис. 14 приведена зависимость у от внешнего давления. Его величина для малых давлений, где зависимость близка к линейной, совпадает с экспериментальным значением коэффициента поверхностного натяжения для расплава

Рис. 11. Зависимость отношения энергии атомов поверхности к полной энергии кристалла

Рис. 13. Зависимость изменения поверхностной энергии сферы (10-15 мкДж) от относительного изменения площади поверхности А5 = (5 - 50)/50, где 50 — начальная площадь поверхности

Рис. 14. Зависимость коэффициента поверхностного натяжения меди

меди. В различных источниках это значение варьируется от 0.5 до 1.2 Дж/м2.

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

K = -V

дР_

dV

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

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

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

Рис. 15. Зависимость коэффициента объемного сжатия кристалла меди

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

Работа выполнена при поддержке Российского фонда фундаментальных исследований (грант № 05-01-00211а) и гранта Президента Российской Федерации для поддержки ведущих научных школ № НШ-2282.2003. 1.

Литература

1. Нанотехнология в ближайшем десятилетии. Прогноз направления

исследований / Под ред. М.К. Роко, Р.С. Уильямса, П. Аливисато-са. - М.: Мир, 2002. - 292 с.

2. Johnson R.A. Alloy models with the embedded-atom method // Phys. Rev. B. - 1989. - V. 39. - P. 12554-12559.

3. Головнев И. Ф., КоневаЕ.И., Фомин В.М. Численное моделирование

разрушения бездефектных кристаллов при динамических нагрузках // Физ. мезомех. - 2001. - Т. 4. - № 5. - С. 1-7.

4. Головнева Е.И., Головнев И.Ф., Фомин В.М. Молекулярно-динамический анализ динамического разрушения наноструктур // Физ. мезомех. - 2003. - Т. 6. - № 2. - С. 37-46.

5. Головнева Е.И., Головнев И. Ф., Фомин В.М. Моделирование квази-

статических процессов в кристаллах методом молекулярной динамики // Физ. мезомех. - 2003. - Т. 6. - № 2. - С. 37-46.

6. ЗубаревД.Н. Неравновесная статистическая термодинамика. - М.:

Наука, 1971. - 415 с.

T