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

УДК 536.1; 539.2

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

Н.Д. Няшина, П.В. Трусов (Пермь)

Abstract

The phase-field model described by system of nonlinear equations is investigated. The steady-state solutions of these equations and the sets of initial conditions leading to these solutions are determined through consideration of the homogeneous problem. The analysis of the linearized stability problem makes it possible to find a spectrum of wave numbers of unstable perturbations, their evolution being investigated numerically for tested problem.

Введение

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

Классическая модель Стефана не учитывает влияние таких физических факторов, как поверхностное натяжение расплава и концентрационное

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

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

Введенная переменная определяется из рассмотрения задачи микроуровня с применением статистической термодинамики [1], допускает существование жидких участков, переохлажденных или перенасыщенных, и учет влияния кривизны фронта.

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

Идеи МФП впервые описываются для затвердевания чистых металлов в работах А.Р. Уманцева, А.Л. Ройтбурга [2], G. Caginalp [4] и для близкой по сути антифазной

граничной миграции - в работах S.M. Allen, J.W. Cahn [5]. O. Penrose, P. Fife [1] сформулировали условия, по которым уравнения для фазового и температурного полей выводятся из термодинамически обоснованного для фазового перехода I-го рода функционала энтропии.

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

Теоретические основы метода

Рассмотрим закрытую систему объема V (мезоуровень), в котором чистый металл претерпевает фазовый переход I-го рода (L-S). Параметр порядка ф(г, t), называемый фазовым полем, определяет фазу материала в точке r ; ф=0 в твердой фазе, ф=1 - в жидкой. Узкая область (с шириной, стремящейся к нулю), в которой 0<ф<1, соответствует межфазной поверхности раздела. Для простоты предположим однородность плотности среды по всей системе, отсутствие конвекции в жидкости; полагаем, что система находится в механическом равновесии.

Внутренняя энергия Е любого подобъема Q^V представляется соотношением

Е = | edQ, (1)

Q

где е - удельная (на единицу объема) плотность внутренней энергии.

Согласно I закону термодинамики изменение Е со временем подчиняется соотношению:

Е + | q • n d (dQ) = 0,

5Q

(2)

где дО - поверхность О с внешней нормалью п, q - поток внутренней энергии. Используя (1) и теорему Остроградского-Гаусса, соотношение (2) преобразуем (с учетом произвольности выбора объема О и гладкости поля q) к виду

ё + У^ = 0. (3)

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

S =

dQ,

(4)

где s(e,ф) - плотность энтропии (расширение классической термодинамики, где ^ = s(e,p), но плотность здесь полагается постоянной); е=еопв1 (для случая изотропии поверхностного натяжения). Такой вид функционала предлагается по аналогии с теорией Ландау-Гинзбурга для фазового перехода 11-го рода [7]. Второе слагаемое в (4) отражает зависимость термодинамического потенциала не только от значений термодинамических переменных, но и от их градиентов (в частности, Уф ).

Дифференцируя (4) по времени и подставляя выражение (3) для ё, получим

*=1

e +

д s .дф

+ є2 V2ф

ф

dQ -

дО

дs

q+ є2 фV2ф

nd (дО).

(5)

Q

Ф

e

J

e

II закон термодинамики утверждает, что в любом 0,<^¥ производство энтропии положительно. Это производство может быть вычислено выделением из £ потока энтропии через границу дО. Тогда производство энтропии имеет вид

і

+

( д ^ дф.

Е2 V2ф

ф

йО. > 0,

(6)

I

а

е

где Т - абсолютная температура ( у =

д я д е,

для потенциала я = я(е,р,ф )).

р,ф

Таким образом можно гарантировать локально положительность производства энтропии, выбирая

q = Мт V

5 я

—(1

— I = Мт VI —

тср =

5 я ,5фУ

+ є2 V2 ф,

где МТ и т - положительны. Тогда управляющие уравнения примут вид:

для плотности внутренней энергии е = -V •

1

для фазового поля тср = - у

—(1

Мт V

т.

д е

.дф

+ є2 V2 ф .

(7)

(8) (9)

Для определения производной в (9) используется плотность свободной энергии

( д е ^ ( д е ^

Гельмгольца: df = -ъйТ + 1---I ёф . Тогда, учитывая, что ёе = Тй8 + ——I ёф, при

чдф^ „ ^дФ^ _

постоянном я получим

ґдД = ( д е I д " I" I = I д е I = —

чдф) т ^ 5 " Vд т _ т _ ) ф Vд Г — - я _ т _ ) т2 ф

Последнее можно проинтегрировать при ф =С0ПБІ

*ҐТ \ т( т е(С,ф)

1(т,ф) = т - '

V тм

с

йс, + С(ф)

(10)

(11)

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

е = е<:(т) + Р(ф)Ь(Г) = еь (Г) + [(ф) - 1]і(Г) = Єь (Ги) + с(т- ти) + [/»(ф) - і]і(Г), (12)

где е$(т) и еЬ(т) - классические плотности внутренней энергии твердой и жидкой фаз соответственно; Ь(т) = е5(т)-еЬ(т) - их разность. Здесь р(0)=0, р(1)=1; Ь0(т)=Ь(тМ) -скрытая теплота плавления.

С учетом (12) плотность свободной энергии в (11) примет вид

(

I (т, ф) = т

| (С, ф) - - 1| £(т) + С(ф) , где д(т) = I .

Л

V ти

с

с

Заметим, что Q(TM) = 0. Из (10), (13) получим д f ( д е ^ ёр dG(ф)

-ф] = — | =-Т/Q(T) + Т = -ТQ(T)р'(ф) + ТО’(ф).

. дфуТ Чдф|^ ёф ёф

Таким образом, управляющие уравнения (8), (9) примут вид

(13)

(14)

ф

е

и

д т д

{с + [р(ф) -1] ’(Т)} + ЦТ) р '(ф) = Ю2 т, (15)

тф = Q(T)р'Ы - с,(ф) +е2 V2ф , (1б)

где к= Мт / Т 2 =свт1 коэффициент теплопроводности.

Выберем функцию р(ф) в виде

ф

| g(Q ) ^

Р(ф) = 7-------= ф3(10 - 15ф + 6ф2), где ^(ф) = ф2(1 -ф)2. (17)

| 8(? ) ^

о

Эта функция удовлетворяет условиям нормировки: р(0)=0, р(1)=1, а также условиям минимума потенциала свободной энергии независимо от Q(T), так как р '(ф) = р " (ф) = 0 при ф=0 и ф=1.

Потенциал / должен иметь локальные минимумы по ф при ф=0 и ф=1 для любых значений температур Т; кроме того / непрерывен по Т в точке плавления. Эти ограничения на /и выбранный вид функции р(ф) (17) требуют, чтобы С(ф) также имел минимумы при ф=0 и ф=1; поэтому С(ф) выбирается симметричным относительно ф=1/2, - образным» потенциалом в форме

1

С(ф) -— 8(ф), а > 0 , а=сотг. (18)

4а

Выбор функции в виде полинома четвертой степени обусловлен простотой аппроксимации функции, имеющей два минимума.

В [6] показана связь параметров а и е с материальными параметрами:

+” (

5 Г Г 2 Т d Ф1 5 (19)

а= Jе T«jx) • (19)

— то

где 5 - толщина фронта (является подгоночным параметром модели); а - поверхностная энергия на единицу площади (для одномерной задачи), которая является отклонением свободной энергии от _ДТМ,ф).

Получим уравнения в безразмерном виде. Для этого введем масштаб длины: w -

размер расчетной области, (w2/к) = (w2/к)с - характерное время распространения температуры по области, где к/с = к = const - температуропроводность в жидкости. Тогда

X ~ t Т — Тм

X - ; t - / - , \ ; Ы - _ __ „

™ ( /к) ТМ - Т0

где Т0 - некоторая температура на внешних границах, задаваемая в условиях Дирихле. В безразмерных величинах управляющие уравнения примут вид

[ г лЬ'(-) | ды Ь(ы) дф ~2

I1+[р(ф) - 1]Нг}аг+4г р'^зГ-~2«- (20)

е2 дф ~2 ~2 ~ ~ 1

—=е V ф+ еа8<~(— р,(ф)-48,(ф^

т д t 4

где Ь(-) = Ььо, (~(‘-) = Ь—:--------—------------------—, (21)

0 {1 + [(Тм - Т0)/Тм]}

>/ЧЬ0] о с (тм Т0) ~ 5 6^/28о

а = ^----, Ь ----------------, е= — , т ---------—. (22)

12сотм Ь0 ™ тк тм

Для простоты выберем теплоемкость на единицу объема твердой фазы, равной постоянной величине теплоемкости жидкости; Ь(Т)=Ь0 , тогда Ь(-) = 1. Также

предположим, что \ТМ - Т0| << ТМ (это соответствует малому переохлаждению

расплава), тогда согласно (21) Q (-) - -.

Таким образом, уравнения (20) примут вид

~(-) = 0 , р'(ф) - 30ф2(1 - 2ф + ф2) - 308(ф) ,

д- + 308(ф) дф ^ е2 дф е2 ^ ч~ о 1 „ч

+ —Ь—"дт =V -, т~ё~ =е V ф+ 8( ф) еа Ь--4 8 (ф) • (23)

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

Рассмотрим эту модель более подробно.

Пространственно однородная задача

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

Итак, постановка пространственно однородной задачи имеет вид

ё- 308( ф ) ёф е2 ёф 1

—~ +----- —= 0, ——~ = 308(ф) еаЬ--48,(ф);

dt Ь dt т dt 4

ф(0)-Фо, -(0) = -0; 8(ф) = ф2(1 -ф)2 . (24)

Решение системы получено методом Рунге-Кутта. На рис.1 приведены зависимости ф(^, полученные при различных начальных условиях. На рис.2 изображены соответствующие фазовые диаграммы.

На рис.2 изображены пунктирные кривые:

Ь - интегральные кривые при различных начальных условиях (241):

-(ф) - Ь[[) - J(ф0)] + -0 , где У(ф) - 10ф3 - 15ф4 + 6ф5;

Г - кривая, соответствующая стационарному решению (242):

дф = 0. -( ) = 8,(ф) = 1 ф(1 -ф)(1 -2ф)

дt ’ ф 120еаЬ8(ф) 60 е аЬ ф 2(1 -ф)2 .

Решения системы изображены сплошными линиями, частично совпадающими с кривыми Ь, начальные условия отмечены точкой на соответствующей кривой.

Рис. 1. Зависимости ф(7) для различных начальных условий

и

1

0.5

□

-0.5

-1

V ч \ >■ ч ф

N ,

0.2 0.4 0.6 0.8

Рис. 2. Фазовые диаграммы системы (24) для различных начальных условий

0.6

о.

-0.6

а б

Рис.3. Области притяжения стационарных решений для различных начальных условий: а - ф0 е[0,0; 1,0],-{) е[-1,0; 1,0]; б - ф0 е[0,3; 0,7],-{) е[-0,35; 0,35];

□ - ф = 0.0; □ - ф = 1.0; ■ - ф = 0,5

б

а

в

Анализ устойчивости линеаризованной задачи

Вернемся к постановке задачи (24), которая для одномерного случая имеет вид:

ди . 30£(ф) дф д2и е2 дф ~2 д2ф . 30 ( )~ о 1 ,( ) (25)

—~ +---- ----= —-г,------------г~ = е —-г + 30е(ф) еаои--е'(ф), (25)

д~ 0 д~ д£2 т дt д£2 4 4 7

где £ - расстояние от фронта. Дополним ее начальными и граничными условиями:

и(0,£) = ф(0,£) = ф0(£); . .

- = (26) ф^ 0, £ ^ -да ; ф^ 1, £ ^ +да ; и ^ и , £ ^ -да ; и ^ и , £ ^ +да .

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

интервале [0, 1] : (30^(ф)) = ^—0{е(ф)^ф = 1. Кроме того, разложим в ряд вблизи

точки ф0 функцию е(ф), сохранив только линейный член:

е,(ф) = е'(ф0) + е"(ф0)(ф - ф0) + °(ф - ф0) .

Рассмотрим устойчивость линеаризованной системы по отношению к малым возмущения стационарных решений: ф(X, £) = ф(Х, £) -ф0(£), и(Х, £) = и(Х, £) - и0(£) .

Функции начальных условий удовлетворяют линеаризованной системе уравнений (25):

д2и0 д2ф ~ 1

-дгг = 0, е ТГТ + еа0и-4£,(ф0) = 0; (27)

д£ д£ 4

вычитая (27) из (25), получим уравнения и начальные условия для пульсаций:

ди дер д2и е2 дер е2 д2ф ~ $ 1 ,

—~ + ^т~ = —-г,-----------=е —г +еаои- — е (ф0)(р ,

д~ д~ д£ т д~ д£2 4 (28)

и(0) = 0, ф(0) = 0,

граничные условия: ф ^ 0, й ^ 0, £ ^-да ; ф ^ 0, и ^ 0, £ ^+да .

Пусть ф(X, £) = Ф ег“£+^ , и(х, £) = иег“£+ц/ , тогда, подставляя эти выражения в (28), получим систему линейных алгебраических уравнений относительно Ф и и :

(д + ш 2)и + 0 ф = 0,

- є а$и +

^~2 + ~2 2 8 " ( ф с) Л

д+ є ш -

Ф = с.

V т 4 у

Условия существования нетривиального решения последней приводят к квадратному уравнению относительно д:

Є2 2 , Г~2 2 1 +т , ~ 8''(фо)) (~2 4 2 8" (ф С)

— д +1 є ш --------+ єа----------Ід + |є ш -ш ------------1= 0

т V т 4 У V 4 У

или а д2 + Ь д + с = 0, (29)

решая которое получим зависимость декремента затухания возмущений д от волнового числа ш. Неустойчивыми будут возмущения с теми волновыми числами, при которых действительная часть корней уравнения (29) больше нуля: ^е(д)>0.

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

Пусть Б = Ь2- 4 ас>0, тогда корни уравнения д1 и д2 действительны и ^е(д)= д.

і ^ ^ 2 Е" (Ф с)

1. Д1 <с< д2 , если с< С , то есть ш <

4

2

2. д1, д2<С, если

В > С Ъ_ 2а с > С;

< С; о

В > С:

2 2 ^ ~2 2 ~ Е"(ФС) Л 2 Е"(Ф С)

є ш +—є ш + єа- -----------> С;о ш >

4 є

2

4 є2 ш 2 - Е"(Ф с ) >С;

3. д1 , д2 >С , если

В > С Ъ_ 2а с > С;

> С; о

В > С;

~2 2 1 ~2 2 ~ Е (фС) ^ 0

є ш +—є ш +єа-----------------<С;о шє0.

4 є2 ш 2 - Е"(Ф с ) >С;

Дискриминант В уравнения (29) является биквадратным полиномом относительно ш с положительным коэффициентом при ш4.

После того, как установлены общие зависимости знаков корней уравнения (29), рассмотрим частный случай на примере чистого железа, где параметры а, 8, т, е определены согласно (22). Необходимые для их вычисления физические константы приведены в таблице. Для таких констант дискриминант уравнения (29) всегда положителен и уравнение имеет только действительные корни д,1,2 = д (ш2; фо). На рис.4 приведены зависимости д=д(ш) при различных значениях ф0 .

Таблица

Физические параметры чистого железа

Ь , «ж кг Сь , Дж кг 0 К Дж а , м 2 Тм , 0 К к , м! с Т

272 1С3 922 2С,4 18С8 3,7 1С-6 1 1

а б в

Рис. 4. Зависимости декремента затухания д от волнового числа возмущения ш: (а) - фо = 0,0; фо = 1,0; (б) - фо = 0,2; фо = 0,8; (в) - фо = 0,5

Итак, система устойчива к возмущениям с очень маленьким (порядка 1С-5 - 1С-4) и достаточно большим (порядка 15 - 75) волновым числом. При ФС = С,5 система устойчива к любым возмущениям.

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

определить, рассмотрим тестовую одномерную задачу (24). В качестве начальных условий возьмем функцию

ФС(£)=

т

2уІ2а

+1

(3С)

которая является аналитическим решением (242) при равновесных условиях и с условиями на границах ф ^ 0 , £ ^ -да ; ф ^ 1, £ ^ +да, с наложенными на нее неустойчивыми возмущениями в виде суммы синусов и косинусов. Результаты расчетов приведены на рис.5.

1.2

1

0.0

0.6

0.4

0.2

0

: Р

£

Рис. 5. Эволюция неустойчивых возмущений:

(а) - начальное распределение;

(б - г) - развитие возмущения определенной частоты с течением времени

а

в

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

Выводы

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

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

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

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

Работа выполнена при поддержке гранта Министерства образования в области

металлургии

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

1. Penrose O., Fife P. Thermodynamically consistent models of phase-field type for the kinetics of phase transitions // Physica D. - 1990. - Vol. 43. P.44-62.

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

3. Warren J.A., Boettinger W.J. Prediction of dendritic growth and microsegregation patterns in a binary alloy using the phase-field method // Acta metall mater. - 1995. - Vol.43. - № 2.

- P. 689-703.

4. Caginalp G. An analysis of a phase field model of a free boundary // Arch. Rat. Mech. Anal. - 1986. - Vol. 92. - P.205-245.

5. Allen S.M., Cahn J.W. A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening // Acta. Metall. - 1979. - Vol. - 27. - P.1085-1095.

6. Wang S.-L., Sekerka R.F., Wheller A.A., Muray B.T., Coriel S.R., Braun R.J., McFadden G.B. Thermodynamically-consistent phase-field models for solidification // Physica D. -1993. - Vol.69. - P.189-200.

7. Ландау Л.Д., Гинзбург В. Л. К теории сверхпроводимости // ЖЭТФ. - 1950. - Т. 20. -С. 1064-1083.

Получено 15.05.99.