• уменьшить загрязняемость трубного пучка по воздушной стороне и увеличить длительность межремонтных периодов. Л И Т Е Р А Т У Р А 1. К е р н, Д. Развитые поверхности теплообмена: пер. с англ. / Д. Керн, А. Краус. – М.: Энергия, 1977. – 464 с. 2. О с н о в ы расчета и проектирования теплообменников воздушного охлаждения: справ. / под общ. ред. В. Б. Кунтыша, А. Н. Бессонного. – СПб.: Недра, 1996. – 512 с. 3. Т е п л о о т д а ч а и аэродинамическое сопротивление шахматных пучков из круг- лых труб с подогнутыми спиральными KLM-ребрами / В. Б. Кунтыш [и др.] // Химическое и нефтегазовое машиностроение. – 2003. – № 11. – С. 10–14. 4. К у н т ы ш, В. Б. Анализ тепловой, объемной и массовой характеристик теплооб- менных секций аппаратов воздушного охлаждения / В. Б. Кунтыш, А. Э. Пиир // Химиче- ское и нефтегазовое машиностроение. – 2009. – № 5. – С. 3–6. 5. К у н т ы ш, В. Б. Основные способы энергетического совершенствования аппаратов воздушного охлаждения / В. Б. Кунтыш, А. Н. Бессонный, А. А. Брилль // Химическое и нефтегазовое машиностроение. – 1997. – № 4. – С. 41–44. 6. Т е п л о о б м е н н а я секция: пат. 2213920 России, МПК С2, F 28 D 3/02 / В. П. Му- лин, И. И. Кочетов, Р. Ф. Теляев, В. Б. Кунтыш, В. И. Мелехов, А. В. Самородов; заявитель ЗАО «Октябрьскхиммаш». – № 2001119695; заявл. 16.07.2001; опубл. 15.07.2003 // Бюл. изобрет. / Роспатент. – 2003. – № 28. – С. 70. 7. У с т р о й с т в о для изготовления теплообменной трубы со спирально-навивными ребрами: заявка на полезную модель Респ. Беларусь, МПК В 21 D 11/06 / В. Б. Кунтыш, В. П. Мулин, Е. С. Санкович, А. Ш. Миннигалеев; заявитель Белорус. гос. технол. ун-т. – № и 20120381; заявл. 05.04.12. 8. Т е п л о о б м е н н а я ребристая труба: пат. 4814 Респ. Беларусь, МПК F 28 F 1/00 / В. Б. Кунтыш, В. И. Володин, Е. С. Санкович, В. П. Мулин, А. Э. Пиир, А. Ш. Миннигалеев, Г. Г. Баранов; заявитель Белорус. гос. технол. ун-т. – № и 20080322; заявл. 17.04.08; опубл. 30.10.08 // Афiцыйны бюл. / Нац. цэнтр iнтэлектуал. уласнасцi. – 2008. – № 5. – С. 36. 9. Т е п л о о б м е н н а я биметаллическая ребристая труба: пат. 14907 Респ. Бела- русь, МПК F 28 F 1/00 / В. Б. Кунтыш, Е. С. Санкович, В. П. Мулин, А. Ш. Миннигалеев, А. Э. Пиир, И. Р. Гаязов, А. Л. Соловьев; заявитель Белорус. гос. технол. ун-т. – № и 20091539; заявл. 28.10.09; опубл. 30.10.11 // Афiцыйны бюл. / Нац. цэнтр iнтэлектуал. уласнасцi. – 2011. – № 4. – С. 58. 10. С п о с о б и устройство для изготовления теплообменной трубы с KLM-ребрами: пат. 16177 Респ. Беларусь, МПК В 21 С 37/15 / В. Б. Кунтыш, В. П. Мулин, Е. С. Санко- вич, А. Э. Пиир, А. Ш. Миннигалеев, А. Л. Соловьев, О. В. Петрович; заявитель Белорус. гос. технол. ун-т. – № а 2010366; заявл. 11.03.10; опубл. 30.08.12 // Афiцыйны бюл. / Нац. цэнтр iнтэлектуал. уласнасцi. – 2012. – № 4. – С. 44. Представлена кафедрой энергосбережения, гидравлики и теплотехники Поступила 12.11.2012 УДК 536.2 (075) ОРТОГОНАЛЬНЫЕ МЕТОДЫ В ЗАДАЧАХ ТЕПЛОПРОВОДНОСТИ ДЛЯ МНОГОСЛОЙНЫХ КОНСТРУКЦИЙ Докт. физ.-мат. наук КУДИНОВ В. А., инженеры КОТОВА Е. В., ЕРЕМИН А. В., КУЗНЕЦОВА А. Э. Самарский государственный технический университет Для определения требуемого сочетания свойств многослойных конст- рукций наилучшим образом подходят аналитические (приближенные ана- 44 литические) решения. Трудности получения аналитических решений крае- вых задач для многослойных конструкций связаны с необходимостью вы- полнения условий сопряжения в виде равенства температур и тепловых потоков в точках контакта слоев. Выполнение этих условий в случае при- менения точных аналитических методов приводит к необходимости реше- ния характеристической системы в виде цепочного трансцендентного уравнения относительно собственных чисел краевой задачи, решение кото- рого может быть получено лишь численными методами [1, 2]. Из приближенных аналитических методов решения достаточно широ- кое распространение получили ортогональные методы взвешенных невязок (Л. В. Канторовича, Бубнова – Галеркина и др.) [2–4]. Однако их примене- ние к расчетам многослойных конструкций существенно сдерживается не- обходимостью построения систем координатных функций, удовлетворяю- щих граничным условиям и условиям сопряжения. В настоящей работе рассматривается метод, позволяющий строить системы координатных функций, в любом приближении точно удовлетворяющих граничным ус- ловиям и условиям сопряжения, независимо от числа контактирующих тел. Построение таких систем оказалось возможным благодаря принятию гло- бальной системы неизвестных функций времени (одинаковой для всех кон- тактирующих тел), введение которой позволяет приводить многослойную конструкцию к однослойной с переменными (кусочно-однородными) свой- ствами среды. Для получения решений при малых и сверхмалых значениях временной переменной необходимо использовать большое число приближений, что в ортогональных методах приводит к необходимости решения алгебраиче- ских уравнений высоких степеней (относительно собственных чисел крае- вой задачи), точность решения которых даже при использовании современ- ной компьютерной техники не всегда может удовлетворять потребностям практики. Наиболее эффективным в области малых значений временной пе- ременной является метод, основанный на определении фронта температурно- го возмущения и дополнительных граничных условий [3, 6]. В настоящей статье он применен для расчета температурного состояния последнего слоя многослойной системы, так как при малых и особенно сверхмалых значе- ниях времени теплообмен протекает лишь в этом слое. В качестве конкретного примера найдем решение задачи теплопровод- ности для трехслойной пластины в следующей математической постановке (рис. 1): 2 2 ( , τ) ( , τ) τ i i i T x T xa x ∂ ∂ = ∂ ∂ (1) ( )1 0 3τ 0; ; 1, 2, 3; 0; δ ;i ix x x i x x−> < < = = = 0( ,0) ;iT x T= (2) 1(0, τ) 0;T x ∂ = ∂ (3) ст3 (δ, τ) ;T T= (4) 45 1( , τ) ( , τ)i i i iT x T x+= ( 1,2);i = (5) 1 1λ ( , τ) λ ( , τ)i i i i i iT x T x x x + +∂ ∂= ∂ ∂ ( 1,2),i = (6) где Ti – температура i-го слоя; x – координата; τ – время; δ1, δ2, δ3 – тол- щины слоев; δ – суммарная толщина трехслойной пластины; ai – коэффи- циент температуропроводности i-го слоя; T0 – начальная температура; Tст – температура стенки при х = δ; λi – коэффициент теплопроводности i-го слоя. Рис. 1. Расчетная схема трехслойной пластины Введем следующие безразмерные переменные: ст ст0 ;ii T T T T − Θ = − ξ ; δ i i x = 2 τFo , δ a = (7) где Θ – относительная избыточная температура; ξ – безразмерная коорди- ната; Fo – число Фурье; a – наименьший из коэффициентов температуро- проводности ( 1, 2, 3).ia i = Задача (1)–(6) с учетом (7) приводится к виду: 2 2 (ξ, Fo) (ξ, Fo) Fo ξ i i ia a ∂Θ ∂ Θ = ∂ ∂ (8) ( )1 0 3Fo 0; ξ ξ ξ ; 1, 2, 3; ξ 0; ξ 1 ;i i i−> < < = = = (ξ, 0) 1;iΘ = (9) 1(0,Fo) 0; ξ ∂Θ = ∂ (10) x 2x 2δ T 1 1 1 δ λ а 1T 0 CTT 3x 1x 2 2 λ а 2T 3δ 3 3 λ а 3T δ1 2 3 1 2 3 46 3 (1,Fo) 0;Θ = (11) 1(ξ ,Fo) (ξ ,Fo) ( 1,2);i i i i i+Θ = Θ = (12) 1 1λ (ξ , Fo) λ (ξ , Fo) ( 1,2). ξ ξ i i i i i i i+ +∂Θ ∂Θ= = ∂ ∂ (13) Решение задачи (8)–(13), следуя ортогональному методу Л. В. Канторо- вича, принимается в виде 1 (ξ, Fo) (Fo) (ξ) n i k ki k f = Θ = ϕ∑ ( 1, 2, 3),i = (14) где (Fo)kf – неизвестные функции; (ξ) ( 1, 2, 3)ki iϕ = – координатные функции, определяемые таким образом, чтобы заранее точно выполнялись граничные условия (10), (11) и условия сопряжения (12), (13). Отличительной особенностью используемого здесь метода является принятие одинаковой системы неизвестных функций времени (Fo)kf для каждого из контактирующих тел, что позволяет получать наиболее просто- го вида координатные функции (ξ) ( 1, 2, 3; 1, ).ki i k nϕ = = Способ построения систем координатных функций, точно удовлет- воряющих граничным условиям и условиям сопряжения, заключается в следующем. Вначале строится координатная функция для третьего слоя. Она принимается в таком виде, чтобы точно выполнялось граничное условие (11) 2 3 ( ) 1 ( 1, ). k k k nϕ ξ = − ξ = (15) Координатная функция для второго слоя принимается в виде 2 2 1 2( ) , k k A Aϕ ξ = + ξ (16) где неизвестные коэффициенты A1 и A2 находятся из условий сопряжения (12), (13) между вторым и третьим слоями. После их определения соотно- шение (16) принимает вид 2 23 3 2 2 2 2 ( ) 1 1 k kk  λ λ ϕ ξ = − − ξ − ξ λ λ  ( 1, ).k n= (17) Координатная функция для первого слоя принимается в виде 2 1 1 2( ) , k k B Bϕ ξ = + ξ (18) где В1, В2 – неизвестные коэффициенты. Очевидно, что соотношение (14) при использовании в качестве коорди- натной функции (18) точно удовлетворяет граничному условию (10). Неиз- вестные коэффициенты В1, В2 находятся из условий сопряжения (12), (13) между первым и вторым слоями. После их определения соотношение (18) будет 47 2 2 23 3 3 3 1 1 2 2 1 2 1 ( ) 1 1k k kk    λ λ λ λ ϕ ξ = − − ξ − − ξ − ξ   λ λ λ λ    ( 1, ).k n= (19) Непосредственной подстановкой можно убедиться, что при использо- вании в качестве координатных функций соотношений (15), (17), (19) со- отношение (14) в любом приближении ( 1, )k n= точно удовлетворяет гра- ничным условиям (10), (11) и условиям сопряжения (12), (13) при любых значениях неизвестных функций (Fo).kf Найдем решение задачи (8)–(13) в нулевом приближении. Для этого по- требуем, чтобы соотношение (14), в котором ограничимся одним членом ряда, удовлетворяло не уравнению (8), а некоторым осредненным уравне- ниям (проинтегрированным в пределах толщины каждого слоя) 1 23 2 1 ( ,Fo) ( ,Fo) 0. Fo i i i i i i a d a − ξ = ξ  ∂Θ ξ ∂ Θ ξ − ξ = ∂ ∂ξ  ∑ ∫ (20) Подставляя (14) в (20), находим 1 2 1 2 1 2 1 2 1 2 2 2 1 1 2 3 0 1 31 2 1 1 2 0 (Fo) ( ) ( ) (1 ) 2 (Fo) 0, f D D d D D d d aa af D d D d d a a a ξ ξ ξ ξ ξ ξ ξ ξ   ′ − ξ ξ + − ξ ξ + − ξ ξ +       + ξ + ξ + ξ =     ∫ ∫ ∫ ∫ ∫ ∫ (21) где 1f ′ – первая производная от функции 1(Fo)f по переменной Fo. Определяя интегралы в (21), относительно неизвестной функции 1(Fo)f получаем следующее обыкновенное дифференциальное уравнение: 1 1 (Fo) (Fo) 0, Fo df f d + ν = (22) где 2 1 ;µν = µ 2 3 3 331 1 1 1 2 2 1 2 1 2 2 1 2( ) ( ) (1 ) (1 ) ; 3 3 3 3 DDD D µ = ξ − ξ + ξ − ξ − ξ − ξ + − ξ − − ξ +    ( )2 1 1 2 3 2 1 3 22 ( ) (1 ) / ;a D a D a aµ = ξ + ξ − ξ + − ξ 2 23 3 3 1 2 2 1 2 1 1 ;k kD    λ λ λ = − − ξ − − ξ   λ λ λ    3 1 1 ;D λ= λ 232 2 2 1 1 ;kD  λ = − − ξ λ  33 2 .D λ= λ Интегрируя уравнение (22), находим 48 1(Fo) exp( Fo),f C= ν (23) где С – постоянная интегрирования. Подставляя (23) в (14), получаем 1( ,Fo) exp( Fo) ( ) ( 1, 2, 3).i iC iΘ ξ = ν ϕ ξ = (24) Для определения постоянной интегрирования С составим невязку на- чального условия (9) и проинтегрируем ее в пределах толщины каждо- го слоя [ ] 1 3 1 ( ,0) 1 0. i i i i d − ξ = ξ Θ ξ − ξ =∑ ∫ (25) Подставляя (24) в (25), находим 1 2 1 2 1 2 2 2 1 2 3 0 ( ) 1 ( ) 1 (1 ) 1 0.С D D d С D D d С d ξ ξ ξ ξ      − ξ − ξ + − ξ − ξ + − ξ − ξ =     ∫ ∫ ∫ Определяя интегралы, относительно постоянной интегрирования С по- лучаем алгебраическое линейное уравнение. Его решение С = 1/µ1. После определения постоянной интегрирования решение задачи (8)–(13) в нулевом приближении принимает вид 1 1 exp( Fo) ( )( , Fo) ii ν ϕ ξ Θ ξ = µ ( 1, 2, 3).i = (26) Если положить 1 2 3 ;λ = λ = λ = λ 1 2 3a a a a= = = (однослойная пласти- на), то соотношение (26) будет 2( , Fo) 1,5exp( 3 Fo)(1 ).Θ ξ = − − ξ (27) Последняя формула совпадает с решением в нулевом приближении для однослойной пластины, приведенным в [2, 3]. Результаты расчетов по формуле (27) представлены на рис. 2. Используя формулу (26), найдем решение конкретной задачи для трех- слойной пластины при следующих исходных данных: 1 0,00086 м;x = 2 0,00261м;x = 3 0,00502 м;x = λ1 = 11 Вт/(м⋅К); 2 2 Вт/(м К);λ = ⋅ 3 1,1 Вт/(м К);λ = ⋅ 6 2 1 3,624 10 м /с;a −= ⋅ 6 2 2 1,02 10 м /с;a −= ⋅ 6 23 0,94 10 м /с.a −= ⋅ (28) Анализ результатов расчетов по формуле (26) в сравнении расчетом численным методом прогонки позволяет заключить, что в диапазоне чисел Фурье 0,187 ≤ Fo < ∞ расхождение составляет около 12 % (рис. 3). Для повышения точности найдем решение задачи (8)–(13) в первом приближении. Для этого составим невязки уравнений (8) и потребуем ор- тогональности невязок к координатным функциям первого приближения 49 1 ( 1, 2, 3)i iϕ = (учитываем, что согласно (27) наименьшим является коэф- фициент температуропроводности а3) 1 ξ 23 11 1 1 12 1 3ξ ( )(Fo) ( ) (Fo) ( ) 0. Fo i i i i i i i af f d a − =  ∂ ϕ ξ∂ ϕ ξ − ϕ ξ ξ = ∂ ∂ξ  ∑ ∫ 0 0,2 0,4 0,6 0,8 1,0 0,2 0,4 0,6 0,8 1,0 Fo = 0,93 Θ 0,5 0,187 0,05 0,03 0,0187 0,005 0,00187 0,0001 ξ 0,3 Рис. 2. Распределение температуры в трехслойной пластине, приведенной к однослойной:  – по формуле (26) (нулевое приближение); × − по формуле (33) (первое приближение); —— – по формуле (14) (восьмое приближение);  – точное решение [5] 50 1,0 Θ 0,8 0,6 0,4 0,2 0 ξ1 0,2 ξ2 0,6 0,8 ξ 1,0 1,87 0,93 Fo= 0,187 0,5 0,005 0,0003 0.5 0.7 0,015 0,0015 Рис. 3. Распределение температуры в трехслойной пластине: × – по формуле (26) (нулевое приближение);  – по формуле (31) (первое приближение); —— – по формуле (14) (восьмое приближение);  – метод прогонки;  – точное решение [5] Последнее соотношение можно переписать в виде 1 2 1 2 1 2 1 2 1 2 2 2 1 1 11 2 3 12 13 0 1 1 2 1 1 11 2 12 13 3 30 (Fo) ( ) ( ) (1 ) 2 (Fo) 0. f D D d D D d d a af D d D d d a a ξ ξ ξ ξ ξ ξ ξ ξ   ′ − ξ ϕ ξ + − ξ ϕ ξ + − ξ ϕ ξ +       + ϕ ξ + ϕ ξ + ϕ ξ =      ∫ ∫ ∫ ∫ ∫ ∫ Определяя интегралы в последнем соотношении, относительно неиз- вестной функции 1(Fo)f получаем обыкновенное дифференциальное урав- нение вида 1 1 1(Fo) / Fo (Fo) 0,df d f+ ν = (29) где 1 4 3/ ;ν = µ µ 3 0,516004 ;µ = 4 1,042163.µ = Интегрируя уравнение (29), находим 1 1 1(Fo) exp( Fo),f C= ν (30) где C1 – постоянная интегрирования. Подставляя (30) в (14), получаем 1 1 1( ,Fo) exp( Fo) ( )i iCΘ ξ = ν ϕ ξ ( 1, 2, 3).i = (31) Для определения постоянной интегрирования C1 составим невязку на- чального условия (9) и потребуем ортогональности невязки к координат- ным функциям 1 ( 1, 2, 3)i iϕ = 51 1 2 1 2 2 2 1 1 11 1 2 3 12 0 1 2 1 13 ( ) 1 ( ) 1 (1 ) 1 0. С D D d С D D d С d ξ ξ ξ ξ    − ξ − ϕ ξ + − ξ − ϕ ξ +     + − ξ − ϕ ξ =  ∫ ∫ ∫ (32) Определяя интегралы в (32), относительно постоянной интегрирования C1 получаем алгебраическое линейное уравнение. Его решение 1 31/ .C = µ После определения постоянной интегрирования C1 решение задачи (8)–(13) в первом приближении находится из (31). Если положить 1 2 3 ;λ = λ = λ = λ 1 2 3 ,a a a a= = = то соотношение (31) приводится к виду 2( ,Fo) 1,25 exp( 2,5 Fo)(1 ).Θ ξ = − − ξ (33) Формула (33) совпадает с решениями в первом приближении для одно- слойной пластины, приведенными в [2–4]. Результаты расчетов по формулам (26), (27), (31), (33) даны на рис. 2, 3. Их анализ позволяет заключить, что точность решения в первом прибли- жении по сравнению с нулевым существенно повышается. И в частности, расхождение с точным решением (рис. 2) не превышает 3 %. Для дальнейшего повышения точности найдем решение задачи (8)–(13) во втором приближении. Составляя невязку уравнений (8) и требуя орто- гональности невязки к координатным функциям 1 ( )iϕ ξ и 2 ( ) ( 1, 2, 3),i iϕ ξ = находим 1 ξ3 1 1 2 2 1 1 2 2 1 3ξ ( ) 0 ( 1, 2), i i i i i i i ji i af f f f d j a − =   ′ ′ ′′ ′′ϕ + ϕ − ϕ + ϕ ϕ ξ = =    ∑ ∫ (34) где 1 2,i i′′ ′′ϕ ϕ – вторые производные от функций 1 ( )iϕ ξ и 2 ( )iϕ ξ по перемен- ной ξ. Соотношение (34) можно переписать в виде: 1 2 1 2 1 2 1 2 1 1 11 2 21 11 1 12 2 22 12 1 13 2 23 13 0 1 1 2 1 11 2 21 11 1 12 2 22 12 1 13 2 23 13 3 30 ( ) ( ) ( ) ( ) ( ) ( ) 0; f f d f f d f f d a af f d f f d f f d a a ξ ξ ξ ξ ξ ξ ξ ξ ′ ′ ′ ′ ′ ′ϕ + ϕ ϕ ξ + ϕ + ϕ ϕ ξ + ϕ + ϕ ϕ ξ − ′′ ′′ ′′ ′′ ′′ ′′− ϕ + ϕ ϕ ξ − ϕ + ϕ ϕ ξ − ϕ + ϕ ϕ ξ = ∫ ∫ ∫ ∫ ∫ ∫ 1 2 1 2 1 2 1 2 1 1 11 2 21 21 1 12 2 22 22 1 13 2 23 23 0 1 1 2 1 11 2 21 21 1 12 2 22 22 1 13 2 23 23 3 30 ( ) ( ) ( ) ( ) ( ) ( ) 0. f f d f f d f f d a af f d f f d f f d a a ξ ξ ξ ξ ξ ξ ξ ξ ′ ′ ′ ′ ′ ′ϕ + ϕ ϕ ξ + ϕ + ϕ ϕ ξ + ϕ + ϕ ϕ ξ − ′′ ′′ ′′ ′′ ′′ ′′− ϕ + ϕ ϕ ξ − ϕ + ϕ ϕ ξ − ϕ + ϕ ϕ ξ = ∫ ∫ ∫ ∫ ∫ ∫ Определяя интегралы в последних соотношениях, относительно неиз- вестных функций 1(Fo)f и 2 (Fo)f с учетом исходных данных (27) получаем следующую систему двух обыкновенных дифференциальных уравнений: 52 1 1 2 2 3 1 4 2 1 1 2 2 3 1 4 2 ' ' 0; ' ' 0, b f b f b f b f l f l f l f l f + + + =  + + + =  (35) где b1 = 0,494399; b2 = 0,585074; b3 = 0,958027; b4 = 1,49522; l1 = 0,585074; l2 = 0,702963; l3 = 1,22521; l4 = 2,17489. Частные решения системы уравнений (35) разыскиваются в виде: 1 2(Fo) exp( Fo); (Fo) exp( Fo),f A r f D r= = (36) где A, D, r – некоторые постоянные. Подставляя (36) в (35), находим: (0,494399 0,958027) (0,585074 1,49522) 0; (0,585074 1,22521) (0,702963 2,17489) 0. A r D r A r D r + + + =  + + + =  (37) Система однородных алгебраических уравнений (37) имеет нетриви- альное решение в случае, если определитель ее равен нулю: 0,494399 0,958027 0,585074 1,4952 0. 0,585074 1,22521 0,702963 2,17489 r r r r + + = + + (38) Из (38) получаем следующее характеристическое уравнение: 0,005232r2 + 0,157067 0,25164 0.r+ + = Его корни: 1 1,698214;r = − 2 28,320207.r = − Подставляя величину корня r1 в систему уравнений (37), получаем: 1 1 1 1 0,113714 0,475086 0; 0,230316 0,962233 0. A D A D + =   + =  (39) Для однородной системы (39) можно положить 1 const 1.A = = Тогда D1 = –0,239355. Аналогично, подставляя r2 в (37), находим: A2 = 1; D2 = –0,881653. С учетом найденных значений постоянных , , ( 1, 2, 3)i i iA D r i = соотно- шения (36) примут вид: 1 2 1 2Fo Fo Fo Fo 1 1 2 2 1 2(Fo) ; (Fo) . r r r rf A e A e f D e D e= + = + (40) Чтобы найти общее решение системы уравнений (35), умножим частное решение, включающее корень r1, на произвольную постоянную C1, а ре- шение, включающее корень r2, – на постоянную C2. Подставляя получен- ные общие решения в (14) (при n = 2), находим 1 2 1 2Fo Fo Fo Fo 1 1 2 2 1 1 1 2 2 2( ,Fo) ( ) ( ) r r r r i i iC A e C A e C D e C D eΘ ξ = + ϕ + + ϕ ( 1, 2, 3).i = (41) Для определения постоянных C1 и C2 составляется невязка начального условия (9) и требуется ортогональность невязки к координатным функци- ям 1 ( )iϕ ξ и 2 ( )iϕ ξ ( 1, 2, 3)i = [ ] 1 2 1 1 2 2 1 1 1 2 2 2 1 ( ) ( ) 1 0 i i i i ji i C A C A C D C D d − ξ = ξ + ϕ + + ϕ − ϕ ξ =∑ ∫ ( 1, 2, 3).j = (42) 53 Соотношение (42) относительно неизвестных C1 и C2 приводится к следующей системе двух алгебраических линейных уравнений: 1 2 1 2 2 1 2 1 11 2 21 11 1 12 2 22 12 0 1 1 13 2 23 13 1 1 11 2 21 21 1 12 2 22 22 0 1 1 13 2 23 23 ( 1) ( 1) ( 1) 0; ( 1) ( 1) ( 1) 0, E E d E E d E E d E E d E E d E E d ξ ξ ξ ξ ξ ξ ξ  ϕ + ϕ − ϕ ξ + ϕ + ϕ − ϕ ξ +   + ϕ + ϕ − ϕ ξ =   ϕ + ϕ − ϕ ξ + ϕ + ϕ − ϕ ξ +   + ϕ + ϕ − ϕ ξ =  ∫ ∫ ∫ ∫ ∫ ∫ где 1 1 1 2 2 2 1 1 2 2; .E C A C A E C D C D= + = + Определяя интегралы, находим: 1 2 1 2 0,356267 0,214339 0,645874 0; 0,4191102 0,346957 0,795890 0. C C C C + + =   + + =  (43) Из решения системы уравнений (43) получаем: C1 = 1,642377; C2 = −3,631466. После определения постоянных C1 и C2 решение задачи (8)–(13) во втором приближении находится из (41). Если положить 1 2 3 ;λ = λ = λ = λ 1 2 3a a a a= = = (однослойная пластина), то соотношение (41) приводится к виду 2,47Fo 25,52Fo 2 2,47Fo 25,52Fo 4 2 ( ,Fo) (1,55 3,3 )(1 ) (0,28 2,9048 )(1 ). e e e A e − − − − Θ ξ = + − ξ − − + − ξ (44) Соотношение (44) совпадает с решением аналогичной задачи во втором приближении, приведенном в [2–4]. Аналогично можно получить решение задачи (8)–(13) и в любом после- дующем приближении. Результаты расчетов по формуле (14) в нулевом, первом и восьмом приближениях в сравнении с точным решением [5] даны на рис. 2. Их ана- лиз позволяет заключить, что в диапазоне числа Фурье 0,001 Fo 0,01≤ < полученные по формуле (14) значения температур в восьмом приближении отличаются от их точных величин не более чем на 2 %, а для чисел Фурье 0,01 Fo≤ < ∞ они практически совпадают с точными. При этом из решения характеристического уравнения были получены следующие собственные числа: 1 2,4674011027;λ = − 2 22,2066099025;λ = − 3 61,6850275455;λ = − 4 120,903753244;λ = − 5 200,477960523;λ = − 6 322,185335578;λ = − 7 647,223264386;λ = − 8 2422,85064772.λ = − 54 Точные значения собственных чисел: 1 2,467401100272;λ = λ2 = 22,2066099024;= 3 61,685027506;λ = 4 120,902653;λ = 5 199,859;λ = 6 298,555;λ = 7 416,99;λ = 8 555,16.λ = Результаты расчетов по формуле (14) для трехслойной конструкции в сравнении с расчетом по методу прогонки даны на рис. 3. Их анализ по- зволяет заключить, что отличие температур в восьмом приближении от их значений, полученных по методу прогонки, в диапазоне чисел Фурье 0,187 Fo≤ < ∞ не превышает 0,5 %. Применительно к решению задачи для трехслойной конструкции в восьмом приближении получены следующие значения собственных чи- сел: 1 1,614384;λ = 2 23,29944;λ = 3 61,66503;λ = 4 121,1786;λ = λ5 = 223,15143;= 6 352,645262;λ = 7 688,46335;λ = 8 2552,72532.λ = При малых и сверхмалых значениях временной переменной до момента времени, пока фронт температурного возмущения, двигаясь от внешней поверхности последнего слоя ξ = 1 (точка задания граничного условия пер- вого рода), не достигает поверхности контакта между предпоследним и последним слоями, теплообмен протекает как бы в однослойной пласти- не. Для получения аналитического решения применительно к этому (по- следнему) слою многослойной системы применим метод, основанный на определении фронта температурного возмущения и дополнительных гра- ничных условий [6]. При его использовании процесс теплообмена разделя- ется на две стадии по времени 10 Fo Fo< ≤ и 1Fo Fo .≤ < ∞ Для этого вво- дится движущаяся во времени граница (фронт температурного возмуще- ния), разделяющая исходную область 0 1≤ ξ ≤ (однослойная пластина) на две подобласти 10 (Fo)q≤ ξ ≤ и 1(Fo) 1,q ≤ ξ ≤ где 1(Fo)q − функция, опре- деляющая продвижение границы раздела во времени. При этом в области, расположенной за пределами фронта температурного возмущения, сохра- няется начальная температура. Первая стадия процесса заканчивается при достижении движущейся границей поверхности контакта между предпо- следним и последним слоями, т. е. когда 1Fo Fo .= Во второй стадии (ста- дия регулярного режима) изменение температуры происходит по всему объему тела ( )0 1 .≤ ξ ≤ Отметим, что для этой стадии процесса решение было получено выше в виде (14). Следовательно, с использованием данно- го метода будем находить решение лишь для первой стадии процесса (для малых и сверхмалых значений времени). Математическая постановка задачи в данном случае имеет вид: 2 2 ( , τ) ( , τ) τ T x T xc x ∂ ∂ γ = λ ∂ ∂ ( )τ 0; 0 ;x> < < δ (45) 0( ,0) ;T x T= (46) (0, τ) 0;T x ∂ = ∂ (47) ст(δ, τ) ,T T= (48) 55 где с − теплоемкость; γ − плотность; Т0 − начальная температура; Тст − температура стенки при х = δ; δ − толщина пластины; λ − коэффициент те- плопроводности. Задача (45)–(48) в безразмерных переменных будет: 2 2 (ξ,Fo) (ξ,Fo) Fo ξ ∂Θ ∂ Θ = ∂ ∂ ( )Fo 0; 0 ξ 1 ;> < < (49) (ξ, 0) 0;Θ = (50) (0,Fo) 0; ξ ∂Θ = ∂ (51) (1,Fo) 1,Θ = (52) где ( ) ( )0 ст 0/ ;T T T TΘ = − − / ;xξ = δ 2Fo / ;a= τ δ а − коэффициент темпе- ратуропроводности. После введения фронта температурного возмущения 1(Fo)q и новой независимой переменной 1ρ = − ξ математическая постановка задачи (49)–(52) для первой стадии процесса приводится к виду (рис. 4): 2 2 ( ,Fo) ( ,Fo) Fo ∂Θ ρ ∂ Θ ρ = ∂ ∂ρ ( )Fo 0; 0 1 ;> < ρ < (53) (0, Fo) 1;Θ = (54) 1( ,Fo) 0;qΘ = (55) 1( ,Fo) 0.q∂Θ = ∂ρ (56) Отметим, что граничное условие (51) в задаче (53)–(56) не требуется, так как оно не влияет на процесс теплообмена в первой стадии. Соотно- шения (55), (56) представляют условия тепловой изоляции подвижной границы. Решение задачи (53)–(56) принимается в виде ( ) ( )1 0 ,Fo , n k k k a q = Θ ρ = ρ∑ (57) где ( )1ka q − неизвестные коэффициенты, определяемые из граничных ус- ловий (54)–(56). После их нахождения соотношение (57) будет ( ) 2 1 ,Fo 1 . q  ρ Θ ρ = −    (58) 56 1,0 Θ ρ 1,0 0 q1(Fo) Рис. 4. Расчетная схема теплообмена Для определения неизвестной функции времени 1(Fo)q составим не- вязку уравнения (53) и проинтегрируем ее в пределах толщины термиче- ского слоя 1 1 2 2 0 0 ( ,Fo) ( ,Fo) . Fo q q d d∂Θ ρ ∂ Θ ρρ = ρ ∂ ∂ρ∫ ∫ (59) Определяя интегралы в (59), находим 1 1 6 Fo.q dq d= (60) Интегрируя уравнение (46), при начальном условии 1(0) 0q = получаем 1(Fo) 12Fo.q = (61) Положив 1 1(Fo ) 1,q = находим время достижения фронтом температур- ного возмущения координаты ρ = 1, которое будет 1Fo 0,08333.= Соотношения (58), (61) представляют решение задачи (53)–(56) в пер- вом приближении. Анализ результатов расчетов безразмерной температу- ры по формуле (58) в сравнении с точным решением для однослойной пла- стины в диапазоне чисел Фурье 0,0001 < Fo < 0,05 позволяет заключить, что расхождение составляет 3–4 %. Для повышения точности решения задачи (53)–(56) необходимо увели- чивать число членов ряда 1( )ka q (57). Появляющиеся при этом дополни- тельные неизвестные коэффициенты находятся из дополнительных гра- ничных условий, определяемых с использованием заданных граничных условий (54)–(56) и дифференциального уравнения (53). Их физический смысл состоит в выполнении решением (57) уравнения (53) на границе об- ласти ( 0)ρ = и на фронте температурного возмущения 1( ).q ρ Так как область перемещения фронта температурного возмущения включает весь диапазон изменения пространственной переменной 0 1,≤ ρ ≤ то, следо- вательно, чем большее число дополнительных граничных условий будет использовано, тем лучше будет выполняться уравнение (53) в этом же диа- пазоне изменения координаты ρ. Отметим, что в каждом последующем приближении вводятся три новых дополнительных граничных условия. Общие формулы их получения при любом числе приближений имеют вид: 57 (0, Fo) 0; i i ∂ Θ = ∂ρ 1( , Fo) 0; i i q∂ Θ = ∂ρ 1 1 1 ( , Fo) 0 i i q+ + ∂ Θ = ∂ρ ( 2, 4, 6, ...),i = (62) где i – число приближений. Найдем решение задачи (53)–(56) во втором приближении. Формулы (62) в этом случае принимают вид: 2 2 (0,Fo) 0;∂ Θ = ∂ρ 2 1 2 ( ,Fo) 0;q∂ Θ = ∂ρ 3 1 3 ( ,Fo) 0.q∂ Θ = ∂ρ (63) Используя основные (54)–(56) и дополнительные (63) граничные усло- вия, можно найти уже шесть неизвестных коэффициентов 1( ) ( 0, 5)ka q k = ряда (57). Для их определения будем иметь цепочную систему шести ал- гебраических линейных уравнений, из которой они легко могут быть най- дены. После определения коэффициентов 1( )ka q соотношение (57) во вто- ром приближении принимает вид 4 1 1 3( , Fo) 1 1 . 2 q q   ρ ρ Θ ρ = − −      (64) Подставляя (64) в (53) и определяя интеграл в пределах от ρ = 0 до 1(Fo),qρ = относительно неизвестной функции 1(Fo)q получаем следую- щее обыкновенное дифференциальное уравнение: 1 1 10 Fo.q dq d= (65) Интегрируя уравнение (65), при начальном условии 1(0) 0q = находим 1(Fo) 20Fo.q = (66) Положив в (66) 1 1(Fo ) 1,q = находим время окончания первой стадии процесса во втором приближении 1Fo 0,05.= Результаты расчетов пере- мещения фронта температурного возмущения по координате ρ во времени Fo позволяют заклю- чить, что с увеличением числа приближений величина времени 1Fo , при котором фронт темпе- ратурного возмущения достигает координаты 1,ρ = уменьшается. И в пределе при n →∞ 1Fo 0.= Этот результат находится в пол- ном соответствии с гипотезой о бесконечной скорости распро- странения теплоты, лежащей в основе вывода параболического уравнения (53) (рис. 5). Рис. 5. Кривые перемещения фронта температурного возмущения: 1, 2, 3, 4, 5, 7, 10, 14 – номер приближения Fo⋅102 q1(Fo) 58 Соотношения (64), (66) представляют решение задачи (53)–(56) во втором приближении. Анализ результатов расчетов по формуле (64) позволяет заключить, что во втором приближении происходит значительное повыше- ние точности решения задачи (отклонение от точного не превышает 1–2 %). Для оценки точности полу- чаемых решений при сверх- малых значениях времени по формуле (57) были выполне- ны расчеты в третьем, седьмом и четырнадцатом приближениях при 8Fo 3 10 .−= ⋅ Результаты рас- четов в сравнении с точным решением приведены на рис. 6. Их анализ позволяет заключить, что значения температур, полу- ченных по формуле (57), отли- чаются от точных их значений в третьем приближении на 0,31 %, в седьмом – на 0,03 % и в че- тырнадцатом – на 0,0004 %. При необходимости число приближений можно увеличить и, следова- тельно, имеется возможность получения решений практически с заданной степенью точности, без каких-либо ограничений на величину числа Фурье в области малых и сверхмалых его значений. Можно заметить, что соотношения (58), (64) представлены в виде про- изведения искомой функции времени 1(Fo)q в отрицательной степени на координатную функцию, зависящую лишь от пространственной перемен- ной, т. е. вид решения формально такой же, как и в методе Л. В. Канторо- вича (соотношения (33), (44)). Дальнейшие выкладки отличаются лишь тем, что верхним пределом интегрирования в соотношении (59) является переменная величина 1(Fo),q тогда как в соотношениях (20), (34) этот пре- дел для каждого из контактирующих тел представляется константой. Таким образом, применяя, по существу, один и тот же подход, можно по- лучать простые аналитические выражения для описания температурного состояния одно- и многослойных конструкций во всем диапазоне времени нестационарного процесса (включая малые и сверхмалые его значения) практически с заданной степенью точности. В Ы В О Д Ы 1. На основе использования ортогонального метода Л. В. Канторовича даны теоретические положения метода построения аналитических решений краевых задач теплопроводности для многослойных конструкций. Приме- няя глобальную систему неизвестных функций времени, многослойная конструкция приводится к однослойной с переменными (кусочно-одно- родными) свойствами среды, что позволяет строить аналитические реше- ния путем использования систем координатных функций, точно удовле- творяющих граничным условиям и условиям сопряжения. 0,996 1 0,992 0,988 0,984 5 6 410⋅ρ 8 - 1 - 2 - 3 - 4 Θ−1 Рис. 6. Изменение температуры в пластине (Fo = 3·10−8): 1 – третье приближение; 2 – седьмое; 3 – четырнадцатое; 4 – точное решение [5] 1,000 59 2. Разработана методика построения систем координатных функций, в любом приближении точно удовлетворяющих граничным условиям и условиям сопряжения при любом числе контактирующих тел. Для их по- строения используется рекуррентный метод, при котором координатные функции строятся последовательно, переходя от слоя к слою, начиная с последнего, при использовании всякий раз метода неопределенных ко- эффициентов. Реализация такого метода построения координатных функ- ций возможна лишь при использовании глобальной системы неизвестных функций времени. 3. На основе введения фронта температурного возмущения и дополни- тельных граничных условий применительно к последнему слою много- слойной системы получено аналитическое решение задачи теплопроводно- сти, позволяющее выполнять оценку температурного состояния конструк- ции для малых и сверхмалых значений времени. Л И Т Е Р А Т У Р А 1. Б е л я е в, Н. М. Методы нестационарной теплопроводности / Н. М. Беляев, А. А. Рядно. – М.: Высш. шк., 1982. – 304 с. 2. К у д и н о в, В. А. Аналитические решения задач тепломассопереноса и термоупру- гости для многослойных конструкций / В. А. Кудинов, Э. М. Карташов, В. В. Калашников. – М.: Высш. шк., 2005. – 429 с. 3. К у д и н о в, В. А. Теплообмен при течении Куэтта с учетом теплоты трения / В. А. Кудинов, И. В. Кудинов // Энергетика… (Изв. высш. учеб. заведений и энерг. объеди- нений СНГ). – 2011. – № 2. – С. 43–51. 4. Ц о й, П. В. Методы расчета задач тепломассопереноса / П. В. Цой. – М.: Энерго- атомиздат, 1984. – 432 с. 5. Л ы к о в, А. В. Теория теплопроводности / А. В. Лыков. – М.: Высш. шк., 1967. – 600 с. 6. К у д и н о в, В. А. Методы решения параболических и гиперболических уравнений теп- лопроводности / В. А. Кудинов, И. В. Кудинов. – М.: Книжный дом «Либроком», 2012. – 280 с. Представлена кафедрой теоретических основ теплотехники и гидромеханики Поступила 19.09.2012 60