Известия РАН. Энергетика, 2019, № 6, стр. 152-158

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

Ю. В. Видин¹, В. С. Злобин^1, *

¹ Сибирский федеральный университет
Красноярск, Россия

Поступила в редакцию 03.08.2019
После доработки 22.11.2019
Принята к публикации 28.11.2019

DOI: 10.1134/S000233101906013X

Аннотация

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

Ключевые слова: температурное поле, коэффициент теплопроводности, характеристическое уравнение, собственные числа, функции Бесселя дробного порядка, аналитическое решение

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

(1)

$\frac{{\partial \vartheta }}{{\partial {\text{Fo}}}} = \frac{1}{\psi }\frac{\partial }{{\partial \psi }}\left[ {\psi \lambda \left( \psi \right)\frac{{\partial \vartheta }}{{\partial \psi }}} \right],\,\,\,\,{\text{0}} \leqslant {\text{Fo}} < \infty ;\,\,\,\,{{\psi }_{{\text{0}}}} \leqslant \psi \leqslant {\text{1}},$

(2)

$\frac{{\partial \vartheta }}{{\partial \psi }} = {\text{0}}\,\,\,\,{\text{при}}\,\,\,\,\psi = {{\psi }_{0}},$

(3)

$\frac{{\partial \vartheta }}{{\partial \psi }} = - {\text{Bi}}\vartheta \,\,\,\,{\text{при}}\,\,\,\,\psi = 1,$

(4)

$\vartheta \left( {\psi ,0} \right) = 1.$

Здесь использованы общепринятые обозначения [1], которые позволяют наиболее удобным аналитическим способом представить изучаемый тепловой процесс. С целью упрощения теоретических выкладок принято, что внутренняя поверхность трубы $\left( {\psi = {{\psi }_{0}}} \right)$ является адиабатической. Очевидно, что в том случае, когда ${{\psi }_{0}} = 0$, полый цилиндр преобразуется в сплошной.

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

(5)

$\lambda \left( \psi \right) = \frac{1}{\psi }$

на интервале изменения радиальной координаты $\psi $ в пределах ${{\psi }_{0}} \leqslant \psi \leqslant 1.$ Тогда исходное дифференциальное уравнение (1) принимает более простой вид:

(6)

$\frac{{\partial \vartheta }}{{\partial {\text{Fo}}}} = \frac{1}{\psi }\frac{{{{\partial }^{2}}\vartheta }}{{\partial {{\psi }^{2}}}}.$

Его аналитическое решение совместно с краевыми условиями (2)–(4) можно записать как бесконечный ряд

(7)

$\vartheta \left( {\psi ,{\text{Fo}}} \right) = \sum\limits_{n = 1}^\infty {{{A}_{n}}{{K}_{n}}\left( \psi \right)\exp \left( { - \mu _{n}^{2}{\text{Fo}}} \right)} ,$

где ${{K}_{n}}\left( \psi \right)$ – некоторые собственные функции поставленной задачи, а ${{\mu }_{n}}$ – ее собственные значения. Для их определения необходимо провести исследование соответствующей задачи Штурма–Лиувилля, которая, в данном случае, будет иметь вид

(8)

$\frac{{{{d}^{2}}K}}{{d{{\psi }^{2}}}} + {{\mu }^{2}}\psi K = 0,$

(9)

$\frac{{dK}}{{d\psi }} = 0\,\,\,\,{\text{при}}\,\,\,\,\psi = {{\psi }_{0}},$

(10)

$\frac{{dK}}{{d\psi }} = - {\text{Bi}}K\,\,\,\,{\text{при}}\,\,\,\,\psi = 1.$

Для нахождения интеграла уравнения (8) целесообразно использовать дробные функции Бесселя [3]. Тогда его решение может быть записано в следующем виде:

(11)

$K = \sqrt \psi \left[ {B{{J}_{{\frac{1}{3}}}}\left( {\frac{{2\mu }}{3}\sqrt {{{\psi }^{3}}} } \right) + {{J}_{{ - \frac{1}{3}}}}\left( {\frac{{2\mu }}{3}\sqrt {{{\psi }^{3}}} } \right)} \right],$

где ${{J}_{{\frac{1}{3}}}}\left( {\frac{{2\mu }}{3}\sqrt {{{\psi }^{3}}} } \right)$ и ${{J}_{{ - \frac{1}{3}}}}\left( {\frac{{2\mu }}{3}\sqrt {{{\psi }^{3}}} } \right)$ – функции Бесселя первого рода соответственно дробного порядка $\frac{1}{3}$ и $ - \frac{1}{3}.$ В справочнике [3] приведены таблицы значений этих функций в пределах аргумента $0 \leqslant \frac{{2\mu }}{3}\sqrt {{{\psi }^{3}}} \leqslant 10.0$ с относительно крупным шагом, равным 0.2. В более ранней работе [4] имеются более подробные таблицы этих же функций при изменении величины $\frac{{2\mu }}{3}\sqrt {{{\psi }^{3}}} $ от 0 до 15 через интервал 0.1. Причем в [4] указанны числовые значения названных функций с четырьмя цифрами после запятой.

Рекомендуемые специальные функции могут быть представлены бесконечными степенными рядами [3]

(12)

$\begin{gathered} {{J}_{{\frac{1}{3}}}}\left( {2Z} \right) = \frac{{{{Z}^{{\frac{1}{3}}}}}}{{\Gamma \left( {\frac{4}{3}} \right)}} - \frac{{{{Z}^{{\frac{7}{3}}}}}}{{\Gamma \left( {\frac{7}{3}} \right)}} + \frac{{{{Z}^{{\frac{{13}}{3}}}}}}{{\Gamma \left( {\frac{{13}}{3}} \right)}} - \frac{{{{Z}^{{\frac{{19}}{3}}}}}}{{\Gamma \left( {\frac{{19}}{3}} \right)}} + ... = \\ = \frac{{{{Z}^{{\frac{1}{3}}}}}}{{\Gamma \left( {\frac{4}{3}} \right)}}\left( {1 - \frac{3}{4}{{Z}^{2}} + \frac{9}{{56}}{{Z}^{4}} - \frac{9}{{560}}{{Z}^{6}} + ...} \right), \\ \end{gathered} $

(13)

$\begin{gathered} {{J}_{{ - \frac{1}{3}}}}\left( {2Z} \right) = \frac{{{{Z}^{{ - \frac{1}{3}}}}}}{{\Gamma \left( {\frac{2}{3}} \right)}} - \frac{{{{Z}^{{\frac{5}{3}}}}}}{{\Gamma \left( {\frac{5}{3}} \right)}} + \frac{{{{Z}^{{\frac{{11}}{3}}}}}}{{\Gamma \left( {\frac{8}{3}} \right)}} - \frac{{{{Z}^{{\frac{{17}}{3}}}}}}{{\Gamma \left( {\frac{{11}}{3}} \right)}} + ... = \\ = \frac{{{{Z}^{{ - \frac{1}{3}}}}}}{{\Gamma \left( {\frac{2}{3}} \right)}}\left( {1 - \frac{3}{2}{{Z}^{2}} + \frac{9}{{20}}{{Z}^{4}} - \frac{9}{{160}}{{Z}^{6}} + ...} \right), \\ \end{gathered} $

где

(14)

$Z = \frac{\mu }{3}\sqrt {{{\psi }^{3}}} .$

Известно, что функции $\Gamma \left( {\frac{2}{3}} \right),$ $\Gamma \left( {\frac{4}{3}} \right),$ $\Gamma \left( {\frac{5}{3}} \right),$ $\Gamma \left( {\frac{7}{3}} \right),$ $\Gamma \left( {\frac{8}{3}} \right),$ $\Gamma \left( {\frac{{10}}{3}} \right),$ $\Gamma \left( {\frac{{11}}{3}} \right),$ $\Gamma \left( {\frac{{13}}{3}} \right)$ и т.д. являются гамма-функциями [3], для которых справедлива следующая рекуррентная формула [3]:

(15)

$\Gamma \left( {X + 1} \right) = X\Gamma \left( X \right),$

причем $\Gamma \left( {\frac{2}{3}} \right)$ = 1.354118 и $\Gamma \left( {\frac{4}{3}} \right)$ = 0.892980. Ограниченные ряды (13) и (14) с инженерной точки зрения обладают приемлемой точностью, если аргумент $Z \leqslant 1.$

Подставляя в (12) и (13) зависимость (14), получим

(16)

${{J}_{{\frac{1}{3}}}}\left( {\frac{{2\mu }}{3}\sqrt {{{\psi }^{3}}} } \right) = \frac{{{{{\left( {\frac{\mu }{3}} \right)}}^{{\frac{1}{3}}}}\sqrt \psi }}{{\Gamma \left( {\frac{4}{3}} \right)}}\left( {1 - \frac{{{{\mu }^{2}}{{\psi }^{3}}}}{{12}} + \frac{{{{\mu }^{4}}{{\psi }^{6}}}}{{504}} - \frac{{{{\mu }^{6}}{{\psi }^{9}}}}{{45360}} + ...} \right),$

(17)

${{J}_{{ - \frac{1}{3}}}}\left( {\frac{{2\mu }}{3}\sqrt {{{\psi }^{3}}} } \right) = \frac{1}{{\Gamma \left( {\frac{2}{3}} \right)\sqrt[3]{{\frac{\mu }{3}}}\sqrt \psi }}\left( {1 - \frac{{{{\mu }^{2}}{{\psi }^{3}}}}{6} + \frac{{{{\mu }^{4}}{{\psi }^{6}}}}{{180}} - \frac{{{{\mu }^{6}}{{\psi }^{9}}}}{{12960}} + ...} \right).$

Таким образом, решение дифференциального уравнения (8) при умеренных значениях аргумента $\frac{{2\mu }}{3}\sqrt {{{\psi }^{3}}} $ можно записать в виде

(18)

$\begin{gathered} \frac{1}{{\Gamma \left( {\frac{2}{3}} \right)\sqrt[3]{{\frac{\mu }{3}}}}}\,K = C\left( {\psi - \frac{{{{\mu }^{2}}{{\psi }^{4}}}}{{12}} + \frac{{{{\mu }^{4}}{{\psi }^{7}}}}{{504}} - \frac{{{{\mu }^{6}}{{\psi }^{{10}}}}}{{45360}} + ...} \right) + \\ + \left( {1 - \frac{{{{\mu }^{2}}{{\psi }^{3}}}}{6} + \frac{{{{\mu }^{4}}{{\psi }^{6}}}}{{180}} - \frac{{{{\mu }^{6}}{{\psi }^{9}}}}{{12960}} + ...} \right), \\ \end{gathered} $

где коэффициент $C$ равен

(19)

$C = \frac{{\Gamma \left( {\frac{2}{3}} \right)\sqrt[3]{{{{{\left( {\frac{\mu }{3}} \right)}}^{2}}}}}}{{\Gamma \left( {\frac{4}{3}} \right)}}\,B.$

Объединяя коэффициенты ряда (7) ${{A}_{n}}$ с сомножителем $\frac{1}{{\Gamma \left( {\frac{2}{3}} \right)\sqrt[3]{{\frac{\mu }{3}}}}},$ окончательно для собственной функции $K$ получим

(20)

$K = C\left( {\psi - \frac{{{{\mu }^{2}}{{\psi }^{4}}}}{{12}} + \frac{{{{\mu }^{4}}{{\psi }^{7}}}}{{504}} - \frac{{{{\mu }^{6}}{{\psi }^{{10}}}}}{{45360}} + ...} \right) + \left( {1 - \frac{{{{\mu }^{2}}{{\psi }^{3}}}}{6} + \frac{{{{\mu }^{4}}{{\psi }^{6}}}}{{180}} - \frac{{{{\mu }^{6}}{{\psi }^{9}}}}{{12960}} + ...} \right).$

Постоянная $C$ находится из граничного условия задачи (9), на основе которого имеем

(21)

$C = \frac{{\frac{{{{\mu }^{2}}\psi _{0}^{2}}}{2} - \frac{{{{\mu }^{4}}\psi _{0}^{5}}}{{30}} + \frac{{{{\mu }^{6}}\psi _{0}^{8}}}{{1440}} - ...}}{{1 - \frac{{{{\mu }^{2}}\psi _{0}^{3}}}{3} + \frac{{{{\mu }^{4}}\psi _{0}^{6}}}{{72}} - \frac{{{{\mu }^{6}}\psi _{0}^{9}}}{{4536}} + ...}}.$

Далее, используя граничное условие третьего рода на внешней поверхности полого цилиндрического тела (10), получим характеристическое уравнение для определения собственных значений поставленной задачи ${{\mu }_{n}}$

(22)

$\begin{gathered} C\left( {1 - \frac{{{{\mu }^{2}}}}{3} + \frac{{{{\mu }^{4}}}}{{72}} - \frac{{{{\mu }^{6}}}}{{45360}} + ...} \right) - \frac{{{{\mu }^{2}}}}{2}\left( {1 - \frac{{{{\mu }^{2}}}}{{15}} + \frac{{{{\mu }^{4}}}}{{720}} - ...} \right) = \\ = - {\text{Bi}}\left[ {C\left( {1 - \frac{{{{\mu }^{2}}}}{{12}} + \frac{{{{\mu }^{4}}}}{{504}} - \frac{{{{\mu }^{6}}}}{{45360}} + ...} \right) + \left( {1 - \frac{{{{\mu }^{2}}}}{6} + \frac{{{{\mu }^{4}}}}{{180}} - \frac{{{{\mu }^{6}}}}{{12960}} + ...} \right)} \right]. \\ \end{gathered} $

Отсюда, в частном случае, а именно, когда ${\text{Bi}} \to \infty $, легко получить

(23)

$C\left( {1 - \frac{{{{\mu }^{2}}}}{{12}} + \frac{{{{\mu }^{4}}}}{{504}} - \frac{{{{\mu }^{6}}}}{{45360}} + ...} \right) + \left( {1 - \frac{{{{\mu }^{2}}}}{6} + \frac{{{{\mu }^{4}}}}{{180}} - \frac{{{{\mu }^{6}}}}{{12960}} + ...} \right) = 0.$

Очевидно, что при ${{\psi }_{0}} = 0$ (сплошное цилиндрическое тело) полученные расчетные формулы существенно упрощаются. Выражение (22) совместно с соотношением (21) позволяет с достаточной точностью определять первое собственное число ${{\mu }_{1}}.$ Для расчета последующих чисел ${{\mu }_{n}}$ необходимо значительно расширить степенные ряды в уравнении (22), что усложняет процедуры нахождения ${{\mu }_{n}}$ ($n > 1$). Поэтому для расширения вычислительных возможностей предлагаемых аналитических зависимостей целесообразно в том случае, когда имеет место условие $\left( {\frac{{2\mu }}{3}\sqrt {{{\psi }^{3}}} } \right) \geqslant 3,$ воспользоваться несколько другими математическими приближениями для дробных функций Бесселя ${{J}_{{\frac{1}{3}}}}\left( {\frac{{2\mu }}{3}\sqrt {{{\psi }^{3}}} } \right)$ и ${{J}_{{ - \frac{1}{3}}}}\left( {\frac{{2\mu }}{3}\sqrt {{{\psi }^{3}}} } \right)$ [3], а конкретно

(24)

${{J}_{{\frac{1}{3}}}}\left( {\frac{{2\mu }}{3}\sqrt {{{\psi }^{3}}} } \right) = \sqrt {\frac{3}{{\pi \mu \sqrt {{{\psi }^{3}}} }}} \cos \left( {\frac{{2\mu }}{3}\sqrt {{{\psi }^{3}}} - \frac{5}{{12}}\pi } \right),$

(25)

${{J}_{{ - \frac{1}{3}}}}\left( {\frac{{2\mu }}{3}\sqrt {{{\psi }^{3}}} } \right) = \sqrt {\frac{3}{{\pi \mu \sqrt {{{\psi }^{3}}} }}} \cos \left( {\frac{{2\mu }}{3}\sqrt {{{\psi }^{3}}} - \frac{1}{{12}}\pi } \right).$

Данные соотношения являются сравнительно простыми и обладают при указанном выше ограничении достаточной точностью. Например, если $\frac{{2\mu }}{3}\sqrt {{{\psi }^{3}}} = 4,$ то согласно (24) и (25) получим ${{J}_{{\frac{1}{3}}}}\left( 4 \right) = - 0.3591;$ ${{J}_{{ - \frac{1}{3}}}}\left( 4 \right) = - 0.3300,$ а табличные значения [3] будут соответственно равны ${{J}_{{\frac{1}{3}}}}\left( 4 \right) = - 0.3554;$ ${{J}_{{ - \frac{1}{3}}}}\left( 4 \right) = - 0.3331.$

Таким образом, для вычисления ${{\mu }_{n}}$ при $n > 1$ целесообразно использовать в качестве собственной функции рассматриваемой задачи математическое соотношение вида

(26)

$K\left( \psi \right) = \frac{1}{{\sqrt[4]{\psi }}}\left[ {C\,\cos \left( {\frac{{2\mu }}{3}\sqrt {{{\psi }^{3}}} - \frac{5}{{12}}\pi } \right) + \cos \left( {\frac{{2\mu }}{3}\sqrt {{{\psi }^{3}}} - \frac{1}{{12}}\pi } \right)} \right].$

Производя подстановку (26) в граничное условие (8), можно вывести характеристическое уравнение для определения собственных значений ${{\mu }_{n}}$($n = 2,\,3$ и т.д.). После проведения такой операции получим

(27)

$\begin{gathered} \mu \left[ {С\sin \left( {\frac{{2\mu }}{3} - \frac{5}{{12}}\pi } \right) + \sin \left( {\frac{{2\mu }}{3} - \frac{1}{{12}}\pi } \right)} \right] = \\ = \left( {{\text{Bi}} - \frac{{\text{1}}}{{\text{4}}}} \right)\left[ {C\,\cos \left( {\frac{{2\mu }}{3} - \frac{5}{{12}}\pi } \right) + \cos \left( {\frac{{2\mu }}{3} - \frac{1}{{12}}\pi } \right)} \right]. \\ \end{gathered} $

Для сплошного цилиндрического тела ($C = 0$) это соотношение существенно упрощается [5, 6]. Из зависимости (27) также следует, что в частном случае, а именно, когда ${\text{Bi}} = \frac{{\text{1}}}{{\text{4}}},$ эта формула также становится проще

(28)

$С\sin \left( {\frac{{2\mu }}{3} - \frac{5}{{12}}\pi } \right) + \sin \left( {\frac{{2\mu }}{3} - \frac{1}{{12}}\pi } \right) = 0.$

При использовании решений (27) и (28) для определения ${{\mu }_{n}}$ ($n = 2,\,3$ и далее) допустимо для комплекса $C$ ограничиться формулой (21), если геометрический параметр ${{\psi }_{0}}$ относительно невелик. Для нахождения коэффициентов ${{A}_{n}}$ решения (7) может быть использовано начальное условие (4), на основе которого будет справедливо выражение

(29)

${{A}_{n}} = \frac{{\int\limits_{{{\psi }_{0}}}^1 {\psi {{K}_{n}}\left( \psi \right)\,d\psi } }}{{\int\limits_{{{\psi }_{0}}}^1 {\psi K_{n}^{2}\left( \psi \right)\,d\psi } }}.$

В табл. 1 приведены значения нескольких первых корней характеристических уравнений (22) и (27) для некоторых величин параметра ${{\psi }_{0}}$ и чисел Био (${\text{Bi}}$).

Таблица 1.

Корни характеристических уравнений (22) и (27)

${\text{Bi}}$	${{\psi }_{0}} = 0$			${{\psi }_{0}} = 0.1$			${{\psi }_{0}} = 0.3$
${\text{Bi}}$	${{\mu }_{1}}$	${{\mu }_{2}}$	${{\mu }_{3}}$	${{\mu }_{1}}$	${{\mu }_{2}}$	${{\mu }_{3}}$	${{\mu }_{1}}$	${{\mu }_{2}}$	${{\mu }_{3}}$
0	0.0000	5.0306	9.7791	0.0000	5.1978	10.3158	0.0000	6.1785	11.3617
1.0	1.2848	5.3154	9.9306	1.2927	5.4913	10.4713	1.3597	6.4538	11.4985
5.0	2.1029	6.0978	10.4570	2.1185	6.3086	11.0162	2.2441	7.2535	11.9875
10.0	2.3463	6.5717	10.9113	2.3651	6.8165	11.4967	2.5160	7.7868	12.4367
50.0	2.6590	7.2444	11.8237	2.6861	7.5591	12.5019	2.9056	8.6155	13.4363
100.0	2.7204	7.3509	11.9942	2.7503	7.6789	12.6950	2.9958	8.7524	13.6350
$\infty $	2.7955	7.4613	12.1737	2.8298	7.8034	12.8993	3.1216	8.8947	13.8459

ВЫВОДЫ

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

Список литературы

Видин Ю.В., Злобин В.С., Иванов Д.И. Нестационарный теплоперенос в неоднородных конструкциях криволинейной конфигурации. Красноярск, СФУ, 2016. 167 с.
Видин Ю.В., Злобин В.С., Федяев А.А. Аналитический метод расчета нестационарного температурного поля при переменном коэффициенте теплопроводности // Системы. Методы. Технологии. 2019. БрГУ № 1(41). С. 57–60.
Янке Е., Эмде Ф., Лёш Ф. Специальные функции. М.: Наука, 1977. 342 с.
Динник А.Н. Таблицы бесселевых функций дробного порядка. Киев, АН УССР, 1933. 29 с.
Видин Ю.В., Злобин В.С., Федяев А.А. Нелинейные процессы переноса тепла в многослойных системах. Братск: Изд-во Братского государственного университета, 2019. 196 с.
Видин Ю.В., Казаков Р.В., Злобин В.С. К расчету собственных значений в задаче нестационарной теплопроводности массивного полого цилиндра // Изв. РАН. Энергетика. 2018. № 5. С. 124–130.

Дополнительные материалы отсутствуют.

Инструменты

предыдущая статья выпуска содержание выпуска

Известия РАН. Энергетика

Архивы выпусков Информация о журнале Отправить рукопись в журнал