УДК 539.3
Вестник СПбГУ. Сер. 1. 2012. Вып. 3
ИССЛЕДОВАНИЕ УСТОЙЧИВОСТИ И СТРУКТУРНОГО ПЕРЕХОДА В ГЦК-РЕШЕТКЕ ПРИ БОЛЬШИХ ДЕФОРМАЦИЯХ*
Е. А. Подольская1, А. М. Кривцов2, А. Ю. Панченко3
1. Институт проблем машиноведения РАН (ИПМаш РАН), магистр механики, мл. научн. сотр., katepodolskaya@gmail.com
2. Институт проблем машиноведения РАН (ИПМаш РАН), д-р физ.-мат. наук, профессор, akrivtsov@bk.ru
3. С.-Петербургский государственный политехнический университет, магистр физики, аспирант, ArtemQT@yandex.ru
Введение. В настоящее время актуальность приобретают задачи, связанные с оценкой прочностных свойств объектов, в силу своей малости обладающих бездефектной кристаллической структурой. Прочность материала тесно связана с его устойчивостью при конечных деформациях [1].
Целью данной работы является исследование устойчивости при больших деформациях идеальной (без дефектов) гранецентрированной кубической (ГЦК) кристаллической решетки. ГЦК-решеткой обладают многие металлы, например, медь, железо (при определенных условиях), серебро, золото, платина. Вопрос устойчивости ГЦК-решетки рассматривался, например, в [2, 3]. В работе [2] показано, что плотноупа-кованные ГЦК- и гексагональная плотноупакованная (ГПУ) структуры устойчивы в малом при любом парном центральном силовом взаимодействии частиц. В работе [3] рассмотрена устойчивость ГЦК- и объемно-центрированных кубических (ОЦК) кристаллов с применением псевдопотенциалов; показано, что потеря устойчивости при сжатии и фазовые переходы связаны с обращением в ноль модуля сдвига. В настоящей работе применяется динамический критерий (вещественность частот упругих волн); с его помощью, например, в работе [4] исследовалась устойчивость нового на-номатериала графена (один слой графита) с использованием парного моментного потенциала [5].
Настоящая работа является продолжением [6], где была исследована устойчивость двумерной треугольной кристаллической решетки — атомного слоя в ГЦК- и ГПУ-структурах — и было показано, что существуют дополнительные области устойчивости, связанные со структурным переходом в материале.
Взаимодействие частиц, составляющих решетку, описывается при помощи парного центрального силового взаимодействия. Это удобная и простая модель для построения теории, проведения аналитических расчетов и вычислительных экспериментов.
В ходе исследования обнаружена возможность структурного перехода ГЦК-решетки в ОЦК-решетку при диагональном тензоре деформации, главные оси которого совпадают с осями кубической симметрии. Кроме того, при наличии сдвиговых деформаций наблюдаются дополнительные области устойчивости, имеющие ту же природу, что и описанные в [6].
*Доклад на Международной конференции по механике «Шестые Поляховские чтения» 31 января— 3 февраля 2012 г., Санкт-Петербург, Россия.
© Е. А. Подольская, А. М. Кривцов, А. Ю. Панченко, 2012
Критерий устойчивости. В работе применяется прямое тензорное исчисление [1]. В качестве закона взаимодействия используются потенциалы Морзе и Леннард-Джонса:
Пм (г) = В
,-2в(т/а-1) _ 2е-в(т/а-1)
Пы (г) = В
12
- 2
(1)
Здесь В — глубина потенциальной ямы, величина параметра в обратно пропорциональна ширине ямы. Вблизи положения равновесия при в = 6 потенциал Морзе эквивалентен потенциалу Леннард-Джонса с теми же значениями глубины потенциальной ямы и равновесного расстояния а [7]. Одно из отличий Пм (г) от П^/(г) состоит в том, что при сжатии материала в точку (г = 0) возникает конечная сила отталкивания, при в = 6 порядка 106В/а. Это позволяет проводить молекулярно-динамическое моделирование при сильном сжатии. Кроме того, быстрое затухание экспонент дает возможность ограничиться меньшим числом частиц. С помощью именно потенциала Морзе можно описать фазовый переход ОЦК^ГЦК, при этом значения параметра в лежат в пределах от 3 до 5 для ряда металлов, у которых наблюдается данный переход [8].
Для описания материала в отсчетной конфигурации введена система координат и положение каждой частицы связано с началом отсчета. Так как ГЦК-решетка простая, любую частицу можно выбрать в качестве нулевой. Частицам присваиваются номера к = ±1... ± N; частицы, расположенные симметрично относительно нулевой, имеют противоположные по знаку номера. Обозначим радиус-векторы, определяющие положения частиц относительно нулевой, через а^, причем а_к = —ак.
Далее будем рассматривать однородную деформацию, описываемую деформаци-
о о
онным градиентом И.У, где V — оператор Гамильтона в отсчетной конфигурации.
Для получения макроскопических уравнений используется длинноволновое приближение [9]. Это означает, что рассматриваются лишь функции, слабо изменяющиеся на расстояниях, сравнимых с длинами основных векторов решетки, то есть волны, длины которых много больше межатомных расстояний.
Связь между векторами Ак (радиус-векторы частицы в актуальной конфигурации) и а к с использованием длинноволнового приближения имеет вид
А к = И (г — а к) — И (г) « а к • VR = ^ • а к.
Запишем уравнение движения сплошной среды в форме Пиолы [1]:
о
рои = V- Р,
(2)
(3)
где ро — плотность материала в отсчетной конфигурации, и — вектор перемещений, Р —тензор напряжений Пиолы, определяющийся формулой (см. [7])
Р
1
— Ф к
П'(Ак) Ак '
(4)
Уо — объем элементарной ячейки
Далее исследуется первая вариация уравнения (3) вблизи деформированного состояния кристаллической решетки. В результате преобразований, более подробно описанных в [10], получаем волновое уравнение
V = 4д • • • vvv, (5)
6
где V = ¿и — первая вариация вектора перемещений, V — оператор Гамильтона в актуальной конфигурации, 4Р — тензор четвертого ранга, зависящий от первой и второй производных потенциала взаимодействия (усилий в связях и жесткостей связей), а также от геометрии окружения частицы [1, 7].
Решением уравнения (5) является V = voeгшíeгKR, К —волновой вектор, ш — частота. Если акустический тензор О = 4Р • •КК положительно определенный, то актуальная конфигурация устойчива [1]. Поскольку ш2 —собственные числа тензора О, то для любого вещественного К должно выполняться условие
ш2 > 0.
Для трехмерной задачи характеристическое уравнение ёе^Б — ш2Е) (Е — единичный тензор) имеет вид
аш6 + Ьш4 + сш2 + й = 0,
а = 1, Ь = — ^Б = —/ь с = 1/2^г2Б — ^Б2) = /2, й = — det Б = —/3,
где /1, /2, /3 —инварианты тензора Б.
Положительность корней кубического уравнения (7) равносильна:
/1 > 0, /2 > 0, /3 > 0, Д > 0,
Д = —4Ь3й + Ь2с2 — 4ас3 + 18аЬсй — 27а2 й2.
(6) 0
(7)
(8)
/ 1+ £х Ч<Рух 0 \
О 0 1+ £у
\ tg^xz 0 1 + ez /
Аффинная деформация ГЦК-решетки. Рассмотрим аффинную деформацию
(9)
Часть внедиагональных элементов полагаются нулевыми, чтобы исключить из рассмотрения эквивалентные деформированные состояния. Пусть оси ж, у, г совпадают с осями кубической симметрии. Предположим также, что tg^>yx = tg^>xz = tg^>zy = 0. Тогда левые части неравенств (8) будут однородными многочленами различных степеней от квадратов компонент волнового вектора К, коэффициенты — функции деформаций ех, £у, е2. Инвариант /1 есть диагональная квадратичная форма, поэтому необходимым условием устойчивости будет положительность коэффициентов этой формы.
Напишем достаточные условия, разбив левые части оставшихся неравенств на квадратичные формы в первых квадрантах. Будем применять метод Монте-Карло там, где, с одной стороны, достаточные условия показывают неустойчивость, с другой, необходимое условие показывает устойчивость. Метод Монте-Карло заключается в том, что для каждой точки из трехмерного пространства начальных деформаций £х, £у, £z составляются неравенства (8) и проверяется их смысл при различных К.
Утверждение. Имеем однородное неравенство Р(ж, у, г) > 0 на положительном октанте. Сделаем замену г = 1 — ж — у и получим неоднородный многочлен ^(ж, у) > 0 при ж> 0, у > 0, х + у< 1.
Рис. 1. Область устойчивости ГЦК-решетки в пространстве £ж, £у, при диагональном тензоре деформации (потенциал Морзе, 0 = 6).
Это утверждение используется для ускорения расчетов методом Монте-Карло с целью обеспечения минимального перебора лучей параметров.
На рис. 1 приведена область устойчивости ГЦК-решетки, частицы которой взаимодействуют посредством потенциала Морзе (1) с параметром в = 6. Рассматриваются три координационные сферы. Видим, что область устойчивости невыпуклая: имеется основная часть, вытянутая вдоль оси £х = £у = ^, и три небольшие области, примыкающие к ней в зоне сильного сжатия. Установлено, что при деформации, реализующейся в этих дополнительных областях, возникает структурный переход из сжатой ГЦК-решетки в ОЦК, сжатой вдоль осей кубической симметрии и повернутой вокруг одной из осей. Это происходит вследствие «выдавливания» некоторых частиц с «п»-й координационной сферы на «п + 1»-ю (рис. 2).
На рис.2 показан пример перехода ГЦК —> ОЦК при ех = еу = \/2/3 — 1, ег = 2/\/3 — 1; частица, принадлежавшая грани ГЦК-решетки, оказывается в центре ОЦК-решетки. На рисунке увеличены частицы, которые присутствуют и на рис. 2, а, и на рис. 2, б. На рис. 1 квадратами показаны три точки неустойчивого равновесия ОЦК-решетки. Неустойчивость связана с выбором потенциала взаимодействия, полученный результат согласуется с [8]. При уменьшении параметра в дополнительные области увеличиваются.
Рис. 2. Переход ГЦК (а) ^ ОЦК (б).
Отдельное исследование показало, что, как и в двумерном случае [6], использование потенциала Леннард-Джонса обеспечивает устойчивость материала при его гидростатическом сжатии, т. е. при его деформировании по линии £х = £у = £г, вплоть до деформаций, сколь угодно близких к точке £х = £у = £г = —1. При этом допо-нительные области, соответствующие ОЦК-решетке, не возникают, что согласуется с[8].
Также решена аналогичная задача деформирования ГЦК-решетки с учетом сдвиговых деформаций. Рассмотрим аффинную деформацию (9). Пусть теперь tg^>yx = 0, tg^xz = 0, tg^>zy = 0. Определитель аффинного преобразования (9) имеет вид
Det = (1 + еж)(1 + £у )(1 + ez) + tg^yxtg^xz tg^zy. (10)
Необходимо ввести ограничения на допустимые значения параметров деформации для того, чтобы аффинное преобразование было собственным (Det > 0), то есть каждая актуальная конфигурация была связана с отсчетной посредством непрерывной деформации. По аналогии с [6] положим все три тангенса положительными.
В ходе исследования областей устойчивости, построенных в шестимерном пространстве деформаций, обнаружены структурные переходы, в результате которых одна из осей кубической симметрии в отсчетной конфигурации становится осью [1,1,1] в актуальной конфигурации.
Молекулярно-динамическое (МД) моделирование. Для проверки полученных результатов использовался метод динамики частиц. Техника моделирования описана в [7]. Для ряда деформированных конфигураций проводился следующий вычислительный эксперимент. В качестве начального условия строилась ГЦК-решетка в деформированном состоянии с периодическими граничными условиями. Взаимодействие частиц осуществлялось посредством потенциала Морзе (1). Начальная кинетическая энергия частиц не превышала 0.0002Д. Эволюция системы описывалась при помощи численного интегрирования уравнений Ньютона методом Верле. Если в процессе эволюции системы наблюдались ограниченные по амплитуде колебания кинетической энергии вокруг некоторого значения, не превышающего 0.0002Д, то делался
Рис. 3. Область устойчивости ГЦК -решетки в пространстве вх, ву, в г при диагональном тензоре деформации, (потенциал Морзе, 0 = 6). Серым цветом показаны точки, полученные при помощи МД-моделирования. Максимальное гидростатическое сжатие составляет 60% при теоретическом подходе и 75% при МД-моделиро-вании.
вывод об устойчивости данной конфигурации. Если наблюдался резкий рост кинетической энергии, то деформированная конфигурация считалась неустойчивой. Области, полученные в результате молекулярно-динамического моделирования (рис.3), совпали с приведенными на рис. 1 в пределах точности компьютерных вычислений. Наблюдающееся расхождение при сильном сжатии связано с тем, что при аналитических расчетах использовалось фиксированное количество координационных сфер, тогда как при МД-моделировании с увеличением сжатия в рассмотрение включалось все больше координационных сфер. Временные затраты на построение зон устойчивости оказались несоизмеримо больше, чем при теоретическом подходе.
Выводы. В ходе исследования устойчивости ГЦК-решетки при больших деформациях выявлено два типа структурных переходов: ГЦК ^ ОЦК при сжатии вдоль осей кубической симметрии, ГЦК ^ ГЦК при добавлении сдвиговых деформаций. При использовании потенциала Леннард-Джонса не наблюдается переход ГЦК ^ ОЦК, материал не теряет устойчивость при сколь угодно большом гидростатическом сжатии. Молекулярно-динамическое моделирование дополняет результаты проведенного исследования в области сильного (60-75%) сжатия.
Литература
1. Лурье А. И. Нелинейная теория упругости. М.: Наука, 1980. 512 с.
2. Wallace D. C., Patrick J. L. Stability of crystal lattices // Phys. Rev. 1965. Vol. 137. N 1A. P. 152-160.
3. Milstein F., Rasky D. Theoretical study of shear-modulus instabilities in the alkali metals under hydrostatic pressure // Phys. Rev. B. 1996. Vol.54. N10. P. 7016-7025.
4. Товстик П. Е., Товстик Т.П. Модель двумерного графитового слоя // Вестн. С.-Петерб. ун-та. 2009. Вып. 3. C. 1-11.
5. Беринский И. Е., Иванова Е. А., Кривцов А. М., Морозов Н. Ф. Применение момент-ного взаимодействия к построению устойчивой модели кристаллической решетки графена // Изв. РАН. МТТ. 2007. №5. С. 6-16.
6. Подольская Е.А., Кривцов А.М., Панченко А.Ю., Ткачев П. В. Устойчивость идеальной бесконечной двумерной кристаллической решетки // Докл. РАН. 2012. Т. 442, №6. С. 755-758.
7. Кривцов А. М. Деформирование и разрушение твердых тел с микроструктурой. М.: Физматлит, 2007. 304 с.
8. Теоретическая механика. Упругие и тепловые свойства идеальных кристаллов: учеб. пособие / И. Е. Беринский и др.; под ред. А.М.Кривцова. СПб.: Изд-во Политехн. ун-та, 2009. 144 с.
9. Born M., Huang K. Dynamical Theory of Crystal Lattices. Oxford: Clarendon Press, 1954. 420 p.
10. Podolskaya E. A., Panchenko A. Yu., Krivtsov A. M. Stability of 2D triangular lattice under finite biaxial strain // Nanosystems: Physics, Chemistry, Mathematics, 2011, 2 (2). P. 8490.
Статья поступила в редакцию 26 апреля 2012 г.