Моделирование кристаллизации метастабильного расплава

( Numerical Modeling of the Solidification Process of a Metastable Melt
Preprint, Inst. Appl. Math., the Russian Academy of Science)

Содержание

Введение……………………………….……………………………………….…..3

§1. Уравнения модели фазового поля...................................................………….6

§2. Постановка задачи для плоского фронта кристаллизации...........…….......11

§3. Разностная схема и алгоритм адаптации сетки к решению...........……….12

§4. Результаты тестовых расчетов......……………………………..………..….15

§5. Анализ квазистационарного режима кристаллизации..……..………...….17

Заключение……...……………………..……………………………………….…29

Литература...…………………………..……………………………………….….29

Введение

Модель процесса кристаллизации расплава, в которой фазовый переход (фронт кристаллизации) рассматривается в приближении подвижной межфазной поверхности, известна как задача Стефана [1]. Задача Стефана состоит из уравнений теплопроводности в жидкой  и твердой  фазах (уравнений сохранения энергии)

                                                                       (1)

и двух граничных условий на подвижной поверхности  разделяющей фазы: уравнения Стефана (уравнения баланса тепла)

                                              (2)

и уравнения, определяющего отклонение температуры межфазной границы от равновесной температуры фазового превращения  – граничное условие Гиббса-Томсона [2, 3]

,        (3)

где  – температура (расплава или твёрдой фазы, соответственно);  – скорость движения межфазной поверхности;  – температуропроводность среды;  – коэффициент теплопроводности;  – теплоёмкость и скрытая теплота фазового превращения единицы объёма; – капилярная длина (капилярная постоянная); – плотность энергии поверхностного натяжения; ;  – кинетический коэффициент роста;  – средняя кривизна межфазной поверхности;  – единичный вектор, нормальный к межфазной поверхности и направленный из твёрдой фазы в жидкую. В классическом варианте задачи Стефана отклонение температуры межфазной поверхности от равновесной температуры фазового превращения из-за кривизны и кинетическое (динамическое) переохлаждение поверхности не учитываются, т.е. , и на фазовом переходе выполняется условие локального термодинамического равновесия

.                                                                          (4)

Решение Задачи Стефана, в постановке (1) – (3), состоит в определении температуры расплава и твердой фазы, а так же координат межфазной поверхности и её скорости перемещения. При кристаллизации переохлаждённого расплава, начальная температура расплава  ниже равновесной температуры , т.е. расплав изначально находится в метастабильном состоянии. В одномерном случае (плоский фронт кристаллизации) для неограниченной области  для задачи (1) – (3) с начальными условиями (задачи Коши) получено интегрального уравнения для скорости движения межфазной поверхности [4].       При невысокой степени переохлаждения расплава  и начальных условиях

          (5)

показано [4], что решение задачи Коши (1) – (3), (5) при  асимптотически стремиться к хорошо известному автомодельному решению классической задачи Коши-Стефана (1), (2), (4), (5) (одностороннее автомодельное решение) кристаллизации в переохлажденный расплав

                   (6)

где 

а  является решением трансцендентного уравнения

           (7)

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

При высокой степени переохлаждения расплава  (гиперохлаждение), теплоты фазового перехода уже недостаточно, чтобы разогреть расплав до равновесной температуры фазового превращения  и, как следствие, одномерная задача Стефана в классической постановке (1), (2), (4) решения не имеет. Твердая фаза, образующаяся в таких условиях кристаллизации, имеет температуру ниже равновесной и граничное условие (4) модифицируется в условие (3) введением поправки, учитывающей кинетическое переохлаждение. Величина кинетического переохлаждения межфазной поверхности связана со скоростью движения этой поверхности в уравнении (3) линейным образом посредством кинетического коэффициента . Линейная зависимость справедлива при нормальном механизме роста кристаллов из расплава [5, 6, 8, 9] и наблюдалась экспериментально [7]. Одномерная модифицированная модель Стефана (1) – (3) кристаллизации в гиперохлажденный расплав с граничными условиями (задача без начальных условий или задача об установлении)

имеет квазистационарное решение

        (8)

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

            Сравнение результатов исследований морфологической устойчивости фронта кристаллизации, проведенных на основе модели (1) – (3), с имеющимися экспериментальными наблюдениями показало, что для адекватного описания роста дендритной структуры при кристаллизации метастабильного расплава необходимо учитывать анизотропию энергии поверхностного натяжения и кинетического коэффициента, т.е. зависимость этих параметров от пространственной ориентации осей симметрии растущей кристаллической фазы. Другая трудность, возникающая при использование моделей типа (1) – (3), связана с необходимостью явного вычисления координат и скорости движения межфазной поверхности, поскольку только величины отклонения температуры от равновесной температуры фазового превращения уже недостаточно для идентификации фаз, как в классической постановке задачи Стефана (1), (2), (4). В случае пересечения межфазных поверхностей, например при явлении коалесценции, возникают затруднения с определением кривизны и нормали в точке пересечения поверхностей.

Альтернативной моделью процесса кристаллизации, появившейся и интенсивно развивающейся в последнее время, является модель фазового поля [10 – 34, 52]. В рамках этой модели фазовый переход рассматривается как переходная область конечной протяженности, разделяющая фазы. Для описания термодинамического состояния такой двухфазной системы, кроме температуры и плотности используется дополнительный внутренний пространственно-неоднородный параметр системы – параметр порядка. Пространственное распределение  параметра порядка называется фазовым полем.

Большое количество работ посвящено сравнению модели фазового поля кристаллизации расплава с моделями Стефановского типа [14 – 34]. Общий вывод, который следует из сравнительного количественного анализа моделей, заключается в том, что динамика фронта кристаллизации и распределение температуры вне фазового перехода, полученные на основе модели фазового поля, практически совпадают с решением задачи Стефана, если длина теплового поля много больше ширины фазового перехода. Таким образом, модель фазового поля, по-видимому, является более общей моделью кристаллизации и включает модель Стефана как предельный случай, когда ширина межфазной области стремиться к нулю. В силу большей общности, на основе модели фазового поля могут быть получены режимы кристаллизации, которые невозможно реализовать в рамках модели Стефана. Анализ устойчивости квазистационарного решения одномерной задачи кристаллизации переохлаждённого расплава () [20, 27, 29 – 32] показал возможность существования режима кристаллизации, при котором образующаяся твёрдая фаза находится в перегретом (метастабильном) состоянии, т.е. её температура выше равновесной температуры фазового превращения.

В данной работе воспроизводится подробный вывод термодинамически-согласованных уравнений модели фазового поля для кристаллизации расплава чистого металла. Формулируются граничные условия для плоского фронта кристаллизации метастабильного (переохлажденного) расплава и рассматривается разностная схема для численного решения уравнений модели на сетке, динамически адаптирующейся к решению. Результаты тестовых расчетов сравниваются с решениями (6), (7), (8) модели (1) – (3). Анализируются квазистационарные решения одномерной задачи кристаллизации расплава различной степени переохлаждения. На основе прямого численного моделирования, определён диапазон параметров модели, при которых возможно формирования твёрдой фазы в перегретом состоянии и проведено сравнение этих результатов с имеющимися теоретическими оценками.

§1. Уравнения модели фазового поля

В континуальных моделях [10 – 15, 19, 21, 26, 29] неравновесных термодинамических систем предполагается, что состояние системы определяется не только плотностью  и температурой , но и пространственно-неоднородным полем дополнительного внутреннего безразмерного параметра системы  (или нескольких параметров), характеризующих отклонение системы от равновесия [35, 36]. Процесс релаксации в таких моделях рассматривается как эволюция этого параметра. В феноменологической теории фазовых переходов [37 – 39] таким релаксационным параметром является параметр порядка. Оставаясь на термодинамическом уровне рассмотрения фазового перехода как неравновесной системы, можно дать только физическую интерпретацию параметра порядка, а точное его определение возможно лишь в рамках микроскопических теорий фазовых переходов [45 – 48]. Фазовый переход в континуальной модели уже не является поверхностью, как в модели Стефановского типа, а представляет собой межфазную область конечной протяженности. Параметр порядка меняется непрерывным образом по ширине этой области, описывая внутреннюю структуру фазового перехода, а вдали от неё имеет постоянную величину, соответствующую фазе.

Предполагается так же, что термодинамический потенциал (свободная энергия Гиббса) всей системы в неравновесном состоянии определен [36]

                                                                       (9)

и плотность термодинамического потенциала задается функционалом Ландау – Гинзбурга [37 – 39] Кана – Хиллиарда [41 – 44]

                     (10)

где объём, занимаемый системой; вариационная производная функционала . Для адиабатически-изолированной системы, т.е. при отсутствие  градиентов температуры и параметра порядка на поверхности , ограничивающей объём

.                                                           (11)

Система рассматривается при постоянном давление. Изменение плотности на фазовом переходе не учитывается, хотя рассматриваемая модель обобщается на этот случай [21]. Для двухфазной системы, функция  при постоянной температуре представляет собой модельный потенциал с двумя “ямами” (рис. 1), иногда его называют синергетическим потенциалом. Точки экстремума  потенциала

определяют устойчивые  однородные (гомогенные) состояния системы при постоянной температуре и неустойчивое состояние .

Как это видно из рис.1, при  фаза  (твердая) находится в стабильном состоянии, фаза  (расплав) находится в метастабильном состоянии. При равновесной температуре фазового превращения , и как расплав, так и твёрдая фаза являются стабильными состояниями системы.

                                                                                      Рис. 1

Плотность энтропии  и энтальпии  системы определяются потенциалом (9) и для них справедливы следующие термодинамические соотношения

              (12)

       (13)

                   (14)

             (15)

где  теплоёмкость единицы объёма среды. Из (10) – (12), (15) следует

   (16)

и изменение плотности энтальпии , учитывая (13), можно представить в виде

.   (`17)

В соответствии с первым началом термодинамики, количество тепла, полученного единицей объёма системы , где плотность внутренней энергии системы и  при постоянной плотности. Предполагая кондуктивный механизм теплопереноса, т.е. закон Фурье для теплового потока

,                                                                                    (18)

 тогда  и, принимая во внимание (17), можно получить уравнение сохранения энергии (уравнение для температуры) модели фазового поля

.                       (19)

Для вывода уравнения эволюции параметра порядка рассматривается изменение энтропиив произвольном объёме  двухфазной системы с течением времени [19, 21, 26]

.   (20)

Из соотношений (17), (19) следует

,

и, поскольку

,

уравнение (20) приводится к виду

,                              (21)

где  поверхность объёма . Используя (18), уравнению (21) окончательно можно придать более наглядную форму

.               (22)

Интеграл в левой части (22) представляет собой поток энтропии через границу  объёма , обусловленный теплопроводностью (плотность потока энтропии ). Интеграл в правой части является источником энтропии в объёме . Производство энтропии обусловлено теплопереносом

и процессом релаксации системы к равновесию. Мощность релаксационного источника энтропии должна быть, в соответствии со вторым началом термодинамики, так же неотрицательна

.                                                      (23)

Наиболее простое уравнение, удовлетворяющее условию (23), которым описывается процесс релаксации двухфазной системы к равновесию, будет

                                                     (24)

где  – коэффициент Онсагера. Уравнение эволюции поля параметра порядка типа (24) известно как временное уравнение Ландау-Гинзбурга [39 – 40], хотя подобное уравнение использовалось и раннее для описания эволюции неравновесных систем [35].

Изменение свободной энергии рассматриваемой неравновесной системы с течением времени можно описать следующим уравнением

.                          (25)

В изотермическом случае  и из (25) следует

,

что находится в соответствии с известной теоремой об убывании свободной энергии адиабатически-изолированной системы при необратимом изотермическом процессе [36]. В состоянии равновесия свободная энергия системы достигает минимального значения и .

Неизвестные параметры модели в уравнение (24) связаны с физически-измеряемыми свойствами системы. Полагая теплоёмкость единицы объёма среды  величиной постоянной, независящей от , для плотности энтальпии , энтропии  и термодинамического потенциала  справедливы соотношения

                                                                    (26)

                                                         (27)

        (28)

где температурно-независимые части плотности энтальпии и энтропии, соответственно, скрытая теплота фазового перехода. Тогда, уравнения модели (19), (24) для безразмерной температуры  и параметра порядка приобретают вид

                                                 (29)

                                                      (30)

где       ,

.

При равновесной температуре фазового превращения предполагается, что

       (31)

где  – расплав,  – твёрдая фаза. Функция  – температурно-независимая часть энтропии может быть выбрана различными способами [18, 20, 21, 26], в дальнейшем используется

.                                                                    (32)

Для плоского, неподвижного фазового перехода в состоянии равновесия, уравнение  параметра порядка (30) для потенциала (31) будет

.                                  (33)

Решение этого уравнения с граничными условиями

             (34)

известно

                                                                   (35)

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

.                                           (36)

Поскольку, как это следует из (31), (33), (35),

интеграл в (36) вычисляется

                                             (37)

Уравнение параметра порядка (30) для плоского, изотермического , квазистационарного  фазового перехода заменой переменных  преобразуется к виду

,                                                                  (38)

                                       (39)

Решение этого уравнения с граничными условиями (34) так же известно [11]

                                                    (40)

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

.                                                                         (41)

Из соотношений (37), (41) параметры модели вычисляются

.                                             (42)

Теплофизические характеристики чистого никеля, расплав которого допускает сильное переохлаждение [49], достаточно хорошо известны [50]

Из требования термодинамической устойчивости [20]  и параметры модели для Ni

§2. Постановка задачи для плоского фронта кристаллизации.

В одномерном случае (плоский фронт) кристаллизации переохлаждённого расплава в неограниченной области  уравнения модели (29), (30) в безразмерных переменных  можно записать как

                                                (43)

                                                           (44)

где .

При невысокой степени переохлаждения  для уравнений (43), (44) ставится задача Коши. Начальные условия для температуры выбираются в соответствии с автомодельным решением (6), (7) задачи Стефана

                     (45)

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

                                   (46)

В этом случае можно ожидать, что решение задачи (43) – (46) будет близко к решению задачи Коши для модели Стефана (1) – (3), и при  выходить на автомодельное решение (6), (7) классической задачи Коши-Стефана кристаллизации переохлаждённого расплава.

Постановка задачи об установлении квазистационарного режима кристаллизации гиперохлаждённого расплава помимо уравнений (43), (44) включает граничные условия для температуры и параметра порядка

,                                   (47)

.                                       (48)

Результатом решения задачи (43) – (44), (47), (48) является определение квазистационарных распределений температуры , параметра порядка  и скорости движения фазового перехода . Первоначально, распределение  температуры задается в соответствии с квазистационарным решением задачи Стефана (8)

                                        (49)

а распределение поля параметра порядка предполагается  равновесным (35)

.                                                    (50)

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

§3. Разностная схема и алгоритм адаптации сетки к решению.

Основная трудность, связанная с численным решением уравнений (43), (44), заключается в том, что ширина межфазной области, где параметр порядка заметно меняется, может быть много меньше длины теплового поля. Для аппроксимации пространственных производных поля параметра порядка в этой области с необходимой точностью, шаг сетки  должен быть достаточно малым. Поскольку межфазная область является подвижной, шаги пространственной сетки должны изменятся с течением времени, адаптируясь к этой особенности решения задачи [52]. Эффективным алгоритмом построения адаптивных сеток для нестационарных задач с большими градиентами является метод динамической адаптации [51].

В основе метода лежит переход от исходной системы координат  физического пространства к нестационарной системе координат  расчетного пространства, в котором область больших градиентов решения неподвижна. Полагая, что соотношение перехода от физического пространства к расчётному и что существует обратное невырожденное преобразование , для частных производных справедливы соотношения [51]

                                  (51)

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

.                                                                            (52)

Уравнения модели (43), (44) для энтальпии  и параметра порядка  и граничные условия (47), (48) в новых независимых переменных  расчетного пространства будут

,  (53)

,                                                       (54)

           (55)

,                 (56)

 – координаты границ расчетной области

.                                       (57)

Дополнительное уравнение для функции скорости узлов , которое определяет механизм адаптации сетки к решению и замыкает систему уравнений (53) – (55), определяется из требования квазистационарности энтальпии в расчетном пространстве. Нестационарная система координат  выбирается таким образом, чтобы . Из уравнения (53) следует требуемое уравнение для функции

.    (58)

Для конечно-разностной аппроксимации используется дивергентная форма уравне-ний (53), (54)

,                                     (59)

,                                                  (60)

,                                                                 (61)

.

Сетка в расчетном пространстве выбирается равномерной по пространственной координате ,– число шагов сетки,   – координаты целых и полуцелых узлов. Шаг интегрирования по времени  может меняться в процессе вычислений . К целым узлам относятся функции, к полуцелым узлам относятся     .

Для конечно-разностной аппроксимации уравнений модели (58) – (61) была выбрана схема с центральными разностями второго порядка точности [53]

   (62)

            (63)

,                            (64)

,  (65)

где – регуляризирующий параметр [51],

.

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

Система (62) – (65) разностных нелинейных алгебраических уравнений решалась методом Ньютона. На каждой итерации лианеризованная система уравнений (62) – (65) “расщеплялась” по физическим процессам. В начале, с помощью прогонки, решалась “энтальпийная” часть задачи, т.е. уравнение (62) при заданном поле параметра порядка. Затем, так же прогонкой, решалось “релаксационное” уравнение (63) параметра порядка при заданной энтальпии. Итерация заканчивается вычислением значений функций  из уравнений (64), (65). Такая процедура решения линеаризованных уравнений с помощью раздельных прогонок оказывается более экономичной, чем с использованием матричной прогонки. Шаг интегрирования по времени выбирался автоматически, изменяясь в зависимости от количества итераций.

§4. Результаты тестовых расчетов.

При невысокой скорости движения фронта кристаллизации , длина теплового поля для типичных металлов много больше ширины фазового перехода, т.е.  . В этом случае, использование модели Стефана для описания процесса кристаллизации является вполне оправданным и результаты расчета динамики фазового перехода на основе модели фазового поля должны быть близки к решению задачи Стефана (1) – (3). Это означает, что ширина фазового перехода (ширина профиля параметра порядка) должна слабо меняться со временем, а скорость его движения и температурное поле вне области фазового перехода должны совпадать, с достаточно высокой точностью, со скоростью движения межфазной поверхности и распределением температуры задачи Стефана. Для проверки модели фазового поля (43), (44) и алгоритма динамической адаптации, выбранного для её численной реализации, было проведено два тестовых расчета процесса кристаллизации (плоский фронт) расплава разной степени переохлаждения  . Результаты моделирования сравниваются с автомодельным (6), (7) и квазистационарным (8) решениями задачи Стефана.

На рис. 2 – 4 представлены результаты расчета процесса кристаллизации переохлаждённого расплава  . Начальные условия для параметра порядка и температуры задавались в виде (45, 46)). Были выбраны следующие параметры, определяющие модель и начальные условия, . На рис. 2 изображены графики зависимостей параметра порядка , энтальпии , и температуры  от координаты  в двух разных масштабах (а, б) для момента времени , а на рис. 3; 4 представлены графики тех же функций в более поздние моменты времени. Маркерами на рисунках обозначены положения узлов адаптивной сетки. Число узлов сетки  в расчётах не превышало 80.

Из приведенных результатов видно, что профиль параметра порядка практически не меняется с течением времени, т.е. его движение носит волновой характер, а температура слабо меняется по ширине фазового перехода. Сравнительный анализ полученных результатов показал, что скорость движения фазового перехода и температурное поле вне фазового перехода отличаются от автомодельного решения задачи Стефана (5), (6) не более чем на 2% при  . С уменьшением величины параметра  отклонение получаемого решения от автомодельного сокращается.

На рис. 5 представлены результаты расчета процесса кристаллизации гиперохлажденного расплав . Граничные и начальные условия для поля параметра порядка и температуры задавались в виде (47) – (50) и определялись параметрами    . Как показали расчеты, решение выходит на квазистационарный режим, т.е. представляет собой волну кристаллизации, движущуюся с постоянной скоростью . Графики зависимостей параметра порядка , температуры  и энтальпии от переменной , после установления квазистационарного режима кристаллизации, представлены на рис. 5 в различном координатном масштабе (а, б). Скорость движения фазового перехода и распределение температуры вне фазового перехода совпадают с квазистационарным решением задачи Стефана с точностью до 2%.

                                              а                                                                                               б

    Рис. 2

                                             а                                                                                             б

   Рис. 3

                                                а                                                                                           б

      Рис. 4

                                                а                                                                                          б

       Рис.  5

§5. Анализ квазистационарного режима кристаллизации.

Задача одномерной кристаллизации переохлаждённого расплава (43), (44) с граничными условиями (47), (48) может иметь квазистационарное решение [20, 24, 25], т.е. фазовый переход движется с постоянной скоростью . Причём, квазистационарный режим кристаллизации реализуется не только при гиперохлаждении , как в модели Стефана (1) – (3), (8), но так же при .

Уравнения модели (43), (44), в случае квазистационарного режима кристаллизации, заменой переменных , , преобразуются к системе обыкновенных дифференциальных уравнений

                                           (66)

                                              (67)

Первый интеграл уравнения (66)

                                             (68)

с учётом граничных условий

выражает закон сохранения энтальпии между гомогенными состояниями вдали от фазового перехода

и температура твёрдой фазы  определяется температурой расплава

.                                                     (69)

Из общего решения уравнения (68)

производную температуры можно представить в интегральном виде

.                                 (70)

При условие  из (70) следует

      (71)

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

.                              (72)

Принимая во внимание, что

соотношение (72) можно преобразовать

.                                    (73)

Для равновесного распределения параметра порядка (35)

              (74)

и в случае  из (71) – (73) следует с точностью до членов  малости

                                                    (75)

При  уравнение (75) упрощается

.                                                     (76)

Из (76) следует, что квазистационарное решение  существует только при . В этом случае, скорость движения фронта кристаллизации совпадает со скоростью движения межфазной поверхности  квазистационарного решения задачи Стефана кристаллизации гиперохлаждённого расплава (8).

            Только одно решение уравнения (75)

                   (77)

удовлетворяет условию кинетической устойчивости  [20]. Из соотношения (77) следует возможность существования квазистационарного решения  при , если число  меньше критической величины . В этом случае, образующаяся твёрдая фаза находится в перегретом (метастабильном) состоянии , где  и максимальная величина перегрева твёрдой фазы .

            На основе модели фазового поля (43), (44) были выполнены расчёты методом динамической адаптации процесса установления квазистационарного режима кристаллизации переохлаждённого расплава для плоского фронта. Граничные и начальные условия задавались в виде (47) – (50). Модель определялась параметрами . Вычисления выполнялись до установления квазистационарного решения. Количество узлов адаптивной сетки не превышало 100. Результаты расчётов представлены на рис. 6 – 15.

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

Рис. 6

– область определения функции . Из результатов, представленных на рисунке, следует, что если величина числа  меньше критического значения , то возможно установление квазистационарного режима кристаллизации расплава с образованием перегретой твёрдой фазы . Когда число  превышает критическое значение , квазистационарный режим реализуется только при гиперохлаждении расплава .

Критическое значение  можно оценить из результатов расчётов, представленных на рис.7 (а, б), где изображены графики зависимости числа  от величины числа  при различной степени переохлаждения расплава .

                                                 а                                                                                         б

     Рис. 7

Из рис. 7 следует, что квазистационарное решение устанавливается при любых значениях числа  в случае гиперохлаждения расплава , а для  только при .   Графики зависимостей параметра порядка , температуры , энтальпии  от переменной и величина числа  после выхода на квазистационарный режим представлены на рис. 8 (а – е) для различных докритических значений числа  при переохлаждение расплава , т. е. режимы с образованием максимально-перегретой твёрдой фазой. Величина  в результатах, представленных на рисунке, определялась путём постепенного уменьшения величины  и проведением тестовых расчетов процесса установления квазистационарного режима. Следует отметить, что по мере приближения величины переохлаждения расплава  к  задаваемые начальные условия задачи должны быть достаточно близки к искомому квазистационарному решению. В противном случае, квазистационарный режим не устанавливается с течением времени. Скорость фазового перехода и величина перегрева твёрдой фазы монотонно уменьшаются, длина теплового поля возрастает и решение, в этом случае, выходит на автомодельное решение задачи Стефана. Расчётные значения критической величины числа  и максимальной величины перегрева твёрдой фазы отличаются от теоретических оценок, основанных на соотношении (77). Более точные теоретические значения этих параметров были получены [20] с помощью бифуркационного анализа системы обыкновенных дифференциальных уравнений (66), (67), а для случая  , на основе приближённого решения уравнения (73) с точностью до членов  малости [31].

            На рис. 9 – 15 представлены графики зависимостей параметра порядка , температуры  и энтальпии  от переменной , а так же величина числа Пекле , квазистационарного решения задачи кристаллизации расплава разной степени переохлаждения  для нескольких значений числа .

       а                                                                                          б

   в                                                                                           г

      д                                                                                         е

Рис. 8

  а                                                                                         б

   в                                                                                           г

  д                                                                                       е

Рис. 9

  а                                                                                          б

в                                                                                        г

д                                                                                         е

Рис. 10

а                                                                                         б

в                                                                                       г

д                                                                                         е

Рис. 11

а                                                                                       б

в                                                                                        г

д                                                                                        е

Рис. 12

а                                                                                         б

в                                                                                        г

д                                                                                         е

Рис. 13

  а                                                                                          б

  в                                                                                          г

д                                                                                         е

Рис. 14

а                                                                                        б

 в                                                                                          г

д                                                                                         е

Рис. 15

Из приведённых результатов следует, что с увеличением величины переохлаждения расплава возрастает скорость квазистационарного движения фазового перехода, длина теплового поля значительно сокращается, а распределение параметра порядка изменяется слабо, оставаясь практически равновесным.

Заключение.

            Проведённое численное исследование процесса  кристаллизации метастабильного расплава на основе модели фазового поля  методом динамической адаптации показало:

1.            результаты моделирования с достаточно высокой точностью согласуются с решениями задачи Стефана в условиях, когда ширина фазового перехода много меньше длины теплового поля, т.е. ;

2.            возможность существования квазистационарного режима кристаллизации переохлаждённого расплава с образованием метастабильной (перегретой) твёрдой фазы при докритической величине числа  – отношения коэффициента термодиффузии среды к коэффициенту “релаксационной” диффузии.

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

Литература.

1.            Мейерманов А.М. Задача Стефана. – Новосибирск: Наука, 1986, 240 с.

2.            Langer J.S. Instabilities and pattern formation in crystal growth. – Rev. Mod. Phys. 1980, v. 52, pp. 1 – 28.

3.            Kessler D.A., Koplik J., Levine H. Pattern selection in fingered growth phenomena. – Adv. Phys. 1988, v. 37, pp. 255 – 339.

4.            Уманцев А.Р.  Движение плоского фронта при кристаллизации. – Кристаллография. 1985, т. 30, с. 153 – 160.

5.            Борисов В.Т. О механизме нормального роста кристаллов. – ДАН СССР. 1963, т. 151,     с. 1311 – 1314.

6.            Coriell S.R., Parker R.L. Interface kinetics and the shape of a solid sphere growing from the melt. – In “Crystal growth”. Ed. Peiser H.S. Pergamon, Oxford, 1967, pp. 703 – 708.

7.            Glicksman M.E., Shaefer R.J. Investigation of solidification interfaces temperatures via isenthalpic solidification. – J. Cryst. Grow. 1967, v. 1, pp. 297 – 310.

8.            Shaefer R.J., Glicksman M.E. Fully time-dependent theory for the growth of spherical crystal nuclei. – J. Cryst. Grow. 1969, v. 5, pp. 44 – 58.

9.            Mikheev L.M., Chernove A.A. Mobility of a diffuse simple crystal-melt interface – J. Cryst. Grow. 1991, v. 112, pp. 591 – 596.

10.    Langer J.S. Model of pattern formation in first-order phase transition. – In “Directions in condense matter physics”. Eds. Grinstein G., Mazenco G. Word Sc. Singapore, 1986, pp. 164 – 186.

11.    Chan S. -K. Steady-state kinetics of diffusionless first order phase transitions. – J.  Chem. Phys. 1977, v.67, pp. 6755 – 6762.

12.    Gunton J.D., San Miguel M., Sahni P.S. The dynamics of first-order phase transition. – In “Phase transition and critical phenomena”. Eds. Domb C., Lebowitz J.L. Academic Press, London, 1983, v. 8, pp. 267 – 482.

13.    Fix G. Phase field models for free boundary problems. – In “Free boundary problems: theory and applications”. Eds. Fasano A., Primicerio M. Pittman, London, 1983, pp. 580 – 589.

14.    Caginalp G. Surface tension and supercooling in solidification theory. – In “Application of field theory to statistical mechanics”. Lecture notes in physics, № 216. Ed. Garrido L. Springer, Berlin, 1984, pp. 216 – 226.

15.    Collins J.B., Levine H.I. Diffuse interface model of diffusion-limited crystal growth. – Phys. Rev. B. 1985, v. 31, pp. 6119 – 6122.

16.    Collins J.B., Chakrabari A., Gunton J. D. Dynamics of phase separation in a model for diffusion-limited crystal growth. – Phys. Rev. B. 1989, v. 39, pp. 1506 – 1511.

17.    Caginalp G. An analysis of phase field model of a free boundary. – Arch. Ration. Mech. Anal. 1986, v. 92, pp. 205 – 245.

18.    Caginalp G., Jones J. A derivation and analysis of phase field models of thermal alloys. –Annal. Phys. 1995, v.237, pp. 66 – 107.

19.    Уманцев А.Р., Ройтбурд А.Л. Неизотермическая релаксация в нелокальной среде. –  ФТТ. 1988, т. 30, с. 1124 – 1131.

20.    Umantsev A. Thermodynamic stability of phases and transition kinetics under adiabatic conditions. – J. Chem. Phys. 1992, v. 96, pp. 605 – 617.

21.    Penrose O., Fife P.C. Thermodynamically consistent models of phase-field type for the kinetic of phase transitions - Physica D. 1990, v. 43, pp. 44 – 62.

22.    Penrose O., Fife P. On the relation between the standard phase-field model and a “thermodynamically consistent” phase-field model. – Physica D. 1993, v. 69, pp. 107 – 113.

23.    Fife P.C., Penrose O. Interfacial dynamics for thermodynamically consistent phase-field models with nonconserved order parameter. – Electron. J. Diff. Equ. 1995, v. 1995, pp. 1 – 49.

24.    Bates P.W., Fife P.C., Gardner R.A., Jones C.K.R.T. Phase field model for hypercooled solidification. – Physica D. 1997, v. 104, pp. 1 – 31.

25.    Bates P.W., Fife P.C., Gardner R.A., Jones C.K.R.T. The existence of travelling wave solution of a generalized phase-field model. – SIAM J. Math. Anal. 1997, v. 96, pp. 60 – 93.

26.    Wang S. -L., Sekerka R.F., Wheeler A.A., Murray B.T., Coriell S.R., Braun R.J., McFadden G.B. Thermodynamically consistent phase-field models for solidification. – Physica D. 1993,    v. 69, pp. 189 – 200.

27.    Braun R.J., McFadden G.B., Coriell S.R. Morphological instability in phase-field models of solidification. – Phys. Rev. E. 1994, v. 49, pp. 4336 – 4352.

28.    Wang S. -L., Sekerka R. Computational of the dendritic operating state at large supercoolings by the phase-field model. – Phys. Rev. E. 1996, v. 53, 3750 – 3776.

29.    Schofield S.A., Oxtoby D.W. Diffusion disallowed crystal growth. I. Landau-Ginzburg model. – J. Chem. Phys. 1991, v. 94, pp. 2176 – 2186.

30.    Löwen H., Schofield S.A., Oxtoby D.W. Diffusion disallowed crystal growth. II. A parabolic model. – J. Chem. Phys. 1991,v. 94, pp. 5685 – 5692.

31.    Löwen H., Bechhofer J. Critical behavior of crystal growth velocity. – Europhys. Lett. 1991,v.16, pp. 195 – 200.

32.    Löwen H., Bechhoefer J., Tuckerman L. Crystal growth at long time: Critical behavior at the crossover from diffusion to kinetic-limited regimes. – Phys. Rev. A. 1992, v.45, pp. 2399 – 2415.

33.    Karma A., Rappel W. -J. Quantitative phase-field modeling of dendritic growth in two and three dimensions. – Phys. Rev. E. 1998, v. 57, pp. 4323 – 4349.

34.    Fabbri M., Voller V.R. The phase-field method in the sharp-interface limit: a comparison between model potentials. – J. Comput. Phys. 1997, v. 130, pp. 256 – 265.

35.    Мандельштам Л.И. , Леонтович М.А. К теории поглощения звука в жидкостях. –ЖЭТФ. 1937, т. 7, с. 438 – 449.

36.    Леонтович М.А. Введение в термодинамику. Статистическая физика. – М. : Наука, 1983, 416с.

37.    Ландау Л. К теории фазовых переходов. – ЖЭТФ. 1937, т. 7, с. 19 – 32.

38.    Ландау Л.Д. К теории фазовых переходов. II – ЖЭТФ. 1937, т. 7, с. 627 – 632.

39.    Гинзбург В.Л., Ландау Л.Д. К теории сверхпроводимости. – ЖЭТФ. 1950, т. 20, с. 1064 – 1082.

40.    Ландау Л.Д., Халатников И.М. Об аномальном поглощении звука вблизи точек фазового перехода второго рода. – ДАН СССР. 1954, т. 96, с. 469 – 473.

41.    Cahn J.W., Hilliard J.E. Free energy of nonuniform system. I. Interfacial free energy. – J. Chem. Phys. 1958, v. 28, pp. 258 – 267.

42.    Cahn J.W., Hilliard J.E. Free energy of nonuniform system. III. Nucleation in a two-component incompressible fluid. – J. Chem. Phys. 1959, v. 31, pp. 688 – 699.

43.    Cahn J.W. Theory of crystal growth and interface motion in crystalline materials. – Acta Metall. 1960, v. 8, pp. 554 – 562.

44.    Cahn J.W. On spinodal decomposition. – Acta Metall. 1961, v. 9, pp. 795 – 801.

45.    Уайт Р., Джабелл Т. Дальний порядок в твердых телах. – М.: “Мир”, 1982, 448с.

46.    Halperin B.I., Hohenberg P.C., Ma S.-K. Renormalization group method for critical dynamics. I. Recursion relations and effects of energy conservation. – Phys. Rev. B. 1974, v.10, pp.139 – 153.

47.    Hohenberg P.C., Halperin B.I. Theory of dynamic critical phenomena. – Rev. Mod. Phys. 1977, v. 49, pp. 435 – 479.

48.    Metiu H., Kitahara K., Ross J. Statistical mechanical theory of the kinetics of phase transitions. - In “Fluctuation phenomena. Studies in statistical mechanics”. Eds. Montroll E.W., Lebowitz J.L. North-Holland Publish. Comp. , Amsterdam, 1979, v. 7, pp. 229 – 291.

49.    Bassler B.T., Hofmeister W.H., Carro G., Bayuzick J. The velocity of solidification of highly undercooled nickel. – Metall. Mater. Trans. 1994, v. 25A, pp. 1301 – 1308.

50.    Desai P.D. Thermodynamic properties of nickel. – Int. J. Thermophys. 1987, v.8, pp. 763 – 780.

51.    Мажукин В.И., Самарский А.А., Кастельянос О., Шапранов А.В. Метод динамической адаптации для нестационарных задач с большими градиентами. – Мат. Моделирование. 1993, т. 5, с. 33 – 56.

52.    McCarthy J.F. One-dimensional phase field models with adaptive grids. – J. Heat Trans.  1998, v. 120, pp. 956 – 964.

53.    Самарский А.А. Теория разностных схем. – М. : Наука, 1989, 616 с.