On one method for solving transient heat conduction problems with asymmetric boundary conditions

Full Text

В теории теплопроводности известны методы, в которых используется понятие фронта температурного возмущения - глубины термического (прогретого) слоя [1-9]. То есть несмотря на то, что решению подлежит параболическое уравнение теплопроводности, в котором заложена бесконечная скорость распространения теплоты, вводится допущение о её конечности. При использовании данного метода процесс теплопроводности формально разделяется на две стадии по времени. Первая стадия заканчивается после достижения фронтом температурного возмущения противоположной стенки пластины. Во второй стадии изменение температуры происходит по всему объёму тела. Эта модель теплопроводности используется во многих методах: Т. Гудмена [3], М. Био [4], А. И. Вейника [5], М. Е. Швеца [6], Ю. С. Постольника [7] и др. Их важным преимуществом является возможность получения простых по форме приближённых аналитических решений многих сложных краевых задач - нелинейных, с переменными коэффициентами, динамических и тепловых задач пограничного слоя и др. Однако к их существенным недостаткам относится малая точность получаемых решений. С целью повышения точности требуется увеличивать число членов аппроксимационного ряда, для определения неизвестных коэффициентов которого используются различные методы. В частности, в работах [1, 10-14] используются дополнительные граничные условия, задаваемые в граничных точках и на фронте температурного возмущения, позволяющие существенно повысить точность получаемых решений. В работах [1, 12-14] показано, что с увеличением числа приближений n наряду с возрастанием точности решений происходит уменьшение времени продвижения фронта температурного возмущения Fo1 от поверхности до центра пластины (симметричная задача). И в пределе при n → ∞, Fo1 → 0, что подтверждает бесконечную скорость распространения теплоты. При этом было отмечено, что с увеличением числа приближений в связи с уменьшением времени первой стадии процесса её значение в оценке температурного состояния конструкции уменьшается, а второй - возрастает. В связи с этим в настоящей работе рассматривается метод решения, в котором исключается первая стадия процесса. Основную его идею рассмотрим на примере решения следующей краевой задачи (рис. 1): ∂ 2 Θ(ξ, Fo) ∂Θ(ξ, Fo) = (Fo > 0, 0 < ξ < 1); ∂Fo ∂ξ 2 Θ(ξ, 0) = 0; Θ(0, Fo) = 1; Θ(1, Fo) = 0; (1) (2) (3) (4) где Θ = (T - T0 )/(Tc1 - T0 ) - безразмерная температура; T0 - начальная температура; Tc1 - температура стенки при x = 0; Tc2 = T0 - температура стенки 343 К у д и н о в И. В., С т е ф а н ю к Е. В., С к в о р ц о в а М. П., К о т о в а Е. В., С и н я е в Г. М. Рис. 1. Расчетная схема для задачи теплообмена [Figure 1. A calculation scheme for a heat transfer problem] при x = δ; x - координата; δ - толщина пластины; ξ = x/δ - безразмерная координата; Fo = at/δ 2 - число Фурье (безразмерное время); t - время; a - коэффициент температуропроводности. В точке ξ = 1 введем дополнительную искомую функцию вида q(Fo) = ∂Θ(ξ, Fo) ∂ξ ξ=1 = tg α, (5) где α - зависящий от времени угол наклона температурной кривой к оси ξ. Минимальное значение этот угол (α = 0) будет иметь при Fo = 0, а максимальное (α = π/2) - при Fo → ∞. Ввиду бесконечной скорости распространения теплоты, описываемой параболическим уравнением теплопроводности (1), угол α будет увеличиваться сразу после приложения граничного условия первого рода в точке ξ = 0 и далее будет возрастать во времени, устремляясь к конечному значению α = π/2 при Fo → ∞, то есть при установлении стационарного режима теплообмена. Следовательно, процесс изменения производной от искомой функции по координате ξ в точке ξ = 1 будет включать весь диапазон времени нестационарного процесса 0 < Fo < ∞. Решение задачи (1)-(4) принимается в виде n Θ(ξ, Fo) = 1 + bk (q)ϕk (ξ), (6) k=1 где bk (q), k = 1, n, - неизвестные коэффициенты; ϕk (ξ) = cos r π ξ 1- 2 2 - координатные функции. Очевидно, что соотношение (6) благодаря принятой системе координатных функций удовлетворяет граничному условию (3). Для определения неизвестных коэффициентов bk будем использовать основное граничное условие (4), соотношение (5) и некоторые дополнительные граничные условия, задаваемые в точках ξ = 0 и ξ = 1 и определяемые таким образом, чтобы их 344 Об одном методе решения нестационарных задач теплопроводности . . . выполнение искомым решением (6) было эквивалентно выполнению уравнения (1) в этих точках [1, 12-14]. Найдём дополнительные граничные условия применительно к точке ξ = 0. Для этого продифференцируем граничное условие (3) по переменной Fo: ∂Θ(0, Fo) = 0. ∂Fo Сравнивания полученное соотношение с уравнением (1) применительно к точке ξ = 0, получаем дополнительное граничное условие вида ∂ 2 Θ(ξ, Fo) ∂ξ 2 ξ=0 = 0. (7) Продифференцируем соотношение (7) по переменной Fo и запишем полученное в виде ∂ 2 ∂Θ(ξ, Fo) = 0. (8) ∂ξ 2 ∂Fo ξ=0 Соотношение (8) с учётом уравнения (1) приводится к следующему дополнительному граничному условию: ∂ 4 Θ(ξ, Fo) ∂ξ 4 ξ=0 = 0. (9) Аналогично дифференцируя соотношение (9) по переменной Fo, с учётом уравнения (1) получаем ещё одно дополнительное граничное условие: ∂ 6 Θ(ξ, Fo) ∂ξ 6 ξ=0 = 0. (10) На основе соотношений (7), (9), (10) можно записать общую формулу дополнительных граничных условий в точке ξ = 0: ∂ i Θ(ξ, Fo) ∂ξ i ξ=0 = 0, i = 2, 4, 6, . . . . (11) Отметим, что благодаря принятой системе координатных функций все дополнительные граничные условия, определяемые формулой (11), решением (6) выполняются при любом числе приближений. Для каких-либо других, более сложных краевых задач (нелинейных, с переменными коэффициентами и др.), дополнительные граничные условия в точке ξ = 0 принятым решением могут и не выполняться заранее. Поэтому их выполнение следует проводить через определение неизвестных коэффициентов bk (q) искомого решения вида (6). Для нахождения дополнительных граничных условий в точке ξ = 1 продифференцируем соотношения (4), (5) по переменной Fo: ∂Θ(1, Fo) = 0; ∂Fo (12) 345 К у д и н о в И. В., С т е ф а н ю к Е. В., С к в о р ц о в а М. П., К о т о в а Е. В., С и н я е в Г. М. d q(Fo) ∂ 2 Θ(ξ, Fo) = d Fo ∂ξ∂Fo ξ=1 (13) . Сравнивая соотношение (12) с уравнением (1), находим первое дополнительное граничное условие вида ∂ 2 Θ(ξ, Fo) ∂ξ 2 ξ=1 (14) = 0. Перепишем соотношение (13) в виде dq(Fo) ∂ ∂Θ(ξ, Fo) = dFo ∂ξ ∂Fo ξ=1 . (15) Соотношение (15) с учётом уравнения (1) приводится к следующему дополнительному граничному условию: ∂ 3 Θ(ξ, Fo) ∂ξ 3 ξ=1 = d q(Fo) . d Fo (16) Продифференцируем соотношения (14), (16) по переменной Fo и запишем полученное в виде ∂ 2 ∂Θ(ξ, Fo) ∂ξ 2 ∂Fo 3 ∂Θ(ξ, Fo) ∂ 3 ∂ξ ∂Fo ξ=1 ξ=1 = 0; (17) = 0. (18) Соотношения (17), (18) с учётом уравнения (1) приводятся к следующим дополнительным граничным условиям: ∂ 4 Θ(ξ, Fo) = 0; ∂ξ 4 ξ=1 ∂ 5 Θ(ξ, Fo) d2 q(Fo) = . ∂ξ 5 d Fo2 ξ=1 (19) (20) Аналогично, дифференцируя соотношения (19) и (20) по переменной Fo, с учётом уравнения (1) получаем дополнительные граничные условия вида ∂ 6 Θ(ξ, Fo) = 0; ∂ξ 6 ξ=1 d3 q(Fo) ∂ 7 Θ(ξ, Fo) = . ∂ξ 7 dFo3 ξ=1 (21) (22) Анализируя соотношения (14), (19), (21), а также (16), (20), (22), можно записать общие формулы дополнительных граничных условий применительно к точке ξ = 1: ∂ i Θ(ξ, Fo) ∂ξ i 346 ξ=1 = 0, i = 2, 4, 6, . . . , (23) Об одном методе решения нестационарных задач теплопроводности . . . ∂ 2i+1 Θ(ξ, Fo) ∂ξ 2i+1 ξ=1 = di q(Fo) , di Fo i = 1, 2, 3, . . . . (24) Очевидно, что граничное условие (4) и условие (5), а также все дополнительные граничные условия, определяемые по общим формулам (23), (24), решением (6) не выполняются. Все эти условия будут использованы для нахождения неизвестных коэффициентов bk (q) решения (6). Для получения решения в первом приближении подставим (6) (ограничиваясь двумя членами суммы) в соотношения (4) и (5): 1 + b1 (q)ϕ1 (1) + b2 (q)ϕ2 (1) = 0; q1 (Fo) - b1 (q)∂ϕ1 (1)/∂ξ - b2 (q)∂ϕ2 (1)/∂ξ = 0. (25) Соотношения (25) относительно неизвестных коэффициентов b1 (q) и b2 (q) представляют собой систему двух алгебраических линейных уравнений, из решения которой находим √ √ 2(4q - 3π) 2(4q - π) b1 (q) = , b2 (q) = . 4π 4π Подставляя найденные значения b1 (q) и b2 (q) в (6), получаем √ 2 π ξ 3π ξ (4q - 3π) cos 1- + (4q - π) cos 1- Θ(ξ, Fo) = 1 + 4π 2 2 2 2 . (26) Потребуем, чтобы соотношение (26) удовлетворяло осредненному в пределах толщины пластины уравнению (1), то есть интегралу теплового баланса вида 1 2 1 ∂ Θ(ξ, Fo) ∂Θ(ξ, Fo) dξ = dξ. (27) ∂Fo ∂ξ 2 0 0 Подставляя (26) в (27), относительно неизвестной функции q(Fo) получаем следующее обыкновенное дифференциальное уравнение: dq(Fo) + A2 q(Fo) + A3 = 0, (28) dFo √ √ √ где A1 = 128 - 64 2, A2 = 12π 2 (2 + 2), A3 = 9π 3 2. Интегрируя уравнение (28), находим q(Fo) = C1 exp(νFo) - µ, (29) A1 где C1 - постоянная интегрирования, √ 3π 2 (2 + 2) √ ν= , 16( 2 - 2) √ 3π 2 µ= √ . 4( 2 + 2) Подставляя (29) в (26), получаем √ Θ(ξ, Fo) = 1 + √ 2 3 2 π ξ ν Fo (C1 e - µ) + cos 1- π 4 2 2 + 347 К у д и н о в И. В., С т е ф а н ю к Е. В., С к в о р ц о в а М. П., К о т о в а Е. В., С и н я е в Г. М. √ √ 2 2 3π ν Fo (C1 e - µ) + cos 1 - ξ/2 π 4 2 + . (30) Для определения постоянной интегрирования C1 составим невязку начального условия (2) и потребуем ортогональности невязки к координатной функции ϕ1 (ξ): 1 (31) Θ(ξ, 0)ϕ1 (ξ)dξ = 0. 0 Подставляя (30) в (31), получаем 1 0 √ √ 2 3 2 π ξ 1+ (C1 - µ) + cos 1- + π 4 2 2 √ √ 2 2 π ξ π ξ + (C1 - µ) + cos 1- + cos 1- π 4 2 2 2 2 dξ = 0. (32) Соотношение (32) ввиду ортогональности косинусов приводится к виду 1 0 √ π 3 2 ξ 2 (C1 - µ) + cos2 1- π 4 2 2 √ dξ = 1 =- cos 0 π ξ 1- 2 2 dξ. (33) После вычисления интегралов в (33) относительно неизвестного коэффициента C1 получаем алгебраическое уравнение, из решения которого находим √ C1 = π(3π(2 + 2) - 16)/16. После определения постоянной интегрирования решение задачи (1)-(4) в первом приближении находится из соотношения (30). Это решение точно удовлетворяет граничным условиям (3), (4) и интегралу теплового баланса (27). Уравнение (1) и начальное условие (2) в данном случае выполняются лишь приближенно в (первом приближении). Отметим, что число приближений будем определять не числом членов суммы решения (6), а порядком дифференциального уравнения относительно неизвестной функции q(Fo) или, что то же самое, числом констант интегрирования этого уравнения. Анализ результатов расчётов по формуле (30) позволяет заключить, что при 0.1 Fo < ∞ их расхождение с точным решением [15] не превышает 4 %. Повышение точности решения связано с увеличением числа членов ряда (6), неизвестные коэффициенты которого будем определять из условий (4), (5), (23), (24). Решение задачи (1)-(4) во втором приближении принимается в виде Θ(ξ, Fo) = 1 + b1 cos π ξ 1- 2 2 3π ξ 1- + 2 2 5π ξ 7π ξ + b3 cos 1- + b4 cos 1- 2 2 2 2 + b2 cos . (34) Подставляя (34) в условия (4), (5), (23) (при i = 2), (24) (при i = 1), относительно неизвестных коэффициентов bk будем иметь систему четырёх 348 Об одном методе решения нестационарных задач теплопроводности . . . алгебраических уравнений, из решения которой получаем b1 = B1 q + 29B2 q - 105B3 , b2 = B1 q + 37B2 q + 35B3 , (35) 11 1 13 1 b3 = B1 q - B2 q - 7B3 , b4 = B1 q + B2 q + 5B3 . 3 3 3 3 √ √ √ Здесь q = dq/dFo; B1 = 2/(2π 3 ); B2 = 2/(32π); B3 = 2/128. Подставляя (34) (с учётом (35)) в интеграл теплового баланса (27), относительно неизвестной функции q(Fo) будем иметь обыкновенное дифференциальное уравнение второго порядка (36) D1 q + D2 q - B3 q - D4 = 0, в котором √ √ D1 = 92/6 2 - 16384, D2 = π 2 (10624 2 - 34816), √ √ D3 = π 4 (6860 2 + 6720), D4 = 3675π 5 2, q = d2 q/dFo2 . Интегрируя уравнение (36), находим q(Fo) = C1 exp(z1 Fo) + C2 exp(z2 Fo) + D, (37) где z1,2 √ √ π 2 (2 22938 - 14368 2 ∓ 83 2 ± 272) √ , =± 16(9 2 - 16) √ 105π 2 √ D= ; 196 2 + 192 C1 , C2 - константы интегрирования, определяемые из начального условия (2). Для этого составляется невязка начального условия и требуется ортогональность невязки к координатным функциям ϕ1 (ξ) и ϕ2 (ξ): 1 Θ(ξ, 0)ϕj (ξ)dξ = 0, j = 1, 2. (38) 0 Подставляя (34) (с учётом (35), (37)) в (38), после определения интегралов относительно C1 и C2 получаем систему двух алгебраических линейных уравнений, из решения которой находим C1 = -1.490559; C2 = 2.022993. Результаты расчётов по формуле (34) в сравнении с точным решением [15] приведены на рис. 2. Их анализ позволяет заключить, что во втором приближении решение существенно уточняется. И, в частности, при 0.09 Fo < ∞ расхождение с точным решением уменьшается с 6 % (в первом приближении) до 2 % - во втором. Результаты расчётов в третьем и четвёртом приближениях в сравнении с точным решением приведены на рис. 2. Их анализ позволяет заключить, что при 0.04 Fo < ∞ расхождение четвёртого приближения с точным решением не превышает 1.5 %. Невязки уравнения (1) и начального условия (2) с увеличением числа приближений уменьшаются, что свидетельствует о сходимости предлагаемого метода решения. 349 К у д и н о в И. В., С т е ф а н ю к Е. В., С к в о р ц о в а М. П., К о т о в а Е. В., С и н я е в Г. М. Рис. 2. Графики изменения температуры: сплошная линия - точное решение; - второе приближение; - третье приближение; ◦ - четвёртое приближение; метки: 1 - Fo = 1.0; 2 - Fo = 0.2; 3 - Fo = 0.09; 4 - Fo = 0.04 [Figure 2. Graphs of the temperature change: the exact solution (solid curve), the second approximation of the solution ( ), the third approximation of the solution ( ), and the fourth approximation of the solution (◦ ); Labels: 1 - Fo = 1.0; 2 - Fo = 0.2; 3 - Fo = 0.09; 4 - Fo = 0.04] Таким образом, основным преимуществом рассмотренного метода является возможность сведения решения исходного дифференциального уравнения в частных производных к интегрированию обыкновенного дифференциального уравнения относительно дополнительной искомой функции, зависящей лишь от времени, что позволяет применять данный метод к решению сложных краевых задач (нелинейных, с переменными коэффициентами и др.), дифференциальные операторы которых не допускают разделения переменных, а также для задач с переменными во времени граничными условиями, с переменными начальными условиями, с источниками теплоты и др. Выводы. Получено приближённое аналитическое решение нестационарной задачи теплопроводности для бесконечной пластины с несимметричными граничными условиями первого рода. Введение дополнительной искомой функции q(Fo), представляющей изменение во времени угла наклона температурной кривой к оси пространственной переменной, позволяет свести решение уравнения в частных производных к интегрированию обыкновенного дифференциального уравнения. Использование функции q(Fo) обосновывается описываемой параболическим уравнением теплопроводности бесконечной скорости распространения теплоты, согласно которой угол наклона температурной кривой к оси ξ в точке ξ = 1 начинает возрастать сразу после приложения граничного условия первого рода в точке ξ = 0, и, следовательно, диапазон изменения функции q(Fo) включает весь диапазон времени нестационарного процесса 0 < Fo < ∞. Дополнительные граничные условия находятся в таком виде, чтобы их выполнение искомым решением было эквивалентно выполнению исходного дифференциального уравнения (в частных производных) в граничных точках 350 Об одном методе решения нестационарных задач теплопроводности . . . области. Показано, что выполнение уравнения на границах приводит к его выполнению и внутри рассматриваемой области с точностью, зависящей от числа приближений принятого решения.

About the authors

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

Ekaterina V Stefanyuk

Samara State Technical University

Email: stef-kate@yandex.ru
(Dr. Techn. Sci.; stef-kate@yandex.ru), Professor, Dept. of Theoretical Fundamentals of Heat-Engineering and Hydromechanics 244, Molodogvardeyskaya st., Samara, 443100, Russian Federation

Marina P Skvortsova

Samara State Technical University

Email: marina.dorozhkina.88@mail.ru
(Postgraduate Student; marina.dorozhkina.88@mail.ru), Assistant, Dept. of Theoretical Fundamentals of Heat-Engineering and Hydromechanics 244, Molodogvardeyskaya st., Samara, 443100, Russian Federation

Eugeniya V Kotova

Samara State Technical University

Email: larginaevgenya@mail.ru
(Cand. Techn. Sci.; larginaevgenya@mail.ru; Corresponding Author), Associate Professor, Dept. of Theoretical Fundamentals of Heat-Engineering and Hydromechanics 244, Molodogvardeyskaya st., Samara, 443100, Russian Federation

Gennady M Sinyaev

Samara State Technical University

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

References

Supplementary files