ОПРЕДЕЛЕНИЕ МЕХАНИЧЕСКИХ СВОЙСТВ МАТЕРИАЛА, МОДЕЛИРУЕМОГО С ПОМОЩЬЮ МЕТОДА МОЛЕКУЛЯРНОЙ ДИНАМИКИ

Полный текст

1. Определение механических свойств материалов, моделируемых с помощью метода молекулярной динамики

   Определение тензора упругих модулей является неотъемлемой задачей любого моделирования, проводимого методом молекулярной динамики. В статье приведена процедура вычисления упругих констант материала, моделируемого с помощью метода молекулярной динамики. Целью исследования является вычисление матрицы упругих постоянных Cij или коэффициентов податливости Sij материалов. Зная эти величины, рассчитывают модули упругости вдоль определенных направлений

   (n) и сдвига в системах сдвига (m, n), а также коэффициент Пуассона. Затем можно рассчитать экстремальные значения величин E(n), G(m, n) и µ(m, n). Кроме того, могут быть определены

    

   

   1Работа выполнена при поддержке РФФИ в рамках гранта № 20-31-90082

   Вестник Самарского университета. Естественнонаучная серия. 2021. Том 27, № 4. С. 14–21

   Vestnik of Samara University. Natural Science Series. 2021, vol. 27, no. 4, pp. 14–21 15

    

   линейные комбинации упругих констант, имеющие физический смысл. В ряде случаев вычисляются технические упругие модули (модули, которые измеряются в опытах на поликристаллах) с помощью усреднения по Фойгту, Ройсу или Хиллу и параметры упругой анизотропии.

   В качестве примеров для моделирования выбраны fcc монокристаллы меди и алюминия. Известно, что закон Гука в тензорном виде имеет следующий вид:

   σij = cijklεkl, (1.1)

   где σij и εkl – тензоры напряжения и деформации. Так как эти тензоры симметричны, то тензор модулей упругости будет обладать симметрией и его можно записать в матричном виде. Закон Гука можно записать в виде

    σ11 

    C1111 C1122 C1133 C1123 C1113 C1112  

   ε11 

    σ22  

   

    

   

    

   

    

    σ33  

   · C2222 C2233 C2223 C2213 C2212  

    

    

   · · C3333 C3323 C3313 C3312  

   ε22 

   

    

   ε33 

    σ23

    = 

     2ε

    . (1.2)

     

   

    

   

    

   

    

    σ13  

   σ12

   · · · C2323 C2313 C2312

   · · · · C1313 C1312

   · · · · · C1212

     23 

     

   2ε

    

     13 

   2ε12

   Закон Гука может быть записан в эквивалентной тензорной форме, через тензор модулей податливости

   Sijkl

    

   Получим обратное матричное уравнение

   εij = sijklσkl. (1.3)

    ε11 

    S1111 S1122 S1133 S1123 S1113 S1112

     σ11 

    ε22  

    ε33  

   · S2222 S2233 S2223 S2213 S2212

     σ22 

     σ 

     

   · · S3333 S3323 S3313 S3312  

   33 

    2ε23

    = 

     σ

    . (1.4)

     

   

    

   

    

   

    

    2ε13  

   2ε12

   · · · S2323 S2313 S2312

   · · · · S1313 S1312

   · · · · · S1212

     23 

     

   σ

    

     13 

   σ12

   В данной статье вычисляют упругие свойства монокристалла меди с определенным потенциалом взаимодействия с помощью метода молекулярной динамики (МД). Механические свойства зависят от кристаллической структуры. Вычислим упругие свойства для гранецентрированной меди. Получить тензор упругих модулей можно с помощью потенциальной энергии (ПЭ), вычисленной в ходе МД расчета. ПЭ зависит от тензора деформаций и может быть разложена в ряд Тейлора:

   6

    

   E(ε) = E(0) + ∑ ∂E ε

   1

    

   ∂

    

   E

    

   6 2

   + ∑ ε ε , (1.5)

   

   

   

   i

    

   ∂ε i

   i

   2

   i,j

   ∂εiεj i j

   где E(0) – энергия начального состояния равновесия, εiεj – деформации в нотации Фойгта.

   Компоненты тензора упругих констант материала можно рассчитать по следующей формуле:

   1 ∂2E

   

   Ci,j = V ∂ε ε

   . (1.6)

   i j

   После минимизации энергии деформирования второй член в разложении 1.5 устраняется.

   Для определения первой компоненты тензора упругих модулей C11 необходимо провести моделирование одноосного растяжения образца. С этой целью зададим тензор деформаций в следующем виде:

    δ 0 0 

   ε =  0 0 0  . (1.7)

   0 0 0

   Выражение для энергии примет следующий вид:

   1 ∂2E 2

   

   1

    

   E(ε1) = E(0) + 2 ∂ε2 δ

   + . . . = α + βδ2

   + γδ3

   + . . . . (1.8)

   Тогда константу C11 определим по формуле

    

   1 ∂2E 2β

   

   C11 = V

   1

    

   ∂ε2

   

   = . (1.9)

   V

   Компоненту C12 не представляется возможным вычислить напрямую, но можно найти разность констант C11 − C12. Зададим компоненты тензора деформаций

    δ 0 0 

   ε =  0 −δ 0  , (1.10)

   0 0 0

   О.Н. Белова Определение механических свойств материала, моделируемого с помощью метода молекулярной...

   16 Belova O.N. Determination of mechanical properties of a material modeled using the molecular dynamics method

    

   тогда уравнение разложения потенциальной энергии примет вид

   1 ∂2E 2

   1 ∂2E 2

   1 ∂2E 2 2 3

   

   E(ε1, ε2) = E(0) + 2 ∂ε2 δ

   

   + δ

   2 ∂ε2

   

   + δ

   2 ∂ε ∂ε

   + . . . = α + βδ

   + γδ

   + . . . (1.11)

   1 2 1 2

   Константу C12 вычисляем из следующего выражения:

   β

   

   C12 = C11 − V . (1.12)

   Для вычисления компоненты тензора упругих модулей C44 зададим тензор деформаций в следующем виде:

    0 0 0

   ε =  0 0 η

   0 η 0

   

    , (1.13)

   где δ = 2η. Выражение для энергии примет вид

   1 ∂2E 2

   

   4

    

   E(ε4) = E(0) + 2 ∂ε2 δ

   ∂2E

   

   + . . . = E(0) + 2 η2

   4

    

   ∂ε2

   + . . . = α + βη2

   + γη3

   + . . . . (1.14)

   Тогда C44 найдем по формуле

   1 ∂2E β

   

   C44 = V

   4

    

   ∂ε2

   

   = (1.15)

   2V

    

2. Определение механических свойств монокристаллов меди и алюминия

МД моделирование выполнено с применением потенциала внедренного атома (embedded atom model — EAM), который является наиболее часто используемым потенциалом для моделирования меди. Данный потенциал реализован в файле Cu_u3.eam, содержащемся в библиотеке потенциалов [1]. Монокристаллическая медь имеет гранецентрированную кубическую структуру (face-centered cubic — fcc) с постоянной решетки 3.615 ˚A. Для определения тензора упругих модулей гранецентрированной кубической меди был проведен ряд расчетов в пакете LAMMPS [2]. Моделировался образец кубической формы, имеющий размеры 36.15х36.15х36.15˚A и содержащий 4000 атомов, к которому прикладывались различные деформации. Температура системы равна 0.1 К и поддерживается с помощью термостата Нойза-Гувера в каноническом NVT ансамбле. В результате расчета сохранялись значения потенциальной энергии и деформаций. Потенциальная энергия в LAMMPS вычисляется в эВ, а объем в ˚A3. Выполнен пересчет энергии в ГПа, а объема в м3. Затем аппроксимировали график зависимости энергии от деформации полиномиальной кривой. Используя коэффициенты данной кривой, проведено вычисление компонентов тензора упругих модулей. На рис. 2.1 приведен график зависимости потенциальной энергии от деформации, полученный для вычисления модуля C11. Точками отмечены результаты молекулярно-динамического расчета, прямая линия – это аппроксимирующая полиномиальная кривая.

На рис. 2.1 приведено уравнение полиномиальной кривой, из которой получим β = 24645. Подставив данное число в уравнение 1.9, получим C11 = 167.25 ГПа Затем прикладываем к данному образцу растягивающую деформацию вдоль оси х и сжимающую вдоль у. Получаем график зависимости

потенциальной энергии от деформации. Строим новую аппроксимирующую кривую и вычисляем

компоненту тензора упругих модулей C12. На рис. 2.2 приведено уравнение полиномиальной кривой, из которой найдем β = 12709. Подставив β в уравнение 1.12, получим C12 = 123.81 ГПа. Для вычисления компоненты тензора упругих модулей C44 к образцу прикладывается сдвиговая нагрузка. Аналогично строим график зависимости потенциальной энергии от деформации и ее аппроксимирующую функцию. По формуле 1.15 вычисляем коэффициент C44 = 76.46 ГПа. Таким образом, получили тензор упругих модулей, содержащий 3 независимых упругих константы, что соответствует кубической кристаллической

решетке.

 

	167.25	123.81	123.81	0	0	0	
	123.81	167.25	123.81	0	0	0	
 	123.81	123.81	167.25	0	0	0	 
 	0	0	0	76.46	0	0	 . 
 	0	0	0	0	76.46	0	 
	0	0	0	0	0	76.46

 

(2.1)

 

Для определения компонент тензора упругих модулей существует код в LAMMPS, написанный Эйданом Томпсоном, который может быть модифицирован для других типов решетки и потенциалов. Автор создает блок моделирования размерами 20х20х20 ячеек моделирования при нулевой температуре

Вестник Самарского университета. Естественнонаучная серия. 2021. Том 27, № 4. С. 14–21

Vestnik of Samara University. Natural Science Series. 2021, vol. 27, no. 4, pp. 14–21 17

 

Рис. 2.1. График зависимости потенциальной энергии системы от деформации и ее аппроксимация

полиномиальной кривой, построенный для вычисления компоненты тензора упругости C11

Fig. 2.1. Graph of the dependence of the potential energy of the system on deformation and its approximation of a polynomial curve constructed to calculate the elasticity tensor component C11

 

 

Рис. 2.2. График зависимости потенциальной энергии системы от деформации и ее аппроксимация

полиномиальной кривой, построенный для вычисления компоненты тензора упругости C12

Fig. 2.2. Graph of the dependence of the potential energy of the system on deformation and its approximation of a polynomial curve constructed to calculate the elasticity tensor component C12

 

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

Известны несколько схем определения упругих констант: Фойгта и Ройса. Общая формула, выражающая зависимость модуля упругости и модуля сдвига по Фойгту:

9BV = (C11 + C22 + C33) + 2(C12 + C23 + C31) ,

15GV = (C11 + C22 + C33) − (C12 + C23 + C31) + 3(C44 + C55 + C66). (2.2)

По Ройсу:

 

15

1

 

BR = (S11 + S22 + S33) + 2(S12 + S23 + S31) ,

GR = 4(S11 + S22 + S33) − 4(S12 + S23 + S31) + 3(S44 + S55 + S66). (2.3)

Часто используют аппроксимацию Фойгта — Ройса — Хилла, где модуль сдвига, модуль упругости и изотермический коэффициент всестороннего сжатия имеют простой вид

1

G = 2 (GV + GR), B =

1

2 (BV + BR), τ =

1

. (2.4)

B

О.Н. Белова Определение механических свойств материала, моделируемого с помощью метода молекулярной...

18 Belova O.N. Determination of mechanical properties of a material modeled using the molecular dynamics method

 

Рис. 2.3. График зависимости потенциальной энергии системы от деформации и ее аппроксимация

полиномиальной кривой, построенный для вычисления компоненты тензора упругости C44

Fig. 2.3. Graph of the dependence of the potential energy of the system on deformation and its approximation of a polynomial curve constructed to calculate the elasticity tensor component C44

 

Средние значения коэффициента Пуассона ν и модуля Юнга E определяются через G и B

ν = 3B − 2G

2(3B + G)

, E =

9BG

3B + G

 

. (2.5)

Для различных кристаллических решеток количество независимых констант различается, поэтому для каждой решетки эти формулы будут уникальны. В кубической решетке три линейно независимые константы: C11, C12, C44. Модуль всестороннего сжатия и модуль сдвига по Фойгту и Ройсу вычисляются по следующим формулам:

C11 + 2C12

C11 − C12 + 3C44

5(C11 − C12)C44

BV = BR =

, GV =

3

, GR =

5 4C44

+ 3(C11

− C12

. (2.6)

)

Для рассматриваемого материала получим модуль всестороннего сжатия BV = BR = 138.29 ГПа и модуль сдвига по Фойгту и Ройсу GV = 54.56 ГПа, GR = 38.07 ГПа, соответственно. По формулам Хилла получим модуль сдвига G = 46.32 ГПа, модуль всестороннего сжатия B = 138.29 ГПа, коэффициент Пуассона ν = 0.35 и модуль Юнга E = 125 ГПа.

Для визуализации упругих свойств можно применить расчетно-графический пакет ELATE – Elastic

tensor analysis. Данная оболочка представлена на сайте [3]. Пользователь вводит матрицу упругих модулей и получает средние значения свойств, двух- и трехмерную визуализацию упругих свойств и др.

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

C11 + C12

−C12

1

(C

 

S11 =

11

− C12

)(C11

+ 2C12

, S12 =

) (C11

− C12

)(C11

+ 2C12

. S44 =

) C44

(2.7)

Таким образом, получили тензор податливости в следующем виде:

16.15 · 10−3	−6.86 · 10−3	−6.86 · 10−3	0	0	0
−6.86 · 10−3	16.15 · 10−3	−6.86 · 10−3	0	0	0
−6.86 · 10−3	−6.86 · 10−3	16.15 · 10−3	0	0	0
0	0	0	13.07 · 10−3	0	0
0	0	0	0	13.07 · 10−3	0
0	0	0	0	0 13.07 ·

 

















 

10−3

 









 

 . (2.8)







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

1  sinθcosφ 

E =

(aiaj ak alSijkl)

, a =  sinθsinφ  , (2.9)

cosθ

Вестник Самарского университета. Естественнонаучная серия. 2021. Том 27, № 4. С. 14–21

Vestnik of Samara University. Natural Science Series. 2021, vol. 27, no. 4, pp. 14–21 19

 

1

G =

(4aibj ak blSijkl)

 cosθcosφcosγ − sinθsinγ

, b =  cosθsinφcosγ + cosθsinγ

−sinθcosγ

aiaj bk blSijkl



 , (2.10)

G = − (a a a a S

. (2.11)

)

i j k l

ijkl

Задавая координаты x, y, z радиус-вектора в любом направлении, согласно вычисленным компонентам тензора податливости откладывается относительное удлинение упругих модулей. Полученные картины иллюстрируют относительное удлинение модулей в каждом направлении для гранецентрированной кубической меди. На рис. 2.4 приведены трехмерные представления упругих констант, полученные из результатов работы электронного ресурса ELATE [3].

 

  

Рис. 2.4. Визуализация модуля Юнга, модуля сдвига, коэффициента Пуассона монокристалла меди

Fig. 2.4. Visualization of Young’s modulus, shear modulus, Poisson’s ratio of copper single crystal

 

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

 

 

	107.87	81.44	81.44	0	0	0	
	81.44	107.87	81.44	0	0	0	
 	81.44	81.44	107.87	0	0	0	 
 	0	0	0	46.16	0	0	 . 
 	0	0	0	0	46.16	0	 
	0	0	0	0	0	46.16

 

(2.12)

 

По формулам Хилла, используя найденные компоненты тензора напряжений, получили модуль Юнга

E = 76.24, коэффициент Пуассона ν = 0.36, модуль всестороннего сжатия B = 90.25, модуль сдвига G =

= 28.05. На рис. 2.5 приведены трехмерные представления упругих констант алюминия, полученные из результатов работы электронного ресурса ELATE.

 

  

Рис. 2.5. Визуализация модуля Юнга, модуля сдвига, коэффициента Пуассона монокристалла алюминия

Fig. 2.5. Visualization of Young’s modulus, shear modulus, Poisson’s ratio of aluminum single crystal

 

Найденные значения упругих модулей меди были использованы в работах [5; 6], посвященных определению коэффициентов разложения Уильямса в окрестности вершины трещины.

О.Н. Белова Определение механических свойств материала, моделируемого с помощью метода молекулярной...

20 Belova O.N. Determination of mechanical properties of a material modeled using the molecular dynamics method

 

Выводы

В статье приведена процедура определения упругих констант материала, моделируемого в пакете LAMMPS, реализующем метод молекулярной динамики. В результате комплекса расчетов найдены матрицы упругих постоянных Cij или коэффициентов податливости Sij fcc монокристаллов меди и алюминия. Были рассчитаны такие свойства материалов, как модуль Юнга, модуль сдвига, модуль всестороннего сжатия и коэффициент Пуассона. Для визуализации упругих свойств применен расчетно-графический пакет ELATE.

 

Об авторах

О. Н. Белова

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

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

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

Россия

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

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