Решение задачи Стефана сведением ее к задаче теплопроводности с движущимся источником тепла Текст научной статьи по специальности «Математика»

УДК 532.546

РЕШЕНИЕ ЗАДАЧИ СТЕФАНА СВЕДЕНИЕМ

ЕЕ К ЗАДАЧЕ ТЕПЛОПРОВОДНОСТИ С ДВИЖУЩИМСЯ ИСТОЧНИКОМ ТЕПЛА

А. Р. Павлов, Е, А. Слепцова

Существующие методы численного решения задачи Стефана можно разбить на две группы: методы с явным выделением границы фазового перехода (см., например, [1-3]) и методы, основанные на замене эквивалентной задачей теплопроводности с источником (стоком) тепла, связанным с границей фазового перехода. Характерная особенность второй группы методов состоит в том, что в них решение задачи ищется как решение однородной разностной схемы. Впервые такие схемы были предложены одновременно в [4,5]. Пространственное распределение источника тепла в уравнении теплопроводности определяется дельта-функцией Дирака. При построении разностных схем она аппроксимируется дельтаобразной функцией, что означает «размазывание» энтальпии фазового перехода на некоторую окрестность границы раздела фаз. В работе [6] вместо дельта-функции аппроксимируется ее первообразная — единичная функция Хевисайда — и тем самым энтальпия фазового перехода «размазывается» на некотором интервале температур, содержащем температуру фазового перехода. В [7] задача Стефана рассматривается как задача теплопроводности с движущимся источником тепла, связанным с поверхностью фазового перехода.

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

© 2008 Павлов А. Р., Слепцова Е. А.

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

1. Постановка задачи. Пусть имеются две фазы 1 и 2 и распределение температуры в них описывается уравнениями

дТ

c!—= div(A1gradT) + /1, Т<Т», (1)

дТ

c2 — = div(A2 gradT) + /, T > Т*, (2)

dt

c c A A

//

через Ф(г, t) = 0 уравнение границы раздела фаз, r = r(xi, x2, • ••,xp). Если точка лежит на поверхности фазового перехода, то r = R и $(R(t), t) = 0. Условие теплового баланса (условие Стефана) на поверхности фазового перехода задается в виде [4,5]

д

((AgradT)i - (AgradT)2,grad$) + L— = 0, Т = Т*, (3)

dt

где L — энтальпия фазового перехода.

На внешних поверхностях областей, занятых соответственно фа-

t>

t

2. Эквивалентная задача теплопроводности. Введем кусочно непрерывную неотрицательную функцию д(Т), удовлетворяющую следующим условиям:

1) #(Т) ф 0 в промежутке [Т*,Т* + А), а вне его тождественно равняется нулю;

2) 2(Т*) = 1;

3) dg/dT < 0 для Т е (Т*, Т + А).

Тогда справедливо следующее утверждение: сформулированная выше задача Стефана эквивалентна задаче теплопроводности, определяемой уравнением

д т д

c—= div(AgradT) + / + L—, (4)

с граничными и начальным условиями исходной задачи.

Для доказательства воспользуемся методикой работы [7]. Рассмотрим некоторый объем V, ограниченный поверхностью А и содержащий внутри себя поверхность раздела фаз £ (рис. 1).

В момент времени £ поверхность раздела фаз £ делит область па две части V, V с соответствующими внешними поверхностями А1, А2. За малое приращение времени поверхность раздела фаз проходит через элементарное приращение объема ¡V и займет положение £'. При этом увеличится объем, занятый фазой 2, и соответственно уменьшится объем фазы 1.

Проинтегрируем уравнение (4) по объему V и воспользуемся формулой Гаусса — Остроградского

Рис. 1.

/

(5)

V

А

V

V

где п — внешняя нормаль к поверхности А.

Рассмотрим второй интеграл правой части (5):

V

V2

V -¿V

У2+ЙУ

Здесь интеграл по области V — ¡V равен нулю по определению функции д(Т), а второй интеграл по V + ¡V представим в виде суммы интегралов по V и ¡V'.

Г 1ддзу= [Ь^ЗУ+Ш / ь(д)1 — (д)2г аV. } Ы } Ы Дг^о у М

V У2 ¿V

Из определения функции д(Т) следует, что (д)1 = 0 и д(Т) = 1 на фазовой поверхности. При ^ 0 отношение dУ^Дí стремится к «П • СЕ, где «П — локальная скорость элемента СЕ фазовой поверхности по нормали к ней в сторону первой фазы. Кроме того, объем ¡V при этом стягивается к поверхности Е, так что областью интегрирования становится Е. Тогда предыдущее выражение примет вид

I Ь^;dУ = I Ь^;dУ — | Ь< СЕ. (6)

V У2 Г

Напишем соотношение (5) для каждой фазы:

/ ^га<1Т,П)1 СА + J Ьдg■dУ + J f1dУ, (7)

V А1+Е V V

У (Аёга<1Т,П)2СА + J Ь^-С^ + J ^СУ.

А2+Е У2

Второй член правой части в (7) равен нулю по определению функ-дТ

нием из (6). Затем, суммируя полученные соотношения, будем иметь

= !(А®:дАТ,П)аА + У^га<1Т,П)1 СЕ

V А Е

+ !(\&аАТ,П)2<^ + У Ь ддд ¿V + у Ь< ¿Е + У fdУ.

Е V £ V

Теперь, вычитая (5) из последнего соотношения, придем к равенству J(А&adТ,n)1^ + У (Аёга<1Т,П)2СЕ + ^ Ьу*п СЕ = 0.

ЕЕ £

(8)

Если локальную нормаль к £, направленную в сторону первой фазы, обозначить через п*, то в первом интеграле п = — п*, а во втором п = п*:

J[((Лёга<1Т)2 — (Лёга<1 Т)ьп*) + Ьу*п] ¿£ = 0. £

В силу произвольности объема V и соответствующей поверхности £ отсюда следует, что

((А^ас1Т)1 — (Лёга<1Т)2,п*) — Ьи*п = 0.

С учетом соотношений

* 1 (¿Я , ж \ _ кга(1 Ф ( ¿Я , ,Д д Ф ^

V* = I-т^г -КгааФ , п* = --г^г, кгааФ + —— = 0

п ^га<1Ф|\, Л'6 У' * ^га<1Ф| V л ) дt

из последнего равенства получаем условие (3). Утверждение доказано.

3. Расчет температурного поля при сварочном нагреве.

Разработанная методика применена при численном моделировании тем-| пературного поля электродуговой сварки встык тонких пластин. Математическая модель процесса теплопереноса сформулирована в виде следующей двухфазной задачи Стефана [8]:

Е £ Ю — > — Т<) + /, Т<Т*. (И

= Е дХк №) — > — Т^ /. Т>Т*. (Ю)

дТ

Л1— = а(Т — Тс), %\= 0, (11)

дж!

дТ

—Л1— = а(Т — Тс), х1=1ь (12)

дж дТ

— = 0, х2 = 0, (13)

дх

—= — Тс), х2=12 (14)

дх

T(xb x2,0) = Tc = const, (15)

дФ

((AgradT)i - + L— = 0, (16)

at

T(P) = T* = const, P(xbx2,t) еФ(xbx2,t) = 0, (17)

где система координат выбрана таким образом, что ее начало находится

xx x

x

задачи введены следующие обозначения: Tc, Т* — температуры среды и фазового перехода, a — коэффициент теплоотдачи, вторые члены правых частей уравнений (9), (10) учитывают поверхностную теплоотдачу, f — объемный источник тепла, вносимого сварочной дугой, который выражается через параметры дуги в следующем виде [9]:

f = qk/n£exp k [(xi — vc t)2 + x|],

где q — эффективная мощность дуги, q = qUI, n — эффективный КПД процесса нагрева изделия; U, I — напряжение и сила тока соответственно; k — коэффициент сосредоточенности теплового потока

vc

Сформулированная задача согласно предыдущему пункту сводится к решению уравнения

дТ * д ( дТ \ 2 a T,MfMTdg пя,

c^ = kL ax;!, Aaxx-k)— t(t — Tat (18)

с граничными и начальным условиями (11)—(15). Коэффициенты c, A имеют разрывы при T = Т*. Уравнение (18) с учетом определения функции g(T) приводится к виду

дТ А д (- дТ \ 2a.,

{c\, T < T*,

c2 — L dT, T* <T < T* + A,

C1+C2 T t — t

2 L dT, T — T *,

c2, T > T* + A,

и имеет конечные скачки в точках Т* и Т* + А. В [Т*, — Д, Т* + Д] проводится сглаживание функции Л, например, линейной функцией. Тогда в (19) коэффициент Л имеет следующую структуру:

Г Л, Т < т* — д,

Л=1 + Л — ^, Т* — д < Т < Т* + д, 1 Лъ Т > Т* + д.

Полученная задача решена локально-одномерным методом. Расчеты по данному алгоритму проведены для случая сварки двух тонких пластин с размерами 10 см х 8 см, изготовленных из низкоуглеродистой стали. Принимались следующие значения параметров:

<}«<}1 КДЖ егми КДЖ С1 = 3634—-, с2 = 5964

м3 • град' м3 • град'

^ КДЖ , Л2 = 123.408- КДЖ

ч • м • град ч • м • град

д = 5400^3^ ус = 96 = ОШм, Т = 1530° С, ч ч

Тс = 20°С, Ь= 596.4 * , & = 3000Дт,

м3 м2

а = 400 ^^ при Т > 500°С и а = 80^1 К для Т < 500° С.

ЧМ^Л ч м^

Расчеты температурного поля выполнены для двух различных видов

функции источника:

д{Т) = 1 — ^Т*-, (20)

д(Т) = ехр{—0 ,69(Т — Т*)/Д}. (21)

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

Т

д Т Т д Т

Т

сглаживания (4,5).

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

Таблица. Распределение температуры в точке B по времени

Т,° C t = 9 сек t = 18 сек t = 27 сек t = 36 сек t = 45 сек t = 54 сек

Т 119.006 832,332 1547,538 1440,138 806,254 575,195

Т2 119,006 832,332 1547,538 1440,138 807,037 575,798

Т3 119,006 832,332 1547,538 1306,591 763,339 659,846

ЛИТЕРАТУРА

1. ЕЪгИсЪ L. W. A numerical method of solving heat flow problem with mouving boundary //J. Assoc. Comput. Machinery. 1958. V. 5, N 2. P. 161-176.

2. Васильев Ф. П., Успенский А. В. Разностный метод решения двухфазной задачи Стефана // Журн. вычислит, математики и мат. физики. 1963. Т. 3, № 5. С. 874-886.

3. Вудак В. М., Васильев Ф. П., Успенский А. В. Разностные методы решения некоторых краевых задач типа Стефана // Численные методы в газовой динамике / Сб. работ ВЦ МГУ. М., 1965. Вып. 4. С. 139-183.

4. Самарский А. А., Моисеенко В. Д. Экономичная схема сквозного счета для многомерной задачи Стефана // Журн. вычислит, математики и мат. физики. 1965. Т. 5, № 5. С. 816-827.

5. Вудак В. М., Соловьева Е. П., Успенский А. В. Разностный метод со сглаживанием коэффициентов для решения задач Стефана // Журн. вычислит, математики и мат. физики. 1965. Т. 5, № 5. С. 828-840.

6. Мажукип В. П., Повещепко Ю. А., Попов С. В., Попов Ю. П. Об однородных алгоритмах численного решения задачи Стефана. М., 1985. 23 с. (Препринт / Ин-т прикл. математики им. М. В. Келдыша АН СССР; № 122).

7. Шамсупдар П., Спэрроу Е. М. Применение метода энтальпии к анализу многомерной задачи теплопроводности при наличии фазового перехода / / Теплопередача. 1975. № 3. С. 14-22.

8. Ларионов В. П., Павлов А. Р., Тихонов А. Г., Слепцов О. И. Применение ЭВМ для численного определения температурного поля при сварке встык тонких пластин // Автомат, сварка. 1979. № 11. С. 19-22.

9. Рыкалпн П. П. Расчет тепловых процессов при сварке. М.: Машгиз, 1951.

г. Якутск

20 декабря 2004 г.