Численное и аналитическое исследование некоторых процессов, описываемых нелинейным уравнением теплопроводности Текст научной статьи по специальности «Математика»

____________УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО УНИВЕРСИТЕТА

Том 157, кн. 4 Физико-математические науки

2015

УДК 563.2:517.927.4

ЧИСЛЕННОЕ И АНАЛИТИЧЕСКОЕ ИССЛЕДОВАНИЕ НЕКОТОРЫХ ПРОЦЕССОВ, ОПИСЫВАЕМЫХ НЕЛИНЕЙНЫМ УРАВНЕНИЕМ ТЕПЛОПРОВОДНОСТИ

А.Л. Казаков, Л.Ф. Спевак

Аннотация

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

Ключевые слова: уравнения с частными производными, нелинейное уравнение теп-лопроводости, метод граничных элементов, степенной ряд, тепловая волна.

Введение

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

Tt = di v(k VT). (1)

Здесь T - искомая функция (температура), t - время, div, и V - дивергенция и градиент по пространственным координатам соответственно.

Наиболее распространенным и хорошо изученным является линейный вариант уравнения (1), когда k = const > 0 (и тогда div(kVT) = kAT). Однако линейная модель в ряде случаев является недостаточно точной для адекватного описания реальных процессов. В этой связи актуальным является исследование нелинейного уравнения теплопроводности, в котором коэффициент теплопроводности k зависит от искомой функции T. Одним из наиболее распространенных вариантов такой зависимости является степенная k(T) = aTа. В этом случае уравнение (1) путем замены искомой функции u = Tа и растяжения по пространственным координатам может быть записано в виде

ut = uAu 4—(Vu)2. (2)

а

42

ИССЛЕДОВАНИЕ НЕЛИНЕЙНОГО УРАВНЕНИЯ ТЕПЛОПРОВОДНОСТИ 43

Уравнение (2), в частности, описывает фильтрацию идеального политропного газа в пористой среде, и тогда u - плотность, а а > 0 - показатель политропы (адиабаты) газа [3]. В англоязычной литературе уравнение (2) обычно называют “the porous medium equation” [4].

В настоящей работе представлены одномерные краевые задачи для нелинейного уравнения теплопроводности в случае степенной зависимости коэффициента теплопроводности от температуры. Проведено обобщение результатов авторов, полученных ранее [5-8]. Представлены и верифицированы новые точные решения в виде кратных степенных рядов, а также численные решения с помощью метода граничных элементов (МГЭ).

1. Основной результат

Уравнение (2) в случае одной пространственной координаты имеет вид

1 2 , v

ut = uurr +-u +— ur

a r

Здесь u(t, r) - искомая функция; t, r - независимые переменные, r > 0; a > 0, v £ No - константы.

Физический смысл переменной r следующий: расстояние до начала координат в пространстве переменных xi, Х2,..., xv+i в евклидовой метрике, то есть r = (х2 + + Х2 + • • • + xv+i )1/2 . При этом наиболее содержательными являются случаи v = 0 (плоскосимметричный), v =1 (цилиндрически симметричный) и v = 2 (сферически симметричный). Поскольку случай плоской симметрии ранее был подробно изучен [5, 6], далее считаем, что v > 0.

Для уравнения (3) рассматриваются следующие краевые условия:

u |r=a(t) u0 (t, r) |r=a(t) g (t) ,

предполагается, что

a(0) = R> 0, g(0) = 0, g'(0) > 0, [g'(0)j2 + K(0)]2 > 0. (5)

При выполнении краевых условиях (4), (5) в момент времени t = 0 при r = R в уравнении (2) коэффициент перед старшей производной обращается в нуль. Данная особенность приводит к тому, что решение задачи (3), (4) может иметь вид тепловой волны, распространяющейся по невозмущенному фону с конечной скоростью.

Для задачи (3), (4) построено новое точное решение в виде сходящегося кратного степенного ряда вида

(4)

(3)

u

Е

k,l =0

uk,itk [r - a(t)]1 kill

(6)

Процедура построения коэффициентов ряда (6) представляет собой перенос на несколько более общий случай подходов, использованных ранее в работах авторов [5, 9], и здесь излагается кратко.

Из краевых условий имеем, что

ui,o = 9 (0), uo,i

aa'(0)

2

±

[aa/(0)]2

4

+ ag/(0).

Условие (5) обеспечивает то, что ui,o, uo,i одновременно в нуль не обращаются. Выбор знака перед корнем определяет направление движения тепловой волны

44

А.Л. КАЗАКОВ, Л.Ф. СПЕВАК

Рис. 1. Сравнение решений МГЭ и в виде ряда

(внутрь области, ограниченной поверхностью r = R, или наружу). Последующие коэффициенты ряда (6), сумма индексов которых n = k + l одинакова, определяются индукцией по n при решении трехдиагональных систем линейных алгебраических уравнений, размерность которых с увеличением n неограниченно возрастает и для которых не выполнено условие диагонального преобладания. Сходимость рядов (6) доказывается методом мажорант. Показано, что построенное точное решение задачи (2), (4) имеет вид тепловой волны (см. выше).

Различные варианты задачи (2)-(4) ранее встречались в литературе [3], однако в данной постановке, насколько известно авторам, она впервые рассмотрена в [5] (см. также [6]) для случая v = 0. Таким образом, настоящее исследование обобщает ранее полученные результаты авторов.

2. Вычислительный эксперимент

Численные решения МГЭ для v = 0 представлены в работе [6]. Алгоритмы решения МГЭ задачи (2), (4) для случаев v =1, 2 при условии a (t) = R описаны в [7, 8]. В настоящей работе построены аналогичные алгоритмы решения при g (t) = = 0. Построенные в работе аналитические и численные решения при различных краевых режимах сравниваются друг с другом и с точными решениями [10].

На рис. 1 показано сравнение решения в виде отрезка степенного ряда (до 4-й степени включительно) и решения МГЭ при v =1, а = 3, a (t) = 1 +t, g (t) = 0.

( 16 y/8

На рис. 2 решение МГЭ при v = 1, а = 3, a (t) = I 1 + — t I , g (t) = 0

сравнивается с точным решением

и

- 16

1 + Yt

+

2

r

1

„ 16 1 + Yt

3/4 '

На рис. 3 решение МГЭ при v = 2, а = 3, a (t) сравнивается с точным решением

г2 1

1 11 \1/111 ^

1 + уч , g(t)

л 22 + ( 22 \ 9/11 •

4 + Тt 42/11 4 + 22 Л

и = —

0

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

ИССЛЕДОВАНИЕ НЕЛИНЕЙНОГО УРАВНЕНИЯ ТЕПЛОПРОВОДНОСТИ 45

Рис. 2. Сравнение решения МГЭ с точным решением, цилиндрическая симметрия

3. Построение фронта тепловой волны

Ряд (6) удобен для исследования (в том числе численного) решения задачи (2), (4) вблизи линии r = a(t). Вместе с тем он не позволяет эффективно строить фронт тепловой волны, то есть линию r = b(t), разделяющую на плоскости Ort области возмущенного решения и фонового решения, которое в данном случае является тривиальным. В самом деле уравнение и = 0 в случае, когда функция u(t, r) представлена в виде ряда с рекуррентно вычисляемыми коэффициентами, решать затруднительно. В данном разделе предлагается процедура, которая позволяет построить функцию b(t) в виде степенного ряда.

Сделаем в уравнении (3) замену переменных

s = t,

У = r - b(t).

(7)

Замена (7) интересна тем, что функция b(t) пока неизвестна, однако из (4) следует, что b(0) = 0.

Задача (3), (4) в результате замены примет вид

{. , / \ 1 2 UUy

Us = b (s)uy + UUyy + ~Uy + y + b(s) ,

U\y=0 0, U\y=a(s)-b(s) f (s), b(0) 0*

(8)

46

А.Л. КАЗАКОВ, Л.Ф. СПЕВАК

Рис. 4. Движение нулевого фронта

Решение задачи (8) будем искать в виде

ж „,1

‘(*,у)= Е Vk’1 ТТ тт> Vk’1

k,l = о

gk+it

dsk dyl

s=0

y=0

dnb

b(s) = V bn ^r, bn = d^b w ^ n.r n dsn

n= 0

co

О

n

0

s

(9)

Коэффициенты рядов (9) определяются индукцией по суммарному порядку дифференцирования n = k + I.

Из граничных условии имеем, что b0 = 0, Vk,0 =0, k = 0,1,..., в том числе V0,0 = 0, vi,0 = 0. Для определения V0,i, bi, положив в (8) s = у = 0, получим систему алгебраических уравнений

biV0,i + - V2 1 = 0, а

V0,i(ai - bi) = fi,

(10)

Решение (10) имеет вид

V0,i

aai 1 2

(а ai)2 4

+ аfl,

bi

V0,i

а

При этом выбор знака в последней формуле соответствует выбору направления движения тепловой волны (к особой точке r = 0 или во внешнюю область).

Для определения коэффициентов

bn, Vn-i,b Vn-2,2, • • • , V0,n (11)

при каждом n решается система линейных алгебраических уравнений с двухдиагональной матрицей и отличным от нуля определителем. Отыскание коэффициентов (11) по описанной процедуре соответствует построению рядов (9), в том числе функции b(t).

Замечание. Из (10) следует, что кривые r = a(t) и r = b(t) при условии f '(0) = 0 имеют общую касательную.

Отметим, что разработанный авторами алгоритм решения задачи (3), (4) при условии a(t) = R методом граничных элементов также включает в себя построение

ИССЛЕДОВАНИЕ НЕЛИНЕЙНОГО УРАВНЕНИЯ ТЕПЛОПРОВОДНОСТИ 47

фронта тепловой волны на каждом шаге по времени [7, 8]. На рис. 4 представлено сравнение фронта волны, полученного при помощи алгоритма на основе МГЭ, с точным решением из примера 2 (рис. 2). Соответствие точного решения и отрезка ряда до слагаемых 4-го порядка также хорошее (график здесь не приводится, так так практически совпадает с графиком МГЭ).

Работа выполнена при частичной поддержке Программы фундаментальных научных исследований УрО РАН (проект № 15-7-1-17).

Summary

A.L. Kazakov, L.F. Spevak. Numerical and Analytical Study of Processes Described by the Nonlinear Heat Equation.

The paper deals with one-dimensional boundary value problems for the nonlinear heat conduction (porous medium) equation in the case of a power law relationship between the heat conduction coefficient and temperature, the solution of which is a heat wave spreading along the cold background at the finite rate. A numerical algorithm of the solution obtained by the boundary element method is developed. Its software implementation is performed. Numerical calculations are carried out. Their results are compared to the known exact solutions of the considered equation, as well as to the multiple power series that are solutions of the investigated problem in case of the analyticity of the input data. The results reported here extend our previous work in this area.

Keywords: partial differential equations, nonlinear heat conduction equation, boundary element method, power series, heat wave.

Литература

1. Тихонов А.Н., Самарский А.А. Уравнения математической физики. - М.: Наука, 1977. - 735 с.

2. Баренблатт Г.И., Ентов В.М., Рыжик В.М. Теория нестационарной фильтрации жидкости и газа. - М.: Недра, 1972. - 288 с.

3. Сидоров А.Ф. Избранные труды: Математика. Механика. - М.: Физматлит, 2001. -576 c.

4. Vazquez J.L. The Porous Medium Equation: Mathematical Theory. - Oxford: Clarendon Press, 2007. - 648 p.

5. Казаков А.Л., Лемперт А.А. О существовании и единственности решения одной краевой задачи для параболического уравнения типа нестационарной фильтрации // Прикл. механика и техн. физика. - 2013. - Т. 54, № 2. - С. 97-105.

6. Kazakov A.L., Spevak L.F. Numerical and analytical studies of a nonlinear parabolic equation with boundary conditions of a special form // Appl. Math. Model. - 2013. -V. 37, No 10-11. - P. 6918-6928.

7. Казаков А.Л., Спевак Л.Ф. Численное исследование нелинейного уравнения теплопроводности в случаях круговой и сферической симметрии // Вестн. КГТУ им. А.Н. Туполева. - 2013. - № 3. - С. 105-110.

8. Казаков А.Л., Кузнецов П.А., Спевак Л.Ф. Об одной краевой задаче с вырождением для нелинейного уравнения теплопроводности в сферических координатах // Труды Ин-та математики и механики УрО РАН. - 2014. - № 1. - С. 119-129.

9. Казаков А.Л., Кузнецов П.А., Лемперт А.А. О построении тепловой волны для нелинейного уравнения теплопроводности в симметричном случае // Изв. Иркут. гос. ун-та. Сер. Матем. - 2015. - Т. 11. - С. 39-53.

48

А.Л. КАЗАКОВ, Л.Ф. СПЕВАК

10. Полянин А.Д., Зайцев В. Ф. Справочник по нелинейным уравнениям математической физики: Точные решения. - М.: Физматлит, 2002. - 432 с.

Поступила в редакцию 29.06.15

Казаков Александр Леонидович - доктор физико-математических наук, заведующий лабораторией математических методов анализа свойств динамических систем, Институт динамики систем и теории управления имени В.М. Матросова СО РАН, г. Иркутск, Россия.

E-mail: kazakov@icc.ru

Спевак Лев Фридрихович - кандидат технических наук, заведующий лабораторией прикладной механики, Институт машиноведения УрО РАН, г. Екатеринбург, Россия.

E-mail: lfs@imach.uran.ru