ЖЭТФ, 2020, том 157, вып. 2, стр. 255-271

ВЗАИМОДЕЙСТВИЕ МЕЖДУ АТОМАМИ АЛЮМИНИЯ

И УГЛЕРОДА НА ГРАНИЦЕ РАЗДЕЛА ФАЗ

АЛЮМИНИЙ-ГРАФЕН И АЛЮМИНИЙ-ГРАФИТ

В. В. Решетнякa*, А. В. Аборкин^b

^a ГНЦ РФ Троицкий институт инновационных и термоядерных исследований

108840, Троицк, Москва, Россия

^b Владимирский государственный университет им. А. Г. и Н. Г. Столетовых

600000, Владимир, Россия

Поступила в редакцию 29 мая 2019 г.,

после переработки 17 сентября 2019 г.

Принята к публикации 19 сентября 2019 г.

С использованием теории функционала плотности и метода классической молекулярной динамики ис-

следовано взаимодействие алюминия, находящегося в жидком и кристаллическом состояниях, с по-

верхностями графена и графита. Для описания взаимодействия между атомами алюминия и углерода

предложена параметризация потенциала Морзе, основанная на результатах расчетов ab initio. Потенциал

использован для теоретического исследования взаимодействия капли расплава с поверхностью графита

(0001). Рассчитаны свойства свободной поверхности расплава алюминия, а также контактной поверхно-

сти, образованной при смачивании графита каплей расплава. Результаты вычислений показали хорошее

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

DOI: 10.31857/S0044451020020066

нию карбидной фазы путем отжига материала при

высокой температуре [15, 16].

Характерным масштабом при расчете взаимо-

1. ВВЕДЕНИЕ

действия атомов вещества является радиус r_c обрез-

В последние десятилетия углеродные наночасти-

ки потенциала, который обычно составляет около

1 нм. Если поверхность частиц образована атомами

цы различной структуры активно используются для

создания гибридных и композиционных материалов

углерода в sp²-состоянии, расположенными в вер-

шинах шестиугольников, а радиус кривизны этой

с заданными свойствами [1-4]. В частности, множе-

ство работ посвящено разработке прочных и лег-

поверхности R ≫ r_c, то с точки зрения локальной

структуры и межатомных взаимодействий она близ-

ких алюмоматричных композитов, использующих в

качестве наполнителя наночастицы, образованные

ка поверхности графена. Кроме графита и графе-

на такой структурой, в частности, характеризует-

атомами углерода в sp²-состоянии [1-3,5-8]. Иссле-

ся поверхность многослойных углеродных нанотру-

дуются возможности применения гибридных струк-

тур, в которых наночастицы или пленки алюминия

бок, луковичных структур и др. Данные, получен-

ные при исследовании взаимодействия различных

синтезируются на углеродной поверхности [9, 10].

Ключевым в этих задачах является вопрос о взаимо-

веществ с такими поверхностями, могут быть обоб-

щены. В этом смысле в настоящей работе локальную

действии алюминия и углерода на границе раздела

фаз, которое, в свою очередь, определяется струк-

структуру поверхности наночастиц будем называть

графеноподобной, не имея при этом в виду наличия

турой межфазной границы и энергией межатомных

взаимодействий. Для повышения адгезии на прак-

двумерной кристаллической решетки.

тике часто прибегают к химической модификации

Исследованию взаимодействия между алюмини-

поверхности наночастиц [11-14] либо к формирова-

ем и графеноподобной поверхностью углерода по-

священо множество экспериментальных и теорети-

^* E-mail: viktor.reshetnyak84@gmail.com

ческих работ. В работах [17-19] было изучено сма-

255

В. В. Решетняк, А. В. Аборкин

ЖЭТФ, том 157, вып. 2, 2020

чивание поверхности графита, пиролитического и

границе диапазона известных из литературы экспе-

стекловидного углерода расплавом алюминия. При

риментальных данных. При этом в работе [26], в ко-

этом капля расплава помещалась на поверхность уг-

торой используется эмпирическая формула, связы-

лерода, после чего измерялось значение контактного

вающая коэффициент поверхностного натяжения с

угла. Существенное влияние на результаты измере-

температурой плавления [27], было получено значе-

ний оказывает окисление поверхности капли [19]. Во

ние γ_Al = 1.032 Дж/м², что ближе к верхней границе

избежание этого эффекта измерения в работе [17]

диапазона.

выполнялись в вакууме, а в [18] — в атмосфере во-

Теоретически взаимодействие алюминия и гра-

дорода. Проведение эксперимента при этом услож-

феноподобной поверхности углерода исследовалось

няется неравновесностью системы: согласно резуль-

с использованием квантовомеханических методов

татам работ [17,18], значение контактного угла в хо-

расчета для кристаллических пленок [28, 29], клас-

де измерений меняется примерно от 140^◦ до 55^◦, что

теров [30-32] и отдельных атомов алюминия [33-35].

авторы связывают с формированием карбида Al₄C₃.

Заметим, что возможность непосредственного ис-

В работе [17] также отмечено влияние на резуль-

пользования результатов таких расчетов для оценки

тат шероховатости поверхности графита, наличия

энергии металлоуглеродной границы при конечных

на ней дефектов и примесей. Тем не менее результа-

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

ты работ [17,18] находятся в неплохом соответствии

ва алюминия вызывает вопросы.

друг с другом. Анализируя экспериментальные дан-

Учет конечной температуры и реальной струк-

ные, авторы пришли к выводу, что алюминий пло-

туры межфазной границы возможен при исполь-

хо смачивает поверхность графита: значение кон-

зовании методов молекулярной динамики (МД) и

тактного угла при этом составляет около 140^◦. Со

Монте-Карло (МК). При этом квантовомеханиче-

временем значение контактного угла уменьшается,

ское описание состояния границы, предполагающее

что, вероятно, связано с нарушением графеноподоб-

учет электронных степеней свободы, связано с суще-

ной структуры поверхности углерода и образовани-

ственными вычислительными трудностями. Поэто-

ем карбидной фазы.

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

Зная значения контактного угла θ и коэффи-

атомов не принимается во внимание, а межатомное

циента поверхностного натяжения γ_Al алюминия,

взаимодействие описывается аналитическими функ-

удельную энергию γ_i адгезии можно вычислить по

циями, параметры которых подбираются для наи-

формуле Юнга - Дюпре

лучшего описания свойств системы, известных из

экспериментов или расчетов ab initio. Анализ ли-

γ_i

= 1 + cosθ.

(1)

тературы позволяет сделать вывод, что для описа-

γ_Al

ния свойств межфазной границы в системе алюми-

Измерению коэффициента поверхностного натя-

ний-углерод межатомные взаимодействия изучены

жения алюминия посвящены работы [20-24]. Одна-

недостаточно.

ко из-за принципиальных экспериментальных труд-

Расчетам ab initio взаимодействия идеальных

ностей, связанных с чувствительностью результатов

кристаллических пленок алюминия и графена по-

к условиям эксперимента, разброс значений в этих

священы работы [28, 29]. В работе [29] также от-

работах оказывается весьма существенным: в рабо-

мечено, что атомы углерода и алюминия в таких

тах [23, 24] в зависимости от условий эксперимента

структурах образуют металлические связи, для опи-

были получены значения от 0.626 до 0.920 Дж/м²,

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

в то время как в работах [20-22] предлагаются зна-

или Ми, чем потенциал Леннарда-Джонса, однако

чения от 1.07 до 1.18 Дж/м². Некоторую определен-

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

ность вносят результаты теоретической работы [25],

потенциалов взаимодействия.

где с использованием метода молекулярной динами-

Исследованию механических свойств алюмомат-

ки были вычислены значения коэффициентов по-

ричных композитов на основе углеродных наноча-

верхностного натяжения для расплавов алюминия

стиц графеноподобной структуры (углеродных на-

и других металлов. Хотя результаты, полученные в

нотрубок и графена) методом МД посвящены рабо-

работе [25] для различных потенциалов взаимодей-

ты [36-42]. Правомерность использования выбран-

ствия между атомами алюминия, заметно отлича-

ных в этих работах параметров потенциалов меж-

ются друг от друга (диапазон примерно от 0.25 до

атомного взаимодействия Al-C вызывает вопро-

0.9 Дж/м²), вычисленные значения γ_Al не превы-

сы. В работах [36-41] был использован потенциал

шают 0.9 Дж/м², т. е. находятся ближе к нижней

Леннарда-Джонса, параметры которого были вы-

256

ЖЭТФ, том 157, вып. 2, 2020

Взаимодействие между атомами алюминия и углерода. . .

числены по правилу смеси Лоренца - Бертло. Потен-

энергии межфазной границы, образованной каплей

циал Леннарда-Джонса имеет вид

расплава алюминия на поверхности графита. Вы-

]

числения выполнялись методом МД. Взаимодей-

[(

)₁₂

(

)₆

ствие между атомами алюминия описывалось в рам-

u_ij = 4ε

(2)

rij

ках модели погруженного атома (embedded atom

model, EAM) с параметрами, предложенными в ра-

где u_ij — энергия взаимодействия пары атомов с ин-

боте [48] для расчета свойств алюминия в кристал-

дексами i и j, r_ij — расстояние между этими атома-

лическом и жидком состояниях. Ковалентное взаи-

ми, ε и σ — параметры.

модействие между атомами углерода описывалось

Согласно правилу Лоренца - Бертло, значения

потенциалом Терсоффа [49]. Для расчета свобод-

параметров для смеси компонент 1 и 2 определя-

ной энергии межфазной границы на начальном эта-

ются формулой

пе вычислялся коэффициент поверхностного натя-

жения расплава алюминия, а затем — контактный

ε₁₂ =

√ε₁ε₂,

угол, образованный каплей расплава на поверхности

(3)

σ₁₂ =

(σ₁ + σ₂).

графита. Значение свободной энергии определялось

путем подстановки полученных данных в формулу

Хорошо известно, что правомерность применения

(1). Для проверки модели использовались известные

параметров, вычисленных таким образом, в каж-

из литературы экспериментальные данные.

дом конкретном случае требует тщательной провер-

ки [43-45].

2. ОПИСАНИЕ

Авторы работы

[42] использовали потенциал

КВАНТОВОМЕХАНИЧЕСКОЙ МОДЕЛИ

Морзе, параметры для которого были взяты из ра-

ВЗАИМОДЕЙСТВИЯ КРИСТАЛЛИЧЕСКИХ

боты [46], где потенциал был использован для описа-

ПЛЕНОК ГРАФЕНА И АЛЮМИНИЯ

ния нековалентных взаимодействий между алюми-

нием и карбидом кремния. Заметим, что парамет-

Вычисления выполнялись с использованием про-

ризация этого потенциала была выполнена в рабо-

граммного пакета CP2K [50, 51]. Расчет электрон-

те [47] по данным расчетов ab initio для идеальных

ных состояний осуществлялся в псевдопотенциаль-

кристаллических пленок алюминия и карбида крем-

ном приближении, методом гауссовых орбиталей и

ния, а его применимость для описания взаимодейст-

плоских волн (GPW) [50,52]. При этом для представ-

вия алюминия с поверхностью углерода в sp²-сос-

ления волновых функций электронов используется

тоянии не доказана.

базисный набор гауссового типа, а для описания рас-

Настоящая работа посвящена теоретическому

пределения электронной плотности — вспомогатель-

исследованию межатомного взаимодействия Al-C

ный базис плоских волн. В настоящей работе для

на границе раздела фаз, между графеноподобной

атомов всех типов использовался дважды валент-

поверхностью углерода и алюминием в кристалли-

но-расщепленный базис с поляризационной функци-

ческом и жидком состояниях. В ходе исследования

ей DZVP [53] и сохраняющий норму псевдопотенци-

были привлечены различные теоретические мето-

ал Гоэдеккера - Тетера - Хуттера (GTH) [54,55], оп-

ды. Сначала с использованием теории функциона-

тимизированный для градиентно-обобщенного при-

ла плотности (density functional theory, DFT) была

ближения [56]. Взаимодействие валентных электро-

рассчитана поверхность потенциальной энергии для

нов вычислялось в рамках градиентно-обобщенного

взаимодействующих между собой идеальных крис-

приближения с использованием функционала Пер-

таллических пленок графена и алюминия, структу-

дью - Бурке - Эрнцерхофа (PBE) [57].

ра и относительное расположение которых соответ-

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

ствуют минимуму потенциальной энергии, а также

в пакете CP2K используется многосеточный метод.

сдвинутых друг относительно друга по нормали к

При этом точность вычисления характеризуется ко-

поверхности раздела фаз. По результатам расчетов

личеством сеток, а также параметрами их детализа-

были определены параметры потенциала Морзе, ко-

ции (параметры CUTOFF и REL_CUTOFF в про-

торый использовался нами для приближенного опи-

грамме CP2K). В настоящей работе было использо-

сания взаимодействия между атомами алюминия и

вано 4 сетки, а значения параметров, которые опре-

углерода.

делялись тестированием сходимости решения, со-

Затем вычисленные параметры потенциала вза-

ставили соответственно 350Ry и 60Ry. Для вычисле-

имодействия использовались при расчете свободной

ния интегралов по зоне Бриллюэна использовалась

257

ЖЭТФ, вып. 2

В. В. Решетняк, А. В. Аборкин

ЖЭТФ, том 157, вып. 2, 2020

сетка k-точек, построенная по алгоритму Монхорс-

ки a =

2.8519Å, параметры суперъячейки алюми-

та - Пака [58]. Число k-точек определялось путем ис-

ния размером 6 × 6 (a =

17.1114Å) и гексагональ-

следования сходимости решения и было разным для

ной суперъячейки графена 7 × 7 (a = 17.2767Å) ма-

различных задач.

ло отличаются друг от друга. При построении плен-

Для предварительной проверки DFT-модели бы-

ки использовались полученные в настоящей работе

ла выполнена серия тестовых расчетов, в кото-

теоретические значения параметров решетки и ко-

рых вычислялись значения параметров решетки и

ординат атомов. Построенная таким образом гете-

констант жесткости графита и алюминия. Расчет

роструктура содержит в ячейке 98 атомов углерода

свойств алюминия выполнялся для кубической эле-

и 180 атомов алюминия.

ментарной ячейки с использованием сетки k-точек

Аналогично, тетрагональная суперъячейка

размером 18 × 18 × 18. Расчет свойств графита вы-

Al(100) размером 6 × 6 (a =

17.1114Å) близка

полнялся для гексагональной элементарной ячейки.

по значению параметра решетки суперъячейке

Для учета нековалентного взаимодействия между

графена 7 × 4 (a =

17.2774Å, b = 17.0996Å). При

соседними слоями графита использовалась диспер-

построении суперъячейки графена использована

сионная поправка Гримме DFT-D3 [59]. Интегри-

орторомбическая ячейка, содержащая

атома

рование по зоне Бриллюэна выполнялось на сетке

углерода [60], параметры которой после релаксации

Монхорста - Пака размером 18 × 18 × 6. Значения

составили a =

2.4682Å, b

= 4.2749Å. Количество

остаточных сил после релаксации атомных позиций

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

и параметров решетки не превышали 5 · 10-3 эВ/Å.

венно 112 и 180. В расчетах использовалась сетка

Рассчитанные свойства алюминия и графита

Монхорста - Пака размером 3 × 3 × 1.

приведены в табл. 1 вместе с экспериментальными

данными работ [60-69] и результатами вычислений

других авторов [28,61,70-80]. Сравнение свидетель-

2.1. Модель взаимодействия пленки Al(111)

ствует о высокой точности модели, использованной

с поверхностью графена

в настоящей работе.

На начальном этапе была выполнена релаксация

При исследовании взаимодействия алюминия с

позиций атомов алюминия с неизменными парамет-

поверхностью графена были рассмотрены кристал-

рами ячейки. Удельная энергия ε_surf поверхности

лические пленки алюминия, вырезанные из крис-

была вычислена по формуле

талла поверхностями (111) и (100). Поверхности

алюминиевой и углеродной структур при этом рас-

(

)

N_slab

полагались на небольшом расстоянии параллельно

ε_surf =

E_slab -

E_bulk

(4)

2S_slab

N_bulk

друг другу. Толщина пленок алюминия составляла

5 атомных слоев.

где E_slab и N_slab (E_bulk и N_bulk) — энергия и ко-

Возникновение некогерентной границы раздела

личество атомов в пленке (в объемном кристалле),

фаз неизбежно ведет к деформации кристалличе-

S_slab — площадь одной из двух поверхностей пленки.

ских пленок. Вычисленная методом ab initio энергия

Полученное для Al(111) значение ε_surf

деформации может отличаться от соответствующе-

= 0.84 Дж/м² находится в разумном соответствии с

го значения, полученного с использованием эмпири-

результатами расчетов ab initio других авторов (см.

ческих потенциалов межатомного взаимодействия.

табл. 1).

В результате этого несоответствия при параметриза-

ции потенциала взаимодействия Al-C могут возни-

Затем вычислялась энергия взаимодействия пле-

кать ошибки. Ввиду хорошего соответствия вычис-

нок. При этом начальное расстояние от крайних ато-

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

мов до верхней и нижней границ ячейки выбиралось

кристаллов экспериментальным данным, характер-

равным 10Å.

ного как ab initio, так и эмпирическим расчетам

Поскольку значения констант жесткости графи-

[48, 49], мы предполагаем, что указанный эффект

та намного больше, чем у алюминия, при построе-

не играет принципиальной роли. Однако для сни-

нии гетероструктуры деформацией ячейки графена

жения вклада вероятных ошибок желательна мини-

пренебрегали, а параметры ячейки алюминия при-

мизация напряжений, возникающих на границе раз-

водились в соответствие параметрам ячейки графе-

дела фаз из-за структурного несоответствия алюми-

на. Обратные координаты атомов при этом сохра-

ниевой пленки и графена. Поскольку пленка Al(111)

нялись неизменными. Затем при неизменных значе-

имеет гексагональную ячейку с параметром решет-

ниях параметров решетки выполнялась релаксация

258

ЖЭТФ, том 157, вып. 2, 2020

Взаимодействие между атомами алюминия и углерода. . .

Таблица 1. Сравнение рассчитанных методом ab initio свойств алюминия и графита с экспериментальными дан-

ными и расчетами других авторов

Графит

Алюминий

Насто-

Расчеты

Насто-

Расчеты

ящая

других

Эксперимент ящая

других

Эксперимент

работа

авторов

работа

авторов

2.4395-2.4658

3.960-4.060

a,Å

2.4666

2.456 [60]

4.0331

4.0469 [60]

[70, 71]

[61, 72-74]

1056-1420

94.0-121.8

107.6-123.0

c₁₁, ГПа

1048.3

1211.3 [75]

117.2

[62-65]

[73, 76, 77]

[66-68]

140-180

54.5-69.6

60.4-70.8

c₁₂, ГПа

190.1

275.5 [75]

56.9

[62, 63, 65]

[73, 76, 77]

[66-68]

c₁₁+c₁₂, ГПа

1238.4

976-1283 [71]

174.1

183 [78]

11.6-37.4

28.3-30.9

c₄₄, ГПа

30.6

[73, 76-78]

[66-68]

74.7-80.2

78.1-80.1

B, ГПа

76.9

[28, 61, 72-74, 77, 78]

[61, 69]

ε_surf (111),

0.67-1.47

0.84

Дж/м²

[28, 73, 74, 77, 79, 80]

ε_surf (100),

0.86-1.63

0.95

Дж/м²

[28, 73, 74, 77, 79, 80]

ε_surf (110),

0.93-1.70

0.97

Дж/м²

[74, 79, 80]

положений атомов. Расчет энергии взаимодействия

горитму для системы Al(111)-графен. Возможная,

выполнялся по формуле

согласно работе [81], реконструкция поверхности

алюминия не учитывалась. Структура системы

ε_i =

[E_Al-C - (E_Al + E_C)] ,

(5)

алюминий-графен представлена на рис. 1б.

2S_slab

где E_Al-C энергия системы алюминий-углерод, E_Al

и E_C — энергии отдельно взятых алюминиевой и уг-

3. ОПИСАНИЕ МОДЕЛИ

леродной пленок.

ВЗАИМОДЕЙСТВИЯ РАСПЛАВА

Затем графен сдвигался относительно равновес-

АЛЮМИНИЯ С ПОВЕРХНОСТЬЮ

ного положения по оси z, перпендикулярной грани-

ГРАФИТА

це раздела фаз (рис. 1a). Расстояние Δz между по-

верхностями изменялось в диапазоне от 3.1 до 4.9Å

Во всех расчетах методом МД, выполненных в

с шагом 0.2Å. Значение энергии взаимодействия вы-

настоящей работе, интегрирование уравнений дви-

числялось по формуле (5) для каждой конфигура-

жения выполнялось с шагом Δt = 5 · 10-16 c. Тести-

ции.

рование свидетельствует о том, что в рассмотренном

диапазоне температур (от 1000 до 2000 К) для алю-

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

2.2. Модель взаимодействия пленки Al(100)

сервативности разностной схемы: при интегрирова-

с поверхностью графена

нии консервативных уравнений движения (nV E-ан-

Алгоритм

исследования

взаимодействия

самбль) значение полной механической энергии си-

Al(100)-графен аналогичен описанному выше ал- стемы флуктуирует относительно среднего значе-

259

В. В. Решетняк, А. В. Аборкин

ЖЭТФ, том 157, вып. 2, 2020

использованные для расчета параметров потенциа-

ла в работе [48], позволяют ожидать высокую точ-

ность при описании свойств алюминия как в жид-

кой, так и в кристаллической фазах.

Для расчета ковалентного взаимодействия меж-

ду атомами углерода был использован потенциал

Терсоффа [49]. Как и EAM, потенциал Терсоффа

учитывает влияние локального окружения на энер-

гию взаимодействия пары атомов. Однако, в отли-

чие от EAM, данный потенциал предназначен для

описания межатомных взаимодействий в атомных

системах с ковалентными связями [49, 86].

Взаимодействие между атомами алюминия и уг-

лерода описывалось парным потенциалом Морзе,

который имеет вид [87]

u_ij = D₀ {exp[-2α(r_ij - r₀)] -

Рис. 1. Периодическая ячейка для расчета ab initio энер-

- 2exp[-α(r_ij - r₀)]},

(6)

гии взаимодействия пленки алюминия с графеном в зави-

где D₀, α, r₀ — параметры потенциала.

симости от межплоскостного расстояния Δz: a — гекса-

гональная ячейка, поверхность Al(111); б — орторомбиче-

Потенциал наиболее пригоден для описания

ская ячейка, поверхность Al(100)

взаимодействия атомов в двухатомных молекулах

(для чего он и был изначально предложен), а его

применимость к описанию состояния конденсиро-

ния, которое, в свою очередь, остается неизменным

ванного вещества ограничена ввиду пренебрежения

при времени интегрирования вплоть до 10-8 с.

многочастичными эффектами. Тем не менее извест-

Для расчетов, предполагающих управление

но, что использование потенциала Морзе позволяет

температурой системы (nV T- и nPT-ансамбли), ис-

предсказывать некоторые свойства металлов в

пользовался термостат Нозе - Гувера. Вычисления,

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

предполагающие контроль давления (nP T -ан-

влетворительной точностью [88-92]. Несомненным

самбль), выполнялись с использованием алгоритма

плюсом при этом является сравнительная простота

работы [82].

потенциала: для выполнения вычислений требуется

Расчет средних по времени выполнялся на раз-

определить всего

параметра для пары атомов

личных временных интервалах в зависимости от

каждого сорта, в то время как при описании вза-

характерного времени задачи. Все результаты для

имодействий с учетом многочастичных эффектов

каждой из рассмотренных задач получены усредне-

в многокомпонентных системах число параметров

нием по 10 статистически независимым вычислени-

может достигать нескольких десятков.

ям.

Параметры потенциала вычислялись методом

наименьших квадратов с использованием результа-

тов расчетов ab initio для кристаллических пленок

3.1. Потенциалы межатомных

Al(111)-графен и Al(100)-графен, после чего ре-

взаимодействий

зультаты усреднялись.

Взаимодействие между атомами алюминия опи-

Необходимо обратить внимание на тот факт, что,

сывалось в рамках EAM с использованием парамет-

поскольку потенциал не учитывает влияния локаль-

ров, предложенных в работе [48]. Модель учитывает

ного окружения пары атомов на энергию взаимодей-

зависимость энергии взаимодействия пары атомов

ствия между ними, модель не может быть использо-

от их локального окружения и успешно применя-

вана для исследования таких эффектов как, к при-

ется при описании межатомного взаимодействия в

меру, влияние алюминия на термическую и меха-

металлах и сплавах [83, 84]. Подробное физическое

ническую стабильность поверхности графеноподоб-

обоснование возможности использования потенциа-

ного углерода. Неправомерно также обобщение мо-

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

дели для исследования взаимодействия алюминия

ном состоянии можно найти в работе [85]. Данные,

с углеродной поверхностью алмазоподобной струк-

260

ЖЭТФ, том 157, вып. 2, 2020

Взаимодействие между атомами алюминия и углерода. . .

туры, межатомного взаимодействия Al-C в карби-

3.2. Расчет коэффициента поверхностного

де Al₄C₃, примесных атомов, учета изменения хи-

натяжения расплава алюминия

мической активности атомов углерода, связанной с

Коэффициент поверхностного натяжения вычис-

неидеальностью структуры графена.

лялся для плоской пленки расплава алюминия, пе-

Радиус обрезки взаимодействия, r_cut, задавался

риодической в направлениях y и z. Структура плен-

равным 12Å. Значения параметров потенциала, вы-

ки задавалась с использованием следующей проце-

численные по методу наименьших квадратов, соста-

дуры. Вначале методом МД вычислялась равновес-

вили D₀ = 4.86 · 10-3 эВ, α = 1.29Å-1, r₀ = 4.15Å.

ная плотность n расплава алюминия при заданной

Для вычисления величины ошибки, возникаю-

температуре и нулевом давлении. Затем строилась

щей из-за обрезки потенциала взаимодействия, был

пленка алюминия, начальные положения атомов в

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

которой задавались соответственно позициям в иде-

хвостовой части потенциала U_cut [93], слегка моди-

альной кристаллической пленке Al(100), периодиче-

фицированный для случая взаимодействия атома с

ской в направлениях y и z, а плотность — равной

поверхностью. Потенциальная энергия, приходяща-

предварительно рассчитанной плотности расплава

яся на 1 атом вещества с парным потенциалом вза-

при заданной температуре. Размер расчетной облас-

имодействия, описывается выражением [94]

ти в каждом из направлений задавался равным 30

∫

размерам ГЦК-ячейки алюминия. Толщина пленки

U =

u(r)g(r) dr,

(7)

считалась равной 21 атомному слою. Полное коли-

чество атомов составило 37800. Расстояние от по-

верхностей пленки до соответствующих границ рас-

где n — плотность атомов вещества.

четной области задавалось равным 50Å; на грани-

Рассмотрим атом алюминия, находящийся над

цах ячейки, перпендикулярных оси x, использова-

поверхностью графита на расстоянии R₀. Графит

лись отражающие граничные условия. Объем меж-

разобьем с помощью сферических поверхностей,

ду границами расчетной области и поверхностями

центр которых совпадает с положением атома алю-

алюминиевой пленки оставлялся незаполненным. В

миния. Пересечение сфер с графитом образует сфе-

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

рические сегменты, радиус которых R = R₀/ cos φ,

Релаксация пленки выполнялась в два этапа.

φ — угол между нормалью к поверхности графи-

Сначала пленка стабилизировалась при температу-

та и отрезком, проведенным из центра сферы в лю-

ре 2000 К в течение 15 пс, а затем охлаждалась до

бую из точек окружности, образованной пересече-

нужной температуры, при которой релаксировалась

нием сферы с поверхностью. Если изменение угла

еще 15 пс. Были рассмотрены значения температур

dφ для любой пары соседних сфер мало, то каждый

в диапазоне от 1000 до 2000 К, включая границы

из атомов углерода в элементарном объеме между

интервала, с шагом 250 К. В ходе стабилизации рас-

ними расположен на расстоянии R от атома алюми-

плава контролировалась эволюция температуры и

ния и взаимодействует с ним с энергией u(R). Ради-

потенциальной энергии системы, а также значений

ус R ≥ r_cut предполагается достаточно большим для

параметров ориентационного порядка связей q₆ и q₁₂

пренебрежения корреляциями: g(r) = 1. Интегрируя

[95, 96]. Анализ эволюции параметров свидетельст-

по всем элементарным объемам, запишем выраже-

вует о достижении равновесного состояния системы.

ние для потенциальной энергии атома в виде

После стабилизации вычислялся коэффициент

поверхностного натяжения γ по формуле [97, 98]

∫

(

R₀

) 1 - cosφ

U = πnR³⁰

dφ.

(8)

∫

cosφ

cos³ φ

φ0

γ =

[P_N (x) - P_T (x)] dx.

(9)

x_l

Для оценки хвостовой части потенциала взаимо-

действия положим R₀ = 3.5Å (это значение при-

Здесь x_l, x_r — границы расчетной области, парал-

мерно соответствует межплоскостному расстоянию,

лельные поверхностям пленки, P_N и P_T — локаль-

вычисленному ab initio), a φ₀ = arccos(R₀/r_cut). Тог-

ные значения нормальной и тангенциальной компо-

да численное интегрирование выражения (8) дает

нент тензора напряжений. Для вычисления компо-

U_cut ∼ 10-7 эВ, что примерно в 10⁴ раз меньше ве-

нент тензора напряжений был применен подход ра-

личины D₀. Малость величины U_cut свидетельствует

боты [99], основанный на использовании вириально-

в пользу выбранного значения радиуса обрезки.

го уравнения состояния. Реализованный в програм-

261

В. В. Решетняк, А. В. Аборкин

ЖЭТФ, том 157, вып. 2, 2020

ме LAMMPS алгоритм расчета тензора напряже-

для вычисления U2∞(R₀), где φ₀ = 0, а расстоя-

ний для систем атомов, взаимодействующих посред-

ние от атома алюминия до поверхности R₀ изме-

ством многочастичных потенциалов, описан в рабо-

няется в диапазоне от 7.5 до 12Å, и выполнив ап-

те [100].

проксимацию по методу наименьших квадратов, мы

Для усреднения по времени был выбран интер-

вычислили параметры потенциала: ε_Al = 7.05 эВ,

вал 2 пс (4000 шагов), для получения окончатель-

σ_Al = 3.45Å. Заметим, что хотя ε_Al ≫ D₀, значение

ных значений n и γ для каждой температуры было

энергии взаимодействия атома с границей, U2∞, не

выполнено 10 статистически независимых вычисле-

превышает 3.5 · 10-4 эВ. Эта величина более чем в

ний.

10 раз меньше значения D₀, а значит, доминирую-

щая роль во взаимодействии с расплавом алюминия

принадлежит атомам верхнего графенового слоя.

3.3. Взаимодействие расплава алюминия с

Это позволяет предположить, что ошибки, внесен-

поверхностью графита (0001)

ные при использовании приближения среднего поля

При моделировании поверхности графита рас-

для описания взаимодействия алюминия со вторым

сматривалось движение только тех атомов углеро-

и более глубокими слоями графита, будут иметь вто-

да, которые относятся к верхнему атомному слою

рой порядок малости и не окажут принципиального

(графеновому листу). Взаимодействие атомов алю-

влияния на результаты вычислений.

миния и углерода с атомами второго и более глубо-

В начальный момент времени лист графена раз-

ких слоев учитывалось в приближении среднего по-

мером 21.99 × 25.27 нм²

(20000

атомов) и сфе-

ля. При этом считалось, что одна из границ расчет-

рическая наночастица алюминия диаметром около

ной области создает постоянное поле, действующее

10 нм (32581 атом), вырезанная из кристалла алю-

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

миния, располагались в орторомбической ячейке,

от сорта атома и его положения в ячейке. В настоя-

ограничивающей расчетную область, с параметра-

щей работе для этого поля был выбран вид (2), где

ми 21.99 × 25.27 × 20.0 нм³. На границах ячейки в

вместо расстояния r_ij между парой атомов задава-

направлениях x и y, параллельных плоскости гра-

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

фена, были заданы периодические граничные усло-

Значения параметров поля, ε_C и σ_C, для ато-

вия. Параметры a и b ячейки графена были вычис-

мов углерода были выбраны согласно известным

лены предварительно путем релаксации при темпе-

экспериментальным данным [60, 101]. Здесь и да-

ратуре 1000 К и нулевом давлении. Лист графе-

лее индексом «C» обозначены значения параметров

на располагался на расстоянии 3.348Å от нижней

для атомов углерода, а индексом Al — для алюми-

границы расчетной области и взаимодействовал с

ния. Параметр σ_C определялся в соответствии с из-

ней посредством потенциала (2) согласно изложен-

вестным для графита межплоскостным расстоянием

ному выше. Наночастица алюминия была помещена

rsep = 3.348Å: σC = rsep · 2-1/6 ≈ 2.983. Значение εC

в центр ячейки, при этом расстояния между атома-

было вычислено с использованием измеренного в ра-

ми алюминия и углерода, а также между атомами

боте [101] значения энергии взаимодействия плоско-

алюминия и границами расчетной области превы-

стей, γ_g ≈ 0.37 Дж/м². Подставляя значение поверх-

шали 4 нм, что намного больше радиусов обрезки

ностной плотности n_surf атомов углерода и перево-

всех использованных в работе потенциалов взаимо-

дя единицы, получим ε_C = γ_g/n_surf ≈ 6.03 · 10-2 эВ.

действия. Система стабилизировалась при темпера-

Для вычисления значений параметров ε_Al и σ_Al

туре 2000 К, а затем при 1000 К, согласно алгорит-

была рассчитана зависимость от расстояния энер-

му, описанному в разд. 3.2, для получения расплава

гии U2∞ взаимодействия атома алюминия со все-

алюминиевой пленки.

ми плоскостями графита, начиная со второй. Ми-

После стабилизации капля алюминия сдвигалась

нимальное возможное расстояние между атомами

по оси z и приводилась в контакт с поверхностью

алюминия и углерода в этом случае можно оценить

графена. Система релаксировалась при температу-

как r_sep + r₀ ≈ 7.5Å, что более чем в 2.5 раза боль-

ре 1000 К в течение 1 нс. Затем был выполнен расчет

ше среднего расстояния между атомами алюминия

траекторий атомов, который использовался для вы-

в расплаве, δ = n-1/3 ≈ 2.7Å. Следовательно, ра-

числения пространственного распределения плотно-

зумным является предположение о малости парных

сти капли алюминия.

корреляций между атомами алюминия и атомами

Плотность вычислялась путем разбиения рас-

углерода второго и более глубоких слоев графена:

четной области и нахождения среднего количества

g2∞(r_Al-C) ≈ 1. Воспользовавшись формулой (8)

атомов в каждом элементе объема. Для разбиения

262

ЖЭТФ, том 157, вып. 2, 2020

Взаимодействие между атомами алюминия и углерода. . .

использовались цилиндрические поверхности, шаг

E_int, Дж/м²

разбиения выбирался равным 2Å. Вычисления про-

водились через каждые 0.4 пс и усреднялись по вре-

мени счета, которое было выбрано равным 1 нс.

-0.05

В ходе расчетов постоянными задавались коли-

чество атомов, температура и объем. Для анализа

-0.10

структуры графенового листа вычислялось значе-

ние параметра трансляционного порядка [102]

-0.15

∑

ρ_k =

exp(ik · r_j ),

(10)

C j=1

-0.20

где k — любая линейная комбинация векторов об-

ратной решетки, rj — радиус-вектор j-го атома,

-0.25

2.5

3.0

3.5

4.0

4.5

5.0

N_C — полное число атомов углерода. В настоя-

z, Å

щей работе в качестве k был выбран вектор k_a =

= 2πb × n/|a × b|, где a и b — векторы решетки

Рис. 2. Зависимости энергии взаимодействия пленок алю-

графена, n — вектор нормали. При таком выборе

миния и графена от межплоскостного расстояния, рассчи-

k невозмущенной идеальной решетке графена соот-

танные методом ab initio (1,3) и с использованием потен-

ветствует значение ρ_k = 0.25, а в отсутствие даль-

циала Морзе (2,4): 1 и 2 — результаты вычислений для

него порядка ρ_k ∼ N-1/2C.

пленки Al(100); 3 и 4 — для Al(111)

4. РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ

величина которой 0.1422Å-2. Таким образом, при

повышении плотности упаковки абсолютное зна-

4.1. Взаимодействие графена и алюминия:

чение энергии взаимодействия уменьшается, что

расчет ab initio и параметризация

может быть связано с различной химической актив-

потенциала Морзе

ностью атомов алюминия, относящихся к разным

Зависимости потенциала взаимодействия графе-

поверхностям.

новой и алюминиевой пленок от межплоскостно-

Этот эффект, обусловленный влиянием локаль-

го расстояния, вычисленные разными методами

ного окружения атомов на их химическую актив-

(ab initio и с использованием потенциала Морзе),

ность, не может быть учтен в рамках модели пар-

представлены на рис. 2. Значения параметров потен-

ного взаимодействия. С данным эффектом связано

циала Морзе, рассчитанные по методу наименьших

качественное расхождение результатов, полученных

квадратов для поверхностей Al(100) и Al(111), а за-

в рамках различных приближений: с использовани-

тем усредненные составили D₀ = 4.86 · 10-3 эВ, α =

ем потенциала Морзе и DFT. Как уже было заме-

= 1.29Å-1, r₀ = 4.15Å. Радиус обрезки взаимодей-

чено, расчет ab initio предсказывает большее зна-

ствия r_c задавался равным 12Å. Отклонение значе-

чение энергии межфазной границы для Al(111), в

ний энергии взаимодействия пленок, вычисленных с

то время как при использовании потенциала Мор-

использованием потенциала Морзе, от результатов,

зе больше (меньше по модулю) оказывается энергия

полученных с использованием DFT, не превышает

для Al(100) (см. рис. 2). Действительно, поскольку

10 %.

плотность упаковки атомов алюминия больше для

Согласно результатам вычислений ab initio

Al(111), число соседних атомов алюминия у каж-

для поверхности Al(111), межплоскостное рас-

дого из атомов углерода для этой поверхности ока-

стояние между пленками после релаксации со-

зывается больше, чем для Al(100). Поэтому парный

ставило d₀

= 3.556Å, а энергия взаимодействия

потенциал Морзе, не учитывающий влияния локаль-

ε_i = -0.201 Дж/м². Для Al(100) значение межплос-

ного окружения на химическую активность ато-

костного расстояния мало отличается и составляет

мов, предсказывает более сильное взаимодействие

d₀

= 3.548Å, в то время как энергия взаимо-

для Al(111).

действия несколько ниже: εi = -0.232 Дж/м2.

Заметим также, что вычисленные в настоящей

При этом плотность упаковки атомов поверхности

работе значения межплоскостного расстояния и

Al(100) составляет примерно 0.1244Å-2, что мень-

энергии межфазной границы существенно отлича-

ше плотности упаковки атомов поверхности Al(111),

ются от полученных для Al(111) значений 3.18Å и

263

В. В. Решетняк, А. В. Аборкин

ЖЭТФ, том 157, вып. 2, 2020

n, кг/м³

Для описания границы раздела жидкость-газ

2400

широко используется следующая аппроксимация:

2(x - x₀)

2350

n_x =

(n_l + n_g)-

(n_l - n_g) th

(12)

— плотность жидкости и газа, а па-

где n_l и n_g

2300

раметры x₀ и ζ, определенные аппроксимацией ре-

зультатов МД-вычислений формулой (12), харак-

2250

теризуют толщину жидкой пленки и ширину гра-

ницы жидкость-вакуум. Ширина границы раздела

фаз «10-90», соответствующая расстоянию, на ко-

2200

тором плотность жидкости убывает от 90 % до 10 %

, равна ζ с точностью до постоян-

относительно n₀

2150

ного множителя: d = 2.1972ζ. Эта характеристика

1000

1200

1400

1600

1800

2000

поверхности может быть измерена эксперименталь-

T, K

но [108-110]. Температурная зависимость величины

Рис. 3. Температурная зависимость плотности расплава

d представлена на рис. 5.

алюминия: 1 — результат вычислений методом МД; 2 —

Температурная зависимость коэффициента по-

линейная аппроксимация

верхностного натяжения расплава алюминия, вы-

численная в настоящей работе, представлена на

рис. 6. Видно, что она близка к линейной и моно-

-0.11 Дж/м² [28]. Вероятно, расхождение обуслов-

тонно убывает в широком диапазоне температур:

лено использованием в расчетах [28] ячейки мало-

γ = aT + b. Значения коэффициента поверхност-

го размера, не обеспечивающей отсутствия напря-

ного натяжения γ_ml, соответствующего температу-

жений в системе.

ре плавления, а также углового коэффициента a

представлены в табл. 2. Из табл. 2 видно, что зна-

чение γ_ml попадает в диапазон экспериментальных

4.2. Свойства расплава алюминия

данных [21-24], находясь при этом ближе к нижней

Температурная зависимость плотности расплава

границе диапазона. Рассчитанное нами значение γ_ml

алюминия, вычисленная в настоящей работе мето-

также неплохо согласуется с теоретическими резуль-

дом МД, а также линейная аппроксимация резуль-

татами работы [25], полученными с использованием

татов представлены на рис. 3.

потенциалов работ [111,112].

Коэффициент объемного теплового расширения

Значение коэффициента a = -1.54·10-4 Дж/м²·К

связан с плотностью выражением

несколько отличается от экспериментальных

данных

[23],

где было получено значение

)

⁽∂n

−1.85 · 10-4 Дж/м²·К, но хорошо согласуется с

β =-

(11)

∂Tp=0

данными работы

[107], согласно которым γ

= 0.868 - 1.52 · 10-4(T - T_melt), где T_melt = 933 К —

Теоретические результаты, полученные для распла-

температура плавления алюминия.

ва алюминия в настоящей работе методом МД, при-

водятся в табл. 2 в сравнении с известными из лите-

4.3. Взаимодействие капли расплава

ратуры экспериментальными данными и результа-

алюминия с поверхностью графита

тами вычислений других авторов. Из табл. 2 видно,

что рассчитанные нами значения плотности распла-

В ходе предварительной релаксации капля алю-

ва и коэффициента теплового расширения находят-

миния была полностью расплавлена. Об этом свиде-

ся в хорошем соответствии с литературными данны-

тельствуют вычисленные значения параметров ло-

ми. Далее в работе также используется безразмер-

кального порядка q₆ и q₁₂ для алюминия. В началь-

ная плотность расплава ñ = n/n₀, где n₀ — теоре-

ный момент значения этих параметров, усредненные

тическое значение плотности в объеме расплава при

по всем атомам алюминия, были равны q₆ = 0.541

заданной температуре.

и q₁₂ = 0.579, что соответствует значениям для иде-

Структура и распределение плотности в алюми-

альной ГЦК-решетки 0.575 и 0.600 [96]. Таким обра-

ниевой пленке представлены на рис. 4.

зом, значения параметров локального порядка для

264

ЖЭТФ, том 157, вып. 2, 2020

Взаимодействие между атомами алюминия и углерода. . .

Таблица 2. Сравнение вычисленных методом МД свойств расплава алюминия с экспериментальными данными и

расчетами других авторов

МД,

настоящая

расчеты других

Эксперимент

работа

авторов

ρ_ml, кг/м³ (933 К)

2409

2338-2397 [48,103,104]

2368 [105]

β, К-1

0.98 · 10-4

0.76 · 10-4 [104, 106]^a

1.2 · 10-4 [105]^b

γ_ml, Дж/м² (933 К)

0.73

0.25-0.9 [25]

0.626-1.18 [20-24, 107]

-1.85 · 10-4 [23]

a, Дж/м²·К

-1.54 · 10-4

-1.52 · 10-4 [107]

◦

136.4

140^◦ [17,18]

137.5^◦c

Примечание. Значения β вычислены^a с использованием представленных в цитируемых работах значений

плотности при разных температурах;^b путем оцифровки представленного в цитируемой работе графика тем-

пературной зависимости плотности;^c с использованием толменовской поправки для поверхностного натяже-

ния капли

1.2

1.0

0.8

0.6

0.4

0.2

22.1

x, Å

Рис. 4. Структура пленки расплава алюминия при T = 1000 K: a — расположение атомов в расчетной ячейке; б —

распределение плотности по координате: 1 — результаты вычислений методом МД; 2 — аппроксимация выражением (12)

вырезанного из кристалла кластера несколько ниже

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

соответствующих значений для идеального кристал-

вали флуктуации относительно средних значений,

ла. Это объясняется наличием поверхности у клас-

которые не менялись со временем.

тера и различием значений параметров локального

Графен сохранял периодическую структуру на

порядка для атомов, расположенных на поверхно-

протяжении всего расчета. Об этом свидетельствует

сти и в объеме кристалла. В результате отжига при

тот факт, что значение параметра трансляционного

температуре 2000 К значения параметров убывали

порядка, определенного согласно выражению (10),

со временем, достигнув равновесных значений q₆ =

на протяжении всего расчета флуктуировало около

= 0.346 и q₁₂ = 0.287 спустя примерно 3 пс. Затем в

среднего, 〈ρ_k〉 ≈ 0.246, и мало отличалось от зна-

265

В. В. Решетняк, А. В. Аборкин

ЖЭТФ, том 157, вып. 2, 2020

d, Å

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

дований методом МД, ограничен. Реакция карбидо-

образования — активационный процесс, характер-

ное время которого связано с энергией активации и

температурой уравнением Арениуса. Для описания

такого процесса методом МД необходимо рассмот-

рение временного интервала, намного большего, чем

характерное время реакции. В настоящей работе вы-

числения были выполнены на временном интервале

2 · 10-9 с, что, по-видимому, не достаточно для мо-

делирования данной химической реакции. Следует

заметить, что протекание реакции в эксперименте

может быть ускорено наличием дефектов структу-

ры графеноподобной поверхности [16, 114]. Данный

эффект не был учтен в моделировании. Кроме того,

1000

1200

1400

1600

1800

2000

в набор данных, предложенный в настоящей рабо-

T, K

те при параметризации межатомного потенциала, не

были включены ни равновесные карбидные струк-

Рис.

5. Температурная зависимость ширины границы

туры, ни неравновесные активированные состояния,

жидкость-вакуум для расплава алюминия: 1 — результа-

соответствующие реакции карбидообразования. По-

ты, полученные при обработке рассчитанных методом МД

этому обобщение потенциала для описания как ука-

траекторий атомов с использованием формулы (12); 2 —

занной химической реакции, так и равновесного со-

линейная аппроксимация: d = 5.27 · 10-3T - 0.41

стояния карбида алюминия не представляется воз-

можным.

, Дж/м²

После релаксации структуры капля была приве-

0.72

дена в контакт с поверхностью графена. Расстояние

z_cm от поверхности подложки до центра масс капли

контролировалось в ходе вычислений. График вре-

0.68

менной зависимости z_cm(t) представлен на рис. 7,

где в качестве начального принят момент приведе-

ния капли в контакт с подложкой. Видно, что центр

0.64

масс капли вначале совершает затухающие колеба-

ния с периодом τ ≈ 100 пс. Возникновение колеба-

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

0.60

веденная в контакт с поверхностью графита капля

изначально была несколько смещена относительно

положения равновесия (расстояние от центра масс

0.56

капли до поверхности графита было несколько боль-

1000

1200

1400

1600

1800

2000

ше равновесного). В результате взаимодействия с

T, K

поверхностью графита капля приобретает отличную

от нуля среднюю скорость, направленную к поверх-

Рис. 6. Температурная зависимость коэффициента поверх-

ности. Затем, в результате столкновения с поверх-

ностного натяжения: 1 — результаты вычислений методом

МД; 2 — линейная аппроксимация: γ = -1.54·10-4T +0.87

ностью графита, в капле возникает волна сжатия,

которая движется в обратном направлении. Таким

образом, причина колебательного движения центра

чения 0.25, вычисленного для идеальной решетки.

масс — распространение упругой волны и ее отраже-

Этот результат противоречит экспериментальным

ние от границ капли. Из-за вязкости расплава энер-

данным [15, 113-118], согласно которым уже при

гия колебательного движения со временем перехо-

температуре 400-450^◦C углеродные наноструктуры

дит в тепловую, которая, в свою очередь, отводится

могут вступать в химическую реакцию с алюминие-

термостатом. Примерно за 300 пс колебания затуха-

вой матрицей с образованием карбидной фазы. Ве-

ют, после чего центр масс движется стохастически,

роятно, указанное расхождение связано с тем фак-

а амплитуда смещений относительно среднего при-

266

ЖЭТФ, том 157, вып. 2, 2020

Взаимодействие между атомами алюминия и углерода. . .

z_cm, Å

z_cm

500

1000

1500

2000

t, пс

Рис. 7. Релаксация капли на поверхности графита: a — атомы алюминия и углерода, для которых вычислялись траекто-

рии; б — временная зависимость положения центра масс капли относительно поверхности графита

мерно соответствует среднему межатомному рассто-

ласти, включающей первые ближайшие к зоне кон-

янию δ = n-1/3.

такта 10 значений z_i, для которых хотя бы в одной

Характерное время τ смещения центра масс (пе-

из точек расчетной области с координатами (r, z_i)

риод колебаний) определяется упругими свойствами

значение плотности приближается к равновесной

расплава и размером капли. Поскольку амплитуда

плотности расплава: ñ(r, z_i) ≥ 0.98. Значение кон-

стохастических смещений центра масс не затухает,

тактного угла было вычислено аналитически путем

время моделирования, необходимое для вычисления

дифференцирования зависимости r(z) в точке кон-

среднего значения контактного угла, должно быть

такта, в качестве которой была выбрана точка с ми-

существенно больше τ. В настоящей работе усред-

нимальным z = z₀. Окончательное выражение для

нение проводилось за время 1 нс, что примерно со-

контактного угла имеет вид

ответствует 10τ.

θ = 180 - arctg

(13)

При вычислении зависимости плотности алюми-

2az₀ + b

ния от координат форма капли считалась осесим-

Пространственное распределение плотности алюми-

метричной, а разбиение расчетной области в коор-

ния, а также результаты аппроксимации представ-

динатах r-z выполнялось с использованием цилинд-

лены на рис. 8.

рических поверхностей, шаг по осям r и z задавался

Вычисленное значение контактного угла θ соста-

равным 2Å.

вило 136.4^◦, что неплохо согласуется с эксперимен-

Для определения положения и ширины поверх-

том [17,18] (см. табл. 2). Существующее несоответ-

ности капли на границе с вакуумом использовалась

ствие результатов вычислений и экспериментальных

формула (12). Анализ выполнялся для каждого ин-

данных может быть связано с приближенным описа-

тервала [z, z+dz], где радиальная зависимость плот-

нием межатомных взаимодействий в настоящей ра-

ности аппроксимировалась выражением (12). По-

боте для пар атомов как Al-C, так и Al-Al. Кро-

ложение границы считалось равным вычисленному

ме того, расхождение может быть связано с измене-

значению x₀, а ширина «10-90» d = 2.1792ζ. Усред-

нием поверхностного натяжения капли жидкости с

ненное по z значение ширины границы расплав-

уменьшением ее размера. Как известно (см. работу

вакуум составило 7.1Å. Это значение несколько

[97] и цитируемую в ней литературу), поверхностное

больше вычисленного для плоской поверхности при

натяжение сферической капли радиуса R равно

T = 1000 К, однако намного меньше радиуса капли,

(

)

δ_∞

величина которого составляет около 50Å.

γ(R) = γ_∞

(14)

Анализ положения границы свидетельствует о

том, что вблизи зоны контакта зависимость r(z)

где γ_∞ — значение коэффициента поверхностного

хорошо аппроксимируется квадратичной функцией:

натяжения для поверхности с бесконечным радиу-

r = az²+bz+c. Аппроксимация выполнялась для об-

сом кривизны, а толменовская длина δ_∞ может быть

267

В. В. Решетняк, А. В. Аборкин

ЖЭТФ, том 157, вып. 2, 2020

z, Å

1.12

100

0.96

0.80

0.64

0.48

0.32

0.16

r, Å

Рис. 8. Пространственное распределение плотности расплава алюминия в переменных r-z (а) и положение границы

жидкость-вакуум (б): 1 — результаты аппроксимации распределения плотности с использованием выражения (12); 2 —

аппроксимация положения границы полиномом второго порядка; 3 — касательная, проведенная к аппроксимирующей

кривой в точке контакта

рассчитана как половина ширины границы раздела

5. ВЫВОДЫ

фаз: δ_∞ = d/2. Используя определенную ранее тем-

пературную зависимость ширины d, получим δ_∞ =

Предложена модель взаимодействия между ато-

= 2.64·10-3T -0.21. Температуре 1000 К соответст-

мами алюминия, находящегося в конденсированном

вует значение δ_∞ ≈ 2.43Å.

состоянии, и углерода, локальная структура поверх-

Подставляя значения радиуса капли, коэффици-

ности которого близка к структуре графена. Модель

ента поверхностного натяжения γ_∞ и толменовской

основана на использовании потенциала Морзе для

длины δ_∞, получим для капли γ_drop = 0.68 Дж/м².

пар атомов Al-C. Параметризация потенциала осу-

В экспериментах [17, 18] рассматривалась капля,

ществлялась по данным расчетов ab initio, выпол-

диаметр которой достаточно велик по сравнению с

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

δ_∞, чтобы зависимостью (14) можно было прене-

и графена с идеальной кристаллической структу-

бречь. Пересчет значения контактного угла для кап-

рой. Потенциал был использован при вычислении

ли бесконечного радиуса с учетом данного эффекта

взаимодействия капли расплава алюминия с поверх-

в соответствии с выражением (1) можно выполнить

ностью графита методом МД. Кроме того, метод

по формуле

МД был применен для исследования плоской по-

верхности расплава чистого алюминия.

[

(

)

]

δ_∞

Результаты настоящей работы сильно отличают-

θ_∞ = arccos (1 + cosθ)

-1 .

(15)

ся от данных, полученных в теоретических исследо-

ваниях другими авторами [28,36-42]. Так, результа-

ты вычислений ab initio, выполненных в настоящей

Однако полученная при пересчете поправка доста-

работе для пленки Al(111), существенно отличают-

точно мала: значение θ_∞ = 137.5^◦ отличается от θ =

ся от данных работы [28]: значение рассчитанной на-

= 136.4^◦ менее чем на 1 %, а от экспериментально

ми потенциальной энергии границы раздела фаз при

измеренной величины θ = 140^◦ — менее чем на 2 %.

нулевой температуре составило -0.201 Дж/Å2, в то

Расчет удельной свободной энергии межфазной

время как в [28] получено значение -0.11 Дж/Å2, а

границы в системе алюминий-углерод по форму-

межплоскостное расстояние, вычисленное в настоя-

ле (1) с использованием вычисленного значения

щей работе, составило 3.556Å, в то время как в [28]

γ_drop = 0.68 Дж/м² дает значение γ_i = 0.18 Дж/м².

было получено значение 3.18Å.

268

ЖЭТФ, том 157, вып. 2, 2020

Взаимодействие между атомами алюминия и углерода. . .

Сравнение вычисленных в настоящей работе зна-

К. В. Кремлев, А. М. Объедков, Н. М. Семенов и

чений параметров потенциала Морзе (глубины по-

др., Письма в ЖТФ 44(19), 24 (2018).

тенциальной ямы D₀ = 4.86·10-3 эВ и ее положения

10.

A. V. Aborkin, K. V. Kremlev, A. M. Obiedkov et

r₀ = 4.15Å) с соответствующими параметрами по-

al., J. Phys.: Conf. Ser. 1164, 012020 (2019).

тенциала Леннарда-Джонса (глубина ямы ε_Al-C =

11.

J. Nie, C. Jia, N. Shi et al., Int. J. Miner. Metall.

= 3.5 · 10-2 эВ и ее положение σ · 21/6 ≈ 3.4Å),

Mater. 18, 695 (2011).

рассчитанными в работах [36-42] с использованием

правила смеси Лоренца - Бертло, также показало

12.

K. P. So, I. H. Lee, D. L. Duong et al., Acta Mater.

существенные различия. В то же время результаты

59(9), 3313 (2011).

вычислений методом МД совпадают с эксперимен-

13.

J. Wang, Z. Li, G. Fan et al., Scripta Mater. 66, 594

тальными данными с высокой точностью: рассчи-

(2012).

танное в настоящей работе значение контактного

угла при взаимодействии капли расплава с поверх-

14.

S. W. Ip, R. Sridhar, J. M. Toguri et al., Mater. Sci.

Eng. A 244, 31 (1998).

ностью графита составило 136.4^◦, в то время как в

экспериментах [17-19] было получено 140^◦. Коррек-

15.

J. G. Park, D. H. Keum, and Y. H. Lee, Carbon 95,

тировка вычисленного для капли малого размера

690 (2015).

значения контактного угла с использованием зави-

16.

W. Zhou, T. Yamaguchi, K. Kikuchi et al., Acta

симости коэффициента поверхностного натяжения

Mater. 125, 369 (2017).

от радиуса кривизны поверхности дает 137.5^◦ и поз-

воляет незначительно улучшить соответствие тео-

17.

K. Landry, S. Kalogeropoulou, and N. Eustathopou-

рии и эксперимента. Вычисленное таким образом

los, Mater. Sci. Eng. A 254, 99 (1998).

значение контактного угла совпадает с эксперимен-

18.

S.-I. Oh, J.-Y. Lim, Y.-C. Kim et al., J. Alloys Comp.

тальным в пределах 2 %. Значение удельной свобод-

542, 111 (2012).

ной энергии межфазной границы, рассчитанное по

19.

K. Landry, S. Kalogeropoulou, N. Eustathopoulos et

формуле Юнга - Дюпре, составило γ_i = 0.18 Дж/м².

al., Scripta Mater. 34(6), 841 (1996).

20.

W. R. Tyson and W. A. Miller, Surf. Sci. 62, 267

Финансирование. Работа выполнена при

(1977).

поддержке Российского научного фонда (проект

21.

C. Garcia-Cordovilla, E. Louis, and A. Pamies, J.

№16-12-10424-П).

Mater. Sci. 21, 2787 (1986).

22.

Х. Х. Калажоков, З. Х. Калажоков, Х. Б. Хоконов,

ЖТФ 73(2), 141 (2003).

ЛИТЕРАТУРА

23.

J. M. Molina, R. Voytovych, E. Louis et al., Int. J.

1. S. R. Bakshi, D. Lahiri, and A. Agarwal, Int. Mater.

Adhesion Adhesives 27, 394 (2007).

Rev. 55, 41 (2010).

24.

I. F. Bainbridge and J. A. Taylor, Metall. Mat. Trans.

2. H. G. P. Kumar and M. A. Xavior, Proc. Eng. 97,

A 44, 3901 (2013).

10331 (2014).

25.

E. B. Webb and G. S. Grest, Phys. Rev. Lett. 86,

2066 (2001).

3. S. C. Tjong, Mater. Sci. Eng. Rep. 74(10), 281 (2013).

26.

K. K. Nanda, Phys. Lett. A 376(19), 1647 (2012).

4. N. Saheb, Z. Iqbal, A. Khalil et al., J. Nanomater.

2012, 983470 (2012).

27.

B. Dayal, Nature 169, 1010 (1952).

5. H. J. Choi, G. B. Kwon, G. Y. Lee et al., Scripta

28.

Y. Qi, L. G. Hector, N. Ooi et al., Surf. Sci. 581, 155

(2005).

Mater. 59, 360 (2008).

29.

W. Lee, S. Jang, M. J. Kim et al., Nanotechnology

6. S. R. Bakshi and A. Agarwal, Carbon 49, 533 (2011).

19, 285701 (2008).

7. I. A. Evdokimov, S. A. Perfilov, A. A. Pozdnyakov et

30.

D.-H. Lim, A. S. Negreira, and J. Wilcox, J. Phys.

al., Inorg. Mater. Appl. Res. 9, 472 (2018).

Chem. C 115, 8961 (2011).

8. I. A. Evdokimov, T. A. Chernyshova, G. I. Pivovarov

31.

I. Garg, H. Sharma, K. Dharamvir et al., J. Phys.

et al., Inorg. Mater. Appl. Res. 5, 255 (2014).

Chem. C 114, 18762 (2010).

269

В. В. Решетняк, А. В. Аборкин

ЖЭТФ, том 157, вып. 2, 2020

32.

I. Moullet, Surf. Sci. 331-333, 697 (1995).

54.

S. Goedecker, M. Teter, and J. Hutter, Phys. Rev.

B 54, 1703 (1996).

33.

A. Ishii, M. Yamamoto, H. Asano et al., J. Phys.:

Conf. Ser. 100, 052087 (2008).

55.

C. Hartwigsen, S. Goedecker, and J. Hutter, Phys.

Rev. B 58, 3641 (1998).

34.

M. T. Baei, Fullerenes, Nanotubes and Carbon Nano-

struct. 20, 681 (2012).

56.

M. Krack, Theor. Chem. Accounts 114, 145 (2005).

35.

M. Hamadanian and F. K. Fotooh, Comput. Mater.

57.

J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev.

Sci. 82, 497 (2014).

Lett. 77, 3865 (1996).

36.

D. M. Park, J. H. Kim, S. J. Lee et al., J. Mech. Sci.

58.

H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13,

Technol. 32, 5845 (2018).

5188 (1976).

37.

A. Mokhalingam, D. Kumar, and A. Srivastava,

59.

S. Grimme, J. Antony, S. Ehrlich et al., J. Chem.

Materials Today: Proc. 4(2A), 3952 (2017).

Phys. 132, 154104 (2010).

60.

R. W. G. Wyckoff, Crystal Structures, Vol. 1, Intersci.

38.

B. K. Choi, G. H. Yoon, and S. Lee, Composites B 91,

119 (2016).

Publ., New York (1963).

61.

R. Gaudoin and W. M. C. Foulkes, Phys. Rev. B 66,

39.

N. Silvestre, B. Faria, and J. N. Canongia Lopes,

Composites Sci. Technol. 90, 16 (2014).

052104 (2002).

62.

O. L. Blakslee, D. G. Proctor, E. J. Seldin et al., J.

40.

S. Xiao and W. Hou, Int. J. Multiscale Comput. Eng.

Appl. Phys. 41, 3373 (1970).

5, 447 (2007).

63.

E. J. Seldin and C. W. Nezbeda, J. Appl. Phys. 41,

41.

H.-Y. Song and X.-W. Zha, Comput. Mater. Sci. 49,

3389 (1970).

899 (2010).

64.

R. Nicklow, N. Wakabayashi, and H. G. Smith, Phys.

42.

J. Xiang, L. Xie, S. A. Meguid et al., Comput. Mater.

Rev. B 5, 4951 (1972).

Sci. 128, 359 (2017).

65.

A. Bosak, M. Krisch, M. Mohr et al., Phys. Rev. B 75,

43.

D. Boda and D. Henderson, Mol. Phys. 106, 2367

153408 (2007).

(2008).

66.

P. M. Sutton, Phys Rev. 91, 816 (1953).

44.

M. Rouha and I. Nezbeda, Fluid Phase Equil. 277,

42 (2009).

67.

J. Vallin, M. Mongy, K. Salama et al., J. Appl. Phys.

35, 1825 (1964).

45.

J. Delhommelle and P. Millié, Mol. Phys. 99, 619

(2001).

68.

J. F. Thomas, Phys. Rev. 175, 955 (1968).

46.

C. R. Dandekar and Y. C. Shin, Composites A 42,

69.

L. Gerward, J. Phys. Chem. Sol. 46, 925 (1985).

355 (2011).

70.

I. V. Lebedeva, A. S. Minkin, A. M. Popov et al.,

47.

H. Zhao and N. Chen, Inverse Probl. 24, 035019

Physica E 108, 326 (2019).

(2008).

71.

N. Mounet and N. Marzari, Phys. Rev. B 71, 205214

48.

H. W. Sheng, M. J. Kramer, A. Cadien et al., Phys.

(2005).

Rev. B 83, 134118 (2011).

72.

R. Gaudoin, W. M. C. Foulkes, and G. Rajagopal, J.

49.

J. Tersoff, Phys. Rev. B 39, 5566 (1989).

Phys.: Condens. Matter. 14, 8787 (2002).

50.

J. VandeVondele, M. Krack, F. Mohamed et al., Com-

73.

P. Jacobs, Y. F. Zhukovskii, Y. Mastrikov et al.,

puter Phys. Comm. 167, 103 (2005).

Comput. Model. New Technol. 6(1), 7 (2002).

51.

J. Hutter, M. Iannuzzi, F. Schiffmann et al., Comput.

74.

N. E. Singh-Miller and N. Marzari, Phys. Rev. B 80,

Mol. Sci. 4, 15 (2014).

235407 (2009).

52.

B. G. Lippert, J. Hutter, and M. Parrinello, Mol.

75.

K. H. Michel and B. Verberck, Phys. Stat. Sol. (b)

Phys. 92, 477 (1997).

245, 2177 (2008).

53.

J. VandeVondele and J. Hutter, J. Chem. Phys. 127,

76.

C. Woodward, D. R. Trinkle, L. G. Hector et al.,

114105 (2007).

Phys. Rev. Lett. 100, 045507 (2008).

270

ЖЭТФ, том 157, вып. 2, 2020

Взаимодействие между атомами алюминия и углерода. . .

77.

Y. Qi and L. G. Hector, Phys. Rev. B 69, 235401

98.

J. S. Rowlinson and B. Widom, Molecular Theory of

(2004).

Capillarity, Courier Corporation (2013).

78.

W. Li and T. Wang, J. Phys.: Condens. Matter. 10,

99.

E. Salomons and M. Mareschal, J. Phys.: Condens.

9889 (1998).

Matter 3, 3645 (1991).

79.

J. Schöchlin, K. P. Bohnen, and K. M. Ho, Surf. Sci.

100.

A. P. Thompson, S. J. Plimpton, and W. Mattson, J.

324, 113 (1995).

Chem. Phys. 131, 154107 (2009).

80.

C. Fiolhais, L. M. Almeida, and C. Henriques, Progr.

101.

W. Wang, S. Dai, X. Li et al., Nature Comm. 6, 7853

Surf. Sci. 74, 209 (2003).

(2015).

81.

R. Ludeke and G. Landgren, Phys. Rev. Lett. 47, 875

102.

D. C. Rapaport, The Art of Molecular Dynamics Si-

(1981).

mulation, Cambridge Univ. Press (2004).

103.

M. I. Mendelev, M. J. Kramer, C. A. Becker et al.,

82.

W. Shinoda, M. Shiga, and M. Mikami, Phys. Rev.

B 69, 134103 (2004).

Phil. Mag. 88, 1723 (2008).

104.

Д. К. Белащенко, УФН 183, 1281 (2013).

83.

J. Lao, T. M. Naghdi, D. Pinisetty et al., JOM 65(2),

175 (2013).

105.

G. E. Totten and D. S. MacKenzie, Handbook of Alu-

minum, Vol. 1, CRC Press, Boca Raton, USA (2003).

84.

P. Brommer, A. Kiselev, D. Schopf et al., Modell.

Simul. Mater. Sci. Eng. 23, 074002 (2015).

106.

Д. К. Белащенко, А. В. Воротягин, Б. Р. Гельчинс-

кий, ТВТ 49, 676 (2011).

85.

M. S. Daw and M. I. Baskes, Phys. Rev. B 29, 6443

(1984).

107.

J. E. Hatch, Aluminum: Properties and Physical Me-

tallurgy, ASM Internat. (1984).

86.

J. Tersoff, Phys. Rev. B 37, 6991 (1988).

108.

M. L. Schlossman, Current Opinion in Colloid & In-

87.

P. M. Morse, Phys. Rev. 34, 57 (1929).

terface Sci. 7, 235 (2002).

88.

R. C. Lincoln, K. M. Koliwad, and P. B. Ghate, Phys.

109.

S. J. Roser, R. Felici, and A. Eaglesham, Langmuir

Rev. 157, 463 (1967).

10, 3853 (1994).

89.

L. A. Girifalco and V. G. Weizer, Phys. Rev. 114,

110.

J. Lekner, Theory of Reflection: Reflection and Trans-

687 (1959).

mission of Electromagnetic, Particle and Acoustic

Waves, Springer, Berlin (2016).

90.

T. H. K. Barron and C. Domb, Proc. Roy. Soc.

London A 227, 447 (1955).

111.

A. Voter, S. Chen, R. Siegel et al., in MRS Symposia

Proc., Materials Research Society Pittsburgh, PA

91.

D. Weaire, M. F. Ashby, J. Logan et al., Acta Metall.

(1987), p. 175.

19, 779 (1971).

112.

Y. Mishin, D. Farkas, M. J. Mehl et al., Phys. Rev.

92.

M. Doyama and J. S. Koehler, Acta Metall. 24, 871

B 59, 3393 (1999).

(1976).

113.

А. В. Аборкин, К. С. Хорьков, А. М. Объедков и

93.

D. Frenkel and B. Smit, Understanding Molecular Si-

др., Письма в ЖТФ 45(2), 22 (2019).

mulation: from Algorithms to Applications, Elsevier,

Amsterdam (2001).

114.

B. Chen, J. Shen, X. Ye et al., Carbon 114, 198

(2017).

94.

Р. Балеску, Равновесная и неравновесная статис-

тическая механика, т. 1, Мир, Москва (1978).

115.

M. H. Vidal-Sétif, M. Lancin, C. Marhic et al., Mater.

Sci. Eng. A 272(2), 321 (1999).

95.

P. J. Steinhardt, D. R. Nelson, and M. Ronchetti,

Phys. Rev. B 28, 784 (1983).

116.

D. Poirier, R. Gauvin, and R. A. L. Drew, Composites

A 40, 1482 (2009).

96.

W. Mickel, S. C. Kapfer, G. E. Schröder-Turk et al.,

J. Chem. Phys. 138, 044501 (2013).

117.

И. А. Евдокимов, Дисс

канд. техн. наук, ВлГУ,

Владимир (2013).

97.

С. Оно, С. Кондо, Молекулярная теория поверх-

ностного натяжения в жидкостях, Рипол Клас-

118.

А. В. Аборкин, М. И. Алымов, А. В. Собольков и

сик, Москва (2013).

др., Металлы 2018(4), 27 (2018).

271