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

УДК 517.958+519.633

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

И. В. Кудинов, В. П. Радченко

Самарский государственный технический университет,

443100, Самара, ул. Молодогвардейская, 244.

E-mail: igor_koudinov@mail.ru; radch@samgtu.ru

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

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

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

Ниже изложен метод, относящийся к группе интегральных методов теплового баланса и позволяющий получать аналитические решения нелинейных задач нестационарной теплопроводности с достаточно высокой точностью практически во всем диапазоне времени нестационарного процесса (0^Ро<то) без каких-либо ограничений на величину числа Фурье в области малых его значений. Основную идею метода рассмотрим на примере решения нелинейной задачи при линейной зависимости коэффициента температуропроводности от температуры а(Т) = ао(1 + /ЗТ). Математическая постановка задачи для бесконечной пластины при граничных условиях первого рода имеет вид:

9Є(Є, Fo) <9Fo

д_

1 + в (0 + Тс)

д0(С, Fo)

д£

0(С, О) = О; 0(О, Fo) = AT; д0(1, Fo)

дС

= О,

(Fo > О; О <С < 1); (1)

(2)

(3)

(4)

Игорь Васильевич Кудинов, аспирант, каф. прикладной математики и информатики. Владимир Павлович Радченко (д.ф.-м.н., проф.), зав. кафедрой, каф. прикладной математики и информатики.

«2 (Го)

Рис. 1. Расчётная схема теплооб-

где Т — температура; Т0 — начальная температура; 0 = (Т — То) — избыточная температура; Ро = ат/Л2 — число Фурье; £ = ж/Л — безразмерная координата; ж — координата; Я — половина толщины пластины; ао — коэффициент температуропроводности; Тст — температура стенки; в — коэффициент; АТ = Тст — Т0.

Процесс нагрева разделим на две стадии по времени: 0 < Ро ^ Рої и Рої ^ Ро < то. Для этого введём движущуюся во времени границу (фронт температурного возмущения), разделяющую исходную область 0 ^ ^ 1 на две

подобласти: 0 ^ ^ ?1(Ро) и ^1(Ро) ^ ^ 1,

где 51 (Ро) —функция, определяющая продви- мена

жение границы раздела во времени (см. рис. 1).

При этом в области, расположенной за фронтом температурного возмущения, сохраняется начальная температура. Первая стадия процесса заканчивается при достижении движущейся границей центра пластины (£ = 1), т. е. когда Ро = Рої. Во второй стадии изменение температуры происходит по всему объему тела (0 ^ ^ 1). Здесь вводится в рассмотрение допол-

нительная искомая функция 52(Ро) = 0(1, Ро), характеризующая изменение температуры во времени в центре пластины.

Математическая постановка краевой задачи для первой стадии процесса в результате введения фронта температурного возмущения 5ї(Ро) будет включать уравнение (1) с граничным условием (2), а также следующие граничные условия, выполняющиеся на фронте температурного возмущения:

0 (5ї, Ро) = 0; д0 (5ї, Ро)

д£

(5)

(6)

где соотношения (5), (6) представляют условия сопряжения прогретой и непро-гретой зон. Соотношение (5) устанавливает равенство температуры тела в точке £ = 5ї(Ро) его начальной температуре. Согласно условию (6) тепловой поток не распространяется за пределы фронта температурного возмущения (условие адиабатной стенки).

Обратим внимание на тот очевидный факт, что на первой стадии процесса задача (1), (2), (5), (6) за пределами фронта температурного возмущения, т. е. на отрезке 5ї(Ро) ^ ^ 1, вообще не определена. В связи с этим здесь

нет необходимости выполнения начального условия вида (2) по всей толщине пластины (поэтому такое условие отсутствует в математической постановке задачи (1), (3), (5), (6)). В данном случае вполне достаточным является выполнение граничного условия (5), согласно которому для всех £ = 5ї(Ро) температура тела равна начальной температуре. Кроме того, в данной задаче отсутствует также граничное условие вида (4) ввиду того, что оно не влияет на процесс теплообмена в первой его стадии.

Потребуем, чтобы искомое решение задачи (1), (3), (5), (6) удовлетворяло не первоначальному уравнению (1), а осредненному по толщине термического

слоя (интегралу теплового баланса)

Г дй(£ Ро) 1

е (7)

Г1 90(£,Ро) ГЯ1 д

/о №0 ^ “ Л Щ

(1 + /3(0 + То))90(е’Р0)

Вычисляя интеграл в правой части соотношения (7), с учетом (6) находим д Г91 0(£, Ро)с£ = - (1 + /ЗТст) 90(°’ Ро). (8)

дРоУо У У 7

Решение задачи (1), (3), (5), (6) разыскивается в виде следующего полинома

П

©(£, Ро) = £ ьк ы £к, (9)

к=0

где неизвестные коэффициенты Ьк (51) (к = 0,1, 2) находятся из граничных условий (3), (5), (6). После их определения соотношение (9) примет вид

0(С> Ро) = АТ ^1 - . (10)

Подставляя (10) в интеграл теплового баланса (8), будем иметь

А_ Г (л ОУ 2(1+#Гст)

йРо Уо V 91/ 91 (ЬЬ) '

Отсюда для определения неизвестной функции 51 (Ро) приходим к следующему обыкновенному дифференциальному уравнению:

51(Ро)^^1(Ро) = 6 (1 + вТст) йРо, (12)

интегрируя которое при начальном условии 51(0) = 0, находим

91 (Ро) = \/12 (1 + РТст) Ро. (13)

Соотношения (10), (13) определяют решение задачи (1), (3), (5), (6) в первом приближении. Положив 51(Ро1) = 1, из (13) найдем время окончания первой стадии процесса

Ро = Ро1 = ——-----------. (14)

1 12 (1 + /ЗТст) К ;

Введём относительную избыточную температуру по формуле

ё(£, Ро) = (15)

Т ст Т 0

С учётом принятого обозначения соотношение (11) запишется так:

©(£, Ро) = (1-^)2- (16)

При в = 0 (линейная задача) соотношение (14) приводится к виду

= л/12БЪ.

(17)

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

Результаты расчетов нелинейной задачи по формуле (16) при в = 0,001, То = 0 °С, Тст = 100 °С в сравнении с расчётом по конечно-разностному методу прогонки представлены на графиках рис. 2. Их анализ позволяет сделать вывод, что максимальное расхождение результатов в диапазоне числа Фурье 1,5 ■ 10-6 ^ Ро ^ 0,05 не превышает 6%.

Для повышения точности аналитического решения привлечем дополнительные граничные условия. Для этого будем последовательно дифференцировать уравнение (1) по переменной £, а граничные условия (3), (5),

(6) —по переменной Ро.

Для получения первого дополнительного граничного условия продифференцируем (3) по Ро:

д0(0, Ро)

дРо

0.

(18)

Запишем уравнение (1) для точки £ = 0:

90(0, №) _

д0(£, №)1

дЕо

_ д_

(1 + /3(0 + То))

де

(19)

5=0

Рис. 2. Распределение безразмерной температуры 0 = : сплош-

Тст-То

, , ная линия по формуле (27) (второе

Сравнивая (18) и (19), получаем первое приближение); о — метод прогонки; х

дополнительное граничное условие

2

по формуле (16) (первое приближение); в = 0,001

в

д0 (0, Ро) \

де

) + (1 + вТст)

920(О, Ро) д?

0.

(20)

Для получения второго дополнительного граничного условия продифференцируем соотношение (5) по Ро, учитывая, что £ = 51 (Ро) является функцией Ро и, следовательно, 0 = (£, Ро) также будет сложной функцией:

<90 (ді, Ро) с£(?і <90 ((/і, Ро)

<9£ сІРо <9Рэ

Запишем уравнение (1) для Є = ^і(Ро):

0.

(21)

д0(Є, Ро)

дРо

5=91

д_

щ

(1 + в (0 + Т0))

де^о)

9£

(22)

5=91

Сравнивая (21) и (22), с учётом соотношений (5), (6) получим второе дополнительное граничное условие

,23,

Для получения третьего дополнительного граничного условия продифференцируем соотношение (6) по Ро с учётом того, что переменная £ является функцией Ро:

д2@(дъ Ро) <1д 1 д2@(дъ Ро) =

д£2 (ЙЪ д£д¥о ' 1 ;

Продифференцируем уравнение (1) по переменной £ и применим полученное соотношение для £ = 51 (Ро):

д20(£, Ро)

дРод£

_ д_ Г д_

(1+/3(0 + То))9е(^Ро)

. (25)

5=91

Сравнивая (24) и (25), с учётом соотношений (6) и (23) получим третье дополнительное граничное условие

д30(51, Ро) д£3 = 0. (26)

Подставляя (9) (ограничиваясь шестью членами ряда) в граничные условия (3), (5), (6), (20), (23), (26), для определения коэффициентов (51) (к =

0, 1, . . . , 5) получим систему шести алгебраических линейных уравнений. Подставляя найденные из решения этой системы коэффициенты (51) в (9), получаем

£ £2 £3 £4 £5

0(£, Ро) = АГ + вЛ + 52\-Б3^з + Б4^ + В5(27) 51 52 53 54 55

где В\ = 2г — л/Б, Б2 = 2 ^2\[~В — БАТ — 4г^, Б3 = —2 (з\[~В — 10АТ — 6г^,

Б4 = 4\/Б —15АТ—8г, Б5 = 4ДТ+2г-\/Б, Б = 10гАТ+4г2, г = | (1 + /ЗТст).

Подставляя (27) в интеграл теплового баланса (8), приходим к следующему обыкновенному дифференциальному уравнению относительно неизвестной функции 51 (Ро):

5^5 = к1^Ро, (28)

где

= -3о2(1+2^);/2^;_<^-^); к = £ (1 + №,) (5ДГ/3 + 4 + 4№г).

Интегрируя уравнение (27), при начальном условии 51(0) = 0 получим

51 (Ро) = \Z2k\Fo. (29)

Положив в (29) 51 (Ро1) = 1, найдем время окончания первой стадии процесса во втором приближении: Ро1 = (2&1 )-1 или с учетом выше приведенных данных —Ро1 = 0,04713.

Соотношения (27), (29) определяют решение задачи (1), (3), (5), (6) во втором приближении. Непосредственной подстановкой можно убедиться, что все

основные и дополнительные граничные условия, а также интеграл теплового баланса удовлетворяются точно. Дифференциальное уравнение (1) благодаря использованию дополнительных граничных условий точно выполняется во всех точках по координате £, в которых в данный момент времени находится фронт температурного возмущения, включая точки £ = 0 и £ = 1.

Результаты расчетов по формуле (27) и по методу прогонки представлены на рис. 2. Из анализа графиков следует, что в диапазоне чисел 0 ^ Ро ^ ^ 0,0000015 расхождение результатов находится в пределах 1 %. Расчёты по формуле (27) выполнялись для функции 0({, Ро) = 0({, Ро)/ДТ.

Вторая стадия теплового процесса, соответствующая времени Ро ^ Ро1, характеризуется изменением температуры уже по всему сечению пластины, вплоть до наступления стационарного состояния. На этой стадии понятие термического слоя теряет смысл, и в качестве дополнительной искомой функции принимается температура в центре пластины 0(1, Ро) = 52(Ро) (рис. 1).

Математическая постановка задачи в данном случае имеет вид

Начальным условием задачи (30), (31) будет соотношение (10) при Ро = = Ро1 и при 51 (Ро1) = 1, т. е.

Отметим, что специального выполнения начального условия (32) не требуется. Оно будет выполнено в процессе получения решения задачи (30), (31). Это связано с тем, что при Ро = Ро1, т. е. когда 51 (Ро1) = 1 и 52 (Ро1) = 0, математические постановки задач (1), (3), (5), (6) и (30), (31) полностью совпадают. Следовательно, и решения этих задач при Ро = Ро1 будут идентичными. Таким образом, происходит плавный переход от задачи (1), (3), (5), (6) к задаче (30), (31).

Осредняя дифференциальное уравнение (30) по всему объёму тела (0 ^ ^ ^ 1), интеграл теплового баланса с учетом третьего соотношения из (31)

приводится к виду

Как и для первой стадии процесса, решение, удовлетворяющее интегралу теплового баланса (33) и краевым условиям (31), будем искать в виде полинома

к=0

После определения неизвестных коэффициентов С к (52) (к = 0, 1, 2) из граничных условий (31) соотношение (34) примет вид

Э0(£, Ро)

д¥о

0(0, Ро) = ДТ; 0(1, Ро) = (?2(Ро); 90(^РО) = 0.

(31)

0 (Є, Рої) = ДТ (1 - е2).

(32)

(33)

п

0(е, Ро) = £ ок ы ек-

(34)

0(е, Ро) = ДТ - 2 (ДТ - 52) е + (ДТ - 52) е2.

(35)

Подставляя (35) в интеграл теплового баланса (33), получаем уравнение d С1

— J [АТ — 2 (АТ — <?2) С + (АТ — (?2) С2] = 2 (1 + (ЗТст) (АТ — д2) •

Вычисляя интеграл, приходим к следующему обыкновенному дифференциальному уравнению относительно неизвестной функции 92 (Fo):

d (AT + 2g2)

AT — q2

Интегрируя уравнение (36), получим

q2(Fo) = AT + r exp [-3 (1 + вТСт) Fo] , (37)

где r — постоянная интегрирования, определяемая из начального условия. Так как при Fo = Foi температура в центре пластины ({ = 1) равна начальной температуре, то начальным условием для функции 92 будет условие

92 (Foi) = 0.

После определения постоянной интегрирования соотношение (37) примет

вид

92(Fo) = AT {1 - exp [-3 (1 + вТст) (Fo - Foi)]} • (38)

Подставляя (38) в (35), получим решение задачи (30), (31) в первом приближении

0(£, Fo) = AT {1 - {(2 - {) exp [-3 (1 + вТст) (Fo - Foi)]} • (39)

Нахождение решения во втором приближении второй стадии аналогично тому, как это было сделано во втором приближении первой стадии процесса.

Результаты расчетов температур во второй стадии процесса в сравнении с расчетом по методу прогонки и с данными из [1] представлены на рис. 3.

Найдём решение нестационарной нелинейной задачи теплопроводности при степенной зависимости коэффициента температуропроводности от температуры вида а(Т) = аоТв. Математическая постановка задачи для бесконечной пластины при симметричных граничных условиях первого рода в данном случае имеет вид

6 (1 + вТст) dFo.

(36)

дТ(£, Fo) _ д_ <9Fo д(

ЭТ(е, Fo)

Т(£, 0)=То; ^— = 0; T(1,Fo)=Tct.

(40)

(41)

Следуя методу, изложенному выше, решение задачи (40), (41) в первом приближении для первой стадии процесса запишем в виде

0(2/, Fo) = ( 1 — — 9i

где 01^0) = 2 л/ЗРо ехр (/? 1п Тст).

Результаты расчётов безразмерной температуры @(у, Ро) = 0(у, Ро)/АТ по формуле (42) при в = 0,1, Т0 = 0 °С, Тст = 100 С представлены на рис. 4.

Решение задачи (40), (41) во втором приближении запишется следующим образом:

в(у, Fo) = АТ - ^4ДТ + У- + ^6АТ +

2^\ y

Ю 92

2:'V; ■^

V в / q3 V Ю 91 2в q5

где Г] = 8 ((ЗАТ + Тст) - 4\/Тст (5f3АТ + 4Тст), ^(Fo) = 2v/3gS^F°^ % = = —(45Тст + 75/ЗАТ) ехр (/31пТст).

Результаты расчетов безразмерной температуры 0(у, Fo) = 0(у, Fo)/AT по формуле (43) в сравнении с решением в первом приближении даны на рис. 4.

0 0,2 0,4 0,6 0,8 1 -у

1,0

0

0,8

0,6

0,4

0,2

0

1 1 1 J. .

k n\\ V- CN 1 1

V

T \\ IV \V2 \\ \\ \ Fo = 0,03

\ \ \ 4^A°.°2 0,00iv0,0054

0,2 0,4

0,6 У 0,f

Рис. 3. Распределение безразмерной тем- Рис. 4. Распределение температуры в

Т-То . сплошная ли- пластине при степенной зависимости ко-

Тст-То эффициента температуропроводности от

пературы ©

ния —по данным [1]; ◦ — метод прог°нкИ5 X - температуры: 1-по формуле (42); 2-по первое приближение; △ — второе приближение; Формуле (45)

в = 0,001 ф р у ( )

Математическая постановка задачи для второй стадии процесса в данном случае имеет вид

dQ(y, Fo) д

<9Fo ду

(0 + To)

р д@(у, Fo) ду

0(0, Fo) = АТ; 0(1, Fo) = q2(Fo); d0(l,F°) = 0

(Fo ^ Fo1; 0 ^ y ^ 1) ; (44)

(45)

ду

Решение задачи (44), (45) в первом приближении записывается следующим образом:

0(у, Fo)=ДT - у (АТ - д2)(2 - у), (46)

где q2 = AT j 1 - exp — 3TCT (Fo - Foi) |.

2

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

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Кудинов В. А., Карташов Э. М., Калашников В. В. Аналитические решения задач теп-ломассопереноса и термоупругости для многослойных конструкций: Учебн. пособие для вузов. — М.: Высш. шк., 2005. — 430 с.

2. Кудинов В. А. Метод координатных функций в нестационарных задачах теплопроводности. Обзор // Извест. РАН. Энергетика, 2004. — №3. — С. 84-107.

3. Гудмен Т. применение интегральных методов в нелинейных задачах нестационарного теплообмена / В сб.: Проблемы теплообмена. — М.: Атомиздат, 1967. — С. 41-96.

4. Беляев Н. М., Рядно А. А. Методы нестационарной теплопроводности: Учебн. пос. для вузов. — М.: Высш. шк., 1978. — 328 с.

Поступила в редакцию 11/Х/2009; в окончательном варианте — 15/111/2010.

MSC: 80A17, 80M25

RECEPTION OF ANALYTICAL DECISIONS NONLINEAR PROBLEMS OF HEAT CONDUCTIVITY ON THE BASIS OF INTRODUCTION OF ADDITIONAL BOUNDARY CONDITIONS

I. V. Kudinov, V. P. Radchenko

Samara State Technical University,

244, Molodogvardeyskaya str., Samara, 443100.

E-mail: igor_koudinov@mail.ru; radch@samgtu.ru

With the use of an integrated method of thermal balance on the basis of introduction of front of temperature indignation and additional boundary conditions the technique of a finding of analytical decisions of nonlinear problems of the n,on,-sta,tion,a,ry heat conductivity is studied, allowing to receive satisfactory accuracy of the decision in all range of change of Fourier number. Decisions have a simple appearance of the sedate algebraic polynoms which do not contain special functions.

Key words: nonlinear problems, analytical methods, front of temperature indignation, additional boundary conditions, integral of thermal balance.

Original article submitted 11/X/2009; revision submitted 15/III/2010.

Igor V. Kudinov, Postgraduate Student, Dept. of Applied Mathematics & Computer Science. Vladimir P. Radchenko (Dr. Sci. (Phys. & Math.)), Head of Dept., Dept. of Applied Mathematics & Computer Science.