Experience in modeling inclined cracks in materials with cubic crystal structure

Full Text

Введение

Появление заметной трещины на макроскопическом уровне является следствием распространения трещины в нескольких масштабах, сильно отличающихся по длине. Наименьший масштаб — это разрыв связей между атомами по мере распространения трещины. Движущей силой, стоящей за разрывом соединения, является поле напряжений, обусловленное нагрузкой, приложенной к образцу. Поскольку материал разрушается, когда концентрация напряжений превышает предел прочности материалов на разрушение, анализ разрушения напрямую связан с концентрацией напряжений или деформаций вокруг вершины трещины. Механическая реакция материалов, подверженных экстремальному приложенному напряжению, контролируется атомистическими механизмами вблизи концентраторов напряжений, таких как вершины трещин. В настоящее время атомистическое моделирование оказалось полезным инструментом для изучения свойств материалов, при разрушении путем исследования атомистической конфигурации и поля напряжений вокруг вершины трещины с атомистической точки зрения. В последнее время пристальное внимание уделяется атомистическим напряжениям, обусловленным разрушением кристаллографической решетки материала вдоль границ зерен. Проведенные исследования показали, что усилия, возникающие вдоль границ зерен, вычисленные с помощью компонент тензора атомистических напряжений, могут служить количественными дескрипторами границ зерен [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.

Для восполнения имеющегося пробела в настоящей статье изучена пластина из монокристаллической меди, ослабленная наклонной трещиной. Трещина наклонена под различными углами по отношению к главным осям симметрии кристаллической решетки материала с кубической симметрией его упругих свойств. Вычисления выполнены с помощью метода молекулярной динамики, реализованного в программе Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS). Ведущей целью вычислений было сравнение угловых распределений компонент тензора напряжений, ассоциированных с вершиной трещины, полученных методом молекулярной динамики, и полученных с помощью аналитического решения задачи, построенного с использованием методов теории упругости анизотропного тела [27–29].

1. Моделирование наклонного дефекта в материалах с кубической кристаллической решеткой

Для исследования полей напряжений в окрестности вершины трещины на наноразмерном уровне был написан скрипт для LAMMPS. В рамках поставленной задачи построения пластинки с центральным дефектом на наноскопическом уровне удобно оперировать с размерностями Å для длины, пс для времени и бар для напряжения, поэтому применялась команда units metal. Несмотря на то что исследумый объект — пластинка, все же она будет обладать некоторой толщиной в несколько атомов, поэтому в коде прописывалась размерность пространства dimension 3. Задавались периодические граничные условия во всех направлениях boundary p p p. Определение стиля представления частиц atom style atomic означает, что будет рассматриваться простейший случай, когда каждая частица в системе считается точечной без внутренней структуры или ориентации. Задание решетки (в данном случае гранецентрированная кубическая решетка с параметром решетки 3.615 Å) определяется командой lattice fcc 3.615. Создание бокса для размещения атомов размерами 1201203 в единицах параметра решетки 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 формировался файл, позволяющий получить данные временной эволюции полей напряжений, воспринимаемые пакетами для визуализации. Кадры, иллюстрирущие реакцию пластинки на растяжение для компонент тензора напряжений ${\sigma }_{11}$, ${\sigma }_{22}$, ${\sigma }_{12}$на временных шагах 4 пс, 7 пс, 9 пс, полученные в пакете OVITO, представлены на рис. 1.2.

 

Рис. 1.1. Геометрия пластинки с центральной трещиной

Fig. 1.1. Geometry of a plate with a central crack

 

Рис. 1.2. Распределение полей напряжений ${\sigma }_{11}$— a) , ${\sigma }_{22}$— b), ${\sigma }_{12}$— c) на временных шагах 4 пс, 7 пс, 9 пс

Fig. 1.2. Distribution of stress fields ${\sigma }_{11}$— a), ${\sigma }_{22}$— b), ${\sigma }_{12}$— c) at time steps of 4 ps, 7 ps, 9 ps

 

Компоненты тензора напряжений в LAMMPS определяются согласно формуле для тензора вириальных напряжений, которая была получена путем определения тензона напряжений Коши в терминах атомистических величин [30]:

${\sigma }_{\text{av}}=-\frac{1}{V}\left[\frac{1}{2}\underset{\alpha \ne \beta }{\sum _{\alpha \mathrm{,}\beta }}\overline{f_{\mathrm{\alpha \beta }}\otimes \left(r_{\alpha }-r_{\beta }\right)}+\sum _{\alpha }{m}_{\alpha }\overline{v_{\alpha }^{\mathrm{rel}}\otimes v_{\alpha }^{\mathrm{rel}}}\right]\mathrm{,}$

где ${\mathbf{r}}_{\mathbf{\alpha }}\mathbf{,}{\mathbf{r}}_{\mathbf{\beta }}$ — радиус-векторы частиц $\alpha \mathrm{,}\beta$; $v_{\mathbf{\alpha }}^{\mathbf{rel}}$— вектор скорости частицы $\alpha$, ${\mathrm{f}}_{\mathbf{\alpha \beta }}$— сила, действующая на частицу $\alpha$со стороны частицы $\beta$.

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

В статье [27] приведено асимптотическое решение, учитывающее ортотропную симметрию упругих свойств материала:

$\left[\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ {\sigma }_{12}\end{array}\right]=2Re\left(\sum _k^{\infty }{A}_k\frac{i^{{\left(k+1\right)}^2}}{{\mu }_1-{\mu }_2}{r}^{\frac{k}{2}-1}\left[\begin{array}{c}{\mu }_2^2{\mu }_1^{\frac{{\left(-1\right)}^{n+1}+1}{2}}{\left(\cos \theta +{\mu }_2\sin \theta \right)}^{\frac{n}{2}-1}-{\mu }_1^2{\mu }_2^{\frac{{\left(-1\right)}^{n+1}+1}{2}}{\left(\cos \theta +{\mu }_1\sin \theta \right)}^{\frac{n}{2}-1}\\ {\mu }_1^{\frac{{\left(-1\right)}^{n+1}+1}{2}}{\left(\cos \theta +{\mu }_2\sin \theta \right)}^{\frac{n}{2}-1}-{\mu }_2^{\frac{{\left(-1\right)}^{n+1}+1}{2}}{\left(\cos \theta +{\mu }_1\sin \theta \right)}^{\frac{n}{2}-1}\\ -\left({\mu }_2{\mu }_1^{\frac{{\left(-1\right)}^{n+1}+1}{2}}{\left(\cos \theta +{\mu }_2\sin \theta \right)}^{\frac{n}{2}-1}-{\mu }_1{\mu }_2^{\frac{{\left(-1\right)}^{n+1}+1}{2}}{\left(\cos \theta +{\mu }_1\sin \theta \right)}^{\frac{n}{2}-1}\right)\end{array}\right]\right)+$

$+2Re\left(\sum _k^{\infty }{B}_k\frac{i^{{\left(k+1\right)}^2}}{{\mu }_1-{\mu }_2}{r}^{\frac{k}{2}-1}\left[\begin{array}{c}{\mu }_2^2{\mu }_1^{\frac{{\left(-1\right)}^n+1}{2}}{\left(\cos \theta +{\mu }_2\sin \theta \right)}^{\frac{n}{2}-1}-{\mu }_1^2{\mu }_2^{\frac{{\left(-1\right)}^n+1}{2}}{\left(\cos \theta +{\mu }_1\sin \theta \right)}^{\frac{n}{2}-1}\\ {\mu }_1^{\frac{{\left(-1\right)}^n+1}{2}}{\left(\cos \theta +{\mu }_2\sin \theta \right)}^{\frac{n}{2}-1}-{\mu }_2^{\frac{{\left(-1\right)}^n+1}{2}}{\left(\cos \theta +{\mu }_1\sin \theta \right)}^{\frac{n}{2}-1}\\ -\left({\mu }_2{\mu }_1^{\frac{{\left(-1\right)}^n+1}{2}}{\left(\cos \theta +{\mu }_2\sin \theta \right)}^{\frac{n}{2}-1}-{\mu }_1{\mu }_2^{\frac{{\left(-1\right)}^n+1}{2}}{\left(\cos \theta +{\mu }_1\sin \theta \right)}^{\frac{n}{2}-1}\right)\\ \end{array}\right]\right)\mathrm{,}$ (2.1)

где $A_i\mathrm{,}{B}_i$ — амплитудные множители поля напряжений, ${\mu }_1\mathrm{,}{\mu }_2$ — константы, описывающие материал (корни характеристического уравнения), $r\mathrm{,}\theta$ — полярные координаты с полюсом в вершине трещины.

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

$A_1=\frac{\sqrt{2a}}{4}\sigma si{n}^2\alpha \mathrm{,}{B}_1=\frac{\sqrt{2a}}{4}\sigma sin\alpha cos\alpha \mathrm{,}$

$A_2=\frac{\sigma }{2Im\left({\mu }_1+{\mu }_2\right)}\left[co{s}^2\alpha +Re\left({\mu }_1{\mu }_2\right)si{n}^2\alpha +\lambda Re\left({\mu }_1+{\mu }_2\right)\right]\mathrm{,}{B}_2=\frac{\sigma Re\left({\mu }_1+{\mu }_2\right)}{2Im\left({\mu }_1{\mu }_2\right)}\left[sin\alpha cos\alpha -\lambda \right]\mathrm{,}$

$A_3=\frac{3}{8\sqrt{2a}}\sigma si{n}^2\alpha \mathrm{,}{B}_3=\frac{3}{8\sqrt{2a}}\sigma sin\alpha cos\alpha \mathrm{,}$

$A_{2n+2}=B_{2n+2}=0,n=1, 2, 3, ...$

$A_{2n+3}=\frac{{\left(-1\right)}^{n+1}\sigma si{n}^2\alpha }{8{\left(2a\right)}^{n+\frac{1}{2}}}\left[-4\times \frac{1\times 3\times 5\times ...\left(2n-1\right)}{2\times 4\times 6\times \mathrm{...2}n}+\frac{3\times 5\times 7\times ...\left(2n+1\right)}{4\times 6\times 8\times ...\left(2n+2\right)}\right]\mathrm{,}n=1, 2, 3, ...$

$B_{2n+3}=\frac{{\left(-1\right)}^{n+1}\sigma sin\alpha cos\alpha }{8{\left(2a\right)}^{n+\frac{1}{2}}}\left[-4\times \frac{1\times 3\times 5\times ...\left(2n-1\right)}{2\times 4\times 6\times \mathrm{...2}n}+\frac{3\times 5\times 7\times ...\left(2n+1\right)}{4\times 6\times 8\times ...\left(2n+2\right)}\right]\mathrm{,}n=1, 2, 3, ...$

где $\lambda$— действительный параметр, характеризуемый величиной смещения тела как абсолютно твердого, $\alpha$— угол наклона трещины, $a$— половина длины трещины, $\sigma$— прикладываемая к пластине нормальная нагрузка.

Поскольку кубическая симметрия упругих свойств является частным случаем ортотропной, приведенное решение можно применить к поставленной задаче для определения полей напряжений монокристалла ГЦК-меди. Константы исследуемого материала ${\mu }_1=0.708728+0.705482i\mathrm{,} {\mu }_2=-0.708728+0.705482i$ были получены из тензора упругих модулей в статье [23]. Для сравнения асимптотического решения и результатов МД-моделирования для компонент тезора напряжения ${\sigma }_{11}$, ${\sigma }_{22}$и ${\sigma }_{12}$, полученных аналитически, был выполнен поворот на 45 градусов.

3. Сравнение аналитического решения и результата МД-моделирования

На рисунках 3.1–3.3 слева изображены контурные графики, построенные по асимпторическому решению, обобщающему решение Уильямса, при удержании 20 слагаемых, справа представлены результаты моделирования, выполненные в программном пакете LAMMPS, где визуализируются вириальные напряжения. Справа приведены кадры из OVITO, показывающие распределение компонент тензора напряжений рис. 3.1 –– ${\sigma }_{11}$, рис. 3.2 – ${\sigma }_{22}$, рис. 3.3 — ${\sigma }_{12}$. На графиках, полученных методом молекулярной динамики, четко видно как расположена трещина относительно осей симметрии рассматриемого материала с ГЦК-решеткой. На всех графиках напряжения имеют размерность ГПа, а трещина располагается в 3 четверти под углом 45 градусов к главным осям материала.

 

Рис. 3.1. Контурный график, построенный согласно аналитическому решению (слева), и распределение полей напряжений, полученное путем МД-моделирования (справа), для компоненты тензона напряжений ${\sigma }_{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 ${\sigma }_{11}$

 

Рис. 3.2. Контурный график, построенный согласно аналитическому решению (слева), и распределение полей напряжений, полученное путем МД-моделирования (справа), для компоненты тензона напряжений${\sigma }_{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${\sigma }_{22}$

 

Рис. 3.3. Контурный график, построенный согласно аналитическому решению (слева), и распределение полей напряжений, полученное путем МД-моделирования (справа), для компоненты тензона напряжений${\sigma }_{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${\sigma }_{12}$

 

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

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

Сравнение результатов асимптотического решения и молекулярно-динамического для касательной компоненты тензора напряжений ${\sigma }_{12}$ представлено на рис. 3.3. Так же, как и для компонент тензора напряжений ${\sigma }_{11}$, ${\sigma }_{22}$, между результатами, полученными различными методами, для ${\sigma }_{12}$ наблюдается схожий характер распределения напряжений. Максимальные напряжения приходятся на вершину трещины, а затем петлеобразно распространяются вправо на уменьшение. Минимальные напряжения тоже имеют петлееобразный характер и распространяются влево с угасанием по модулю.

Для молекулярно-динамического решения выделялась кольцевая область в окрестности вершины трещины, и далее было проведено сравнение угловых распределений для компонент тезора напряжений ${\sigma }_{11}$, ${\sigma }_{22}$ и ${\sigma }_{12}$ на расстоянии $r/a=0.75$ (рис. 3.4).

 

Рис. 3.4. Угловые распределения для компонент тензора напряжений ${\sigma }_{11}$, ${\sigma }_{22}$и ${\sigma }_{12}$, где линия — результат аналитического решения, точки — данные, полученные МД-методом

Fig. 3.4. Angular distributions for stress tensor components ${\sigma }_{11}$, ${\sigma }_{22}$и ${\sigma }_{12}$, where the line is the result of the analytical solution, the points are data obtained by the MD method

 

Выводы

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

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

About the authors

Karina A. Mushankova

Samara National Research University

Author for correspondence.
Email: karinamushankova@yandex.ru
ORCID iD: 0000-0002-6081-1169

Master Degree student of the Department of Mathematical Modeling in Mechanics

Russian Federation, 34, Moskovskoye shosse, 443086

Larisa V. Stepanova

Samara National Research University

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

Doctor of Physical and Mathematical Sciences, associate professor, head of the Department of Mathematical Modeling in Mechanics

Russian Federation, 34, Moskovskoye shosse, 443086

References

Supplementary files