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

МОДЕЛЬ ДВУХМЕРНОГО ГРАФИТОВОГО СЛОЯ*

П. Е. Товстик1, Т. П. Товстик2

1. С.-Петербургский государственный университет, д-р физ.-мат. наук, профессор, Peter.tovstik@mail.ru

2. Институт Проблем машиноведения РАН,

канд. физ.-мат. наук, ст. научн. сотр., tovstik_t@mail.ru

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

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

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

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

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

В нашей работе проводятся исследования свойств плоского графитового слоя, построенного по модели [4]. Построена плотность упругой энергии как функция деформаций при аффинном преобразовании решетки. Рассматриваются деформации, напряже-

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

© П. Е. Товстик, Т. П. Товстик, 2009

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

1. Плоская решетка с парным моментным взаимодействием между частицами. Рассмотрим бесконечную двухмерную решетку частиц, каждая из которых имеет три степени свободы: декартовы координаты х, у и угол поворота у>. Парное мо-ментное взаимодействие между частицами описывается потенциалом (см. [4])

3к\

Щ{х2,У2, .г’ъ Уъ Ч>\) = и0 (»’, 7, к) = П0(г) + П^г) эт) — \ ,

По (г) = Б1

-2(е.у , П,И = С!

г г г

(1.1)

где (см. рис. 1а)

г = \](ж2 - Ж1)2 + (г/2 - ,У 1)2, ж2 - Х\ = г соэ в, г/2 - У\ = г вт в,

7 = в - (^1 + ^2)/2, к = ^2 - ^1-

(1.2)

©

<3-------------------В>

Рис. 1. Две частицы: а) углы, Ъ) устойчивое положение равновесия.

Первое слагаемое в (1.1) есть потенциал Леннарда—Джонса, а второе слагаемое описывает моментное взаимодействие между частицами. При П1 =0 моментное взаимодействие отсутствует и устойчивое равновесное положение двух изолированных частиц будет при г = р. Предполагая провести исследование в безразмерных переменных, считаем р = ^1 = 1. Примем значения В2 и т, наиболее близко соответствующие графитовому слою [4]:

Во =--------, т = 12. (1-3)

~ 1.26’ 1 ;

В этом случае существует множество устойчивых равновесных положений двух частиц при г = 0.76872, отличающихся друг от друга тем, что углы у>1 и ^>2 определены

с точностью до слагаемых, кратных 2п/3. При в = 0 (см. рис. 1Ь) будет

у>1 = п/3 + 2кп/3, <^>2 = 2кп/3, к = 0,1, 2. (1-4)

Остальные равновесные положения (у>1 = ^2 =0; у>1 = ^2 = п/3; у>1 = 0, ^2 = п/3) неустойчивы.

Частицы с потенциалом (1.1) могут образовывать устойчивую бесконечную решетку из правильных шестиугольников с расстоянием между частицами го = 0.75831 (в размерных величинах го* = 0.1395нм). Целью нашего анализа является рассмотрение аффинных деформаций этой решетки, определение возникающих при этом внутренних усилий и исследование устойчивости положений равновесия после деформации.

2. Аффинная деформация решетки. После аффинной деформации правильные шестиугольники превратятся в одинаковые шестиугольники с параллельными противоположными сторонами (см. рис. 2). Найдем возможные равновесные формы этих шестиугольников .

Рис. 2. Повторяющаяся ячейка: недеформированное (а) и деформированное (Ъ) положение.

Все частицы находятся в одинаковых условиях. Поэтому как в недеформированном, так и в деформированном положении будем вычислять потенциальную энергию V взаимодействия одной из частиц (например, частицы А с координатами х = у = 0, у>о) с остальными частицами:

V = ^ ио(х*,у*,^*; 0, 0, у>о), (2.1)

г

где суммирование распространяется на те частицы, взаимодействие которых с частицей А целесообразно учитывать (в силу (1.1) и (1.2) с ростом расстояния г между частицами силы взаимодействия быстро убывают).

После деформации все ячейки имеют форму шестиугольников с параллельными противоположными сторонами и согласованными с <^>о углами <^г (срг = ^о+кп/3). Пусть х, у —декартовы координаты частицы, имевшей до деформации координаты £, ц, и одна из сторон шестиугольника параллельна оси х. Форма ячейки описывается шестью параметрами

аг, г = 1, .. ., 5, о, (2.2)

(где размеры аг показаны на рис. 26). Имеют место формулы

дх 2а\ + а2 + а5 ду аз — а4

1 + £1 = ^ = Згп ’ ^1 = ^ = ^Г’

ду аз + а4 дх а2 — 0,5

связывающие размеры аг ячейки с компонентами градиента деформации Т:

Т = ( 1 + £1 1+2 ) . (2.4)

\ «1 1 + £2 у

(2.3)

Нелинейные деформации решетки єіз- считаем заданными:

1

1

£ц — £1 + ^(е1 + а,1)> е22 — £2 + 2 (е2 + ^2)1 £12 — ^1(1 + £2) + ^2(1 + £1)- (2-5)

В силу формул (2.3) деформации £г^ являются функциями размеров ячейки аг:

£^'к = £^'к (аг). (2.6)

Для вычисления координат хг,уг,фг частиц в потенциале (2.1) необходимо найти величины аг и угол фо. Величины аг в силу формул (2.3) и (2.5) удовлетворяют трем уравнениям (2.6). Не любая форма деформированной ячейки является равновесной. Три других ограничения, накладываемые на шесть величин (2.2), вытекают из трех условий равновесия частицы А. Для записи этих условий введем малые проекции перемещения и поворот (х1, у1, <£>1) частицы А. Тогда потенциал (2.1) примет вид

и = ^ Цз(жі, уі, фі, жі, уі, фо + фі). і=і

(2.7)

Условия равновесия частицы А имеют вид

= —

дії

джі

М = -

аг/

ду>і

(2.8)

где и М — проекции силы и момент, действующие на частицу А.

Используем следующую последовательность вычислений. Три из шести величин (2.2), а именно аі, а2, фо, будем определять из уравнений (2.8). Считая заданными величины аі, а2, єіз-, семь величин аз, а4, а5, єі, ^і находим из системы уравнений (2.3), (2.5). Эта система легко решается методом итераций. Перепишем ее в виде итерационной схемы:

(п+і)

(п+і) _

1 + 2єіі - ^и) - 1,

3го(1 + є і )) — 2аі — а2,

.("+і)

= (у0^5(1 + є£°) + 1-5и[п)) го,

1 + 2Є22 — — 1

(п)

є , ,(и)(і + є(п))

(п+1) _ £12 — Iі + £1 )

Лп)\

го

-,(п+і)

1 + є2п)

= (у0^5(1 + £2П)) - 1.5с^п)) го,

(о)

(о)

причем итерации следует начинать со значений єі ' = о>>'' =0.

После определения величин аі координаты вершин шестиугольника на рис. 2Ь получаем непосредственно, а координаты остальных точек деформированной решетки находим из условий периодичности. В результате потенциал (2.1) оказывается функцией шести величин и = и (є із, аі,а2,фо), причем величины єіз- считаем заданными, а величины аі,а2,фо находим из условий (2.8) равновесия точки А. После определения величин аі,а2,фо потенциал и в деформированном состоянии становится функцией трех параметров деформации єіз-, т. е. и = и(єіз-).

3. Деформационные свойства решетки. Будем рассматривать решетку как двухмерную сплошную среду, упругие свойства которой определяются потенциалом

0

0

і

2

5

2

і

з

4

и = и(е^-). В связи с тем, что упругая энергия и отнесена к одной частице, плотность энергии и*, отнесенная к единице площади до деформации, пропорциональна и:

и* = В и. (3.1)

В размерных переменных коэффициент пропорциональности В имеет вид

В=\В^ 8=^г1, (3.2)

где Б — площадь, приходящаяся на одну частицу до деформации, и р — те же, что

и в формуле (1.1). Множитель 1/2 связан с тем, что в потенциале (1.1) участвуют две

частицы. Ниже нас будут интересовать относительные величины, поэтому множитель В опускаем.

Производные

Кц = д.7Гз (3'3)

дают компоненты энергетического тензора напряжений Кирхгофа [7], а вторые производные

Сцы = т,—т,— (3.4)

06^ 06^1

определяют коэффициенты жесткости среды.

Для вычисления производных функции и (ж) используется численное дифференцирование типа

<ш

йж

сРи

(1х2

_и{К)-и{-К) , и{-2К)-2и{-К) + 2и{К)-и{2К) х=0 21 + 12Н :

и(й) + и(-й) - 2и(0) и(-2й) - 4и(-й) + би(0) - 4и(й) + и(2й)

0 й2 12й2

(3.5)

При отсутствии начальных деформаций (е^- = 0) среда изотропна. Это подтверждается тем, что модуль Юнга Е = 475 и коэффициент Пуассона V = 0.241 удалось подобрать таким образом, чтобы одновременно были выполнены соотношения

Е Ev Е

Сцц = С*2222 = ~л------2 = ^4, С1122 = ~л---й = ^2, С1212 = Г = О = 191.

1 - V2 1 - V2 2(1 + V)

(3.6)

Если в сумме (2.1) удерживать достаточно большое число слагаемых, соотношения (3.6) выполняются с высокой степенью точности (а не только с тремя десятичными знаками, которые удержаны в (3.6)). Таблица 1 иллюстрирует нарушение изотропии, если в сумме (2.1) брать N =12 или N = 3 частиц, ближайших к частице А.

Таблица 1. Зависимость коэффициентов жесткости от числа учитываемых частиц

N Сцц = С2222 СЦ22 О ю ю Е V С

оо 504 122 191 475 0.241 191

12 510 124 192 480 0.243 193

3 538 134 194 504 0.250 202

Таблица 2. Зависимость коэффициентов жесткости от начальной деформации

ЕЦ £22 £12 Сцц С2222 О ю ю СЦ22 Сц 12 О ю ю ю

0 0 0 504 504 191 122 0 0

0.1 0 0 27 415 59 26 0 0

0 0.1 0 301 42 93 46 0 0

0 0 0.1 643 535 207 193 -237 -117

В первой строке таблицы 1 повторена информация, содержащаяся в формулах (3.6). В последнем столбце таблицы 1 приведен модуль сдвига О, вычисленный по формуле О = Е/(2(1 + V)). Для изотропного материала должно быть О = С1212, что при N =12 и при N = 3 не имеет места.

При наличии начальной деформации среда становится анизотропной. Степень анизотропии иллюстрируется в таблице 2 и на рис. 3.

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

При построении таблицы 2 и далее расчеты проводились при достаточно большом числе N удерживаемых в сумме (2.7) слагаемых. В связи с использованием численного дифференцирования значения потенциала V в (3.5) нужно знать с высокой точностью для получения нужного числа верных знаков в результате. Например, при построении области устойчивости в п. 5 для получения трех верных десятичных знаков было наложено ограничение \/(х.,, — .г’1)2 + (г/* — г/1)2 < 40, при котором сумма (2.7) содержит N = 6729 слагаемых.

^ \ уС2222^—) \ Л 1500 ■ \ \ \ \ \ \ ^■1111 (Бц) с2222^22)

Сии<еи> Чч\ 1000 \ч NN.

5о9~ N..

0 ж 611 У622

1 1 1 —1 =г" а

-0.1 -0.05 0.05 0.1

Рис. 3. Жесткости при больших деформациях.

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

С1111 (ец) = 0 при ец =0.110; С2222 (е22)=0 при е22 = 0.132, (3.7)

а при дальнейшем росте деформаций становятся отрицательными, что связано с особенностями потенциала По (г) в (1.1).

4. Устойчивость решетки после аффинной деформации. Рассматривается

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

шетке частиц с парным потенциалом взаимодействия (1.1). Упругие свойства определяются потенциалом и(ег,). Вопрос о выходе частиц из плоскости не рассматривается. Уравнения равновесия в проекциях на оси до деформации имеют вид [7]

др - -

^ • Р + Г = О ИЛИ "^■+Л'=0’ І = 1;2, (4.1)

где рг, — компоненты тензора напряжений Пиолы, связанные с напряжениями Кирхгофа (3.3) формулами

_ _т т, дх, ди / ди, N

к‘ = к''дГг=Щ,{,р + ^)' ( ’

причем (^1,^2) = (£, п), (х1,х2) = (х, у) = (£1 + М1, £2 + м2), а иьи2 —проекции пере-

мещения частицы, 5,р — символ Кронекера, Т — градиент деформации (2.4), а верхний индекс Т означает транспонирование. Через fj в (4.1) обозначены проекции на оси интенсивности внешней нагрузки до деформации. Здесь и далее предполагается суммирование по немым индексам.

Для исследования устойчивости среды рассмотрим уравнения в вариациях для системы (4.1), положив иг = и0 + V*. В тензорной записи эти уравнения имеют вид [8, 9]

где Е — тензор деформации Коши—Грина, индекс 0 у производных энергии по тензору Е и у тензора Т отмечает, что величины вычисляются при V* = 0. Вариации Е и Т тензоров Е и Т вычисляются по формулам

Ё = 1- (ТГТ0 + Т0тТ) , Т={|^}. (4.4)

В качестве вариации внешней нагрузки возьмем силы инерции

А - (4.5)

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

В развернутой записи система уравнений (4.3) имеет вид

где — акустический тензор (четвертого ранга) [9].

Имея в виду использовать динамический критерий устойчивости, рассмотрим решение системы (4.6) в виде распространяющейся волны

V]. = г^ехр{*(рі£і +Р2& - ^)}, к= 1,2, і = \Гл. (4.7)

Если при некоторых волновых числах рі, р2 величина Л = ро<^2, определяемая из уравнения

Й(Л)

іі Л і2

а — Л а

2і

а21 а

2і 22 Л

а

ы = Акірі + 2А2рір2 + Ак2 р2, (4.8)

0

Рис. 4- Область устойчивости.

оказывается неположительной, то среда неустойчива, ибо распространение волны сопровождается ростом амплитуды. Критерием устойчивости в силу а12 = а21 является положительная определенность трех функций

gi(pi,p2)= а11, g2(pi,p2)= а22, g3(pi,p2)= а11а22 - (а12)2 . (4.9)

Приведем выражения для компонентов акустического тензора, необходимых в силу (4.8) для проверки выполнения условий устойчивости (4.9):

А11 = ^ 1C1111 + 2£1^2С1112 + ^2C1212 + K1b

АЦ = (t 1^2 + W1W2)C! 112 + t1^1Cnn + t2^2^1212,

A11 = t2Cn12 + 11^2(^1122 + C1212) + ^2^2212 + K12,

A12 = (t1t2 + <^1^2)(C1122 + C1212)/2 + t1^1C1 112 + £2^2^2212,

A22 = t2C1212 + 2t1^2C2212 + ^2^2222 + K22, (4.10)

A22 = (t1t2 + <^1^2)622 12 + £1^1^1212 + t2^2^2222,

A11 = £ 2 C1212 + 2t2^1C1112 + ^2C1111 + K1b

A22 = t2C2212 + £2^1(^1122 + C1212) + ^2C1112 + K12,

A22 = £2^2222 + 2t2^1C2212 + ^2C1212 + K22,

Здесь для краткости введены обозначения £1 = 1 + £1, £2 = 1 + £2.

5. Область устойчивости. В пространстве деформаций £ц, £22, £12 область устойчивости симметрична относительно плоскости £12 = 0. Приведем сначала размеры этой области:

-0.143 < £ц < 0.118 при £22 = £12 = 0,

-0.140 < £22 < 0.152 при £ц = £12 = 0, (5.1)

-0.154 < £12 < 0.154 при £ц = £22 = 0,

В частности, сравнение со значениями (3.7) показывает, что положительная определенность матрицы d2U/dE2 = {Cjki} не может служить критерием устойчивости, ибо нарушение ее положительной определенности может предшествовать потере устойчивости.

На рис. 4 представлены сечения области устойчивости четырьмя плоскостями £12 = 0, £12 = 0.05, £12 = 0.1 и £12 = 0.15. Последнее значение £12 близко к максимально возможному значению £12 = 0.154.

Границу области устойчивости определяют линии, на которых нарушается положительная определенность хотя бы одной из функций (4.9). Угловые точки линий на рис. 4 связаны с тем, что в этих точках положительная определенность нарушается сразу для двух функций (4.9).

При £12 = 0 граничной точкой области устойчивости является, по-видимому, точка £11 = £22 = -0.5, соответствующая неограниченному сближению частиц. Однако в связи с потерей точности приблизиться к этой точке не удалось.

Авторы благодарят Е. А. Иванову, А. М. Кривцова и Н. Ф. Морозова за весьма полезное обсуждение представленных здесь результатов.

Литература

1. Иванова Е.А., Кривцов А. М., Морозов Н. Ф., Фирсова А. Д. Теоретическая механика. Описание механических свойств кристаллических твердых тел на микро- и макроуровне. СПб.: Изд-во СПбГТУ. 2003. 32 с.

2. Иванова Е. А., Кривцов А. М., Морозов Н. Ф. Получение макроскопических соотношений упругости сложных кристаллических решеток с учетом моментных взаимодействий на микроуровне // Прикл. матем. и механ. 2007. Т. 71. Вып. 4. С. 595-615.

3. Бызов А. П., Иванова Е. А. Математическое моделирование моментных взаимодействий частиц с вращательными степенями свободы // Научно-технические ведомости СПбГПУ. 2007. №2. С. 260-268.

4. Беринский И. Е., Иванова Е. А., Кривцов А. М., Морозов Н. Ф. Применение моментного взаимодействия к построению устойчивой модели кристаллической решетки графита // Изв. РАН, МТТ. 2007. №5. С. 6-16.

5. Кривцов А. М. Деформирование и разрушение твердых тел с микроструктурой. М.: ФИЗМАТЛИТ, 2007. 304 с.

6. Ivanova E.A., Krivtsov A.M., Zhilin P. A. Description of rotational molecular spectra by means of an approach based on rational mechanics // ZAMM. Z Angew. Math. Mech. 2007. Vol. 87, №2. P. 139-149.

7. Кошелев А. И., Нарбут М. А. Лекции по механике деформируемого твердого тела. СПб.: Изд-во С.-Петерб. ун-та, 2003. 276 с.

8. Елисеев В. В. Механика упругих тел. СПб.: Изд-во СПбГТУ, 1999. 340 с.

9. Лурье А. И. Нелинейная теория упругости. М.: Наука, 1980. 512 с.

Статья поступила в редакцию 19 марта 2009 г.