Об одном методе решения задачи кристаллизации многокомпонентного раствора в цилиндрической ампуле Текст научной статьи по специальности «Физика»

УДК 519.6

DOI: 10.18698/1812-3368-2017-4-118-138

ОБ ОДНОМ МЕТОДЕ РЕШЕНИЯ ЗАДАЧИ КРИСТАЛЛИЗАЦИИ МНОГОКОМПОНЕНТНОГО РАСТВОРА В ЦИЛИНДРИЧЕСКОЙ АМПУЛЕ

О.В. Щерица1' 2 [email protected]

А.О. Гусев2 [email protected]

О.С. Мажорова1' 2' 3 [email protected]

1 Институт прикладной математики им. М.В. Келдыша РАН, Москва, Российская Федерация

2 МГТУ им. Н.Э. Баумана, Москва, Российская Федерация

3 МГУ им. М.В. Ломоносова, Москва, Российская Федерация

Аннотация Ключевые слова

Предложен численный метод решения термодиффузион- Задача Стефана, фазовый пере-ной задачи Стефана. Рассмотрена самосогласованная ход, математическое моделиро-модель квазиравновесного процесса кристаллизации вание трехкомпонентного раствора в тонкой цилиндрической ампуле. Состояние системы описывается средней по радиусу температурой и концентрациями. Соответствующая одномерная задача о фазовом переходе в многокомпонентной системе аппроксимирована неявной консервативной разностной схемой и решена совместно. Разработанный алгоритм относится к классу методов, в явном виде отслеживающих положение фронта кристаллизации. Предложенный метод использован для численного моделирования процесса кристаллизации трехком- Поступила в редакцию 24.10.2016 понентного соединения © МГТУ им. Н.Э. Баумана, 2017

Работа выполнена при поддержке гранта РФФИ № 15-01-03741 и программы фундаментальных исследований Президиума РАН П-5

Введение. Численные методы решения задач типа Стефана развиваются преимущественно в двух направлениях. Первое составляют алгоритмы сквозного счета, в которых граница фазового перехода явно не выделяется [1, 2], второе — методы, в явном виде отслеживающие положение фронта кристаллизации [2-7]. Последние более эффективны для решения задач о фазовом переходе в многокомпонентном растворе. Этот класс задач отличается от классической задачи Стефана тем, что температура фазового перехода зависит от состава твердой и жидкой фазы. Контроль за положением границы раздела фаз осуществляется либо за счет использования подвижных сеток, согласованных в исходных переменных с формой фронта кристаллизации [3-6, 8, 9], либо с помощью динамической замены переменных [1, 5-7, 10], которая выбирается так, чтобы в новых координатах расчетная область была регулярной, с фиксированными границами, совпадающими с координатными линиями — метод выпрямления фронта [11].

Задача Стефана относится к классу задач с подвижной внутренней границей. Фазовый переход определяется комплексом физических процессов: динамикой изменения температуры; динамикой изменения концентрации каждого компонента соединения; движением границы раздела фаз. На фронте кристаллизации температура и концентрации удовлетворяют балансным соотношениям, вытекающим из законов сохранения энергии и массы, и фазовой диаграмме системы: как правило, достаточно сложным нелинейным соотношениям, связывающим равновесные концентрации всех компонентов и температуру. Опыт решения подобных задач [5-7, 12, 13] свидетельствует о том, что использование алгоритмов, в основе которых лежит расщепление по физическим процессам, может приводить к искажению физического явления: например, по результатам расчета в системе должно происходить растворение, а в физическом процессе происходит рост, и наоборот [14]. В связи с этим наиболее надежными методами решения термодиффузинной задачи Стефана являются алгоритмы, в которых система уравнений, описывающая процессы тепломассопереноса и движение фронта кристаллизации, решается совместно [5, 6, 13].

Численное исследование кристаллизации трехкомпонентного соединения было проведено в работах [5, 15, 16]. Используемая математическая модель учитывала диффузию в твердой и жидкой фазах и движение фронта кристаллизации, при этом распределение температуры в системе предполагалось однородным по пространству, но изменяющимся во времени. Однако существует класс задач, в которых изменение температуры на границе раздела фаз играет важную роль [6]. В настоящей работе рассмотрен процесс кристаллизации трехкомпонентного соединения в тонкой ампуле. Математическая модель процесса учитывает диффузионный тепломассоперенос в твердой и жидкой фазах, движение фронта кристаллизации, выделение теплоты в результате фазового перехода, а также теплообмен с окружающей средой. Состояние системы можно описать средней в сечении температурой и концентрациями. В результате осесиммет-ричная задача сводится к одномерной, в которой распределение температуры описывается уравнением температуропроводности с источником, возникающим вследствие учета контактного теплообмена с окружающей средой, а распределение массы — уравнениями диффузии. Для решения соответствующей одномерной нестационарной задачи построен чисто неявный совместный алгоритм. Учет движения фронта кристаллизации осуществляется с помощью замены переменных типа Ландау, в результате которой граница раздела фаз остается неподвижной на протяжении всего процесса. В расчетах использована неявная консервативная разностная схема, предложенная в работе [6], модифицированная для учета источника теплоты в уравнении теплопроводности. Соответствующая система нелинейных сеточных уравнений решается методом Ньютона относительно вектора неизвестных, компонентами которого являются все искомые функции: концентрации в жидкой и твердой фазах; поле температуры; скорость движения фронта кристаллизации. Тем самым реализуется одновре-

менное определение положения границы раздела фаз, распределения температуры и массы в системе. Предложенный метод применяется для численного исследования процессов роста и растворения соединений А х БХ_ХС.

Постановка задачи. Рассмотрим задачу о кристаллизации тройного соединения Ах Бх_х С. Раствор компонентов А и Б в расплаве С приводится в контакт с затравкой, расположенной на дне тонкой цилиндрической ампулы (рис. 1). Затем температура системы понижается, жидкая фаза становится пересыщенной, и начинается процесс кристаллизации. Математическое моделирование проводится в рамках следующих предположений. На фронте кристаллизации состав жидкой и твердой фазы находится в квазиравновесии. Кинетические процессы на поверхности протекают достаточно быстро и не нарушают квазиравновесия. Фронт кристаллизации плоский. Конвекция в жидкой фазе незначительна, тепломассопере-нос в системе осуществляется механизмом диффузии.

А„В

А, В

Т*

( > ( у -«(J-.-T) ^ Г > /

* ') к 5 / ] < L

Твердая фаза Жидкая фаза

Рис. 1. Цилиндрическая ампула

Подобласть О/ = (£, Ь] х (0, Л], где Ь — длина ампулы; = £(£) — положение фронта кристаллизации; Я — радиус ампулы (см. рис. 1), занимает жидкая фаза, которая представляет собой раствор компонентов А и Б в расплаве С. Состав жидкой фазы определяется концентрациями компонентов А и Б, а концентрация вещества С вычисляется из условия

l(A) , VZ(B) , 1 (C) _

-х

X

= 1,

где х1 (^), ] = А, Б, С — мольная доля соответствующего компонента в жидкой фазе. Твердая фаза занимает область = [0, £) х (0, Я] и является твердым раствором Ах Б1_х С. Состав твердой фазы определяется одной переменной х, при этом мольные доли каждого ее компонента х5(] = АС, БС, удовлетворяют соотношениям

х5 (AC) = X,

,s (BC) _

= 1 - X.

(1)

Распределение температуры и состава в системе описывается уравнениями

dT 1 d ( ,дТ

d (,S дТ

cpр5— =--1 rkS— 1+—I kS—

dt r dr v дг ) dz v dz

(z, r) eQs;

(2)

clp Pl

T ~öt

l i.( k dT

r dr V dr

—{k-dT

dz V dz

(z, r)eQ/.

(3)

6C

s(A)

1 6

6t r бг

6Cl( j) 1 б

rD

s(A)

6C

s(A) Л

бг

(

6z

D

s(A)

6C

s(A) Л

6z

6t

r 6r

rDl( j)

6Cl(j) Л 6

6r

6z

Dl( j)

6cl( jj) Л

6z

, (z, r)eQs;

(z, r)eQl.

(4)

(5)

Здесь г, г — радиальная и вертикальная координаты цилиндрической системы координат; 1 — время; С1(1) — объемная концентрация 1-го компонента в твердой и жидкой фазе; ^^ — коэффициент диффузии 1-го компонента в твердой и жидкой фазе ) = А, В; Т — температура; Ср — удельная теплоемкость; к1 — коэффициент теплопроводности; р1 — плотность; 1 = 5 при г < % (твердая фаза) и 1 = I при г > % (жидкая фаза). Содержание компонента В в твердой фазе определяется условием (1), записанным в объемных единицах

(aa + flo)Cs(A) + (ав + ac)Cs(B)= ps,

(6)

где аА, аВ, ас — атомные веса соответствующих компонентов.

На границе раздела фаз (г = )) выполнены приведенные ниже условия. 1. Условие фазового равновесия: уравнения, описывающие связь концентраций компонентов А и В в твердой и жидкой фазах

Fm(Cs(A), Cs(B), Cl(A), Cl(B), T) = 0, m = 1, 2, 3.

(7)

Конкретный вид уравнений фазовой диаграммы определяется свойствами соединения А х В1_х С.

2. Закон сохранения энергии (условие Стефана):

k-6T

6z

z = ВД-

- f 6T

6z

(8)

z = ВД+

Здесь X — скрытая теплота плавления; %% — скорость движения фронта кристаллизации.

3. Закон сохранения массы:

Ds( j)

6Cs( j)

6z

- Dl( j)

6Cl( j)

z = ВД-

6z

= _(С5( 1) _ С1( 1%, ) = А, В. (9)

г = %(!)+

Стенки цилиндра непроницаемы, поэтому граничные условия для поля концентрации имеют вид

DKj) Cj)

6n

= 0, j = A, B, у = s, l,

(10)

6Q

где 6Q — граница области Q = Qs Ц| Ql.

На поверхности цилиндра происходит теплообмен по закону Ньютона

6T

k у-

6n

6Q

= ((T* - Tl6Q).

(11)

Здесь а — коэффициент контактного теплообмена; Т* — температура окружающей среды.

Состояние системы можно описать функциями Т(г,1) и С^^(г, 1), ] = А, В, 1 = 5,I, где

1) Т (г ,1) — средняя температура в сечении г в момент времени 1:

2

R

T(z, t) = —fT(r, z, t)rdr;

R2 0

(12)

;у( j),

2) С (г, 1) — средняя концентрация компонента ) (] = А, В, 1 = 5,1) в сечении г в момент времени 1:

~1()) 2 Я

С (г,1) = -2- Г С1( %, г, 1)тйт. Я 2 0

Для того чтобы получить уравнение для температуры Т (г, I), проинтегрируем (2), (3) по радиусу г:

- Rdpy—rdr = -2- R1— Гkуr—1 rdr+-2- R — f k'1 — 1 rdr, у =

0 P 6t R 2j0 r 6r ^ 6r J R 2j0 6z l 6z J

с учетом (12) получим

R 20

c Ppy

T=А

6t R2

kу r

_6T 6r

- kу r

r=R

_6T 6r

> Я f

r=0 J

6z

kу

"6z

, у = s, /.

(13)

Поскольку значение температуры на боковой поверхности цилиндра не зависит от его радиуса, верно соотношение

2 Я

— Г Т(Я, г, 1)г^г = Т(Я, г, 1) = Т(г, 1), Я 0

что позволяет переписать (13) с учетом граничных условий на боковой поверхности (11) в виде

у у6Т 6

cp p у — = —

p 6t 6z

f

f 6T

6z

+(T* - T), у = s, l. R

(14)

Граничные условия на торцах ампулы

k у —

6z

z=0

* ~ 6T = -а(T* -Т), kу —

6z

= а (T* -т), у = s, l.

(15)

z=L

Аналогично получим уравнение для средней концентрации С ^ ^(г ,1), проинтегрировав (4), (5) по радиусу г. С учетом граничных условий запишем

8C

у( j)

8t

8z

Dу

8C

у(j) Л

8z

у = 5, l, j = A,B.

(16)

Граничные условия (10) примут вид

:у( j)

D у

8C

8z

= 0, у = 5, l, j = A, B.

(17)

z=0, L

На границе раздела фаз z = ^(t) выполнены условия Стефана:

= ^Р5 Ь;

8z

z = S(t)-

- ki8L

8z

z = 5(t) +

D5 (j)

8C

5 ( j)

8z

- Dl( j)

8С

i (j)

8z

z=ад-

~5( j) С (j) ?

= -(C - С ,

j = A, B.

(18)

(19)

z = §(t) +

В результате осесимметричная задача (2)-(11) сведена к одномерной задаче (14)-(19). Контактный теплообмен на боковой поверхности ампулы теперь учитывается в уравнении температуропроводности в виде источника теплоты, зависящего от температуры.

Замена переменных. Для решения задачи (14)—(19) с внутренней подвижной границей, положение которой необходимо определять в ходе решения задачи, применим метод выпрямления фронта. Для этого выполним замену независимых переменных (1, г) так, чтобы в новой системе координат , у) положение границы раздела фаз было фиксировано и совпадало с точкой у = 1, а границы области, точки г = 0 и г = Ь, перешли в точки у = 0 и у = 2 [11]. Связь между системами координат задается равенствами

У;(0, у е [0,1];

t = t, z = ф(у, t) = <

S(t) + (L-!#))(у -1), у e [1,2].

(20)

Якобиан преобразования (20) равен ] = д(у *) =1, где I = Щ) — длина со-

д(г, *) I

ответствующей фазы в исходной системе координат: в твердой фазе I = 15 = ), в жидкой фазе I = I1 = Ь ). В новых переменных области Оу ={0 < у <1} и Оу ={1< у < 2} заняты твердой и жидкой фазами, Оу = Оу и Оу.

В новой системе координат задача (14)-(19) записывается в виде [6]

8 (lу T У8 (ФtС)

8t 8у

f

иу 8T

\

8у { lу 8у

+ l(ayT+ßy), у e (0,1)u(1,2), у = 5, l. (21)

Здесь ку = ку / р у — кусочно-постоянная функция, принимающая значения к при 7 е[0,1) и к1 при 7 е (1,2]; ау = 2< /( руЯ); 3у =2(3 /( руЯ). На границе раздела фаз (у = 1) условие Стефана примет вид

жs дГ_

~ ну

у=1-

Ж- dT Т" "dy

= х

у=1+

d; dt

(22)

Уравнения для концентраций в переменных (t, у):

— (iу с1 (A) dt \

Ы (*с

У (A)

_d_

dy

f d^ dC!^ ^

ly dy

у e (0,1) u (1,2), у = s, l; (23)

,_ S(A) _ s(B)

(üa + flc) C + (йв + ac) C = ps

dt

(Ti K')

_d_

dy

' D«B) dCl(B) ^

dy

У e (0,1); , у e (1,2).

(24)

(25)

На границе раздела фаз имеем

Ds (j) dC

S ( j)

ls

dy

Dl(j) dC

i (j)

у=1-

= -(_s(j) - С (j))^, j = A, B. dt J

у=1+

ll dy

Pm(Cs(A), Cs(B), Cl(A), Cl(B),T)| =0, m = 1,2,3.

1у=1

Граничные условия (11), (17) в переменных (у, t):

AI dT

ls dy

у=0

=-а T -T

I у=0

Ai dT

ll dy

у =2

= а T - T

ly=2

Ds(A) dC

S (A)

lS

dy

= 0,

Dl(j) dC

l (j)

y=0

ll

dy

= 0, j = A, B.

(26)

(27)

(28) (29)

у=2

Решить систему уравнений (21)—(29) — означает определить положение фронта кристаллизации, распределение теплоты и концентраций каждого компонента в твердой и жидкой фазах. Эта система нелинейна: от искомой скорости движения фронта зависят длины фаз 15 и 11, а также функция замены переменных ф(у, t). Для аппроксимации уравнений (21)-(29) использовалась разностная неявная консервативная схема, предложенная в работе [6], но модифицированная для рассматриваемой задачи.

Разностная схема. В области Оу введем разностную сетку юу = {уг, 0 < г < N, у0 = 0, у * =1, уN =2}, Ь+У2 = уг+1 -уг — шаг сетки. Пусть граница раздела фаз, точка у = 1 совпадает с одним из узлов этой сетки; обозначим его у *. Наряду с узлами уу, введем потоковые точки уг+1/2 = (уг + уг+1)/2, г = 0,1,.., N -1, и шаги

сетки Йу = у1+1/2 - Уг-ш- Разностную сетку по времени обозначим как ={£о=0, tk+1 = tk + х, к = 0,1,...}, х — шаг сетки по времени. Сеточные функции температуры и концентрации будем относить к узлам

сетки. Пусть / = Т, СУ<у = 5, I, ] = А, В, тогда /(уг,tj) = /. Доопределим эти

функции в интервалах между узлами: /(у, tj ) = /(уг, tj) при у е (уу-1/2, уг+1/2),

г Ф г . В силу непрерывности температуры на границе раздела фаз, Т ¡* - о = Т ¡* + о.

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

^0) ^0) ^0) два ее значения С * = С и С» = С , С , С — равновесные концен-

г - 0 г +0

трации, удовлетворяющие фазовой диаграмме системы; С ' (у, tk) = С0, если у е (у *-1/2, у.*) и С0)( у, tk ) = С+ 0, если у е (у.*, у*+1/2). Для единообразия в

регулярных точках (г Ф г*) также будем использовать обозначения С (-о и С(+о,

С(-) = С(+о = С(j). Далее тильду над Т, СУ<^ опустим. Длина зоны, удельная теплоемкость, коэффициенты температуропроводности, коэффициенты диффузии относятся к полуцелым узлам сетки. Слева и справа от точки г эти функции постоянны и принимают значения, соответствующие твердой и жидкой фазам.

Разностную производную по времени обозначим /г = / = (- /)/ х, где

/ = / (у1, 1 к + х) [17]. Определим операторы разностного дифференцирования по пространству / (г) = [ / (уг+1, tk) - / (уг Л к)] / к+ш, /у (г) = /у (г -1).

Кратко приведем этапы построения разностной схемы для уравнения теплопроводности (21). Детальное описание аппроксимации задачи для концентраций приведено в работе [6]. Проинтегрируем уравнение (21) по разностной ячейке [у{-Ш,у{+ш] х ^,^ +х]:

tkУ+1/2

J I

tk yi-1/2

сp ^(ZyT) - с(9tT)

öt öy

^ yi+1/2

= J J

tk yi-1/2

ö Г Xy ÖT ^

ö X -V(a'T+ßy)

öy ^ Zу öy

dydt =

dydt, у = 5, Z.

(30)

Подставим в равенство (30) сеточные функции и вычислим последовательно интегралы, стоящие в левой и правой части равенства:

tk+х уг я уг+1/2 я

71= | Л[ | Сг-1/2 —(УТ)у + | с,+1/2 —((УТ)у] =

^ уг-1/2 ^ уг ^ (31)

нУ —-— ^у —-—

= С,-1/2 [(¿г-1/2 Т ) - (1г-1/2 Т )] + С,+1/2[(/«- + 1/2 Т ) - (¡,+ 1/2 Т )], У = 5, I.

Если отрезок интегрирования [ /¿-1/2, Уг+1/2] не содержит границу раздела фаз

(г Ф г ), то с,-1/2 = с,+1/2, /г-1/2 = ¡¿+т, и при вычислении интеграла нет необходимости разбивать область интегрирования на отрезки, расположенные слева и справа от точки уг. Представление (31) обеспечивает однородную, не зависящую от номера узла, форму записи разностной аппроксимации производной по времени.

Интеграл от второго слагаемого в левой части (30) также разобьем на две

части и при его вычислении учтем, что (фТ)(у, -0) = (фТ)(y¡ + 0) и ф.( у,) =

tk +т

h = J dt

tk

yi Q У1+1/2 Q

J Ci-1/2(pT)dy + J Ci+1/^— (pT)dy .yi-1/2 Qy yi Qy

tk+x

= J [-1/2 [(фtT)(yi)-(<PtT)(yi-1/2)] + Ci+1/2 [((ptT)(yi+1/2)-(<PtT)(yi)]]dt « (32)

tk

(Ci-1/2 - Ci+1/2 )^t Ti* + Ci+1/2 (pt T)(yi+1/2) - Ci-1/2 (pt T)(yi-1/2)

Произведение (фТ) в полуцелых точках вычисляется по формуле

(фТ)(У1 ±1/2 ) = ф (У1 ±1/2)(Т + Т±1) / 2.

Перейдем к вычислению интеграла, стоящего в правой части равенства (30). Тепловой поток непрерывен в регулярной точке, а на границе раздела фаз испытывает скачок, величина которого определяется из условия Стефана (22), поэтому приближенное выражение для интеграла 13 можно записать в виде

tk+x yi Q (

1з= J dt[ f —

tJ У1 -J1/2 Qy l l y Qy )

кy QT V yi+1/2 Q( к у QT ^ dy + f---

y!i Qyl iy Qy,

■ x ^ A,8 ..* + x[ Ki+1/2 () - к-1/2 ()],

dt Ii+1/2 Ii

dy] =

(33)

Н+1/2 Н-1/2

где 8 » — символ Кронекера.

Для интеграла источника теплоты получим следующую аппроксимацию:

= x

tk+x Х4 = J dt

tk

hi-1/2

yi yi+1/2

J ly (aT + ßy )dy + J ly (ayT + ßy )dy

L yi-1/2 yi

h

(34)

2 -1/2 (ai-1/2 Ti + ßi-1/2 ) + ' +21/2 Ii+1/2 (ai+1/2 Ti + ßi+1/2 )

Заменив в (30) интегралы приближенными формулами (31)-(34), получим умноженную на площадь ячейки дивергентную разностную схему для уравнений (21):

hy 1/2 hy+m

-Ci-1/2 (li-1/2Ti )t +~L--Ci+1/2 (li+ 1/2Ti )t -

2 2

- [Ci+1/2 (pt T)(yi+1/2 ) - Ci-1/2 (pt T)(yi-1/2 )] =

K i+1/2 hi-1/2

( T ^

Ty

v Ii+1/2 )

- Ki-1/2

i T-

1 y

Ii-1/2

+ + (Ci -1/2 Ci+1/2 )T /

(35)

2 /г-1/2(а,_1/2 Тг + Рг-1/2 ) + ' '^^ 1г+1/2 («г +1/2 Тг +Рг+1/2 )

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

Для аппроксимации граничных условий в точках у0 и уы уравнение (21) необходимо проинтегрировать с учетом (28) по отрезкам [ у0, у1/2] и [ уы-1/2, уЫ ]. В результате получим

Y C1/2 (I1/2T0 )t - C1/2 (pt T)(y1/2 ) = K1/2

( t ^

—L

v l1/2 )

-T0) + ^Ma^o +ß1/2); P1/2 2

hy

"n-1/2

CN-1/2 (lN-1/2TN )t + CN-1/2 (ptT)(Ln-1/2 ) =

a N-1/2

pN-1/2

(T* -TN ) - Kn.

-1/2

/ У

lN-1/2

hNy

N-1/2

lN-1/2 (aN-1/2 TN + ßN-1/2 ).

Разностная аппроксимация уравнений диффузии (23) и (25):

hy1/2 (l C(j)) , hi+1/2 (i r(j)) —-—(li-1/2Ci-0)t +---— (li+1/2Ci+0)t -

pt с(j)) (Li+1/2) - (pt C(j)) (Li-1/2)

= D,

( t (j) ^

C y

i+1/2

li+1/2 V )

- D..

i-1/2

( T(-j) ^

C y

li-1/2 V )

(37)

Аппроксимация граничных условий (29) получается, если в (37) при г = 0 и г = N формально принять равными нулю члены с индексами 0-1/2 и N +1/2. В этом нетрудно убедиться, проинтегрировав уравнения (23) и (25) по приграничным ячейкам с учетом условий (29).

Для решения нелинейной системы уравнений (35)-(37) применим метод Ньютона. Введем обозначения ^ = (С(А), Сг(в), Тг )т, г = 0,1,..., г *-1, г * +1,..., N, =(С>(А), С^*(Б), Т>, ^, С?А), С>(в))т. Запишем схему (35)-(37) в операторной форме

ЯО = 0, (38)

2

где переменные упорядочены следующим образом: С = (С0, Сь—, -1, С/*, С*+1, ..., См)т. Определим итерационный процесс по методу Ньютона для системы (38):

Т' (Сп )5С = -Т (Сп ), (39)

где Т'(С) — матрица Якоби системы (38); 5С = Сп+1 ~Сп; п — номер итерации.

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

А-ц 5Сг-1 + А^ 5Сг + А+1,г- 5С,г+1 + в/ 5^ = -Т, г = 0, ..., г * -1, г * +1, ..., N;

А-1,0 = А+1м =0; (40)

А , 5С-* + А..*5С* + А , 5С-* + = -Т*.

-1,г ^г - 0,г ^г + 1,г ^г + г

Здесь А -1,г, А0, ¡, А+и — матрицы 3 х 3; Т — трехкомпонентный вектор;

А-1г>, А+1 > — матрицы 3 х 6; А0г* — матрицы 6 х 6; Т* — шестикомпонент-

ный вектор. Систему уравнений (40) можно решить с помощью модифицированного метода матричной прогонки [18]. Пусть искомые 5С связаны следующими соотношениями:

5Сг = «5+15С+1 +Р5+1 +У5+15^, / = 0,1,...,/* -1;

(41)

5Сг = а^С-! +р?-1 + у I, г = N, N -1,..., С +1.

Здесь ау — матрицы 3х 3; Ру, уу — трехкомпонентные векторы (г = 1,..., N). С помощью соотношений (41) задача (40) сводится к системе уравнений на границе раздела фаз. Для этого последовательно вычислялись коэффициенты а5, Р5, у 5, г = 0,1,..., Г, и а\, р1, у1, г = N, N - 1,...,Г. Затем с использованием коэффициентов а^, р5*, у^, а1*, р1*, у1 * из уравнений, соответствующих границе раздела фаз ( г = Г), исключим неизвестные 5С г*-1, 5С>+1. В результате получим систему из шести уравнений с шестью неизвестными 5С* , которая решается методом Гаусса. По найденным неизвестным 5С* с помощью прогоночных соотношений (41) вычислим приращения , 0< ¿<Г, и г <г<N, 15, /г, на (п + 1)-й итерации по методу Ньютона.

В качестве критерия окончания итерационного процесса используем условие

||5С|| < Сп\ + В2, (42)

где 81, в2 — величины, определяющие точность сходимости итераций.

Алгоритм решения задачи (38) на каждом временном слое состоит из приведенных ниже этапов.

1. Вычисление элементов матрицы Якоби Т '(С) и правых частей Т (С). В качестве начального приближения для искомых величин используется решение с предыдущего временного слоя.

2. Определение прогоночных коэффициентов ау, ру, уу. Сведение исходной задачи к системе из шести уравнений на границе раздела фаз.

3. Расчет 5с0(А), 5с0(В), 5Т0, 5^, 5с0(А), 5с0(В). Решение шести линейных уравнений на границе раздела фаз.

4. Вычисление искомых функций 5С//,0 < г < Г и г*< г < N, I5, 11, ф(.

5. Проверка условия сходимости итераций по методу Ньютона.

6. Переход на следующий шаг по времени.

Результаты расчетов. В качестве примера использования разработанного алгоритма приведем результаты моделирования процессов роста и растворения соединений АХВ1-ХС, фазовая диаграмма системы подробно описана в работе [19]. Совместное определение полей температуры, концентраций и положения фронта кристаллизации позволяет задавать начальные данные естественным образом, без предварительного согласования.

Изотермический рост. Для получения кристалла в режиме изотермического роста, насыщенный при температуре Т1 раствор компонентов А и В в расплаве С приводится в контакт с затравкой состава А ЖВ1_ЖС, и выдерживается при постоянной температуре Т0 = ТЦ-ДТ. Величина ДТ >0 называется степенью пересыщенности расплава. Поскольку жидкая фаза пересыщена, начинается процесс кристаллизации. Основная особенность такого режима — рост кристалла постоянного состава.

Рассматривалась затравка состава А Х0В1-Х0С, х0=0,2, толщиной I5 =

= 30 мкм, длина жидкой фазы 11 =15 см, радиус ампулы 0,5 см, ДТ = 10 К, Т0 =783 К, Г = 500, N - г" =5104, х = 0,1В, Ь =5с, В5 = 5-10-12 см2/с, В1 = = 5-10-5 см2/с, с5р = 441-10-3 Дж/(г-К), с1р = 393-10-3 Дж/(г ■ К), ж5 =

= 2 -10-2 (см2-Вт)/(г-К), ж1 =8-10-2 (см2-Вт)/(г-К). В этом расчете боковая поверхность ампулы считалась теплоизолированной, температура на торцах ампулы постоянна и равна Т0. Распределения концентрации компонентов АС в твердой фазе, компонентов А и В в жидкой фазе приведены на рис. 2. Поскольку раствор является пересыщенным, в жидкой фазе выстраивается концентрационный градиент, способствующий росту кристалла постоянного состава = 0,184. Температура на фронте кристаллизации не изменяется, поэтому состав на границе жидкой фазы также остается постоянным.

Известны экспериментальные данные о составе выросшего кристалла. Зависимость состава твердой фазы от степени первоначального пересыщения раствора ДТ приведена в работах [20-22]. Расчеты проводились для ДТ = 2,5,7,10 К.

Рис. 2. Содержание АС в твердой фазе и состав жидкой фазы при изотермическом росте (Х0=0,2, tl=0tD, t2=30tD, t3=100tD, г4 =300^в)

Падение содержания АС на один градус пересыщения Ах = (х0 - хе)/АТ, а также аппроксимация полученных результатов линейной функцией представлены на рис. 3. Мольная доля АС в кристалле убывает с ростом температуры АТ. Зависимость хе (АТ) с хорошей точностью описывается функцией хе = 0,2 - 0,0016 АТ, что согласуется с данными, приведенными в работах [20-22].

Отметим следующие достоинства используемого вычислительного алгоритма. Во-первых, возможность задания начальных данных для расчетов естественным образом. При изотермическом режиме выращивания в начальный момент состав расплава не равновесен затравке, т. е. концентрации веществ в твердой и жидкой фазах не удовлетворяют уравнениям фазовой диаграммы. Однако для предложенного алгоритма в рассмотренном диапазоне параметров не требуется специального согласования начальных данных, о котором, например, упоминается в работах [7, 23, 24]. Во-вторых, численный метод является консервативным. Выполнение закона сохранения массы контролировалось на каждом шаге по времени. Масса каждого компонента оставалась постоянной с точностью, определяемой точностью решения системы уравнений.

Рис. 3. Зависимость содержания АС в выросшем кристалле от степени пересыщения раствора АТ (х0 =0,2):

1 — аппроксимация; 2 — результаты расчета

Рост в условиях принудительного охлаждения. В начальный момент вре- 0,207 мени затравка ЛХБ1-ХС приводится в о,206 контакт с равновесным ей раствором Л-Б-С. Температура на торцах ампулы изменяется по закону: Т(£) = Т0 - ц, где Т0 — начальная температура соединения; ц — скорость охлаждения; I —

0,205 0,204 0,203 0,202

Затравка Кристалл

> | V \\ \\ V5 \ \ \ \ \ 1 1 1

время; боковая поверхность ампулы

г „ 7 0,99950 0,99975 1,00000 1,00025 z, см

теплоизолированная. Вследствие понижения температуры раствор становится Рис. 4. Распределение состава твердой пересыщенным, начинается процесс кри- фазы на момент времени t = 10 мин сталлизации. Распределение состава твер- (*о =0,207, £ = 10 мин, Г0 =783 К) при дой фазы на момент времени при раз- скорости охлаждения 0,1 (1), 0,25 (2), личных значениях скорости t = 10 мин 0,5 (3) К/мин

охлаждения ц приведено на рис. 4,

начальная толщина твердой фазы составляет 1 см. В расчетах число узлов в твердой фазе составляло г = 5-104, остальные параметры приняты такими же, как и в предыдущих расчетах. С увеличением скорости ц рост кристалла происходит быстрее. Так, при ц = 0,1 К/мин за 10 мин толщина твердой фазы увеличивается на 1 мкм, а при ц = 0,5 К/мин — на 4,9 мкм. Кроме того, содержание АС в твердой фазе снижается, поскольку кристаллизация происходит при более низкой температуре [20, 21].

Проведено исследование влияния теплообмена между боковой поверхностью ампулы и окружающей средой на процесс кристаллизации. На торцах ампулы задавали условия охлаждения, а на боковой поверхности — условия контактного теплообмена (11). Результаты расчетов показали, что если температура окружающей среды ниже, чем температура внутри ампулы, то теплообмен выступает в качестве дополнительного механизма охлаждения, поэтому рост кристалла происходит 0,05 быстрее (рис. 5). Таким образом, темпера-

туру окружающей среды можно исполь-0 2 4 6 8 мин зовать как дополнительный регулятор

Рис. 5. Зависимость скорости движения скорости роста кристалла.

границы раздела фаз от времени Рост кРисталла из недосыЩенного

Вт/(К ■ см2), К/мин, раствора. При температуре Т0 нед°сы-

К) при температуре окружающей щенный раств°р к°мп°нентов А и Б

среды = 780 (1), 773 (2) К и при в расплаве С, равновесный тверд°й фазе

изоляции (3) Х0 при температуре Т = Т -АТ, поме-

щается в ампулу, на дне которой расположена затравка состава Л^Б^С. Торцы ампулы охлаждаются по закону Т(р) = Т0 до температуры Т2, затем система выдерживается при постоянной температуре. Здесь ЛТ >0 характеризует степень недосыщения жидкой фазы. Боковая поверхность ампулы теплоизолированная1. При таком режиме выращивания насыщение раствора происходит как за счет понижения температуры, так и за счет растворения затравки, в результате которого в жидкую фазу поступает дополнительный материал (рис. 6). С течением времени раствор в окрестности границы раздела фаз становится насыщенным, и продолжающееся понижение температуры приводит к смене растворения ростом кристалла при I = 10,5 мин. Вблизи фронта кристаллизации выстраивается концентрационный градиент, способствующий росту кристалла ( = 14,5 мин).

0,2 0,4 0,6 0,8 1,0 z, см

Рис. 6. Содержание ЛС в твердой фазе и состав жидкой фазы соединения

К/мин, при t = 1 с (1), 1 мин (2), 10,5 мин (3) и 14,5 мин (4)

K,

Эволюция состава в системе на этапе, когда охлаждение ампулы прекращено, представлена на рис. 7. Расчеты показывают, что происходит рост кристалла постоянного состава. При этом предварительное подрастворение подложки приводит к тому, что в твердой фазе присутствует зона, в которой концентрация ЛС заметно отличается от начального состава подложки. Высота и ширина концентрационного «всплеска» уменьшаются с течением времени за счет диффузии в твердой фазе. Вследствие большой длины ампулы и достаточно малого коэффициента диффузии в жидкой фазе в окрестности фронта кристаллизации наблюдается зона локального пересыщения раствора.

1 В расчетах использованы те же параметры, что и для изотермических режимов.

АС 0,212

0,211

0,210

0,209

0,218

0,207

Затравка I

1 Кристалл

1 л А

/ А

А

А

10 15

В

0,1770 0,1767

0,1764

; А 0,01112

0,01104

0,01096

20 г, мкм 0,5 1,0 1,5 2,0 г, см

Рис. 7. Содержание АС в твердой фазе и состав жидкой фазы соединения на этапе роста эпитаксиального слоя постоянного состава (ДТ = 10 К, ц = 0,5 К/мин, х0 =0,207)

при t = 23 (1), 35 (2) и 52 (3) мин

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

ЛИТЕРАТУРА

1. Самарский А.А., Вабищевич П.Н. Вычислительная теплопередача. М.: УРСС, 2003. 784 с.

2. Samarskii A.A., Vabishchevich P.N., Iliev O.P., Churbanov A.G. Numerical simulation of convection/diffusion phase change problems — a review // Journal of Heat Mass Transfer. 1993. Vol. 36. No. 17. P. 4095-4106.

3. Muray W.D., Landis F. Numerical and machine solutions of the transient heat conduction problems involving melting or freezing // Journal of Heat Transfer. 1959. Vol. 81. P. 106-112.

4. Vermolen F.J., Vuik C. A mathematical model for the dissolution of particles in multicomponent alloys // J. of Computational and Applied Math. 2000. Vol. 126. No. 1-2. P. 233-254.

DOI: 10.1016/S0377-0427(99)00355-6

URL: http://www.sciencedirect.com/science/article/pii/S0377042799003556

5. Мажорова О.С., Попов Ю.П., Щерица О.В. Алгоритм расчета задачи о фазовом переходе в многокомпонентной системе // Дифференциальные уравнения. 2004. Т. 40. № 7. С. 1051-1060.

6. Мажорова О.С., Попов Ю.П., Щерица О.В. Консервативные разностные схемы для термодиффузионной задачи Стефана // Дифференциальные уравнения. 2013. Т. 49. № 7. С. 897-905.

7. Illingworth T.C., Golosnoy I.O. Numerical solutions of diffusion-controled moving boundary problems which conserve solute // Journal of Computation Physics. 2005. Vol. 209. No. 1. P. 207-225. DOI: 10.1016/j.jcp.2005.02.031

URL: http://www.sciencedirect.com/science/article/pii/S0021999105000859

8. Бакирова О.И. Численное моделирование процесса зонной плавки на основе решения задачи о фазовом переходе в бинарной системе // Математическое моделирование. Получение металлов и полупроводниковых структур. М.: Наука, 1986. С. 142-158.

9. Дегтярев Л.М., Дроздов В.В., Иванова Т.С. Метод адаптивных к решению сеток в сингулярно-возмущенных одномерных краевых задачах // Дифференциальные уравнения. 1987. Т. 23. № 7. С. 1160-1169.

10. Pandelaers L, Verhaeghe F., Wollants P., Blanpain B. An implicit conservative scheme for coupled heat and mass transfer problems with multiple moving interfaces // Int. J. of Heat and Mass Transfer. 2011. Vol. 54. No. 5-6. P. 1039-1045.

11. Landau H.G. Heat conduction in a melting solid // J. App. Math. 1950. Vol. 8. P. 81-94.

12. Chtcheritsa O.V., Mazhorova O.S., Popov Yu.P. Implicit numerical algorithm for the solution of phase transition problems in multi-component alloys // Mathematical Modelling and Analysis. 2004. Vol. 9. No. 4. P. 253-266.

URL: http://www.tandfonline.com/doi/abs/10.1080/13926292.2004.9637258

13. Мажорова О.С., Попов Ю.П., Похилко В.И. Матричный алгоритм численного решения нестационарных задач концентрационной конвекции для многокомпонентных сред // Получение монокристаллов и полупроводниковых структур / под ред. А.А. Самарского, Ю.П. Попова, О.С. Мажоровой. М.: Наука, 1986. C. 19-31.

14. Ghez R., Small M.B. Growth and dissolution kinetics of ternary alloys of ternary III-V hetero-structures formed by liquid phase epitaxy. III. Effect of temperature programming // Journal of Applied Physics. 1982. Vol. 53. No. 7. P. 4907-4918. DOI: 10.1063/1.331324

15. Shcheritsa O.V., Mazhorova O.S., Popov Yu.P. Numerical study for diffusion processes in dissolution and growth of CdHgTe/CdTe heterostructures formed by LPE. Part I. Isothermal conditions // Journal of Crystal Growth. 2006. Vol. 290. No. 2. P. 357-362.

16. Мажорова О.С., Попов Ю.П., Щерица О.В. Чисто неявный метод решения задач о фазовом переходе // Препринт института прикладной математики им. М.В. Келдыша РАН. 2004. № 29. 42 с. URL: http://library.keldysh.ru/preprint.asp?id=2004-29

17. Самарский А.А. Введение в теорию разностных схем. М.: Наука, 1971. 552 с.

18. Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. М.: Наука, 1978. 592 с.

19. Денисов И.А., Лакеенков В.М., Мажорова О.С., Попов Ю.П. Математическое моделирование эпитаксиального выращивания твердых растворов CdYHg1-YTe из жидкой фазы // Препринт института прикладной математики им. М.В. Келдыша РАН. 1992. № 65. 42 с.

20. Denisov I.A., Lakeenkov V.M., Mazhorova O.S., Popov Yu.P. Numerical modelling for liquid phase epitaxy of CdxHg1-xTe solid solution // Journal of Crystal Growth. 2002. Vol. 245. No. 1-2. P. 21-30.

21. Denisov I.A., Mazhorova O.S., Popov Yu.P., Smirnova N.A. Numerical modelling for convection in growth/dissolution of solid solution CdxHg1-xTe by liquid phase epitaxy // Journal of Crystal Growth. 2004. Vol. 269. No. 2-4. P. 284-291.

22. Sanz-Maudes J., Sangador J., Rodriguez T., et al. Numerical simulation of the growth of HgCdTe layers by liquid phase epitaxy from Te-rich solutions: The effect of liquid dimensions and mercury loss // Journal of Crystal Growth. 1990. Vol. 106. No. 2-3. P. 303-317.

23. Dost S., Qin Z., Kimura M. A model for convective mass transport in liquid phase epitaxial growth of semiconductors // Journal of Heat Mass Transfer. 1997. Vol. 40. No. 13. P. 3039-3047.

24. Kimura M., Qin Z., Dost S. A solid — liquid diffusion model for growth and dissolution of ternary alloys by liquid phase epitaxy // Journal of Crystal Growth. 1996. Vol. 158. P. 231-240.

Щерица Ольга Владимировна — канд. физ.-мат. наук, научный сотрудник Института прикладной математики им. М.В. Келдыша РАН (Российская Федерация, 125047, Москва, Миусская пл., д. 4), доцент кафедры «Прикладная математика» МГТУ им. Н.Э. Баумана (Российская Федерация, 105005, Москва, 2-я Бауманская ул., д. 5, стр. 1).

Гусев Андрей Олегович — бакалавр кафедры «Прикладная математика» МГТУ им. Н.Э. Баумана (Российская Федерация, 105005, Москва, 2-я Бауманская ул., д. 5, стр. 1).

Мажорова Ольга Семеновна — д-р физ.-мат. наук, главный научный сотрудник Института прикладной математики им. М.В. Келдыша РАН (Российская Федерация, 125047, Москва, Миусская пл., д. 4), профессор МГУ им. М.В. Ломоносова (Российская Федерация, 119991, Москва, ул. Ленинские Горы, д. 1), профессор кафедры «Прикладная математика» МГТУ им. Н.Э. Баумана (Российская Федерация, 105005, Москва, 2-я Бауманская ул., д. 5, стр. 1).

Просьба ссылаться на эту статью следующим образом:

Щерица О.В., Гусев А.О., Мажорова О.С. Об одном методе решения задачи кристаллизации многокомпонентного раствора в цилиндрической ампуле // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. 2017. № 5. C. 118-138.

DOI: 10.18698/1812-3368-2017-5-118-138

ON SOLUTION TO PHASE TRANSITION PROBLEM IN MULTICOMPONENT ALLOY IN THE CYLINDRICAL AMPULE

O.V. Shcheritsa1, 2 A.O. Gusev2 O.S. Mazhorova1, 2 3

[email protected]

[email protected]

[email protected]

1 Keldysh Institute of Applied Mathematics, Russian Academy of Sciences, Moscow, Russian Federation

2 Bauman Moscow State Technical University, Moscow, Russian Federation

3 Lomonosov Moscow State University, Moscow, Russian Federation

Abstract

The paper presents a self-consistent model of multicompo-nent alloy crystallization in a cylindrical ampule. The mathematical model accounts for the heat and matter transfer in both solid and liquid phases. We described the system by the interface position and radially average temperature and concentrations. Special efforts are required to solve a corresponding one dimensional phase transition problem in multicomponent alloy. To handle evolution of solid/liquid interface, the moving boundary problem is mapped to a new coordinate system. We obtained a conservative and implicit finite difference scheme in a new coordinate system with control volume technique and constructed a fully implicit coupled approach. Furthermore, we solved a corresponding set of nonlinear equations by Newton method for the unknown vector, whose components are concentrations of all species, interface rate and temperature. The proposed method was used for numerical simulation of the crystallization process of A-B—C solution

Keywords

Stefan problem, phase transition, mathematical simulation

Received 24.10.2016 © BMSTU, 2017

REFERENCES

[1] Samarskii A.A., Vabishchevich P.N. Vychislitel'naya teploperedacha [Computational heat transfer]. Moscow, URSS Publ., 2003. 784 p.

[2] Samarskii A.A., Vabishchevich P.N., Iliev O.P., Churbanov A.G. Numerical simulation of convection/diffusion phase change problems — a review. Journal of Heat Mass Transfer, 1993, vol. 36, no. 17, pp. 4095-4106.

[3] Muray W.D., Landis F. Numerical and machine solutions of the transient heat conduction problems involving melting or freezing. Journal of Heat Transfer, 1959, vol. 81, pp. 106-112.

[4] Vermolen F.J., Vuik C. A mathematical model for the dissolution of particles in multicomponent alloys. J. of Computational and Applied Math., 2000, vol. 126, no. 1-2, pp. 233-254. DOI: 10.1016/S0377-0427(99)00355-6

Available at: http://www.sciencedirect.com/science/article/pii/S0377042799003556

[5] Mazhorova O.S., Popov Yu.P., Shcheritsa O.V. An algorithm for solving a phase transition problem in a multicomponent system. Differential Equations, 2004, vol. 40, no. 7, pp. 10511059. DOI: 10.1023/B:DIEQ.0000047035.96793.be

Available at: https://link.springer.com/article/10.1023/B%3ADIEQ.0000047035.96793.be

[6] Mazhorova O.S., Popov Yu.P., Shcheritsa O.V. Conservative scheme for the thermodiffusion Stefan problem. Differential Equations, 2013, vol. 49, no. 7, pp. 869-882.

DOI: 10.1134/S0012266113070094

Available at: https://link.springer.com/article/10.1134/S0012266113070094

[7] Illingworth T.C., Golosnoy I.O. Numerical solutions of diffusion-controled moving boundary problems which conserve solute. Journal of Computation Physics, 2005, vol. 209, no. 1, pp. 207-225. DOI: 10.1016/j.jcp.2005.02.031

Available at: http://www.sciencedirect.com/science/article/pii/S0021999105000859

[8] Bakirova O.I. Chislennoe modelirovanie protsessa zonnoy plavki na osnove resheniya zadachi o fazovom perekhode v binarnoy sisteme. Matematicheskoe modelirovanie. Poluchenie metallov i poluprovodnikovykh struktur [Numerical simulation of zone melting process based on solution of problem of phase transition in binary system. In: Math. modeling. Metals and semiconductor structures production]. Moscow, Nauka Publ., 1986, pp. 142-158 (in Russ.).

[9] Degtyarev L.M., Drozdov V.V., Ivanova T.S. The method of nets adapted to the solution in singularly perturbed one-dimensional boundary value problems. Differentsial'nye uravne-niya, 1987, vol. 23, no. 7, pp. 1160-1169 (in Russ.).

[10] Pandelaers L., Verhaeghe F., Wollants P., Blanpain B. An implicit conservative scheme for coupled heat and mass transfer problems with multiple moving interfaces. Int. J. of Heat and Mass Transfer, 2011, vol. 54, no. 5-6, pp. 1039-1045.

[11] Landau H.G. Heat conduction in a melting solid. J. App. Math., 1950, vol. 8, pp. 81-94.

[12] Chtcheritsa O.V., Mazhorova O.S., Popov Yu.P. Implicit numerical algorithm for the solution of phase transition problems in multi-component alloys. Mathematical Modelling and Analysis, 2004, vol. 9, no. 4, pp. 253-266.

Available at: http://www.tandfonline.com/doi/abs/10.1080/13926292.2004.9637258

[13] Mazhorova O.S., Popov Yu.P., Pokhilko V.I. Matrichnyy algoritm chislennogo resheniya nestatsionarnykh zadach kontsentratsionnoy konvektsii dlya mnogokomponentnykh sred. Matematicheskoe modelirovanie. Poluchenie monokristallov i poluprovodnikovykh struktur [Matrix algorithm of numerical solution of non-stationary concentration-induced convection problems in multicomponent medium. Monocrystal and semiconductor structure production]. Moscow, Nauka Publ., 1986, pp. 19-31 (in Russ.).

[14] Ghez R., Small M.B. Growth and dissolution kinetics of ternary alloys of ternary III-V heterostructures formed by liquid phase epitaxy. III. Effect of temperature programming. Journal of Applied Physics, 1982, vol. 53, no. 7, pp. 4907-4918. DOI: 10.1063/1.331324

[15] Shcheritsa O.V., Mazhorova O.S., Popov Yu.P. Numerical study for diffusion processes in dissolution and growth of CdHgTe/CdTe heterostructures formed by LPE. Part I. Isothermal conditions. Journal of Crystal Growth, 2006, vol. 290, no. 2, pp. 357-362.

[16] Mazhorova O.S., Popov Yu.P., Shcheritsa O.V. Implicit numerical algorithm for solution of phase transition problems. Preprint instituta prikladnoy matematiki im. M. У. Keldysha RAN [KIAM Preprint], 2004, no. 29, 42 p.

Available at: http://library.keldysh.ru/preprint.asp?id=2004-29

[17] Samarskii A.A. Vvedenie v teoriyu raznostnykh skhem [Introduction to the theory of difference schemes]. Moscow, Nauka Publ., 1971. 552 p.

[18] Samarskiy A.A., Nikolaev E.S. Metody resheniya setochnykh uravneniy [Finite-difference equation solution methods]. Moscow, Nauka Publ., 1978. 592 p.

[19] Denisov I.A., Lakeenkov V.M., Mazhorova O.S., Popov Yu.P. Mathematical simulation of epitaxial growing of solid solutions CdyHg1-rTe liquid phase. Preprint instituta prikladnoy matematiki im. M.V. Keldysha RAN [KIAM Preprint], 1992, no. 65, 42 p. (in Russ.).

[20] Denisov I.A., Lakeenkov V.M., Mazhorova O.S., Popov Yu.P. Numerical modelling for liquid phase epitaxy of Cd^Hg^Te solid solution. Journal of Crystal Growth, 2002, vol. 245, no. 1-2, pp. 21-30.

[21] Denisov I. A., Mazhorova O.S., Popov Yu.P., Smirnova N.A. Numerical modelling for convection in growth/dissolution of solid solution Cd^Hg^Te by liquid phase epitaxy. Journal of Crystal Growth, 2004, vol. 269, no. 2-4, pp. 284-291.

[22] Sanz-Maudes J., Sangador J., Rodriguez T., et al. Numerical simulation of the growth of HgCdTe layers by liquid phase epitaxy from Te-rich solutions: The effect of liquid dimensions and mercury loss. Journal of Crystal Growth, 1990, vol. 106, no. 2-3, pp. 303-317.

[23] Dost S., Qin Z., Kimura M. A model for convective mass transport in liquid phase epitaxial growth of semiconductors. Journal of Heat Mass Transfer, 1997, vol. 40, no. 13, pp. 30393047.

[24] Kimura M., Qin Z., Dost S. A solid — liquid diffusion model for growth and dissolution of ternary alloys by liquid phase epitaxy. Journal of Crystal Growth, 1996, vol. 158, pp. 231-240.

Shcheritsa O.V. — Cand. Sc. (Phys.-Math.), Researcher of Keldysh Institute of Applied Mathematics, Russian Academy of Sciences (Miusskaya ploshchad' 4, Moscow, 125047 Russian Federation), Assoc. Professor of Applied Mathematics Department, Bauman Moscow State Technical University (2-ya Baumanskaya ul. 5, str. 1, Moscow, 105005 Russian Federation).

Gusev A.O. — Bachelor' Degree student of Applied Mathematics Department, Bauman Moscow State Technical University (2-ya Baumanskaya ul. 5, str. 1, Moscow, 105005 Russian Federation).

Mazhorova O.S. — Dr. Sc. (Phys.-Math.), Senior Researcher of Keldysh Institute of Applied Mathematics, Russian Academy of Sciences (Miusskaya ploshchad' 4, Moscow, 125047 Russian Federation), Professor of Lomonosov Moscow State University (Leninskie gory 1, Moscow, 119991 Russian Federation), Professor of Applied Mathematics Department, Bauman Moscow State Technical University (2-ya Baumanskaya ul. 5, str. 1, Moscow, 105005 Russian Federation).

Please cite this article in English as:

Shcheritsa O.V., Gusev A.O., Mazhorova O.S. On Solution to Phase Transition Problem in Multicomponent Alloy in the Cylindrical Ampule. Vestn. Mosk. Gos. Tekh. Univ. im. N.E. Baumana, Estestv. Nauki [Herald of the Bauman Moscow State Tech. Univ., Nat. Sci.], 2017, no. 5, pp. 118-138. DOI: 10.18698/1812-3368-2017-5-118-138

Издательство МГТУ им. Н.Э. Баумана 105005, Москва, 2-я Бауманская ул., д. 5, стр. 1 [email protected] www.baumanpress.ru

Подписано в печать 02.10.2017 Формат 70 х 108/16 Усл.-печ. л. 12,0 Отпечатано в ПАО «Т8 Издательские Технологии» 109316, Москва, Волгоградский пр-т, д. 42, корп. 5 Тираж 100 экз.