COMPUTER SIMULATION OF CRACK GROWTH. MOLECULAR DYNAMICS METHOD

Full Text

Введение

Большинство исследований разрушения основано на энергии Гриффитса, рассчитанной по поверхностной энергии, и коэффициенте интенсивности напряжений, полученном из механических формулировок континуума. Оба подхода, однако, отражают механику сплошных сред и не имеют прямой связи с атомистическими процессами вблизи кончика трещины. В то же время они могут быть значимыми при разрушении материала [1]. Множество исследований посвящено изучению процессов на различных масштабных уровнях и сравнению полученных результатов. В работе [2] анализируют простой сдвиг и кручение монокристаллической меди, используя эксперименты, моделирование молекулярной динамикой и конечных элементов с тем, чтобы сравнить процессы на разных масштабах. Основным сходством, наблюдаемым во всех моделированиях и экспериментах, была колебательная картина деформации, наблюдавшаяся на внешней поверхности образца. Колебательная картина возникла потому, что кинематика качественно одинакова в молекулярной динамике (МД), конечно-элементном моделировании и экспериментах. Коэффициенты интенсивности напряжений (КИН) используются в механике сплошного разрушения для количественной оценки полей напряжений, окружающих трещину в однородном материале в линейно-упругом режиме. Критические значения КИН определяют внутреннюю меру сопротивления материала распространению трещины. Вычислениям коэффициентов интенсивности напряжений посвящено много исследований, в том числе с помощью метода молекулярной динамики. Способность характеризовать разрушения на атомном масштабе дает возможность исследовать, как микроструктура материала влияет на точность определения КИН. Вычислению коэффициентов интенсивности напряжений и высших приближений разложения Вильямса посвящен ряд работ [3–14] В работе [15] предложена методика вычисления КИН с помощью МД. Точность этого подхода сначала исследуется для простой модели, а затем применяется к атомистическому моделированию разрушения в аморфном кремнеземе. Показана зависимость вычисленных КИН от времени и от выбора точек. Работа [16] посвящена исследованию сопряжения атомистического моделирования разрушения в железе с континуальным решением. Предметом исследования является многомасштабный анализ трещины. Рассматривается влияние ограничений (геометрия, размер трещины и режим нагружения) на механизмы разрушения путем изменения T-напряжения. Авторы пришли к выводу, что атомистическое моделирование будет очень полезным инструментом для понимания перехода от хрупкого состояния к пластичному в сталях. В статье [17] представлена оценка вязкости разрушения в смешанном режиме для межфазной трещины с использованием МД.

 

1. Моделирование пластины

   В статье с помощью метода молекулярной динамики моделируется процесс распространения трещины в монокристаллической пластине с трещиной под действием смешанного нагружения. Моделирование проведено в программном комплексе LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) [18]. Моделирование проходит в несколько этапов. На первом этапе описывается ячейка моделирования, в нашем случае ячейка типа BOX, и заполняется атомами с определенной кристаллической решеткой. Пластина моделировалась квадратной формы. Построен ряд моделей с различными размерами

   модели и трещины. Длина стороны пластины W выбиралась от 50 до 200 элементарных ячеек, что соответствует ширине от 180.75 ˚A до 723 ˚A. Толщина пластины составляла от 4 до 10

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

   В работах по моделированию трещин с помощью метода молекулярной динамики встречаются два способа их задания. В основном трещина создается путем удаления нескольких слоев кристаллических ячеек.

   Другой способ — это задание двух областей, одна под трещиной и вторая над ней. При этом связи между атомами этих областей устанавливали нулевыми. Таким образом, формируется надрез.

   В данной работе трещина создавалась путем удаления одного ряда кристаллических ячеек, при этом между областями под трещиной и над ней удалялись связи между атомами. Длина трещины варьировалась от 5 до 10 % от ширины пластины.

   Схематически описание модели приведено на рис. 1.

   На втором этапе проводится описание материала и задание потенциала взаимодействия. Наиболее полный репозиторий межатомных потенциалов создан на сайте [19]. Пользователям предлагается свободно загружать и использовать проверенные исследователями межатомные потенциалы. В ходе исследования моделировалась медная и алюминиевая пластины. Поведение атомов в пластине задавалось с помощью потенциала погруженного атома (англ. Embedded Atom Model (EAM)). Данный потенциал активно используется для моделирования поведения меди и алюминия [2]. Для моделирования

   Белова О.Н., Степанова Л.В., Чаплий Д.В. Компьютерное моделирование роста трещины...

   46Belova O.N., Stepanova L.V., Chapliy D.V. Computer simulation of crack growth. Molecular dynamics method

    

   

   Рис. 1. Описание модели

   Fig. 1. Description of the model

    

   меди в работе используются потенциалы Cu_u3.eam, Cu_u6.eam, для моделирования алюминия — Al_u3.eam.

   Далее производится релаксация системы, приведение системы к нужной температуре. Температура моделирования устанавливалась в разных экспериментах от 10 до 300 K.

   Приложение нагрузки моделируется путем задания растягивающей и сдвиговой деформации ячейки моделирования. Скорость деформации P варьируется от 0.001 ˚A/пс до 0.01 ˚A/пс. Соотношение растягивающей и сдвиговой деформаций характеризуется параметром смешанности нагружения — M p, впервые введенным в работе Ши [20], который вычисляется по следующей формуле:

   

   

   M p = 2 arctan lim σθθ (r, θ = 0) . (1)

   π r→0 σr,θ (r, θ = 0)

   Сдвигу соответствует нулевое значение параметра смешанности, чистому растяжению соответствует единица, значение 0 < M p < 1 принимается для смешанного нагружения образца с трещиной. В данном исследовании моделируется чистое растяжение пластины.

   На последнем этапе осуществляются расчет модели, вывод результатов и их анализ. Расчет растяжения пластины с центральной трещиной производился на суперкомпьютере «Сергей Королев» [21]. Визуализация результатов проводится с помощью программы Open Visualization Tool (OVITO) [22]. Для вычисления коэффициентов интенсивности напряжений была использована программа, написанная авторами в пакете символьной математики Maple. Подробнее с реализацией этой программы можно познакомиться в работах [13; 14].

    

2. Вычисление коэффициентов интенсивности напряжений

   Создана пластина со сторонами W = 200, H = 200 элементарных ячеек, что соответствует размерам 723 ˚A x 723 ˚A. Толщина пластины равна четырем кристаллическим элементарным ячейкам, что соответствует 14.46 ˚A. Длина трещины составляла 10 % от ширины пластины. Во всех трех направлениях были заданы периодические граничные условия. Всего в моделировании участвовало порядка 720 тыс. атомов. После создания элементарной ячейки моделирования система приводилась в равновесное состояние в течение 2 пс (пикосекунд) шагов с ансамблем nve к температуре, равной 10 K. Один шаг по времени равен tstep = 0.001 пс. Приложение нагрузки осуществили путем деформации

   ячейки моделирования. Изменение размера ячейки осуществляется с постоянной скоростью инженерной деформации, вдоль вертикальной оси она равна erate = 2 ∗ 10−3 1/пс. Высота ячейки моделирования,

   а соответственно и пластины L в зависимости от времени будет меняться следующим образом: L(t) =

   = L0(1 + erate ∗ dt), где dt — прошедшее время (в единицах времени), L0 — исходная длина коробки.

   Нагружение пластины моделировалось в течение 31 пс.

   На рисунках далее представлена молекулярная модель трещины в медной пластине в различные моменты времени. На рис. 2 показано распределение компоненты тензора напряжений σ11 через каждые 5 пс.

   На рис. 3 приведено распределение компоненты тензора напряжений σ22 через каждые 5 пс. На рис. 4 приведено распределение компоненты тензора напряжений σ12 через каждые 5 пс.

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

   Вестник Самарского университета. Естественнонаучная серия. 2020. Том 26, № 4. С. 44–55

   Vestnik of Samara University. Natural Science Series. 2020, vol. 26, no. 4, pp. 44–55 47

    

   

   a б

   

   в г

   Рис. 2. Угловое распределение компоненты тензора напряжений σ11 в зависимости от времени:

   а — t = 5 пс, б — t = 10 пс, в — t = 15 пс, г — t = 20 пс

   Fig. 2. Angular distribution of the stress tensor component σ11 depending on time:

   a — t = 5 ps, b — t = 10 ps, c — t = 15 ps, d — t = 20 ps

    

   

   a б

   

   в г

   Рис. 3. Угловое распределение компоненты тензора напряжений σ22 в зависимости от времени:

   а — t = 5 пс, б — t = 10 пс, в — t = 15 пс, г — t = 20 пс

   Fig. 3. Angular distribution of the stress tensor component σ22 as a function of time:

   a — t = 5 ps, b — t = 10 ps, c — t = 15 ps, d — t = 20 ps

   Белова О.Н., Степанова Л.В., Чаплий Д.В. Компьютерное моделирование роста трещины...

   48Belova O.N., Stepanova L.V., Chapliy D.V. Computer simulation of crack growth. Molecular dynamics method

    

   

   

   

   a б

    

   в г

   Рис. 4. Угловое распределение компоненты тензора напряжений σ12 в зависимости от времени:

   а — t = 5 пс, б — t = 10 пс, в — t = 15 пс, г — t = 20 пс

   Fig. 4. Angular distribution of the stress tensor component σ12 depending on time:

   a — t = 5 ps, b — t = 10 ps, c — t = 15 ps, d — t = 20 ps

    

    

   

   

   

   a б

    

   в г

   Рис. 5. Способы выбора областей: а — кольцо радиусом 36 ˚A, б — кольцо радиусом от 15 до 36 ˚A,

   в — кольцо радиусом от 25 до 36 ˚A, г — кольцо радиусом от 34 до 36 ˚A

   Fig. 5. Methods for selecting regions: a — a ring with a radius of 36 ˚A, b — a ring with a radius of 15 to

   36 ˚A, c — ring with radius from 25 to 36 ˚A, g — ring with radius from 34 to 36 ˚A

    

3. Вычисление коэффициентов

 

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

Вестник Самарского университета. Естественнонаучная серия. 2020. Том 26, № 4. С. 44–55

Vestnik of Samara University. Natural Science Series. 2020, vol. 26, no. 4, pp. 44–55 49

 

распределение компоненты тензора напряжений σ22 в момент времени 5 пс. Множество атомов выбиралось в виде круговой или кольцевой области.

На рис. 6–11 приведены графики сравнения распределения компонент напряжений σ11, σ22, σ12, полученные с помощью аналитического решения (пунктирная линия) и с помощью восстановленного решения по коэффициентам разложения Вильямса, вычисленным с использованием метода молекулярной динамики (сплошная линия) в зависимости от выбора множества атомов в момент времени 5 пс.

 

 

a б в

Рис. 6. Графики сравнения решений, для атомов, выбранных из кольцевой области радиусом

от 10 до 15 ˚A: а — σ11, б — σ22, в — σ12

Fig. 6. Graphs of comparison of solutions, for atoms selected from an annular region with a radius of 10 to 15 ˚A:

a — σ11, b — σ22, c — σ12

 

 

a б в

Рис. 7. Графики сравнения решений для атомов, выбранных из кольцевой области радиусом

от 10 до 20 ˚A: а — σ11, б — σ22, в — σ12

Fig. 7. Graphs of comparison of solutions, for atoms selected from an annular region with a radius of 10 to 20 ˚A:

a — σ11, b — σ22, c — σ12

 

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

Белова О.Н., Степанова Л.В., Чаплий Д.В. Компьютерное моделирование роста трещины...

50Belova O.N., Stepanova L.V., Chapliy D.V. Computer simulation of crack growth. Molecular dynamics method

 

 

a б в

Рис. 8. Графики сравнения решений для атомов, выбранных из кольцевой области радиусом

от 15 до 25 ˚A: а — σ11, б — σ22, в — σ12

Fig. 8. Graphs of comparison of solutions, for atoms selected from an annular region with a radius of 15 to 25 ˚A: a — σ11, b — σ22, c — σ12

 

 

a б в

Рис. 9. Графики сравнения решений для атомов, выбранных из кольцевой области радиусом

от 15 до 36 ˚A: а — σ11, б — σ22, в — σ12

Fig. 9. Graphs of comparison of solutions, for atoms selected from the annular region with a radius of 15 to 36 ˚A: а — σ11, b — σ22, c — σ12

 

 

a б в

Рис. 10. Графики сравнения решений для атомов, выбранных из кольцевой области радиусом

от 25 до 36 ˚A: а — σ11, б — σ22, в — σ12

Fig. 10. Graphs of comparison of solutions, for atoms selected from an annular region with a radius of 25 up to 36 ˚A: a — σ11, b — σ22, c — σ12

Вестник Самарского университета. Естественнонаучная серия. 2020. Том 26, № 4. С. 44–55

Vestnik of Samara University. Natural Science Series. 2020, vol. 26, no. 4, pp. 44–55 51

 

 

 

a б в

Рис. 11. Графики сравнения решений для атомов, выбранных из кольцевой области радиусом

от 34 до 36 ˚A: а — σ11, б — σ22, в — σ12

Fig. 11. Graphs of comparison of solutions, for atoms selected from the annular region with a radius of 34 up to 36 ˚A: a — σ11, b — σ22, c — σ12

 

 

 

a б в

Рис. 12. Графики сравнения решений для атомов, выбранных из круговой области радиусом 36 ˚A:

а — σ11, б — σ22, в — σ12

Fig. 12. Graphs of comparison of solutions for atoms selected from a circular region with a radius of 36 ˚A:

a — σ11, b — σ22, c — σ12

 

 

 

a б в

Рис. 13. Графики сравнения решений в момент времени 5 пс: а — σ11, б — σ22, в — σ12

Fig. 13. Graphs of comparison of solutions at the moment of time 5 ps: a — σ11, b — σ22, c — σ12

Белова О.Н., Степанова Л.В., Чаплий Д.В. Компьютерное моделирование роста трещины...

52Belova O.N., Stepanova L.V., Chapliy D.V. Computer simulation of crack growth. Molecular dynamics method

 

 

 

a б в

Рис. 14. Графики сравнения решений в момент времени 10 пс: а — σ11, б — σ22, в — σ12

Fig. 14. Graphs of comparison of solutions at the moment of time 10 ps: a — σ11, b — σ22, c — σ12

 

 

a б в

Рис. 15. Графики сравнения решений в момент времени 15 пс: а — σ11, б — σ22, в — σ12

Fig. 15. Graphs of comparison of solutions at time 15 ps: a — σ11, b — σ22, c — σ12

 

На рис. 13, 14, 15 приведены аналогичные графики сравнения с течением времени для круговой области выбора атомов.

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

 

Результаты вычислений и выводы

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

About the authors

O. N. Belova

Samara National Research University

Author for correspondence.
Email: BelovaONik@yandex.ru
ORCID iD: 0000-0002-4492-223X

postgraduate student of the Department of Mathematical Modelling in Mechanics

34, Moskovskoye shosse, 443086, Russian Federation.

L. V. Stepanova

Samara National Research University

Email: stepanovalv@samsu.ru
ORCID iD: 0000-0002-6693-3132

Doctor of Physical and Mathematical Sciences, professor, Department of Mathematical Modelling in Mechanics

Russian Federation

D. V. Chapliy

Samara National Research University

Email: DCH300189@yandex.ru
ORCID iD: 0000-0001-9510-3659

postgraduate student of the Department of Mathematical Modeling in Mechanics

Russian Federation

References

Supplementary files