CC BY
200
Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки. — 2008. — № 2 (17). — С. 171—184
УДК 536.2(075)
АНАЛИТИЧЕСКИЕ РЕШЕНИЯ ЗАДАЧ ТЕПЛОПРОВОДНОСТИ ПРИ ПЕРЕМЕННЫХ ВО ВРЕМЕНИ КОЭФФИЦИЕНТАХ ТЕПЛООТДАЧИ
Е. В. Стефанюк, В. А. Кудинов
Самарский государственный технический университет, 443100, Самара, ул. Молодогвардейская, 244.
E-mails: [email protected]; [email protected]
С помощью интегрального метода теплового баланса на основе введения фронта температурного возмущения получено аналитическое решение задачи нестационарной теплопроводности при граничных условиях третьего рода с переменным во времени коэффициентом теплоотдачи. Построены графики распределения изотерм и скоростей их движения. Для повышения точности решения вводятся дополнительные граничные условия, определяемые из исходного дифференциального уравнения и основных граничных условий, включая условия, задаваемые на фронте температурного возмущения.
Ключевые слова: интегральные методы, аналитические решения, методы взвешенных невязок, фронт температурного возмущения, дополнительные граничные условия, изотермы, скорости движения изотерм.
При протекании любого естественного процесса теплового взаимодействия наряду с изменением температурного состояния тел происходит также изменение и граничных условий теплообмена. Особенно это характерно для коэффициента теплоотдачи. Переменность коэффициента теплоотдачи в процессах нагрева (охлаждения) тел в условиях конвекции является твёрдо установленным экспериментальным фактом [1—5], в связи с чем допущение о его постоянстве во многих практических расчётах может быть весьма грубым приближением к реальным физическим процессам теплообмена.
Известно, что решения задач теплопроводности, полученные с помощью классических аналитических методов, представляются в форме бесконечных рядов, плохо сходящихся в окрестностях граничных точек и при малых значениях временной координаты. Исследования, выполненные авторами настоящей работы, показывают, что сходимость точного аналитического решения нестационарной задачи теплопроводности для бесконечной пластины при граничных условиях первого рода в диапазоне чисел Фурье 10-12 ^ Fo ^ 10-7 наблюдается лишь при использовании от 1000 (Fo = 10-7) до пятисот тысяч (Fo = 10-12) членов ряда.
Эта проблема еще в большей степени характерна и для вариационных методов (Ритца, Треффтца, Л. В. Канторовича и др.), а также методов взвешенных невязок (ортогональный метод Бубнова—Галёркина, метод моментов, коллокаций и др.). Эти методы для получения решений нестационарных задач теплопроводности при малых значениях временной координаты
Стефанюк Екатерина Васильевна — ассистент кафедры теоретических основ теплотехники и гидромеханики Самарского государственного технического университета; к.т.н. Кудинов Василий Александрович — заведующий кафедрой теоретических основ теплотехники и гидромеханики Самарского государственного технического университета; д.ф.-м.н., профессор.
практически неприменимы ввиду того, что при большом числе приближений относительно неизвестных коэффициентов искомого решения получаются системы алгебраических линейных уравнений большой размерности. Матрицы коэффициентов таких систем, являясь заполненными квадратными матрицами с большим разбросом коэффициентов (по абсолютной величине), как правило, плохо обусловлены. В связи с этим с увеличением числа приближений точность решения может не улучшаться, а ухудшаться [6, 7].
К методам, позволяющим избежать указанных трудностей, относятся интегральные методы теплового баланса. Однако их широкое применение сдерживается недостаточной точностью получаемых решений. Всякие попытки увеличения точности не приводили к существенным результатам [8, 9].
Ниже будет изложен метод, относящийся к группе интегральных методов теплового баланса, позволяющий получать аналитические решения задач нестационарной теплопроводности с достаточно высокой точностью практически во всем диапазоне времени нестационарного процесса (0 ^ Ео < то) без каких-либо ограничений на величину числа Фурье в области малых его значений. Основную идею метода рассмотрим на примере решения задачи теплопроводности для бесконечной пластины при граничных условиях третьего рода в случае линейной зависимости коэффициента теплоотдачи от времени:
а(т) = ао(1 + /ЗТ),
где ао —начальное значение коэффициента теплоотдачи; в — коэффициент.
Математическая постановка задачи для бесконечной пластины при граничных условиях третьего рода имеет вид:
д 6(£, Ео) д2 6(£, Ео) _ А
—дЕо—^ = ——^ (Ео>0<^<1); (1) б(£, °) = °; (2)
^ = (3)
- Б1 (1 + Рё Ео) [6(0, Ео) - 1] =0, (4)
д Q(0, Fo)
где в = (Т — То)/(Тср — То) —относительная избыточная температура; £ = x/R — безразмерная координата; x — координата; R — половина толщины пластины; Fo = ат/R2 —число Фурье; а — коэффициент температуропроводности; т — время; То — начальная температура; Тср —температура среды; Bi = aoR/Л — число Био; Л — коэффициент теплопроводности; Pd = ¡3R2/a — критерий Предводителева.
Процесс нагрева разделим на две стадии по времени: 0 < Fo ^ Foi и Foi ^ Fo < то. Для этого введём движущуюся во времени границу (фронт температурного возмущения), разделяющую исходную область 0 ^ £ ^ 1 на две подобласти: 0 ^ £ ^ qi(Fo) и qi(Fo) ^ £ ^ 1, где qi(Fo) —функция, определяющая продвижение границы раздела во времени (рис. 1). При этом в области, расположенной за фронтом температурного возмущения, сохраняется начальная температура. Первая стадия процесса заканчивается при достижении движущейся границей центра пластины (£ = 1), т.е. когда Fo = Foi.
Во второй стадии изменение температуры происходит по всему объему тела 0 ^ £ ^ 1. Здесь вводится в рассмотрение дополнительная искомая функция 92 (Ро) = ©(1, Ро), характеризующая изменение температуры во времени в центре пластины.
Разделение процесса теплопроводности на два взаимосвязанных процесса позволяет существенно упростить последовательность получения решения задачи, т. к. в данном случае оказывается возможным применение метода аппроксимационного представления приближенного решения с определением любого числа его слагаемых. При определении неизвестных коэффициентов полинома возникает необходимость использования дополнительных граничных условий, нахождение которых связано с включением граничных точек по пространственной координате в область определения исходного дифференциального уравнения. Физический смысл дополнительных граничных условий состоит в том, что их выполнение равносильно выполнению исходного дифференциального уравнения в граничных точках и на фронте температурного возмущения (т. е. внутри области). При этом точность выполнения уравнения полностью определяется количеством дополнительных граничных условий, от которых, в свою очередь, зависит число членов аппроксимирующего полинома (число приближений).
Математическая постановка краевой задачи для первой стадии процесса в результате введения фронта температурного возмущения 91 (Ро) будет включать уравнение (1) с граничным условием (4), а также следующие граничные условия, выполняющиеся на фронте температурного возмущения:
Рис. 1. Расчётная схема теплообмена
6(91, РО) = 0;
д 6(91, Ро) д£
= 0, 0 ^ £ ^ 91 (РО),
(5)
(6)
где соотношения (5), (6) представляют условия сопряжения прогретой и непро-гретой зон. Соотношение (5) устанавливает равенство температуры тела в точке £ = 91 (Ро) его начальной температуре. Согласно условию (6) тепловой поток не распространяется за пределы фронта температурного возмущения (условие адиабатной стенки). Математическое доказательство условий (5), (6) дано в [10].
Обратим внимание на тот очевидный факт, что на первой стадии процесса задача (1), (4)—(6) за пределами фронта температурного возмущения, т.е. на отрезке 91 (Ро) ^ £ ^ 1, вообще не определена. В связи с этим здесь нет никакой необходимости выполнения начального условия вида (2) по всей толщине пластины (поэтому такое условие отсутствует в математической постановке задачи (1), (4)—(6)). В данном случае вполне достаточным является
выполнение граничного условия (5), согласно которому для всех £ = 51 (Ео) температура тела равна начальной температуре. Кроме того, в данной задаче отсутствует также граничное условие вида (3) ввиду того, что оно не влияет на процесс теплообмена в первой его стадии. Все это позволяет существенно упростить как процесс получения аналитического решения задачи (1), (4)— (6), так и окончательное выражение для него по сравнению с решением задачи (1)—(4) классическими аналитическими методами.
При использовании классического аналитического решения, полученного, например, с использованием метода разделения переменных, наибольшие проблемы возникают в случае нахождения температуры для малых значений числа Фурье (Ео ^ 0). Это связано с тем, что искомое решение в данном случае должно описывать две линии изменения температуры ЕО и ОС (см. рис. 1). Причем на линии ЕО температура зависит от координаты £, а на прямой ОС она остается постоянной и равной начальной температуре. Удовлетворить всем этим столь разнородным и к тому же изменяющимся во времени условиям в одном аналитическом выражении (при Ео ^ 0) можно лишь при использовании в нём бесконечно большого числа членов ряда, о чём уже упоминалось выше.
Простота предлагаемого здесь метода состоит в том, что получаемое решение в данном случае не связано с необходимостью аппроксимации температуры на участках ОС, ВС и др., т. е. за фронтом температурного возмущения 51 (Ео) ^ £ ^ 1, так как в этой области на первой стадии процесса задача (1), (4)—(6) не определена. Получаемое здесь аналитическое решение описывает изменение температуры во времени, характеризуемое лишь кривыми вида ЕО, АВ и др. (рис. 1), для практически точного описания которых достаточно всего нескольких членов ряда решения вне зависимости от величины временной координаты.
Отметим, что задача (1), (4)—(6) не относится к классу задач, в которых учитывается конечная скорость продвижения тепловой волны. Получение решений таких задач сводится к интегрированию гиперболического (волнового) уравнения теплопроводности [11]. Введённый в задаче (1), (4)—(6) фронт температурного возмущения по физическому смыслу является аналогом движущейся изотермы (но не тепловой волны). Ввиду того, что на фронте температурного возмущения в процессе его движения по координате поддерживается начальная температура 6(51, Ео) = 0, он является аналогом нулевой изотермы (см. ниже соотношение (39) и график на рис. 5).
Решение задачи (1), (4)—(6) разыскивается в виде следующего полинома:
где ак(51) —неизвестные коэффициенты. Для их определения используются граничные условия (4)- (6). Подставляя (7) (ограничиваясь тремя членами ряда) в (4)-(6), для определения ак(51) (к = 0,1, 2) будем иметь систему трёх алгебраических линейных уравнений. После определения ак(51) соотношение (7) примет вид
п
(7)
к=0
(8)
Для нахождения неизвестной функции qi(Fo) в первом приближении составим невязку уравнения (1) и проинтегрируем её в пределах глубины термического слоя (что равнозначно построению интеграла теплового баланса):
qi(Fo) qi(Fo)
« = / ^* (9)
0 0
Подставляя (8) в (9), после вычисления интегралов относительно функции qi(Fo) будем иметь следующее обыкновенное дифференциальное уравнение:
qiql [Bi qi(i + 2PdFo + Pd2Fo2) +4(PdFo + i)] +2Pdq?-
- 6Bi qi(l + 2PdFo + Pd2Fo2) - 12(l + PdFo) = 0. (10)
Начальное условие для уравнения (10) имеет вид
qi(0) = о. (11)
Ввиду нелинейности нахождение точного решения уравнения (10) весьма затруднительно. Поэтому будем искать его приближённое решение в виде
qi(Fo) = kFo\ (12)
где k, Л — неизвестные постоянные.
Очевидно, что соотношение (12) удовлетворяет начальному условию (11).
Для определения постоянных k и Л составим невязку уравнения (10) и запишем полученное соотношение для Fo = Foi. Подставляя (12) в (10), после некоторых преобразований получим (примем Bi = 1; Pd = 1):
k(1 +Foi)2(k2AFof - 6Fo(11+A)) +
+ 4(1 + Foi)(k2AFo2A - 3Fo(1+2A)) + 2k2Fo((1+2A) =0. (13)
Трансцендентное уравнение (13) содержит три неизвестных величины k, Fo1 и A. Для их определения привлечём уравнение вида
91(Fb1) = 1. (14)
Соотношение (12) с учётом (14) запишется так:
kFoA = 1. (15)
Выражая k из (15), получим
k = Fo-A. (16)
Подставляя (16) в (13), найдём
A(5 + 6Fo1 + Fo1) - 2Fo1 (8 + 12Fo1 + 3Fo1) = 0. (17)
175
Выразим Л из соотношения (17)
2Foi (8 + 12Foi + 3Fo?)
Л =
5 + 6Foi + Fof Подставляя Л в (15), будем иметь
k = Fo
2Fo^8+12Fo1+3Fol) 5+6Fo1+Fo1
i
Выражая из (14) Foi, найдём
Foi = k-1/x.
Подставляя (19) в (18), получим
_ 12+3fe1/A +8fe1/A k = 1+6fc2/A+6fc1/A
(18)
(19)
(20)
Задаваясь величиной Л, из соотношения (20) найдём такие значения к, при которых невязка уравнения (10) в диапазоне чисел 0 ^ Fo ^ Foi (где Foi находится из (19) ) наименее отклонялась бы от нуля. Графики изменения невязки при различных значениях Л даны на рис. 2. Их анализ позволяет заключить, что наименьшая невязка (е = 0,0009) наблюдается при величине Л = 0,511 (кривая 4). При этом к = 2,6079; Foi = 0,1532.
Соотношения (8), (12) с учётом найденных значений к и Л определяют решение задачи (1), (4)—(6) в первом приближении (при Bi = 1; Pd = 1).
На рис. 3 даны результаты расчёта перемещения фронта температурного возмущения по формуле (12) (при Bi = 1; Pd = 1). Результаты расчётов при Bi = 1 и Pd = 0 показали, что график qi(Fo) полностью совпадает с аналогичным графиком для случая граничных условий третьего рода с постоянным во времени коэффициентом теплоотдачи.
Рис. 2. Изменение невязки уравнения (10) при различных значениях Л, k, Foi: 1- Л = 0,505; k = 2,5936; Foi = 0,1515; 2- Л = 0,505; k = 2,6; Foi = 0,1515; 3- Л = = 0,51; k = 2,6055; Foi = 0,1529; 4 - Л = 0,511; k = 2,6079; Foi = 0,1532; 5 - Л = 0,515;
k = 2,6174; Foi = 0,1544
3
О 0,03 0,06 0,09 0,12 К}
Рис. 3. Перемещение фронта температурного возмущения по координате £ во времени Ро: 1-К = 1; Р(1 = то, 2-К = 1; Ра = 0, 3 —К = 1; Р(1 = 1
При Рё ^ то (а ^ то) теплообмен протекает при граничных условиях первого рода с температурой стенки, равной температуре среды. График (Ро) в данном случае совпадает с аналогичным графиком для граничных условий первого рода.
Результаты расчётов температуры по формуле (8) (Рё = 0; Б1 = 1) в сравнении с точным решением [11] даны на рис. 4. Их анализ позволяет сделать вывод, что отклонение температур, полученных по формуле (8), от точных их значений не превышает 2,5%.
Для повышения точности решения необходимо увеличивать число членов ряда (7), однако в случае, когда число коэффициентов а^ (51) больше трёх, для дополнительные граничные условия [12, продифференцируем граничное условие
0
0,4
0,2
0
0 ^0,2
0,4
0,6
Ро
Рис. 4. Изменение температуры в пластине при граничных условиях 3-го рода: сплошная линия — точное решение [11]; а — по формуле (8) (первая стадия); о — по формуле (31) (вторая стадия)
их определения необходимо вводить 13]. Для получения первого из них (4) по переменной Ро:
д2е(0'Ро) - Б.ае(0:Р0) - Б1Рё[в(0, Ро) - 1] - Б Рё РодМРО) =0. (21)
дРод£
дРо
д Ро
Продифференцируем уравнение (1) по переменной £ применительно к точке £ = 0:
д2е(0, Ро) = д 3е(0, Ро) д Род£ = д£з ' ( )
Подставляя (1), (22) в (21), получим первое дополнительное граничное условие вида
д3е(0, Ро) д2е(0, Ро)
д£3
- Б.
д£2
■ (1 + Рё Ро) - Б. Рё [е(0, Ро) - 1=0. (23)
Последовательность получения второго и третьего дополнительных граничных условий дана в [12, 13]. В данном случае они имеют вид
д26%Е0) = 0; (24)
de2
d 3Q(gi, Fo)
de3
0. (25)
Используя основные (4)—(6) и дополнительные (23)—(25) граничные условия, можно определить уже шесть коэффициентов ряда (7) и получить решение задачи (1), (4)—(6) во втором приближении для первой стадии процесса.
Во второй стадии процесса понятие фронта температурного возмущения теряет смысл и в качестве дополнительной искомой функции принимается температура в центре пластины 0(1, Ео) = 52 (Ео). Математическая постановка задачи в данном случае включает уравнение (1) с граничным условием (4), а также условия вида
6(1, Ео) = ^(Ео); (26)
96(1, Ео)
de
= 0 (Fo ^ Foi; 0 < e < 1). (27)
Начальным условием в задаче (1), (4), (26), (27) является решение задачи (1), (4)—(6) вида (8) в конце первой стадии процесса, т.е. при Fo = Foi. Соотношение (8) (ввиду того, что при Fo = Foi, qi(Foi) = 1) принимает вид
ч Bi(l+PdFbi)(l - 2e + С2) , ,
6(e, Foi) = —-(-iA )—^. (28)
u Bi(l+PdFoi) +2 v 7
Соотношение (28) является начальным условием задачи (1), (4), (26). Однако в его специальном выполнении нет необходимости. Это условие будет выполнено в процессе получения решения задачи (1), (4), (26), (27) ввиду того, что при Fo = Foi математические постановки задач (1), (4), (26), (27) и (1), (4)—(6) полностью совпадают, следовательно, совпадают для данного момента времени и их решения. Таким образом, в данном случае происходит плавный переход от первой стадии процесса ко второй.
Осредняя уравнение (1) в пределах всей толщины пластины (0 ^ e ^ 1), для второй стадии процесса получим следующий интеграл теплового баланса:
i i
de = / ^ de. (29)
0 0
Решение задачи (1), (4), (26), (27) разыскивается в виде ряда
n
6(e, Fo) = £ bk(q2)ek. (30)
k=0
Неизвестные коэффициенты Ък в первом приближении находятся из граничных условий (4), (26), (27). После определения Ък соотношение (30) примет вид
6(С Ео) = Б1 [1 + РёЕ)1 + (92 - 1) (1 + РёРо^ (2£ - С2)] + 2^2 (31) ( 0) 2 + Б1(1+РёРо) ' ( )
Подставляя (31) в (29), после вычисления интегралов относительно неизвестной функции 52 (Ро) получим следующее обыкновенное дифференциальное уравнение:
92 [6 + Б1(5 + Б1) + Б1 Рё Ро(5 + 2Б1 + Б1 Рё Ро)] -
- 92 [Б1(Рё - 6 - 3Б1) - 3Б1 Рё Ро(2Б1 + 2 + Б1 Рё Ро)] +
+ Б1(Рё - 6 - 3Б1) - 3Б1 Рё Ро(2Б1 + 2 + Б1 Рё Ро) = 0. (32)
Решение уравнения (32) разыскивается в виде суммы двух функций
92 (Ро) = П + Р, (33)
где п — частное решение неоднородного уравнения (32); р — общее решение соответствующего однородного уравнения, т. е. при
Б1(Рё - 6 - 3Б1) - 3Б1 Рё Ро(2Б1 + 2 + Б1 Рё Ро) = 0.
Частное решение неоднородного уравнения принимается в виде п = С, где константа С находится из выполнения уравнения (32). В данном случае С = 1.
Общее решение соответствующего однородного уравнения имеет вид
Г 2 + Б1 (1 + РёРо) 91п[3 + Б1 (1 + РёРо)] 1 Р = С еХР{1П3 + Б1(1+РёРо) + Б1РЧ-" " ■ (34)
где С1 — константа интегрирования, определяемая из начального условия 92(Ро1) = 0, задаётся формулой:
_ Г . 2 + Б1(1+РёРо1) 91п [3 + Б1(1+РёРо1)] \
С1 = - ехр < 3Ро1 - 1п-;---------—--— > ' (35)
1 \ 1 3 + Б1(1+РёРо1) Б1Рё / у 7
Соотношение (33) после определения частного решения неоднородного и общего решений соответствующего однородного уравнения принимает вид:
92 (Ро) = 1 - ехр
' 3 + Б1 (1 +РёРо1) , 2 + Б1 (1 + РёРо1) 1п----„ , „ „--1п
3 + Б1 (1 + РёРо) 2 + Б1 (1 + РёРо)
(36)
9 - 3 (Ро - Ро1)
Б1Рё 3 + Б1 (1 + РёРо) В частном случае при Рё = 0 значение 92(Ро) не определено.
Переходя к пределу при Рё ^ 0, соотношение (36) приводим к виду
3Б1Рс 3Б1РС1 ( 3Б1Ео\ 52(Ро) = ехр — - ехр ех^-—) . (37)
Непосредственной подстановкой можно убедиться, что соотношение (31) с учётом (36) точно удовлетворяет граничным условиям (4), (26), (27) и начальному условию (28), а соотношения (36), (37) удовлетворяют начальному условию 52(0) = 0. Уравнение (1), как это следует из (29), в данном случае удовлетворяется лишь в среднем.
Результаты расчётов по формуле (31) при Рё = 0, Б1 = 1 в сравнении с точным решением [11] представлены на рис. 4. Их анализ позволяет заключить, что в диапазоне чисел Фурье 0,25 ^ Ро < то полученное здесь решение практически совпадает с точным. Расхождение результатов расчётов при Ро = Р01 = 0,1532 составляет 2,5%.
Для увеличения точности решения необходимо увеличивать число членов ряда (30) и, следовательно, вводить дополнительные граничные условия [12, 13].
Первое дополнительное граничное условие в данном случае будет совпадать с аналогичным условием для первой стадии процесса, т. е. с соотношением (23). Второе и третье дополнительные граничные условия имеют вид:
д2 е(1, Ро) = ¿52 д3 е(1, Ро) = 0
д£ 2 = ¿Ро' д£ 3 =°- (38)
Основные граничные условия (4), (26), (27) и дополнительные (23), (38) позволяют определить уже шесть неизвестных коэффициентов (51) (к = = 0, 1, . . . , 5) ряда (30) и задать температурную функцию в виде полинома пятой степени.
Характерной особенностью полученных здесь аналитических решений является полиномиальная зависимость температуры от координаты £ в отличие от классических точных аналитических решений, где такая зависимость выражается через тригонометрические функции. Полиномиальная зависимость позволяет получить решение в виде поля изотермических линий. Отметим, что получение изотерм на основе классических аналитических решений рассмотрено в [14].
Для построения изотерм, используя формулы (8), (31), выразим координату £ как функцию температуры е(£, Ро) и времени. Выражая координату £ из формулы (8), получим
./251(1 + Ро) + 2е^1 [51 + 2 + Ро (251 + 2 + ^Ро)] £ = --2(РЪ+Т)-• (39)
Аналогично, выражая координату £ из (31), будем иметь
у/(3 + 4Ро + Ро2) (52 - 52е - 52 + е)
£ = ^-(1+Ро)(52 - 1)-• (40)
Соотношения (39), (40) позволяют для любых конкретных 6(£, Fo) = = const построить графики зависимости температур от £ и Fo (графики изотерм). Изотермические кривые, найденные по формулам (39), (40), приведены на рис. 5.
Первые производные по времени от соотношений (39), (40) позволяют определить безразмерные скорости движения изотерм и = d£/dFo по координате £ в зависимости от времени, а вторые производные — ускорений и = = d2£/dFo2. Формулы скоростей для первой и второй стадий процесса соответственно будут:
U1
0,000354k [10226k (Fo2'51 + 2F01'51 + Fo0'51)] Fo0'48^ ■ Fo0'5 [1 + Fo + 36kFo0'51 + 26 (1 + Fo)]
0,000354k [489Fo (1 + 26) + 10226]
u2 = — 2 2
Fo0'48 y kFo0'5 [1 + Fo + 36kFo0'51 + 26 (1 + Fo)] 1 2q2(92 - 1 - 6) + q2(Fo2 - Fo26 - 4Fo6 + 4Fo - 36 + 3) + 26
(1 + Fo)(92 - 1)^/(3 + 4Fo + Fo2)(q2 - 1)(q2 - 6)
(41)
(42)
Графики скоростей движения изотерм, построенных на основании формул (41), (42), даны на рис. 6.
Анализ распределения изотерм и скоростей их движения позволяет сделать вывод, что каждая изотерма появляется на поверхности пластины в строго определённый момент времени, имея при этом некоторую начальную скорость (см. рис. 5, 6). При этом для изотерм малого потенциала (0 ^ 6 < 0,1) наблюдаются высокие начальные скорости, устремляющиеся к бесконечности при 6 ^ 0. При увеличении 6 начальные скорости уменьшаются и, начиная с 6 = 0,4, незначительно отличаются друг от друга, стабилизируясь в пределах 0,7 ^ ио ^ 0,8. Скорости изотерм при их приближении к значению координаты £ = 1 устремляются к бесконечным величинам.
Первая
зона
Вторая зона
£ 0,8
0,4
е=:о/ |0,04|0Д /о,2
/0,3 /о,4 /0,5 /0,6
So,7
0,1 Foi 0,2
0,4
0,8
Fo
Рис. 5. Распределение изотерм в пластине при переменном во времени коэффициенте теплоотдачи (первое приближение, Foi = 0,1532; Bi = 1; Pd = 1)
Первая Вторая зона зона
Рис. 6. Распределение скоростей изотерм в пластине при переменном во времени коэффициенте теплоотдачи (первое приближение Foi = 0,1532; Bi = 1;
Pd = 1)
Отметим, что нулевая изотерма в(£, Fo) = 0 совпадает с графиком движения фронта температурного возмущения по координате £ в зависимости от времени Fo (см. рис. 3). В самом деле, при в(£, Fo) = 0 выражение (39) принимает вид соотношения (12). Отсюда следует, что по физическому смыслу фронт температурного возмущения является аналогом изотермы, движущейся по координате £ во времени. В данном случае это есть нулевая изотерма — изотерма начального условия.
Отметим, что графики ускорений изотерм по форме (качественно) практически соответствуют графикам скоростей и отличаются от них лишь количественно.
Ввиду невысокой точности решений в первом приближении изотермы, определяемые по формулам (39), (40), имеют незначительный излом при Fo = = Foi, т. е. в точке сопряжения решений для первой и второй стадий процесса (см. рис. 5). В связи с этим на графиках рис. 6 в этой точке имеет место некоторый скачок в эпюрах скоростей движения изотерм, который с увеличением точности решения не наблюдается, так же как и излом в изотермах.
Выводы
1. На основе введения фронта температурного возмущения с использованием дополнительных граничных условий разработана методика получения приближенных аналитических решений нестационарных задач теплопроводности во всем диапазоне изменения числа Фурье. Решения имеют вид степенных алгебраических полиномов, не содержащих специальных функций. Получение таких решений оказалось возможным благодаря разделению процесса теплопроводности на два взаимосвязанных процесса и введению соответственно этому двух дополнительных искомых функций — qi(Fo) для первой и q2(Fo) для второй стадии процесса. Введение функции qi(Fo), характеризующей перемещение фронта температурного возмущения по координате £ во времени, позволяет отказаться от выполнения начального условия вида ©(£, 0) = 0 по всей ширине пластины, заменив его начальным условием
(0) = 0, выполняемым лишь в точке £ = 0. Такой подход физически обосновывается тем, что на первой стадии процесса за пределами фронта температурного возмущения краевая задача не определена.
2. Дополнительные искомые функции 51(^0) и 52(^0) вводятся в полном соответствии с физическим смыслом краевой задачи. Обе они являются определяемыми величинами при получении решений любыми другими методами. Необходимость специального выделения этих функций связана с появляющейся при этом возможностью разделения исходной краевой задачи на два взаимосвязанных процесса — инерционный (нерегулярный) и установившийся (регулярный), закономерности изменения температуры во времени в которых настолько различны, что объединение их в рамках единой краевой задачи приводит к значительным трудностям получения ее аналитического решения.
3. Физический смысл применения дополнительных граничных условий заключается в возможности как можно более точного (в зависимости от числа граничных условий — числа приближений) выполнения исходного дифференциального уравнения внутри области и в ее граничных точках. Это свойство уже заложено в их выводе, основанном на требовании точного выполнения дифференциального уравнения и производных от него в граничных точках и на фронте температурного возмущения. Подчинение искомого решения дополнительным граничным условиям приводит к точному выполнению дифференциального уравнения в тех точках по координате £, в которых в данный момент времени находится фронт температурного возмущения. Так как область изменения фронта температурного возмущения охватывает весь диапазон изменения пространственной координаты 0 ^ 51(^0) ^ 1, то, следовательно, чем большее количество дополнительных граничных условий будет использовано, тем более точно будет выполняться уравнение внутри области.
4. Простота выражений для аналитических решений позволяет проводить исследования краевых задач в полях изотерм и определять скорости их перемещения во времени, что затруднительно выполнить, если использовать для этих целей классические точные аналитические решения.
При наличии экспериментальных данных по изменению температуры в пластине методика может быть использована для определения закономерности изменения коэффициента теплоотдачи во времени.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Шумаков, Н. В. Метод последовательных интервалов в теплометрии нестационарных процессов [Текст] / Н. В. Шумаков. —М.: Атомиздат, 1979. —212 с.
2. Тихонов, А. Н. Методы решения некорректных задач [Текст] / А.Н. Тихонов, В. Я. Ар-сенин. — М.: Наука, 1974.—288 с.
3. Тихонов, А.Н. Уравнения математической физики [Текст] / А.Н. Тихонов, А. А. Самарский.— М.-Л.: Гостехтеориздат, 1951. —724 с.
4. Смольский, Б. М. Нестационарный теплообмен [Текст] / Б. М. Смольский, Л. А. Сергеева, В. Л. Сергеев. — Мн.: Наука и техника, 1974. — 293 с.
5. Сергеев, В. Л. Влияние свойств калориметрического элемента на результаты измерения теплового потока [Текст] / В. Л. Сергеев, А. Г. Шашков, Л. А. Сергеева // Инженерно-физический ж. — 1968. — Т. 15, № 4. — С. 660-668.
6. Кудимов, В. А. Аналитические решения задач тепломассопереноса и термоупругости для многослойных конструкций [Текст]: Учебн. пособ. для вузов / В. А. Кудинов, Э. М. Карташов, В. В Калашников.—М.: Высш. шк., 2005.—430 с.
7. Кудинов, В. А. Метод координатных функций в нестационарных задачах теплопровод-
ности. Обзор [Текст] / В. А. Кудинов // Изв. АН. Энергетика. — 2004. — № 3. — С. 84-107.
8. Гудмен, Т. Применение интегральных методов в нелинейных задачах нестационарного теплообмена [Текст] / Т. Гудмен / В кн.: Проблемы теплообмена. — М.: Атомиздат, 1967.— С. 41-96.
9. Беляев, Н. М. Методы нестационарной теплопроводности [Текст]: Учебн. пособ. для вузов / Н.М. Беляев, А. А. Рядно. —М.: Высш. шк., 1978. — 328 с.
10. Карташов, Э. М. Аналитические методы в теории теплопроводности твердых тел [Текст] / Э. М. Карташов. — М.: Высш. шк., 2001. —550 с.
11. Лыков, А. В. Теория теплопроводности [Текст] / А. В. Лыков. —М.: Высш. шк., 1967. — 600 с.
12. Аналитические решения задач теплопроводности с переменным начальным условием на основе определения фронта температурного возмущения [Текст] / В. А. Кудинов, Б. В. Аверин, Е. В. Стефанюк и др. // Инженерно-физический ж. — 2007. — Т. 80, № 3. — С. 27-35.
13. Кудинов, В. А. Решения задач теплопроводности при переменных во времени граничных условиях на основе определения фронта температурного возмущения [Текст] / В. А. Кудинов, Б. В. Аверин, Е. В. Стефанюк // Извест. АН. Энергетика. — 2007. — № 1. — С. 55-68.
14. Цирельман, Н. М. Прямые и обратные задачи тепломассопереноса [Текст] / Н. М. Ци-рельман. — М.: Энергоатомиздат, 2005. — 392 с.
Поступила в редакцию 19/1У/2008; в окончательном варианте — 26/1Х/2008.
MSC: 80A17, 80M25
HEAT CONDUCTION PROBLEM ANALYTICAL SOLUTION AT TIME DEPENDENT HEAT TRANSFER COEFFICIENTS
E. V. Stefanjuk, V. A. Kudinov
Samara State Technical University,
443100, Samara, Molodogvardeyskaya str., 244.
E-mails: [email protected]; [email protected]
With the help of thermal balance integrated method and on the basis of introduction of temperature indignation front the analytical decision of a problem of non-stationary heat conductivity is obtained under boundary conditions of the third kind with variable in time factor of heat return. Graphs of isotherms distribution and velocity of their movement are constructed. In order to provide for the accuracy of the solution the additional boundary conditions defined from initial differential equation and basic boundary conditions, including the data obtained, at the front of temperature indignation are entered.
Key words: integral methods, analytical solutions, front of temperature perturbation, additional boundary conditions, isotherms, speeds of movement of isotherms.
Original article submitted 19/IV/2008; revision submitted 26/IX/2008.
Stefanjuk Ekaterina Vasil'evna, Ph.D. (Tech.), Dept. of the Theoretical Basics of Heattechnics and Hydromechanics of Samara State Technical University.
Kudinov Vasiliy Alexandrovich, Dr. Sci. (Phis. & Math.) Prof., Head of Dept. of the Theoretical Basics of Heattechnics and Hydromechanics of Samara State Technical University.
