Двухкомпонентный параметр порядка в модели фазового поля для описания кристаллизации расплава Текст научной статьи по специальности «Физика»

УДК 669.017

П.В. Меленев, Н.Д. Няшина, П.В. Трусов Пермский государственный технический университет

ДВУХКОМПОНЕНТНЫЙ ПАРАМЕТР ПОРЯДКА В МОДЕЛИ ФАЗОВОГО ПОЛЯ ДЛЯ ОПИСАНИЯ КРИСТАЛЛИЗАЦИИ

РАСПЛАВА

Abstract

The two-component order parameter is introduced to derive constitutive relations based on the principle of entropy production positiveness. The application of this principle instead of the free energy one (as in isothermal case) allows obtaining thermodynamically valid evolutionary equations for the components of the order parameter and thermal field and taking into account the nonisothermal conditions of solidification. The analysis of the simplified version of the equations enables displaying the physical sense of the parameters in the evolutionary equations. The test calculations for the one-dimensional problem show that the model describes qualitatively physics of solidification and grain formation.

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

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

Постановка задачи

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

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

должна эволюционировать при условии положительности производства энтропии (в том числе и в неизотермических условиях).

Воспользуемся для вывода уравнений модели формализмом, предложенным в работах [2, 3].

Рассмотрим закрытую область й, заполненную чистым металлом, способным претерпевать фазовый переход «жидкость - кристалл». В последующих выводах рассматривается двухмерный вариант задачи (используется декартова система координат с ортонормированным базисом {/, у'}). Введем характеристику фазового состояния материала - параметр порядка (1111), состоящий из двух переменных: {(,г), 9(,г)} [3]. Переменную ф можно интерпретировать как степень

упорядоченности материала (доля заполненных узлов кристаллической решетки) в точке (точнее, в микрообъеме) с радиус-вектором г в момент времени 7; ф = 0 соответствует жидкому состоянию, а ф = 1 - кристаллическому. Переменная 9 определяет угол между одним из главных кристаллографических направлений и осью абсцисс в выбранной системе координат; очевидно, что, учитывая порядок оси симметрии Л0 рассматриваемого типа решетки, можно рассматривать 9 в пределах от 0 до 2п/Л0 (для двухмерного случая).

Предположим, что для энтропии произвольного подобъема рассматриваемой области V еП можно записать:

5 = |{<е,<р)--2(г|Уф|)2 -1 (пИ)2, (1)

V

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

Для учета анизотропии поверхностной энергией кристалла предположим, что коэффициент 8 является функцией ориентации поверхности кристалла, которую будем описывать углом у = у-9, где у - угол между осью абсцисс и вектором -Уф, направленным, как предполагается, по нормали к границе раздела фаз (см. рис. 1).

Третье слагаемое в (1) отражает влияние разориентации соседних кристаллов. Предположим, что влияние это зависит от упорядоченности материала, тогда п = п(ф).

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

е + V • Ц = 0, (2)

тогда, выделяя выражение для потока энтропии через границу подобъема V, изменение энтропии по времени можно записать в следующем виде:

5=|(ц Ч т (¥1 + хф]ф+х99]л' -Iп-1® ц+

V

где обозначено

дф

де

зу о уф

х 9 = 0(уМу)|уф| 2 + У • (п(ф)у9);

У =ф [г(у)2 Уф + г(у)г'(у) J • Уф]+9п(ф)У9. Здесь тензор I = И + у , тензор J - единичный.

Рис. 1. Схематичное изображение границы раздела фаз

Согласно второму закону термодинамики производство энтропии в любом V еП положительно. Производство энтропии можно найти, если вычесть из 5 поток энтропии через границу (V :

[Ц ' т

В подынтегральном выражении

£ +/п+ г\с1 (ЗУ)> 0. ч

ЗУ

Т

поток энтропии за счет теплового потока,

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

= Нч Ч Т] +

(де ^ .дф.

ф + Xе 9 \ё¥ > 0

5

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

Ц = Мт У , (3)

Т V

1 (де

тф = - т|д-| +V

Т1д^£

- [є(М/)(є(М/)^ + є'(ц/)./)- V-]-n^(ф)|vе2, (4)

т9 = 0(у)0,(у)уф 2 + У • (n(ф)v9), (5)

где т, Мт - положительные параметры. Это наиболее простой, но не единственно возможный выбор вида соотношений (3) - (5). Воспользовавшись (2), можно записать уравнение для изменения плотности внутренней энергии:

1

Є = -V -

МТ VI

Т

(6)

Уравнения (4)-(6) представляют собой систему управляющих уравнений для параметра порядка и плотности внутренней энергии соответственно.

Для определения производной от плотности внутренней энергии в (4) целесообразно использовать плотность свободной энергии Гельмгольца [2].

Параметр Мт можно задать в виде Мт =ХТ2, что обеспечит пропорциональность теплового потока градиенту температур. Тогда коэффициент X имеет смысл теплопроводности. Зададим вид приградиентных функций:

г(У) = 0 00Ы> П(ф) = П0 П(ф),

где 00, п0 - положительные постоянные. Опуская подробные выкладки, запишем эволюционные уравнения модели:

тф = 02у • [0(у)(0(у)~ + 0'(у) ~)-уф]-П0 2 - д(т)р,(ф) - 00g,(ф); (7)

т9 = 0 0 0(у)0'(у)уф 2 + П0у • (п(ф)у9); (8)

(с - р(ф)Ь'(т))т - р'(ф)Щ)ф = ГУ2т. (9)

где функции Ь(т), р(ф), 00g(ф) появились при выводе уравнения для поля температуры в предположении, что плотность внутренней энергии имеет вид

е =е(т) +[1 - р(ф)] Щ) = е1(т) - р(ф) Щ),

где е,, е1 - классические плотности внутренней энергии твердой и жидкой фаз соответственно, а ^ )= е, (т)- е, (т), причем 1>(тм,) = 10 - скрытая теплота плавления; функция р(ф) определяет зависимость плотности внутренней энергии от упорядоченности материала, то есть характер зависимости е(ф, т) в двухфазной области, где 0 < ф < 1, и в соответствии с выбранными значениями ф для объемных фаз: р(0) = 0, р(1) = 1.

Функция Q(T) в соотношении (7) связана с частью свободной энергии, учитывающей выделение теплоты кристаллизации; она появляется после интегрирования по температуре при переходе от плотности внутренней энергии к плотности свободной энергии и имеет вид

Q(т )= 1Ж.

тм

Вид функции 0(ф) выбирается таким образом, чтобы потенциал свободной энергии имел минимум при ф=1 и ф=0; для простоты берется аппроксимация этой функции полиномом 4-го порядка:

0(ф) = 00g(ф) = °0ф (1 - ф)2.

Функция р(ф) удовлетворяет указанным выше условиям нормировки, вытекающим из смысла представления плотности свободной энергии (р(0) = 0, р(1) = 1 ) и имеет вид

р(ф) = |g(^)^и |g(С)^С = ф3 (10 - 15ф + 6ф2).

0 / 0

Полученная система (7)-(9) содержит неизвестные константы: 80, п0, т, 00.

Физический смысл параметров, входящих в уравнения

Рассмотрим одномерную изотропную формулировку задачи (7)-(9) при равновесных условиях: т = тм , ф = 0 . Предположим, что ориентация всего объема

одинакова - 9 = 90, неизменна и ориентация границы раздела у = у0. Тогда, рассматривая направление £, по углу 90, получим модификацию уравнения (7):

8 2 (У 0 - ° £ ,(ф) = °. (10)

д-

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

ф ^+00 > 0, ф ^-о >1. (11)

Решение (10)-(11) имеет вид

ф(-) = -1

ік

О,

рв(у 0 )'

+1

(12)

Если оценить размер области 5 существенного изменения ф (в диапазоне значений от 0,05 до 0,95), то выяснится, что (учитывая произвольность выбора ориентации кристалла)

3^2"0 0

5

л/О0

(13)

Величину 5 можно интерпретировать как ширину размытой границы раздела

фаз.

Далее, найдем свободную поверхностную энергию (энергию вновь образующейся размытой границы раздела фаз), соответствующую полученному решению (12). Для этого рассмотрим потенциал свободной энергии Гельмгольца при температуре плавления, учитывая выбранный вид энтропии (1),

^ = Е- тм5 = |{/(,ф)+ 2ти0-(уХУф)!^V.

Соответствующая поверхностная энергия на единицу площади, определяемая как величина отклонения свободной энергии от /(тм, ф) = 0 (предполагается, что это значение соответствует устойчивым фазам вещества), то есть как свободная энергия двухфазной области (где 0 < ф < 1), равна в одномерном случае:

а =

| Тм 8 2 (У 0

дф

В результате получаем соотношение, связывающее параметры модели 00 и 80 (через 5) с а,

а = ^°0тм8 . (14)

2

-о

Далее оценим параметр п0. Для этого воспользуемся оценкой времени поворота 1г кристаллического зародыша в поле вязких сил, предложенной в [4, 5]:

~ г3

, (15)

г кт

где ~ - динамическая вязкость, г - радиус зародыша.

Предположим, что микрокристалл встраивается (то есть при слиянии образуется кристалл без дефектов) в уже существующий кристалл (для которого ф=1 и Уф = 0). Тогда из уравнения (8) для этого случая

т9 = п0 V2 9

можно найти оценку для времени встраивания (если предположить, что радиус кластера соизмерим с размером элементарного объема):

кТ • (16)

По кТ

Надо отметить, что (16) является достаточно грубой оценкой.

Коэффициент т можно оценить, рассматривая гомогенное зародышеобразование при постоянном переохлаждении на начальной стадии кристаллизации. В этом случае доля закристаллизовавшегося вещества невелика и зародыши развиваются независимо. Зная размер критического зародыша г* и частоту зародышеобразования J (обе эти характеристики являются функциями переохлаждения АТ ), можно оценить изменение величины ф при образовании в объеме одного критического зародыша: ф* (АТ). С учетом этого эволюционное уравнение (7) примет вид

тф = -0(Ты - АТ)Р'(Ф*) - ^оЯ'(Ф*) .

В последнем выражении первое слагаемое больше нуля в силу вида функции Q(T)• Частоту зародышеобразования и скорость роста твердой фазы можно связать соотношением ф|ф =ф*/1* = АУ = Jv, где и, V - время ожидания и объем

критического зародыша, соответственно. Тогда т можно оценить из следующего соотношения:

т = - ^Тм - АТ)Р'(Ф*) - ^оЯ'(Ф*) (17)

Jv

Безразмерные уравнения

Для удобства расчетов запишем эволюционные уравнения (7)-(9) в безразмерной форме, используя следующие масштабы:

• координаты (Г - характерный размер рассматриваемой области);

• время (Г 2/х = (г2/ ^)рс - характерное время распространения температуры по области Г2 (х - коэффициент температуропроводности));

• безразмерная температура (и = (Т — Тм )/ АТ, где АТ - некоторый

характерный интервал температур (например, величина максимального переохлаждения для рассматриваемого металла)).

Используя эти масштабы и умножив уравнения (7)-(8) на Г2 /тх, а уравнение (9) на Г 2 Дат , запишем уравнения в безразмерном виде

|~- = • [в(ф)(в(ф)~ + в ' (ф) 0 )• V ф]- 0 ДМ|у0

2 ~

— Усиб(и)р (ф) -уО0g (ф); (18)

д0 = єє(ф)Б ' (ф)уф 2 + 0У • (п(ф)у 0); (19)

где

50 ~ , ч,, 2

д?

(1 - р(ф)&1'(и))-дд0 - р'(Ф)иДи)-^0 = у2и, (20)

дї д 1і

У = г /тх, и = Ь0/сАТ, ? = вІ/XX, 0 = По/ТХ, (21)

ь{и) = Ь(тм + АТи')/Ь0 , (22)

0(«)=} ( ,1(~) '2 <*. (23)

о (АТ + х)2

Таким образом, получена модель, описываемая уравнениями (7)-(9) (и их безразмерными аналогами (18)-(20) соответственно). Так как данная модель была построена на основе функционала энтропии и закона ее возрастания (справедливого и в неизотермических условиях, которые, как было отмечено ранее, могут иметь место при кристаллизации металла), то она согласуется с принципами термодинамики необратимых процессов.

Результаты решения одномерной задачи

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

= _0пН|1г) ~уш(~(и)Р,(Ф)-у°оя,(Ф); (24)

д г дх2 2 1дх)

^0 = 0 дї

,ґ \ дф 50 ґ \д2 0

п ^дГдг+^фЬгт

дх дх дх2

(25)

(1 - р(ф)иі'(и )) - р' (ф) иь(и) д?=-д-и, (2б)

д ї дї дх2

где х - безразмерная координата.

Предположим, что теплота фазового перехода слабо зависит от температуры: Ь (и ) = 1, следовательно,

е(« )=Т +). (27)

Тм \Тм /АТ + и)

Перепишем систему (24)-(26) с учетом сделанных допущений:

дф ? д2ф ,(ф)Гд0>12 АTug(ф)

^7 = вТТ~0Л21I I _усШ ( АТ + )-уСоg,(ф); (28)

дї дх 2 ) Тм (Тм /АТ + и)

*0 = 0

дї

\дф 50 і \д20

пИ^^+тфЬт

дх дх дх

(29)

| - 30?(ф)о§Ф = 0. (30)

дг дг дх

В тестовой задаче рассматривалась кристаллизация в толще расплава при постоянном переохлаждении -Аи на границах рассматриваемой области:

и(0, х) = и(г,-да) = и(г,+да) = -Аи. (31)

Для фазового поля в качестве начальных условий были выбраны случайно

распределенные возмущения с амплитудой 0,1 (угол ориентации 9 распределен случайным образом) и периодические граничные условия

ф(, - да) = ф(г, + да), (32)

9(г, -да) = 9(г, +да). (33)

Уравнения (28), (29) содержат нелинейные дифференциальные операторы, поэтому для решения использовалась явная схема с постоянными шагами по времени и координате. При этом использовались разностные аналоги производных, обеспечивающие аппроксимацию первого порядка по времени и второго порядка по координате. Для обеспечения устойчивости численного решения величина шага интегрирования по времени выбиралась в соответствии с условиями Куранта [6], при

этом величина шага по координате диктовалась размером области и количеством узлов в расчетной сетке.

ї = 0,05

і

1 1.4

0.8 1.2

1

0.6 А А 0.8

0.4 \ \ 0.6

1 \ 1 \ Г\ г\ Ґ 0.4

0.2 0.2

20

ї = 0,07 і

40

60

80

100

ї = 0,10 і

ї = 0,35 і

ї = 0,50 і

Рис. 2. Вид решения одномерной изотропной задачи для нескольких моментов времени

я

q

я

я

На рис. 2 показано распределение упорядоченности и угла ориентации для нескольких моментов времени, отражающие развитие процесса кристаллизации. Видно, что на начальном этапе (t = 0,05; 0,07) кристаллизация начинается в нескольких центрах, за счет роста которых, в основном, и происходит увеличение доли твердой фазы. В области существуют несколько участков закристаллизовавшегося вещества (ф близок к 1), разделенных менее упорядоченными промежутками. На более поздних стадиях процесса (t = 0,10; 0,35) кристаллы начинают взаимодействовать: прирастают один за счет другого, выравнивая свою ориентацию. Уменьшение количества кристаллов идет, во-первых, из-за расплавления малых кристаллов за счет скрытой теплоты кристаллизации, выделяющейся при росте больших; во-вторых, за счет выравнивания близких ориентаций у соседних кристаллов. В итоге количество зерен в области уменьшается, а их размер возрастает. Более подробный анализ полученных эффектов будет проведен для двухмерных задач (результаты решения которых готовятся к публикации).

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

Выводы

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

Работа выполнена при финансовой поддержке из средств гранта "Урал-2004" № 04-02-96030 Российского фонда фундаментальных исследований.

Библиографический список

1. Влияние переохлаждение на структуру слитка никеля /Д.Е. Овсиенко, В.П. Костюченко, В.В. Маслов, Г.А. Алфинцев // Кинетика и механизм кристаллизации. - Минск: Наука и техника, 1973. - С. 75-81.

2. 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 // Phys. D. - 1993. - Vol. 69. - P. 189-200.

3. Kobayashi R., Warren J.A., Carter W.C. Vector-valued phase field for crystallization and grain boundary formation // Phys. D. - 1998. - Vol. 119. - P. 415-423.

4. Смирнов Ю.М. Физика кристаллизации. - Калинин: КГУ, 1986. - 81 с.

5. Ландау Л.Д., Лифшиц М.Ф. Механика сплошных сред. - М.: Наука, 1953. - 788 с.

6. Бояршинов М.Г. Численные методы. Часть 2. - Пермь: ПГТУ, 1999. - 200 с.

Получено 15.07.2004.