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

Полный текст

Введение
Появление заметной трещины на макроскопическом уровне является следствием распространения
трещины в нескольких масштабах, сильно отличающихся по длине. Наименьший масштаб — это раз-
рыв связей между атомами по мере распространения трещины. Движущей силой, стоящей за разрывом
соединения, является поле напряжений, обусловленное нагрузкой, приложенной к образцу. Поскольку
материал разрушается, когда концентрация напряжений превышает предел прочности материалов на
разрушение, анализ разрушения напрямую связан с концентрацией напряжений или деформаций во-
круг вершины трещины. Механическая реакция материалов, подверженных экстремальному приложен-
ному напряжению, контролируется атомистическими механизмами вблизи концентраторов напряжений,
таких как вершины трещин. В настоящее время атомистическое моделирование оказалось полезным
инструментом для изучения свойств материалов, при разрушении путем исследования атомистической
конфигурации и поля напряжений вокруг вершины трещины с атомистической точки зрения. В послед-
нее время пристальное внимание уделяется атомистическим напряжениям, обусловленным разрушени-
ем кристаллографической решетки материала вдоль границ зерен. Проведенные исследования показали,
что усилия, возникающие вдоль границ зерен, вычисленные с помощью компонент тензора атомистиче-
ских напряжений, могут служить количественными дескрипторами границ зерен [1]. Дескрипторы да-
ют возможность предсказать напряжение, соответствующее критической эмиссии дислокаций и свойства
абсорбции или пропускания непосредственно из атомистической структуры границы зерен [1–3]. Кон-
цепция атомистических напряжений в настоящее время широко применяется для изучения разрушения
на атомистическом масштабе, а также для того, чтобы пролить свет на поведение таких материалов
как металлическое стекло, гранулированные материалы на атомном уровне [4–6]. Это связано с тем,
что атомистическое напряжение может быть широко интерпретировано как напряжение, эквивалент-
ное напряжению Коши, которое позволяет срастить между собой атомистические явления и процессы,
описываемые механикой сплошных сред [4–6]. В моделировании классической молекулярной динамики
(MD) взаимодействие между атомами определяется межатомным потенциалом, и атомистические напря-
жения могут быть непосредственно вычислены из теоремы вириала [7]. Следует отметить, что понятие
атомистического напряжения все еще актуально за пределами области MD, однако его нелегко опреде-
лить количественно. Например, в расчетах теории функционала плотности (DFT) глобальное (макроско-
пическое) квантово-механическое напряжение может быть получено из изменения полной энергии как
функции равномерной деформации [8; 9]. Однако локальное изменение напряжения, или напряжение,
приходящееся на атом, внутри суперячейки DFT не может быть однозначно определено количественно,
несмотря на то что межатомные потенциалы в MD часто калибруются по DFT [10–13]. Определение ло-
кального квантово-механического напряжения, приходящееся на атом в DFT, с помощью аналитических
или численных методов в настоящее время является активной областью исследований [13–15]. Аналогич-
ным образом понятие атомистического напряжения недостаточно четко определено в экспериментах по
просвечивающей электронной микроскопии высокого разрешения (HRTEM). Однако достижение субна-
нометрового разрешения измерений внутренней деформации может стать реальностью благодаря новым
методам, таким как наноразмерная голографическая интерферометрия [16–18] или использование дви-
жущихся дислокаций в качестве механических зондов [19].
В настоящее время текущие результаты [20–26] получены в специальной атомистической конфигура-
ции, когда трещина лежит в плоскости, параллельной одной из плоскостей симметрии [21–26]. Эволюция
микроструктуры и процесс разрушения монокристаллических материалов могут быть иными, если изме-
нить ориентацию кристалла. Кроме того, атомный потенциал и ограничивающие условия также могут
влиять на процесс атомарного разрушения и эволюцию внутренней микроструктуры. Чтобы понять бо-
лее общий процесс разрушения материала и микромеханизмы монокристаллов на атомарном уровне,
необходима более подробная и расширенная модель MD.
Для восполнения имеющегося пробела в настоящей статье изучена пластина из монокристалличе-
ской меди, ослабленная наклонной трещиной. Трещина наклонена под различными углами по отно-
108
Мушанкова К.А., Степанова Л.В. Опыт моделирования наклонных трещин в материалах с кубической ...
Mushankova K.A., Stepanova L.V. Experience in modeling inclined cracks in materials with cubic crystal ...
шению к главным осям симметрии кристаллической решетки материала с кубической симметрией его
упругих свойств. Вычисления выполнены с помощью метода молекулярной динамики, реализованного
в программе Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS). Ведущей целью вы-
числений было сравнение угловых распределений компонент тензора напряжений, ассоциированных с
вершиной трещины, полученных методом молекулярной динамики, и полученных с помощью аналити-
ческого решения задачи, построенного с использованием методов теории упругости анизотропного тела
[27–29].
1. Моделирование наклонного дефекта в материалах с кубической
кристаллической решеткой
Для исследования полей напряжений в окрестности вершины трещины на наноразмерном уровне был
написан скрипт для LAMMPS. В рамках поставленной задачи построения пластинки с центральным
дефектом на наноскопическом уровне удобно оперировать с размерностями ˚A для длины, пс для време-
ни и бар для напряжения, поэтому применялась команда units metal. Несмотря на то что исследумый
объект — пластинка, все же она будет обладать некоторой толщиной в несколько атомов, поэтому в
коде прописывалась размерность пространства dimension 3. Задавались периодические граничные усло-
вия во всех направлениях boundary p p p. Определение стиля представления частиц atom_style atomic
означает, что будет рассматриваться простейший случай, когда каждая частица в системе считается
точечной без внутренней структуры или ориентации. Задание решетки (в данном случае гранецентри-
рованная кубическая решетка с параметром решетки 3.615 ˚A) определяется командой lattice fcc 3.615.
Создание бокса для размещения атомов размерами 120×120×3 в единицах параметра решетки region
box prism 0 120 0 120 0 3 0 0 0 (часто бывает удобрым моделировать именно призму, поскольку в слу-
чаях, когда к пластине прикладывается сдвиговая нагрузка, происходит наклон ее боковых границ, в
таких ситуациях region box block выдает ошибку). В результате было смоделировано 172 800 атомов.
Для создания области с 3 типами атомов в указанном регионе используется команда create_box 3 box.
Команда create_atoms 1 region box создает атомы типа 1 внутри региона box. Для создания трещи-
ны длины 10 (в единицах параметра решетки) определялось два региона, представляющих верхнюю и
нижнюю области, где были определены атомы типов 2 и 3. На рис. 1.1 атомы первого типа обозначе-
ны цифрой 1, второго — цифрой 2, третьего — 3. Регионы выделялись из соображений, чтобы трещина
располагалась под углом 45 градусов как к осям симметрии рассматриевомого материала с ГЦК-решет-
кой, так и по отношению к прикладываемой далее нормальной нагрузке. Потенциал взаимодействия
между атомами устанавливался командой pair_style eam, где eam — embedded atom method, то есть
использовался потенциал внедренного атома. Командой pair_coeff * * Cu_u3.eam задавались коэфици-
енты потенциала, определенные в файле Cu_u3.eam. Чтобы исключить взаимное влияние друг на друга
атомов типа 2 и 3, использовалась команда neigh_modify exclude type 2 3. Напряжения на атом рассчи-
тывались командой compute peratom all stress/atom. Для моделирования прикладываемой нормальной
нагрузки в 80 ГПа использовалась команда fix npt, тепмература поддерживалась равной 10 К для умень-
шения влияния пластических деформаций. Командой dump формировался файл, позволяющий получить
данные временной эволюции полей напряжений, воспринимаемые пакетами для визуализации. Кадры,
иллюстрирущие реакцию пластинки на растяжение для компонент тензора напряжений 11, 22, 12 на
временных шагах 4 пс, 7 пс, 9 пс, полученные в пакете OVITO, представлены на рис. 1.2.
Рис. 1.1. Геометрия пластинки с центральной трещиной
Fig. 1.1. Geometry of a plate with a central crack
Вестник Самарского университета. Естественнонаучная серия 2023. Том 29, № 4. С. 106–116
Vestnik of Samara University. Natural Science Series 2023, vol. 29, no. 4, pp. 106–116 109
Компоненты тензора напряжений в LAMMPS определяются согласно формуле для тензора вириаль-
ных напряжений, которая была получена путем определения тензона напряжений Коши в терминах
атомистических величин [30]:
av = − 1
V
[
1
2
Σ
 ; 
 ̸= 
f   ⊗ (r  − r ) +
Σ
 
m vrel
 
⊗ vrel
 
]
;
где r ; r  — радиус-векторы частиц  ;  ; vrel
  — вектор скорости частицы  , f   — сила, действующая
на частицу   со стороны частицы  .
Рис. 1.2. Распределение полей напряжений 11 — a) , 22 — b), 12 — c) на временных шагах 4 пс,
7 пс, 9 пс
Fig. 1.2. Distribution of stress fields 11 — a), 22 — b), 12 — c) at time steps of 4 ps, 7 ps, 9 ps
2. Аналитическое представление полей напряжений
в анизотропных материалах с кубической сингонией
В статье [27] приведено асимптотическое решение, учитывающее ортотропную симметрию упругих
свойств материала:


11
22
12

 = 2Re


1Σ
k=1
Ak
i(k+1)2
1 − 2
r
k
2
????1


22

(−1)n+1+1
2
1 (cos + 2sin)
n
2
????1 − 21

(−1)n+1+1
2
2 (cos + 1sin)
n
2
????1

(−1)n+1+1
2
1 (cos + 2sin)
n
2
????1 −
(−1)n+1+1
2
2 (cos + 1sin)
n
2
????1
−
(
2
(−1)n+1+1
2
1 (cos + 2sin)
n
2
????1 − 1
(−1)n+1+1
2
2 (cos + 1sin)
n
2
????1
)




+
+2Re


1Σ
k=1
Bk
i(k+1)2
1 − 2
r
k
2
????1


22

(−1)n+1
2
1 (cos + 2sin)
n
2
????1 − 21

(−1)n+1
2
2 (cos + 1sin)
n
2
????1

(−1)n+1
2
1 (cos + 2sin)
n
2
????1 −
(−1)n+1
2
2 (cos + 1sin)
n
2
????1
−
(
2
(−1)n+1
2
1 (cos + 2sin)
n
2
????1 − 1
(−1)n+1
2
2 (cos + 1sin)
n
2
????1
)




; (2.1)
где Ai;Bi — амплитудные множители поля напряжений, 1; 2 — константы, описывающие материал
(корни характеристического уравнения), r; — полярные координаты с полюсом в вершине трещины.
Амплитудные коэфициенты для бесконечной анизотропной пластины с центральной трещиной при
одноосном нагружении:
A1 =
√
2a
4
sin2 ; B1 =
√
2a
4
sin  cos ;
110
Мушанкова К.А., Степанова Л.В. Опыт моделирования наклонных трещин в материалах с кубической ...
Mushankova K.A., Stepanova L.V. Experience in modeling inclined cracks in materials with cubic crystal ...
A2 =

2Im(1 + 2)
[
cos2  + Re(12)sin2  + Re(1 + 2)
]
; B2 =
Re(1 + 2)
2Im(12)
[sin  cos  − ] ;
A3 =
3
8
√
2a
sin2 ; B3 =
3
8
√
2a
sin  cos ;
A2n+2 = B2n+2 = 0; n = 1; 2; 3; :::
A2n+3 =
(−1)n+1sin2 
8(2a)n+1
2
[
−4 × 1 × 3 × 5 × :::(2n − 1)
2 × 4 × 6 × :::2n
+
3 × 5 × 7 × :::(2n + 1)
4 × 6 × 8 × :::(2n + 2)
]
; n = 1; 2; 3; :::
B2n+3 =
(−1)n+1sin  cos 
8(2a)n+1
2
[
−4 × 1 × 3 × 5 × :::(2n − 1)
2 × 4 × 6 × :::2n
+
3 × 5 × 7 × :::(2n + 1)
4 × 6 × 8 × :::(2n + 2)
]
; n = 1; 2; 3; :::
где — действительный параметр, характеризуемый величиной смещения тела как абсолютно твердого,
  — угол наклона трещины, a — половина длины трещины, — прикладываемая к пластине нормальная
нагрузка.
Поскольку кубическая симметрия упругих свойств является частным случаем ортотропной, приведен-
ное решение можно применить к поставленной задаче для определения полей напряжений монокристал-
ла ГЦК-меди. Константы исследуемого материала 1 = 0:708728 + 0:705482i; 2 = −0:708728 + 0:705482i
были получены из тензора упругих модулей в статье [23]. Для сравнения асимптотического решения и
результатов МД-моделирования для компонент тезора напряжения 11, 22 и 12, полученных анали-
тически, был выполнен поворот на 45 градусов.
3. Сравнение аналитического решения и результата
МД-моделирования
На рисунках 3.1–3.3 слева изображены контурные графики, построенные по асимпторическому ре-
шению, обобщающему решение Уильямса, при удержании 20 слагаемых, справа представлены резуль-
таты моделирования, выполненные в программном пакете LAMMPS, где визуализируются вириальные
напряжения. Справа приведены кадры из OVITO, показывающие распределение компонент тензора на-
пряжений рис. 3.1 –– 11, рис. 3.2–22, рис. 3.3 — 12. На графиках, полученных методом молекулярной
динамики, четко видно как расположена трещина относительно осей симметрии рассматриемого матери-
ала с ГЦК-решеткой. На всех графиках напряжения имеют размерность ГПа, а трещина располагается
в 3 четверти под углом 45 градусов к главным осям материала.
Рис. 3.1. Контурный график, построенный согласно аналитическому решению (слева),
и распределение полей напряжений, полученное путем МД-моделирования (справа), для компоненты
тензона напряжений 11
Fig. 3.1. Contour plot plotted according to the analytical solution (left) and the distribution of stress fields
obtained by MD modeling (right) for the stress tensor component 11
На рис. 3.1 как слева, так и справа можно видеть характерные области минимальных и максималь-
ных напряжений, представленных различными оттенками в соответствии со шкалой. Область макси-
мальных напряжений как на контурном графике, так и на визулизации МД-решения расположена под
трещиной, область минимальных напряжений на обоих графиках образуют характерные петли, прихо-
дящиеся на центр между 1 и 2 четвертями окружностей.
Вестник Самарского университета. Естественнонаучная серия 2023. Том 29, № 4. С. 106–116
Vestnik of Samara University. Natural Science Series 2023, vol. 29, no. 4, pp. 106–116 111
На рис. 3.2 приведены визуализации нормальной компоненты тензора напряжений. Как на левом,
так и на правом графиках концентрация напряжений приходится на вершину дефекта, а берега трещи-
ны обладают сравнительно меньшими значениями напряжений. Максимальные напряжения формируют
некие петлеобразные области, распространяющиеся вправо от вершины трещины, и наблюдаются как
для аналитического решения, так и для численного.
Рис. 3.2. Контурный график, построенный согласно аналитическому решению (слева),
и распределение полей напряжений, полученное путем МД-моделирования (справа), для компоненты
тензона напряжений 22
Fig. 3.2. Contour plot plotted according to the analytical solution (left) and the distribution of stress fields
obtained by MD modeling (right) for the stress tensor component 22
Сравнение результатов асимптотического решения и молекулярно-динамического для касательной
компоненты тензора напряжений 12 представлено на рис. 3.3. Так же, как и для компонент тензора
напряжений 11, 22, между результатами, полученными различными методами, для 12 наблюдается
схожий характер распределения напряжений. Максимальные напряжения приходятся на вершину тре-
щины, а затем петлеобразно распространяются вправо на уменьшение. Минимальные напряжения тоже
имеют петлееобразный характер и распространяются влево с угасанием по модулю.
Рис. 3.3. Контурный график, построенный согласно аналитическому решению (слева),
и распределение полей напряжений, полученное путем МД-моделирования (справа), для компоненты
тензона напряжений 12
Fig. 3.3. Contour plot plotted according to the analytical solution (left) and the distribution of stress fields
obtained by MD modeling (right) for the stress tensor component 12
Для молекулярно-динамического решения выделялась кольцевая область в окрестности вершины тре-
шины, и далее было проведено сравнение угловых распределений для компонент тезора напряжений 11,
22 и 12 на расстоянии r=a = 0:75 (рис. 3.4).
112
Мушанкова К.А., Степанова Л.В. Опыт моделирования наклонных трещин в материалах с кубической ...
Mushankova K.A., Stepanova L.V. Experience in modeling inclined cracks in materials with cubic crystal ...
Рис. 3.4. Угловые распределения для компонент тензора напряжений 11, 22 и 12, где линия —
результат аналитического решения, точки — данные, полученные МД-методом
Fig. 3.4. Angular distributions for stress tensor components 11, 22 и 12, where the line is the result of the
analytical solution, the points are data obtained by the MD method
Выводы
В данной статье выявлено хорошее совпадение атомистических и континуальных полей напряжений у вершины трещины в условиях смешанного нагружения в анизотропной среде с кубической симметрией упругих свойств. Атомистические распределения напряжений, ассоциированные с вершиной трещины, получены с помощью метода молекулярной динамики. Континуальные распределения получены из теоретического решения задачи определения напряженно-деформированного состояния у вершины трещины, базирующегося на методах теории упругости анизотропных сред и последующего разложения комплексных потенциалов по собственным функциям. В рамках молекулярно-динамического вычислительного эксперимента рассмотрена монокристаллическая гранецентрированная медь при низкой температуре, с тем чтобы выделить упругий режим деформирования монокристалла, и был использован потенциал внедренного атома. Отличительной особенностью проведенного молекулярно-динамического моделирования является рассмотрение трещины, которая составляет различные углы с плоскостями симметриии кристалла. В окрестности вершины трещины выбирались точки, лежащие в кольцевых областях, на различном расстоянии от вершины трещины и различной толщины, и строились зависимости компонент тензора напряжений в зависимости от полярного угла. Сравнение угловых зависимостей, полученных посредством атомистического расчета и с помощью теоретического решения, показало их хорошую согласованность. Обнаружено, что сходство угловых зависимостей компонент тензора напряжений наблюдается при всех изученных значениях двух углов: угла между осью симметрии кристаллической решетки (в плоскости пластины) и направлением трещины и угла между направлением действия растягивающей нагрузки и линией трещины.
В силу указанного свойства решений можно заключить, что решения континуальной механики разрушения могут служить для описания полей напряжений на атомистических расстояниях от вершины дефекта.

Об авторах

К. А. Мушанкова

Самарский национальный исследовательский университет имени академика С.П. Королева

Автор, ответственный за переписку.
Email: karinamushankova@yandex.ru
ORCID iD: 0000-0002-6081-1169

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

Россия, 443086, Российская Федерация, г. Самара, Московское шоссе, 34

Л. В. Степанова

Самарский национальный исследовательский университет имени академика С.П. Королева

Email: stepanovalv2015@yandex.ru
ORCID iD: 0000-0002-6693-3132

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

Россия, 443086, Российская Федерация, г. Самара, Московское шоссе, 34

Список литературы

Дополнительные файлы