Non-stationary heat exchange in cylindrical channel at laminar flow of fluids

Full Text

В процессе решения краевой задачи нестационарной теплопроводности при ламинарном течении жидкости в цилиндрическом канале происходит её физически обоснованное разделение на две задачи — нестационарную и стационарную, границы применимости которых определяются соотношениями между временной и продольной пространственной переменными. Дальнейшее решение каждой из полученных задач, определённых в строго фиксированных границах указанных переменных, выполняется раздельно. Нестационарная часть задачи представляет собой задачу теплопроводности для сплошного цилиндра (теплообмен протекает как бы в неподвижной жидкости), точные аналитические решения которой известны. Наибольшую трудность в данной случае представляет решение стационарной задачи с учётом движения жидкости (задача Гретца—Нуссельта). Эта задача впервые была решена Гретцем в 1885 г. [1]. Независимо от Гретца в 1910 г. её также решил Нуссельт [2]. Уточнённое решение задачи Гретца—Нуссельта дано в [3]. Отметим, что приведённое в [3] решение представляет собой бесконечный функциональный ряд, плохо сходящийся при малых значениях продольной координаты (x < 0,01), то есть в области, где исследования изменения температуры представляют наибольший интерес. Отметим также, что в зависимости от 122 Нестационарный теплообмен в цилиндрическом канале . . . величины поперечной координаты решение содержит функции Бесселя различного (в том числе и дробного) порядка. Такое решение малопригодно для инженерных приложений. В связи с этим разработка приближённых аналитических методов решения подобных задач имеет как научную ценность, так и практическое значение. Одно из таких решений путём совместного использования интегральных преобразований Лапласа и ортогонального метода Бубнова—Галёркина получено П. В. Цоем [4]. Решение найдено лишь в третьем приближении. Получение решения в последующих приближениях ограничивается проблемой, связанной с трудностями решения систем алгебраических линейных уравнений в общем виде, а также перехода от изображения искомой функции к её оригиналу. В настоящей работе излагается метод получения решения задачи Гретца— Нуссельта, основанный на совместном использовании метода Фурье и ортогонального метода Бубнова—Галёркина. Его отличительными особенностями являются простота получаемых аналитических выражений и возможность нахождения решения для малых значений продольной пространственной переменной [5–7]. Рассмотрим полную последовательность получения аналитического решения нестационарной задачи теплообмена при течении жидкости в цилиндрическом канале. Введём следующие допущения: течение жидкости и процесс теплообмена стационарны; жидкость несжимаема, её физические свойства постоянны; профиль скорости не изменяется по длине трубы; температура жидкости на входе в трубу неизменна по сечению и равна t0 ; темпера- Рис. 1. Схема стабилизированного ламинартура внутренней поверхности стенки ного течения жидкости в круглой трубе трубы постоянна и равна tc , причём tc = t0 ; начальная температура жидкости равна tн ; внутренние источники тепла и диссипацию энергии не учитываем; переносом теплоты за счёт теплопроводности вдоль оси трубы пренебрегаем (рис. 1). Математическая постановка задачи имеет следующий вид [3, 4]: ∂t(ξ, η, τ ) a ∂t(ξ, η, τ ) ∂ 2 t(ξ, η, τ ) ∂t(ξ, η, τ ) + ωx = +a , ∂τ ∂η ξ ∂ξ ∂ξ 2 τ > 0, η > 0, 0 < ξ < r; ∂t(0, η, τ ) t(ξ, η, 0) = tн ; t(ξ, 0, τ ) = t0 ; t(r, η, τ ) = tc ; = 0. ∂ξ (1) (2) Здесь t — температура; η, ξ — продольная и поперечная координаты; r — радиус трубы; t0 — температура жидкости на входе в трубу (при η = 0); tc — температура стенки; tн — начальная температура; a — коэффициент температуропроводности; ωx = 2ωср (1 − ξ 2 /r2 ) — распределение скорости по координате ξ; ωср = ωmax /2 — средняя скорость; ωmax — максимальная скорость ламинарного течения (рис. 1). 123 А. В. Е р е м и н, Е. В. С т е ф а н ю к, А. Ю. Р а с с ы п н о в, А. Э. К у з н е ц о в а Введём следующие безразмерные переменные и параметры: Θ= t − tc ; t0 − tc ξ y= ; r x= 1 η ; 2 Pe r Pe = ωcp r ; a Fo = aτ . r2 Здесь Θ — относительная избыточная температура; y — безразмерная поперечная координата; x — безразмерная продольная координата; Pe — число Пекле; Fo — число Фурье. С учётом принятых обозначений задача (1), (2) записывается так: ∂Θ(y, x, Fo) ∂Θ(y, x, Fo) 1 ∂Θ(y, η, Fo) ∂ 2 Θ(y, η, Fo) + (1 − y 2 ) = + , (3) ∂Fo ∂x y ∂y ∂y 2 Fo > 0, x > 0, 0 < y < 1; ∂Θ(0, x, Fo) Θ(y, x, 0) = Tн ; Θ(y, 0, Fo) = T0 = 1; Θ(1, x, Fo) = 0; = 0, (4) ∂y где Tн = (tн − tc )/(t0 − tc ). Для решения задачи (3), (4) применим двукратное преобразование Лапласа—Карсона по безразмерному времени Fo и относительной координате x [4]: ∞ T ∗ (y, S, Fo) = S Θ(y, x, Fo) exp(−Sx)dx; 0 ∞ T (y, S, P ) = P (5) T ∗ (y, S, Fo) exp(−P Fo)dFo, (6) 0 где P , S — параметры двукратного преобразования Лапласа—Карсона. Соотношение (6) с учётом (5) примет вид ∞ ∞ Θ(y, x, Fo) exp(−Sx − P Fo)dxdFo. T (y, S, P ) = P S 0 0 Задача (3), (4) в области изображений запишется следующим образом: P [T (y, S, P ) − Tн ] + S(1 − y 2 )[T (y, S, P ) − T0 ]− − T (1, S, P ) = 0; 1 ∂T (y, S, P ) ∂ 2 T (y, S, P ) − = 0; (7) y ∂y ∂y 2 ∂T (0, S, P ) = 0. ∂y (8) Решение задачи (7), (8) в соответствии с ортогональным методом Бубнова—Галёркина в первом приближении принимается в виде T (y, S, P ) = ν1 (S, P )ϕ1 (y), (9) где ν1 (S, P ) — неизвестный коэффициент-изображение; ϕ1 (y) — координатная функция, определяемая соотношением ϕ1 (y) = 1 − y 2 . 124 (10) Нестационарный теплообмен в цилиндрическом канале . . . Очевидно, что соотношение (9) при использовании координатной функции вида (10) удовлетворяет граничным условиям (8). Для определения неизвестного коэффициента ν1 (S, P ) составим невязку уравнения (7) и потребуем ортогональности невязки и координатной функции ϕ1 (y): 1 yP [T (y, S, P ) − Tн ] + yS(1 − y 2 )[T (y, S, P ) − T0 ]− 0 − ∂T (y, S, P ) ∂ 2 T (y, S, P ) −y ϕ1 (y)dy = 0. ∂y ∂y 2 Подставляя (10) в (9), преобразуем полученное выражение: 1 yP ν1 (S, P )(1 − y 2 ) − Tн + yS(1 − y 2 )[ν1 (S, P ) − T0 ]+ 0 + 2ν1 (S, P )y + 2ν1 (S, P )y (1 − y 2 )dy = 0. (11) Вычисляя интегралы в (11), относительно неизвестного коэффициента ν1 (S, P ) получаем P F1 + SF2 , ν1 (S, P ) = P F3 + SF4 + F5 где 1 1 y(1 − y 2 )dy, F1 = Tн 0 1 y(1 − y 2 )2 dy, F2 = T0 0 0 1 1 y(1 − y 2 )3 dy, F4 = 0 y(1 − y 2 )2 dy, F3 = y(1 − y 2 )dy. F5 = 4 0 В пространстве оригиналов соотношение (9) запишется в виде   F1 exp − F5 Fo ϕ1 (y) при x > F4 Fo,  F F3 F3 Θ(y, x, Fo) = F3 F5  2  exp − x ϕ1 (y) при x < F4 Fo. F3 F4 F3 (12) Верхняя строка формулы (12) совпадает с решением в первом приближении задачи (3), (4) при равенстве нулю конвективного члена (второе слагаемое в левой части) уравнения (3), то есть с решением нестационарной задачи теплопроводности для бесконечного цилиндра. Нижняя строка формулы (12) совпадает с решением в первом приближении стационарной задачи конвективного теплообмена, то есть при равенстве нулю первого члена в левой части уравнения (3). Таким образом, для областей теплообменника, которых не достигло возмущение, вызванное начальной температурой tн жидкости на входе в канал (при x = 0), нестационарный теплообмен происходит как бы в неподвижной жидкости, то есть перенос теплоты происходит только за счёт теплопроводности. Данная задача представляет собой нестационарную задачу теплопроводности для бесконечного цилиндра, классические точные аналитические решения которой известны. 125 А. В. Е р е м и н, Е. В. С т е ф а н ю к, А. Ю. Р а с с ы п н о в, А. Э. К у з н е ц о в а Для областей теплообменника, подверженных влиянию температурных условий на входе в трубу (жидкость при течении достигла этих областей), теплообмен не зависит от начальной температуры t0 . Теплообмен в этих областях не зависит от времени и полностью определяется течением среды, то есть задача становится стационарной с учётом конвективного переноса теплоты по оси x. Подобное разделение решения нестационарной задачи конвективного теплообмена при течении жидкостей в трубах и плоских каналах на стационарную и нестационарную составляющие физически обосновано и в других работах [3, 4]. В данном случае наибольшую трудность представляет нахождение решения стационарной задачи (задачи Гретца—Нуссельта). Её математическая постановка в данном случае имеет вид (1 − y 2 ) ∂Θ(y, x) ∂ 2 Θ(y, x) 1 ∂Θ(y, x) = + , x > 0, 0 ∂ ∂y 2 y ∂y ∂Θ(0, x) = 0; Θ(1, x) = 1. Θ(y, 0) = 0; ∂y y < 1; (13) (14) Следуя методу Фурье, решение задачи (13), (14) будем искать в виде Θ(y, x) = ϕ(x)ψ(y). (15) Подставляя (15) в (13), находим dϕ(x) + µϕ(x) = 0, dx d2 ψ(y) 1 dψ(y) + + µ(1 − y 2 )ψ(y) = 0, dy 2 y dy (16) (17) где µ — некоторая постоянная. Решение уравнения (16) известно и имеет вид ϕ(x) = A exp(−µx), (18) где A — неизвестный коэффициент. Удовлетворяя (15) двум последним условиям из (14), получаем ∂ψ(0) = 0, ∂y ψ(1) = 0. (19) Уравнение Бесселя (17) и граничные условия (19) представляют собой задачу Штурма—Лиувилля. Её точное аналитическое решение, полученное классическими методами, содержит функции Бесселя различного порядка (нулевого и 1/3), и к тому же решение выражается различными формулами для различных участков жидкости по радиальной координате y. В частности, рассматриваются три участка: вблизи оси трубы (для малых y), вблизи стенки (для больших y), для средних y [3]. Использование таких решений для инженерных приложений затруднительно, так как область определения 126 Нестационарный теплообмен в цилиндрическом канале . . . каждого из этих трёх решений по радиальной координате не задана, поэтому возникают трудности определения координат точек их сопряжения. Найдём приближённое аналитическое решение задачи (17), (19). Искомое решение принимается в виде n Bi y i+1 , ψ(y) = B0 + (20) i=1 где Bi — неизвестные коэффициенты, определяемые из основных граничных условий (19) и некоторых дополнительных условий. Далее будем полагать, что n = 48. Первое дополнительное условие следует из первого условия (19): ψ(0) = const := 1. (21) Подставляя (20) в (21), находим B0 = 1. В качестве других дополнительных условий будем использовать соотношения, получаемые путём удовлетворения решением (20) уравнения (17) в отдельных произвольно выбранных точках переменной y (они называются дополнительными условиями, так как выполняются не на границах, а в промежутке 0 < y < 1). Для получения дополнительных условий потребуем, чтобы соотношение (20) удовлетворяло уравнению (17) в следующих точках переменной y: y = 1/48, 2/48, . . . , 47/48. Ещё одно уравнение получим подстановкой (20) во второе граничное условие (19) (отметим, что первое граничное условие (19) соотношением (20) удовлетворяется), и, таким образом, всего будем иметь 49 алгебраических линейных уравнений относительно такого же количества неизвестных коэффициентов Bi . После определения коэффициентов Bi из решения системы, включающей 49 алгебраических линейных уравнений, и подстановки их в (20) это соотношение будет удовлетворять основным граничным условиям (19) и уравнению (17) в отмеченных выше 47 точках переменной y. Отметим, что в соотношении (20) неопределёнными оказываются собственные числа µk . Для их нахождения потребуем, чтобы соотношение (20) удовлетворяло уравнению (17) во всей области изменения переменной y (0 y 1). Для этого найдём интеграл взвешенной невязки уравнения (17): 1 0 d2 1+ dy 2 n Bi y i=1 i+1 1 d + 1+ y dy n Bi y i+1 + i=1 n 2 Bi y i+1 + µ(−y ) 1 + dy = 0. (22) i=1 После вычисления интегралов в (22) относительно собственных чисел µk получаем алгебраическое уравнение десятой степени, из решения которого находим µ1 = 7,3135869; µ2 = 44,609461; µ3 = 113,92103; µ4 = 215,24054; µ5 = 348,564115; µ6 = 513,8899; µ7 = 711,215; µ8 = 940,70; µ9 = 1202,9; µ10 = 1465,5. Точные значения собственных чисел [3]: µ1 = 7,3135868; µ2 = = 44,609460; µ3 = 113,92104; µ4 = 215,24054; µ5 = 348,56412; µ6 = 513,89; µ7 = 711,217; µ8 = 940,54; µ9 = 1201,8; µ10 = 1495,2. 127 А. В. Е р е м и н, Е. В. С т е ф а н ю к, А. Ю. Р а с с ы п н о в, А. Э. К у з н е ц о в а Высокая точность определения собственных чисел объясняется тем, что соотношение (20) удовлетворяет дифференциальному уравнению (17) в 47 точках переменной y. Собственные функции по найденным значениям собственных чисел находятся из (20). Подставляя (18), (20) в (7) находим n Bi y i+1 . Θ(y, x) = A exp(−µx) B0 + (23) i=1 Каждое частное решение (23) удовлетворяет двум последним граничным условиям (14) и уравнению (13) во всех 47 точках пространственной переменной y, в которых было выполнено уравнение (17). Однако ни одно из этих частных решений не удовлетворяет первому условию (14). Для его выполнения составим сумму частных решений n m Bi y i+1 , Ak exp(−µk x) B0 + Θ(y, x) = (24) i=1 k=1 где m = 10 — количество собственных чисел (десятое приближение). Неизвестные коэффициенты Ak найдём так, чтобы соотношение (24) удовлетворяло первому условию (14) в десяти точках (по числу собственных чисел) пространственной переменной y: y1 = 0, y2 = 1/10, y3 = 2/10, . . ., y10 = 9/10. Подставляя (24) в первое условие (14) и записывая полученное соотношение в десяти точках переменной y, относительно неизвестных коэффициентов Ak получаем систему алгебраических линейных уравнений 10 48 i+1 Bi yj Ak exp(−µk x) B0 + k=1 = 0, j = 1, 2, . . . , 10. (25) i=1 Из решения системы уравнений (25) находим A1 = 1,47205; A2 = −0,82047; A3 = 0,56376; A4 = −0,53858; A5 = 0,28605; A6 = −0,47515; A7 = 0,79096; A8 = 0,91147; A9 = 0,49147; A10 = −1,68158. После определения неизвестных коэффициентов Ak решение задачи (13), (14) в замкнутом виде находится из (24). Это решение в 47 точках пространственной переменной y удовлетворяет уравнению (13), в десяти точках переменной y — первому граничному условию из (14) и также удовлетворяет двум последним граничным условиям (14). На рис. 2 сплошными линиями приведены результаты расчётов по формуле (24) при m = 10, штриховыми линиями — точное решение [3], пунктирными линиями — результаты расчётов по методу конечных разностей с использованием метода прогонки. Анализ представленных результатов позволяет сделать вывод, что полученные по формуле (24) значения температур для x 0,01 практически совпадают с точными, а при x = 0,005 расхождение между полученным и точным решением не превышает 2 %. Выводы. Получено приближённое аналитическое решение задачи теплообмена при течении жидкости в круглой трубе, позволяющее для инженерных приложений с достаточной точностью определять температурное состояние 128 Нестационарный теплообмен в цилиндрическом канале . . . Рис. 2. Распределение температуры в цилиндрическом канале: сплошная линия — расчёт по формуле (24); штриховая линия — точное решение; пунктирная линия — метод конечных разностей; 1) x = 0,005, 2) x = 0,01, 3) x = 0,025, 4) x = 0,05, 5) x = 0,1, 6) x = 0,2, 7) x = 0,3 жидкости от продольной пространственной переменной x 0,005. В отличие от точного решения, полученное решение не содержит специальных функций и имеет вид алгебраических полиномов с коэффициентами, экспоненциально стабилизирующимися во времени. Получение решения для достаточно малых значений продольной переменной x оказалось возможным благодаря использованию дополнительных условий. Для получения этих условий используется искомое в виде полиномов решение дифференциального уравнения (17) задачи Штурма—Лиувилля, которое должно удовлетворять уравнению в заданном числе точек поперечной пространственной переменной. Использование этих условий также приводит к удовлетворению искомым решением исходного дифференциального уравнения (13) краевой задачи Гретца—Нуссельта в указанных точках.

About the authors

Anton V Eremin

Samara State Technical University

Email: a.v.eremin@list.ru
Postgraduate Student, Dept. of Theoretical Basis of Heat Engineering & Flow Mechanics 244, Molodogvardeyskaya st., Samara, 443100, Russia

Ekaterina V Stefanyuk

Samara State Technical University

Email: stef-kate@yandex.ru
(Dr. Techn. Sci.), Associate Professor, Dept. of Theoretical Basis of Heat Engineering & Flow Mechanics 244, Molodogvardeyskaya st., Samara, 443100, Russia

Alexander Yu Rassypnov

Samara State Technical University

Postgraduate Student, Dept. of Theoretical Basis of Heat Engineering & Flow Mechanics 244, Molodogvardeyskaya st., Samara, 443100, Russia

Anastasiya E Kuznetsova

Samara State Technical University

Email: kuznetsovaae@rambler.ru
Postgraduate Student, Dept. of Theoretical Basis of Heat Engineering & Flow Mechanics 244, Molodogvardeyskaya st., Samara, 443100, Russia

References

Supplementary files