Математическое моделирование нестационарной теплопроводности при кристаллизации частицы сферической формы в расплаве с движущимися

границами

Сергеев С.А. (netserg@mail.biysk.ru)

Бийский технологический институт (филиал) Алтайского государственного технического университета им. И.И. Ползунова

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

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

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

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

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

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

лофизических характеристик каждой из фаз. Конвекционные тепловые по-

Рисунок 1. Схема температурного поля.

Задача нестационарной теплопроводности для зоны прогрева расплава формулируется с использованием уравнения переноса энергии

токи отсутствуют.

А

Т

--►

5(т) r

-д r2 - = -, [R(t) < r < 5(т)] (1)

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

T(r,0) = Т0; 7Щт),т) = Tc; Т(5(т),т) = Тз, (2)

где Т0 - температура окружающей среды; Tc - температура кристаллизации; т - время; a - коэффициент температуропроводности; r - текущий радиус; 8(т) - радиус границы прогрева в окружающей среде; гс(т) - радиус частицы.

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

dT V dr & r=8(т)

где X - коэффициент теплопроводности окружающей среды.

Количество теплоты dQ характеризуется известным уравнением

dQ = сТ0 dm, (4)

где dm=Spd8(т), dm - масса прироста зоны прогрева, с - удельная теплоемкость окружающей среды.

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

f = <Р^. (5)

где p - плотность окружающей среды.

Приравнивая выражения (3) и (5), и учитывая, что X=cpa, получим

d8( т ) * dT

0-7- = -a

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

dR( т) T dR( т ) , dT dT

) dr & r-R(т)

В результате получим систему уравнений (1), (2), (6), (7), которая будет описывать зависимость температуры зоны прогрева от текущего радиуса частицы и радиуса границы прогрева в окружающей среде, где зона прогрева - это часть окружающей среды, в которой происходит изменение температуры от Tc до T0.

Решая уравнение (1) с учетом граничных условий (2), найдем уравнение распределения температуры в расплаве с учетом границы прогрева в окружающей среде

Merf( -) - erf( n))+- eeX

T = T0 +(TC-T0)-f-r,-Vr> (8)

4i(erf( -) - erf( в))+fSidkl - f3t±±

где e- Ж, L=,

2л[ах 2л[ат 2л[ат Интеграл

J e x dx - erf (x)

Vtc 0

называется функцией ошибок Гаусса.

При 8(т) -> оо уравнение (8) преобразуется в зависимость, аналогичную уравнению, которое встречается в работе Любова [3]

Т = Т0 +(Тс -Т0)ГввГ]- (9)

в

Для того чтобы выразить связь величины радиуса сферы R(x) и величины зоны прогрева 8(т), для начала используем тепловой баланс на границе раздела фаз (7). В результате совместного решения (7) и (8) имеем

Л

. (10)

= 2в2 *1 - в ехр(в2 {Merf (Z) - erf (в)) + --Х^

cc+q V V z jj

Если в уравнении (10) 8(x) - оо и пренебречь количеством теплоты, которое поступает в частицу за счет движения границы фазового перехода, то получается выражение [3]

(Т -Т){ = 2в2(1 - л/Пфехр(в2)erfc(e)). (11)

Ч

При решении уравнения теплового баланса на границе прогрева в окружающей среде (6) и уравнения (8) получим следующую формулу

(12)

Т0 в

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

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

Л 0.003Н

0.0005 j

0.00251

0.0015]

0.002Н

0.001 Н

Рисунок 2. Зависимость радиуса сферы R от времени т (1 - зона прогрева безгранична 8(т) - да, 2 - зона прогрева конечной длины 8(т) ).

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

влиянии на рост кристалла динамики поведения границы прогрева в окру-

Рисунок 3. График зависимости величины в сферы от значений переохлаждения ДТ (1 - зона прогрева безгранична 8(т) - да, 2 - зона прогрева конечной длины 8(т) ).

жающей среде.

в

На рис. 4 приведены значения величин Р и При возрастании значения переохлаждения AT увеличиваются величины Р и Для нормального решения задачи необходимо потребовать, чтобы Р < иными словами, величина роста зоны прогрева всегда должна превосходить радиус частицы.

Рисунок 4. График зависимостей величин Р и - сферы от переохлаждения AT для задачи с зоной прогрева конечной длины.

R т-.

ООО? □.□015 □ 001

0.0005

AT = 20 К 15 К

10 К 5 К

Рисунок 5. Зависимость радиуса сферы R от времени т для разных значений переохлаждения AT.

т

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

Обозначим Ф,

и Ф 2 =

T - T T

Перепишем уравнения (10) и (12) в следующем виде

2в2 (1 - в expie2 i(erf(z) - erf (в)) +

* * ехр(-в2)

-2Z2 ( 1 - Z exp(z2) (Z)- erf (в))

в

(13)

(14)

Разделим выражение (13) на (14), получим

ф =в3 ехр(\в2)

Из уравнения (15) выразим величину Z

(15)

ф

2 exp°2

(16)

В уравнении (16) W - это функция, известная как функция Ламберта [5], которая удовлетворяет следующему выражению

W (z )eW (z) = z,

где z - комплексное число.

Подставив выражение (16) в (14), окончательно находим

( ехр(в2 )Г3

27Лв3 ехр(в2 )erf (в)- erf

Данные задачи были реализованы в среде MATLAB и Maple.

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

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

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

4. Проведено сравнение решения задачи, в которой учитывается величина зоны прогрева с решением задачи, в которой зона прогрева считается бесконечной.

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

1. Карташов Э. М. Аналитические методы решения краевых задач нестационарной теплопроводности в областях с движущимися границами. ИФЖ. 2001. т. 74. № 2. С. 171-195.

2. Лыков А. В. Теория теплопроводности. - М.: Высшая школа, 1967. - 599 с.

3. Любов Б. Я. Теория кристаллизации в больших объемах. - М.: Наука, 1975. - 256 с., ил.

4. Самойлович Ю. А., Тимошпольский В. И., Трусова И. А. Анализ кристаллизации переохлажденного расплава методом интегрального баланса. ИФЖ. 2001. т. 74. № 1. С. 139-144.

5. Corless R.M., Gonnet G.H., Hare D.E., Jeffrey D.J., Knuth D.E. On the Lambert W function. Advances Computational Maths. 1996. Vol. 5, p.

329-359.