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

®

www.volsu.ru

DOI: https://doi.Org/10.15688/jvolsu1.2017.2.6

УДК 536.2.02 ББК 31.3

ЧИСЛЕННОЕ РЕШЕНИЕ НАЧАЛЬНО-КРАЕВЫХ ЗАДАЧ ДЛЯ УРАВНЕНИЯ ТЕПЛОПРОВОДНОСТИ МЕТОДОМ ИНТЕГРАЛЬНЫХ УРАВНЕНИЙ1

Анатолий Михайлович Афанасьев

Доктор технических наук, профессор кафедры информационной безопасности, Волгоградский государственный университет a.m. afanasiev@yandex. ru, infsec@volsu. ru

просп. Университетский, 100, 400062 г. Волгоград, Российская Федерация

Андрей Юрьевич Глухов

Кандидат технических наук, старший преподаватель кафедры радиофизики, Волгоградский государственный университет aglukhov@live.ru, rf@volsu.ru

просп. Университетский, 100, 400062 г. Волгоград, Российская Федерация

Борис Николаевич Сипливый

Доктор технических наук, профессор кафедры теоретической физики и волновых процессов, Волгоградский государственный университет Siplivyj@list.ru, tf@volsu.ru

просп. Университетский, 100, 400062 г. Волгоград, Российская Федерация

Аннотация. В статье рассмотрен метод решения начально-краевых задач для

уравнения теплопроводности в областях произвольной формы, основанный на сведе-

г- нии исходной задачи к решению интегро-дифференциального уравнения. Интегральный

оператор, порождающий это уравнение, обладает свойствами, позволяющими пред-

^ ставить решение в виде разложения по собственным функциям оператора. Такой под-

w ход позволяет естественным образом учесть краевые условия Дирихле, Неймана тре->0

й тьего рода.

| Ключевые слова: уравнение теплопроводности, функция Грина оператора Лап-

ласа, интегро-дифференциальное уравнение, собственные функции и числа, начальные

. краевые задачи.

<

« о

При численном исследовании процессов удаления влаги из влагосодержащих материалов ^ электромагнитным излучением возникает необходимость решения начально-краевых задач для S уравнения теплопроводности при различных краевых условиях. Для некоторых тел канонической геометрии удается получить аналитическое решение [1-3]. Для областей произвольной формы предлагается численный алгоритм, основанный на использовании функции Грина опера-« тора Лапласа.

^ Рассмотрим твердое тело с постоянным коэффициентом температуропроводности a и дей-© ствующими внутри источниками тепла с заданной плотностью f(M, t).

Температурное поле внутри тела определяется решением начально-краевой задачи [7]:

1 дТ

АТ---= f,

а д

Т|5 =ц(М,X), М е 5, (1)

ТМ,0) = Ч>(М), М е V,

здесь 5 - поверхность тела; V - объем. Запишем уравнение в виде:

1 лт

АТ = f +--= F (2)

а дх

Пусть известна функция Грина G(M, X) оператора Лапласа для области V (способ приближенного построения этой функции будет указан ниже). Тогда решение уравнения представимо в виде [7]:

1 1

Т (М, X ) = ~— ( F (М, X ЩМ, N ^- — Ш, X)— (М, N )05„. 4% V 4я 5 дпм

Подставляя сюда вместо F(M, X) его значение из (2), получаем: 1 1

Т(М,X)+ —(М,XЩМ,N=-f (М,XЩ(М,N-

дX 4я V

1 Л/^

(ф, X)— (М, N ^ = g(M, X).

1тг Ип

(3)

4я5 дПИ

Уравнение (3) есть интегро-дифференциальное уравнение, описывающее температурное поле Т(М, X) в объеме V. К нему необходимо добавить начальное условие Т(М, 0) = у(М).

Особенностью уравнения является независимость ядра интегрального оператора от времени и его симметрия, что позволяет применить для построения решения теорему Гильберта -Шмидта [4]. В соответствии с этой теоремой решение и правую часть уравнения представим в виде разложения по собственным функцияму^(М) интегрального оператора, порождающего уравнение (3):

Т(М, X )=X ^ (X )/к (М ) (4)

к '

g (М, X ) = Х ак (X )/к (М )

к '

где/к(М) - нормированное решение уравнения:

/(М ) = ^{ /^ЩМ, N №.

V

Подставляя (4) в (3), получим:

х Ск т (м)+х ^ (РЩМ, N)dVN =

= Х ак (X )/к (М ),

к

где .) «V' Г' N -

V

ак (' )= | Я (N, ?) / (N )с^ = ак (')+ а^ (?),

V

а? (0 = ^ ^ ,? СSN

кл'

ОК (2) .

к л V

а(2)

В силу линейной независимости системы собственных функций/(М имеем:

ск (0 + ^ (х) = ак (0 или

4акла а?

Сс

—к (?)+ к ласк (О = к л аак (0,

с

Ск(0) = Ьк = ^)/к^% .

V

Решение этой задачи Коши представимо в виде:

Ск (?) = 4^к ла 1 ак (тЩ(? - т)Ст + ^е^'

'•■к

0

Удовлетворяя начальному условию, имеем Ск = Ък. Здесь Щ( ') - импульсная реакция дифференциального оператора, определяемая решением задачи Коши:

—(')+ 4л^ каЩ = 0 Ж

щ(о)=1

Щ(? ) = е .

Подставляя найденные ск(') в (4), получим окончательный вид решения:

Т (М,' ) = 4ла X / (М)

ч[ак (т)е-4"кла -т)Ст + Ъке л

0

(5)

где

V

(1) (2)

(') = 1 g N, О/к Ис^ = ак (0+ ак (0, Ък =1^)/ (N % .

V

Для построения решения по этой формуле необходимо знание функции Грина G(M, N1, собственных функций /(М и собственных чисел А^.

Приближенное построение функции Грина

По определению функции Грина задача Дирихле для уравнения Лапласа в области V, ограниченной поверхностью имеет вид [7]:

G (М, N ) = — + и (М, N )

к

G (М, N ) = 0 при N е 5, АЩ (М, N )= 0.

Задача построения Щ(М, Щ заключается в отыскании гармонической в V функции и(М, Ж), обеспечивающей выполнение граничного условия Щ(М, N = 0, при N е 5.

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

Будем искать функцию Грина в виде:

Щ(М, N ) = —

(

\

■ +... + -

МЫ у М1Ы М гЫ МпА у

здесь М1 - точки расположения источников вне V; - мощности источников.

Совмещая точку N последовательно с граничными точками Q i =1, 2 ... т, потребуем выполнения граничного условия:

Щ(М, Q¡ ) = 0, i = 1,2... т. Эти соотношения задают систему алгебраических уравнений относительно с. :

где

а11 а12

а21 а21

м м

ат1 ат1

1

а.. = -

¡1 Г ^М

Л Л

2п

м

Л а„

С1 Ь

С2 Ь2

м м

Сп Ь т

Ь =--

М е¥

MQ¡

Заметим, что при т > п система не совместна, тогда ее можно решить приближенно, например, методом наименьших квадратов [6]. Для исключения такого случая целесообразно положить п = т, в результате получим систему с квадратной матрицей. Такой способ вычисления функции Грина можно рассматривать как распространение метода электростатических изображений [8] на области произвольной формы.

Второй способ приближенного построения функции Грина оператора Лапласа основан на использовании потенциала двойного слоя. Гармоническую в области функцию ЩМ, N), входящую в выражение для функции Грина:

Щ (М, N ) = — + Щ (М, N )

можно искать в виде потенциала двойного слоя:

Щ (М, N ) = )

д 1

4- 5

05,

Q .

Удовлетворяя граничному условию:

Щм {М, N } М% =-±

^MеV Г

• Лл

АеЯ MеV

с

с

с

2

п

+

а

1

и используя свойства потенциала двойного слоя, получим интегральное уравнение:

8 1

т^)-4"1т(б)8-—dSв = —, N е s, М е v,

4л х 8пв Гив М

здесь п - внешняя к S нормаль.

Для каждой точки М е V, необходимой для дальнейшего использования функции Грина G(M, N1, решается это ИУ и вычисляется функция Грина.

Если решение искать в виде потенциала простого слоя:

и (м, n ) = ±- ^ dSв, 4 л г

то приходим к интегральному уравнению первого рода:

dSв =——, N е s, М Еv

г в г

Я 'ив МЫ

Вычисляя приближенно интеграл по формуле прямоугольников, получим:

С*в СЯв -I» 1 ^ = - к = 1.2...П .

я N • АЯ, ^кв г АЯ, Хв 'шк

Эти соотношения задают систему линейных алгебраических уравнений (СЛАУ) относительно о:

г

1 г СЯв

1СТгагк =--, а1к = 1 к = 1,2...П

^^ Г * V

АЯ, ^кв

Из-за наличия особенностей у^^ диагональные элементы матрицы этой СЛАУ превосходят остальные, поэтому определитель СЛАУ не равен нулю и она разрешима и при этом единственным способом.

Синусоидальный режим

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

Т(М )+л[ ^(М, N )С¥м = ^М ), (6)

V

. гга

где л =-; ю - частота,

4ла

1 1 л/^

СМ) = - —1 /№(М, N)dVN - ^Ф^)— (М, N)dSN.

Представляя Т(М) и СМ) в виде разложений по собственным функциям ядра и подставляя эти разложения в (6), получим:

ТМ)=1С/ (М), ^М)=£ Щ/ (М),

кк

I Ск/к (М )+£А-С/ (М )-1 Щ/к (М )=0.

к к л к к

Отсюда имеем:

( \ \

411+ ХЧ = ],

У х к)

4

к

1 +

х

х^

где

к)

= 4 + 4

■'к ~ I и к V" ' М ~ ^к т к 1

( 4м ) / (М )0УМ

V

4 = ^^ )£/к N ^,

4 = ^ /¥) /к ^ )цуя.

Подставляя найденные 4 в решение, получим:

4М)/ (М )=Ц 4е318 С4/к (М). к 1 + _ к

1к

хк

Соответствующая временная функция имеет вид:

Т (М, X) = XI4/ (М >ш (ш + 3Г8 4). (7)

к

Решение задачи при краевых условиях Неймана

Пусть в задаче (1) краевые условия имеют вид:

дТ (м , X)

дп

= ф(М, X), М е 5 .

Решение этой задачи может быть выражено через функцию Грина задачи Неймана [8]: Т(М,X) = -±1 F(N, X)GN (М,N)dVN - 4-(дТNXУЩ" (М,N)dVN +

+ 4-

¡-(ф^, X Щ (М, N )05м

где

(М, N) = — + Щ (М, N),

дЩА

АЩ (М, N ) = 0, (М, N ) =--—, N е 5.

дпн 4-5

Эта функция определяется с точностью до слагаемого с0(М), зависящего от положения точки М.

Нормирующее условие ( GN (М, N )d5N = 0 определяет функцию Щ^(М, N однозначно и де-

5

лает ее симметричной [8]. Интегро-дифференциальное уравнение, аналогично (3), имеет вид:

Т (М, ?)+ 4а 18Т N, ? ^ (М, N )СУи =- 4-1 / (N, ? (М, N

1-fф(N, ? )GN (М, N)СSN = g (М, ?).

4л о

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

Т (М,' ) = 1 / (М)

4лаЛ к 1 ак (т)е^к 4-('-т)Ст + Ъке

Ък = 1^ )/к N )СУн ,

V

ак = ак(')+ ак('), (? )= 4Л!г1ф(N, ? )/к (N ^ ,

!(?) = - 4-^1 /^, ?)/к (N^ .

/ N), л к -

решение уравнения

нормированное условием

/ (М )=л1 / N (М, N )СУн ,

1 / ^ ^ = 0.

Приближенное построение функции Грина задачи Неймана

Как и ранее, функцию Грина GN(M, Щ будем искать в виде:

GN(М,N) = + + ^ +... + ^ + с0(М), N е V г г Г г

'ММ ' М1И ' М2И ' МпМ

(8)

Здесь с0 - постоянная, зависящая от точки М, с точностью до которой определяется GN(M, N1. Эта постоянная однозначно определяется нормирующим условием:

GN (М, N= 0,

,(М ) = - ^ 1

^+1-

'МУ к=1 'М^

Подставляя GN(M, N в граничное условие:

8п„

(М, в, )=--Ц, в, Е Я, , = 1,2К

4-Я

да,

получим СЛАУ относительно ск. В матричном виде эта СЛАУ имеет вид:

с

к

с

а. = -

У

а11 а12

а21 а21

м м

ат1 ат1

д ( 1

дпа г {'а, МУ

л

Л

л

ь. = -

а1п с1 ь

а2п С2 ь2

м м м

а тп Сп ь т

д ( 1 \

дпг

ш

4^

Подстановка найденных ск и с0 в (8) дает приближенное значение функции Грина задачи Неймана для уравнения Лапласа.

Приближенные значения собственных функций/(М и собственных чисел X к интегральных

операторов с ядрами G(M, N и GN(M, N можно построить методом Келлога [4]. Другой способ вычисления собственных функций и чисел заключается в замене в (5) интеграла квадратурной формой (например, прямоугольников). При этом задача отыскания собственных функций и чисел сводится к алгебраической проблеме вычисления спектра симметричной матрицы, которая может быть решена методом спектрального разложения матрицы [6].

Из формул (6) и (7) следует, что искомое температурное поле Т(М, X) определяется только собственными функциями/(М и числами \к удовлетворяющими уравнению (5). Известно [8], что решение (5) эквивалентно решению задачи Штурма - Лиувилля:

А/ + \/ = 0,

(8)

/ (М)

= 0.

Эту задачу можно приближенно решить методом, аналогичным методу построения функции Грина G(M, N1, изложенному выше.

При положительных X фундаментальное решение уравнения (8) имеет вид [7]:

R(M, N ) =

¡ш^М)

Вне области решения Vвыберем произвольно точки N 1, . = 1,2...п и расположим в них источники мощностью с.. Решение уравнения (8) представим в виде наложения полей этих источников:

/ М, N ) = X

¡ш^М, )

Последовательно совмещая точку М с граничными точками а.., . = 1,2...п и приравнивая результат к нулю, получим однородную СЛАУ с симметричной матрицей А = (а.), где

а у =■

¡ш )

ненулевое решение этой СЛАУ существует при условии det(a У) = 0. Корни этого уравнения есть искомые собственные числа Хк Метод решения этого уравнения и вычисления собственных функций изложен в [5].

Заметим, что при таком способе определения собственных чисел и функций не требуется предварительного вычисления функции Грина G(M, И), что сокращает объем вычислений.

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

1

< V

г

г

1=1

г

МЛ

^^^^^^^^^^^^^^^^^^^^ КОМПЬТЕРНОЕ МОДЕЛИРОВАНИЕ

ПРИМЕЧАНИЕ

1 Исследование выполнено при финансовой поддержке РФФИ и Администрации Волгоградской области в рамках научного проекта «16-48-340527 р_а».

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

1. Афанасьев, А. М. Задача о сушке шара электромагнитным полем / А. М. Афанасьев, Б. Н. Сипливый // Инженерно-физический журнал. - 2013. - N° 5. - С. 3-8.

2. Афанасьев, А. М. Теория электромагнитной сушки: асимптотическое решение начально-краевой задачи для прямоугольной области / А. М. Афанасьев, Б. Н. Сипливый // Физика волновых процессов и радиотехнические системы. - 2012. - Т. 15, № 1. - С. 77-83.

3. Афанасьев, А. М. Теория электромагнитной сушки: асимптотическое решение начально-краевой задачи для цилиндра / А. М. Афанасьев, Б. Н. Сипливый // Теоретические основы химических технологий. -2014. - Т. 48. - С. 222.

4. Васильева, А. Б. Интегральные уравнения / А. Б. Васильева, Н. А. Тихонов. - М. : Изд-во МГУ, 1989.- 157 с.

5. Князев, С. Ю. Применение метода точечных источников поля при численном решении задач на собственные значения для уравнения Гельмгольца / C. Ю. Князев, Е. Е. Щербакова // Известия высших учебных заведений. Электромеханика. - 2016. - № 3. - С. 11-17.

6. Малышев, А. Н. Введение в вычислительную линейную алгебру / А. Н. Малышев. - Новосибирск : Наука. Сиб. отд-ние, 1991. - 229 с.

7. Тихонов, А. Н. Уравнения математической физики / А. Н. Тихонов, А. А. Самарский. - 7-е изд. - М. : Наука, 2004. - 798 с.

8. Функция Грина оператора Лапласа / А. Н. Боголюбов, Н. Т. Левашова, И. Е. Могилевский, Ю. В. Мухта-рова, Н. Е. Шапкина. - М. : Физ. фак. МГУ 2012. - 130 с.

REFERENCES

1. Afanasyev A.M., Siplivyy B.N. Zadacha o sushke shara elektromagnitnym polem [The Problem of Ball Drying in Electromagnetic Field]. Inzhenerno-fizicheskiyzhurnal, 2013, no. 5, pp. 3-8.

2. Afanasyev A.M., Siplivyy B.N. Teoriya elektromagnitnoy sushki: asimptoticheskoe reshenie nachalno-kraevoy zadachi dlya pryamougolnoy oblasti [The Theory of Electromagnetic Drying: Asymptotic Solution of the Initial Boundary Value Problem for the Rectangular Area]. Fizika volnovykh protsessov i radiotekhnicheskie sistemy, 2012, vol. 15, no. 1, pp. 77-83.

3. Afanasyev A.M., Siplivyy B.N. Teoriya elektromagnitnoy sushki: asimptoticheskoe reshenie nachalno-kraevoy zadachi dlya tsilindra [The Theory of Electromagnetic Drying: Asymptotic Solution of the Initial Boundary Value Problem for a Cylinder]. Teoreticheskie osnovy khimicheskikh tekhnologiy, 2014, vol. 48, p. 222.

4. Vasilyeva A.B., Tikhonov N.A. Integralnye uravneniya [Integral Equations]. Moscow, MGU Publ., 1989.

157 p.

5. Knyazev S.Yu., Shcherbakova E.E. Primenenie metoda tochechnykh istochnikov polya pri chislennom reshenii zadach na sobstvennye znacheniya dlya uravneniya Gelmgoltsa [Applying the Point Sources Method of the Field in the Numerical Solution of Eigenvalue Problems for the Helmholtz Equation]. Izvestiya vysshikh uchebnykh zavedeniy. Elektromekhanika, 2016, no. 3, pp. 11-17.

6. Malyshev A.N. Vvedenie v vychislitelnuyu lineynuyu algebru [Introduction to Numerical Linear Algebra]. Novosibirsk, Nauka. Sibirskoe otdelenie Publ., 1991. 229 p.

7. Tikhonov A.N., Samarskiy A.A. Uravneniyamatematicheskoy fiziki [Equations of Mathematical Physics]. Moscow, Nauka Publ., 2004. 798 p.

8. Bogolyubov A.N., Levashova N.T., Mogilevskiy I.E., Mukhtarova Yu.V., Shapkina N.E. Funktsiya Grina operatoraLapplasa [Green's Function of the Laplace Operator]. Moscow, Fiz. fak. MGU Publ., 2012. 130 p.

NUMERICAL SOLUTION OF INITIAL BOUNDARY VALUE PROBLEMS FOR THE HEAT EQUATION BY THE METHOD OF INTEGRAL EQUATIONS

Anatoliy Mikhaylovich Afanasyev

Doctor of Technical Sciences, Professor, Department of Information Security,

Volgograd State University

a.m. afanasiev@yandex. ru, infsec@volsu. ru

Prosp. Universitetsky, 100, 400062 Volgograd, Russian Federation

Andrey Yuryevich Glukhov

Candidate of Technical Sciences, Senior Lecturer, Department of Radio Physics, Volgograd State University aglukhov@live.ru, rf@volsu.ru

Prosp. Universitetsky, 100, 400062 Volgograd, Russian Federation

Boris Nikolaevich Siplivyy

Doctor of Technical Sciences, Professor, Department of Theoretical Physics and Wave Processes, Volgograd State University Siplivyj@list.ru, tf@volsu.ru

Prosp. Universitetsky, 100, 400062 Volgograd, Russian Federation

Abstract. Sometimes it becomes necessary to solve initial boundary value problems for the heat equation under various boundary conditions in the numerical study of the processes of moisture removal from moisture-containing materials by electromagnetic radiation. For some bodies of canonical geometry, one can obtain an analytic solution. For domains of arbitrary shape, a numerical algorithm based on the use of the Green's function of the Laplace operator is proposed.

The study considers a solid body with a constant coefficient of thermal diffusivity a and internal heat sources with a given density fM, t). Consequently, the solution of the initial boundary value problem for a body with a surface S and volume V determines the temperature field inside the body.

If the Green's function G(M, N) of the Laplace operator is known for one of the listed boundary conditions, then the desired solution can be represented in the form (for boundary conditions of the first kind a = 1, P = 0).

As a result, the relation is an integral-differential equation with respect to T(M, t). The integral operator generating this equation is self-adjoint in the class of quadratically-integrable functions. Therefore, the Hilbert-Schmidts theorem is applicable for the solution. Expanding the functions T(M, t), f(M, t) into series in eigenfunctions of the integral operator obtains the Cauchy problem for the coefficient ck(t) of the expansion T(M, t).

Carrying out the necessary transformations, the solution of the starting initial boundary value problem is expressed in terms of eigenfunctions and the number of the kernel G(M, N) which can be calculated by the Kellogg method. The authors propose the certain form of the kernel G(M, N) which is necessary for the realization of this method.

Thus, putting N consecutively points into the boundary points and equating the results to zero it obtaines a system of linear algebraic equations for c..

As a result, the proposed method for constructing G(M, N) can be treated as a generalization of electrostatic images method to regions of arbitrary shape.

In addition, the study obtains expressions for boundary conditions of the second and third kind. Thus, the proposed method makes it possible to take into account boundary conditions of various types in a natural way.

Key words: heat equation, Green's function of the Laplace operator, integral-differential equation, proper functions and numbers, initial boundary value problems.