CC BY
19
https://doi.Org/10.15350/17270529.2021.3.28
УДК 519.688:538.91
Градиентно-устойчивый алгоритм для динамики фазового поля при моделировании направленного затвердевания в Si-As
В. Г. Лебедев1,2, А. А. Обухов1, В. П. Бовин1, В. И. Ладьянов2
1 Удмуртский государственный университет, Россия, 426034, Ижевск, ул. Университетская, 1
2 Удмуртский федеральный исследовательский центр УрО РАН, Россия, 426067, Ижевск, ул. Т. Барамзиной, 34
Аннотация. Рассмотрены проблемы численного моделирования термодинамически согласованной модели динамики параметра порядка для фазовых превращений первого рода в рамках фазово-полевого подхода. Модель основана на расширенной необратимой термодинамике, использующей для описания состояния системы потенциалы Гиббса реальных систем. Численные алгоритмы для моделирования уравнений фазового поля как для локально-равновесных, так и локально-неравновесных моделей релаксационной динамики получены из условия градиентной устойчивости. Выполненные оценки позволяют построить градиентно-устойчивый алгоритм для моделирования уравнений фазового поля в реальных термодинамических системах. В работе проведён анализ градиентно-устойчивых численных схем на устойчивость по малым отклонениям для уравнений фазового поля в Фурье-приближении. Выполненные численные расчеты для раствора Si-As демонстрируют устойчивость полученного алгоритма как при увеличении шага по времени на несколько порядков, так и при росте движущих сил из-за сильного переохлаждения системы.
Ключевые слова: быстрые процессы затвердевания, локально-неравновесная модель фазового поля, потенциалы Гиббса, градиентно-устойчивые алгоритмы, численное моделирование.
И Владимир Лебедев, e-mail: lvg@Mdsu.ru
Gradient-Stable Algorithm for the Dynamics of the Phase Field at Modeling of Directional Solidification in Si-As
1,2 1 1 2 Vladimir G. Lebedev , , Alexander A. Obukhov , Vladimir P. Bovin , Vladimir I. Ladyanov
1 Udmurt State University (Universitetskaya Str., 1, Izhevsk, 426034, Izhevsk, Russian Federation)
2 Udmurt Federal Research Center UB RAS (T. Baramzina Str., 34, Izhevsk, 426067, Russian Federation)
Summary. The main provisions of the phase field method for concentrated solutions and problems of numerical simulation of obtained thermodynamically consistent equations have been considered briefly in the paper. The phase field is an order parameter for phase transformations of the first kind and takes certain values in the phase areas, and phase field quickly changes inside the diffuse interface. The model is based on the extended irreversible thermodynamics, which uses Gibbs potentials of real systems to describe state of the system. The complete statement of the problem includes directly the equation of the phase field and the diffusion equations in each of phases. The purpose of this study is construction of the numerical algorithm that allows to effectively simulate processes of the phase boundary movement. For the sake of simplicity at the first stage, the problem of the phase field modeling is separated from the diffusion tasks and takes into account only the driving force for interface migration, which is a large thermodynamic potential calculated with a constant distribution of concentrations. Numerical algorithms for modeling phase field equations for both locally equilibrium and locally-nonequilibrium models of relaxation dynamics are derived from the gradient stability condition, which guarantees decrease of free energy of the system during relaxation. The fulfilled estimates of eigenvalues of the Hessian matrix in the lattice space allow to construct a gradient-stable algorithm for modeling the time evolution of the phase field equation in real thermodynamic systems, despite existence of differences in the approximation of Gibbs energies of different materials. The paper analyzes gradient-stable numerical schemes for stability over small deviations for phase field equations under Fourier approximation. The numerical calculations performed for the Si-As solution demonstrate stability of the obtained algorithm both with an increase in the time step by several orders of magnitude, and with an increase in driving forces due to the severe undercooling of the system.
Keywords: fast solidification processes, locally non-equilibrium phase-field model, Gibbs potentials, gradient-stable algorithms, numerical modeling.
И Vladimir Lebedev, e-mail: lvg@udsu.ru
312
CHEMICAL PHYSICS AND MESOSCOPY, 2021, vol. 23, no. 3
ВВЕДЕНИЕ
Численное моделирование [1, 2] уравнений критической динамики для фазовых превращений в растворах является важным инструментом исследования процессов структурообразования реальных систем [3]. Одной из особенностей фазово-полевых моделей критической динамики является их нелинейность, необходимая для описания физической неустойчивости, связанной с зарождением и развитием микроструктур. Другая особенность -отсутствие локального равновесия, что приводит к сингулярному возмущению классических уравнений структурообразования в виде появления второй производной по времени с малым коэффициентом. Стандартные вычислительные алгоритмы для таких задач, как правило, неустойчивы [4], или в лучшем случае неэффективны [5]. К еще одной особенности можно отнести достаточно широкое разнообразие потенциалов Гиббса [6] для различных систем, несмотря на которое должно выполняться условие устойчивости численных алгоритмов.
Градиентно-устойчивые методы [5, 7, 8], гарантирующие монотонное понижение свободной энергии при моделировании процессов структурообразования с произвольным шагом по времени, представляют собой важное исключение среди других алгоритмов и получили в последнее время достаточно широкое распространение [9 - 12]. Обобщение градиентно-устойчивых методов на случай задач с отсутствием локального равновесия выполнено в работах [13, 14]. Появление термодинамически согласованных подходов для описания фазовых превращений в растворах [15 - 19] и повсеместное использование потенциалов Гиббса реальных систем [6, 20] приводит к необходимости дополнительных исследований, позволяющим соотнести ограничения, накладываемые градиентной устойчивостью, с устойчивостью по фон Нейману.
Целью данной работы является построение градиентно-устойчивого численного алгоритма для локально-неравновесной динамики фазового поля в рамках термодинамически согласованной модели затвердевания [19] и его анализ на устойчивость по фон Нейману.
УРАВНЕНИЯ ФАЗОВОГО ПОЛЯ ДЛЯ ЗАТВЕРДЕВАНИЯ БИНАРНОГО РАСТВОРА
Объектом исследования в данной работе являются численные алгоритмы для моделирования уравнения Аллена-Кана [21], описывающего процесс быстрой кристаллизации бинарного раствора в рамках идеологии фазового поля [1, 2]. Бинарный раствор - система, состоящая из растворителя и некоторого количества примеси с мольной концентрацией л при постоянной температуре T (по объему системы) и постоянном давлении P. Рассмотрим сосуществование двух различных термодинамически равновесных фаз а (а = L, S), ассоциированных с жидкой (L) и твердой (S) фазами. Мольные
концентрации примеси в каждой из фаз обозначим как ха. Пренебрегая внутренними напряжениями и, соответственно, давлением в системе, для описания термодинамически равновесных состояний используем обобщение объемной плотности фазовых потенциалов Гиббса Ga(xa, T) [22]:
G,, =Ъа Р“ (?)Ga (ха, T) + Wg (р) (1)
Предполагается, что температура T задана значительно ниже ликвидуса при заданной концентрации, так что твердая фаза при данной температуре устойчива, а жидкая -метастабильна. Начальное состояние соответствует сосуществанию жидкой и твердой фаз с границей, заданной фазовым полем.
Модельные функции
g(P) = р2(1 -р)2, Р(Р) = р2(3 - 2Р\ (2)
позволяют описывать термодинамические характеристики многофазных систем с помощью переменной фазового поля р, являющейся параметром порядка фазового превращения.
ХИМИЧЕСКАЯ ФИЗИКА И МЕЗОСКОПИЯ. 2021. Том 23, № 3
313
В общем случае, для каждой фазы а вводится свой параметр порядка - фазовое поле фа, такое, что '^фа = 1. Параметр порядка фа принимает значение 1 в фазе а, и значение 0 в остальных фазах. Для двухфазной системы достаточно одного поля фх= ф и дуальной переменной ф2 = 1 — ф. Функция Wg (ф) определяет потенциальный барьер между состояниями, выбор функции р(ф) связан с требованием монотонности на интервале [0,1] и граничными условиями р (0) = 0 и р (1) = 0, необходимыми для устойчивости термодинамических фаз [22].
Уравнения критической динамики термодинамических переменных и фазового поля могут быть получены [23] из условия убывания интегральной энергии Гиббса локальноравновесной системы вида (1). В рамках модели фазового поля для определённости примем, что жидкой фазе (L) соответствует ф = 0, а твёрдая (S) фаза опеределяется условием ф = 1. Считая для простоты, что концентрации фаз являются постоянными, остановимся более конкретно на уравнении Аллена-Кана для фазового поля, которое имеет вид:
Тфф + ф = Мф{е2У2ф — Wg' (ф) — тр ф)) (3)
где градиентное слагаемое описывает эффекты межфазной поверхности [24]. Величина АЛ (разность больших термодинамических потенциалов фаз) определена как:
АЛ = GS (xS ) — Gl (xL ) — Ах^арава( р а), (4)
где т - время релаксации скорости изменения фазового поля, Мр > 0 - коэффициент
мобильности фазового поля (кинетический коэффициент роста фазы), ва(ра) - функции Хэвисайда, которые равны единице при положительном значении аргумента и нулю в остальных случаях. В соотношениях (3) - (4) использованы следующие обозначения:
Ра (Ф)
дРа(ф) дф ’
g '(ф)
dg (ф) дф ’
ф =
дф
~dt'
ха(ф)
дха dt '
Уравнение (3) должно быть дополнено граничными условиями, которые также определяются из условия убывания интегральной энергии Гиббса локально-равновесной системы вида (1). Наиболее естественными условиями для фазового поля являются однородные условия Неймана (однородные граничные условия второго рода), поскольку они не навязывают системе дополнительного (кроме действия термодинамических движущих сил) поведения.
Используя стационарное одномерное решение уравнения (3) вида
ф0 = 2 (1 — tanh(x / 28)), (5)
где 8 - полуширина диффузной границы, выразим неизвестные параметры в уравнении (3) с помощью соотношений [22]:
- 6
S2 = 3о8, W = — (6)
08
где о - коэффициент поверхностного натяжения.
Для перехода к безразмерным переменным в уравнении (3), сделаем замену z ^ z8 и t ^ t82 / v, выразив мобильность фазового поля через параметр у [25, 26]:
М
ф
v
12о8
(7)
В этом случае уравнение для фазового поля в одномерном случае будет иметь вид
тфф + ф
ф’’ — 1 g ,ф) — АЛР(фХ
(8)
314
CHEMICAL PHYSICS AND MESOSCOPY, 2021, vol. 23, no. 3
в котором AQ = ADA /12а выступает в роли безразмерной движущей силы фазового перехода. Предполагая использование градиентно-устойчивых алгоритмов [11, 12] для моделирования динамики межфазной границы, запишем уравнение (8) для фазового поля в
одномерном случае как: Г Ф = р, (9)
\тФр + P = Ф" — ФФ,
где введены обозначения ФФ = - g , + A^p\ (10)
у” = д 2ф/ dz2, и Тф = тчу / 82 - безразмерное время релаксации фазового поля. Отметим, что
функция Ф в уравнении (10) является однородной частью уравнения, хотя выглядит, как
вклад источника. Эта функция, из-за нелинейности, как раз и является проблемой при численном моделировании.
ТЕОРЕМА ЭИРА О ГРАДИЕНТНОЙ УСТОЙЧИВОСТИ
Рассмотрим произвольную, дважды дифференцируемую функцию f (x) в евклидовом пространстве Е(и). Для f (x) можно показать (см. приложение в работе [14]), что:
f (x+y ^Е у
Е у
'удс. JХ+у
+
\дсг Jx
\dxi J
+ -4шп| y| ,
<я
<
2
(11)
где Ях, Яп - значения, ограничивающие спектр собственных значений матрицы вторых производных Mtj =д2 f (x)/ dxt dxj.
Перенося эти соотношения на функциональное пространство сеточных функций s, где компоненты s. определяют значения функции s в j -м узле сетки, и полагая, что вектора x, у определены как:
x = s (t + At), y = s (t) — s (t + At),
из (11) находим, что соответствущий функционал F[s ], понимаемый как интеграл (сумма) по области определения, будет удовлетворять условию
F[s(t + At )]< F[s(t)] (12)
Эйр показал [7, 8], что разделение энергетического функционала на две части:
F[?] = FC И + fe PL
где Fc [s ] - сжимающая (contractive), а FE[s] - растягивающая (expansive) часть полной энергии Г иббса, приводит к теореме (Эйра), согласно которой, при условии
ЯЯ < - Я
max ^ 1
где Я^ - максимальное собственное значение гессиана НЕ = д2FE /dstds. растягивающей части, а Яп - минимальное собственное значение гессиана для полного функционала свободной энергии Н =д2 F / ds. ds j, аппроксимация релаксационных уравнений
ХИМИЧЕСКАЯ ФИЗИКА И МЕЗОСКОПИЯ. 2021. Том 23, № 3
315
критической динамики для несохраняющегося параметра порядка s (t), к которым относится фазовое поле,
? SF[s ]
s =-----—
Ss
(13)
в виде
? л F И s + At—CLJ
U+At Ss
Н Л SFe Г?] = s -At—ELA
t+At
Ss
(14)
гарантирует невозрастание полной энергии Гиббса (12).
Под сжимающей частью функционала свободной энергии Fc [s] понимается выпуклая часть функционала F[s ]. Последнее означает, что ее гессиан неотрицательно определен (HC > 0). Соответственно, растягивающуя часть свободной энергии FE[s] определим как вклад с неположительным гессианом (HE < 0 ). Собственные значения Лс матрицы гессиана для выпуклой функции Fc [s ] строго неотрицательны, а собственные значения ЛЕ матрицы для вогнутой функции Fe [s ] строго неположительны для любых возможных конфигураций поля.
УСЛОВИЕ ЭИРА ДЛЯ РЕАЛЬНЫХ МАТЕРИАЛОВ
Задача (9) отличается от уравнений Эйра (13) в двух аспектах:
• В локально-неравновесных моделях появляется сингулярное возмущение локальноравновесных уравнений (14). Обобщение градиентной устойчивости на локальнонеравновесный случай было выполнено в работах [13, 14].
• Рассмотренные в работах [13, 14] примеры носят модельный характер и предполагают постоянство коэффициентов уравнений, в отличие от уравнения (8). Рассмотрим эти отличия более подробно.
Как показано в [13, 14], разностная схема для локально-неравновесной динамики несохраняющегося параметра порядка (9) (фазового поля) с шагом At по времени может быть записана в виде отображения с n -го момента времени на (n +1):
p{n+X) +
At2
At + 2г„
SFC [p] Sp
(n+1)
= p{n)
At
p v
2 S [pfi(n) 2 At г
At + 2r
p v
Sp
+ -
At + 2r
p n(n)
p ,
p
(n+1) = -p(n) +2 p
At
(p( n+1) -p{n))
где свободная энергия F[p] в нашем случае может быть представлена как
F [p] = J (j (V p)2 + | g (p) + P(p)AD)dV =Fc [p] + Fe [p].
(15)
(16)
p
Представленная свободная энергия должна быть разложена на сжимающую Fc [p] и растягивающую части F[p] в соответствие с идеологией работ [7, 8]. Выберем
Fe [p]=ZakFk[p], Fc [p]=Z (1 -ak)Fk[p] (17)
k k y ;
где вклады Fk [p] выбраны в виде:
Fp] = I\(Vp)2 dV, F2 [p] = 2 j(1 + 6AD)p2dV,
F3 [p] = -j (1 + 2AD)p3dV, F4 [p] = 2 J p4dV.
316
CHEMICAL PHYSICS AND MESOSCOPY, 2021, vol. 23, no. 3
Тогда оператор гессиана свободной энергии (17) определен второй вариацией по фазовому полю и равен:
И = 1
к
52
р2
Fk [И
V2 + (1 + 6AD) — 6(1 + 2AD )р + 6р2.
(18)
Поскольку собственное значение оператора Лапласа V2, с учетом однородных условий Неймана на границе, равно — /л1, нетрудно получить оценку собственных значений оператора H :
Л(И) = м2 + (1 + 6 AD) — 6(1 + 2 AD )р + 6р2.
(19)
Эта оценка достигает минимума при значениях м = 0 и р0 = 1 + AD.. Предполагая, что безразмерная движущая сила |AD| << 1, имеем:
3
Лтп (И) = (1 + 6AD) — -(1 + 2AD )2. (20)
В целях простоты, выберем а3 = 1, a = 1, что позволяет вычислять нелинейные вклады в растягивающую часть свободной энергии (19) явным образом. Тогда оценка максимального собственного значения гессиана для растягивающей части свободной энергии дает
ЛтхИ ) = «М* + «2(1 + 6 AD) — 6(1 + 2AD)p + 6р2. (21)
Для процесса затвердевания движущая сила отрицательна: AD < 0, поэтому
квадратичное выражение р2 — (1 + 2AD )р будет максимально при (р = 1, что в итоге дает
Лтах(ИE ) = «1 Мта* + «2(1 + 6 AD) — 12 AD.
Условие градиентной устойчивости требует
Лах( НЕ ) < 1 Лп( И)
(22)
(23)
или
«ML + «(1 + 6 AD) < 12 AD — 3AD2 — -.
(24)
Максимальное собственное значение на d -мерной разностной сетке с характерным пространственным шагом Az можно оценить как
М
ж2 d
max '“А2
>> 1,
что позволяет, выбрав а = —1 получить ограничение на а2 :
2 1 л
« (1 + 6AD) <^d + 12AD — 3AD2 —. 2 Az2 4
(25)
(26)
Поскольку правую часть (26), за счет соответствующего выбора Az, можно считать заведомо положительной, приравняем параметр « нулю, что переводит вклад F2 в неявную (сжимающую) часть схемы (18). Явно выписывая соотношения, в итоге получаем численный алгоритм нахождения р(и+1) :
ХИМИЧЕСКАЯ ФИЗИКА И МЕЗОСКОПИЯ. 2021. Том 23, № 3
317
P{n+X) +
At2
At + 2t,
-( 2V2p- (1 + 6AQ)p)
( n+1)
( \ At2 / о 0 \(n) 2AtT , .
= p(n) - t ^ (V2p - 3(1 + 2AQ)p2 + 2p3) + t J p{n),
At + 2 t
At + 2 t
(27)
p(n+1) =-p(n) +2 (p(n+1) -p(n)).
At v '
Отметим, что хотя вычислительный алгоритм (18) имеет первый порядок точности по времени, дальнейшее повышение точности алгоритма по времени не имеет особого смысла из-за нелинейности задачи, которая сама по себе ограничивает шаг по времени.
СПЕКТРАЛЬНЫЕ СВОЙСТВА ЧИСЛЕННЫХ СХЕМ
Градиентная устойчивочть не гарантирует устойчивости по малым отклонениям. Поэтому рассмотрим этот вопрос более подробно. Для анализа устойчивости полученных отображений по времени рассмотрим отображение малых отклонений от точного решения. Для этого будем полагать, что фазовое поле и концентрации в каждый момент t(n) представимы в виде точного решения (обозначенного индексом "0") и малого отклонения от него, которое представим в виде Фурье-компоненты:
P(n) = Po') + A(n) exp(iqz),
P
(n) =P0n + A)exp(iqz),
(28)
где q - волновой вектор возмущения.
Для амплитуд Фурье - компонент отклонений фазового поля от точного решения имеем
A
(n+1)
At2
At + 2t
-(2q2 +1 + 6AQ)
= A(n)
At2
At + 2 t
-(q2 + 6(1 + 2AQ)Pn) -6(Pn))2)
A(n+1) = -A(n) +-2 (A(n+1) - A(n)).
2AtTp At + 2
A(n),
(29)
Собственные значения отображения (29) приведены на рис. 1.
На рис. 1 (слева), видно, что длинноволновая часть собственного значения А+ для отображения (29) по величине несколько превосходит единицу, что формально соответствует неустойчивости численного расчета. Однако, физически, эта неустойчивость является неотъемлемой частью уравнения фазового поля, которое описывает переход из термодинамически метастабильного состояния в устойчивое состояние под действием движущей силы AD. Поэтому при численном расчете нужно ожидать малых флуктуаций в области чистых фаз и устойчивое решение на границе между ними. Второе собственное значение А (рис. 1, справа) всегда остается по модулю меньше единицы.
Общий вывод по последнему разделу, очевидно состоит в том, требование устойчивости по фон Нейману для малых отклонений от точного решения не дает адекватной проверки устойчивости численных алгоритмов, связанных с процессами возникновения структур. По сути, можно говорить, что построение таких алгоритмов должно быть основано на теореме Эйра и "a posteriori" проверяться на вычислительном эксперименте.
318
CHEMICAL PHYSICS AND MESOSCOPY, 2021, vol. 23, no. 3
Рис. 1. Зависимость собственных значений Л+ (слева) и Л (справа) для отображения (29)
в зависимости от значения фазового поля р0 и волнового числа q
Fig. 1. Dependence of the eigenvalues Л+ (left) and Л_ (right) for the mapping (29) depending on the value of the phase field and the wave number q
ЧИСЛЕННАЯ ПОСТАНОВКА ЗАДАЧИ
Для численного моделирования динамики фазового поля по алгоритму (27) рассмотрим одномерную постановку задачи, описывающую движение диффузной границы между фазами при затвердевании бинарного раствора на бесконечном интервале. Для простоты будем полагать, что концентрации примеси в фазах постоянны и не меняются со временем. Поэтому распределения концентраций примеси в фазах переохлажденного раствора считаем однородным x^ = C0 и x^ = C ■ Значение C0 выбирается произвольно (при расчете полагалось С0 = 0.1), а концентрация примеси в твердой фазе С определялась из условия равновесия на поверхности ликвидуса раствора при заданном C ■
В качестве начального условия для фазового поля используем стационарное решение (5). Материальные параметры модели взяты из работы [27] и приведены в таблице.
Таблица - Параметры моделирования для расплава Si-As
Table - Simulation parameters for Si-As melt
Параметр Parameter Величина The quantity Значение Meaning
Молярный объем Molar volume v m 1.2 -10 5 м3/моль [m3/mol]
Поверхностное натяжение Surface tension a 0.477 Дж/м2 [J/m2]
Мобильность р Mobility р Мр 8.777 м3/(Дж-с) [m3/(J-s)]
Кинетический коэффициент для р Kinetic coefficient for р V 1.22 -10 ~8 м2/с [m2/s]
Ширина диффузной границы Diffuse border width 5 1.875 -10 9 м [m]
Газовая постоянная Gas constant R 8.314 ДжАмоль-Ю [J/(mol-K)]
Время релаксации для р Relaxation time for р тр О 1 о сл
ХИМИЧЕСКАЯ ФИЗИКА И МЕЗОСКОПИЯ. 2021. Том 23, № 3
319
Для перехода к расчетам на конечном интервале, чтобы избежать искажений решения вблизи границы, правая граница области, к которой со временем приближается фронт кристаллизации, должна смещаться с течением времени. В силу однородности системы, при численном решении задачи удобнее сделать обратную операцию: решение транслируется в направлении, противоположном направлению движения фронта, когда фронт приближается к правой границе на определенное расстояние. Используемые граничные условия для задачи на конечном интервале имеют вид:
1,
ди
dz
z=L
0.
(30)
Для расчета методом конечных разностей использовалась пространственная область размером 100ё разбиваемая на N = 1000 интервалов (пространственных шагов). Данное разбиение выбрано как оптимальное по точности и затратам времени, что определялось из численных экспериментов с различными разбиениями.
Поскольку нам не удалось найти исследований устойчивости для алгоритмов моделирования фазового поля в литературе, для сравнения используем аналитическое исследовании устойчивости явной разностной схемы для уравнения локально-неравновесной диффузии, проводившееся в работе [28], где было получено достаточное условие устойчивости в виде At < qAx2 /8, где ё2 / DL Itd « 9.77, DL и zD - коэффициент диффузии
и диффузионное время релаксации в жидкости. В результате полученное в [28] ограничение может быть представлено в виде At < 1.22Ax2.
По недавней работе [29], где для моделирования полной изотермической задачи с перераспределением примеси и напряжениями на решетке 100 х100 х100 понадобилось около 1170 часов, можно утверждать, что оценка At «Ax2 для устойчивости явных схем останется адекватной оценкой. Поэтому для сравнения на рис. 2 показаны профили фазового поля при At = 10Ax2 и At = 107Ax2. Визуально графики практически идентичны, хотя численное сравнение показывает слабые отличия в нижней половине фронта. Графики равновесных потенциалов Гиббса и химпотенциалов для Si-As [27] приведены на рис. 3 при T = 1330 K.
Рис. 2. Профили фазового поля при различных шагах по времени при переохлаждении от температуры равновесия при данном составе на величину AT = 100 K
Fig. 2. Phase field profiles at different time steps under undercooling from the equilibrium temperature at a given composition by AT = 100 K
Рис. 3. Потенциалы Гиббса и химпотенциалы фаз при AT = 1330 K Fig. 3. Gibbs potentials and chemical potentials of phases at AT = 1330 K
320
CHEMICAL PHYSICS AND MESOSCOPY, 2021, vol. 23, no. 3
На рис. 4 показана зависимость скорости движения фронта от переохлаждения (показано, как меняется скорость межфазной границы в зависимости от понижения температуры на величину AT по отношению к температуре равновесия при заданных концентрациях раствора). Зависимость представлена в двойном логарифмическом масштабе. Она показывает слабую нелинейность от переохлаждения, вида v = AT1 045, связанную с нелинейность потенциалов Гиббса. Соответствующая гладкая зависимость демонстрирует, что устойчивость выполняется даже при достаточно больших движущих силах.
Рис. 4. Зависимость скорости движения фазовой границы при различных переохлаждениях
(концентрации фаз фиксированы)
Fig. 4. Dependence of the phase interface velocity at various undercoolings (concentrations of the phases are fixed)
ОБСУЖДЕНИЕ И ВЫВОДЫ
В работе рассмотрен вопрос построения градиентно-устойчивых схем для нелинейного уравнения фазового поля, описывающего диссипативную динамику движения границы раздела между твердой и жидкой фазами на основе равновесных потенциалов Гиббса. Для вывода численного алгоритма использована теорема Эйра, определяющая разбиение операторов уравнения на явную и неявную части в соответствие с их собственными значениями. В рассмотренной постановке (моделирование динамики фазового поля при постоянных концентрациях) задача о движении фронта оторвана от проблемы перераспределения примеси на межфазной границе, что было сделано с целью упрощения. Исследование численных алгоритмов для транспорта примеси на фоне заданного движения фронта, а также их совместной динамики является объектом дальнейших исследований.
Исследование устойчивости полученного алгоритма для численного решения уравнения фазового поля по малым отклонениям показывает наличие слабой неустойчивости при малых значениях волнового вектора. Последнее связано с явлением самоорганизации движущейся границы и не проявляется в виде неустойчивости расчетов при численном моделировании направленного процесса затвердевания Si-As, демонстрирующим устойчивость расчетов в широком диапазоне значений выбранного шага по времени. Представленный алгоритм может быть использован при разработке программных комплексов по прогнозированию микроструктуры материалов в неравновесных условиях.
ХИМИЧЕСКАЯ ФИЗИКА И МЕЗОСКОПИЯ. 2021. Том 23, № 3
321
Исследование выполнено при частичной финансовой поддержке РФФИ и Правительства Удмуртской Республики в рамках научного проекта № 18-42-180002.
The study was carried out with partial financial support from the Russian Foundation for Basic Research and the Government of the Udmurt Republic within the framework of scientific project No. 18-42-180002.
Авторы благодарят А.В. Обухова за ряд ценных замечаний по численному моделированию.
СПИСОК ЛИТЕРАТУРЫ
1. Singer-Loginova I., Singer H. M. The phase field technique for modeling multiphase materials // Reports on Progress in Physics, 2008, vol. 71, iss. 10, 106501. http://dx.doi.org/10.1088/0034-4885/71/10/106501
2. Emmerich H. Advances of and by phase-field modelling in condensed-matter physics // Advances in Physics, 2008, vol. 57, iss. 1. pp. 1-87. https://doi.org/10.1080/00018730701822522
3. Херлах Д., Галенко П., Холланд-Мориц Д. Метастабильные материалы из переохлажденных расплавов / пер. с англ. Н. В. Худойкиной, под ред. П. К. Галенко М.-Ижевск: НИЦ "РХД", ИКИ, 2010. 496 с.
4. Rogers T. M., Elder K. R., Desai R. C. Numerical study of the late stages of spinodal decomposition // Physical Review B, 1988, vol. 37, iss. 16, pp. 9638-9649. https://doi.org/10.1103/PhysRevB.37.9638
5. Vollmayer-Lee B. P., Rutenberg A. D. Fast and accurate coarsening simulation with an unconditionally stable time step // Physical Review E, 2003, vol. 68, iss. 6, 066703. https://doi.org/10.1103/physreve.68.066703
6. Computational Phase Diagram Database (CPDDB). URL: http://cpddb.nims.go.ip/cpddb/periodic.htm (дата обращения:25.03.2021).
7. Eyre D. J. An Unconditionally Stable One-Step Scheme for Gradient Systems, preprint. Published 1997. https://citeseerx. ist.psu. edu/viewdoc/download?doi= 10.1.1.38.313 &rep=rep 1 &type=pdf
8. Eyre D. J. Unconditionally gradient stable time marching in the Cahn-Hilliard equation // MRS Online Proceedings Library, 1998, vol. 529, pp. 39-46. https://doi.org/10.1557/PROC-529-39
9. Cheng M., Warren J. A. Controlling the accuracy of unconditionally stable algorithms in the Cahn-Hilliard equation // Physical Review E, 2007, vol. 75, iss. 1, 017702. https://doi.org/10.1103/PhysRevE.75.017702
10. Cheng M., Warren J. A. An efficient algorithm for solving the phase field crystal model // Journal of Computational Physics, 2008, vol. 227, iss. 12, pp. 6241-6248. https://doi.org/10.1016/i.icp.2008.03.012
11. Cheng M., Cottenier S., Emmerich H. Thermodynamic consistency and fast dynamics in phase field crystal modeling // Materials Science, 6 Aug 2009, 4 p. https://arxiv.org/abs/0807.3579
12. Choi J. W., Lee H. G., Jeong D., Kim J. An unconditionally gradient stable numerical method for solving the Allen-Cahn equation // Physica A: Statistical Mechanics and its Applications, 2009, vol. 388, iss. 9, pp. 1791-1803. https://doi.org/10.1016/i.physa.2009.01.026
13. Lebedev V., Sysoeva A., Galenko P. Unconditionally gradient-stable computational schemes
in problems of fast phase transitions // Physical Review E, 2011, vol. 83, iss. 2, 026705.
https://doi.org/10.1103/PhysRevE.83.026705
14. Галенко П. К., Лебедев В. Г., Сысоева А. А. Градиентная устойчивость численных алгоритмов в локально-неравновесных задачах критической динамики // Журнал вычислительной математики и математической физики. 2011. Т. 51, № 6. C. 1148-1165.
15. Steinbach I., Zhang L., Plapp M. Phase-field model with finite interface dissipation // Acta Materialia 2012, vol. 60, iss. 6-7, pp. 2689-2701. https://doi.org/10.1016/i.actamat.2012.01.035
16. Wang H., Galenko P. K., Zhang X., Kuang W., Liu F., Herlach D. M. Phase-field modeling of an abrupt disappearance of solute drag in rapid solidification // Acta Materialia, 2015, vol. 90, pp. 282-291. https://doi.org/10.1016/i.actamat.2015.02.021
17. Wang H., Kuang W., Zhang X., Liu F. A hyperbolic phase-field model for rapid solidification of a binary alloy // Journal of Materials Science, 2015, vol. 50, pp. 1277-1286. https://doi.org/10.1007/s10853-014-8686-1
18. Zhang X., Wang H., Kuang W., Zhang J.. Application of the thermodynamic extremal principle to phase-field modeling of non-equilibrium solidification in multi-component alloys // Acta Materialia, 2017, vol. 128, pp. 258-269. https://doi.org/10.1016/i.actamat.2017.02.026
19. Лебедев В. Г., Галенко П. К. Высокоскоростное затвердевание и плавление концентрированных
растворов и конструкция параллельных Хиллерта // Расплавы. 2016, № 5. С. 422-433.
https://doi.org/10.1134/S0036029516080097
20. The Microstructure Evolution Simulation Software. URL: http://web.micress.de/ (дата
обращения:10.03.2021).
21. Allen S. M., Cahn J. W. A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening // Acta Metallurgica, 1979, vol. 27, iss. 6, pp. 1085-1095. https://doi.org/10.1016/0001-6160(79)90196-2
322
CHEMICAL PHYSICS AND MESOSCOPY, 2021, vol. 23, no. 3
22. Kessler D. Sharp interface limits of a thermodynamially consistent solutal phase field model // Journal of Crystal Growth, 2001, vol. 224, iss. 1-2, pp. 175-186. https://doi.org/10.1016/S0022-0248(01)00814-4
23. Galenko P., Jou D. Diffuse-interface model for rapid phase transformations in nonequilibrium systems // Physical Review E, 2005, vol. 71, iss. 4, 046125. https://doi.org/10.1103/PhysRevE.71.046125
24. Cahn J. W., Hilliard J. E. Free Energy of a Nonuniform System. I. Interfacial Free Energy // The Journal of Chemical Physics, 1958, vol. 28, iss. 2, 258. https://doi.org/10.1063/1.1744102
25. Lebedev V. G., Abramova E. V., Danilov D. A., Galenko P. K. Phase-field modeling of solute trapping: comparative analysis of parabolic and hyperbolic models // International Journal of Materials Research, 2010, vol. 101, iss. 4, pp. 473-479. https://doi.ora/10.3139/146.110297
26. Galenko P. K., Abramova E. V., Jou D., Danilov D. A., Lebedev V. G., Herlach D. M. Solute trapping in rapid solidification of a binary dilute system: A phase-field study // Physical Review E, 2011, vol. 84, iss. 4, 041143. https://doi.ore/10.1103/PhysRevE.84.041143
27. Li S., Zhang J., Wu P. A comparative study on migration of a planar interface during solidification of nondilute alloys // Journal of Crystal Growth, 2010, vol. 312, iss. 7, pp. 982-988. https://doi.ore/10.1016/i.icrvsgro.2009.12.070
28. Galenko P. K., Krivilyov M. D. Model for isothermal pattern formation of growing crystals in undercooled binary alloys // Modelling and Simulation in Materials Science and Engineering, 2000, vol. 8, pp. 67-79. https://doi.ore/10.1088/0965-0393/8/1/306
29. Bottger B., Apel M, Budnitzki M., Eiken J., Laschet G., Zhou B. Calphad coupled phase-field model with mechano-chemical contributions and its application to rafting of y’ in CMSX-4 // Computational Materials Science, 2020, vol. 184, 109909. https://doi.org/10.1016/i.commatsci.2020.109909
REFERENCES
1. Singer-Loginova I., Singer H. M. The phase field technique for modeling multiphase materials. Reports on Progress in Physics, 2008, vol. 71, iss. 10, 106501. http://dx.doi.org/10.1088/0034-4885/71/10/106501
2. Emmerich H. Advances of and by phase-field modelling in condensed-matter physics. Advances in Physics, 2008, vol. 57, iss. 1. pp. 1-87. https://doi.org/10.1080/00018730701822522
3. Herlach D., Galenko P., Holland-Moritz D. Metastable solids from undercooled Melts. Amsterdam: Elsevier, 2007; 432 p.
4. Rogers T. M., Elder K. R., Desai R. C. Numerical study of the late stages of spinodal decomposition. Physical Review B, 1988, vol. 37, iss. 16, pp. 9638-9649. https://doi.org/10.1103/PhysRevB.37.9638
5. Vollmayer-Lee B. P., Rutenberg A. D. Fast and accurate coarsening simulation with an unconditionally stable time step. Physical Review E, 2003, vol. 68, iss. 6, 066703. https://doi.org/10.1103/physreve.68.066703
6. Computational Phase Diagram Database (CPDDB). URL: http://cpddb.nims.go.jp/cpddb/periodic.htm (accessed March 25, 2021).
7. Eyre D. J. An Unconditionally Stable One-Step Scheme for Gradient Systems, preprint. Published 1997. https ://citeseerx.ist.psu. edu/viewdoc/download?doi= 10.1.1.38.313 &rep=rep 1 &type=pdf
8. Eyre D. J. Unconditionally gradient stable time marching in the Cahn-Hilliard equation. MRS Online Proceedings Library, 1998, vol. 529, pp. 39-46. https://doi.org/10.1557/PROC-529-39
9. Cheng M., Warren J. A. Controlling the accuracy of unconditionally stable algorithms in the Cahn-Hilliard equation. Physical Review E, 2007, vol. 75, iss. 1, 017702. https://doi.org/10.1103/PhysRevE.75.017702
10. Cheng M., Warren J. A. An efficient algorithm for solving the phase field crystal model. Journal of Computational Physics, 2008, vol. 227, iss. 12, pp. 6241-6248. https://doi.org/10.1016/i.icp.2008.03.012
11. Cheng M., Cottenier S., Emmerich H. Thermodynamic consistency and fast dynamics in phase field crystal modeling. Materials Science, 6 Aug 2009, 4 p. https://arxiv.org/abs/0807.3579
12. Choi J. W., Lee H. G., Jeong D., Kim J. An unconditionally gradient stable numeric al method for solving the Allen-Cahn equation. Physica A: Statistical Mechanics and its Applications, 2009, vol. 388, iss. 9, pp. 1791-1803. https://doi.org/10.1016/i.physa.2009.01.026
13. Lebedev V., Sysoeva A., Galenko P. Unconditionally gradient-stable computational schemes in problems of fast phase transitions. Physical Review E, 2011, vol. 83, iss. 2, 026705. https://doi.org/10.1103/PhysRevE.83.026705
14. Galenko P. K., Lebedev V. G., Sysoeva A. A. Gradient stability of numerical algorithms in local nonequilibrium problems of critical dynamics. Computational Mathematics and Mathematical Physics, 2011, vol. 51, no. 6, pp. 1074-1090. https://doi.org/10.1134/S0965542511060078 (In Russian).
15. Steinbach I., Zhang L., Plapp M. Phase-field model with finite interface dissipation. Acta Materialia, 2012, vol. 60, iss. 6-7, pp. 2689-2701. https://doi.org/10.1016/i.actamat.2012.01.035
16. Wang H. F., Galenko P. K., Zhang X., Kuang W. W., Liu F., Herlach D. M. Phase-field modeling of an abrupt disappearance of solute drag in rapid solidification. Acta Materialia, 2015, vol. 90, pp. 282-291. https://doi.org/10.1016/i.actamat.2015.02.021
17. Wang H. F., Kuang W. W., Xiao Z., Liu F. A hyperbolic phase-field model for rapid solidification of a binary alloy. Journal of Materials Science, 2015, vol. 50, pp. 1277-1286. https://doi.org/10.1007/s10853-014-8686-1
ХИМИЧЕСКАЯ ФИЗИКА И МЕЗОСКОПИЯ. 2021. Том 23, № 3
323
18. Zhang X., Wang H., Kuang W., Zhang J.. Application of the thermodynamic extremal principle to phase-field modeling of non-equilibrium solidification in multi-component alloys. Acta Materialia, 2017, vol. 128, pp. 258-269. https://doi.org/10.1016/j.actamat.2017.02.026
19. Lebedev V. G., Galenko P. K. High-rate solidification and melting of concentrated solutions and the
Hillert parallel construction. Russian Metallurgy (Metally), 2016, vol. 8, pp. 785-792. (In Russian).
https://doi.org/10.1134/S0036029516080097
20. The Microstructure Evolution Simulation Software. URL: http://web.micress.de/ (accessed March 10, 2021).
21. Allen S. M., Cahn J. W. A microscopic theory for antiphase boundary motion and its application to antiphase
domain coarsening. Acta Metallurgica, 1979, vol. 27, iss. 6, pp. 1085-1095. https://doi.org/10.1016/0001-
6160(79)90196-2
22. Kessler D. Sharp interface limits of a thermodynamially consistent solutal phase field model. Journal of Crystal Growth, 2001, vol. 224, iss. 1-2, pp. 175-186. https://doi.org/10.1016/S0022-0248(01)00814-4
23. Galenko P., Jou D. Diffuse-interface model for rapid phase transformations in nonequilibrium systems. Physical Review E, 2005, vol. 71, iss. 4, 046125. https://doi.org/10.1103/PhysRevE.71.046125
24. Cahn J. W., Hilliard J. E. Free Energy of a Nonuniform System. I. Interfacial Free Energy. The Journal of Chemical Physics, 1958, vol. 28, iss. 2, 258. https://doi.org/10.1063/1.1744102
25. Lebedev V. G., Abramova E. V., Danilov D. A., Galenko P. K. Phase-field modeling of solute trapping: comparative analysis of parabolic and hyperbolic models. International Journal of Materials Research, 2010, vol. 101, iss. 4, pp. 473-479. https://doi.org/10.3139/146.110297
26. Galenko P. K., Abramova E. V., Jou D., Danilov D. A., Lebedev V. G., Herlach D. M. Solute trapping in rapid solidification of a binary dilute system: A phase-field study. Physical Review E, 2011, vol. 84, iss. 4, 041143. https://doi.org/10.1103/PhysRevE.84.041143
27. Li S., Zhang J., Wu P. A comparative study on migration of a planar interface during solidification of nondilute alloys. Journal of Crystal Growth, 2010, vol. 312, iss. 7, pp. 982-988.
https://doi.org/10.1016/i.icrvsgro.2009.12.070
28. Galenko P. K., Krivilyov M. D. Model for isothermal pattern formation of growing crystals in undercooled binary alloys. Modelling and Simulation in Materials Science and Engineering, 2000, vol. 8, pp. 67-79. https://doi.org/10.1088/0965-0393/8/1/306
29. Bottger B., Apel M, Budnitzki M., Eiken J., Laschet G., Zhou B. Calphad coupled phase-field model with mechano-chemical contributions and its application to rafting of f in CMSX-4. Computational Materials Science, 2020, vol. 184, 109909. https://doi.org/10.1016/i.commatsci.2020.109909
Поступила 27.07.2021; принята к опубликованию 23.09.2021 Received 27 July 2021; accepted 23 September 2021
Лебедев Владимир Геннадьевич, кандидат физикоматематических наук, заведующий кафедрой теоретической физики ИМИТИФ УдГУ, старший научный сотрудник НЦ МФМ УдмФИЦ УрО РАН, Ижевск, Российская Федерация, e-mail: lvs@udsu.ru
Обухов Александр Андреевич, младший научный сотрудник НЦ МФМ УдмФИЦ УрО РАН, Ижевск, Российская Федерация
Бовин Владимир Павлович, кандидат физикоматематических наук, доцент кафедры теоретической физики ИМИТИФ УдГУ, Ижевск, Российская Федерация
Ладьянов Владимир Иванович, доктор физикоматематических наук, руководитель Научного центра МФМ УдмФИЦ УрО РАН, Ижевск, Российская Федерация
Vladimir G. Lebedev, Cand. Sci. (Phys.-Math.), Head of the Department of Theoretical Physics, IMITIF Udmurt State University, Senior Researcher, Scientific Center of the MPM, Udmurt Federal Research Center UB RAS, Izhevsk, Russian Federation, e-mail: lvg@udsu.ru
Alexander A. Obukhov, Junior Researcher, Scientific Center of the MPM, Udmurt Federal Research Center UB RAS, Izhevsk, Russian Federation
Vladimir P. Bovin, Cand. Sci. (Phys.-Math.), Associate Professor of the Department of Theoretical Physics, IMITIF Udmurt State University, Izhevsk, Russian Federation
Vladimir I. Ladyanov, Dr. Sci. (Phys.-Math.), Head of the Scientific Center of the MPM Udmurt Federal Research Center UB RAS, Izhevsk, Russian Federation
324
CHEMICAL PHYSICS AND MESOSCOPY, 2021, vol. 23, no. 3
