CC BY
36
Челябинский физико-математический журнал. 2018. Т. 3, вып. 3. С. 344-352.
УДК 51-73 БОТ: 10.24411/2500-0101-2018-13308
АЛГОРИТМ ДЛЯ АНАЛИЗА РАСПРЕДЕЛЕНИЯ
ПОР ПО РАЗМЕРАМ, ОСНОВАННЫЙ НА РЕЗУЛЬТАТАХ
МОЛЕКУЛЯРНО-ДИНАМИЧЕСКОГО
МОДЕЛИРОВАНИЯ
А. Е. Майер", П. Н. Майерь
Челябинский государственный университет, Челябинск, Россия "mayer@csu.ru, ьpolina.nik@mail.ru
Разработан алгоритм для поиска и определения размеров пор в конденсированной фазе по набору координат центров атомов, получаемому в результате молекулярно-динамического моделирования. Алгоритм основан на разбиении расчётной ячейки на подъячейки с размерами меньше среднего межатомного расстояния, которые считаются пустыми или заполненными в зависимости от расстояний до центров ближайших атомов. Проведены тестирование алгоритма и анализ с его помощью распределения полостей по размерам в расплаве алюминия при высокоскоростном растяжении со скоростью деформации 100/нс при температуре 1100 К. Показано, что на ранних стадиях образования и роста пор их распределение по размерам близко к экспоненциальному; отклонение от экспоненциального распределения имеет место для пор наибольших размеров.
Ключевые слова: молекулярно-динамическое моделирование, конденсированное вещество, полости, 'распределение по размерам, высокоскоростное 'растяжение.
Введение
Образование и рост пор в конденсированном веществе сопровождает процессы вязкого разрушения металлов [1] и кавитации в жидкостях [2], включая расплавы металлов. Образование пор чаще всего происходит при отрицательных давлениях, когда межатомные расстояния больше равновесных значений. Такие состояния ме-тастабильны, т. е. неустойчивы по отношению к образованию полостей, появление которых позволяет окружающим объёмам вещества сжаться до равновесной плотности. В ходе распада растянутого метастабильного состояния идут параллельные процессы образования и роста пор, что приводит к появлению распределения пор по размерам. Информация о характере такого распределения и его эволюции важна для построения моделей разрушения твёрдого и жидкого вещества в растягивающих напряжениях.
Метод классической молекулярной динамики (МД) часто используется для исследования динамического разрушения твёрдых металлов [3-6] и расплавов [7-9] при высоких скоростях деформации. Непосредственной информацией, получаемой при МД-моделировании, являются положения центров атомов, их скорости и силы, действующие между ними. По этим данным могут быть восстановлены напряжения, в частности, предельные напряжения, приводящие к образованию полостей. Для определения размеров пор, возникающих в МД-системе при растяжении,
Исследование выполнено за счёт гранта Российского научного фонда (проект 17-71-10205).
необходимо по положениям центров атомов находить занятые и пустые области пространства. Можно выделить два подхода к этому. Первый состоит в разбиении пространства на тетраэдры Делоне, в вершинах которых находятся центры атомов. Пустыми считаются тетраэдры, не вписывающиеся в сферу порогового радиуса [10], т.е. такие, расстояние между атомами в вершинах которых достаточно велико. В качестве порогового радиуса обычно берётся положение первого максимума парно-корреляционной функции. Данный подход реализован в алгоритме Construct Surface Mesh программного комплекса OVITO [11]. Второй подход состоит в разбиении всего пространства, занятого МД-системой, на регулярно расположенные прямоугольные ячейки (воксели) [7]; заполненными считаются воксели, находящиеся от центров атомов на расстоянии не далее порогового радиуса.
В данной статье описана реализация алгоритма на основе разбиения объёма МД-системы на воксели, его тестирование и применение для анализа распределения полостей по размерам в расплаве алюминия при высокоскоростном растяжении.
1. Алгоритм поиска пор
Для исследования распределения пор по размерам выбираются такие условия растяжения (скорость деформации и размер системы), при которых в МД-системе образуется и одновременно растёт сразу множество пор. Поэтому задачей алгоритма обработки является поиск всех изолированных друг от друга пор в МД-системе и вычисление их объёмов, каждой поры по отдельности. Входными данными являются положения центров всех атомов. Атомы находятся внутри области пространства в виде прямоугольного параллелепипеда (расчётной ячейки), на гранях которого заданы периодические граничные условия. Рёбра расчётной ячейки параллельны осям Ox, Oy и Oz и имеют длины Lx, Ly и Lz соответственно.
Расчётная ячейка разбивается сеткой подъячеек (вокселей) также в форме прямоугольных параллелепипедов с рёбрами hx, hy и hz. Длины рёбер выбираются так, что отношение La/ha является целым и ha < а/и, где параметр а представляет собой пороговый радиус, целое число и контролирует степень подробности сетки вокселей, а = x,..., z. В качестве а используется значение, близкое к положению первого максимума парно-корреляционной функции для рассматриваемой МД-системы. Тестовые расчёты показывают, что значение и =1 (близко к тому, что использовалось в [7]) даёт недостаточно подробную сетку, а при и > 1 число вокселей существенно превышает число атомов в МД-системе и растёт как и3. Поэтому нами используется и = 2. Состояние сетки вокселей хранится в трёхмерном массиве с минимальным типом данных: значение 1 соответствует заполненному вокселю, значение 0 — пустому вокселю, ещё не приписанному ни к одной поре, значение 2 — пустому вокселю, уже приписанному к одной из пор.
Воксель считается заполненным, если на расстоянии r < а от одной из его вершин находится центр хотя бы одного атома. Для заполнения вокселей организуется цикл по всем атомам МД-системы. Вначале по координатам атома {x; y; z} вычисляются индексы {i0; j0; k0} вокселя, в котором лежит его центр, например i0 = [(x — Xi)/hx] + 1, где {xi; yi; zi} — координаты начальной точки расчётной ячейки, квадратные скобки означают взятие целой части. Далее пробегаются все воксели в диапазоне {i = i0 — и,..., i0 + и; j = j0 — и,... ,j0 + и; k = k0 — и,..., k0 + и} и проверяется расстояние r от атома до ближайшей вершины вокселя, положение которой вычисляется по индексам {i,j,k}. Если выполняется условие r < а, состояние вокселя менятеся на 1 (заполнен). При этом учитываются периодические граничные условия.
Далее происходит идентификация отдельных полостей, для чего в цикле ищутся пустые воксели, объединённые в полость. Условием выхода из цикла является то, что все пустые воксели связаны в полости. Два пустых вокселя считаются соединёнными друг с другом, если у них есть общая грань. Все соединённые друг с другом пустые воксели формируют отдельную полость, объём которой Vm равен сумме объёмов вокселей Vm = (hxhyhz)Nm, где Nm — число вокселей в m-й полости. Поиск каждой поры начинается с поиска пустого вокселя, не приписанного ни к одной поре. Далее пробегаются воксели, граничащие гранями с данным, и среди них обнаруживаются пустые. Процедура повторяется рекурсивно, пока среди соседей, входящих в пору вокселей, не останется пустых вокселей. Поиск соседей проводится с учётом периодических граничных условий. Рекурсия реализуется при помощи стека, хранящего индексы входящих в пору вокселей с ещё не обработанными соседями. Стек пополняется при обнаружении новых вокселей, входящих в пору, и уменьшается при обработке соседей одного из вокселей стека. Пустой стек означает конец обработки поры. Такой алгоритм при простоте реализации обеспечивает приемлемую производительность, в том числе при обработке МД-систем с большой долей пор. В предположении сферичности поры эффективный радиус может быть вычислен как Rm = ((3Vm)/(4n))1/3. Поры с малым радиусом Rm < a не рассматриваются, поскольку соответствуют отдельным вакансиям или тепловым флуктуациям.
2. Тестирование алгоритма
Тестирование алгоритма проводилось на МД-образцах с заранее вырезанными порами. Атомы размещались в узлах правильной ГЦК-решётки с параметром 4.05 А, что соответствует алюминию; первый максимум парно-корреляционной функции в этом случае находится на расстоянии a = 2.88 А_. Затем все атомы, центры которых расположены внутри сферы радиуса R, удалялись, полученный набор координат анализировался разработанным алгоритмом. Результаты для единичной поры приведены в таблице. Видно, что оценённый алгоритмом радиус поры R1 систематически больше заданного на величину AR ^ 1.1 ^ 1.2 А, что составляет порядка 0.4a вне зависимости от R. Можно было бы добавить в алгоритм простую корректировку, но разность SR является не ошибкой, а возникает из-за различий в определении того, что называть пустым объёмом. При подготовке образца с порой удаляются все атомы, центры которых попали внутрь сферы; т. е. пустым считается объём, в котором нет центров атомов. В алгоритме обработки пустыми считаются области, расположенные далее a от центра любого атома. Последнее в большей степени соответствует физическому определению пустого объёма, как области, в которую может быть помещён пробный атом, который из-за действия отталкивающих сил не может подойти близко к центру другого атома. Поэтому мы не вносили в алгоритм дополнительную корректировку. Проблема невозможности точного определения объёмов и площадей групп атомов хорошо известна в физике наноструктур. При росте размера поры относительная разность AR/R ^ 0.
Помимо единичных пор тестирование проводилось на системах с несколькими
Результаты тестирования алгоритма: Я — заданный радиус поры, Я\ — вычисленный радиус поры, ДЯ = Я —
R A Ri A AR A
8.10 6.93 1.17
12.15 11.10 1.05
16.2 15.00 1.20
20.25 19.15 1.10
24.30 23.07 1.23
28.35 27.16 1.19
32.40 31.29 1.11
36.45 35.35 1.10
40.50 39.38 1.12
44.55 43.36 1.19
48.60 47.44 1.16
52.65 51.48 1.17
56.70 55.61 1.09
порами. Также рассматривались системы с порами, пересекающими периодическую границу, т. е. когда части поры расположены вблизи противоположных граней расчётной ячейки. В этих случаях алгоритм определял размер каждой поры в МД-системе с точностью, соответствующей таблице.
3. Распределение по размерам пор в расплаве алюминия при высокоскоростном растяжении
В качестве примера применения алгоритма анализировалось распределение по размерам полостей в расплаве алюминия при растяжении со сверхвысокой скоростью деформации 100/нс. МД-моделирование проводилось при помощи пакета LAMMPS [12] и межатомного потенциала [13], основанного на модели погружённого атома (EAM) [14]. Результаты расчётов визуализировались при помощи программы OVITO [11] и анализировались разработанным нами алгоритмом.
Рассматривался образец алюминия в форме куба, содержащий 4 миллиона атомов. Атомы изначально помещались в узлы правильной ГЦК-решётки, затем образец переводился в состояние расплава нагревом до температуры 2500 К с постепенным охлаждением до 1100 К в течение 10 пс. Интегрирование по времени проводилось с шагом 0.001 пс. На начальном этапе подготовки расплава заданная температура поддерживалась термостатом Нозье — Гувера, а нулевое давление в системе поддерживалось баростатом Нозье — Гувера. После этапа подготовки моделировалось однородное всестороннее растяжение расплава масштабированием координат атомов при помощи стандартной команды deform пакета LAMMPS. Масштабирование координат физически соответствует растяжению расплава по инерции. Рассматривалась скорость деформации 100/нс. В ходе растяжения температура системы поддерживалась термостатом на уровне 1100 К. Постановка задачи соответствует статье [9], но здесь при помощи разработанного алгоритма дополнительно анализировалось распределение по размерам возникающих в расплаве полостей, в то время как в [9] только лишь приближённо оценивалось их количество.
(а) (Ь)
3 4 5 6 7 8
Pore size (angstrom)
Рис. 1. Расплав алюминия на стадии роста числа полостей: общий вид МД-системы (а) с раскраской атомов по значению энергии, шкала в эВ, и гистограмма распределения полостей по размерам (Ъ), общее число пор в системе равно 8540. Скорость деформации 100/нс, температура 300 К, момент времени 2пс с начала растяжения
(а) (Ь)
5 10 15 20 25 30 35 Pore size (angstrom)
Рис. 2. Расплав алюминия на стадии уменьшения числа полостей: общий вид МД-системы (а) с раскраской атомов по значению энергии, шкала в эВ, и гистограмма распределения полостей по размерам (Ь), общее число пор в системе равно 3071. Скорость деформации 100/нс, температура 300 К, момент времени 5пс с начала растяжения
Как следует из [9] и проведённых здесь расчётов, в эволюции растягиваемого расплава можно выделить стадию активной нуклеации полостей, когда после достижения некоторого порога по давлению их число быстро увеличивается на фоне продолжающегося уменьшения давления. Результаты расчётов для данной стадии представлены на рис. 1. Согласно рис. 1 (а) распределение полостей по размерам близко к экспоненциальному в диапазоне радиусов 3-6 А. В то же время для самых крупных полостей, появившихся в начале процесса нуклеации, имеет место существенное отклонение: их количество меньше, чем следовало бы из аппроксимации экспоненциальным распределением.
(а) (Ь)
Pore size (angstrom)
Рис. 3. Расплав алюминия на стадии объединения полостей: общий вид МД-системы (а) с раскраской атомов по значению энергии, шкала в эВ, и гистограмма распределения полостей
по размерам (Ь), общее число пор в системе равно 490. На (Ь) появление большой полости с эффективным радиусом 311А означает разрушение перегородок между полостями, что можно видеть на (а). Скорость деформации 100/нс, температура 300 К, момент времени 16пс с начала
растяжения
Падение давления в МД-системе останавливается, когда скорость роста объёма полостей начинает превышать скорость роста объёма системы, что приводит к остановке процесса нуклеации новых полостей. Далее количество полостей сокращается за счёт схлопывания мелких полостей при одновременном росте крупных. Результаты расчётов для этой стадии представлены на рис. 2. За счёт схлопывания мелких полостей распределение выполаживается, за счёт роста крупных — уширяется, но всё равно является близким к экспоненциальному распределению в диапазоне 3-25 A. Общее количество полостей в момент времени, представленный на рис. 2, более чем вдвое сократилось по сравнению с рис. 1.
При достижении объёмной доли полостей величины порядка 0.9 начинается их объединение за счёт разрушения перегородок между полостями. Распределение полостей по размерам для стадии объединения представлено на рис. 3. В распределении (рис. 3 (b)) появляется отдельно стоящая крупная пора, которая соответствует полостям с разрушенными перегородками. Форма пор на этой стадии уже далека от сферической (рис. 3(a)), но эффективный радиус всё равно является удобной характеристикой их размера.
Заключение
Разработан алгоритм для поиска и определения размеров полостей в конденсированной фазе по набору координат центров атомов, получаемому в результате молекулярно-динамического моделирования. Алгоритм основан на разбиении расчётной ячейки на подъячейки с размерами меньше среднего межатомного расстояния, которые считаются пустыми или заполненными в зависимости от расстояний до центров ближайших атомов. Проведены тестирование алгоритма и анализ с его помощью распределения полостей по размерам в расплаве алюминия при высокоскоростном растяжении со скоростью деформации 100/нс при температуре 1100 К. Показано, что на ранних стадиях образования и роста полостей их распределение по размерам близко к экспоненциальному; отклонение от экспоненциального распределения имеет место для полостей наибольших размеров.
Список литературы
1. Tang, Y. Ductile tensile failure in metals through initiation and growth of nanosized voids / Y. Tang, E. M. Bringa, M. A. Meyers // Acta Materialia. — 2012. — Vol. 60. — P. 4856-4865.
2. Theory and molecular dynamics modeling of spall fracture in liquids / A. Yu. Kuksin, G. E. Norman, V. V. Pisarev, V. V. Stegailov, A. V. Yanilkin // Physical Review B. — 2010. — Vol. 82. — P. 174101.
3. Dynamic fracture kinetics, influence of temperature and microstructure in the atomistic model of aluminum / A. Kuksin, G.Norman, V. Stegailov, A. Yanilkin, P. Zhilyaev // International Journal of Fracture. — 2010. — Vol. 162. — P. 127-136.
4. Pogorelko, V. V. Influence of titanium and magnesium nanoinclusions on the strength of aluminum at high-rate tension: Molecular dynamics simulations / V. V. Pogorelko, A. E. Mayer // Materials Science and Engineering: A. — 2016. — Vol. 662. — P. 227-240.
5. Mayer, A. E. Influence of free surface nanorelief on the rear spallation threshold: Molecular dynamics investigation / A. E. Mayer, A. A. Ebel // Journal of Applied Physics. — 2016. — Vol. 120, no. 16. — P. 165903.
6. On the ultimate tensile strength of tantalum / E. N. Hahn, T. C. Germann, R. Ravelo, J. E. Hammerberg, M.A.Meyers // Acta Materialia. — 2017. — Vol. 126. — P. 313-328.
7. Cai, Y. Cavitation in a metallic liquid: Homogeneous nucleation and growth of nanovoids / Y. Cai, H.A.Wu, S. N. Luo // The Journal of Chemical Physics. — 2014. — Vol. 140. — P. 214317.
8. Mayer, A. E. Continuum model of tensile fracture of metal melts and its application to a problem of high-current electron irradiation of metals / A. E. Mayer, P. N. Mayer // Journal of Applied Physics. — 2015. — Vol. 118, no. 3. — P. 035903.
9. Mayer, P. N. Late stages of high rate tension of aluminum melt: Molecular dynamic simulation / P. N. Mayer, A. E. Mayer // Journal of Applied Physics. — 2016. — Vol. 120, no. 7. — P. 075901.
10. Stukowski, A. Computational analysis methods in atomistic modeling of crystals /
A. Stukowski // JOM. — 2014. — Vol. 66, no. 3. — P. 399-407.
11. Stukowski, A. Visualization and analysis of atomistic simulation data with OVITO — the Open Visualization Tool / A. Stukowski // Modelling and Simulation in Materials Science and Engineering. — 2010. — Vol. 18. — P. 015012.
12. Plimpton, S. Fast parallel algorithms for short-range molecular dynamics / S.Plimpton // Journal of Computational Physics. — 1995. — Vol. 117, no. 1. — P. 1-19.
13. Zope, R. R. Interatomic potentials for atomistic simulations of the Ti-Al system / R. R. Zope, Y. Mishin // Physical Review B. — 2003. — Vol. 68. — P. 024102.
14. Daw, M. S. Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals / M. S. Daw, M. I. Baskes // Physical Review
B. — 1984. — Vol. 29. — P. 6443.
Поступила в редакцию 12.04-2018 После переработки 03.08.2018
Сведения об авторах
Майер Александр Евгеньевич, доктор физико-математических наук, заведующий кафедрой общей и прикладной физики, Челябинский государственный университет, Челябинск, Россия; e-mail: mayer@csu.ru.
Майер Полина Николаевна, кандидат физико-математических наук, старший научный сотрудник кафедры общей и прикладной физики, Челябинский государственный университет, Челябинск, Россия; e-mail: polina.nik@mail.ru.
Chelyabinsk Physical and Mathematical Journal. 2018. Vol. 3, iss. 3. P. 344-352.
DOI: 10.24411/2500-0101-2018-13308
ALGORITHM FOR ANALYSIS OF PORE SIZE DISTRIBUTION BASED ON RESULTS OF MOLECULAR DYNAMIC SIMULATIONS
A.E. Mayer", P.N. Mayerb
Chelyabinsk State University, Chelyabinsk, Russia "mayer@csu.ru; bpolina.nik@mail.ru
We develop an algorithm for searching of pores in the condensed phase and determination of their sizes from the set of coordinates of the atomic centers obtained as a result of the molecular dynamics simulation. The algorithm is based on dividing the simulation cell into subcells with dimensions less than the average interatomic distance; the subcells are considered empty or filled depending on the distances to the centers of the nearest atoms. The algorithm is tested and applied for analyzing the size distribution of cavities in aluminum melt under high rate tension at the strain rate of 100/ns at the temperature of 1100 K. It is shown that at early stages of pore formation and growth, the size distribution of pores is close to the exponential one; the deviation from the exponential distribution occurs for largest pores.
Keywords: molecular-dynamic simulations, condensed matter, cavities, size distribution, high rate tension.
References
1. TangY., BringaE.M., Meyers M.A. Ductile tensile failure in metals through initiation and growth of nanosized voids. Acta Materialia, 2012, vol. 60, pp. 4856-4865.
2. KuksinA.Yu., Norman G.E., PisarevV.V., Stegailov V.V., YanilkinA.V.
Theory and molecular dynamics modeling of spall fracture in liquids. Physical Review B, 2010, vol. 82, p. 174101.
3. KuksinA., Norman G., Stegailov V., YanilkinA., ZhilyaevP. Dynamic fracture kinetics, influence of temperature and microstructure in the atomistic model of aluminum. International Journal of Fracture, 2010, vol. 162, pp. 127-136.
4. Pogorelko V.V., Mayer A.E. Influence of titanium and magnesium nanoinclusions on the strength of aluminum at high-rate tension: Molecular dynamics simulations. Materials Science and Engineering: A, 2016, vol. 662, pp. 227-240.
5. Mayer A.E., EbelA.A. Influence of free surface nanorelief on the rear spallation threshold: Molecular dynamics investigation. Journal of Applied Physics, 2016, vol. 120, no. 16, p. 165903.
6. HahnE.N., GermannT.C., RaveloR., Hammerberg J.E., Meyers M.A. On the
ultimate tensile strength of tantalum. Acta Materialia, 2017, vol. 126, pp. 313-328.
7. CaiY., WuH.A., LuoS.N. Cavitation in a metallic liquid: Homogeneous nucleation and growth of nanovoids. The Journal of Chemical Physics, 2014, vol. 140, p. 214317.
8. Mayer A.E., Mayer P.N. Continuum model of tensile fracture of metal melts and its application to a problem of high-current electron irradiation of metals. Journal of Applied Physics, 2015, vol. 118, no. 3, p. 035903.
9. Mayer P.N., Mayer A.E. Late stages of high rate tension of aluminum melt: Molecular dynamic simulation. Journal of Applied Physics, 2016, vol. 120, no. 7, p. 075901.
10. Stukowski A. Computational analysis methods in atomistic modeling of crystals. JOM, 2014, vol. 66, no. 3, pp. 399-407.
11. Stukowski A. Visualization and analysis of atomistic simulation data with OVITO — the Open Visualization Tool. Modelling and Simulation in Materials Science and Engineering, 2010, vol. 18, p. 015012.
12. Plimpton S. Fast parallel algorithms for short-range molecular dynamics. Journal of Computational Physics, 1995, vol. 117, no. 1, pp. 1-19.
13. ZopeR.R., MishinY. Interatomic potentials for atomistic simulations of the Ti-Al system. Physical Review B, 2003, vol. 68, p. 024102.
14. Daw M.S., BaskesM.I. Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals. Physical Review B, 1984, vol. 29, p. 6443.
Accepted article received 12.04-2018 Corrections received 03.08.2018
