Определение изгибной жесткости графенового листа Текст научной статьи по специальности «Физика»

УДК 539.2 + 539.3

Определение изгибной жесткости графенового листа

И.Е. Беринский1,2, A.M. Кривцов1,2, A.M. Кударова3

1 Институт проблем машиноведения РАН, Санкт-Петербург, 199178, Россия 2 Санкт-Петербургский государственный политехнический университет, Санкт-Петербург, 195251, Россия 3 Делфтский технический университет, Делфт, 2628 CN, Нидерланды

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

Ключевые слова: изгибная жесткость, графен, моментное взаимодействие, дискретные модели, стержневые модели

Bending stiffness of a graphene sheet

I.E. Berinskii1-2, A.M. Krivtsov12, and A.M. Kudarova3

1 Institute of Problems of Mechanical Engineering, RAS, St. Petersburg, 199178, Russia 2 St. Petersburg State Polytechnical University, St. Petersburg, 195251, Russia 3 Delft University of Technology, Delft, 2628 CN, Netherlands

The paper proposes a discrete mechanical model of graphene sheets. A relationship between parameters of the model and elastic characteristics of its equivalent continuum is found by comparing the energy of small strains on micro- and macroscales. This makes it possible to determine the microscale interaction parameters from known experimental data and, knowing the microscale parameters, to determine the mechanical properties of graphene. The main aim of the work is to estimate the bending stiffness of a graphene sheet. With the proposed discrete model, an analytical dependence of the graphene sheet bending stiffness on the microscale interaction parameters is derived.

Keywords: bending stiffness, graphene, moment interaction, discrete models, rod models

1. Введение

В последнее время активно развивается область исследований и разработок наноэлектромеханических систем. В частности, ведется поиск методов создания нано-резонаторов нового поколения. Существующие кварцевые резонаторы обладают высокой частотой и добротностью (до 400 МГц и до 2500 соответственно), однако при уменьшении их толщины до нанометровых размеров поверхностные эффекты ухудшают их свойства. Кроме того, существуют трудности в изготовлении тонких кварцевых пластин с высокой степенью параллельности рабочих сторон, что сильно влияет на появление резонанса на гармониках вблизи рабочей частоты [1]. Поэтому внимание исследователей привлекают попытки создания резонаторов на базе углеродных нанострук-

тур, в частности графена. Графен — монослой графита — является самым тонким известным материалом, т.к. содержит лишь один слой атомов. Поскольку работы с графеном начаты недавно [2], графеновые резонаторы находятся на начальной стадии разработки [3, 4] и существует ряд актуальных проблем, связанных с их проектированием и производством.

Механические свойства графеновых резонаторов определяются свойством листа графена сопротивляться как растяжению, так и изгибу. Наличие изгибной жесткости (способности упругой оболочки сопротивляться выходу из плоскости) — это принципиальное свойство графена. Р. Пайерлс [5] и Л.Д. Ландау [6] показали, что бесконечный идеальный двумерный кристалл теряет устойчивость под действием тепловых флуктуаций. Од-

© Беринский И.Е., Кривцов A.M., Кударова A.M., 2014

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

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

М = Dк. (1)

Жесткость D выражается через модуль Юнга пластины Е, коэффициент Пуассона V и толщину h [7]:

Однако для однослойного графена невозможно однозначно определить его толщину, что не позволяет использовать напрямую формулу (2). Используем для графена параметры Eh = 340 Н/м и v = 0.17, которые могут быть получены из экспериментальных данных, приведенных в [8]. Также, следуя наиболее распространенной идее [9-14], примем толщину графена равной расстоянию между слоями графена в графите h = 0.34 нм. Тогда получим D = 3.37 нН-нм. Однако это значение на порядок превышает оценки изгибной жесткости однослойного графена (табл. 1), полученные другими методами, не использующими понятие толщины, а именно ab initio (квантово-механические расчеты из первых принципов) и с помощью эмпирических потенциалов взаимодействия.

Как видно из табл. 1, теоретические оценки изгибной жесткости имеют один порядок, хотя различаются до двух раз, меняясь от 0.13 до 0.26 нН-нм (0.8-1.6 эВ). Попытки подтвердить эти данные экспериментально также сталкиваются с проблемой невозможности использования классических теорий сплошной среды. Для многослойного графена такие теории используются, и значения изгибной жесткости определены экспериментально. Так, в [18] методом «прощелкивания» двухслойной графеновой мембраны (потери устойчивости вследствие действия сосредоточенной силы) получено значение изгибной жесткости 5.7 нН - нм. В работе [19] изгибная жесткость определялась методом индентиро-вания графеновой мембраны в круглое отверстие. Было показано, что экспериментальные данные согласуются

Таблица 1

Изгибная жесткость графена, определенная различными способами

Источник D,нН - нм Способ

[15] 0.24

[12] 0.26 Ab initio

[13] 0.23

[14] 0.13 Эмпирический потенциал

[16] 0.13

[17] 0.22

с кривой, описываемой зависимостью (2). В работе [20] приводится следующая зависимость изгибной жесткости от числа слоев N в многослойной мембране: D = = 9.8N 2 нН • нм. Это значение получено из оценки «разглаживания» податливой ребристой подложки при закреплении на ней графеновой мембраны. Таким образом, в двух последних упомянутых работах экспериментально подтверждается зависимость изгибной жесткости от куба толщины для большого числа слоев. Однако в случае однослойной мембраны или имеющей несколько слоев зависимость (2) неприменима.

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

D = Л1(Ы - 1)Ы (Ы + 1) + Л2(3Ы - 1). (3)

Здесь Л1 и Л2 — постоянные коэффициенты, зависящие от свойств кристаллической решетки. Таким образом, изгибная жесткость нанокристалла складывается из двух слагаемых, первое из которых обращается в ноль для одного слоя и стремится к (2) при увеличении числа слоев. Наличие второго слагаемого вызвано учетом моментного взаимодействия в дополнение к силовому. Эта поправка становится малой по сравнению с первым слагаемым с увеличением толщины нанокрис-талла, но существенна для нанокристаллов, содержащих малое число слоев и однослойных структур.

В данной работе предлагается подход, позволяющий определить изгибную жесткость одного слоя графена. Как и в [21], этот подход основан на использовании моментного взаимодействия в дополнение к силовому. Считается, что атомы углерода в графене взаимодействуют между собой посредством сил и моментов. В разделе 2.1 приведена сводка уравнений, описывающих это взаимодействие. Переход от микроструктуры материала к континуальному уровню рассмотрен в разделе 2.2. Этот переход позволяет связать жесткости силового и моментного взаимодействия с упругими характеристиками материала. Таким образом, изгибная жесткость графенового листа оказывается функцией продольной, поперечной, изгибной и крутильной жесткости углеродной связи. В разделе 2.3 рассмотрена модель угле-

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

2. Определение изгибной жесткости графенового листа

2.1. Взаимодействие на микроуровне

Для описания взаимодействий в кристаллической решетке графена используется подход, предложенный в работах [21, 22]. Атом углерода в графене моделируется телом-точкой, т.е. материальным объектом, занимающим нулевой объем в пространстве, положение которого считается определенным, если задан его вектор положения и тензор поворота. Взаимодействие между телами-точками характеризуется вектором силы и вектором момента.

Рассмотрим систему из двух тел-точек, моделирующих атомы кристаллической решетки. В актуальной конфигурации их положение задается радиус-векторами Г[, г2, а ориентация — векторами поворотов ф1; ф2. В равновесном положении г2 - г = г0, ф1 = 0, ф2 = 0. Вводятся обозначения: т1 — сила и момент, действующие на тело-точку 1 со стороны тела-точки 2; т2 — сила и момент, действующие на тело-точку 2 со стороны тела-точки 1. Для них справедливы соотношения:

f = ^ =^2, (4)

т = т1 +1/2 (г - Г2)х £1 = -т2 -1/2 (Г2 - г )х £2 . (5) Рассматривая случай линейного упругого деформирования, для внутренней энергии можно принять следующую аппроксимацию:

и = 1"0 • £ + т о • к +1/2 £ • А • £ + £ • В • к + + 1/2 к • С • к. (6)

Коэффициенты А, В и С называются тензорами жесткости связей, а векторы ^ и т0 — это начальные усилия. В линейной теории тензоры жесткости — постоянные величины, причем тензоры А и С симметричные, а В произвольный. На векторах е и к, называемых векторами деформации, совершают работу векторы силы и момента:

£ = £0 + А • £ + В к, т = т0 + £ • В + С к. (7)

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

£ = г - Г +1/2 г0 х (Ф1 + Ф2),

(8)

к = Ф2 - Ф1> Г = Г2 - Г1-

Уравнения в виде (7), (8) удобны для определения структуры тензоров жесткости. Характеристики взаимодействия рассматриваемой системы из двух частиц были вычислены относительно центра этой системы, поэтому она имеет две ортогональные плоскости симметрии. Можно показать, что тензоры жесткости в этом случае имеют вид

А = САп + Св jj, В = 0,

С = СТ и + Св (jj + кк), (9)

где j и к — векторы ортонормированного базиса, при этом частицы лежат в плоскости, образованной векторами i и j; к ортогонален этой плоскости.

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

£1 = Г -Г0 + Г0 ХФ2> к1 = Ф2 -Ф^ Г = Г2 -Г1- (10) Векторы силы и момента в этом случае могут быть представлены в виде

£1 = %0 + Та • £1 + тв • к 1,

т1 = т 0 + £1 • Тв + Тс • к1.

Начальные усилия %0 и т 0 отличаются от начальных усилий ^ и т0, а тензоры жесткости ТА, Тв и ТС отличаются от тензоров А, В и С. В дальнейшем будем рассматривать ненапряженные системы, т.е. начальные усилия будем считать равными нулю. Нетрудно убедиться, что верны следующие соотношения:

та = А, Тв = В -1/2 А хГ0,

т (12)

Тс = С +1/2(г0 хВ - Вт хГ0) -1/4г0 х АхГ0.

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

г = ^£Па, (13)

2'0 а

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

(15)

Потенциальная энергия Па может быть представлена как квадратичная форма векторов деформации и тензоров жесткости связи:

Па = V2 га • Аа • га + га • Ва •ка +

+ 1/2 К а • Са- К а, (14)

где А а, В а, Са — тензоры жесткости связи а, которые содержат информацию о деформировании связи в разных направлениях:

Аа = СЛПаПа + ^ (Е - ПаПа X

Са = СТ П ап а + Св (Е - П аП а ), аП а = а а. Здесь Е — единичный тензор; вектор п а = а а/а, где а — длина связи; аа соединяет две соседние частицы. Тензор В равен нулю, если решетка имеет две взаимно перпендикулярные плоскости симметрии [16], что верно в случае решетки графена.

Векторы деформаций можно представить в виде [16]:

га = иа-и + V2ааХ(Фа + фХ ка = Фа-ф. (16) Здесь и а, и — это смещения частицы с индексом а и некоторой отсчетной частицы; фа, ф — их повороты.

Уравнения (14), (15) можно подставить в (13), и таким образом будет получена потенциальная энергия системы частиц в зависимости от их положения и а, и, ориентации ф а, ф, геометрии начальной конфигурации аа и силовых констант взаимодействия СЛ, Св, Ст, СГ1.

2.2. Переход на макроуровень

Сопоставим смещениям и поворотам частицы смещения и повороты элемента сплошной среды: и = и(г), ф = ф (г), где г — радиус-вектор в начальной конфигурации. Тогда смещения и повороты соседних частиц можно представить в виде иа = и(г - аа), фа = ф(г - аа). Воспользуемся длинноволновым приближением [23], которое состоит в том, что рассматриваемые длины волн значительно превосходят начальное межатомное расстояние между ближайшими соседями аа, принимаемое малым параметром. Тогда смещения и повороты можно представить в виде разложения

иа = и + аа ^и + фа= ф + аа^ф + ПУ- (17)

Здесь V — дифференциальный набла-оператор. Параметр п равен единице, если в элементарной ячейке решетки содержится две частицы; если одна частица, он равен нулю. Векторы у обозначают относительные смещения и повороты частиц из разных подре-шеток, составляющих полную решетку. Выражения (17) могут быть подставлены в энергию W, вычисленную в предыдущем параграфе. Затем необходимо найти векторы у из условия обеспечения ими сдвига одной подрешетки относительно другой, соответствующего минимуму энергии деформирования:

После получения выражения для векторов у из условия (18) плотность энергии эквивалентного континуума W будет иметь вид

W = W(Vu, Vф, ф, Сл, CD, Ст, Св). (19)

Также она может быть представлена в виде квадратичной формы тензоров жесткости четвертого ранга и тензоров деформаций:

W = 1/22 в "4 А-2 г + 2 г -4 В-2 к + +1/22 к -4 С ••2 к. (20)

Тензоры деформаций имеют вид

2 г ^и + Е Хф, 2 к ^ф. (21)

Для получения соотношений между силовыми константами взаимодействий на микроуровне и компонентами тензоров жесткости необходимо сравнить выражения для энергий (19) и (20).

Теперь применим описанный выше метод к кристаллической решетке графена. Решетка графена — двухмерная решетка, в которой атомы расположены в узлах шестиугольников. Введем ортонормированный базис (е1; е2, е3}, где векторы е1 и е2 лежат в плоскости решетки, а вектор е3 перпендикулярен плоскости решетки. Обозначим оси х, у и г, сонаправленные с векторами е1; е2, е3 соответственно. Тогда векторы направления связей можно представить как

п1 = е1, п2 =-1/2е1 + Т3/2е2,

п3 =-1/2 е1 -л/3/2е2. ( )

Смещения и повороты также можно разложить по базису:

и = их (х, у)е1 + иу (х, у)е2 + и2 (х, у)е3,

(23)

ф = ф (х, у)е1 (х, у )е2 (х, у)ез.

Тогда нетрудно произвести сравнение выражений для энергий (19) и (20), где энергия (19) получена для решетки графена с учетом (22). Для сравнения нужно приравнять множители перед производными компонент вектора смещений и вектора поворота, а также перед компонентами вектора поворота в выражениях (19) и (20), предварительно представив их в покомпонентной записи. В итоге получим соотношения между коэффициентами тензоров жесткости и параметров взаимодействия на микроуровне:

,л/з Сл (Сл + 3CD)

dW = 0, dW = 0. Эд Эу

(18)

A1111 = A2222

A1122 = A2211

A1212 = A2121

A1221 = A2112

6 ca + Cd V3 Ca (Ca - cd) 6 Ca + Cd '

V3 Cd (Cd + 3CA )

6 Ca + Cd '

V3 Cd (Ca - Cd ) 6 Ca + Cd '

A3131 = A3232 __^CD'

C1111 = C2222

C1122 = C2211

C1212 = C2121

C1221 = C2112

£

3

Уз CT (CT + 3CB) 6 CT + CB '

V? Ct (Ct - Cb ) 6 CT + CB '

л/3 CB (CB + 3CT)

"6 CT + CB ''

V3 CB (cT - cB )

6 CT + CB '

(24)

С3131 = С3232 СВ-

В уравнениях (24) представлены ненулевые компоненты тензора 4 А и тензора 4 С. Эквивалентный решетке графена континуум инвариантен относительно поворотов вокруг нормали е3 и отражений от базисных плоскостей. Тензоры жесткости должны иметь такую же симметрию согласно принципу Кюри [18], что отражено в (24): только такой набор ненулевых компонент и равенство нулю тензора 4 В обеспечивают указанную симметрию. Компоненты из (24) не являются независимыми, поэтому можно выразить некоторые из них через другие:

A1212 + A1221 = A1111 A1122> C1212 + C1221 = C1111 - C1122-

(25)

Еще одно соотношение может быть получено из системы уравнений (24):

A3131

A2 - A2 A1111 A1122

3A1122 + A1111

C1221 =

C1122(C1111 C1122 ) 3C1122 + C1111

(26)

Тензор А связывает тензор силовых напряжений Т и тензор деформаций е, тензор 4 С — тензор моментных напряжений М и тензор деформаций к:

Т11 = А1111е11 + А1122е22, Т22 = А2211е11 + А2222е22>

Т21 = А1221е12 + А1212е21, Т12 = А2121е12 + А2112е21>

(27)

T13 = A3131e13> T23 = A3232e23> M11 = C1111K11 + C1122 K22> M 22 = C2211K11 + C2222 K22> M 21 = C1221K12 + C1212K21> M12 = C2121K12 + C2112 K21, M13 = C3131K13> M 23 = C3232 K23-

Компонентами тензоров деформации e, k являются e11 = Ux> e12 = u,X -Ф2> e13 = u,X +ФУ> e21 = Uy + Ф2

е22 = иУу, £23 = и2у-Ф", Кц = ф", к21 = ф*, (28) к12 =Ф,х, К22 =фУу, К13 =Ф,х, К23 =Ф,2у

Исходя из этих соотношений можно определить физический смысл различных компонент тензоров жесткости. Коэффициент А1111 описывает свойства материала на растяжение, А1122 — эффект Пуассона, А1212 и А1221 связывают сдвиговые деформации и напряжения в плоскости векторов е1; е2, А3131 — в ортогональной плоскости, С1111 описывает свойства материала на кручение, С1221 — аналог эффекта Пуассона в случае кручений.

Определим изгибную жесткость материала D как коэффициент между моментным напряжением в сечении, ортогональном е1; и деформацией к12 (или е2 и к21). Тогда этот коэффициент равен С1212 (С2121):

л/3 Св (Св + ЗСт )

D = -

Cb + Ct

(29)

К сожалению, для определения всех упругих модулей не хватает экспериментальных данных. Для уменьшения количества независимых модулей рассмотрим упрощенную теорию. Предположим, что в плоскости отсутствуют моментные напряжения. Тогда тензор жесткости 4 А должен быть инвариантен по отношению к транспозиции пар индексов 12 и 21. По этой причине введем новую компоненту А1*212, соответствующую упрощенной теории, в виде (А1212 + А1221)/2:

л/3 СаСв

A1212 = A1291 = A?1?1 = A

A1221

2121

A2112

3 CA + CD

(30)

Результаты, полученные для коэффициентов A1111 и A1122, описывающих деформацию в плоскости, совпадают с результатами, полученными в [22], где экспериментальные данные для графита применялись для вычисления силовых констант CA и CD - Выяснилось, что Cd/Ca = 0.55, это говорит о том, что поперечная жесткость сравнима с продольной и ее необходимо привлекать во внимание.

Из выражений (24) видно, что эффект Пуассона исчезает, если CA = CD - Подобный эффект происходит и при CT = CB: нет скручивания материала в направлении, ортогональном направлению приложения крутящего момента.

2.3. Стержневая модель углеродной связи в графене

В настоящее время предложено множество подходов, которые могут успешно использоваться для моделирования динамики графена на микроуровне. Среди них, помимо ab initio и молекулярно-динамических подходов [10-17, 24], можно отметить молекулярно-меха-нические [25, 26].Отличие подхода, предлагаемого в данной работе, заключается в том, что предлагается конкретная механическая модель углеродной связи —

упругий стержень, работающий на сжатие-растяжение, изгиб и кручение. Подобным стержнем моделируется взаимодействие электронных облаков, которые формируют направленную ковалентную химическую связь. Ясно, что упругие характеристики этого стержня должны быть подобраны так, чтобы полученная в результате модель удовлетворяла экспериментальным значениям упругих характеристик исследуемого кристалла. С другой стороны, в работах [22, 27] было показано, что эта же связь может быть описана с помощью моментного подхода, причем характеристики моментного взаимодействия можно определить для графена и алмаза.

Для удобства читателя далее приведена сводка основных уравнений линейной теории прямолинейных стержней, составленная на основе [26], а также приведены несколько задач, решения которых понадобятся в дальнейшем. Представим, что смещения точек стержня описываются следующими уравнениями: и = ^ + w, t • = 0,

у = у + t х 0, 0 • t = 0. (31)

Здесь и — это продольные смещения точек стержня; w — вектор поперечных смещений; у—кручение; t — единичный вектор касательной. Векторы деформации задаются соотношениями

е = Et + у, е = и, у = -те' - 0, ф = уЧ + t х 0'. (32) Здесь е — относительное удлинение стержня; у — вектор деформации поперечного сдвига; У — относительное закручивание стержня; 0 — вектор изгибных деформаций.

В линейной теории стержней без учета естественной крутки соотношения упругости для действующих в стержне усилий и моментов принимают следующий вид: п = А • е, т = С • ф,

п = П + q, t • q = 0, (33)

т = Н +1 х1, t • 1 = 0. Здесь Т— продольное усилие в стержне; Н — крутящий момент; q — вектор поперечных усилий; 1 — вектор изгибающих моментов.

Тензоры напряжений А и С строятся следующим образом:

A = EFtt + kGF(E - tt), С = GJrtt + E 2 c,

(34)

2с = + J2d2d2, k = П2/12.

В приведенных выше формулах Е — это модуль Юнга; F — площадь поперечного сечения стержня; G — модуль сдвига; к—коэффициент поперечного сдвига; d1 и d2 — векторы нормали и бинормали; /1 и /2 — соответствующие моменты инерции поперечного сечения стержня. За /г обозначена геометрическая жесткость на кручение, причем для эллиптического в сечении стержня она равна полярному моменту инерции /р.

С учетом используемых обозначений без учета внешних распределенных сил и моментов уравнения равно-

весия стержня примут следующий вид:

п = 0, т +1Х п(я) = 0, (35)

где 5 — естественная координата стержня.

Рассмотрим систему из двух частиц, одна из которых жестко закреплена, а вторая сместилась относительно первой на трансляционный вектор и* и вектор поворота ф*. Согласно моментному подходу, между частицами начали действовать сила и момент взаимодействия (7), которые характеризуются компонентами Сл, Св, Ст, CD. С другой стороны, можно представить, что рассматриваемые частицы соединены упругим стержнем, левый конец которого закреплен, а правый смещен из положения равновесия на вектор и* и повернут на угол ф*. В результате в стержне возникают усилия, зависящие от и* и ф*. Уравнения статики (35) для них дают

n = n0 = const,

m = m0 -1 xn0

(36)

На основании (36) введем обозначения для силы и момента, вычисленных на конце стержня:

n* = nо = const, m* = m0 -1 xn0l. (37)

С другой стороны, согласно (33) n = A - (u'+1 Хф), m = С - ф'. (38)

Решая совместно системы (36) и (38), получаем выражения для векторов смещения и поворота:

u = A-1 - n 0s -1 x

x(C-1 - (1/2 m0s2 -161 x n0s3) + ф0s) + u0, (39)

ф = C-1 - (m 0s -1/21 x n0s2) + ф0. С учетом граничных условий

"Is=0 = 0 ф|s=0 = 0 "Is=l = "*' ф|s=l = Ф* (40)

находим

"0 = 0, Ф0 = 0,

ф* = C-1 - (m 0l -121 x n0l2), (41)

u* = A-1 -n0l -1 x (C-1 - (12m0l2 -161 xn0l3)). Разрешим системы (37) и (41) относительно u* и ф*. В результате действующую на конце стержня силу можно привести к виду

n* = A -s*, s* = u* +121 xФ*, (42)

A = (A-1 -1/12 (t x С-1 x t )l3 )-1. Сравнивая последние соотношения с (7)-( 10), видим, что нам удалось получить силу взаимодействия в той же форме, в какой она была получена при использовании дискретного подхода. Для определения компонент тензора A используем следующее свойство: пусть некий тензор второго ранга представляется в виде

Л = А111 + A2d1d1 +A1d2d2. (43)

Нетрудно убедиться, что обратный тензор в этом случае имеет вид

Л-1 = 111 +1А2 d1d1 +1/^ d2d2. (44)

Используя это и подставляя выражения (34), найдем интересующие нас компоненты:

Сл = t • А • t = EF/l, CD = ^ • А • ^ =■

12kE/2 F

(45)

(48)

Ш3 + 24/2(1+у ) Уравнения теории стержней, которые были использованы нами для получения (45), учитывают деформацию поперечного сдвига, что соответствует модели балки Тимошенко. В результате использования этой модели коэффициент Пуассона V вошел в выражение для поперечной жесткости CD как независимый параметр. Однако если необходимости учета деформации поперечного сдвига нет, то можно воспользоваться более простой моделью Бернулли-Эйлера. Переход к этой модели может быть произведен, если положить в (45) к ^ Тогда получим

С EF С НЕ/2

СЛ = I > CD = "з-• (46)

Момент, вычисленный на конце стержня, можно представить в виде, аналогичном (11):

т* = г* • В + С • к*, (47)

где

В = 121 х А, (С = 1/1С -141 Х АХЫ2,

г* = и + ^ х ф*)1, к* = ф*. Используя (11), нетрудно убедиться, что изгибная жесткость Св может быть найдена по формуле

Св = 11d2 • С • d2 = (Е/2уI. (49)

В качестве замечания отметим, что из выражений (48) легко также получить жесткость на кручение упругого стержня вокруг своей оси, аналога которой нет в рассматриваемой нами моментной постановке:

Ст = 1/11 • С • t = (в/г )/1. (50)

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

2.4. Вычисление изгибной жесткости графенового листа

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

F =

/2 =

~64~

/ =

т1

(51)

46464

Здесь d — диаметр стержня. Воспользуемся моделью Эйлера для определения упругих характеристик стержня. Из (46) с учетом (51) получим выражения для модуля Юнга и диаметра через значения продольной и поперечной жесткостей связи:

Е = 3Сл

п1СГ

(52)

которые полностью определяют упругие свойства стержня для данной модели.

Подставляя (51), (52) в (49) и (50), можем получить выражения для изгибной и крутильной жесткости стержня:

= Cлd = Сг1 ="

С =. Сл^

СГ2

(53)

16 12 ' ' 16(1 + у) 12(1 + у) Подставив выражения (53) в формулу (29), получим выражение для изгибной жесткости графенового листа

Г = — 2 = —к(у)СГа2,

48 л 36

(54)

кМ =

4 + V

4 + 2v

Изгибная жесткость содержит независимый параметр — коэффициент Пуассона, который, вообще говоря, остается неопределенным. Для этого параметра выполняется ограничение -1.0 < V < 0.5, что дает для коэффициента к^) оценку:

9/10 < к(v) < 3/2. (55)

Параметры Сл и СГ могут быть однозначно определены из упругих характеристик материала. В частности, в [22] показано, как определить такие характеристики для графена на основе экспериментальных данных из [3]:

Сл = 730.2 Н/м, СГ = 401.6 Н/м, 1 = 0.142 нм. (56) Тогда, подставляя (55), (56) в (54), получаем

0.35 < D < 0.58 нН • нм. (57)

3. Обсуждение результатов и выводы

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

Подставляя (56) в (52), получим Е = 8.928 ТПа, d = = 0.122 нм, т.е. диаметр стержня оказывается сопоставим с его длиной, и модель Бернулли-Эйлера может давать ощутимую погрешность при оценке изгиба. Правильнее было бы использовать модель Тимошенко, однако в этом случае появляется еще один неопределенный параметр — коэффициент поперечного сдвига к. Вообще говоря, этот коэффициент зависит от формы сечения. Для стержней и оболочек переход от теории, учитывающей сдвиг, к классической теории происходит посредством предельного перехода к ^ <». В теории оболочек при учете сдвига удается показать, что параметр к лежит в интервале л2/12 < k < 1 [28]. Однако в теории стержней подобный результат не доказан [29], поэтому, с учетом специфики моделируемого материала, не будем ограничивать данный параметр сверху. Из (46) следует, что

= 4 + 2(1 + у)

СГ 3 d2 к '

(58)

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

Отсюда следует, что при V = 0 должно выполняться соотношение k > 2СВ/СА . При этом зависимость d от параметра k — быстро убывающая монотонная функция.

При k = 2 СВ/СА диаметр стержня устремляется в бесконечность, а минимальный диаметр достигается при k ^ ^ и совпадает со значением, даваемым классической теорией d = 0.122 нм. При k = 4.12 длина стержня сравнивается с его диаметром. Значение изгиб-ной жесткости при этом составляет 0.53 нН-нм. При k =10 диаметр стержня отличается от значения, даваемого классической теорией, менее чем на 5 %.

Если воспользоваться классическим значением коэффициента поперечного сдвига и положить k = п2/12, то для материала стержня, моделирующего углеродную связь, условие (59) положительности диаметра дает следующий интервал возможных значений коэффициента Пуассона:

-1 < V < -0.252. (60)

Заметим, что при V ^ -0.252 диаметр стержня неограниченно возрастает, а при V ^ -1 выражение для изгиб-ной жесткости становится равным полученному ранее с помощью модели Бернулли-Эйлера (0.122 нм).

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

микропараметров (жесткостей связи на изгиб и кручение), например из экспериментальных данных либо из данных, используемых для определения параметров распространенных силовых полей для графитовых структур. На данный момент есть трудности с экспериментальным определением значений жесткости ковалент-ной связи в графене на кручение и изгиб. Поэтому в этой работе предлагается использовать дополнительные предположения о связи между атомами с использованием стержневой модели. Используемая модель позволяет выразить изгибную жесткость графена через поперечную и продольную жесткости связи, отвечающие за деформацию в плоскости листа графена. Значения этих жесткостей можно определить из известных экспериментальных данных для графита. В результате вычислена изгибная жесткость графенового листа. Полученные аналитически значения изгибной жесткости имеют тот же порядок, что и результаты компьютерного моделирования, основанные на моделировании с использованием эмпирических потенциалов и ab initio расчетов.

Работа выполнена при поддержке грантов РФФИ (№№ 12-01-31459, 13-01-90718 и 11-01-00809-a). Авторы благодарят Е.А. Иванову, Д.А. Индейцева и

H.Ф. Морозова за ценное обсуждение результатов, полученных в данной работе.

Литература

I. Ламберт Е. Кварцевые и кремниевые генераторы: Golledge и Silicon Labs // Компоненты и технологии. - 2010. - № 7. - С. 80-82.

2. Geim A.K., Novoselov K.S. The rise of graphene // Nature Materials. -

2007. - No. 6. - P. 183-191.

3. Bunch S.J., van der Zande A.M., Verbridge S.S., Frank I.W., Tanen-baum D.M., Parpia J.M., Craighead H.G., McEuen P.L. Electromechanical resonators from graphene sheets // Science. - 2007. -No. 315.- P. 490-493.

4. Chen C., Rosenblatt S, Bolotin K.I., Kalb W, Kim P., Kymisis I. Performance of monolayer graphene // Nature Nanotechnology. - 2009. -V. 4. - P. 861-867.

5. Peierls R. Remarks on transition temperatures // Helv. Phys. Acta. -

1934. - No. 7. - P. 81.

6. Ландау Л.Д., Лифшиц Е.M. Теоретическая физика. Статистическая

физика. Ч. 1. - М.: Наука, 1976. - 584 с.

7. Тимошенко С.П. Сопротивление материалов. - М.: Наука, 1965. -Т. 2. - 480 c.

8. Blakslee O.L., Proctor G.B., Seldin E.J., Spence G.B., Weng T. Elastic constants of compression-annealed pyrolitic graphite // J. Appl. Phys. -1970. - No. 41. - P. 3373-3382.

9. Lee C., Wei X., Kusar J.W., Honel J. Measurement of the elastic properties and intrinsic strength of monolayer graphene // Science. -

2008. - V. 321. - No. 5887. - P. 385-388.

10. Reddy C.D., Rajendran S., Liew K.M. Equilibrium configuration and continuum elastic properties of finite sized graphene // Nanotechnology. - 2006. - V. 17. - No. 3. - P. 864-870.

11. Sakhaee-Pour A., Ahmadian M.T., Naghdabadi R. Vibrational analysis of single-layered graphene sheets // Nanotechnology - 2008. -V. 19. - P. 085702.

12. Zhang D.B., Akatyeva E., Dumitrica T. Bending ultrathin graphene at the margins of continuum mechanics // Phys. Rev. Lett. - 2011. -V. 106. - P. 255503.

13. Wei Y., Wang B., Wu J., Yang R., Dunn M.L. Bending rigidity and Gaussian bending stiffness of single-layered graphene // Nano Lett. -2013. - V. 13. - P. 26-30.

14. Arroyo M., Belytschko T. Finite crystal elasticity of carbon nanotubes based on the exponential cauchy-born rule // Phys. Rev. B. - 2004. -V. 69. - P. 115415.

15. Kudin K.N., Scuseria G.E. C2F, BN, and C nanoshell elasticity from ab initio computations // Phys. Rev. B. - 2001. - V. 64. - No. 23. -P. 235406.

16. Atalaya J., Isacsson A., Kinaret J.M. Continuum elastic modeling of graphene resonators // Nano Lett. - 2008. - V. 8. - No. 12. - P. 41964200.

17. Lu Q., Arroyo M., Huang R. Elastic bending modulus of monolayer graphene // J. Phys. D: Appl. Phys. - 2009. - V. 42. - P. 102002.

18. Lindahl N. Determination of the bending rigidity of graphene via electrostatic actuation of buckled membranes // Nano Lett. - 2012. -V. 12. - P. 3526-3531.

19. Poot M., van der Zant H.S. Nanomechanical properties of few-layer graphene membranes // Appl. Phys. Lett. - 2008. - V. 92. -P. 063111.

20. Scharfenberg S. Probing the mechanical properties of graphene using a corrugated elastic substrate // Appl. Phys. Lett. - 2011. -V. 98.- P. 091908.

21. Иванова E.A., Кривцов A.M., Морозов Н.Ф., Фирсова А.Д. Учет моментного взаимодействия при расчете изгибной жесткости наноструктур // ДАН. - 2003. - Т. 391. - № 6. - С. 764-768.

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

23. Борн М., Кунь X. Динамическая теория кристаллических решеток. - М.: ИЛ, 1958. - 488 с.

24. ГоловневИ.Ф., Головнева Е.И., Фомин В.М. Расчет термодинамических свойств наноструктур методом молекулярной динамики // Физ. мезомех. - 2007. - Т. 10. - № 5. - С. 71-76.

25. Aлехин В.В., Amm Б.Д., БабичевA.В., Коробейников С.Н. Собственные колебания и выпучивание графеновых листов // Изв. РАН. МТТ. - 2013. - № 5. - С. 34-38.

26. Гольдштейн Р.В., Ченцов A.B. Дискретно-континуальная модель нанотрубки // Изв. РАН. МТТ. - 2005. - № 4. - С. 57-74.

27. Кривцов A.M. Упругие свойства одноатомных и двухатомных кристаллов. - СПб.: Изд-во СПбГПУ, 2009. - 127 с.

28. Жилин n.A. Прикладная механика. Теория тонких упругих стержней. - СПб.: Изд-во СПбГПУ, 2007. - 101 с.

29. Жилин n.A. Прикладная механика. Основы теории оболочек. -СПб.: Изд-во СПбГПУ, 2006. - 167 с.

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

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

Беринский Игорь Ефимович, к.ф.-м.н., зав. отд. ИПМаш РАН, доц. СПбГПУ, iberinsk@gmai1.com Кривцов Антон Мирославович, д.ф.-м.н, проф., зав. каф. СПбГПУ, зав. лаб. ИПМаш РАН, akrivtsov@bk.ru Кударова Асия Муратовна, асп. Делфтского технического университета, asya.kudarova@gmai1.com