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

Е.В. Стефанюк, В.П. Радченко ТЕПЛОПРОВОДНОСТЬ В ПЛАСТИНЕ

ПРИ ПЕРЕМЕННЫХ ВО ВРЕМЕНИ ГРАНИЧНЫХ УСЛОВИЯХ ТРЕТЬЕГО РОДА. ТЕМПЕРАТУРА СРЕДЫ - ЭКСПОНЕНЦИАЛЬНАЯ ФУНКЦИЯ ВРЕМЕНИ

Путем совместного использования методов Л. В. Канторовича и ортогонального метода Бубнова-Галеркина получено приближенное аналитическое решение задачи теплопроводности для бесконечно протяженной пластины при переменных во времени граничных условиях 3—го рода.

Математическая постановка задачи при наличии переменных во времени граничных условий существенно усложняется. Применение точных аналитических методов для решения таких задач оказывается затруднительным. Для их решения весьма полезным может быть применение комбинации точных (Фурье, интегральных преобразований и др.) и приближенных аналитических (вариационных, взвешенных невязок, коллокаций и др.) методов [1-5].

Рассмотрим совместное использование метода Л.В. Канторовича и ортогонального метода Бубнова-Галеркина применительно к решению задачи теплопроводности для бесконечнопротяженной пластины при экспоненциальном изменении температуры среды. Математическая постановка задачи имеет вид

(Т > 0; 0£ Ц < «); (1)

дт д Г1

т ад = Т0; (2)

дТ (0,т)

дц

= 0; (3)

Я д7|Цт) + «[т (8 ,т) - Тср (т)]= 0, (4)

где Т - температура; Т0 - начальная температура; Тср (т) = Тм - (Тм - Т0)е~Ьт - температура среды; ц - координата; т - время; 8 - толщина пластины; а - коэффициент теплоотдачи; Я - коэффициент теплопроводности; а - коэффициент температуропроводности; Ь - коэффициент; Тм - максимальная температура внешней среды.

Введем следующие безразмерные переменные и параметры:

0 = (Т (х, Го) - Т0)/(Тм - Т0); Ві = а8/Я; х = ц/8 ; Го = ат / 82, Рё = Ь82/а, (5) где 0 - относительная избыточная температура; х - безразмерная координата; Го - число Фурье; Ві - число Био; Рё - критерий Предводителева.

С учетом принятых обозначений задача (1) - (4) примет вид

д0( х,Гв > = д 2 0( х2 Гв) (Го > 0; 0 £ х < 1); (6)

дГо дх2

0(х,0) = 0; (7)

д0(0,Гв > = 0; (8)

дх

д0(1, Го )

- Ві[і - ехр(-РёГо) - 0(1, Го)] = 0 . (9)

дх

Следуя методу Л.В. Канторовича, решение задачи (6) - (9) в первом приближении [4] принимается в виде

0(х, Го) = 1 — exp(—PdГo)/(Го)р(х), (10)

где / (Го)- неизвестная функция времени; р (х) = Б + 2 — х2 - координатная функция.

Б1

Благодаря принятой конструкции координатной функции р(х) соотношение (10) точно удовлетворяет граничным условиям (8), (9). Для нахождения неизвестной функции времени / (Го) в нулевом приближении составляется интеграл взвешенной невязки дифференциального уравнения (6), т.е.

Рё ехр(-РёГо) -”^ "у р(х) + /(Го)д Р2х)

дх

д/ (Го) дГо

ёх = 0.

(11)

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

Рё + //(Го)

ґВі + 2 1Л

Ві

3

+ 2 /1( Го) = 0 / '(Го) + Б/ (Го) - 0.5БРё сс^( РёГо) = 0, (12)

5 (Б/ )2 +15 Б/

где Б =-------2------------коэффициент.

2 (Б/ )2 + 10Б/ +15

Для решения уравнения (12) воспользуемся методом варьирования произвольной постоянной [1]. Для этого найдем сначала решение однородного уравнения

/ '(Го) + Б/ (Го) = 0, (13)

которое имеет вид

/(Го) = С ехр(—БГо). (14)

Решение для уравнения (12) будем разыскивать в виде

/(Го) = С (Го) ехр(—БГо), (15)

считая С не постоянной, а функцией времени. Подставляя решение (15) в уравнение (12), получим

С'(Го) = 0,5 Б Pd оо8(PdГo) ехр(БГо). (16)

Определяя интеграл в (16), получим

Б Pd ехр(БГо) |

С (Го) =

[Б ос8(РёГо) + Рё 8Іп(РёГо)] + С1.

2(Б2 + Рё2)

где С1 - постоянная интегрирования. Подставляя (17), (15) в (10), получим

Б Рё ехр(БГо)

(17)

0( х, Го) = 8Іп( РёГо) -

2(Б2 + Рё2)

[Б ос8( РёГо) + Рё 8Іп( РёГо)] + ]

ехр(-БГо)р(х). (18)

Для определения постоянной интегрирования С^ составляется интеграл взвешенной невязки начального условия, т.е.

| ^т(Рё • 0) - /(0)р(х))р(х)ёх = 0.

(19)

Б2 Рё

Отсюда С, =-------------—. С учетом найденного значения С, соотношение (18) в первом

1 2(Б2 + Pd2) 1

приближении принимает вид 0( х, Го) = 8т( PdГo) —

Б Рё ехр(БГо) [Б ссъ(РёГо) + Рё біп(РёГо)] - - ° 2

2( Б2 +Рё2)

2(Б2 Р Рё2)

Xехр(-БГо)р(х). (20)

Соотношение (20) полностью совпадает с решением в первом приближении задачи (6)—(9), полученным в работе [1] путем совместного использования интегральных преобразований Лапласа и ортогонального метода Бубнова - Галеркина.

Решение задачи (6) - (9) при любом числе приближений находится в виде

0( х Го) = вІп( РёГо) - £ Тк (Го )р к (х),

(21)

к=1

Б/ + 2к 2к 1 т-г

где рк (х) =-х - координатные функции. При использовании таких координатных

Б/

функций соотношение (21) в любом приближении удовлетворяет граничным условиям (8), (9). Найдем решение задачи во втором приближении. Для нахождения /к (Го) (к = 1,2) из условия ортогональности невязки уравнения (6) к координатным функциям р1 (х) и р 2 (х) получим

Рё - £(/; (то рк (х) - /к (Го р (х))

к=1

р ■ (х)ёх = 0 (у = 1,2).

(22)

Определяя интегралы, относительно неизвестных функций времени приходим к следующей системе двух неоднородных дифференциальных уравнений первого порядка

[ Ап/1+ Au /' + Бп/ + Б12 / + Ъ = 0;

|_ А21/ 7 + А22 /2 + Б2\/\ + Б22 /2 + Ъ2 = 0

1 1 1

где Ак = | Рк (x)PJ(х№; Бк = —| Рк (х)р3 (х)dx; Ък = Pd | Рк(х) ^.(У = к =1,2).

0 0 0

Решение системы уравнений (23) разыскивается в виде

/:(Го) = / Д Го) + /!* (Го); (24)

/2(Го) = /2(Го) + /2* (Го), (25)

где /1(Го), /2(Го)- частные решения неоднородной системы уравнений (23); /1*(Го),

/2* (Го) - общие решения соответствующей однородной системы, т.е. при Ъ1 = Ъ2 = 0 .

Общие решения соответствующей однородной системы находятся в виде

// (Го) = Де ^; /2* (Го) = Б2* ^, (26)

где Б1, Б2, т - некоторые пока неизвестные постоянные. Подставляя (26) в (23) и положив

= Ъ2 = 0, получим следующую систему двух однородных алгебраических линейных урав-

нений:

Б1 (А11т + Б11) + Б2(А12 т + Б12) = 0; (2уч

Б1 (А21т+б21)+б2 ( а22 т+б22) = 0.

Однородная система имеет нетривиальное решение в случае, если ее определитель равен нулю, т.е.

А11т + Б11 А12т + Б12 = 0 А2т + Б21 А22 т + Б22

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

2-ой степени. Так как элементы матрицы положительны и симметричны относительно главной диагонали, то корни тк (к = 1,2) должны быть действительными отрицательными числами. Эти корни являются приближенными значениями собственных чисел краевой задачи Штурма -Лиувилля для уравнения теплопроводности с переменными во времени граничными условиями

3-го рода [2].

Подставляя тк (к = 1,2) в (25), для каждого собственного числа находятся неизвестные коэффициенты Бку (к = у = 1,2). При этом в однородной системе следует положить Бк1 = 1

(к = 1,2). Неизвестные функции времени /к * (Го) находятся из (26).

С учетом найденных значений Бу (к = у = 1,2) соотношения (25) примут вид

А*(Го) = Би ехр(т1Го)+Б21 ехр(т 2Го); (28)

/2* (Го) = Б12 ехр(т 1Го) + Б22 ехр(т2Го). (29)

Умножая в (28), (29) слагаемые, соответствующие корню т1, на произвольную постоянную С1 и корню т 2 - на С 2, а также учитывая, что Бп = Б21 = 1, получим

./1*(Го) = С1 ехр(т 1Го) + С 2 ехр(т 2Го); (30)

/2* (Го) = С1Б12 ехр(т1 Го) + С2Б22 ехр(т2Го), (31)

где С1 и С2 находятся из начального условия (7) (см. ниже формулу (43)). Соотношения (30), (31) представляют общее решение системы уравнений (23) при Ъ1 = Ъ2 = 0 .

Частные решения системы уравнений (23) разыскиваются в виде (30), (31), считая С1 и С2 не постоянными, а функциями времени, т.е.

/1 (Го) = С1 (Го) ехр(т 1 Го) + С2 (Го) ехр(т2Го); (32)

23

/2(Г°) = С1 (Р°)°12 ЄХР(^ 1 ^°) + С2(Р°)°22 ЄХР(т2¥°) •

Подставляя (32), (33) в (23), найдем

' / /

Сі у (-^ц + А12^12) + С2 V2 (Ац + А12^22) + -^1 = 0 ;

/ /

С1 "У1 (-^21 + А22^12) + С2 У 2 (А21 + А22^22) + ^2 = 0 •

Решая систему уравнений (34), (35), получим

С, = _ 7 (33 + 26^22) 1 624 у1(Д2 - Б22)

С =

2

7 (33 + 26Д2) 624 у2(Д2 - Е22)

Интегрируя (36), найдем

7 (33 + 26Е22)

+X;С2^°) = -

7 (33 + 26Е12)

(33)

(34)

(35)

(36)

(37)

624 г1т1(Б12 — Б22) 624 у2т2(Б12 — Б22)

где Х2 - постоянные интегрирования (так как находятся частные решения, то можно при-

нять х = Х2 = 0).

Подставляя (37) в (32), (33), получим

/ 1(Р°) =

7 (33 + 26Е22) 7 (33 + 26Е12) ;

624 ^(А2 -Б22) 624 т2(Д2 -Б22)

4-7^12 (33 + 26^22) 7^22 (33 + 26^)

/ 2 (^ °) •

2 624 т1(^12 -е22) 624 т2(А2 -е22)

Подставляя соотношения (30), (31), (38), (39) в (24), (25), найдем

7 (33 + 26Е22) 7 (33 + 26Е12)

(38)

(39)

/1 (г°) = с1 ехр(т р° ) + с 2 ехр(т 2 р° ) +

624 т1(Р12 -Е22) 624 т2ф12 -е22)

/2 (г°) = д2с ехр(т р° )+е22с 2 ехр(т 2 р° ) +

Подставляя (40), (41) в (21), получим

7Д2 (33 + 26Е22) 7Б22 (33 + 26Д2)

624 ^(А2 -А2) 624 т2(Д2 -Б22)

(40)

• (41)

0(х, ¥°) = Рс№° +

с1 ехр(т р° ) + с 2 ехр(т 2 р° ) +

7 (33 + 26Е22) 7 (33 + 26Е12)

Е12С1 ехр(т 1^°) + Е22С2 ехр(т2¥°) +

624 т 1 (А2 - Е22) 624 т2(а2 - а2 )

7£>12 (33 + 26^22) 7^22 (33 + 26^12)

624 т 1СА2 Е22) 624 т 2СА2 Е22)

91( X) + (42)

Для определения постоянных С1 и С2 составляется невязка начального условия (7) и требуется ортогональность невязки к координатным функциям <р1 (х) и <р2 (х), т.е.

10(х,0) 9 ■ (x)dx = 0 (у = 1,2).

(43)

Из решения системы двух алгебраических линейных уравнений (43) находятся постоянные С и С2. После их определения решение задачи (6) - (9) в замкнутом виде находится из (42).

Результаты расчетов относительной избыточной температуры по формуле (22) в сравнении с точным решением [2] представлены на рис.1. В точном решении было взято 30 членов ряда (формула (9), с. 289, [2]).

Собственные числа для двух и четырех приближений при Ві = 1, а также точные их значения представлены в таблице.

Анализируя изменение невязки Є уравнения (6) по координате х при Г° = 0,01 (рис. 2), можно заметить, что максимальное значение она принимает в точке х =1 . Анализ изменения невязки уравнения (6) во времени в точках х = 0 и х = 1 (рис. 3 и 4) показывает, что при ¥° > 0,01 она становится практически равной нулю. Невязка начального условия во всем диапазоне координаты х равна нулю. Отметим, что нулевая невязка начального условия имеет место уже в первом приближении (см. формулу (20)).

+

Я Число приближений Точные значения [2]

2 4

к 0,86033361831 0,860333589028290 0,86033358901938148

Я2 3,464594611 3,42561892126255 3,425618459481728

Я3 6,44986878307776 6,4372981791719471

І4 10,71473561814056 9,5293344053619636

©(Цо)

РА

0.008 0.006 0.004 0.002

/

/

7

О 0.01 0.02 0.03 0.04

^0

Р и с. 1. Графики распределения относительной избыточной температуры от времени ( Б1 = 1, Го = 0,1, х = 1):-------------расчет по формуле (21) (четвертое приближение);-------точное решение [2]

Р и с. 2. Изменение невязки £ уравнения (6) от координаты х при п = 4 (четыре при -ближения): Го = 0,01, Б1 = 1

Рис. 3. Изменение невязки £ уравнения (6) во времени для п = 4 в точке х = 0 при Б1 = 1

Р и с. 4. Изменение невязки £ уравнения (6) во времени для п = 4 в точке х = 1 при Б1 = 1

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

1. ЦойП.В. Методы решения отдельных задач тепломассопереноса. М.: Энергия. 1971. 384 с.

2. ЛыковА.В. Теория теплопроводности. М.: Высшая школа, 1967. 600 с.

3. Кудинов В.А., Карташов Э.М. и др. Тепломассоперенос и термоупругость в многослойных конструкциях. М.: Энергоатомиздат, 1997, 425 с.

4. Кудинов В.А. Аналитические методы решения краевых задач для многослойных конструкций (Обзор) Изв. АН. Энергетика. №5, 1999. С.86-107.

5. Карташов Э.М. Аналитические методы в теории теплопроводности твердых тел. М.: Высшая школа. 2001. 550 с.

Поступила 2.02.2004 г.