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

Вестник ПНИПУ. Аэрокосмическая техника. 2016. № 44

DOI: 10.15593/2224-9982/2016.44.08 УДК 538.911:539.32

С.С. Стволова, И.Ю. Зубко

Пермский национальный исследовательский политехнический университет, Пермь, Россия

ПОСТРОЕНИЕ ТЕНЗОРА УПРУГИХ СВОЙСТВ АНИЗОТРОПНЫХ МАТЕРИАЛОВ ПРИ ДИСКРЕТНО-АТОМИСТИЧЕСКОМ МОДЕЛИРОВАНИИ

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

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

S.S. Stvolova, I.Yu. Zubko

Perm National Research Polytechnic University, Perm, Russian Federation

CONSTRUCTION OF ELASTIC PROPERTIES TENSOR OF ANISOTROPIC MATERIALS USING DISCRETE-ATOMISTIC

APPROACH

Prediction of mechanical properties of nanostructured materials which are widely used in aerospace industry is based (according to many authors) on discrete atomistic simulation. Such approach often provides a unique way to studying nanomaterials. Interatomic interaction in solids is modeled

by a huge amount of different potentials - pairwise, many-particles potentials, embedded atom method, covalent bond potential etc. It is well known that computed mechanical properties in some cases may differ from experimental data even qualitatively. The paper aims at the demonstration of different potentials ability to explain elastic anisotropy by studying invariant representation of tensor of elastic moduli in the exact form which has been built using different potentials of interatomic interaction. The type of potentials is not specified, their derivatives are expressed through interatomic forces; it makes it possible to study the opportunities of different potentials in order to describe anisotropy and symmetry properties of material elastic response with crystalline microstructure. The authors demonstrated the ability of two-particle or multi-particle potentials of interatomic interaction for description of anisotropy of materials elastic properties using the obtained invariant representation with an example of two-dimensional quasi-crystalline structures. The pairwise potentials in contrast to the many-particle embedded atom potential are shown to be unable to explain elastic anisotropy.

Keywords: discrete-atomistic approach, nonsymmetric elasticity, elastic properties prediction, elastic anisotropy, finite sized crystals, many-particles potentials, embedded atom method, fast sums method.

Введение

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

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

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

радиус-вектор произвольного к-го атома образца в начальной конфигурации, тогда в текущей конфигурации его положение задается вектором гк = К • Як (правило Коши-Борна [1, 2]). Задача определения

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

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

£ = 1(Уи + иУ) малых деформаций, где вектор V - дифференциальный

оператор Гамильтона (градиент); и - вектор перемещений, для которого принимается, что компоненты как самого этого вектора, так и его градиента малы. В этом случае отсчетная К0 и текущая К конфигурации являются достаточно близкими, чтобы отождествить для них векторы локального материального базиса И • =Эк/дХ' е К0

и Г' =дг/дХ' е К, где И = И(XX1,X2, X3) - радиус-вектор частиц среды в отсчетной конфигурации; г = г(х1,X2,X3,1) - радиус-вектор тех

(1 2 3 \ / 1 2 3 \

X1,X2,X3,0) = И(X1,X2,X3). Здесь X' - выбранная (в общем случае криволинейная) система мате-

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

д/дХ' или в текущей V = г' д/дХ' конфигурациях. Основной тензорной характеристикой, связывающей две введенные конфигурации, является аффинор, или деформационный градиент [3]

о

К = г Я' = (V г)т, описывающий аффинное искажение малой окрестности произвольной материальной частицы сплошной среды. Используем

о

вектор перемещений частиц тела и = г - Я, аффинор К = I + (Vи)т зао

писывается через тензор дисторсии в = (V и) , где I - единичный тензор второго ранга. Обратный деформационный градиент К-1 = Я;г' =

= (V Я )т = I - (V и )т. При выполнении условий близости отсчетной и текущей конфигураций можно отождествить векторы локального базиса в этих конфигурациях: г' « Я', векторные операторы Гамиль-

о Л о Л

тона V = V и тензоры дисторсии (Vи)т «(^7и)т = ^и)1'. При этом можно записывать все соотношения для деформируемой среды в терминах внешней лабораторной системы координат или отсчетной конфигурации.

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

а = П: е. (1)

При использовании симметричных мер о и £ тензор П симметричен относительно перестановок внутри первой и последней пар индексов:

П ^ =П = П1)1к. (2)

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

р и = 2 а :е = 2 е :П:е. (3)

Эта величина из термодинамических соображений представляет собой положительно определенную квадратичную форму, т.е. и > 0 при е Ф 0 и и = 0 при е = 0. Отсюда следует дополнительная симметрия П относительно перестановки пар индексов:

П']к1 =П щ. (4)

Компоненты тензора П (упругие модули) в анизотропном случае имеют ясный физический смысл в материальных осях, связанных со структурой материала, например в кристаллографических осях для монокристаллов или в осях, связанных с ориентацией армирующих элементов композиционного материала. Использование тензора малых деформаций подразумевает при замене системы отсчета допустимость наложения только малого поворота системы отсчета, поэтому закон (1) применим лишь для деформируемых тел, не испытывающих (больших) поворотов относительно инерциальной системы отсчета, в которой решается задача. По отношению к малым поворотам тензор £ является инвариантным, поэтому при исследовании больших деформаций, которые могут сопровождаться значительными поворотами осей, задающих материальную симметрию, использование тензора малых деформаций £ и закона (1) даже при малых искажениях материальных волокон некорректно. Для описания больших искажений материального базиса используются тензорные меры, имеющие различный смысл. Например, могут использоваться тензоры, в компонентном представлении описывающие изменение локальной материальной метрики в окрестности частицы среды текущей конфигурации g ~ = г • г;. по отношению к метрике О^ = К ; • И ■ окрестности той же частицы в отсчет-

ной конфигурации: -2 - О ~). В локальном базисе отсчетной конфигурации эта разность дает тензор деформаций Коши-Грина

о 1 / т \

£ = 2(К • К -1), а в текущей конфигурации - тензор деформаций Альманси £ = 2 (I - К - т • К-1). В условиях малых деформаций эти тензоры принимают вид £ = 2 (т + К)-1 и £ = I - 2(-т + К-1), но при том, что в силу малости деформаций справедливо

о о

К-1 = I - (Х7и)т «I - (V и)т = 2! - К, следует равенство £ « £ « £.

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

р и = с: у^7, (5)

где р - локальная плотность среды в текущей конфигурации; V - вектор скорости материальных частиц; V - дифференциальный набла-оператор (градиент) в текущей конфигурации; ^^ = с: vV - мощность напряжений, либо в форме

р и = Р: V V = Р: К, (6)

где р - плотность среды в отсчетной конфигурации; Р - (двухточечный) первый тензор напряжений Пиолы-Кирхгофа, Р = IК-1 • о;

I = ёй К = р/р; V - набла-оператор в отсчетной конфигурации. При выводе выражения (6) мощность напряжений определяется как W(¡) = I _1Р: К. Для скорости изменения внутренней энергии и доказывается ее независимость от выбора системы отсчета. Эта величина может быть представлена с помощью различных пар тензорных мер напряженного и деформированного состояния. Например, (5)-(6) можно записать как

р и = т : уУ,

где т = I а - взвешенный тензор напряжений Кирхгофа. Аналогичным способом может быть введен набор различных пар энергетически сопряженных тензоров напряжений Т и деформаций Е (тензоры правого семейства, т.е. построенные в отсчетной конфигурации с помощью правой меры искажений и), с помощью которых записывается мощность напряжений единицы объема в отсчетной конфигурации Т: ЕЕ = р и = I ). Подобным образом вводятся и энергетически сопряженные пары £ и е (левого семейства, построенные с помощью левой меры искажений V) в текущей конфигурации: р: е = р и = Заметим, что тензор напряжений Коши и тензор деформаций Альманси не являются энергетически сопряженными, и упругий закон в текущей конфигурации не может быть записан в виде линейной связи между ними.

При записи упругого закона в терминах отсчетной конфигурации

0 11 т \

используется тензор деформации Коши-Грина £ = 2(К • К -1), а

энергетически сопряженным ему тензором напряжений оказывается второй тензора Пиолы-Кирхгофа Рп = IК-1 • о • К-т:

о/

I ) = Рп : d £/ ¿1. (7)

Аналогично могут быть получены другие пары энергетически сопряженных тензоров.

Если в качестве кинематического аргумента и рассматривается деформационный градиент К, то и (К) = иК : Кт, следовательно,

( -Рт): Кт = 0, тогда

Р = рйи/№т. (8)

Заметим, что в чистом виде аргумент К входить в представление для энергии деформирования и(К) не может в силу требования инвариантности и по отношению к наложенному жесткому движению. Это замечание вызвано тем, что деформационный градиент преобразуется как К* = О • К, где О - собственно ортогональный тензор, и* = и, т.е. и (О • К) = и(К). Это приводит к требованию, чтобы тензор К

входил в функцию и как Кт • К : и = и(Кт • К), йи/ёТ = = ёи/ё ( К т • К): ((К т • К)/ Ж) . Аналогично при выборе в качестве кинематического аргумента тензора деформаций Коши-Грина получается выражение для второго тензора напряжений Пиолы-Кирхгофа:

Рп =рёи /ё£. (9)

Такая запись упругого закона соответствует определению гиперупругого материала (упругость по Грину). В физически линейном случае (9) дает упругий закон в виде

Рп = Н :£, (10)

где Н - тензор линейно-упругих свойств, определенный в отсчетной конфигурации (совпадает с тензором П в выражении (1)). Для тензора напряжений Коши справедливо

о = I-1 К • рп(£) • Кт. (11)

Подставляя (10) в (11), получим формулировку упругого закона в К:

о = I"1 К • (Н : £) • Кт. (12)

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

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

Форма несимметричного упругого закона для наночастиц

Второй тензор напряжений Пиолы-Кирхгофа Рп, определенный в отсчетной конфигурации, при малых упругих искажениях связан

о

с линейным тензором деформаций £ - £ = ( + Кт )/2-1 обобщенным

законом Гука (10). Построим обобщение (10) на случай несимметричных мер напряженного и деформированного состояния. В частном случае это обобщение должно сводиться к симметричной форме (10). В симметричном случае Н= Н = Н и закон (10) при малых деформациях можно представить в виде

Рп = Н : (( + Кт )/2-1) = Н :(К -1 + А), (10')

где К -1 = в - тензор «малой» дисторсии; А - любой «малый» косо-симметричный тензор второго ранга, описывающий изменение конфигурации образца при наложении на него малого поворота. Для несимметричного случая будем использовать следующую форму закона Гука (здесь подразумевается, что искажения и повороты малы) для несимметричного внутри пар индексов тензора Н:

Рп = Н : (К-1- w#), (13)

где w # - кососимметричный тензор малого поворота кристаллической

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

О -I + ю,

К*-1*= О • К -1 - (I + ю) • К -1 = (I + ю) • (К -1) + ю - (К -1) + ю,

в итоге справедливо соотношение К*-1* « (К -1) + ю, w # = w # + ю. Окончательно получаем К*-1* - w # « К -1 - w #, т.е. введенная несимметричная мера малых деформаций является инвариантной. При дополнительном наложении на тензор линейно-упругих свойств условия симметрии И~ы = Н .]Ы = Н ~1к закон (13) сводится к обычному

виду (10).

Мера деформаций в правой части (13) соответствует мере, предложенной для больших деформаций в работе [5], в которой при учете вращательной степени свободы частиц среды обосновано использование объективной лагранжевой меры Е = Qт • К -1, где тензор Q описывает вращение ориентационного триэдра, который для среды Коссе-ра характеризует поворот материальной частицы. При ортогональном

преобразовании О тензор Е преобразуется как Е* = О • Е • От. В частности, тензор Q может описывать ориентацию кристаллической решетки относительно фиксированного базиса пространственной системы координат: Q = Q#, г = (т • К)• Я. При малых искажениях и поворотах мера Е сводится к используемой в формуле (13) мере К -1 - w#:

Q#т • К -1-(I + w # )т • (I + в) -1 = (1 - w # )• (I + в) -1 = = в - ^^# - ^^# •в - в - ^^# = Г -1 - w#.

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

и = ^ (К -1 - w#): Н : (К -1 - w#). (14)

Плотность внутренней энергии и (упругий потенциал) из термодинамических соображений представляет собой положительно определенную квадратичную форму, и > 0 при любых К -1 - w# Ф 0, и и = 0 только при нулевых искажениях. Отсюда следует симметрия компонент тензора упругих свойств Н относительно перестановки первой и последней пар индексов Н Г]Ш = НЫ}!. При этом не предполагается априорной симметрии тензора напряжений Коши, т.е. согласно структуре соотношения (14) из него не следует симметрия внутри пар индексов

компонент тензора линейно-упругих модулей, т.е. в общем случае

Н1]к1 ^ Н ¡1Ы ^ Н1]1к.

При разложении тензора дисторсии в на сумму в = £ + w симметричной £ и кососимметричной w частей из (13) следует Рп = Н : £ + Н : ( - w#). В тензоре четвертого ранга Н можно выделить

полусимметричную часть Н(4), для которой выполняются условия симметрии в парах индексов Н{щ ) = Н(к) = Н\-к), «полуантисиммет-

Н( а ,а) ут( а,а) уу (а,а) уу (а,а)

, Нж = -Нцк/ = ~Нтк и две смешанные

ijkl jikl ijlk

1 ijkl

части H(s,a) и H(a,s-1, для которых выполняются условия H(kf) = = H = - H^\ Ha) = -H ^ ) = j\ Подставляя разложение H = = H() + H(s,a ) + H(a,s ) + H(a,a) в упругий закон (13), получим

Pn = H(s,s): s + H(a,s>: £ + H(a,a) : (w - w#) + H(s,a) : (w - w#), (15)

следовательно, несимметричный тензор напряжений разделяется на симметричную P^) = H(s,s): £ + H(s,a): (w - w#) и антисимметричную

PI(Ia) = H(a,s) : s + H(a,a): (w - w#) части. Плотность упругой энергии в «несимметричном случае» представляется в виде

u = 2£ : H s,s) : £ +1 (w - w#): H(a,a) : (w - w#). (16)

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

s : H(s ,a): (w - w # ) = 0, (w - w#): H(a,s): s = 0,

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

Из представления (16) следует, что части H(s,s) и H(a,a), составляющие тензор линейно-упругих свойств H, являются положительно определенными. Сумма тензоров H(s,a) и H(a,s) также является положительно определенной. Модули H(s,s) при малых искажениях совпадают с обычными упругими модулями, поэтому идентификация параметров потенциала при дискретно-атомистическом моделировании может быть проведена при сравнении компонент тензора H(s,s) с известными из экспериментов упругими свойствами материала.

Энергетический подход при дискретно-атомистическом

моделировании

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

М М-1 м

Ф = Е ф (Я (•) )+ЕЕ Ф2 (Я (,)) +

^=1 ^=1 ¡=¿+1 (17) М -2 М-1 М /

+ ЕЕ Е ф3 (г к( к^ )+...,

•=1 ¡=1+1 к= ¡+1

где Я(•) - радиус-вектор, задающий положение ¿-го атома, Я(¡) = = Я(]) - Я(•); М - полное число атомов образца; Ф1 (К•)) - часть потенциальной энергии атомов, которая не зависит от их взаимодействия, а определяется полем некоторой внешней силы; Ф2 (Н(•)) - потенциальная энергия парного взаимодействия или двухчастичный потенциал; Ф3 (Я(;7), Я(и), Я^)) - трехчастичный потенциал взаимодействия

атомов. Двухчастичный потенциал межатомного взаимодействия характеризует изменение потенциальной энергии при изменении расстояния между парами атомов. Этот потенциал с помощью некоторой функциональной зависимости описывает, что при сближении два атома начинают отталкиваться, а при удалении притягиваться. Поскольку атомы нельзя сдвинуть бесконечно близко, то в ноле такая функция стремится к бесконечности. При увеличении расстояния между парой атомов эта функция выходит на горизонтальную асимптоту, а сила взаимодействия, модуль которой равен тангенсу угла наклона касательной к графику функции Ф2 (к(¡¡}), стремится к нулю. Трех-

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

потенциалов аддитивности сил взаимодействия нет. В рассматриваемом случае действие внешних сил не рассматривается, поэтому далее

Ф1 ( (•) ) =

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

Для расчета упругих модулей кристаллических систем в рамках энергетического подхода в статической постановке [8-10] принимается, что плотность упругой энергии и плотность потенциальной системы взаимодействующих атомов кристалла в текущей конфигурации совпадают. Тогда производные от плотности потенциальной энергии дискретной системы атомов по мерам деформации дадут выражения для вычисления компонент тензоров напряжений (с помощью первых производных) и компонент тензора линейно-упругих свойств (вторые производные). Пусть Ф(Е) - полная потенциальная энергия однородно деформированного кристаллического образца в текущей конфигурации. Энергия образца в отсчетной конфигурации, которая соответствует системе атомов с заданной кристаллической структурой и минимальным значением потенциальной энергии по параметрам решетки, 0 Л

равна Ф = Ф(1). Относя изменение полной потенциальной энергии де-

0

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

Р = й-1 Э(Ф(Е) -Ф)/ЭЕТ =й-1 ЭФ(Е)/ЭЕт. (18)

Тензор напряжений Коши о находится как о = I 1К • Р1, где

I = ёе К = ££ / й, т.е. тензор напряжений Коши о = й -1 К -ЭФ(К)/ЭК т.

Рассмотрим два вида представления потенциальной энергии в случае учета только двухчастичного взаимодействия, используя для простоты вместо набора векторных аргументов, соединяющих различные пары атомов, обозначение ДК в отсчетной конфигурации или Дг в текущей:

1) Ф(К) = X Ф(Дг) = X Ф(К •ДК), где X (•) - упрощенное обозначение суммы потенциалов (17) в текущей конфигурации;

2) Ф(Е) = X Ф (|Дг|) = X Ф (к • ДК|) = X ДК • Кт • К ДК).

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

Производная по тензору деформационного градиента К в этих двух случаях приводит к следующим выражениям:

1. Производная Ф (Дг) по векторному аргументу Дг дает вектор силы Г (Дг), который определяется конкретным потенциалом межатомного взаимодействия,

о .

й Р = ЭФ(Е)/ЭК т = X Ф '(Дг) •Э(Е •ДК)/ ЭЕТ = X ДКГ (Дг), о = й -1 К •Х ДКГ (Дг) = й -1 X ДгГ (Дг),

Рп = £ -1 X ДК (к-1 Г (Дг)). (19)

Для двухчастичных потенциалов векторы Г (Дг) и Дг коллине-арны, что приводит к симметрии вычисляемых тензоров напряжений (19). Для многочастичного взаимодействия (как в методе погруженного атома для металлов или для ковалентной связи в углеродных материалах) направление силы Г (Дг) может не совпадать с направлением вектора Дг , соединяющего атомы, что приведет к нарушению несимметрии тензора напряжений Коши.

2. Производная Ф'(|Лг|) по скалярному аргументу является скалярной величиной:

& Р = ЭФ(К)/дКт = X ф'(|Дг|) д^ДЯ ■ Кт • К -ДК/ЭКт =

= X Фф '(М )М-1 ДК Дг,

следовательно, для тензора напряжений Коши получим выражение

о = &-1 X Ф'(|Дг| )Дг Дг/|Дг|, Р11 = & -1 X Ф'(|Дг| )ДК ДК/|Дг|, (20)

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

При учете многочастичного взаимодействия в методе погруженного атома [11-13] потенциальная энергия системы атомов определяется выражением, в котором учтено, что отталкивание всех атомов описывается согласно закону парного взаимодействия, а притяжение описывается нелинейной функцией у(.), задающей влияние окружения произвольного атома (функция погружения). Окружение определяется множеством номеров атомов, участвующих в формировании электронного газа вокруг г-го атома:

ад=хг Xф*М+X М г(| „ Ф-(г( я/ (21)

Функции ф± (г( г)) могут зависеть как только лишь от расстояний

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

слагаемое при использовании обозначений (ф- («))) = Г/ дает

Л^(XXх,Ф- () // = *х,Ф- (%) ))Xх, К(«й)•

Р"'{X М-1 X И

'(Ч)

" К(ч)г(ч)+X,=1 (у'X«

/е X,- К(чЛч),

• =а{I» Iи(ф+) |т„>г%л„)+1М( I%Г,)}. (22)

р.=* -'{х:,41 ",+, (Ф+)' +1М, ( I ,, )} (23)

где Г(„) ^ Г 1 • ). Поскольку вектор Г, (или Г(„)) не обязательно

направлен вдоль вектора Г(„) (или «„)), то получаемые тензоры (22),

(23) могут быть несимметричными.

Для вычисления упругих модулей в несимметричном случае при малых деформациях необходимо найти производную от тензора напряжений (23):

й 2ы /йГ 2 = йРт/ йГ =

=°-' л [I м-1 х :=,+' и к,) г1 Г(,)«(,)М (у' I „, Г(-„) к (,) ))

Рассмотрим производные этого выражения по отдельности:

-1

-3 Г/ л'

г( „) (ф+) г( „)

V

Э (ф+)' Г( „)"' Г( „) К („)/Эг =

М ] Г( „)К („ )Г( „)К („) +М'

г( „)

-1

екк („ )е'к („V

где ек - векторы фиксированного пространственного базиса,

1 ,,, к („ )/дГ = = У" (I ы к(„)Г(-)) (I * * к(„)Г(-)) + ^I* * §(„) • екк(„/к(„,

где производные у' и у" функции у имеют аргумент I ф-(г(„)), а тензор второго ранга § {„) =ЭГ(-Л/ дГ( „) =д2ф-(Г( „•))/дг2^. В итоге полу-

чим тензор четвертого ранга:

^дРп/дГ = I ■

1< .<, < М

-3 Г/ +\"

г(„) (ф+) г(„)

V

г( „)К ( „ )г( „)К („) +

+

И

%)

-1

е»Кс-)е»Ксе) \ +

+

ЕМ=1 {^ (Е * Кс-) (Е * Кс-) ) + У Е * 3, (с-) ■ е»КС-)е»

Я

с-)

При малых деформациях К ^ I инвариантное представление тензора линейно-упругих модулей для кристаллического материала в отсчетной конфигурации имеет вид

н = ^-1 е

1<К]< м

-3

К с ч) (ф+) V К с а)

+

И

Я

с-)

-1

екЯс-/Яс-) \ +

с24)

+

+ ^-1 Е М1Е * 3, К(Е * 3, К с-)Е-) )( 3, К с-) Е-)

+ у Е^(Ссе) ■ е»Яс-)екЯс-))},

где С с-) = Э2ф- (К с--) ))ЭК 2„), Е- = ^(К-) 0У)- В общем случае выражение с24) не обладает симметрией внутри пар первых и вторых диад, которая возможна только в случае коллинеарности векторов Е-и Я с-). Без учета погружения ус х) = -1/2 х:

н = ^ Е К«I-3 [(ф+)К)|>+)']КсА)КсА-) +

1<1 < ¡< м

+

Кс«Г1 е»-Кс« \-

И ,

о

"-1 Е {С с -) ■ е»Я с-)е»Я с -)}.

1<1< ] < м

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

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

пользуется метод погруженного атома [11-13], основанный на применении потенциалов многочастичного взаимодействия. Для того чтобы группы атомов взаимодействовали на большом расстоянии согласно экспериментальным законам (описываются потенциалами Морзе или Ми [14]), а на малых расстояниях учитывалась металлическая связь, предлагается модификация метода, основанная на потенциале Ми:

Ф в

1

т - п

М-1 м

ЕЕ

;=1 ]=;+1

т Е 2 Е

а - т (1 - са) а

V Г( ¡0 V Г( ¡0

М

Е с а

( \пр а

]=1 ] *'

г( ¡о

1р'

(25)

где М - число всех атомов образца, т,пе Ъ, с у е {0,1}, р е{1,2}, Су = 1 для атомов, участвующих в образовании электронного газа вблизи положения г-го атома, су = 0 для всех остальных атомов. Раз-

Ч

мер окрестности для учета соседей может быть от одной до нескольких координационных сфер. При р = 1 и любом числе атомов или для двух изолированных атомов М = 2 и любого р потенциал (25) совпадает с потенциалом Ми:

Ф

"в"

1

М -1 м

Е ЕЕ

т - п г=1 у=г+1

а - т а

V Г( ¡0 V Г( ¡0 )

(26)

При условии, что для всех атомов Су = 0 также получается потенциал Ми, при р = 2 получается модификация потенциала Финни-са-Синклэра [12]. При одновременном выполнении условий, что для всех атомов с у = 1 и р = 2 получается модификация [13]. Параметр а

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

Для плоских квазикристаллических структур с осями симметрии различного порядка (примеры - на рисунке) получено, что при использовании парных потенциалов для оси симметрии любого порядка кро-

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

€ ФФФ £ © С» €€€ С

X

а

••У

X б

ссс сссс

СССС <Ь

сссс ссс

Се®®

б

г

Рис. Элементы плоских квазикристаллических структур с осями симметрии: а - 4-го порядка; б - 5-го порядка; в - 6-го порядка; г - 7-го порядка

Потен- Ось симметрии

циал 4-го порядка 5-го порядка 6-го порядка 7-го порядка

С = С = П111 2222 С = С = 1111 2222 С = С = 1111 2222 С = С = 1111 2222

Потен- = 30,06 в/а2 = 26,16 в/а2 = 29,39 в/а2 = 28,32 в/а2

циал Ми С = С = 1122 1212 С = С = 1122 1212 С = С = 1122 1212 С = С = 1122 1212

m = 5, п = 3 = -2,02 в/а2 = 8,72 в/а2 = 9,80 в/а2 = 9,44 в/а2

Е = 29,92 в/а2 Е = 23,25 в/а2 Е = 26,13 в/а2 Е = 25,17 в/а2

V = -0,07 V = 1/3 V = 1/3 V = 1/3

С = С = 4111 "-2222 С = С = 1111 2222 С = С = 1111 2222 С = С = 1111 2222

Потен- = 65,42 в/а2 = 61,74 в/а2 = 76,35 в/а2 = 58,62 в/а2

циал Ми С = С = 1122 1212 С = С = 1122 1212 С = С = 1122 1212 С = С = 1122 1212

m = 12, п = 6 = -3,61 в/а2 = 20,58 в/а2 = 25,45 в/а2 = 19,54 в/а2

Е = 65,22 в/а2 Е = 54,88 в/а2 Е = 67,87 в/а2 Е = 52,10 в/а2

V = -0,05 V = 1/3 V = 1/3 V = 1/3

Потен- С = С = 1111 2222 С = С = 1111 2222 С = С = 1111 2222 С = С = 1111 2222

циал погруженного = 8,85 в/а2 С1122 = 5,92 в/а2 = 8,06 в/а2 С1122 = 7,16 в/а2 = 8,10 в/а2 С1122 = 7,95 в/а2 = 8,60 в/а2 С1122 = 6,40 в/а2

атома (25) т = 12, С1212 =-0,04 в/а2 Е = 4,88 в/а2 С1212 = 0,45 в/а2 Е = 1,70 в/а2 С1212 = 0,08 в/а2 Е = 0,30 в/а2 С1212 = 1,10 в/а2 Е = 3,84 в/а2

п = 6 V = 0,67 V = 0,89 V = 0,98 V = 0,74

Из условия положительной определенности тензора линейно-упругих свойств Н для двумерной среды, как в рассмотренных примерах, следуют ограничения E> 0, уе (-1; 1). Причем модуль Юнга определяется как E = (1^ -Н?122)/НШ1, коэффициент Пуассона

V = H1Ш/H1111. При использовании потенциала погруженного атома

(25) для различных двумерных структур получалось число независимых модулей, которое частично соответствует результатам линейной теории упругости [15]. Отличие от классических результатов теории упругости для двумерной среды дает случай оси симметрии 6-го порядка, который должен описываться только двумя независимыми ненулевыми упругими модулями. Метод погруженного атома при этом дает три независимых компоненты тензора упругих свойств. Также для всех рассмотренных случаев при использовании потенциала (25) полу-

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

Заключение

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

Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (гранты № 14-01-00069 и № 15-01-08678) и в рамках задания № 2014/152 на выполнение государственных работ в сфере научной деятельности в рамках базовой части госзадания Минобрнауки РФ (код проекта - 1911).

Библиографический список

1. Arroyo M., Belytschko T. Finite crystal elasticity of carbon nano-tubes based on the exponential Cauchy-Born rule // Phys. Rev. B. - 2004. -Vol. 69. - Р. 115-126. DOI: 10.1103/PhysRevB.69.115415

2. Reddy C.D., Rajendran S., Liew K.M. Equilibrium configuration and continuum elastic properties of finite sized graphene // Nanotechno-logy. - 2006. - Vol. 17. - P. 864-870.

3. Поздеев А. А., Трусов П.В., Няшин Ю.И. Большие упругопла-стические деформации: теория, алгоритмы, приложения. - М.: Наука, 1986. - 232 с.

4. Zubko I.Yu., Kochurov V.I. Estimation of elastic moduli of graphene monolayer in lattice statics approach at nonzero temperature // AIP Conf. Proc. - 2015. - Vol. 1683. - 020241.

5. Pietraszkiewicz W., Eremeyev V.A. On natural strain measures of the nonlinear micropolar continuum // International Journal of Solids and Structures. - 2009. - Vol. 46. - P. 774-787.

6. Clayton J. Nonlinear Mechanics of Crystals. - London: Springer, 2011. - 715 p.

7. Беринский И.Е., Кривцов А.М. Об использовании многочастичных межатомных потенциалов для расчета упругих характеристик графена и алмаза // Известия РАН. Механика твердого тела. - 2012. -№ 6. - С. 60-85.

8. Симонов М.В., Зубко И.Ю. Определение равновесных параметров решетки различных ГПУ-монокристаллов с помощью потенциала межатомного взаимодействия Ми // Вестник Пермского национального исследовательского политехнического университета. Механика. - 2012. - № 3. - С. 204-217.

9. Зубко И.Ю., Симонов М.В. Энергетический способ расчета упругих модулей образцов конечных размеров с ГПУ-решеткой // Известия Том. политехн. ун-та. - 2013. - Т. 323, № 2. - С. 194-200.

10. Зубко И.Ю. Вычисление упругих модулей монослоя графена в несимметричной постановке с помощью энергетического подхода // Физическая мезомеханика. - 2015. - Т. 18, № 2. - С. 37-50.

11. Daw M.S., Baskes M.I. Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals // Physical Review B. - 1984. - Vol. 29, № 12. - P. 6443-6453.

12. Finnis M.W., Sinclair J.E. A simple empirical N-body potential for transition metals // Philosophical Magazine A. - 1984. - Vol. 50:1. -P. 45-55.

13. Sutton A.P., Chen J. Long-range Finnis-Sinclair potentials // Philosophical Magazine Letters. - 1990. - Vol. 61, iss. 3. - P. 139-146.

14. Israilishvili J.N. Intermolecular and surface forces. - Academic Press: Harcourt Brace and Company, 1998. - 450 p.

15. Черных К.Ф. Введение в анизотропную упругость. - М.: Наука, 1988. - 190 с.

References

1. Arroyo M., Belytschko T. Finite crystal elasticity of carbon nano-tubes based on the exponential Cauchy-Born rule. Phys. Rev. B, 2004, vol. 69, pp. 115-126. DOI: 10.1103/PhysRevB.69.115415

2. Reddy C.D., Rajendran S., Liew K.M. Equilibrium configuration and continuum elastic properties of finite sized graphene. Nanotechnology, 2006, vol. 17, pp. 864-870.

3. Pozdeev A.A., Trusov P.V., Nyashin Yu.I. Bolshie uprugoplas-ticheskie deformatsii: teoriya, algoritmy, prilozheniya [Finite elastic-plastic deformations: theory, algorithms, application]. Moscow: Nauka, 1986. 232 p.

4. Zubko I.Yu., Kochurov V.I. Estimation of elastic moduli of graphene monolayer in lattice statics approach at nonzero temperature. AIP Conf. Proc., 2015, vol. 1683, 020241.

5. Pietraszkiewicz W., Eremeyev V.A. On natural strain measures of the nonlinear micropolar continuum. International Journal of Solids and Structures, 2009, vol. 46, pp. 774-787.

6. Clayton J. Nonlinear Mechanics of Crystals. London: Springer, 2011. 715 p.

7. Berinskiy I.E., Krivtsov A.M. Ob ispolzovanii mnogochastich-nykh mezhatomnykh potentsialov dlya rascheta uprugikh kharakteristik grafena i almaza [On application of multi-particle potentials for computation of elastic characteristics of graphene and diamond]. Izvestiya Rossiyskoy akademii nauk. Mekhanika tverdogo tela, 2012, no. 6, pp. 60-85.

8. Simonov M.V., Zubko I.Yu. Opredelenie ravnovesnykh paramet-rov reshetki razlichnykh GPU-monokristallov s pomoshchyu potentsiala mezhatomnogo vzaimodeystviya Mi [Finding equilibrium lattice parameters of different HCP-monocrystals with use of Mie interatomic potential]. Vest-nic permskogo natsionalnogo issledovatelskogo politekhnicheskogo univer-siteta. Mekhanika, 2012, no. 3, pp. 204-217.

9. Zubko I.Yu., Simonov M.V. Energeticheskiy sposob rascheta uprugikh moduley obraztsov konechnykh razmerov s GPU-reshetkoy [Energy-based approach to estimation of elastic moduli of finite sized specimens with HCP-lattice]. Izvestiya Tomskogo politekhnicheskogo univer-siteta, 2013, vol. 323, no. 2, pp. 194-200.

10. Zubko I.Yu. Vychislenie uprugikh moduley monosloya grafena v nesimmetrichnoy postanovke s pomoshchyu energeticheskogo podkhoda [Computation of elastic moduli of graphene monolayer in non-symmetric formulation using energy-based approach]. Fizicheskaya Mezomekhanica, 2015, vol. 18, no. 2, pp. 37-50.

11. Daw M.S., Baskes M.I. Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals. Physical Review B, 1984, vol. 29, no. 12, pp. 6443-6453.

12. Finnis M.W., Sinclair J.E. A simple empirical N-body potential for transition metals. Philosophical Magazine A, 1984, vol. 50:1, pp. 45-55.

13. Sutton A.P., Chen J. Long-range Finnis-Sinclair potentials. Philosophical Magazine Letters, 1990, vol. 61, iss. 3, pp. 139-146.

14. Israilishvili J.N. Intermolecular and surface forces. Academic Press, Harcourt Brace and Company, 1998. 450 p.

15. Chernykh K.F. Vvedenie v anizotropnuyu uprugost [Introduction in anisotropic elasticity]. Moscow: Nauka, 1988. 190 p.

Об авторах

Стволова Софья Сергеевна (Пермь, Россия) - студентка, ФГБОУ ВПО ПНИПУ (614990, г. Пермь, Комсомольский пр., д. 29, e-mail: [email protected]).

Зубко Иван Юрьевич (Пермь, Россия) - кандидат физико-математических наук, доцент кафедры «Математическое моделирование систем и процессов» ФГБОУ ВПО ПНИПУ (614990, г. Пермь, Комсомольский пр., д. 29, e-mail: [email protected]).

About the authors

Sofya S. Stvolova (Perm, Russian Federation) - Student, Perm National Research Polytechnic University (29, Komsomolsky av., Perm, 614990, Russian Federation, e-mail: [email protected]).

Ivan Yu. Zubko (Perm, Russian Federation) - Ph. D. in Physical and Mathematical Sciences, Associate Professor, Mathematical Modelling of Systems and Processes Department, Perm National Research Polytechnic University (29, Komsomolsky av., Perm, 614990, Russian Federation, e-mail: [email protected]).

Получено 18.11.2015