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

Л И Т Е Р А Т У Р А

1. К е р н, Д. Развитые поверхности теплообмена: пер. с англ. / Д. Керн, А. Краус. -М.: Энергия, 1977. - 464 с.

2. О с н о в ы расчета и проектирования теплообменников воздушного охлаждения: справ. / под общ. ред. В. Б. Кунтыша, А. Н. Бессонного. - СПб.: Недра, 1996. - 512 с.

3. Т е п л о о т д а ч а и аэродинамическое сопротивление шахматных пучков из круглых труб с подогнутыми спиральными КЬМ-ребрами / В. Б. Кунтыш [и др.] // Химическое и нефтегазовое машиностроение. - 2003. - № 11. - С. 10-14.

4. К у н т ы ш, В. Б. Анализ тепловой, объемной и массовой характеристик теплооб-менных секций аппаратов воздушного охлаждения / В. Б. Кунтыш, А. Э. Пиир // Химическое и нефтегазовое машиностроение. - 2009. - № 5. - С. 3-6.

5. К у н т ы ш, В. Б. Основные способы энергетического совершенствования аппаратов воздушного охлаждения / В. Б. Кунтыш, А. Н. Бессонный, А. А. Брилль // Химическое и нефтегазовое машиностроение. - 1997. - № 4. - С. 41-44.

6. Т е п л о о б м е н н а я секция: пат. 2213920 России, МПК С2, Е 28 Б 3/02 / В. П. Му-лин, И. И. Кочетов, Р. Ф. Теляев, В. Б. Кунтыш, В. И. Мелехов, А. В. Самородов; заявитель ЗАО «Октябрьскхиммаш». - № 2001119695; заявл. 16.07.2001; опубл. 15.07.2003 // Бюл. изобрет. / Роспатент. - 2003. - № 28. - С. 70.

7. У с т р о й с т в о для изготовления теплообменной трубы со спирально-навивными ребрами: заявка на полезную модель Респ. Беларусь, МПК В 21 Б 11/06 / В. Б. Кунтыш, В. П. Мулин, Е. С. Санкович, А. Ш. Миннигалеев; заявитель Белорус. гос. технол. ун-т. -№ и 20120381; заявл. 05.04.12.

8. Т е п л о о б м е н н а я ребристая труба: пат. 4814 Респ. Беларусь, МПК Е 28 Е 1/00 / В. Б. Кунтыш, В. И. Володин, Е. С. Санкович, В. П. Мулин, А. Э. Пиир, А. Ш. Миннигалеев, Г. Г. Баранов; заявитель Белорус. гос. технол. ун-т. - № и 20080322; заявл. 17.04.08; опубл. 30.10.08 // Афщыйны бюл. / Нац. цэнтр штэлектуал. уласнасщ. - 2008. - № 5. - С. 36.

9. Т е п л о о б м е н н а я биметаллическая ребристая труба: пат. 14907 Респ. Беларусь, МПК Е 28 Е 1/00 / В. Б. Кунтыш, Е. С. Санкович, В. П. Мулин, А. Ш. Миннигалеев, А. Э. Пиир, И. Р. Гаязов, А. Л. Соловьев; заявитель Белорус. гос. технол. ун-т. -№ и 20091539; заявл. 28.10.09; опубл. 30.10.11 // Афщыйны бюл. / Нац. цэнтр штэлектуал. уласнасщ. - 2011. - № 4. - С. 58.

10. С п о с о б и устройство для изготовления теплообменной трубы с КЬМ-ребрами: пат. 16177 Респ. Беларусь, МПК В 21 С 37/15 / В. Б. Кунтыш, В. П. Мулин, Е. С. Санкович, А. Э. Пиир, А. Ш. Миннигалеев, А. Л. Соловьев, О. В. Петрович; заявитель Белорус. гос. технол. ун-т. - № а 2010366; заявл. 11.03.10; опубл. 30.08.12 // Афщыйны бюл. / Нац. цэнтр штэлектуал. уласнасщ. - 2012. - № 4. - С. 44.

Представлена кафедрой энергосбережения,

гидравлики и теплотехники Поступила 12.11.2012

УДК 536.2 (075)

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

Докт. физ.-мат. наук КУДИНОВ В. А., инженеры КОТОВА Е. В., ЕРЕМИН А. В., КУЗНЕЦОВА А. Э.

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

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

44

выполнения условий сопряжения в виде равенства температур и тепловых потоков в точках контакта слоев. Выполнение этих условий в случае применения точных аналитических методов приводит к необходимости решения характеристической системы в виде цепочного трансцендентного уравнения относительно собственных чисел краевой задачи, решение которого может быть получено лишь численными методами [1, 2].

Из приближенных аналитических методов решения достаточно широкое распространение получили ортогональные методы взвешенных невязок (Л. В. Канторовича, Бубнова - Галеркина и др.) [2-4]. Однако их применение к расчетам многослойных конструкций существенно сдерживается необходимостью построения систем координатных функций, удовлетворяющих граничным условиям и условиям сопряжения. В настоящей работе рассматривается метод, позволяющий строить системы координатных функций, в любом приближении точно удовлетворяющих граничным условиям и условиям сопряжения, независимо от числа контактирующих тел. Построение таких систем оказалось возможным благодаря принятию глобальной системы неизвестных функций времени (одинаковой для всех контактирующих тел), введение которой позволяет приводить многослойную конструкцию к однослойной с переменными (кусочно-однородными) свойствами среды.

Для получения решений при малых и сверхмалых значениях временной переменной необходимо использовать большое число приближений, что в ортогональных методах приводит к необходимости решения алгебраических уравнений высоких степеней (относительно собственных чисел краевой задачи), точность решения которых даже при использовании современной компьютерной техники не всегда может удовлетворять потребностям практики. Наиболее эффективным в области малых значений временной переменной является метод, основанный на определении фронта температурного возмущения и дополнительных граничных условий [3, 6]. В настоящей статье он применен для расчета температурного состояния последнего слоя многослойной системы, так как при малых и особенно сверхмалых значениях времени теплообмен протекает лишь в этом слое.

В качестве конкретного примера найдем решение задачи теплопроводности для трехслойной пластины в следующей математической постановке (рис. 1):

8Т (х, т) _ а 8%(х,т)

— ? 2 8т 8х

(т > 0; х1._1 < х < х(; 1 _ 1,2,3; х0 _ 0; х3 _ 5);

(1)

Т (х, 0) _ Т0; 8Т1(0, т) __

(2)

8х

(3)

Т3(5,т) _ Г;

(4)

Ti (х, т) _ Т+1(х,т) (1 _ 1,2);

(5)

КдТ г (х 1 ,т) +1дТ<+1 (X ,т)

дх

дх

(1 = 1,2),

(6)

где Тг - температура 1-го слоя; х - координата; т - время; 81, 62, 83 - толщины слоев; 8 - суммарная толщина трехслойной пластины; аг - коэффициент температуропроводности 1-го слоя; То - начальная температура; Тст - температура стенки при х = 8; Х - коэффициент теплопроводности 1-го слоя.

Т Х а1

81

Тг

Х а

х1 х2

х

Тз

Хз

а

Рис. 1. Расчетная схема трехслойной пластины Введем следующие безразмерные переменные:

© = Т - ТСт с = х р = ат

М1 Т - Т ' ? 1 5; 0 52,

(7)

где © - относительная избыточная температура; Е - безразмерная координата; Бо - число Фурье; а - наименьший из коэффициентов температуропроводности а1 (1 = 1, 2, 3).

Задача (1)-(6) с учетом (7) приводится к виду:

д©, (^о) = а д2© , (^о) ЗБо а д£,2

(Бо > 0; < ^ < 1 = 1,2,3; ^ = 0; ^ = 1);

©1 &0) = 1; д©1 (0, Бо)

д^

= 0;

© з(1,Ро) = 0; ©,(^о) = ©1+1 (^о) (1 = 1,2);

(8)

(9) (10)

(11) (12)

Т

х

З®,(4,, Fo) = X,+1ЗЭ,+1 (4,, Fo) (, = ^ (13)

84 34

Решение задачи (8)—(13), следуя ортогональному методу Л. В. Канторовича, принимается в виде

®,(4, Fo) = ¿Л^с) фк,(4) (, = 1, 2, 3), (14)

к=1

где /к (Бо) - неизвестные функции; фк.(4) (, = 1,2,3) - координатные функции, определяемые таким образом, чтобы заранее точно выполнялись граничные условия (10), (11) и условия сопряжения (12), (13).

Отличительной особенностью используемого здесь метода является принятие одинаковой системы неизвестных функций времени /к (Бо) для каждого из контактирующих тел, что позволяет получать наиболее простого вида координатные функции фш (4) (, = 1, 2,3; к = 1, п).

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

Фкз(4) = 1 -42к (к = М). (15)

Координатная функция для второго слоя принимается в виде

Фк2(4) = Л + 442к, (16)

где неизвестные коэффициенты Л\ и Л 2 находятся из условий сопряжения (12), (13) между вторым и третьим слоями. После их определения соотношение (16) принимает вид

' О

X

Фк2(4) = 1 - 1 -т^ 42к-т^Г (к = 1,п). (17)

Ч Х2 ) Х2

Координатная функция для первого слоя принимается в виде

Фк1(4) = В + В£2к, (18)

где В\, В2 - неизвестные коэффициенты.

Очевидно, что соотношение (14) при использовании в качестве координатной функции (18) точно удовлетворяет граничному условию (10). Неизвестные коэффициенты В1, В2 находятся из условий сопряжения (12), (13) между первым и вторым слоями. После их определения соотношение (18) будет

^ Х3 Х3 ^

2 Х1 )

Фк1(4) = 1 - ^4?к - 1 42к -^Ч2к (к = 1,п). (19)

' хЛ

Х 2 )

X,

Непосредственной подстановкой можно убедиться, что при использовании в качестве координатных функций соотношений (15), (17), (19) со-

отношение (14) в любом приближении (к = 1, п) точно удовлетворяет граничным условиям (10), (11) и условиям сопряжения (12), (13) при любых значениях неизвестных функций /к (Ро).

Найдем решение задачи (8)—(13) в нулевом приближении. Для этого потребуем, чтобы соотношение (14), в котором ограничимся одним членом ряда, удовлетворяло не уравнению (8), а некоторым осредненным уравнениям (проинтегрированным в пределах толщины каждого слоя)

3 Ч

XI

с&,(Ч,Ро) а, д=©,(Ч,Ро)

5Ро

дЧ2

, = 0.

(20)

Подставляя (14) в (20), находим

"Ч Ч2

/1'(Ро)

| (в - вч2)йч +1 В - +1 (1 - Ч2)йЧ

Ч1

+2 /1(Ро)

В М й Ч+в2 а2( й Ч+аз I й Ч

а * а * а

(21)

= 0,

где / - первая производная от функции /1(Ро) по переменной Ро.

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

й/1(Ро)

йРо

" У/1 (Ро) = 0,

(22)

где

V =

В

В.

1

^ =Ч11 вЧ 1+в2{Ч-Ч) -^-(Ч-Ч3) + (1 -Ч2)--(1-Ч2)+т;

^2 =2 (а1 ВЧ + а2В3 (Ч - Ч) + аз (1 - Ч2)) /а;

В = 1 -

Хз Хз КХ 2 Х1 у

Ч2к -

1 -Хз

Х

Ч2к .

2

2

(

В = В2 = 1 -Х,

1 Хз

Х

Л

Ч2к; Вз

2 У

Интегрируя уравнение (22), находим

/1(Ро) = Cexp(vРo),

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

© , (Ч,Ро) = Сехр^Ро^ДЧ) (, = 1,2, з).

(2з)

(24)

, -1

2

Для определения постоянной интегрирования С составим невязку начального условия (9) и проинтегрируем ее в пределах толщины каждого слоя

3 %

£{[©, -1] ^ = 0. (25)

¿=11

Подставляя (24) в (25), находим

^ 1

| [С(D - 2) -1]d^ +1 [С(D2 - Dз^2) -1]+| [С(1 - ^2) -1]d£ = 0.

0 ^ ^

Определяя интегралы, относительно постоянной интегрирования С получаем алгебраическое линейное уравнение. Его решение С = 1/ць

После определения постоянной интегрирования решение задачи (8)-(13) в нулевом приближении принимает вид

©■ &Р°) = еХР(УР°)ф1-® (■ = 1,2,3). (26)

Если положить Х1 = X2 = Х3 = X; а1 = а2 = а3 = а (однослойная пластина), то соотношение (26) будет

©(£, Б°) = 1,5ехр(-3Б°)(1 -^2). (27)

Последняя формула совпадает с решением в нулевом приближении для однослойной пластины, приведенным в [2, 3]. Результаты расчетов по формуле (27) представлены на рис. 2.

Используя формулу (26), найдем решение конкретной задачи для трехслойной пластины при следующих исходных данных:

X = 0,00086 м; х2 = 0,00261 м; х3 = 0,00502 м; X = 11 Вт/(мК);

X2 = 2 Вт/(мК); Х3 = 1,1 Вт/(мК); а1 = 3,624 • 10-6 м2/с;

а2 = 1,02 • 10-6 м2/с; а3 = 0,94 • 10-6 м2/с. (28)

Анализ результатов расчетов по формуле (26) в сравнении расчетом численным методом прогонки позволяет заключить, что в диапазоне чисел Фурье 0,187 < < да расхождение составляет около 12 % (рис. 3).

Для повышения точности найдем решение задачи (8)—(13) в первом приближении. Для этого составим невязки уравнений (8) и потребуем ортогональности невязок к координатным функциям первого приближения Фи (■ = 1, 2, 3) (учитываем, что согласно (27) наименьшим является коэффициент температуропроводности аъ)

3 %

ЕI

¿=1 £ £ -1

^ <М9 - №) а

а3 ос,

Ф1ДМ £ = 0.

Рис. 2. Распределение температуры в трехслойной пластине, приведенной к однослойной: А - по формуле (26) (нулевое приближение); х - по формуле (33) (первое приближение); --по формуле (14) (восьмое приближение); О - точное решение [5]

Рис. 3. Распределение температуры в трехслойной пластине: х - по формуле (26) (нулевое приближение); А - по формуле (31) (первое приближение); --по формуле (14) (восьмое приближение); О - метод прогонки; □ - точное решение [5]

Последнее соотношение можно переписать в виде

"Е Е2 1

/те

| (В - ДЕ2)ф^Е +1 (- ф^Е + | (1 - Е2)Ф1з^Е

0

+ 2 /(Бо)

Г „ Е1 „ Е2 1 ^

а г а г г

А~ I Фп^Е+В2— I Ф12^Е + I Ф13^Е

а ^ а ^ ^

У "з о "з Е1 Е2

= 0.

Определяя интегралы в последнем соотношении, относительно неизвестной функции Л (Бо) получаем обыкновенное дифференциальное уравнение вида

^(Бо)/ dFo + v1 ¿(Бо) = 0, (29)

где v; =ц4/ц3; ц3 = 0,516004; =1,042163. Интегрируя уравнение (29), находим

Л1(Бо) = С^хр^Бо), (30)

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

© ,(Е,Бо) = С ехр^Бо^ДЕ) (/■ = 1,2,3). (31)

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

I [С1 (D - DlЕ2)-1]ф1^Е +I [С ф2 - DзЕ2) -1]Ф12

0 Е1

+ | [С1(1 -Е2) -1]ф1э^ Е = 0.

(32)

Определяя интегралы в (32), относительно постоянной интегрирования С1 получаем алгебраическое линейное уравнение. Его решение С1 = 1/ ц,3.

После определения постоянной интегрирования С1 решение задачи (8)-(13) в первом приближении находится из (31). Если положить X = X2 = Х3 = X; а1 = а2 = а3 = а, то соотношение (31) приводится к виду

©(Е,Бо) = 1,25 ехр(-2,5 Бо)(1 -Е2). (33)

Формула (33) совпадает с решениями в первом приближении для однослойной пластины, приведенными в [2-4].

Результаты расчетов по формулам (26), (27), (31), (33) даны на рис. 2, 3. Их анализ позволяет заключить, что точность решения в первом приближении по сравнению с нулевым существенно повышается. И в частности, расхождение с точным решением (рис. 2) не превышает 3 %.

Для дальнейшего повышения точности найдем решение задачи (8)-(13) во втором приближении. Составляя невязку уравнений (8) и требуя орто-

гональности невязки к координатным функциям фь. (4) и ф2.. (4) (. = 1,2,3), находим

3 %

е ц

.=1 £

ЛФ + /2ф 2.--^ ( /"Ф". + /2 Ф2'. )

ф „

= 0 (] = 1,2), (34)

где ф", ф'. - вторые производные от функций фь. (4) и ф2. (4) по переменной 4.

Соотношение (34) можно переписать в виде:

4 ъ. "

| (Л"ф1" + У" )ф"1^4 + | (Л"'ф12 + /2'ф22)ф"244 + | (Л"'ф13 + /2'ф23)ф13^4 -

0 4" 4'

4" а 4' 1

11 (Л"ф"1 + ЛЮф^ 4 —-1 (/"Ф"' + /2Ф2'2)Ф12 4 4 -1 ( Л"ФТЭ + /2ф2э)ф13^ 4 = 0;

а3 4" 4'

4" 4' 1

| (/"ф11 + /2'ф21)ф214 4 + | (Л\ф\2 + /2'ф22)ф2244+ | (/"Ф13 + Л'ф'эМэ4 40 4\ 4'

4' 1

11 (/\фп + Лф^ф^-—I (/"ф"' + /2ф2'2)ф22 4 4 - / (/\ф"з + Лф*)^ 4 = 0.

"3 0

•♦3 0

-3 4"

Определяя интегралы в последних соотношениях, относительно неизвестных функций /1 (Ео) и /2 (Ео) с учетом исходных данных (27) получаем следующую систему двух обыкновенных дифференциальных уравнений:

Ь\/1'+ Ъ2 /''+ ¿3/ + ¿4 /' = 0;]

/\/1'+ /'/2'+/3/1 + /4/2 = 0, ]

(35)

где Ъ\ = 0,494399; ¿2 = 0,585074; ¿3 = 0,958027; ¿4 = 1,49522; /\ = 0,585074; /2 = 0,702963; /3 = 1,22521; /4 = 2,17489.

Частные решения системы уравнений (35) разыскиваются в виде:

/"(Ео) = А ехр(гЕо); /'(Ео) = В ехр(гЕо),

где А, В, г - некоторые постоянные. Подставляя (36) в (35), находим:

А(0,494399 г + 0,958027) + В (0,585074 г +1,49522) = 0;) А(0,585074 г +1,22521) + В(0,702963 г + 2,17489) = 0.

(36)

(37)

Система однородных алгебраических уравнений (37) имеет нетривиальное решение в случае, если определитель ее равен нулю:

0,494399 г + 0,958027 0,585074 г +1,4952 0,585074 г +1,22521 0,702963 г + 2,17489

= 0.

(38)

Из (38) получаем следующее характеристическое уравнение: 0,005232г2 + 0,157067г + 0,25164 = 0. Его корни: гх =-1,698214; г2 =-28,320207.

Подставляя величину корня т\ в систему уравнений (37), получаем:

0,113714 A\ + 0,475086 Д = 0;| 0,2303\6 A + 0,962233 Д = 0.J

Для однородной системы (39) можно положить A\ = const = 1. Тогда D\ = -0,239355.

Аналогично, подставляя Г2 в (37), находим: A2 = 1; D2 = -0,881653. С учетом найденных значений постоянных Ai, D, r (i = 1, 2, 3) соотношения (36) примут вид:

/\(Fo) = A^0 + A^; /2 (Fo) = DxerFo + D2erFo.

r\Fo

(40)

Чтобы найти общее решение системы уравнений (35), умножим частное решение, включающее корень п, на произвольную постоянную С1, а решение, включающее корень тг, - на постоянную С2. Подставляя полученные общие решения в (14) (при п = 2), находим

©,№) = (С;А^ + С2А^К + (СДе^ + С2ДОф2, (, = 1, 2,3). (41)

Для определения постоянных С1 и С2 составляется невязка начального условия (9) и требуется ортогональность невязки к координатным функциям Ф1Д) и Ф2 ,(%) 0' = 1,2,3)

2 %

[(СД +С2А2)ф,+(С1Б1 +С2Б2)ф2,-1]ф^% = 0 (у = 1,2,3). (42)

Соотношение (42) относительно неизвестных С1 и С2 приводится к следующей системе двух алгебраических линейных уравнений:

f (Е\Ф\\ + £>21 - 1)Ф\\^% + f (Е\Ф\2 + E2^22 - 1)Ф\2d% -

0 %1 \

+ f (Е\Ф\3 + Е2Ф23 - 1)Ф\3d% = 0;

%2

1 %2

f (Е\Ф\\ + Е2Ф2\ - 1)Ф2\d% + f (Е\Ф\2 + Е2Ф22 1) Ф 22

0 %1 \

+ f (Е\Ф\3 + Е2Ф23 -1)Ф2зd% = 0,

%2

где Е = CA\ + C2A2; Е2 = ClDl + C2D2. Определяя интегралы, находим:

0,356267C\ + 0,2\4339C2 + 0,645874 = 0; | 0,4\91102C\ + 0,346957C2 + 0,795890 = 0.J

Из решения системы уравнений (43) получаем: С = 1,642377; С2 = -3,631466.

После определения постоянных С и С2 решение задачи (8)-(13) во втором приближении находится из (41). Если положить Х1 =Х2 = Х3 = X; а1 = а2 = а3 = а (однослойная пластина), то соотношение (41) приводится к виду

©( Е,Бо) = (1,55 в-2А17° +3,3 е-25'52¥о)(1 -42) -

(44)

- (0,28е~2'47ро +2,9048Л2е"25'52ро)(1 -^4).

Соотношение (44) совпадает с решением аналогичной задачи во втором приближении, приведенном в [2-4].

Аналогично можно получить решение задачи (8)-(13) и в любом последующем приближении.

Результаты расчетов по формуле (14) в нулевом, первом и восьмом приближениях в сравнении с точным решением [5] даны на рис. 2. Их анализ позволяет заключить, что в диапазоне числа Фурье 0,001 < Бо < 0,01 полученные по формуле (14) значения температур в восьмом приближении отличаются от их точных величин не более чем на 2 %, а для чисел Фурье 0,01 < Бо < да они практически совпадают с точными. При этом из решения характеристического уравнения были получены следующие собственные числа: Х1 =-2,4674011027; Х2 =-22,2066099025; Х3 =-61,6850275455; Х4 =-120,903753244; Х5 =-200,477960523; Х6 =-322,185335578; X7 =-647,223264386; X8 =-2422,85064772.

Точные значения собственных чисел: X1 = 2,467401100272; X2 = = 22,2066099024; X3 = 61,685027506; X4 = 120,902653; X, = 199,859; X6 = 298,555; X 7 = 416,99; X8 = 555,16.

Результаты расчетов по формуле (14) для трехслойной конструкции в сравнении с расчетом по методу прогонки даны на рис. 3. Их анализ позволяет заключить, что отличие температур в восьмом приближении от их значений, полученных по методу прогонки, в диапазоне чисел Фурье 0,187 < Бо < да не превышает 0,5 %.

Применительно к решению задачи для трехслойной конструкции в восьмом приближении получены следующие значения собственных чисел: X1 = 1,614384; X2 = 23,29944; X3 = 61,66503; X4 = 121,1786; X5 = = 223,15143; X6 = 352,645262; X7 = 688,46335; X, = 2552,72532.

При малых и сверхмалых значениях временной переменной до момента времени, пока фронт температурного возмущения, двигаясь от внешней поверхности последнего слоя Е, = 1 (точка задания граничного условия первого рода), не достигает поверхности контакта между предпоследним и последним слоями, теплообмен протекает как бы в однослойной пластине. Для получения аналитического решения применительно к этому (последнему) слою многослойной системы применим метод, основанный на определении фронта температурного возмущения и дополнительных граничных условий [6]. При его использовании процесс теплообмена разделяется на две стадии по времени 0 < Бо < Бо1 и Бо1 < Бо < да. Для этого

вводится движущаяся во времени граница (фронт температурного возмущения), разделяющая исходную область 0 < Ч < 1 (однослойная пластина)

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

было получено выше в виде (14). Следовательно, с использованием данного метода будем находить решение лишь для первой стадии процесса (для малых и сверхмалых значений времени).

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

Щхг) =х8^Щт) (т о < х <8); (45)

дт дх '

Т (х, 0) = Т0; (46)

дТ (0,Т) = 0; (47)

дх

Т (5,т) = Тст, (48)

где с — теплоемкость; у — плотность; Т0 — начальная температура; Тст — температура стенки при х = 5; 5 — толщина пластины; Х — коэффициент теплопроводности.

Задача (45)-(48) в безразмерных переменных будет:

«Ш. (Ро > 0. 0 < 5 < 1); (49)

дРо 34

© (4,0) = 0; (50)

д©(Мо) = 0. (51)

55

© (1,Ро) = 1, (52)

где © = (Т -Т0)/ (Тст -Т0); Ч = х/5; Ро = ат /52; а — коэффициент температуропроводности.

После введения фронта температурного возмущения дДРо) и новой независимой переменной р = 1 - Ч математическая постановка задачи (49)-(52) для первой стадии процесса приводится к виду (рис. 4):

д©(р,Ро) = д 2©(р2Ро) (Ро >0; 0 <р<1); (5з)

ЗРо ф2 '

0 (0, Бо) = 1;

(54)

0 (й,Ро) = 0; д0(д1,Бо)

др

= 0.

(55)

(56)

Отметим, что граничное условие (51) в задаче (53)-(56) не требуется, так как оно не влияет на процесс теплообмена в первой стадии. Соотношения (55), (56) представляют условия тепловой изоляции подвижной границы.

Решение задачи (53)-(56) принимается в виде

0(р,Бо) = £ак (% У

(57)

к=0

где ак (%) — неизвестные коэффициенты, определяемые из граничных условий (54)-(56). После их нахождения соотношение (57) будет

0(р,Бо ) =

%

(58)

1,0 0

Р 1,0

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

Для определения неизвестной функции времени %1(Бо) составим невязку уравнения (53) и проинтегрируем ее в пределах толщины термического слоя

■д0(р,Бо) , д20(р, Бо)

дБо

др2

(59)

Определяя интегралы в (59), находим

= 6аТо. (60)

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

%1(Бо) = 412Бо. (61)

Положив % (Бо1) = 1, находим время достижения фронтом температурного возмущения координаты р = 1, которое будет Бо1 = 0,08333.

0

Соотношения (58), (61) представляют решение задачи (53)-(56) в первом приближении. Анализ результатов расчетов безразмерной температуры по формуле (58) в сравнении с точным решением для однослойной пластины в диапазоне чисел Фурье 0,0001 < Бо < 0,05 позволяет заключить, что расхождение составляет 3-4 %.

Для повышения точности решения задачи (53)-(56) необходимо увеличивать число членов ряда ак (57). Появляющиеся при этом дополнительные неизвестные коэффициенты находятся из дополнительных граничных условий, определяемых с использованием заданных граничных условий (54)-(56) и дифференциального уравнения (53). Их физический смысл состоит в выполнении решением (57) уравнения (53) на границе области (р = 0) и на фронте температурного возмущения дДр). Так как область перемещения фронта температурного возмущения включает весь диапазон изменения пространственной переменной 0 < р < 1, то, следовательно, чем большее число дополнительных граничных условий будет использовано, тем лучше будет выполняться уравнение (53) в этом же диапазоне изменения координаты р. Отметим, что в каждом последующем приближении вводятся три новых дополнительных граничных условия. Общие формулы их получения при любом числе приближений имеют вид:

*©(0, Р0) = 0; *©4 Р0) = 0; ^©^ Р0) = 0 (' = 2,4,6,...), (62)

др' др' др"1 ^

где ' - число приближений.

Найдем решение задачи (53)-(56) во втором приближении. Формулы (62) в этом случае принимают вид:

д 2©(02Бо)=0, д 2©(д/о)=0; д3©(дуБо)=0. (63)

др др др

Используя основные (54)-(56) и дополнительные (63) граничные условия, можно найти уже шесть неизвестных коэффициентов ак (41) (к = 0, 5) ряда (57). Для их определения будем иметь цепочную систему шести алгебраических линейных уравнений, из которой они легко могут быть найдены. После определения коэффициентов ак(41) соотношение (57) во втором приближении принимает вид

©(р, Бо) =

'1 - 3 ^

2 41

, 41 .

4

(64)

Подставляя (64) в (53) и определяя интеграл в пределах от р = 0 до р = 41 (Бо), относительно неизвестной функции 41 (Бо) получаем следующее обыкновенное дифференциальное уравнение:

4^4\ = 10^Бо. (65)

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

^(Ро) = у1 20Ро.

(66)

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

И в пределе при п Ро1 = 0. Этот результат находится в полном соответствии с гипотезой о бесконечной скорости распространения теплоты, лежащей в основе вывода параболического уравнения (53) (рис. 5). Соотношения (64), (66) представляют решение задачи (53)-(56) во втором приближении. Анализ результатов расчетов по формуле (64) позволяет заключить, что во втором приближении происходит значительное повышение точности решения задачи (отклонение от точного не превышает 1-2 %).

Для оценки точности получаемых решений при сверхмалых значениях времени по формуле (57) были выполнены расчеты в третьем, седьмом и четырнадцатом приближениях при Ро = 3 • 10-8. Результаты расчетов в сравнении с точным решением приведены на рис. 6. Их анализ позволяет заключить, что значения температур, полученных по формуле (57), отличаются от точных их значений в третьем приближении на 0,31 %, в седьмом - на 0,03 % и в четырнадцатом - на 0,0004 %. При необходимости число приближений можно увеличить и, следовательно, имеется возможность получения решений практически с заданной степенью точности, без каких-либо ограничений на величину числа Фурье в области малых и сверхмалых его значений.

Можно заметить, что соотношения (58), (64) представлены в виде произведения искомой функции времени q1(Fo) в отрицательной степени на координатную функцию, зависящую лишь от пространственной перемен-

Рис. 5. Кривые перемещения фронта температурного возмущения: 1, 2, 3, 4, 5, 7, 10, 14 - номер приближения

1,000 1 -©

0,996 0,992 0,988 0,984

/'

/ / //X

/ / /// ///

' П '///

//// / /// --4

5

р-104

Рис. 6. Изменение температуры в пластине (Ро = 3-10-8): 1 - третье приближение; 2 - седьмое; 3 - четырнадцатое; 4 - точное решение [5]

6

8

ной, т. е. вид решения формально такой же, как и в методе Л. В. Канторовича (соотношения (33), (44)). Дальнейшие выкладки отличаются лишь тем, что верхним пределом интегрирования в соотношении (59) является переменная величина % (Бо), тогда как в соотношениях (20), (34) этот предел для каждого из контактирующих тел представляется константой. Таким образом, применяя, по существу, один и тот же подход, можно получать простые аналитические выражения для описания температурного состояния одно- и многослойных конструкций во всем диапазоне времени нестационарного процесса (включая малые и сверхмалые его значения) практически с заданной степенью точности.

В Ы В О Д Ы

1. На основе использования ортогонального метода Л. В. Канторовича даны теоретические положения метода построения аналитических решений краевых задач теплопроводности для многослойных конструкций. Применяя глобальную систему неизвестных функций времени, многослойная конструкция приводится к однослойной с переменными (кусочно-однородными) свойствами среды, что позволяет строить аналитические решения путем использования систем координатных функций, точно удовлетворяющих граничным условиям и условиям сопряжения.

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

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

Л И Т Е Р А Т У Р А

1. Б е л я е в, Н. М. Методы нестационарной теплопроводности / Н. М. Беляев,

A. А. Рядно. - М.: Высш. шк., 1982. - 304 с.

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

3. К у д и н о в, В. А. Теплообмен при течении Куэтта с учетом теплоты трения /

B. А. Кудинов, И. В. Кудинов // Энергетика... (Изв. высш. учеб. заведений и энерг. объединений СНГ). - 2011. - № 2. - С. 43-51.

4. Ц о й, П. В. Методы расчета задач тепломассопереноса / П. В. Цой. - М.: Энерго-атомиздат, 1984. - 432 с.

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

6. К у д и н о в, В. А. Методы решения параболических и гиперболических уравнений теплопроводности / В. А. Кудинов, И. В. Кудинов. - М.: Книжный дом «Либроком», 2012. - 280 с.

Представлена кафедрой теоретических

основ теплотехники и гидромеханики Поступила 19.09.2012