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

## Abstract

Using additional boundary conditions and additional required function in integral method of heat-transfer we obtain approximate analytical solution of transient heat conduction problem for an infinite plate with asymmetric boundary conditions of the first kind. This solution has a simple form of trigonometric polynomial with coefficients exponentially stabilizing in time. With the increase in the count of terms of a polynomial the obtained solution is approaching the exact solution. The introduction of a time-dependent additional required function, setting in the one (point) of the boundary points, allows to reduce solving of differential equation in partial derivatives to integration of ordinary differential equation. The additional boundary conditions are found in the form that the required solution would implement the additional boundary conditions and that implementation would be equivalent to executing the original differential equation in boundary points. In this article it is noted that the execution of the original equation at the boundaries of the area only (via the implementation of the additional boundary conditions) leads to the execution of the original equation also inside that area. The absence of direct integration of the original equation on the spatial variable allows to apply this method to solving the nonlinear boundary value problems with variable initial conditions and variable physical properties of the environment, etc.

## 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 Об одном методе решения нестационарных задач теплопроводности . . . области. Показано, что выполнение уравнения на границах приводит к его выполнению и внутри рассматриваемой области с точностью, зависящей от числа приближений принятого решения.
×

### 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

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

1. Кудинов В. А., Стефанюк Е. В. Задачи теплопроводности на основе определения фронта температурного возмущения // Известия РАН. Энергетика, 2008. № 5. С. 141-157.
2. Лыков А. В. Методы решения нелинейных уравнений нестационарной теплопроводности // Известия АН СССР. Энергетика и транспорт, 1970. № 5. С. 109-150.
3. Goodman T. R Application of integral methods to transient nonlinear heat transfer // Advances in Heat Transfer, 1964. Т. 1. С. 51-122. doi: 10.1016/S0065-2717(08)70097-2.
4. Biot M. A. Variational Principles in Heat Transfer: A Unified Lagrangian Analysis of Dissipative Phenomena / Oxford Mathematical Monographs. Oxford: Clarendon Press, 1970. x+185 pp.
5. Вейник А. И. Приближенный расчет процессов теплопроводности. М., Л.: Госэнергоиздат, 1959. 184 с.
6. Швец М. Е. О приближенном решении некоторых задач гидродинамики пограничного слоя // ПММ, 1949. Т. 13, № 3. С. 257-266.
7. Тимошпольский В. И., Постольник Ю. С., Андрианов Д. Н. Теоретические основы теплофизики и термомеханики в металлургии. Минск: Белорусская наука, 2005. 560 с.
8. Глазунов Ю. Т. Вариационные методы. М., Ижевск: НИЦ “Регулярная и хаотическая динамика”; Институт компьютерных исследований, 2006. 470 с.
9. Беляев Н. М., Рядно А. А. Методы нестационарной теплопроводности. М.: Высш. шк., 1978. 328 с.
10. Федоров Ф. М. Граничный метод решения прикладных задач математической физики. Новосибирск: Наука, 2000. 220 с.
11. Кудряшев Л. И., Меньших Н. Л. Приближенные решения нелинейных задач теплопроводности. М.: Машиностроение, 1979. 232 с.
12. Кудинов В. А., Стефанюк Е. В. Аналитический метод решения задач теплопроводности на основе введения фронта температурного возмущения и дополнительных граничных условий // Инженерно-физический журнал, 2009. Т. 82, № 3. С. 540-558.
13. Стефанюк Е. В., Кудинов В. А. Получение приближенных аналитических решений при рассогласовании начальных и граничных условий в задачах теории теплопроводности // Изв. вузов. Матем., 2010. № 4. С. 63-71.
14. Кудинов В. А., Кудинов И. В., Скворцова М. П. Обобщенные функции и дополнительные граничные условия в задачах теплопроводности для многослойных тел // Ж. вычисл. матем. и матем. физ., 2015. Т. 55, № 4. С. 669-680. doi: 10.7868/ S0044466915040080.
15. Лыков А. В. Теория теплопроводности. М.: Высш. шк., 1967. 600 с.

Copyright (c) 2016 Samara State Technical University