Численный расчет одной модели, описывающей кристаллизацию металлов I. Одномерный случай

( The numerical calculation of some model, which describes the metals’ crystallization process I. 1-D case
Preprint, Inst. Appl. Math., the Russian Academy of Science)

Авторы выражают благодарность Е.В.Радкевичу и Н.Н.Яковлеву, познакомивших авторов с моделью, за интенсивные и плодотворные обсуждения.

Введение

Целью настоящего препринта является численное исследование недавно появившейся в литературе [5, Гл.7] модели, описывающей процесс кристаллизации в эвтектических сплавах, не образующих твердых растворов. Эта модель является расширением широко известной модели Кана-Хилларда [3], [4] и ее модификаций, см., например, [2].

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

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

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

§1. Постановка задачи о кристаллизации и математическая формулировка модели

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

 — мольная концентрация кристаллизующегося металла в жидкой фазе,

 — усадка, т.е. изменение объема при изменении состава жидкой фазы и фазовом переходе,

 — скорость конвекции в жидкой фазе,

 — осреднённое смещение твёрдой фазы,

 — средняя скорость нарастания твёрдой фазы (осреднённая скорость микрофронтов),

 — осреднённое условное смещение жидкости относительно твёрдой фазы,

 — температура.

За основу модели в [5] принята макромодель, которая опирается на модель Био насыщенной пористой среды [1] и закон сохранения массы жидкой фазы, а также микромодель ликвации Канна-Хилларда [3], [4], которая расширена учетом конвекции и теплопереноса. Таким образом, при моделировании необходимо учитывать несколько масштабных уровней.

          Обозначим . Тогда система одномерных уравнений, моделирующих процесс кристаллизации, может быть записана в следующем виде:

                      (1.1)

Система (1.1) записана как объединение трех подсистем, выделенных квадратными скобками. Эти подсистемы представляют собой описания основных физических процессов, участвующих в создании модели и о которых мы говорили выше. Первая подсистема представляет собой модель распространения волн в пористом «каркасе», заполненным жидкой фазой, – упрощенный вариант так называемой модели Био. Вторая подсистема, состоящая из одного уравнения, представляет собой просто уравнение неразрывности для жидкой фазы и управляет эволюцией усадки , поскольку является функцией и  (см. ниже). Третья подсистема является расширенной моделью Кана-Хилларда образования и роста зародышей, включающей конвективный перенос и теплопередачу.

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

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

Дальнейшие пояснения и интерпретации указанных констант и функций см. в [5].

Также скажем несколько слов относительно краевых условий, используемых нами при решении системы (1.1):

Функции и .

Начальные условия:

 , что соответствует ;

краевые условия:

, что соответствует .

Усадка .

Начально-краевые условия:

.

Концентрация .

; . Константы здесь необходимо выбирать в соответствии с начальными данными для , чтобы обеспечить непрерывность. Начальные данные для  выбираются, исходя из конкретной физической ситуации.

Температура .

; ; ,

где , , является некоторым параметром скорости остывания.

§2. Численная схема для решения системы уравнений, моделирующей процесс кристаллизации

Как можно видеть, система (1.1) содержит подсистемы разного характера: первая подсистема, вообще говоря, представляет собой систему уравнений первого порядка (в области допустимости модели эта подсистема должна быть гиперболической); вторая подсистема – это уравнение переноса; а третья подсистема является подсистемой уравнений высокого порядка (первое уравнение при наших предположениях имеет четвертый порядок, а второе уравнение имеет второй порядок). Математическая форма системы (1.1) подсказывает, что наиболее удобным способом решения будет метод расщепления, когда «гиперболическая» и «переносная» подсистемы решаются с помощью использования явных схем (учитывая возможное появление разрывов), а «диффузионно-дисперсионная» подсистема решается с помощью использования неявной схемы (учитывая возможное расширение методики на многомерный случай). При этом нелинейный конвективный член в уравнении Кана-Хилларда трактуется в рамках технологии TVD-схем, принимая во внимание возможность распространения разрывов концентрации, которые в некоторых режимах могут не сглаживаться «диффузионно-дисперсионными» членами.

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

(2.1)

где

(2.2)

Рассмотрим гиперболическую часть системы (2.1) (а значит первую подсистему системы (1.1)) более подробно. Введем обозначения , . Тогда рассматриваемая гиперболическая часть (без правых частей) будет выглядеть следующим образом

                                  ,                                    (2.3)

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

                  .                     (2.4)

Обозначив эти собственные числа через , мы легко найдем соответствующие собственные направления

           .

Заметим при этом, что наша модель будет иметь смысл только если все четыре собственных значения будут действительными, в противном случае мы будем считать, что модель выходит за рамки применимости.

Итак, напомним, что для решения системы (2.1) (а значит и для (1.1)) мы будем использовать метод расщепления. Для первых двух подсистем системы (1.1) возьмем явную схему, которая будет выглядеть так

(2.5)

Здесь  и  представляют собой шаги по времени и пространству соответственно, а индексы  и  относятся к величинам в -ый момент времени в -ой точке расчетной сетки по пространству. Схема (2.5) является схемой первого порядка, реализующей принципы TVD-схем для последнего уравнения переноса. Естественно она дополняется соответствующими краевыми условиями.

Стандартными методами можно показать, что (2.5) будет устойчивой при выполнении следующих ограничений на временной шаг

                   ,                     (2.6)

где  являются корнями уравнения (2.4).

Разностная схема для последнего уравнения системы (2.1) будет выглядеть следующим образом (неявная схема, в предположении, что  уже найдено)

                  .                    (2.7)

Получающаяся система линейных уравнений для определения температуры на верхнем слое решается стандартным методом прогонки с учетом поставленных начально-краевых условий.

Уравнение Кана-Хилларда для эволюции концентрации  имеет наиболее сложную структуру, и должно быть рассмотрено более подробно. Напомним, как выглядит уравнение Кана-Хилларда (см. (1.1) или (2.1))

.   (2.8)

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

. (2.9)

Производя необходимые дифференцирования и группируя члены при различных производных от  получим из (2.9)

  . (2.10)

          Заметим теперь, что переносный член в (2.10) для правильного расчета распространения разрыва с точки зрения теории TVD схем необходимо записывать в дивергентной форме, т.е. (2.10) окончательно преобразуется к виду

                  (2.11)

где

Для уравнения (2.11) будем использовать следующую разностную схему

  .  (2.12)

Соберем коэффициенты при неизвестных на новом слое

 ,  

где

а

можно считать решением по явной TVD-схеме уравнения

                                 .

После перегруппировки система уравнений (2.12) может быть решена стандартным методом прогонки.

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

Как уже говорилось выше, представленная модель (даже в одномерном случае) с математической и вычислительной точек зрения является конгломератом сложных нелинейных и разномасштабных процессов. Поэтому представляет определенный интерес просто проверить работоспособность модели и возможность производить расчеты по представленной выше численной методологии. (Это является особенно важным при переходе к моделированию в двух- и трехмерном случаях, когда физическая интерпретация полученных результатов приобретет более осмысленный характер.) Поэтому далее мы приведем демонстрационный расчет, где будем следить только за математической эволюцией основной неизвестной функции – концентрации – без детальной привязки к какому-либо физическому опыту. Единственная физическая интерпретация, которую мы будем использовать, состоит в следующем: если выделить некоторое промежуточное значение концентрации , которое мы обозначим как «критическую» концентрацию , то будем считать, что кристаллизующийся металл находится в жидком состоянии при и в твердом состоянии при .

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

Рисунок 1. Начальные функции температуры, концентрации и плотности жидкой фазы для демонстрационного расчета.

Рисунок 2 демонстрирует эволюцию во времени функции концентрации. Для каждого момента времени приведены последовательно две картинки: первая – это расчет по явной схеме, вторая – расчет по неявной схеме. Можно заметить, что обе картинки совпадают с высокой точностью. При этом величина шага по времени для этой сетки, в согласии с теоретическими представлениями, увеличивается для неявной схемы примерно в 70-80 раз, что обуславливает затрату значительно меньшего времени на один расчет.

          Рисунок 3 демонстрирует вид других функций модели, которые расчитаны по неявной схеме на момент .

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

Рисунок 2. Эволюция во времени функции концентрации, расчитанная по явной и неявной схемам .

Рисунок 3. Другие функции модели, расчитанные по неявной схеме к моменту времени .

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