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

Математические заметки СВФУ Январь—март, 2014. Том 21, № 1

УДК 519.63

ЧИСЛЕННОЕ РЕШЕНИЕ ОБРАТНОЙ ЗАДАЧИ СТЕФАНА МЕТОДОМ ВВЕДЕНИЯ РАСПРЕДЕЛЕННОГО ИСТОЧНИКА ТЕПЛОТЫ Е. А. Крылова

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

Ключевые слова: теплота фазового перехода, распределенный источник тепла, задача Стефана, краевая задача для нагруженного уравнения, обратная задача, восстановление правой части уравнения теплопроводности, сеточная задача, разностная схема.

Введение

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

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

В работе рассматривается задача по восстановлению зависимости правой части параболического уравнения от времени при известном распределении по пространству. Такая линейная обратная задача относится к классу некорректных в классическом смысле задач математической физики при специальных предположениях о точках дополнительных измерений — при условии действия источника в точках наблюдения [5]. На базе методики из [7] построен вычислительный алгоритм приближенного решения одномерной по пространству двухфазной обратной задачи Стефана при новом способе учета теплоты фазового перехода — путем введения распределенного в окрестности межфазовой границы источника тепла [8, 9].

© 2014 Крылова Е. А.

1. Постановка задачи

Распределение температуры и положение границы фазового перехода определяется решением системы уравнений:

дТ д ( дТ"

= «г «г )+/(*>*)> 0<Ж<^),

дж у дж

дТ д / дТ\ ,, , ^ , С2^гг = т- Л2— Ф)<х<1.

д£ дж у дж На границе раздела фаз выполнены условия:

Т = Т*, ж = £(£),

х дТ_х дТ_

2 дх 1 дх <М Для уравнений (1.1)—(1.4) зададим граничные условия третьего рода:

дТ

А1 = а(Т-Тс), ж = 0,

^ = 0, х = 1.

дж

1.1) 1.2)

1.3)

1.4)

1.5)

1.6)

Задано также начальное условие

Т(ж, 0) = Тс, 0 < ж < I. (1.7)

В уравнениях (1.1)—(1.7) введены следующие обозначения: с — объемная теплоемкость, Л — коэффициент теплопроводности, Т*,Тс — температуры фазового перехода и окружающей среды соответственно, Ь — скрытая теплота кристаллизации (плавления), а — коэффициент теплоотдачи, ж = £(£) — уравнение границы раздела фаз. Индексы 1 и 2 относятся к фазам, у которых Т>Т* и Т < Т* соответственно.

Прямая задача формулируется в виде (1.1)—(1.7).

В данной работе рассматривается обратная задача, в которой неизвестной помимо Т(ж,£) является и правая часть /(ж, £) уравнений (1.1), (1.2). Будем считать, что /(ж,£) представляется в виде

/ (ж^) = сп(гЖж), (1.8)

где функция ^(ж) задана, а неизвестной является зависимость источника от времени — функция п(£). Эта зависимость восстанавливается по дополнительному наблюдению за Т(ж, £) в некоторой внутренней точке 0 < ж* <1:

Т (ж*,£)= (1.9)

Решение задачи по восстановлению зависимости правой части параболического уравнения от времени при известном распределении по пространству рассматривается при следующих ограничениях:

1) ^(ж*) = 0,

2) ^(ж) — достаточно гладкая функция (^ € С2[0,1]),

3) ^(ж) = 0 на границе расчетной области (для простоты изложения). Особого внимания заслуживает первое предположение, которое связано с

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

3) Ш < 0 для Т е (Т* - Д,Т*).

2. Метод введения распределенного источника теплоты

Для численного решения задачи Стефана (1.1)—(1.7) разработан ряд методов. Наиболее широкое распространение получил метод, описанный в [10,11], где с помощью введения функции энтальпии решение задачи Стефана сводится к решению краевой задачи для уравнения теплопроводности с разрывными коэффициентами. Разностная схема сквозного счета для таких уравнений строится с заменой дельта-функции Дирака дельтаобразной функцией на интервале (Т* — Д,Т* + А). При этом функция энтальпии аппроксимируется различными непрерывными в указанном интервале функциями. В [12] вместо функции энтальпии рассматривается аппроксимация единичной функции Хевисайда. В указанных методах принимается допущение, что фазовый переход при кристаллизации происходит начиная с некоторой температуры, которая выше температуры кристаллизации.

В настоящей работе для численного решения обратной задачи Стефана предлагается модификация учета теплоты фазового перехода, более точно описывающая реальный процесс тепловыделения на границе фазового перехода, путем введения распределенного в окрестности (Т* — А, Т*] (в сторону образующейся фазы) источника тепла [8]. Для этого вводится кусочно непрерывная неотрицательная функция, удовлетворяющая следующим условиям:

1) у(Т) определена во всем диапазоне изменения температур, отлична от нуля в промежутке (Т* — А, Т*), а вне его тождественно равна нулю;

2) у(т*) = 1;

ду дТ

Утверждение. Если г>(Т) удовлетворяет указанным выше условиям, то уравнения (1.1), (1.2), (1.4) заменяются одним уравнением (2.1), заданным во всей области 0 < х < I:

дТ д / дТ\ Эу .

ст=1Гх[х^)+ьт+ПхЛ (2Л)

где с = сх, А = Ах при Т > Т* и с = с2, А = А2 при Т < Т*.

Доказательство данного утверждения по методике [10] приведено в работах [8, 9]: в [8] — в случае одномерной области, когда правые части уравнений (1.1), (1.2) равны нулю (/(х,£) = 0), в [9] — в случае двумерной области при / (х,4) = 0.

Коэффициенты уравнения (2.1) как функции температуры терпят разрывы при Т = Т* и не определены в этой точке. Последний член правой части объединим с левой частью, тогда эффективная теплоемкость с — Ьфр имеет разрывы в точках Т = Т* и Т = Т* + Д. Исходя из сказанного, уравнение (2.1) представим в виде

дТ д дТ

коэффициенты которого определены во всем диапазоне изменения температуры следующим образом:

с2, если Т < Т* — А,

с2 — если Т* - А < Т < Т*,

, ~ (2.3)

0.5(с1 + с2) если Т = Т*,

сх, если Т > Т*,

Лх, если Т < Т* — А,

(Лх + Л2)/2 — (Л2 + Лх)(Т — Т*)/2А, если Т* — А < Т < Т* + А, (2.4)

Л2, если Т > Т* + А.

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

Таким образом, прямая задача (1.1)—(1.7) сводится к решению уравнения (2.2) с граничными условиями (1.5), (1.6) и с начальным условием (1.7).

При выборе функции г>(Т) необходимо учитывать направление хода процесса — фазовый переход происходит в результате понижения температуры или ее повышения. Она должна быть выбрана отличной от нуля в области вновь образующейся фазы.

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

3. Сведение обратной задачи к краевой задаче для нагруженного уравнения

Будем искать решение обратной задачи в виде

т (ж, г) = 0(гЖж) + ш(ж, г), (3.1)

где

t

0(г) = ^ ф) (3.2)

о

Подстановка (1.8), (1.9), (3.2) в (2.2) дает уравнение для ш(ж, г):

д д . . дЬ дх\ дх) дх\ дх)

С учетом (3.2) условие (1.9) приводит к следующему представлению для неизвестной 0(г):

6{1) = щ-^{ч>{1)-и{х*,1)). (3.4)

Подстановка (3.4) в (3.3) дает искомое нагруженное параболическое уравнение

Граничное условие с учетом наших предположений о правой части на границе области имеет вид

А^М = аЦ0,г)+Тс), = 0 < 4 < ¿о- (3.6)

Как следует из (3.2), для вспомогательной функции 0(г) имеем

0(0) = 0, (3.7)

что позволяет использовать начальное условие

ад(х, 0) = Тс, 0 <х<1. (3.8)

Тем самым обратная задача (2.2), (1.5)—(1.7) формулируется как краевая задача для нагруженного уравнения (3.5)—(3.8) с представлением (3.2), (3.4) для неизвестной зависимости источника от времени.

4. Разностная схема

Будем считать, что по переменной х введена равномерная сетка ш с шагом Л. Через XI = ¿Л, г = 0,1,..., N, ЖЛ = I, обозначим узлы сетки, и пусть V = V = г>(х^). Для простоты будем считать, что точка наблюдения х = х* совпадает с внутренним узлом, отвечающим номеру г = к.

При переходе с одного временного слоя , = ^'т, ] = 0,1,...,^°, т > 0, на другой временной слой ,+ 1 будем использовать чисто неявную разностную схему для уравнения (3.5). Во внутренних узлах сетки по пространству имеем

- Г , ,+ -, ч 1

С-

= (ао4+1)а! + —-4+1)(аф*)х. (4.1)

т у х фк

Для задач с достаточно гладким коэффициентом Л положим а = А • 0, 5(х^ + х^_1). Аппроксимируя (3.6), (3.8), получим

А"1 =*Ы+1-Тс),

г,+ 1 - ш,'+1

" 1г ДГ'1 =0, ^ = (4.2)

= Тс, г = 1, 2, ...,Ж - 1. (4.3)

Из решения разностной задачи (4.1)—(4.3) в соответствии с (3.4) определим

= .7 = 0,1,..., .70-1, (4.4)

дополнив эти соотношения условием (см. (3.7)) 0° = 0. Принимая во внимание (3.2), для искомой зависимости правой части от времени используем простейшую процедуру численного дифференцирования:

0,+ 1 — 0,

г]3+1 =-, .7 = 0,1,..., .70-1. (4.5)

т

Есть необходимость отдельно остановиться на вопросе решения сеточной задачи для реализации рассматриваемой неявной схемы.

5. Сеточная нелокальная задача и программная реализация

Никаких особых проблем с вычислительной реализацией схемы (4.1)—(4.3), несмотря на нестандартность (нелокальность) сеточной задачи на новом временном слое, не возникает. Следуя методике [5,7], запишем уравнение (4.1) во внутренних узлах в виде

г , + 1 1

^ - («-Г),, + ТкЫ*)*А+1 = я1 (5.1)

с заданной правой частью д? и граничными условиями (4.2). Решение системы (4.2), (5.1) ищется в виде

= Уг + г = 0,(5.2)

Подстановка (5.2) в (5.1) позволяет сформулировать следующие сеточные задачи для вспомогательных функций уг, 2г:

с--(аух)х,г = д1, г = 1,2,...,М-I, (5.3) т

\У1-Уо , гг, , \VN-VN-1 п /к ^

А—-— = а(у0-Тс), А---= 0, (5.4)

2г 1 с— - (а,г5)х^ + — (афх)хл = 0, г = 1, 2,..., Ж - 1, (5.5)

т фи

А - 20 А ^ - -1 „ /К „ ,

А—-— = а^о, —А---= I). (5-Ь)

После этого с учетом представления (5.2) находится

¿+1_ № к 1 - гк

(5.7)

Корректность алгоритма обеспечивается необращением знаменателя (5.7) в нуль. Для сеточной задачи (5.5), (5.6) на основе принципа максимума для разностных схем устанавливается априорная оценка

тах |.гг| < т тах

0<г<№ 0<г<N

1

фк

(афх)х,г.

Поэтому |2г| < 1 при достаточно малых т = 0(1), т. е. необходимо использовать малые шаги по времени.

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

1. По известному предыдущему распределению температуры выбирается параметр А для сглаживания разрывных коэффициентов си А так, чтобы промежуток сглаживания (Т* — А,Т*] содержал в себе значения температуры хотя бы в двух соседних узлах сетки, между которыми находится граница фазового перехода. Если в указанный интервал не попадает ни одно из узловых значений температуры, то температурное поле будет определено без учета выделения теплоты фазового перехода.

2. Вводится функция распределенного источника г>(Т), удовлетворяющая условиям утверждения, и с помощью нее определяются согласно формулам (2.3), (2.4) значения разрывных коэффициентов А и А.

3. Решаются вспомогательные сеточные задачи: для у — (5.3), (5.4), для 2 — (5.5), (5.6).

4. Далее последовательно определяются значения ^к+1 — по формуле (5.7), — по формуле (5.2), — по формуле (4.4).

5. Восстановление неизвестной зависимости правой части от времени п?^1 производится по формуле (4.5).

6. Распределение температурного поля для данного временного слоя находится по формуле (3.1) и присваивается значению предыдущего распределения температуры.

7. Процесс итерации (пп. 1-6) продолжаем до выполнения условия

Т(ж, г) - Т(ж,г)| < 5,

где 5 — уровень погрешности, Т(ж, г) — решение прямой задачи. Далее переходим на следующий временной слой.

По данному алгоритму проведены численные расчеты. В рамках концепции квазиреального эксперимента рассматривается прямая задача (2.2), (1.5)-(1.7) с некоторой заданной правой частью (1.8), где

( г, если 0 < г < 0.6, /пж

= 1 п ,ЧПЙ Щх) = 8111 Т~

{ 0, если г > 0.6, \ I

Рассмотрены два варианта выбора функции г>(Т):

1МТ) = 1 + ^,

2) и(Т) = ехр(0.69(Т* - Т + А)/А) - 1.

Такая задача численно решается при следующих значениях параметров: С1 = 2814 кДж/м3-К, Ах =6, 3 кДж/м-ч-К, с2 = 2016 кДж/м3-К, А2 = 8.4 кДж/ м-ч-К, Т* = 0° С, Тс = -25° С, Ь = 0.4175406 кДж/м3, а = 83, 5 кДж/м2-ч-К, I = 0.1 м, N = 100, к = 40, т = 0.0125.

Для решения обратной задачи (3.1)-(3.8) по значениям функции Т(ж, г) в точке ж* = 0.6, найденной решением прямой задачи, определяется функция дополнительного наблюдения ^(г).

Для проверки корректности разработанного алгоритма результаты расчетов обратной задачи сопоставлены с результатами прямой задачи.

В табл. 1 приведены результаты расчетов в узлах сетки при ^о = 50 (г = 0.625 ч). В первых двух строках сопоставлены значения температуры, полученные решением прямой задачи (Тх) и обратной задачи (Т2) при первом виде функции г>(Т); в строках 3, 4 — результаты прямой задачи (Тх) и обратной задачи (Т3) при втором виде функции г>(Т); в строках 5, 6 — результаты прямой задачи (Тх) и обратной задачи (Т4) при традиционном методе сглаживания.

Таблица 1

Узлы 0 5 10 20 30 40 50

1 ть °с -0.210 1.481 3.007 5.611 7.445 8.481 8.805

2 т2, °с -0.262 1.589 3.178 5.701 7.376 8.411 8.814

3 Ть °С -0.210 1.584 3.070 5.632 7.450 8.482 8.806

4 Тз, °с -0.261 1.593 3.184 5.711 7.385 8.416 8.816

5 ть °с -0.312 1.404 2.961 5.598 7.441 8.480 8.805

6 Та, °С -0.262 1.587 3.174 5.694 7.368 8.403 8.799

Из таблицы видно, что значения температуры, полученные традиционным методом сглаживания (строки 5 и 6), ниже значений, полученных предложенным методом введения распределенного источника тепла. Это обусловлено выбором величины объемной теплоемкости, которая при традиционном методе сглаживания принимается большей, чем ее действительное значение.

Рис. 1. Распределение температурного поля = 0.5)

На рис. 1 представлено распределение температурного поля, полученное решением обратной задачи при г>(Т) = 1 + в различные моменты времени,

6 — уровень погрешности.

Иллюстрацией корректности рассматриваемой обратной задачи и предлагаемого способа учета теплоты фазового перехода являются данные расчетов при различном уровне погрешностей. На рис. 2 слева дано решение задачи при 6 = 0.75, справа — при 6 = 0.25. При уменьшении уровня погрешности решение восстанавливается более точно.

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

ЛИТЕРАТУРА

1. Алифанов О. М. Обратные задачи теплообмена. М.: Машиностроение, 1988.

2. Бек Дж., Блакуэлл Б., Сент-Клэр Ч. Некорректные обратные задачи теплопроводности. М.: Мир, 1989.

3. Samarskii A. A., Vabishchevich P. N. Computational heat transfer. Chichester: Wiley, 1995.

4. Тихонов А. Н., Арсенин В. Я. Методы решения некорректных задач. М.: Наука, 1986.

5. Самарский А. А., Вабищевич П. Н. Численные методы решения обратных задач математической физики. М.: Изд-во ЛКИ, 2009.

6. Нахушев А. М. Уравнения математической биологии. М.: Высш. шк., 1995.

7. Борухов В. Т., Вабищевич П. Н. Численное решение обратной задачи восстановления источника в параболическом уравнении // Мат. моделирование. 1998. Т. 10, № 11. С. 93—100.

8. Павлов А. Р., Слепцова Е. А. Численное решение задачи Стефана введением распределенного источника // Мат. VI научно-техн. конф. «Современные проблемы теплофизики в условиях Крайнего Севера», посвящ. памяти проф. Н. С. Иванова. Якутск: ???, 2004. С. 74-79.

9. Слепцова Е. А. Математическое моделирование развития напряженно-деформированного состояния тонких пластин при их стыковой сварке: Дис. ... канд. физ.-мат. наук. — Якутск, 2009. — 109 с.

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

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

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

Статья поступила 25 февраля 2014 г. Крылова Екатерина Анатольевна

Северо-Восточный федеральный университет им. М. К. Аммосова, Институт математики и информатики ул. Кулаковского, 48, Якутск 677000 krekan2012@mail.ru