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

УДК 539.3

Е. Л. Кузнецова, И. А. Селин, В. Ф. Формалев

ЗАДАЧА ТИПА СТЕФАНА С ПРОИЗВОЛЬНЫМ ЧИСЛОМ ПОДВИЖНЫХ ГРАНИЦ ФАЗОВЫХ ПРЕВРАЩЕНИЙ

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

E-mail: lareyna@mail.ru

Ключевые слова: теплоперенос, задача Стефана, подвижные границы,

композиционные материалы, разложение связующих, аналитические решения.

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

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

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

Ранее на основе общего подхода и моделирования теплового состояния композиционных материалов при высокотемпературном на-гружении рассматривалось численное решение подобных задач [1-5].

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

границ положительна; тепловые эффекты и теплофизические характеристики постоянны; различные фазы отличаются друг от друга тепло-физическими характеристиками.

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

д 2Тп 1 д Тп

д х2 an dt д2 Т, 1 д Т,

= 0, 0 < x < жП_1 (t), t > 0; = 0, x* (t) < x < x*_ 1 (t),

д х2 as дt

t > 0, s = 2,...,n - 1; д2 Ti 1 д Ti

дх2 a1 дt

= 0, x* (t) < x < то, t > 0;

(1)

(2) (3)

-А,

д Т,

дх

+ А,-1

д Ts-1

L(i)-0

дх

= 1Ps_1x*_ 1 (t)

l(t)+0

x = x*_ 1 (t), t> 0, s = 2, ...,n; (4)

TsL=®*_1(i)_o = Ts_1L=®*_1(i)+o = Ts_1, x = x*_1(t), t> 0, s = 2,..., n; T (0,t)= TW, x = 0, t> 0; T (то, t) = T0, x = то, t ^ 0; T (x, 0) = To, x £ (0, то); t = 0;

x* (0) = 0, s =1,n - 1, t = 0.

(5)

(6)

(7)

(8) (9)

В задаче (1)-(9) условия (4), (5) — граничные условия Стефана на подвижных границах, разделяющих фазы под номерами 8 — 1 и 8, 8 = 2,..., п. Здесь А,с, р, а,Т, — теплопроводность, теплоемкость, плотность, температуропроводность, температура, теплота фазовых превращений соответственно.

Наружная граница с температурой ТШ1 (условие (6)) может также рассматриваться подвижной под действием конвективно-кондуктивных тепловых потоков. Однако, чтобы не загромождать выкладки, будем рассматривать границу (х = 0) неподвижной.

x=x

x=x

Задача (1)-(9) — нелинейная, в чем можно убедиться подставив уравнения (1)-(3) в краевые условия (4), (5).

Решением задач для уравнений (1)-(3) с граничными условиями первого рода (5)-(7) будут функции [2]

1 _ erf f x ^ erf f x \ _ ^ / xi (t) \ T (x t)_ T* еГ V ) + erV 2Vait/ f 2-^alt)

1 (X '^ _ Tl 1 - erf ( 4EJ + T0 1 -erf ( ^ЕJ '

x^ (t) < x < to , t > 0; (10)

Ts (x , t) _

T- erf (ц) erf (щ) + (TS - ^ erf ( 2^)

_ erf (Щ ) - erf (%JS ) '

xS (t) < x < xS-1 (t), t > 0; s _ 2 , n - 1; (11)

Tn (x ) t)

erf (——J

_ TW1 - (TW1 - ТП-i) —VV^ , 0 < x < xn-i (t), t > 0. (12)

erf ( ^^M) \ 2\Jant у

В формулах (10)-(12) координаты х* (£), 8 = 1,п — 1, подвижных границ не известны. Для их определения имеются 2п — 2 краевых условия стефановского типа (4), (5). Из представления решений (10)-(12) ясно, что зависимости х* (£), 8 = 1,п — 1, должны быть пропорциональны т.е.

х* (£) = • 8 = 1,п — 1, (13)

где коэффициенты пропорциональности xs, s _ 1,n - 1, подлежат определению, причем поскольку xS-1 (t) > xS (t), то для любых моментов времени постоянные xs-1 > xs. В этом случае решения (10)-(12) с учетом (13) примут форму

T (x t x ) _ T*1 erk2v/a^) + T erk2v/a!t) erf (Xl)

Ti (x ^ ,Xi) _ Ti-1-77—ч--Г T0-:-77—\->

1 - erf(xi) 1 - erf(xi)

xl (t) < x < to , t > 0; (14)

Ts (x,t,Xs-i,Xs) =

erf (x,) - r; erf (xs—1 + (Г; - Г/-1) erf^^)

erf (xa) - erf i Xs-i ^ 1

ж* (t) < x < x*—1 (t), t > 0; s = 2, n - 1; (15)

erf )

Г (x t x ) = Г - (Г - Г* ) — \2Vantj

±n An—1) wi w 1 n—1/

г» I I an— 1

erf Xn—1

an

0 <x<xn_ 1 (t), t> 0. (16)

Для определения коэффициентов xs, s = 1,n — 1, подставим решения (14)-(16) в краевые условия (4), получим систему n — 1 трансцендентных уравнений относительно xs, s = 1,n — 1. Не теряя общности всего алгоритма, рассмотрим задачу (1)-(9) для n = 3 с двумя нестационарно подвижными границами x = xj (t) и x = x£ (t), для определения координат которых получаем следующую систему двух трансцендентных уравнений (качественная картина процесса представлена на рис. 1):

exP ( — Xl ~ ) / 24

F (Xi ,X2 ) = --F ff- — X! = 0,

fxi — erf (X2) 1 — ^ (Xl)

2 (17)

exp

fii (x1,x2) = c

2 a2 аз

где

— D

exp (-x2)

E =

Л2С2Р2С1 TT - tt nAipi Q1

F =

а1 ) - erf (X2) a2 .

Ci TT - To

- X2 = 0,

П

QT

(18)

С=

A3C3P3C2 TWl - T* ПА2Р2 Q2

D =

T2T - tt

Q2

Поскольку Т/ (х 1, Х2), Рц (х 1, Х2) — дифференцируемые функции, то для нахождения х 1,Х2 можно применить итерационный процесс Ньютона (при условии, что на каждой итерации матрицы Якоби системы (17), (18) — не вырождены). Однако для применения итерационных процедур необходимо найти начальное приближение вектора

неизвестных ^х10),х2°^ . Для системы двух уравнений это можно

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

( (°) (°)

пересечения используем для задания начального вектора х1 , х2 после чего применяем итерационный процесс уточнения х1 и х2.

Подставив далее найденные значения х1, х2 и теплофизические характеристики в решения (13)-(16), получим нестационарное температурное поле в трех областях с двумя нестационарно подвижными границами, координаты которых определены по формуле (13).

Сложность решения системы (17), (18) двух трансцендентных уравнений заключается в том, что элементы матрицы Якоби

дТ/ (х1,х2) дТ/ (х1,х2)

dXi

dX2

(х1,х2) дТи (х1,х2) дх1 дх2

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

Поэтому компоненты начального вектора ^х10), х2°^ необходимо вычислять с высокой точностью, с отклонением от компонент точного вектора (х1,х2)т на несколько процентов.

г, к 2000

1500

1000

500

0

Рис. 2. Распределение температур в трех областях, ограниченных подвижными границами x*(t), x£ (t), для разных моментов времени t при Q = Ql = Ql = 1000 кДж/кг; T22 = 600 K; T£ = 900 K (a), T£ = 1300 K (б)

Анализ результатов. Для варьируемых значений Q* , Q* , Т* , Т2* и входных данных Ai = Л2 = Л3 = Л = 0,001 кВт/(мК); с1 = с2 = с3 = = c = 1,5 кДж/(кгК); р1 = р2 = р3 = 2000 кг/м3; а1 = а2 = а3 = а =

= — = 1 • 10-6 м2/с; То = 300 К; TW1 = 2000 К; х* (0) = 0; х2* (0) = 0 с р 3

получены решения задачи (1)-(9) в виде распределения температур и изменения координат х* (t), х* (t) подвижных границ (рис. 2-4).

На рис. 2 приведены зависимости от времени температурных полей в трех областях для случая Q = Q* = Q* = 1000 кДж/кг при различных температурах Т* и T2*. Четко видны положения границ x1 (t) и x* (t) по разрыву касательных к графикам функций, причем если разно сть ДТ* = T2* - Т* мала (ДТ* = 300К, см. рис.2,а), то профили температур находятся на более близком расстоянии друг от друга по сравнению со случаем ДТ * = 700 К на рис.2, б. При этом скорость движения границ в существенной степени зависит от температур Т* и Т2* и теплот фазовых превращений Q 1 и Q*.

На рис.3, 4 приведены координаты подвижных границ x* (t) и x2 (t), причем рис. 3 соответствует фиксированным Т* и Т2* и варьируемым значениям Q = Q 1 = Q2, а рис. 4 — фиксированным значениям Q = Q* = Q* и варьируемым значениям Т* и Т2*. Из рисунков видно, что при фиксированных значениях Т * и Т2* и монотонно возрастающих значениях Q = Q1 = Q* (рис. 4) скорости движения границ x1 (t) и x* (t) монотонно убывают (с возрастанием Q* расстояния между кривыми 1-3 на рис. 4 убывают).

Наоборот, при фиксированных значениях Q = Q* = Q* и монотонно возрастающих Т*, Т2* (рис. 4) скорость движения границы x* (t) монотонно убывает, а скорость движения границы x* (t) монотонно возрастает (несмотря на то, что Т* = 600К = const); это является неожиданным результатом.

Во всех случаях скорости x* (t), x* (t) движения границ с течением времени уменьшаются, что является следствием уменьшения разности тепловых потоков на подвижных границах из-за увеличения расстояния от источника тепла на наружной границе w 1.

а б

Рис. 3. Изменение координат нестационарно подвижных границ ж*(£), ж* (£) при Т* = 600 X, Т* = 1200 Ю

Я = Я* = Я* = 600 кДж/кг (7); Q = Я* = Я* = 1000 кДж/кг (2); Я = Я* = Я* = 1400 кДж/кг (3)

Рис. 4. Изменение координат нестационарно подвижных границ ж* (£), ж* (£) при д = д* = д* = 1000 кДж/кг:

Т* = 600 К; Т2* = 900 К (7); Т2* = 1100 К (2); Т2* = 1300 К (3)

Выводы. Поставлена и аналитически решена задача типа Стефана со многими нестационарно подвижными границами. Анализ полученных результатов показал, что температурные профили в точках х* (£) и X (£) имеют изломы касательных в соответствии с условиями Стефана, а скорости движения границ в существенной степени влияют друг на друга.

Работа выполнена при финансовой поддержке Гранта Президента Российской Федерации МК-7784.2009.8., поддержке ФЦП "Научные и научно-педагогические кадры инновационной России" на 2009-2073 годы в рамках мероприятия 7.2.7 "Проведение научных исследований научными группами под руководством докторов наук", госконтракт № П887 от 78 августа 2009 г., и ФЦП "Научные и научно-педагогические кадры инновационной России" на 2009-2073 годы в рамках мероприятия 7.7 "Проведение научных исследований коллективами научно-образовательных центров", госконтракт № 02.740.77.0227 от 77 июля 2009г.

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

1. Формалев В. Ф., ФедотенковГ.В., КузнецоваЕк. Л. Общий подход к моделированию теплового состояния композиционных материалов при высокотемпературном нагружении // Механика композиционных материалов и конструкций. - 2006. - Т. 12. - № 1. - С. 141-156.

2. К а р с л о у Г., Е г е р Д. Теплопроводность твердых тел. - М.: Наука, 1964. -487 с.

3. Формалев В. Ф., Колесник С. А., Миканев С. В. Моделирование теплового состояния композиционных материалов // Теплофизика высоких температур. - 2003. - Т. 41. - С. 935-941.

4. Кузнецова Е к. Л. Математическое моделирование тепломассопереноса в композиционных материалах при высокоинтенсивном нагреве: Дис. ... канд. техн. наук. - МАИ, 2006. - 138 с.

5. Ф о р м а л е в В. Ф., К у з н е ц о в а Е к. Л. Многомерный теплоперенос при наличии фазовых переходов в анизотропных композиционных материалах // Механика композиционных материалов и конструкций. - 2007. - Т. 13. - № 4. -С. 129-141.

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

Екатерина Львовна Кузнецова родилась в 1980 г., окончила

МАИ в 2003 г. Канд. фнз.-мат. наук, доцент кафедры "Прикладная математика и механика" Московского авиационного института (государственного технического университета). Автор более 30 научных работ в области математического моделирования задач тепломассопереноса.

Ye.A. Kuznetsova (b. 1980) graduated from the Moscow Aviation Institute in 2003. Ph. D. (Phys.-Math.), assoc. professor of "Applied Mathematics and Mechanics" department of the Moscow Aviation Institute (state technical university). Author of more than -30 publications in the field of mathematical simulation of heat- I mass-transfer problems.

Владимир Федорович Формалев родился в 1938 г., окончил в 1968 г. МАИ им. С. Орджоникидзе. Д-р физ.-мат. наук, профессор кафедры "Вычислительная математика и программирование" Московского авиационного института (государственного технического университета). Заслуженный деятель науки РФ. Автор более 170 научных работ в области математического моделирования задач тепломассопереноса и вычислительной математики.

V.F. Formalev (b. 1938) graduated from the Moscow Aviation Institute n. a. S. Ordzhonikidze in 1968. D. Sc. (Phys.-Math.), professor of "Computational Mathematics and Programming" department of the Moscow Aviation Institute (state technical university). Honored Science Worker of the Russian Federation. Author of more than 170 publications in the field of mathematical simulation of heat-mass-transfer problems and computational mathematics.

Илья Александрович Селин родился в 1983 г., аспирант кафедры "Вычислительная математика и программирование" Московского авиационного института (государственного технического университета). Автор 5 научных работ в области математического моделирования тепловых волн.

I.A. Selin (b. 1983) — post-graduate of "Computational Mathematics and Programming" department of the Moscow Aviation Institute (state technical university). Author of 5 publications in the field of mathematical simulation of heat waves.