Russian Federation
Russian Federation
Russian Federation
Samara, Russian Federation
Russian Federation
Russian Federation
Using additional boundary conditions (ABC) and an additional sought function (ASF), a solution to the heat conduction problem with non-homogeneous boundary conditions has been obtained. ABC allows satisfying the differential equation at the boundaries, leading to its fulfillment both inside the domain and without resorting to integration over Cartesian coordinates. ASF reduces the original partial differential equation to a temporal ordinary equation, from which the eigenvalues of the boundary value problem are determined, as formulated in the classical methods of the Sturm-Liouville problem, stated in spatial variables. Thus, this work considers an alternative way of determining eigenvalues. The integration constants are found from the initial condition using the least squares method, which allows their values to be determined with a given accuracy. The solution obtained based on ABC and ASF approximates n → ∞ the exact analytical solution in the form of an infinite series, including trigonometric coordinate functions with coefficients that stabilize exponentially in time. In this case, the eigenvalues determined from the solution of the temporal ordinary differential equation regarding the additional sought function coincide with their exact values at any approximation. The accuracy of the integration constants, determined by the method of least squares, is controlled by the number of approximation points in the range of the auxiliary variable's variation. It should be noted that the additional boundary conditions considered in this work hold true for any other method of obtaining a solution to the problem under consideration, including the exact method, as can be verified by direct substitution. Therefore, their introduction does not distort the original mathematical formulation of the problem but only significantly simplifies the process of obtaining its analytical solution.
thermal conductivity, boundary conditions, additional sought function, additional boundary conditions, analytical solution, least squares method
Введение
В теории теплопроводности применяются методы, связанные с определением интеграла взвешенной невязки исходного дифференциального уравнения, что приводит к его осреднению по пространственной переменной (Л. В. Канторовича, Бубнова – Галеркина, многие разновидности интегрального метода теплового баланса – Т. Гудмена, А. П. Ваничева, М. Био и др.) [1–7]. Все они относятся к группе приближенных аналитических методов, основной проблемой которых является низкая точность получаемых решений. Причина в том, что осреднение дифференциального уравнения по пространственной переменной приводит к низкой точности нахождения собственных чисел, определяющих выполнение дифференциального уравнения. Следовательно, повышение точности связано с увеличением точности выполнения уравнения в области изменения временной и пространственной координат. Для достижения этих целей в настоящей работе используются дополнительная искомая функция (ДИФ) и дополнительные граничные условия (ДГУ) [5, 7]. Их применение не влияет на математическую постановку исходной задачи и является лишь вспомогательным средством упрощения получения ее решения. Следует отметить, что все эти условия выполняются и при любом другом методе получения аналитических решений (включая и классические методы), несмотря на отсутствие их отдельного рассмотрения в процессах получения решений.
Математическая постановка задачи
Найдем решение следующей задачи
где c – теплоемкость пластины, Дж/(кг·°С); ρ – плотность, кг/м3; T - температура °С ; x – координата, м; τ – время, с; λ – теплопроводность, Вт/(м·К); δ − толщина пластины, м; T0 − начальная температура, К; , − температура пластины, К, при x = 0 и x = δ.
Обозначим
(5)
где Θ, Fo, ξ – безразмерные температура, время (число Фурье), координата соответственно; a = λ / (сρ) – коэффициент температуропроводности тела.
С учетом (5) задача (1)–(4) принимает вид (рис. 1)
Рис. 1. Схема теплообмена: α – угол наклона
температурной кривой к оси ξ в точке ξ = 1
Fig. 1. Heat exchange diagram: α is the inclination angle
of the temperature curve to the axis ξ at the point ξ = 1
Введем ДИФ q(Fo), характеризующую изменение градиента температуры в точке ξ = 1:
где α(Fo) – угол наклона температурной кривой к оси ξ в точке ξ = 1. Очевидно, что при Fo = 0 tgα = 0, ∂Θ(1, Fo) / ∂ξ = 0. При Fo →∞ tgα → –1, ∂Θ(1, Fo) / ∂ξ = –1.
Реализация метода
Рассмотрим получение решения задачи (6)–(9), приняв его в следующем виде:
, (11)
где bk(q) – неизвестные функции; φk(ξ) = sin(kπξ); 
.
Соотношение (11) удовлетворяет условиям (8), (9). Для нахождения коэффициентов bk(q) используется соотношение (10) и ДГУ, определяемые в граничных точках [5–9]. В работах [7, 8] математически доказывается, что выполнение уравнения в частных производных на границах влечет его выполнение и внутри области, исключая его непосредственное интегрирование, по пространственной переменной.
Найдем ДГУ, которые должны быть выполнены решением (11) в точке ξ = 0 [5, 7]. Дифференцируя условие (8) по Fo, находим
. (12)
Сопоставляя (12) с (6), получаем ДГУ вида
. (13)
Продифференцируем (13) по Fo:
. (14)
Сравнивая (14) и (6), получаем еще одно ДГУ в точке ξ = 0:
.
Аналогично могут быть получены и последующие ДГУ в точке ξ = 0. Общая формула для них имеет вид
(15)
Для получения ДГУ в точке ξ = 1 продифференцируем (9), (10) по Fo:
(16)
(17)
Выражения (16), (17) с учетом (6) приводятся к ДГУ
Продифференцируем (18), (19) по Fo:
Соотношения (20), (21), учитывая уравнение (6), приводятся к следующим:
Аналогично можно получить и другие ДГУ в точке ξ = 1. Общие формулы для них
Дополнительные граничные условия, определяемые формулами (15), (22), принятым решением (11) выполняются, что эквивалентно выполнению уравнения (6) в точках
ξ = 0 и ξ = 1 (в предельном смысле – правая и левая его части равны нулю [5, 8]). Следовательно, невыполненными соотношением (11) являются условия (10), (23), которые и будут использованы далее при определении функций bk(q).
В первом приближении, подставляя (11) в (10), находим: bk(q) = (1 + q(Fo)) / π.
С учетом b1(q) соотношение (11) будет
. (24)
В граничных точках ξ = 0 и ξ = 1 соотношение (24) удовлетворяет уравнению (6) в предельном смысле, что не позволяет найти неизвестную функцию q(Fo) из выполнения уравнения (6) в этих точках, в связи с чем потребуем, чтобы соотношение (24) удовлетворяло уравнению (6) в любой из точек переменной (0 < ξ < 1), исключая точки ξ = 0 и ξ = 1. Подставляя (24) в (6), получаем
. (25)
Интеграл уравнения (25) будет
, (26)
где C1 – постоянная интегрирования; v1 = π2 – первое собственное число, равное его точному значению [9, 10].
Подставим (26) в (24):
. (27)
Соотношение (27), независимо от величины C1, точно удовлетворяет граничным условиям (8), (9) и уравнению (6) во всем диапазоне пространственной переменной (0 ≤ ξ ≤ 1), включая точки ξ = 0 и ξ = 1, где оно удовлетворяется в предельном смысле, причем уравнение (6) выполнено, минуя процесс его непосредственного интегрирования по переменной ξ, что можно объяснить выполнением соотношением (27) ДГУ (15), (22). Таким образом, выполнение уравнения внутри области связано с использованием ДГУ и ДИФ. В отличие от классических методов, где собственные числа находятся из решения задачи Штурма – Лиувилля, интегрируемой в области определения пространственной переменной, в данном случае они находятся из решения временно́го обыкновенного дифференциального уравнения, описывающего изменение во времени дополнительной искомой функции.
Постоянная интегрирования C1 находится из выполнения решением (27) начального условия (7) путем использования метода наименьших квадратов. Целью применения этого метода является исключение нахождения сложных интегралов, получающихся в результате выполнения требования ортогональности невязки начального условия ко всем координатным функциям. Использование метода наименьших квадратов позволяет максимально арифметизировать процесс выполнения начального условия, определяя постоянные интегрирования с заданной точностью при сохранении решения в аналитическом виде. При этом процесс выполнения начального условия максимально упрощается, т. к. наиболее трудоемкая часть расчетов, связанных с определением постоянных интегрирования, перекладывается на ЭВМ.
Потребуем выполнения начального условия (7) в следующих 10 точках области: ξi = 0,1; 0,2; 0,3; ...; 0,9; 1,0, (
). Точка ξ = 0 не входит в состав принятых точек ввиду невозможности одновременного выполнения в ней отличающихся друг от друга условий (7) и (8) в начальный момент времени Fo = 0. Подставив (27) в (7) и записав полученное соотношение для 10 выделенных точек, относительно постоянной интегрирования C1 получим систему 10 алгебраических линейных уравнений. Данная система является переопределенной (число уравнений превышает число неизвестных). Найдем ее решение в смысле наименьшего квадратичного уклонения, т. е. используем метод наименьших квадратов. Определяя сумму квадратов начального условия в выделенных точках, получаем следующий функционал:
. (28)
Подставим (27) в (28):
. (29)
Для нахождения минимума функционала (29) необходимо найти производную по неизвестному C1 и приравнять полученное соотношение к нулю:
(30)
(30)Соотношение (30) представляет алгебраическое уравнение относительно C1. Из его решения находим
. 
(31)
Из (31) находим C1 = 1,98352. Точная его величина составляет C1 = 2,0 [11]. Увеличивая число точек аппроксимации начального условия до 100, для C1 получаем значение, отличающееся от его точной величины лишь в четвертом знаке после запятой (C1 = 1,9998). Однако уже полученная точность при десяти точках аппроксимации вполне достаточна для решения задачи в первом приближении.
Подставляя C1 в (27), получаем первое приближение задачи (6)–(9):
. (32)
Решением (32) точно удовлетворяются уравнение (6) и условия (8), (9) и приближенно – начальное условие (7). Расчеты по формуле (32) в диапазоне
0,1 ≤ Fo ≤ ∞ отличаются от точного решения [12] на 4 %.
Во втором приближении, подставляя (11) в (10), (23) (i = 1) для b1(q) и b2(q), получаем следующие формулы:
где 
.
Подставим (33) и (34) в (11):
(35)
Подставив (35) в (6) и применив полученное соотношение к любой точке переменной ξ (исключая точки ξ = 0 и ξ = 1), находим
,
Интеграл уравнения (36) имеет вид
, (37)
где C1, C2 – постоянные интегрирования; v1 = π2, v2 = 4π2 – собственные числа, равные точным их значениям [12].
Подставим (37) в (35):
(38)
Соотношение (38), независимо от величин C1 и C2, точно удовлетворяет граничным условиям (8), (9) и уравнению (6) во всем диапазоне изменения временно́й и пространственной переменных, включая точки ξ = 0 и ξ = 1. Выполнение уравнения (6) во всем диапазоне пространственной переменной, минуя процесс непосредственного интегрирования по ней, является следствием того, что соотношение (38) точно удовлетворяет ДГУ (15), (22), (23). Отметим, что эти условия выполняются n при любом другом методе получения решения задачи (6)–(9), – отличие лишь в том, что они не выделяются в них в качестве условий для отдельного рассмотрения.
Для определения постоянных интегрирования C1 и C2, аппроксимируя начальное условие в 100 точках области (0 < ξ < 1), получаем переопределенную систему 100 алгебраических уравнений с неизвестными C1 и C2. Используя метод наименьших квадратов, получаем: C1 = 1,9998; C2 = –1,9998. Полученные величины C1 и C2 в пятом знаке после запятой отличаются от точных их значений. Точность определения C1 и C2 можно повысить, увеличивая число точек аппроксимации начального условия. В дальнейшем будем считать, что постоянные интегрирования найдены с точностью до пятого знака после запятой (по сравнению с точными их значениями), и записывать в решения точные значения C1 и C2. Соотношение (38) с учетом C1 и C2 будет
(39)
Соотношение (39) представляет второе приближение задачи (6)–(9). Оно точно удовлетворяет уравнению (6), граничным условиям (8), (9) и приближенно – начальному условию (7). Расчеты по формуле (39) в диапазоне 0,1 ≤ Fo ≤ ∞ отличаются от точного решения [12] на 1 %.
В третьем приближении неизвестные b1(q), b2(q), b3(q)) находятся из (10), (23) (при i = 1, 2). Из решения системы 3-х алгебраических уравнений для них получены следующие формулы:
Подставляя (11) (с учетом b1(q), b2(q), b3(q)) в уравнение (6) и применяя полученное соотношение для любой точки переменной ξ (исключая точки ξ = 0 и ξ = 1), получаем
Интеграл уравнения (40) будет
где v3 = 9π2 – третье собственное число, равное его точному значению [12]; C1, C2, C3 – постоянные, определяемые из (7). Аппроксимируя (7) в 100 точках пространственной переменной ξ (исключая точку ξ = 0) и используя метод наименьших квадратов, получаем их значения, совпадающие с точными до четвертого знака после запятой. Учитывая, что точность определения C1, C2, C3 можно повысить (путем увеличения числа точек аппроксимации начального условия), далее будем использовать их точные значения, равные C1 = 2, C2 = –2, C3 = 2.
С учетом полученных bk(q), функции q(Fo) и C1, C2, C3 решение в третьем приближении принимает вид
,(41)
где vk = kπ.
Формула (41) справедлива и в последующих приближениях k = 1, 2, 3, ..., n. Общая формула для них будет
(42)
Соотношение (42) при n → ∞ совпадает с классическим точным аналитическим решением задачи (6)–(9) [12]. Расчеты по формуле (42) в различных приближениях приведены на рис. 2.
Рис. 2. Изменение температуры: – 1-е приближение;
– 2-е приближение;
– точное решение [12]
Fig. 2. Temperature changes: – 1st approximation;
– 2nd approximation;
– exact solution [12]
Приближение к формуле для коэффициентов Ak с заданной точностью выполняется методом наименьших квадратов.
Заключение
Путем использования ДГУ и ДИФ получено решение задачи теплопроводности для пластины с неоднородными граничными условиями первого рода, приближающееся точному решению при n → ∞. Применение ДГУ позволяет выполнить исходное уравнение во всей области пространственной переменной без проведения непосредственного интегрирования по ней, в связи с чем процесс получения аналитического решения сводится лишь к интегрированию обыкновенного уравнения относительно ДИФ (изменяющейся лишь во времени), из которого находятся собственные числа, определяемые в классических методах из задачи Штурма – Лиувилля, заданной в пространственных переменных. Следовательно, в настоящей работе приводится другой способ определения собственных чисел, основанный на решении временного уравнения для ДИФ.
Постоянные интегрирования находятся путем выполнения начального условия методом наименьших квадратов, позволяющим находить их значения практически с заданной точностью. Преимущество такого способа их определения заключается в отсутствии необходимости нахождения интегралов в пределах области изменения пространственной переменной при сохранении аналитического вида решения.
1. Beliaev N. M., Riadno A. A. Metody nestatsionarnoi teploprovodnosti [Methods of non-stationary thermal conductivity]. Moscow, Vysshaia shkola Publ., 1978. 328 p.
2. Bio M. Variatsionnye printsipy v teorii teploobmena [Variational principles in the theory of heat transfer]. Mos-cow, Energiia Publ., 1975. 208 p.
3. Glazunov Iu. T. Variatsionnye metody [Variational methods]. Moscow, Izhevsk, NITs «Reguliarnaia i khaoti-cheskaia dinamika»; In-t komp'iuter. issled. Publ., 2006. 470 p.
4. Gudmen T. Primenenie integral'nykh metodov v nelineinykh zadachakh nestatsionarnogo teploobmena [Integral methods application in nonlinear problems of non-stationary heat transfer]. Problemy teploobmena: sbornik nauchnykh trudov. Moscow, Atomizdat, 1967. Pp. 41-96.
5. Kudinov V. A., Eremin A. V., Trubitsyn K. V., Stefaniuk E. V. Modeli termomekhaniki s konechnoi i besko-nechnoi skorost'iu rasprostraneniia teploty [Thermomechanics models with finite and infinite heat propagation velocity]. Moscow, Prospekt Publ., 2020. 224 p.
6. Kudriashov L. I., Men'shíkh N. L. Priblizhennye resheniia nelineinykh zadach teploprovodnosti [Approximate solutions of nonlinear heat conduction problems]. Moscow, Mashinostroenie Publ., 1979. 232 p.
7. Kudinov I. V., Kotova E. V., Kudinov V. A. Metod polucheniia analiticheskikh reshenii kraevykh zadach na osnove opredeleniia dopolnitel'nykh granichnykh uslovii i dopolnitel'nykh iskomykh funktsii [A method for obtaining analytical solutions to boundary value problems based on the determination of additional boundary conditions and additional sought functions]. Sibirskii zhurnal vychislitel'noi matematiki, 2019, vol. 22, no. 2, pp. 153-165.
8. Fedorov F. M. Granichnyi metod resheniia priklad-nykh zadach matematicheskoi fiziki [Boundary method for solving applied problems of mathematical physics]. Novosibirsk, Nauka Publ., 2000. 220 p.
9. Kantorovich L. V. Ispol'zovanie idei metoda Galerkina v metode privedeniia k obyknovennym differentsial'nym uravneniiam [Using the idea of the Galerkin method in the method of reduction to ordinary differential equations]. Prikladnaia matematika i mekhanika, 1942, vol. 6, no. 1, pp. 31-40.
10. Kantorovich L. V., Krylov V. I. Priblizhennye metody vysshego analiza [Approximate methods of higher analysis]. Leningrad, Fizmatgiz Publ., 1962. 708 p.
11. Kartashov E. M. Analiticheskie metody v teorii tep-loprovodnosti tverdykh tel [Analytical methods in the theory of thermal conductivity of solids]. Moscow, Vysshaia shkola Publ., 2001. 550 p.
12. Lykov A. V. Teoriia teploprovodnosti [Theory of thermal conductivity]. Moscow, Vysshaia shkola Publ., 1967. 600 p.
