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

УДК 532.78

Г. Н. К у в ы р к и н, А. К. Лепешкин

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

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

В настоящее время в технике все большее применение находят материалы с особыми физическими свойствами, обусловленными их внутренней структурой. К ним относятся аморфные металлы [1, 2], а также металлы с особым строением кристаллической решетки [3]. В зависимости от состава и термообработки эти материалы проявляют различные свойства.

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

Методы закалки из жидкого состояния имеют несколько разновидностей. В данной работе для нас представляют интерес методы распыления расплава и кавитации, применяемые для производства порошков. С помощью этих методов получают не только аморфные, но и кристаллические порошки. Возможность получения аморфной или кристаллической структуры зависит от скорости охлаждения расплава. Для аморфных материалов это, как правило, 104-106 К/с и выше. Для кристаллических порошков — порядка 103-104 К/с.

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

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

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

Пусть фазовый переход происходит на границе раздела фаз Б = = Б (*), где t — время. Эта граница разделяет всю расчетную область П на две подобласти. В дальнейшем индексом "+" будем обозначать величины, соответствующие жидкой фазе, а индексом "—" — твердой фазе.

В связи с тем, что фазовый переход сопровождается выделением/поглощением теплоты и на границе фазового перехода происходит скачкообразное изменение теплофизических свойств материала, на этой границе имеет место разрыв термодинамических параметров процесса. Поэтому на границе фазового перехода должны выполняться условия сопряжения, называемые условиями Стефана [4-8]. Это условие неразрывности температуры

[Т (р, ¿)] = о, р е Б

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

-Л^ (Р,*) = Р$*А Р е Б. (1)

Здесь [•] = (•)+ — (-)- означает скачок функции при переходе через поверхность разрыва со стороны твердой фазы на сторону жидкой фазы; п — единичный вектор нормали к поверхности Б, направленный в сторону ее движения; Л(Т) — теплопроводность материала; р(Т) — плотность материала; О — модуль нормальной скорости движения границы фазового перехода.

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

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

В простейшем случае обобщенное уравнение теплопроводности имеет вид [4, 6-8]

^ = V (АУТ), (2)

где с(Т) — эффективная теплоемкость,

с(Т) = с (Т) - (Т — Т*); (3)

с(Т) — удельная массовая теплоемкость материала; 8(Т—Т*) — дельта-функция.

Уравнение теплопроводности (2) приемлемо для описания медленно протекающих процессов. В случае высокоинтенсивного процесса оно нуждается в уточнении.

Поскольку состояние рассматриваемой термодинамической системы быстро изменяется во времени, то рассматриваемый термодинамический процесс является неравновесным. Поэтому для учета неравновесности процесса введем в рассмотрение термодинамическую температуру Ф (М,Ь), М е П, которая совпадает с абсолютной, если скорость изменения Ф стремится к нулю. Термодинамическая температура служит мерой средней кинетической энергии движения молекул в неравновесном процессе.

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

Изменение введенных параметров состояния описывается соответствующими кинетическими уравнениями. В линейном приближении примем следующий вид кинетических уравнений:

для термодинамической температуры [9]

тт Ф = Ф — Ф, (4)

для параметра порядка [10]

Т8Х. + х = П (Т — Т*), (5)

где тт — время релаксации параметра Ф, обратно пропорциональное частоте собственных колебаний молекул в неравновесном процессе; Ф — равновесное значение этого параметра; т^ — время релаксации параметра х; П (Т — Т*) — функция Хевисайда. Решение уравнения (4) имеет вид г

* Ж [ ( 1 — т \ дФ ,

Ф = Ф — ехр--

У V тт / дт

о

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

X =1 - - ехр(-П (Т* - Т), (6)

где г* — время начала фазового перехода, соответствующее достижению температурой значения T*.

Поскольку термодинамическая температура совпадает с абсолютной, если скорость ее изменения стремится к нулю, то в условиях термодинамического равновесия Ф = T. Тогда

г

г-Л дт

Ф = T -J expl(7)

0

Представим объемную плотность свободной энергии в виде

pA (T, Ф) = р (B (T) + Bi (T, Ф)).

При температуре Te = const естественного состояния термодинамической системы B (Te) = 0, B1 (Te,Te) = 0. Тогда массовая плотность энтропии выражается соотношением

h = _ dA = _ dB_dB1 h dT dT dT. (8)

Поскольку затвердевание жидкости сопровождается выделением определенного количества теплоты, то можно считать, что на S имеется поверхностный тепловой источник с некоторой мощностью qs. Тогда для всей расчетной области П можем записать закон сохранения энергии в виде [6]

dh

pT— = -V- q + isqs. (9)

dt

Здесь слагаемое 8Sqs представляет собой объемную мощность поверхностного источника теплоты, а 8S — поверхностная 5-функция, определяемая так, что для любой функции f (M) справедливо соотношение

J 5s f (M) dfi = y f (M) dS. П s

Из (9) с учетом (8) получим

id2B д2Bi \ dT ^ д2Bi дф „

-pT(dT2 + -ST21 Ж - ~dTd$ m = -V-q +5sqs.

Положим с = —Т (й2Б/йТ2 + д2Бх/дТ2) = -Тд2Б1 /(дТдФ). Тогда, используя закон теплопроводности Фурье я = —ХЧТ, а также представление термодинамической температуры (7), получим уравнение теплопроводности в следующем виде:

dT pc dt tt

exp

I \ {Г)Л7~1

-— —dT = V (AVT) + Ssqs. (10)

TT J OT

Вместо условия Стефана (1) с учетом параметра х имеем

A

dT (P,t) dn

= pQ*Dx, P e S.

(11)

Разность нормальных составляющих q- и д+ векторов плотности теплового потока (соответственно отводимого от поверхности фазового перехода Б и подводимого к этой поверхности) равна мощности поверхностного источника теплоты qs:

q

q+ = qs.

Поскольку направления векторов плотности теплового потока я- и противоположны направлению нормали п к Б, то с учетом соотношения (11) из последнего равенства получаем

qs = рЯ*БХ. (12)

Покажем [6], что

T

SsD = S (T - T*) —.

t

(13)

Вблизи границы фронта кристаллизации введем локальную ортогональную систему координат (х',у',г'), метрические коэффициенты которой равны единице. В этих новых координатах поверхностная 5-функция 5S есть 5S = 5 (х' — х0), где уравнение х' = х0 определяет границу Б. Аналогично для скорости движения свободной границы имеем Б = йх'/йЬ. Условие Т (М,Ь) = Т*, М е Б(Ь) соответствует тому, что в новых координатах Т = Т (х', у', г', Ь), Т (х'0, у', г', Ь) = Т*. С учетом этого получаем выражение (13).

Тогда уравнение теплопроводности (10) с учетом (6), (12) и (13) примет вид Г

дТ рс Г ( Ь — т \ дТ

pcs ~дЬ + TT,

exp

TT

—dT = V (AVT), (14)

где cS (T) — эффективная теплоемкость,

Cs (T, t) = c (T) -

- QU 1 - 1 - exp

t -1* ts

n (T* - T)\ S (T - T*). (15)

t

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

Рассмотрим каплю расплава алюминия. Объем капли У0 = = 32п/3 мм3. Материал капли изотропен. Теплофизические характеристики материала (плотность р, удельная массовая теплоемкость с и теплопроводность Л) зависят от температуры Т и приведены в табл. 1. Температура плавления алюминия Т* = 933 К, удельная теплота фазового перехода Q* = 380000 Дж/кг. Температура в начальный момент времени Т0 = 1100 К. Температура подложки Тп = 300 К. Ввиду высокой интенсивности теплообмена с подложкой, теплообменом с окружающей средой можно пренебречь.

Таблица 1

Теплофизические свойства алюминия

T, K р, кг/м3 Л, Вт/(м-К) с, Дж/(кг-К)

200 2750 201 850

300 2703 207 900

400 2670 213 948

500 2645 222 993

600 2620 233 1041

700 2600 251 1088

800 2570 271 1135

900 2540 282 1183

1000 2400 150 1085

1100 2350 120 1050

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

Рис. 1. Равновесная форма капли

Рис. 2. Пространственная расчетная область

свободной энергии [11]. Он заключается в том, что при постоянных температуре и объеме физически наблюдаемой конфигурацией будет та, для которой свободная энергия минимальна. Поскольку данная задача непосредственно не касается цели исследования данной работы, опустим ее рассмотрение. Приведем лишь результат ее решения. На рис. 1 представлена половина осевого сечения капли.

Итак, мы имеем задачу расчета температурного поля в области ^ (рис.2). Здесь $1 — поверхность капли, контактирующая с окружающей средой, $2 — поверхность капли, контактирующая с подложкой. Поскольку область ^ осесимметричная, то можем ввести цилиндрическую систему координат (г, г).

Поставим начальное и граничные условия для уравнений (2) и (14) в области Начальное условие (при I = 0) зададим в виде

Т (М, 0)= То(М), М е П,

где Т0(М) — известная функция.

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

Р е $1 : ^=0;

д п

Р е $2 : Т (Р,Ь) = Тп,

где п — нормаль к поверхности $1.

Температура в любом осевом сечении расчетной области полностью определяется двумя координатами г и г. Поэтому вместо про-

/

-5»

г

странственной расчетной области ^ можно рассматривать симметричную половину Е осевого сечения области ^ (рис. 3).

Начальное и краевые условия для плоской области Е примем в следующем виде:

где 11 состоит из двух частей: части границы области Е, соответствующей части границы Б1 области и части, соответствующей оси симметрии области а 12 есть часть границы области Е, соответствующая части границы Б2 области

Некоторые аспекты численного решения. Для численного анализа построенных математических моделей воспользуемся методом конечных элементов [12-14]. Будем аппроксимировать область Е симплексными треугольными конечными элементами.

Зависимость теплофизических свойств материала от температуры аппроксимируем кусочно-линейными функциями.

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

При использовании традиционной формулировки метода конечных элементов в начальные моменты времени нестационарного процесса теплопроводности могут возникнуть осцилляции решения в точках конечно-элементной сетки. Эти осцилляции могут быть объяснены проявлением неконсервативности схемы метода конечных элементов для одного элемента при консервативности всей системы в целом [16]. С целью избавления от осцилляций решения применяют процедуру диагонализации матрицы теплоемкостей [10].

T (M, 0) = To, M e F,

P e ¡2 : T (P,t)= Tn,

Шаг интегрирования по времени выбирается автоматически [17,18].

В выражения (3) и (15), определяющие соответственно c(T) и cS (T, t), входит S-функция. При численном решении задачи ее необходимо сглаживать [6]. Для этого применяются различные методы сглаживания, однако выбор метода не оказывает существенного влияния на результаты расчетов [5]. Более существенное влияние оказывает ширина интервала сглаживания.

Заменим S (u) функцией S (u, А), отличной от нуля только внутри некоторого интервала сглаживания [—А, А]. Интервал сглаживания и вид функции S (u, А) определим исходя из условия

сю

J S (u, А) du = 1. Аппроксимацию S (u, А) примем в виде

i оА, |u| < А; S (u, А) = I оА

[о, |u| > А. Тогда выражения для с (T) и cS (T, t) запишем в виде

с (T) = c (T) - Q*S (T - T*, А) ;

CS (T,t) = c (T) -

- QN 1 - ' / * - **

n (T * - T U S (T - T *, A).

1 — ехр

тя

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

Для определения влияния параметров состояния вещества на развитие процесса затвердевания при его численном моделировании проведена серия расчетов. Сначала проведен расчет для простейшей модели затвердевания, в результате чего определено время затвердевания капли £з в поставленных условиях: £з = 0,000918083 с. Затем была проведена серия расчетов для модели с учетом параметров состояния при различных значениях отношения времен релаксации тт и т? (соответственно термодинамической температуры Ф и параметра порядка х) к £з.

На рис. 4 приведено распределение температуры по оси капли в близкие моменты времени при различных значениях времен релаксации параметров состояния, а в табл. 2 — время затвердевания капли

Рис. 4. Распределение температуры по оси капли:

1 — модель без учета параметров состояния, Ь = 4,49642 • 10-4 с; 2 — модель с учетом параметров состояния, Ь = 4,50217 • 10-4 с, тт/Ьз = 10-6, т^/Ьз = 10-6; 3 — модель с учетом параметров состояния, Ь = 4,50115 • 10-4 с, тт/Ьз = 10-6, т^/Ьз = 1 4 — модель с учетом параметров состояния, Ь = 4,49924 • 10-4 с, тт/Ьз = 10 т^/Ьз = 10-6; 5 — модель с учетом параметров состояния, Ь = 4, 49014 • 10 тт/Ьз = 10-3, тБ/Ьз = 1; 6 — модель с учетом параметров состояния, Ь = 4, 50125 • 10-4 с, тт/Ьз = 1, т^/Ьз = 10-6; 7 — модель с учетом параметров состояния Ь = 4, 50223 • 10-4 с, тт/Ьз = 1, т5/Ьз = 1

-4

для модели с учетом параметров состояния при различных значениях времен релаксации этих параметров.

Таблица 2

Время затвердевания капли для уточненной модели

TT Дз TS/*з = 10-6 TS Лз = 10-3 TS Лз = 1

10-6 0,000931083 0,000969864 0,001015082

10-3 0,001078104 0,001118245 0,001173212

1 0,001224507 0,001278568 0,00133012

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

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

СПИСОК ЛИТЕРАТУРЫ

1. Аморфные металлические сплавы / Под ред. Ф.Е.Люборского: Пер. с англ. - М.: Металлургия, 1987. - 584 с.

2. Судзуки К., Фудзимори Х., Хасимото К. Аморфные металлы: Пер. с япон. - М.: Металлургия, 1987. - 328 с.

3.Абрамов А. В., Дерягин Н. Г., Третьяков Д. Н. Применение сверхбыстрого (102... 103 °С/с) охлаждения раствора-расплава в жидкофазной эпитаксии полупроводников // Физика и техника полупроводников. - 1999. -Т. 33. - Вып. 9. - С. 1130-1133.

4. А в д о н и н Н. А. Математическое описание процессов кристаллизации. -Рига: Зинатне, 1980. - 180 с.

5. Вабищевич П. Н. Численные методы решения задач со свободной границей. - М.: Изд-во Моск. ун-та, 1987. - 164 с.

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

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

8. Р у б и н ш т е й н Л. И. Проблема Стефана. - Рига: Звайгзне, 1967. - 457 с.

9. З а р у б и н В. С., К у в ы р к и н Г. Н. Математические модели термомеханики. — М.: Физматлит, 2002. - 168 с.

10. Кувыркин Г. Н. Термомеханика деформируемого твердого тела при высокоинтенсивном нагружении. - М.: Изд-во МГТУ им. Н.Э. Баумана, 1993. -142 с.

11. Шехтер Р. Вариационный метод в инженерных расчетах. - М.: Мир, 1971. -291 с.

12. В л а с о в а Е. А., Зарубин В. С., Кувыркин Г. Н. Приближенные методы математической физики. - М.: Изд-во МГТУ им. Н.Э. Баумана, 2001. -700 с.

13. З а р у б и н В. С., С е л и в а н о в В. В. Вариационные и численные методы механики сплошной среды. - М.: Изд-во МГТУ им. Н.Э. Баумана, 1993. - 360 с.

14. Зенкевич О. Метод конечных элементов в технике / Пер. с англ. - М.: Мир, 1975.- 541 с.

15. Кувыркин Г. Н., Шаров С. Н. Решение двумерных задач теплопроводности с подвижной внешней границей методом конечных элементов // Изв. вузов. Машиностроение. - 1985. - Вып. 7. - С. 52-56.

16. Трощиев В. Е., Шагалиев Р. М. Консервативные узловые схемы методов конечных разностей и конечных элементов для двумерного уравнения теплопроводности // Численные методы механики сплошной среды. - 1984. -Т. 15.-Вып.4.-С. 131-157.

17. К у в ы р к и н Г. Н., Панин С. Д., Ц и ц и н А. Г. Особенности численного решения задач нестационарной теплопроводности при высокоинтенсивном нагреве // Изв. АН СССР. Энергетика и транспорт. - 1988. - Вып. 5. - С. 162-165.

18. Особенности численного решения задач нестационарной теплопроводности для составных тел переменного объема / В.С. Зарубин, Г.Н. Кувыркин, А.Г. Цицин, С.М. Шаров // Мат. моделир. в науке и техн. Всес. школа-семинар. Тез. докл. - Пермь, 1986. - С. 144-145.

Статья поступила в редакцию 26.03.2007

Георгий Николаевич Кувыркин родился в 1946 г, окончил МВТУ им. Н.Э. Баумана в 1970 г. Д-р техн наук, профессор кафедры "Прикладная математика" МГТУ им. Н.Э Баумана. Автор более 150 научных работ в области термомеханики и термопрочности элементов конструкций.

G. N. Kuvyrkin (b. 1946) graduated from the Bauman Moscow Higher Technical School in 1970. D. Sc. (Eng), professor of "Applied Mathematics" department of the Bauman Moscow State Technical University. Author of 150 publications in the field of thermal mechanics and thermal structural mechanics of construction elements.

Антон Костантинович Лепешкин родился в 1982 г., окончил МГТУ им. Н.Э. Баумана в 2006 г. Аспирант кафедры "Прикладная математика" МГТУ им. Н.Э. Баумана. Специализируется в области моделирования процессов с фазовыми переходами.

A.K. Lepyoshkin (b. 1982) graduated from the Bauman Moscow State Technical University in 2006. Post-graduate of "Applied Mathematics" department of the Bauman Moscow State Technical University. Specializes in the field of simulation of processes with phase transitions.

"ВЕСТНИК МОСКОВСКОГО ГОСУДАРСТВЕННОГО ТЕХНИЧЕСКОГО УНИВЕРСИТЕТА имени Н.Э. БАУМАНА"

Журнал "Вестник МГТУ имени Н.Э. Баумана" в соответствии с постановлением Высшей аттестационной комиссии Федерального агентства по образованию Российской Федерации включен в перечень периодических и научно-технических изданий, в которых рекомендуется публикация основных результатов диссертаций на соискание ученой степени доктора наук.

Журнал издается в трех сериях: "Приборостроение", "Машиностроение", "Естественные науки" с периодичностью 12 номеров в год. Подписку на журнал "Вестник МГТУ имени Н.Э. Баумана" можно оформить через ОАО "Агентство "Роспечать".

Подписка по каталогу "Газеты, журналы" ОАО "Агентство "Роспечать"

Индекс Наименование серии Объем выпуска Подписная цена (руб.)

Полугодие 3 мес. 6 мес.

72781 "Машиностроение" 2 250 500

72783 "Приборостроение" 2 250 500

79982 "Естественные науки" 2 250 500

Адрес редакции журнала: 105005, Москва, 2-я Бауманская ул., д.5. Тел.: (495) 263-62-60; 263-60-45. Факс: (495) 261-45-97. E-mail: [email protected]