COMPUTER SIMULATION OF CRACK GROWTH. MOLECULAR DYNAMICS METHOD


Cite item

Abstract

The aim of the study is to determine the stress intensity factors using molecular dynamics (MD) method. In the course of the study, a computer simulation of the propagation of a central crack in a copper plate was carried out. The simulation was performed in the LAMMPS (Large-scale Atomic / Molecular Massively Parallel Simulator) software package. A comprehensive study of the influence of geometric characteristics (model dimensions, crack length), temperature, strain rate and loading mixing parameter on the plate strength, crack growth and direction was carried out. The article proposes a method for determining the coefficients of the asymptotic expansion of M. Williams’ stress fields. The analysis of the influence of the choice of points on the calculation of the coefficients and the comparison of the results obtained with the analytical solution are carried out.

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

     

    image

    Рис. 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], который вычисляется по следующей формуле:

    image

    image

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

    π r0 σ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 103 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

     

    image

    a б

    image

    в г

    Рис. 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

     

    image

    a б

    image

    в г

    Рис. 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

     

    image image

    image

    image

    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

     

     

    image image

    image

    image

    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 пс.

 

image image image

 

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

 

image image image

 

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

 

image image image

 

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

 

image image image

 

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

 

image image image

 

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

 

 

image image image

 

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

 

 

image image image

 

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

 

 

image image image

 

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

 

 

image image image

 

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

 

image image image

 

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

  1. Thaulow С. Application of CTOD in atomistic modeling of fracture. Procedia Materials Science, 2014, vol. 3, pp. 1542–1547. DOI: http://dx.doi.org/10.1016/j.mspro.2014.06.249.
  2. Horstemeyer M.F., Lim J., Lu W.Y., Mosher D.A., Baskes M.I., Prantil V.C., Plimpton S.J. Torsion/Simple Shear of Single Crystal Copper. ASME. Journal of Engineering Materials and Technology, 2002, vol. 124 (3), pp. 322—328. DOI: http://dx.doi.org/10.1115/1.1480407.
  3. Malikova L. Multi-parameter fracture criteria for the estimation of crack propagation direction applied to a mixed-mode geometry. Engineering Fracture Mechanics, 2015, vol. 143, pp. 32–46. DOI: http://dx.doi.org/10.1016/j.engfracmech.2015.06.029.
  4. Malikova L., Vesely V., Seitl S. Crack propagation direction in a mixed mode geometry estimated via multi-parameter fracture criteria. International Journal of Fatigue, 2016, vol. 89, pp. 99–107. DOI: http://dx.doi.org/10.1016/j.ijfatigue.2016.01.010. [5] Malikova L., Vesely V. Estimation of the crack propagation direction in a mixed-mode geometry via multi-parameter fracture criteria. Frattura ed Integrita Strutturale, 2015, vol. 33, pp. 25–32. DOI: http://doi.org/10.3221/IGF-ESIS.33.04.
  5. Malikova L., Vesely V. Influence of the elastic mismatch on crack propagation in a silicate-based composite. Theoretical and Applied Fracture Mechanics, 2017, vol. 91, pp. 25–30. DOI: https://doi.org/10.1016/J.TAFMEC.2017.03.004.
  6. Patil P., Vyasarayani C.P., Ramji M. Linear least squares approach for evaluating crack tip fracture parameters using isochromatic and isoclinic data from digital photoelasticity. Optics and Lasers in Engineering, 2017, vol. 93, pp. 182–194. DOI: http://dx.doi.org/10.1016/j.optlaseng.2017.02.003.
  7. Williams M.L. Stress Singularities Resulting From Various Boundary Conditions in Angular Corners of Plates in Extension. Journal of Applied Mechanics, 1952, vol. 74, pp. 526–528. Available at: http://authors.library.caltech.edu/47672/1/382785.pdf.
  8. Williams M.L. On the Stress Distribution at the Base of a Stationary Crack. Trans. ASME. Journal of Applied Mechanics, 1957, vol. 24, pp. 109–114. Available at: https://core.ac.uk/reader/33111630.
  9. Hello G., Tahar M.B., Roelandt J.-M. Analytical determination of coefficients in crack-tip stress expansions for a finite crack in an infinite plane medium. International Journal of Solids and Structures, 2012, vol. 49, pp. 556–566. DOI: http://dx.doi.org/10.1016/j.ijsolstr.2011.10.024.
  10. Stepanova L.V. Asymptotic Analysis of the Crack Tip Stress Field (Consideration of Higher Order Terms). Numerical Analysis and Applications, 2019, vol. 22, no. 3, pp. 284–296. DOI: https://doi.org/10.15372/SJNM20190307. (English; Russian original)
  11. Sobek J., Frantik P., Vesely V. Analysis of accuracy of Williams series approximation of stress field in cracked body - influence of area of interest around crack-tip on multi-parameter regression performance. Frattura ed Integrita Struсturale, 2017, vol. 11, no. 39, pp. 129–142. DOI: http://dx.doi.org/10.3221/IGF-ESIS.39.14.
  12. Stepanova L.V. Experimental determination and finite element analysis of coefficients of the multi-parameter williams series expansion in the vicinity of the crack tip in linear elastic materials. Part i. PNRPU Mechanics Bulletin, 2020, vol. 4, pp. 237–249. DOI: https://doi.org/10.15593/perm.mech/2020.4.20.
  13. Belova O.N., Stepanova L.V. Determination of the coefficients of asymptotic crack — tip stress expansion. Mixed mode loading of the plate. Vestnik Samarskogo universiteta. Estestvennonauchnaia seriia = Vestnik of Samara University. Natural Science Series, 2020, vol. 26, no. 3, pp. 40-–62. Available at: https://journals.ssau.ru/est/article/view/8689.
  14. Wilson M.A., Grutzik S.J., Chandross M. Continuum stress intensity factors from atomistic fracture simulations. Computer Methods in Applied Mechanics and Engineering, 2019, vol. 354, pp. 732—749. DOI: http://dx.doi.org/10.1016/j.cma.2019.05.050.
  15. Ersland C.H., Thaulow C., Vatne I.R., Ostby E. Atomistic modeling of micromechanisms and T-stress effects in fracture of iron. Engineering Fracture Mechanics, 2012, vol. 79, pp. 180-–190. DOI: http://dx.doi.org/10.1016/j.engfracmech.2011.10.012.
  16. Lee G.H., Beom H.G. Mixed-mode fracture toughness testing of a Cu/Ag bimetallic interface via atomistic simulations. Computational Materials Science, 2020, vol. 183, p. 109806. DOI: http://dx.doi.org/10.1016/j.commatsci.2020.109806.
  17. Website of LAMMPS Molecular Dynamics Simulator. Available at: https://lammps.sandia.gov.
  18. Website of NIST Interatomic Potentials Repository Project. Available at: https://www.ctcms.nist.gov/potentials.
  19. Shih C.F. Elastic-plastic analysis of combined mode crack problems: PhD Thesis. Harvard University, 1973.
  20. Website of Supercomputing Center of Samara University. Available at: http://hpc.ssau.ru. (In Russ.)
  21. Website of OVITO – Open Visualization Tool. Available at: https://www.ovito.org.

Copyright (c) 2021 Belova O.N., Stepanova L.V., Chapliy D.V.

Creative Commons License
This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.

This website uses cookies

You consent to our cookies if you continue to use our website.

About Cookies