A method for solving problems of heat transfer during the flow of fluids in a plane channel

Full Text

Введение. В теории теплопроводности известны методы, основанные на определении глубины термического слоя, называемые интегральными методами теплового баланса [1-6]. Их применение весьма эффективно при решении ряда краевых задач, получение аналитических решений которых классическими методами не представляется возможным. Особенностью данной группы методов является использование допущения о конечной скорости распространения теплоты, несмотря на то, что решению подлежит параболическое уравнение теплопроводности, в основе которого лежит бесконечная скорость перемещения температурного импульса. В случае применения данных методов процесс нагрева (охлаждения) формально разделяется на две стадии, первая из которых характеризуется постепенным перемещением фронта температурного возмущения от поверхности тела к его центру. Дополнительной искомой функцией в первой стадии является величина фронта температурного возмущения (глубина погретого слоя). Во второй стадии изменение температуры происходит по всей толщине тела. Здесь в рассмотрение вводится дополнительная искомая функция, характеризующая изменение температуры в центре тела. Многочисленные исследования [1-3] показали, что с увеличением числа приближений (n) первой стадии процесса время (F01 ) достижения фронтом температурного возмущения центра тела уменьшается и в пределе при n → ∞ F01 → ∞. Таким образом, противоречие, связанное с допущением конечного значения скорости распространения теплоты при решении параболического уравнения теплопроводности интегральным методом теплового баланса, снимается ввиду того, что получаемое решение подтверждает бесконечную скорость ее распространения. Отметим, что с увеличением числа приближений первой стадии процесса диапазон времени, в котором она определена, уменьшается, а диапазон второй стадии - возрастает. Следовательно, при большом числе приближений первой стадии процесса из полученного решения могут быть найдены температуры лишь для малых и сверхмалых значений временн´й и пространо ственной переменных. Роль второй стадии в определении температурного состояния при этом возрастает. В настоящей работе рассматривается метод получения приближенного аналитического решения, позволяющий избежать использования первой стадии процесса в интегральном методе теплового баланса. Постановка задачи. Последовательность применения предлагаемого метода рассмотрим на примере задачи теплообмена для жидкости, движущейся в плоскопараллельном канале (рис. 1) в следующей математической постановке: ∂ 2 t(ξ, η) ∂t(ξ, η) =a (η > 0; 0 < ξ < r0 ); ω(ξ) ∂η ∂ξ 2 (1) ∂t(0, η) t(ξ, 0) = t0 ; = 0; t(r0 , η) = tc , ∂ξ 110 Об одном методе решения задач теплообмена . . . где t - температура; ξ - поперечная координата; η - продольная координата; a - коэффициент температуропроводности; r0 - половина ширины плоского канала; t0 - начальная температура; tc - температура стенки; ω(ξ) = 2 = 1.5ωср (1 - ξ 2 /r0 ) - распределение скорости по координате ξ (0 ξ r0 ); ωср - средняя скорость. Рис. 1. Схема стабилизированного ламинарного течения жидкости в плоскопараллельном канале [Figure 1. Scheme of stabilization of a laminar flow in a plane-parallel channel] При решении задачи (1) приняты следующие допущения: 1) течение жидкости ламинарное стабилизированное, процесс теплообмена стационарный; 2) температура жидкости на входе в канал постоянна по сечению и равна t0 ; 3) температура внутренней поверхности стенки трубы во всем диапазоне продольной переменной 0 η < ∞ постоянна и равна tc , причем tc = t0 ; 4) теплопроводностью жидкости в продольном направлении, а также теплопроводностью стенки канала пренебрегаем; 5) внутренние источники тепла и диссипация энергии не учитываем [7, 8]. Введем следующие безразмерные переменные и параметры: Θ= t - t0 ; tc - t0 y= ξ ; r0 x= 8 1 η , 3 Pe h (2) где Θ - относительная избыточная температура; y - безразмерная поперечная координата; x - безразмерная продольная координата; h = 2r0 - ширина канала. С учетом обозначений (2) задача (1) принимает вид (рис. 2) (1 - y 2 ) ∂Θ(y, x) ∂ 2 Θ(y, x) = (x > 0; 0 < y < 1); ∂x ∂y 2 Θ(y, 0) = 0; ∂Θ(0, x) = 0; ∂y Θ(1, x) = 1. (3) (4) (5) (6) Решение задачи. Введем дополнительную искомую функцию q(x) = Θ(0, x), (7) 111 Е р е м и н А. В., К у д и н о в И. В., Ж у к о в В. В. Рис. 2. Схема теплообмена при течении жидкости в плоскопараллельном канале [Figure 2. Scheme of heat transfer of a fluid flow in a plane-parallel channel] представляющую изменение температуры в центре канала (рис. 2). Температура в точке y = 0 будет изменяться во времени тотчас же после воздействия теплового удара (граничного условия первого рода) на стенку канала (y = 1), что связано с бесконечной скоростью распространения теплоты, заложенной в параболическом уравнении (3). Следовательно, процесс изменения функции q(x) включает весь диапазон изменения переменной 0 x < ∞ и весь диапазон изменения искомой функции 0 Θ(0, x) 1. Так как температура в центре пластины является искомой величиной задачи (3)-(6), введение функции q(x) её не изменяет и рассматривается лишь с целью упрощения процесса получения решения данной задачи. Решение задачи (3)-(6) ищется в виде n Θ(y, x) = 1 + bk (q) cos k=1 rπy , 2 r = 2k - 1, (8) где bk (q) - неизвестные коэффициенты. Соотношение (8) удовлетворяет граничным условиям (5), (6). Для определения коэффициентов bk (q) используются соотношение (7) и дополнительные граничные условия, задаваемые в точках y = 0 и y = 1 и определяемые таким образом, чтобы их выполнение решением (8) было эквивалентно выполнению уравнения (3) в этих точках. Для получения дополнительных граничных условий применительно к точкам y = 0 и y = 1 продифференцируем граничные условия (5), (6) по переменной x: ∂ 2 Θ(0, x) = 0; ∂y∂x ∂Θ(1, x) = 0. ∂x 112 (9) (10) Об одном методе решения задач теплообмена . . . Продифференцируем уравнение (3) по переменной y: -2y ∂Θ(y, x) ∂ 2 Θ(y, x) ∂ 3 Θ(y, x) + (1 - y 2 ) = . ∂x ∂x∂y ∂y 3 (11) Сравнивая соотношения (3) и (9), а также (10) и (11), получаем следующие дополнительные граничные условия: ∂ 3 Θ(0, x) = 0; ∂y 3 ∂ 2 Θ(1, x) = 0. ∂y 2 (12) (13) В точке y = 1 уравнение (3) приводится к соотношению (13), которое решением (8) выполняется в любом приближении. Следовательно, уравнение (3) в этой точке всегда выполняется, в связи с чем какие-либо другие дополнительные граничные условия здесь больше не рассматриваются. Для получения следующих дополнительных граничных условий в точке y = 0 продифференцируем условие (12) по x, а уравнение (3) - трижды по переменной y: ∂ 4 Θ(0, x) = 0; ∂x∂y 3 ∂ 2 Θ(y, x) ∂ 3 Θ(y, x) ∂ 4 Θ(y, x) ∂ 5 Θ(y, x) -6 - 6y + (1 - y 2 ) = . ∂x∂y ∂x∂y 2 ∂x∂y 3 ∂y 5 (14) (15) Сравнивая (14) и (15), применительно к точке y = 0 с учетом (9) находим ∂ 5 Θ(0, x) = 0. ∂y 5 Отсюда можно записать общую формулу для дополнительных граничных условий в точке y = 0: ∂ i Θ(0, x) = 0, ∂y i i = 3, 5, 7, . . . . (16) Отметим, что благодаря принятой системе координатных функций условия (16) решением (17) выполняются в любом приближении. Для нахождения неизвестных коэффициентов bk (q) используется соотношение (7) и некоторые получаемые из него дополнительные условия. Для их определения продифференцируем соотношение (7) по переменной x: dq(x) ∂Θ(0, x) = . dx ∂x (17) Сравнивая (17) с уравнением (3), положив y = 0, находим дополнительное условие вида dq(x) ∂ 2 Θ(0, x) = . (18) dx ∂y 2 113 Е р е м и н А. В., К у д и н о в И. В., Ж у к о в В. В. Для получения следующего дополнительного условия продифференцируем соотношение (18) по переменной x: d2 q(x) ∂ 3 Θ(0, x) = . dx2 ∂x∂y 2 (19) Продифференцируем уравнение (3) дважды по переменной y: -2 ∂Θ(y, x) ∂ 2 Θ(y, x) ∂ 3 Θ(y, x) ∂ 4 Θ(y, x) - 4y + (1 - y 2 ) = . ∂x ∂x∂y ∂x∂y 2 ∂y 4 (20) Сравнивая (19) и (20) с учетом (17), находим еще одно дополнительное граничное условие: d2 q(x) dq(x) ∂ 4 Θ(0, x) -2 = . dx2 dx ∂y 4 Таким образом, дифференцируя каждое предыдущее дополнительное граничное условие по переменной x, а уравнение (3) - многократно - по переменной y и сравнивая получаемые соотношения, можно получить какое угодно число дополнительных граничных условий. Например, следующие три таких условия будут d2 q(x) ∂ 6 Θ(0, x) d3 q(x) - 14 = ; dx3 dx2 ∂y 6 d4 q(x) d3 q(x) d2 q(x) ∂ 8 Θ(0, x) - 44 + 60 = ; dx4 dx3 dx2 ∂y 8 d5 q(x) d4 q(x) d3 q(x) ∂ 10 Θ(0, x) - 100 + 8444 = . dx5 dx4 dx3 ∂y 10 (21) Для нахождения решения задачи (3)-(6) в первом приближении, подставляя (8) (ограничиваясь одним членом ряда) в соотношение (7), относительно неизвестного коэффициента b1 (q) получаем алгебраическое линейное уравнение, из решения которого находим b1 (q) = q(x) - 1. Подставляя найденное значение b1 (q) в (8), получаем Θ(y, x) = 1 + (q(x) - 1) cos(πy/2). (22) Потребуем, чтобы соотношение (21) удовлетворяло не уравнению (3), а некоторому осредненному уравнению - интегралу теплового баланса: 1 (1 - y 2 ) 0 ∂Θ(y, x) dy = ∂x 1 0 ∂ 2 Θ(y, x) dy. ∂y 2 (23) Подставляя (22) в (23), после определения интегралов относительно неизвестной функции q(x) получаем обыкновенное дифференциальное уравнение dq(x) π 4 π4 + q(x) - = 0. dx 32 32 114 (24) Об одном методе решения задач теплообмена . . . Интегрируя уравнение (24), находим q(x) = 1 + C exp -π 4 /32 x , (25) где C - постоянная интегрирования. Подставляя (25) в (22), получаем Θ(y, x) = 1 + C exp (-π 4 /32)x cos(πy/2). (26) Для определения постоянной интегрирования C составим невязку граничного условия (4) и потребуем ортогональности невязки к первой координатной функции: 1 (1 + C cos(πy/2)) cos(πy/2)dy = 0. (27) 0 Вычисляя интеграл (27), относительно константы интегрирования C получаем алгебраическое линейное уравнение. Его решение C = -4/π. После определения константы интегрирования решение задачи (3)-(6) в первом приближении находится из (26). Результаты расчетов по формуле (26) в сравнении с численным решением приведены на рис. 3. Из их анализа следует, что в диапазоне 0.2 x ∞ расхождение результатов не превышает 5 %. С уменьшением x расхождение возрастает. Для повышения точности решения необходимо увеличивать число членов ряда (8). При этом следует вводить дополнительные граничные условия. Так, во втором приближении (ограничиваясь двумя членами ряда (8)), для определения неизвестных коэффициентов b1 (q), b2 (q) используются соотношения (7), (18). Рис. 3. (online в цвете) Распределение температуры при x = 0.04 (чёрные кривые), x = = 0.2 (синие кривые), x = 0.4 (зелёные кривые), x = 0.8 (красные кривые); метки 1, 2, 3 - номера приближений; 4 - метод прогонки [Figure 3. (color online) Temperature distribution when x = 0.04 (black curves), x = 0.2 (blue curves), x = 0.4 (green curves), x = 0.8 (red curves); The labels 1, 2, 3 are approximate numbers of solutions for boundary-value problem (3)-(6); the label 4 is a solution by sweep method] 115 Е р е м и н А. В., К у д и н о в И. В., Ж у к о в В. В. Подставляя (8) в (7), (18), получаем  q(x) - b1 (x) - b2 (x) - 1 = 0;  2 2  dq(x) + π b (x) + 9π b (x) = 0.  1 2 dx 4 4 (28) Из решения системы уравнений (28) находим b1 (x) = 9 1 dq(x) 9 + q(x) - , 2π 2 dx 8 8 b2 (x) = - 1 dq(x) 1 1 - q(x) + . 2π 2 dx 8 8 (29) С учетом (29) соотношение (8) принимает вид Θ(y, x) = 1 + πy 1 dq(x) 9 9 cos - + q(x) - 2π 2 dx 8 8 2 1 dq(x) 1 3πy 1 - - 2 cos . (30) - q(x) + 2π dx 8 8 2 Подставляя (30) в (23), после определения интегралов относительно неизвестной функции q(x) получаем обыкновенное дифференциальное уравнение: 896 d2 q(x) dq(x) + (1952π 2 + 108π 4 ) + 81π 6 q(x) - 81π 6 = 0. 2 dx dx (31) Общее решение уравнения (31) имеет вид q(x) = 1 + C1 exp(-30.38233x) + C2 exp(-2.86059x), где C1 , C2 - постоянные интегрирования. Для определения постоянных интегрирования C1 , C2 составим невязку начального условия (4) и потребуем ортогональности невязки к первой и второй координатным функциям: 1 1 Θ(y, 0) cos(πy/2)dy = 0; 0 Θ(y, 0) cos(3πy/2)dy = 0. (32) 0 После подстановки (30) в (32) и определения интегралов соотношения (32) относительно C1 и C2 приводятся к системе линейных уравнений -0.20629C1 + 0.48812C2 + 0.63662 = 0; 0.70707C1 + 0.009959C2 - 0.21221 = 0, из решения которой получаем C1 = 0.31661, C2 = -1.17042. После определения неизвестной функции q(x) и констант интегрирования C1 , C2 решение задачи (3)-(6) во втором приближении находится из (30). Результаты расчетов безразмерной температуры по формуле (30) в сравнении с точным решением представлены на рис. 3. Их анализ позволяет заключить, что в диапазоне изменения продольной координаты 0.1 x < ∞ расхождение результатов не превышает 1 %. 116 Об одном методе решения задач теплообмена . . . На рис. 3 приведены также результаты расчетов температуры в третьем приближении (три члена ряда (8)). Из их анализа следует, что в диапазоне 0.1 x < ∞ полученное решение практически совпадает с численным. На рис. 4, 5 приведены графики изменения невязок уравнения (3) и граничного условия (4) в первом, втором и третьем приближениях. Их анализ позволяет заключить, что с увеличением числа приближений невязки уменьшаются, что свидетельствует о сходимости рассматриваемого метода решения краевой задачи (3)-(6). Рис. 4. Невязка уравнения (3) при x = 0.05: 1, 2, 3 - первое, второе и третье приближения [Figure 4. The residual error of the equation (3) when x = 0.05; The labels 1, 2, 3 are approximate numbers of solutions for boundary-value problem (3)-(6)] Рис. 5. Невязка граничного условия (4): 1, 2, 3 - первое, второе и третье приближения [Figure 5. The residual error of the boundary condition (4); The labels 1, 2, 3 are approximate numbers of solutions for boundary-value problem (3)-(6)] Заключение. Рассмотренный метод получения решения, основанный на использовании дополнительной искомой функции в интегральном методе теплового баланса, позволяет свести решение дифференциального уравнения в частных производных к интегрированию обыкновенного дифференциального уравнения. При этом требуется выполнение не исходного дифференциального уравнения, а этого же уравнения, осредненного по поперечной пространственной переменной. Такой метод позволяет находить решения краевых задач со сложными дифференциальными операторами в уравнениях (нелинейные, с переменными коэффициентами и др.), получение решений которых с помощью классических аналитических методов либо затруднительно, либо вообще 117 Е р е м и н А. В., К у д и н о в И. В., Ж у к о в В. В. не представляется возможным. Показано, что выполнение дифференциального уравнения на границах области с увеличением числа приближений приводит к его выполнению и внутри нее с точностью, зависящей от числа приближений.

About the authors

Anton V Eremin

Samara State Technical University

Email: a.v.eremin@list.ru
(Cand. Tech. Sci.; a.v.eremin@list.ru; Corresponding Author), Associate Professor, Dept. of Theoretical Fundamentals of Heat-Engineering and Hydromechanics 244, Molodogvardeyskaya st., Samara, 443100, Russian Federation

Igor V Kudinov

Samara State Technical University

Email: igor-kudinov@bk.ru
(Cand. Tech. Sci.; igor-kudinov@bk.ru), Associate Professor, Dept. of Theoretical Fundamentals of Heat-Engineering and Hydromechanics 244, Molodogvardeyskaya st., Samara, 443100, Russian Federation

Vitaliy V Zhukov

Samara State Technical University

Email: vrbatacom@mail.ru
Postgraduate Student, Dept. of Theoretical Fundamentals of Heat-Engineering and Hydromechanics 244, Molodogvardeyskaya st., Samara, 443100, Russian Federation

References