АППРОКСИМАЦИЯ ПЕРИОДИЧЕСКИХ РЕШЕНИЙ В ДВУХМОДОВОЙ МОДЕЛИ КРИСТАЛЛИЧЕСКОГО ФАЗОВОГО ПОЛЯ Текст научной статьи по специальности «Физика»

ISSN 2079-3316 ПРОГРАММНЫЕ СИСТЕМЫ: ТЕОРИЯ И ПРИЛОЖЕНИЯ т. 13, №2(53), с. 65-84

научная статья _ математическое моделирование

УДК 519.688:519.63

10.25209/2079-3316-2022-13-2-65-84;

*>

Аппроксимация периодических решений в двухмодовой модели кристаллического фазового

поля

Владимир Евгеньевич Анкудинов1", Илья Олегович Стародумов2 Институт физики высоких давлений им. Л. Ф. Верещагина РАН, Москва, Россия1, Уральский федеральный университет им. Б.Н. Ельцина,

г. Екатеринбург, Россия2,

vladimir@ankudinov.orgM (подробнее об авторах на с. 81)

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

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

(see abstract in English on p. 82)

Ключевые слова и фразы: метод кристаллического фазового поля, численные расчеты, конечные элементы, аппроксимация

Благодарности:

Исследование выполнено при поддержке гранта Российского научного фонда № 21-73-00263, https://rscf.ru/project/21-73-00263/.

Для цитирования: Анкудинов В. Е., Стародумов И. О. Аппроксимация периодических решений в двухмодовой модели кристаллического фазового поля // Программные системы: теория и приложения. 2022. Т. 13, № 2(53). С.65-84. http://psta.psiras.ru/read/psta2022_2_65-84.pdf

© Анкудинов В. Е., Стародумов И. О.

2[юШаи0

Введение

Модель кристаллического фазового поля или Phase-Field Crystal (КФП, PFC) была сформулирована для описания переходов из однородного состояния в периодическое кристаллическое, а также между различными периодическими состояниями на временах, сопоставимых со временами диффузии [1—3]. КФП - континуальная модель, учитывающая симметрию периодической фазы, она основывается на описании свободной энергии в виде функционала поля атомной плотности (концентрации) [2]. Уравнение КФП является консервативной версией модели Свифта-Хоэнберга [4] и используется при моделировании упорядочения структур на микронных масштабах [5], движения фронтов кристаллизации и формирования дендритов [6].

Ранее, в наших работах была рассмотрена модифицированная модель КФП [7], разработана и исследована ее неявная численная постановка, основанная на теории изогеометрического анализа (IGA) [8]. Отличие представленной в данной работе реализации заключается в использовании метода конечных элементов в постановке Галерки-на [9]. Для вычисления временных производных использовался метод дифференцирования назад (BDF, Backward differentiation formula) до 5-го порядка точности.

1. Уравнения модели кристаллического фазового поля

Свободная энергия в модели КФП, позволяющая описывать фазовые переходы и кристаллизацию в периодическую фазу для системы в некоторой области Q, имеет следующий вид [10]:

где п - параметр порядка, мгновенное поле плотности системы; а, V -коэффициенты свободной энергии в форме Ландау, отвечающие за межфазный барьер; £д - дифференциальный оператор, полученный с помощью приближенного градиентного разложения обменного вклада [1,11]. Этот оператор может быть представлен в одномодовом приближении при К = 1 как:

(1)

(2)

Li = ДВс + B0x(V2 + q02)2.

Для описания более сложных структур, таких как ГЦК или гексагональные решетки требуется двухмодовое разложение более высокого порядка, К =2:

(3) £ = ДВо + ВХ(го + (V2 + 92)2)(Г1 + (V2 + д2)2).

Здесь коэффициенты г0, г1, до, отвечают за параметры корреляционной функции, определяющей взаимодействие и формирование периодического решения. Коэффициенты соответствуют значе-

ниям параметров решетки первой и второй координационных сфер. Коэффициенты го, Г1 напрямую влияют на стабильность воспроизводимых периодических структур [12]. Движущая сила ДВо = Вд — Вд определена как разница между модулем объемной упругости жидкости Вд и модулем упругости кристалла Вд. Управляющий параметр ДВо может рассматриваться как переохлаждение ДТ, то есть отклонение от температуры плавления/кристаллизации Тт:

(4) ДВо = ДВ* (1 + ДТ) .

Здесь ДВ* — критическая движущая сила соответствующая температуре плавления Тт, которую можно выразить с помощью коэффициентов уравнения (1) как ДВ* = 2а2/(%) [10,13]. При ДВо = ДВ* достигается равенство энергий жидкой и кристаллической фазы [13].

Первоначально модель КФП была сформулирована в параболической форме [1], которая позволяет описывать динамику системы вблизи критической точки в приближении локального термодинамического равновесия. Для учета высокоскоростных переходов и конечного времени релаксации т потока плотности Л динамическое уравнение КФП было расширено второй производной по времени [4,14]. Данное расширение (МКФП модифицированная КФП) позволяет учитывать «быструю» упругую релаксацию одновременно с «медленной» концентрационной диффузией. Само динамическое уравнение имеет следующий вид:

. . д2п дп , гж_,2 ГдВ(п)

(5) + Ж = ()

дп

здесь М — мобильность фазового поля, которая устанавливает масштабы времени £ и является аналогом диффузионного коэффициента; правая часть уравнения включает функциональную производную

от свободной энергии уравнения (1), имеющую смысл химического потенциала ^(п) = №(п)/6п. Данное уравнение в частных производных учитывает невозрастание свободной энергии во времени. В двухмодовом приближении динамическое уравнение (5) после подстановки свободной энергии уравнение (1) с дифференциальным оператором уравнение (3) примет следующий вид:

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

Ранее для одномодового уравнения КФП с пространственными производными шестого порядка была разработана безусловно устойчивая схема решения [15], в которой, однако, граничное условие на поток после разделения переменных является мягким (требуется непрерывность производной потока J). В нашей работе мы будет использовать понижение пространственной производной, а производную по времени вычислять явно, с использованием предыдущих шагов по времени.

2. Численное решение двухмодового уравнения КФП

Численное решение динамического уравнения (6) с двухмодо-вым оператором уравнения (3) было выполнено с помощью метода конечных элементов (МКЭ) в пакете COMSOL Multiphysics [16] и Fenics (DOLFIN) [17] в прямом пространстве. При дискретизации использовались тетраэдральные (в трехмерном случае) и треугольные (в двухмерном) лагранжевые конечные элементы (КЭ) различной степени. При использовании ПО COMSOL Multiphysics решения выполнялись с помощью прямых решателей MUMPS и PARDISO; в пакете FeniCS (DOLFIN) применялось семейство решателей PetSC [18] с реализацией MUMPS. Основные расчеты выполнялись с применением

(6)

д2 n dn „ г„ „ о о

~тгрг + ^т =МV2 [ДВоп - an2 + vn3+ dt2 dt

+ B0x(ro + (V2 + q2)2)(n + (V2 + q2)2)n].

dt

двухпроцессорного сервера на базе AMD Epyc 7302, 64 ядра, 1024 Gb RAM, COMSOL Multiphysics, MUMPS.

Для понижения степени пространственных производных введем дополнительные переменные P2, P4, P6 и запишем уравнение (6) следующим образом:

д 2n dn „

т ^ + Ж = M V2 "

(7)

^ = (ДВ0 + BQ) n - an2 + vn3+

+BXV2 (2Q2n + Q3P2 + 2Q4P4 + Рб): P2 = V2n,

P4 = V2P2,

Рб

V2P4.

Здесь Q 1..4 - константы, в которые входят параметры q0, q 1, r0, r 1:

(8)

Qi = qo^i + qori + qiro + rori, Q2 = qote + q^ + ri) + q^o

Q3 = q° + + ro + ri, Q4 = q° + q°°.

Вторая производная по времени рассчитывалась при помощи интегрирования назад с точностью до пятого порядка. Для большей части расчетов, если не оговорено специально, в данной работе использовалась двумерная область с периодическими граничными условиями, заданными на протяженных стенках, размеры расчетной области Ьх = 200, Ьу = 100, параметры расчетов устанавливались как т = 0.1, М = 1, Вд =1, а = 0, V = 1. Начальный шаг по времени был задан как д£о = 10-6, который затем менялся адаптивно до д£тах = 10-1. Начальное распределение в расчетной области задавалось в виде постоянного значения поля п=по в основной части домена, при этом, для инициации фазового превращения, вблизи границы х = 0 задавалось возмущение в виде полоски амплитудой 1 шириной 1.

3. Отбор параметра решетки

3.1. Периодические решения в модели КФП

Решением в модели МКФП (6) является периодическая структура (кристалл) с периодом А (в двухмодовом случае при до = периодов

решетки может быть несколько [19]), либо константное распределение параметра порядка п. При этом в динамике можно наблюдать распространение фронта кристаллизации в переохлажденную гомогенную фазу.

На рисунке 1 распространение фронта кристаллизации происходит слева направо в переохлажденную жидкость, параметры расчета: т = 0.1, п0 = -0.45, ДВ0 = -1, Яо = 0.7, 91 = 1, г0 = -0.1, г1 = 0.1, Ьх = 500. После заполнения расчетного домена периодической фазой

1

0.5 0

-0.5 -1

1[ 0.5 0

-0.5

0 100 200 300 400 х

(б) в момент времени £ = 60

Рисунок 1. Распределение параметра порядка п при решении одномерной задачи КФП

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

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

0 100 200 300 400 *

(а) в момент времени t = 25

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

Чрезмерное количество КЭ может приводить к неприемлемой вычислительной сложности. В случае двухмодовой постановки задачи баланс размера КЭ особенно важен, так как после введения дополнительных независимых переменных Р2,4,6 даже в случае линейной С1 аппроксимации каждому узлу ассоциированы 5 степеней свободы.

Отбор периодических решений в моделях КФП и МКФП связан с существованием устойчивых кристаллических модификаций в виде периодических решений [10]. При этом аналитические решения позволяют найти минимальную равновесную энергию существования только конкретного кристалла (полосчатого, треугольного и т.д.). С помощью численных расчетов возможно найти не только равновесную конфигурацию для произвольных граничных условий, но и исследовать динамику структурных превращений в данной системе, например рассчитать скорость фронта или рассмотреть упорядочение на фронте [8, 20,21]. Метод конечных элементов предлагает приближенное решение задачи на ограниченном наборе конечных элементов с базисными функциями Сп.

При численных расчетах в модели МКФП существует два возможных сценария возникновения неравновесных конфигураций:

(1) неполная релаксация в системе, связанная с кинетическими эффектами (например связанная с присутствием в уравнении (6) второй производной по времени),

(2) невозможность интерполяции конкретного периодического распределения на данном наборе КЭ.

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

3.2. Отбор параметра решетки

В сравнительно больших системах, как на рисунке 1, рост периодической фазы не подвергается значительному влиянию граничных условий (размер области Ь ^ Л). Рассмотрим отбор равновесного параметра решетки при ненулевом временах релаксации. Влияние кинетического вклада, определяемого гиперболическим членом уравнения (6) можно обнаружить построив фурье-образ Т(п) или структурный фактор £(к) параметра порядка п в обратном пространстве к. Как видно из расчетов с параметрами п0 = -0.45, ДВ0 = -1, 9о = 0, 7, ^ = 1, го = -0.1, г 1 = 0.1, Ьх = 500 на рисунке 2, при т = 0 в случае

ЗД

160 120 80 40 0

(а) время релаксации потока плотности т = 0

т 120 100 80 60 40 20 0

0 0.5 1 1.5

(б) время релаксации потока плотности т = 0.1

Рисунок 2. Структурный фактор £(к) (фурье-образ) периодического решения п одномерной задачи КФП после частичной релаксации £ = 100

нулевого времени релаксации отбирается единственный период решетки, равный среднему между д0 и Л = (д0 + )/2. В случае т = 0.1 кинетический вклад оказывает влияние на отбор Л и пик структурного фактора промежуточной структуры расщепляется, что связано с неокончательной релаксацией к равновесному параметру решетки

вследствие возникающей неустойчивости Экхауса [20] на фронте, а затем и в объеме. При t = 200 пик остается один и система приходит к равновесному А.

При решении задачи МКФП выбор дискретизации расчетной области должен учитывать возможное существование вторичных гармоник, связанных с кинетическими эффектами, однако, как показано на рисунке 2 их частота в случае неустойчивости Экхауза меньше частоты основной структуры. Таким образом, лимитирующим фактором величины конечного элемента I является период основной структуры А.

4. Аппроксимация численного решения

Целью численных экспериментов было оценить влияние размера конечного элемента на качество аппроксимации решения задачи МКФП. Был подготовлен сценарий кристаллизации с распространением плоского фронта вдоль вытянутого вычислительного домена. Для экспериментов изменялись размеры конечных элементов I в домене, а также порядок аппроксимации элементов Cn, n = 1, 2. Такой набор выбирался исходя из возможности продемонстрировать сходимость решений и относительных соотношений КЭ при минимальном наборе численных расчётов.

Размер расчетной сетки выбирался исходя из параметра кристаллической решетки моделируемой структуры, отбираемого в одномерных расчетах А = 5.3, а именно, вблизи размера А/5 (в предположении, что для корректного описания одной ячейки требуется 5 узлов).

Исходя из рассуждений о качестве аппроксимации пиков атомной плотности на фронте кристаллизации, для случая плоского фронта определяющим параметром I должен быть его линейный размер вдоль направления движения фронта. Для случая одномодового уравнения МКФП при использовании IGA таким размером является I « 5..6 [8], однако, в случае двухмодовой модели МКФП в классической постановке МКЭ размеры элемента должны быть уменьшены, что связано с повышенными требованиями к точности вычисления пространственных производных. Требования к размеру КЭ в постановке IGA также ниже из-за наличия дополнительных степеней свободы в NURBS (Non-uniform rational В^рИпе)-базисе [7].

Расчеты были выполнены при до = 1, ^ = 1.18, Го = -0.22, Г1 = 0.08, что соответствовало равновесной квадратной решетке при начальных параметрах по = -0.55, ДВо = -1.

На рисунке 3 показаны решения двухмодового МКФП уравнения при различных порядках аппроксимации в пределах конечного элемента для величины КЭ I = 2.

С1 линейные элементы момент с2 квадратичные элементы

времени

(а)

(с)

(е)

I = 50 (Ь)

I = 100 ((1)

I = 2000

ЩЁмтштшж

о О О О о „ООО ООО О О О О

0л°0 0°о0л О о о о о О О О о о ОиО О»

°в° °о °о0о0о0о0о0о0о°о°о00 ° °о° < о°лООл°лО о о о о о о о о °л о о о:

"" О О О ОлО ООО о ООиоО о ° оо О О О ° °ЛО » °0 О О Оо О О О( О о° о ооо О о оо.ооуо оо^оо о о ° о О 000 о о о°о!о0о0о°о0о °о0о00?0° о

шт шш

Рисунок 3. Кристаллизация в двухмодовой модели МКФП, расчеты для различной аппроксимации КЭ для уравнения (7)

Отчётливо видно, что линейная аппроксимация С1 при данном размере КЭ приводит не только к завышению скорости фронта (а)-(Ь), но и огрублению самой структуры (а). После окончания релаксации система не достигает равновесной квадратной решетки (е). При использовании квадратичных элементов

С2

качественно наблюдается каскад структурных переходов от полосчатой к треугольной (Ь), затем к сотовидной и наконец к квадратной решетке (£). После окончания релаксации при £ = 2000 наблюдается равновесная поликристаллическая квадратная фаза.

Как видно из сравнения панелей (а) и (Ь) рисунка 4, при величине I =1 различия в динамике распространения фронта для аппроксимации

Рисунок 4. Кристаллизация в двухмодовой модели МКФП, расчеты в единый момент времени £ = 50 для двух типов аппроксимации полиномами в КЭ

С1 и С2 по прежнему сохраняются, хотя в данном случае они уже не так велики как (а)-(Ь) рисунка 3. При этом для I = 1.5 при аппроксимации С1 весь каскад структурных превращений, аналогичный (Ь), (ё) и (£) рисунка 3, наблюдается в полной мере.

Равновесный параметр решетки Л, позволяющий определить сходимость результирующих решений, указан в таблице 1. На основе

Таблица 1. Зависимость равновесных параметров квадратной решетки А от порядка аппроксимации КЭ и его размера

Размер элемента 1 0.5 1 1.5 2

А из линейной С1 -аппроксимации А из квадратичной С2-аппроксимации 6.8 7.2 7.5 6.8 6.8 7.1 7

приведенных данных можно сделать вывод о достижении оптимального размера конечного элемента и сходимости двумерного решения на сетке: I ~ 0.5 и Л/1 = 13.6 для линейной С1 -аппроксимации КЭ и I ~ 1, Л/1 = 6.8 в случае квадратичной

С2

-аппроксимации. Аналогичные расчеты были проведены и для трехмерного случая в расчетной области Ьх = 120, Ьу = 50, Ьг = 50, здесь сходимость решений происходила

быстрее, оптимальный размер КЭ I к 1 с Л/1 = 7.2 для линейной С1-аппроксимации и, соответственно, I к 1.7, Л/1 = 4.23 для квадратичной С2 - аппроксимации.

При размере конечного элемента I =1 оценка скоростей фронтов для линейных элементов даёт скорость V = 2.08 безразмерных единиц, а для квадратичных V = 1.66. При I = 0.5 для линейных элементов V = 1.8, а для квадратичных V = 1.65. Таким образом, мы считаем, что сходимость решения по скорости фронта V для С1 достигается в районе I = 0.3..0.5, а для С2 при I = 0.5..1. Сходимость по равновесному параметру решетки Л достигается для линейных КЭ при I = 0.5..1, для квадратичных элементов при I = 1..1.5.

В таблице 2 приведено сравнение производительности (времени

Таблица 2. Время расчета двухмодовой модели в минутах, в скобках— количество степеней свободы

Размер элемента 1 Количество КЭ 0.5 1 1.5 2 203500 49712 22324 12480

Расчет линейной С^-аппроксимации Расчет квадратичной С2-аппроксимации 16 4 2 1 (511755) (125785) (56820) (31955) 151 24 7 3 (2041005) (500125) (225255) (126305)

расчета в минутах) в зависимости от количества конечных элементов и степеней свободны. Показано, что совокупное влияние уменьшения количества конечных элементов и увеличения сложности базисных функций на вычислительную сложность задачи линейно, при этом количество степеней свободы при повышении порядка аппроксимации КЭ растёт как 4N.

5. Сравнение эффективности расчетов по модели МКФП

Для анализа эффективности решения двухмодового уравнения МКФП (7) была подготовлена трехмерная тестовая задача Lx — Ly — Lz — 20, I — 1.7, т — 0.01, C2, количество степеней свободы DOF — 118000, количество потоков при решении было ограничено N — 8. При расчетах использовались стандартные параметры адаптивного временного шага.

При использовании ПО COMSOL Multiphysics на процессорах AMD была подключена библиотека BLAS (Basic Linear Algebra Subprograms) для оптимизации производительности [16].

Реализация в пакете Fenics (DOLFIN) была выполнена с использованием интерфейса C++, компилятор GCC 9.4.0, Ubuntu 20.04 [17].

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

• AMD Epyc 7302: 3.0 GHz, 16 яд., 32 пот., DDR4-3200;

• Intel i7-9700KF: 4.7 GHz, 8 яд., 8 пот., DDR4-2666;

• AMD TR-2970wx2: 4 GHz, 24 яд., 48 пот., DDR4-2933;

• Xeon X5650: 2.4 GHz, 6 яд., 12 пот., DDR3-1333.

Для всех систем количество памяти составляло более 64 Гб.

Сравнение производительности при решении задачи МКФП, уравнение (7) с использованием различных видов прямых решателей (MUMPS, PARDISO) для пакетов COMSOL Multiphysics и Fenics (DOLFIN) представлено на рисунке 5.

Рисунок 5. Сравнение производительности с использованием различных видов прямых решателей

В приведенных расчетах показано, что для уравнений МКФП наиболее оптимальным является использование решателей MUMPS, при этом лимитирующим фактором для подобных задач является частота отдельного ядра. Как показано в работе по исследованию эффективности при распараллеливании подобных задач [22], оптимальным является N = 8 потоков, так дальнейший рост их количества повышает издержки на распараллеливание. Для случая реализации расчётов на процессорах нами была выполнена проверка этой гипотезы. Результаты производительности при решении с помощью COMSOL MUMPS при различном количестве используемых потоков на процессоре AMD Epyc 7302 приведены на рисунке 6.

Рисунок 6. Производительность решения задачи МКФП, уравнение (7) с использованием СОМЯОЬ МиШрЬувкв

Показано, что заметное изменение скорости расчета происходит до N =8 — 10 потоков, что в целом согласуется с результатами, полученными в [22]. При дальнейшем увеличении количества потоков после N = 8 накладные расходы на распараллеливание становятся значительны, что обуславливает малый прирост производительности при дальнейшем увеличении N.

6. Выводы

В работе представлена двухмодовая модель кристаллического фазового поля (КФП), описывающая микроскопическую эволюцию и упорядочение вещества в процессе кристаллизации из разупорядоченной фазы. К данной модели, представленной нелинейным дифференциальным уравнением десятого порядка по пространству и второго по времени, применен метод конечных элементов. Показано, что в силу периодического вида численных решений и появления дополнительного масштаба структуры, точность аппроксимации решения значительно зависит от параметров дискретизации. Исследовано влияние размеров конечных элементов (КЭ) и порядка их дискретизации на аппроксимацию решения двухмодовой модели КФП для случая движения плоского фронта кристаллизации. Исследована сходимость решения на сетке для описания динамики движения фронта и формирования периодической структуры. Рассчитаны безразмерные критерии сходимости относительно периода Л описываемого периодического решения в двумерных задачах: для линейных КЭ - Л/1 = 13.6; для квадратичных - Л/1 = 6.8. Также рассчитаны критерии для трехмерной постановки задачи КФП: для линейных КЭ - Л/1 = 7.2; для квадратичных - Л/1 = 4.23.

Рассмотрено возникновение неустойчивости решения, связанной с кинетическими эффектами, обусловленными вкладом второй

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

Выполнено сравнение эффективности расчетов по модели КФП на различных аппаратных платформах с помощью программных пакетов COMSOL Multiphysics и Fenics (DOLFIN) для различных решателей и количества доступных потоков. Показано, что лимитирующим фактором является быстродействие отдельного ядра, а также преимущество решателя MUMPS независимо от программной реализации.

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

[1] Elder K.R., Katakowski M., Haataja M., Grant M. Modeling elasticity in crystal growth // Physical Review Letters.- 2002.- Vol. 88.- No. 24.-pp. 245701. ее 6-

[2] Provatas N., Elder K. R. Phase-Field Methods in Materials Science and Engineering.- Weinheim: Wiley-VCH Verlag GmbH & Co. KGaA.- 2010.-isbn 978-3-527-40747-7.- 312 pp. fee

[3] Starodumov I., Ankudinov V., Nizovtseva I. A review of continuous 'modeling of periodic pattern formation with modified phase-field crystal models // The European Physical Journal Special Topics.- 2022,- Vol. 231,- pp. 1135-1145. d 66

[4] Galenko P., Danilov D., Lebedev V. Phase-field-crystal and Swift-Hohenberg equations with fast dynamics // Physical Review E.- 2009.- Vol. 79.- No. 5.051110. 66 6-

[5] Elder K. R., Rossi G., Kanerva P., Sanches F., Ying S. C., Granato E., Achim C. V., Ala-Nissila T. Patterning of heteroepitaxial overlayers from nano to micron scales // Physical Review Letters.- 2012.- Vol. 108.- No. 22.-pp. 226102. 66

[6] Asadi E., Asle Zaeem M. A Review of Quantitative Phase-Field Crystal Modeling of Solid-Liquid Structures jj JOM - 2015,- Vol. 67,- No. 1-pp. 186-201. I '66

[7] Bueno J., Starodumov I., Gomez H., Galenko P., Alexandrov D. Three dimensional structures predicted by the modified phase field crystal equation // Computational Materials Science.- 2016.- Vol. 111.- pp. 310-312. 66 тз

[8] Стародумов И.О., Галенко П. К., Кропотин Н. В., Александров Д. В. Об аппроксимации периодического решения уравнения кристаллического фазового поля при расчетах методом конечных элементов // Программные

системы: теория и приложения.- 2018.- Т. 9.- №4(39).- с. 265-278. t66

71 73

[9] Зенкевич О. Метод конечных элементов в технике.- М.: Мир.- 1975.534 с. fee

[10] Ankudinov V. E., Galenko P. K., Kropotin N. V., Krivilyov M. D. Atomic density functional and diagram of structures in the phase field crystal model // Journal ol Experimental and Theoretical Physics.- 2016.- Vol. 122.- No. 2.-pp. 298-309. tee 67 71

[11] Jaatinen A., Achim C. V., Elder K. R. Thermodynamics of bcc metals in phase-field-crystal models // Physical Review E.- 2009.- Vol. 80.- No. 3.031602,- 10 pp. ее

[12] Emdadi A., Asle Zaeem M., Asadi E. Revisiting phase diagrams of two-mode phase-field crystal models // Computational Materials Science.- 2016.-Vol. 123,- pp. 139-147. 67

[13] Ankudinov V., Elder K. R., Galenko P. K. Traveling waves of the solidification and melting of cubic crystal lattices // Physical Review E.- 2020.- Vol. 102.-No. 6,- 062802. 67

[14] Stefanovic P., Haataja M., Provatas N. Phase-field crystals with elastic interactions // Physical Review Letters.- 2006,- Vol. 96,- No. 22,- 225504.

d j67

[15] Galenko P. K., Gomez H., Kropotin N. V., Elder K. R. Unconditionally stable method and numerical solution of the hyperbolic phase-field crystal equation // Physical Review E.- 2013,- Vol. 88,- No. 1,- 013310. 68

[16] COMSOL Multiphysics® v. 6.0, http://www.comsol.com.- Stockholm: COMSOL AB. fee, 7e

[17] Alnaes M., Blechta J., Hake J., Johansson A., Kehlet B., Logg A., Richardson C., Ring J., Rognes M.E., Wells G.N. The FEniCS Project Version 1.5 // Archive of Numerical Software.- 2015.- Vol. 3.- No. 100.-pp. 9-23. ' 68 76

[18] Balay S., Gropp W. D., Mclnnes L. C., Smith B. F. Efficient management of parallelism in object-oriented numerical software libraries // Modern Software Tools for Scientific Computing, eds. E. Arge, A. M. Bruaset, H. P. Langtangen.- Birkhauser Press.- 1997.- pp. 163-202. fee

[19] Mkhonta S. K., Elder K. R., Huang Z. F. Exploring the complex world of two-dimensional ordering with three modes // Physical Review Letters.- 2013.-Vol. 111.- No. 3,- 035501. d 70

[20] Galenko P. K., Elder K. R. Marginal stability analysis of the phase field crystal model in one spatial dimension // Physical Review B.- 2011.- Vol. 83.-No. 6,- pp. 064113. 7i 73

[21] Ankudinov V., Galenko P. K. Growth of different faces in a body centered cubic lattice: A case of the phase-field-crystal modeling // Journal of Crystal Growth.- 2020,- Vol. 539,- 125608. I ' 71

[22] Стародумов И.О., Павлюк Е. В., Абрамов С.М., Клюев Л. В., Га-ленко П. К., Александров Д. В. Эффективность распараллеливания алгоритма решения уравнения PFC с использованием библиотеки PetIGA // Вестн. Удмуртск. ун-та. Матем. Мех. Компьют. науки.- 2016.- Т. 26.-№ 3.- с. 445-450. d 78

Поступила в редакцию 28.04.2022;

одобрена после рецензирования 25.05.2022; принята к публикации 21.06.2022.

Рекомендовал к публикации к.т.н. С. А. Амелькин

Информация об авторах:

Владимир Евгеньевич Анкудинов К.ф.-м.н., научный сотрудник теоретического отдела ИФВД РАН, г. Москва (г. Троицк). Область научных интересов: тепломассоперенос, формирование упорядоченных структур и кристаллизация, фазовые превращения, метод фазового поля, численные методы.

0000-0001-8563-5862 e-mail: vladimir@ankudinov.org

Илья Олегович Стародумов К.ф.-м.н., научный сотрудник лаборатории многомасштабного математического моделирования УрФУ, г. Екатеринбург. Область научных интересов: структурно-фазовые первращения, кристаллизация, гидродинамика, численные методы.

0000-0001-6397-488X e-mail: ilya.starodumov@urfu.ru

Вклад авторов: В. Е. Анкудинов — 85% (постановка задачи, программирование, формальный анализ, проведение экспериментов, ресурсы, написание черновой версии, доработка и редактирование, визуализация); И. О. Стародумов —15% (методология исследований, верификация, доработка и редактирование).

Авторы заявляют об отсутствии конфликта интересов.

ISSN 2079-3316 PROGRAM SYSTEMS: THEORY AND APPLICATIONS vol. 13, No2(53), pp. 65-84 Research Article _ computational science

UDC 519.688:519.63

10.25209/2079-3316-2022-13-2-65-84;

Approximation of periodic solutions of two-mode phase-field crystal model

Vladimir Evgenyevich Ankudinov1 M, Ilya Olegovich Starodumov2

Vereshchagin Institute for High Pressure Physics of RAS, Moscow, Russia1, Ural

Federal University, Ekaterinburg, Russia2,

vladimir@ankudinov.orgM (learn more about the authors in Russian on p. 81)

Abstract. In present paper we consider a mathematical model of two-mode phase-field crystal (PFC). This model describes the microscopic evolution and ordering of matter during crystallization from the homogeneous phase. The model is represented by a nonlinear partial differential equation of the tenth order in space and second order in time. The solution of PFC-model was performed using the Galerkin finite-element method. Due to the periodic form of the numerical solutions of this model, the additional spatial scale appeared and so this requires an increased discretization accuracy. The mesh convergence criteria and discretization parameters for the numerical solutions is considered, taking into account the computational complexity of two-mode PFC-model. The influence of size of finite elements (FE) and their order of base functions on the approximation of the solution in FE is considered. The correspondent numerical solution is devoted to the motion of planar crystallization front. The optimal sizes of FEs are determined, and the efficiency of numerical simulations using various software packages and solvers is compared. (In Russian).

Key words and phrases: crystal phase field method, numerical calculations, finite elements, approximation

2020 Mathematics Subject Classification: 35Q35; 35Q68; 68N30

For citation: Ankudinov V. E., Starodumov I. O. Approximation of periodic solutions of two-mode phase-field crystal model // Program Systems: Theory and Applications, 2022, 13:2(53), pp. 65-84. (In Russian). http://psta. psiras.ru/read/psta2022_2_65-84.pdf

© Ankudinov V. E., Starodumov I. O.

C«feW.<a»0

References

[1] K. R. Elder, M. Katakowski, M. Haataja, M. Grant. "Modeling elasticity in crystal growth", Physical Review Letters, 88:24 (2002), pp. 245701. 66 67

[2] N. Provatas, K. R. Elder. Phase-Field Methods in Materials Science and Engineering, Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, 2010, isbn 978-3-527-40747-7, 312 pp. fee

[3] I. Starodumov, V. Ankudinov, I. Nizovtseva. "A review of continuous modeling of periodic pattern formation with modified phase-field crystal models", The European Physical Journal Special Topics, 231 (2022), pp. 1135-1145. fee

[4] P. Galenko, D. Danilov, V. Lebedev. "Phase-field-crystal and Swift-Hohenberg equations with fast dynamics", Physical Review E, 79:5 (2009), 051110. t66

67

[5] K. R. Elder, G. Rossi, P. Kanerva, F. Sanches, S.C. Ying, E. Granato, C. V. Achim, T. Ala-Nissila. "Patterning of heteroepitaxial overlayers from nano to micron scales", Physical Review Letters, 108:22 (2012), pp. 226102.

d 66

[6] E. Asadi, Asle Zaeem M.. "A Review of Quantitative Phase-Field Crystal Modeling of Solid-Liquid Structures", JOM, 67:1 (2015), pp. 186-201. tee

[7] J. Bueno, I. Starodumov, H. Gomez, P. Galenko, D. Alexandrov. "Three dimensional structures predicted by the modified phase field crystal equation", Computational Materials Science, 111 (2016), pp. 310-312. I 166 73

[8] I. O. Starodumov, P. K. Galenko, N. V. Kropotin, D. V. Aleksandrov. "On approximation of a periodic solution of the phase field crystal equation in simulations by the finite elements method", Program Systems: Theory and Applications, 9:4(39) (2018), pp. 265-278 (In Russian), d 66 ti 73

[9] O. Zenkevich. Finite element method in engineering, Mir, M., 1975, 534 pp. fee

[10] V. E. Ankudinov, P. K. Galenko, N. V. Kropotin, M. D. Krivilyov. "Atomic density functional and diagram of structures in the phase field crystal model", Journal of Experimental and Theoretical Physics, 122:2 (2016), pp. 298-309.

d 66 67 71

[11] A. Jaatinen, C.V. Achim, K. R. Elder. "Thermodynamics of bcc metals in

phase-field-crystal models", Physical Review E, 80:3 (2009), 031602, 10 pp.

66

[12] A. Emdadi, Asle Zaeem M., E. Asadi. "Revisiting phase diagrams of two-mode phase-field crystal models", Computational Materials Science, 123 (2016), pp. 139-147. ' 67

[13] V. Ankudinov, K.R. Elder, P. K. Galenko. "Traveling waves of the solidification and melting of cubic crystal lattices", Physical Review E, 102:6 (2020), 062802. I te7

[14] P. Stefanovic, M. Haataja, N. Provatas. "Phase-field crystals with elastic

interactions", Physical Review Letters, 96:22 (2006), 225504.

84

VladimirE. Ankudinqy, IlyaO. Starodumqv

[15] P. K. Galenko, H. Gomez, N. V. Kropotin, K. R. Elder. "Unconditionally stable method and numerical solution of the hyperbolic phase-field crystal

equation", Physical Review E, 88:1 (2013), 013310. I -fes

[16] COMSOL Multiphysics® v. 6.0, http://www.comsol.com, COMSOL AB, Stockholm. ^68, 76

[17] M. Alnaes, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg,

C. Richardson, J. Ring, M. E. Rognes, G.N. Wells. "The FEniCS Project Version 1.5", Archive of Numerical Software, 3:100 (2015), pp. 9-23. 68 re

[18] S. Balay, W. D. Gropp, L. C. Mclnnes, B.F. Smith. "Efficient management of parallelism in object-oriented numerical software libraries", Modern Software Tools for Scientific Computing, eds. E. Arge, A. M. Bruaset, H. P. Langtangen, Birkhauser Press, 1997, pp. 163-202. d es

[19] S. K. Mkhonta, K.R. Elder, Z. F. Huang. "Exploring the complex world of two-dimensional ordering with three modes", Physical Review Letters, 111:3 (2013), 035501. I t™

[20] P. K. Galenko, K. R. Elder. "Marginal stability analysis of the phase field crystal model in one spatial dimension", Physical Review B, 83:6 (2011), pp. 064113. 1 171 73

[21] V. Ankudinov, P. K. Galenko. "Growth of different faces in a body centered cubic lattice: A case of the phase-field-crystal modeling", Journal of Crystal Growth, 539 (2020), 125608. fi

[22] I. O. Starodumov, E. V. Pavlyuk, S. M. Abramov, L. V. Klyuev, P. K. Galenko,

D. V. Alexandrov. "The effectiveness of parallelizing an algorithm of the PFC equation solution using PetlGA library", Vestnik Udmurtskogo Universiteta: Matematika, Mekhanika, Komp'yuternye Nauki, 26:3 (2016), pp. 445-450 (In Russian). 77 78