Obtaining exact analytical solutions for nonstationary heat conduction problems using orthogonal methods

Full Text

При получении приближенных аналитических решений большое распространение получили ортогональные методы Л. В. Канторовича, Бубнова- Галеркина и интегральный метод теплового баланса, важное преимущество которых перед классическими аналитическими методами (Фурье, интегральных преобразований, функций Грина и др.) состоит в том, что при их использовании на вид дифференциального оператора краевой задачи не накладывается каких-либо ограничений [1, 2]. И, в частности, они могут быть применены к интегрированию дифференциальных уравнений, не допускающих разделения переменных [3, 4], например, нелинейных, с переменными физическими свойствами среды и др. Решение находится в виде степенных (или тригонометрических) полиномов по координатам текущей точки с экспоненциально стабилизирующимися во времени коэффициентами. В результате применения этих методов для линейных и нелинейных краевых задач можно получить простые по форме аналитические решения, в том числе и точные (для линейных задач), содержащие все основные физические свойства среды в явном виде. В настоящей работе показывается, как на основе ортогональных методов Л. В. Канторовича, Бубнова-Галеркина и интегрального метода теплового баланса можно получать не только приближенные, но и точные аналитические решения задач нестационарной теплопроводности. Рассмотрим идею метода на примере решения краевой задачи для бесконечной пластины толщиной δ в следующей постановке: ∂T (x, t) ∂ 2 T (x, t) =a , t > 0, 0 < x < δ; ∂t ∂x2 T (x, 0) = T0 , ∂T (0, t)/∂x = 0, T (δ, t) = Tст . (1) (2) Здесь T - температура, x - координата, t - время, Tст - температура на стенках пластины, T0 - начальная температура, a - коэффициент температуропроводности. Для приведения задачи (1), (2) к безразмерному виду введем следующие обозначения: Θ = (T - Tст )/(T0 - Tст ), Fo = at/δ 2 , ξ = x/δ, (3) где Θ, Fo, ξ - соответственно, безразмерные температура, время и координата. Задача (1), (2) в обозначениях (3) принимает вид ∂ 2 Θ(ξ, Fo) ∂Θ(ξ, Fo) = , Fo > 0, 0 < ξ < 1; ∂Fo ∂ξ 2 Θ(ξ, 0) = 1; ∂Θ(0, Fo) = 0, Θ(1, Fo) = 0. ∂ξ 198 (4) (5) (6) Получение точных аналитических решений нестационарных задач теплопроводности. . . Следуя методу Л. В. Канторовича, решение задачи (4)-(6) ищется в виде n Θ(ξ, Fo) = fk (Fo)ϕk (ξ), (7) k=1 где здесь и далее fk (Fo) - неизвестные функции времени; ϕk (ξ) - координатные функции. В работах [3, 4] в качестве координатных принимаются алгебраические полиномы. Так как алгебраические координатные функции не обладают свойством ортогональности, то при их использовании могут быть получены лишь приближенные аналитические решения. Найдем решение задачи (4)-(6) при помощи тригонометрических координатных функций вида (2k - 1)π ξ . 2 ϕk (ξ) = cos (8) Очевидно, что (7) с учетом (8) удовлетворяет граничным условиям (6). Потребуем, чтобы уравнение (4) на отрезке пространственной переменной (0 ξ 1) выполнялось в среднем, то есть потребуем выполнения следующего интеграла теплового баланса: 1 0 ∂Θ(ξ, Fo) dξ = ∂Fo 1 0 ∂ 2 Θ(ξ, Fo) dξ. ∂ξ 2 (9) Для нахождения решения задачи (4)-(6) в первом приближении подставим (7), ограничиваясь одним членом ряда, в (9): 1 0 ∂ π f1 (Fo) cos ξ ∂Fo 2 1 dξ - 0 ∂2 π f1 (Fo) cos ξ 2 ∂ξ 2 dξ = 0. (10) Соотношение (10) относительно неизвестной функции f1 (Fo) приводится к обыкновенному дифференциальному уравнению df1 (Fo) π 2 + f1 (Fo) = 0, dFo 4 после интегрирования которого находим f1 (Fo) = C1 exp - π2 Fo , 4 (11) где C1 - постоянная интегрирования. Подставляя (8), (11) в (7), находим Θ(ξ, Fo) = C1 exp - π2 π Fo cos ξ . 4 2 (12) Для определения постоянной интегрирования C1 составим невязку начального условия (5) и потребуем ортогональности невязки к координатной функции ϕ1 (ξ): 1 π π C1 cos ξ - 1 cos ξ dξ = 0. (13) 2 2 0 199 К у д и н о в В. А., К л е б л е е в Р. М., К у к л о в а Е. А. Определяя интегралы в (13), относительно постоянной C1 получаем алгебраическое линейное уравнение, из решения которого находится C1 = 4/π. Соотношение (12) с учетом найденного значения постоянной интегрирования принимает вид Θ(ξ, Fo) = 4 π exp(-µ1 Fo) cos ξ , π 2 (14) где µ1 = π 2 /4. Отметим, что µ1 совпадает с первым собственным числом краевой задачи Штурма-Лиувилля при решении задачи (4)-(6) классическими аналитическими методами [5]. Соотношение (14) представляет решение задачи (4)-(6) в первом приближении. Оно точно удовлетворяет уравнению (4) и граничным условиям (6) и в первом приближении - начальному условию (5). Результаты расчетов по формуле (14) в сравнении с точным решением [5] приведены на рис. 1. Анализ численных данных позволяет заключить, что для Fo 0.1 расхождение результатов не превышает 5 %. Для решения задачи (4)-(6) во втором приближении составим невязку уравнения (4) и потребуем ортогональности невязки к координатным функциям ϕ1 (ξ) и ϕ2 (ξ) (реализация ортогонального метода Бубнова-Галеркина): 1 0 0 1 ∂Θ(ξ, Fo) ∂ 2 Θ(ξ, Fo) - ϕ1 (ξ)dξ = 0, ∂Fo ∂ξ 2 ∂Θ(ξ, Fo) ∂ 2 Θ(ξ, Fo) - ϕ2 (ξ)dξ = 0. ∂Fo ∂ξ 2 (15) Рис. 1. Решение задачи (4)-(6) при различных Fo: сплошные линии - точное решение, - первое приближение, ◦ - второе приближение, - третье приближение; 1 - Fo = 0.02, 2 - Fo = 0.05, 3 - Fo = 0.1, 4 - Fo = 0.2, 5 - Fo = 0.4, 6 - Fo = 0.6 [Figure 1. The solution of the problem (4)-(6) for different Fo: solid lines - the exact solution, - the first approximation, ◦ - the second approximation, - the third approximation; 1 - Fo = 0.02, 2 - Fo = 0.05, 3 - Fo = 0.1, 4 - Fo = 0.2, 5 - Fo = 0.4, 6 - Fo = 0.6 ] 200 Получение точных аналитических решений нестационарных задач теплопроводности. . . Подставляя (7) (ограничиваясь двумя членами ряда) в (15), находим 1 f1 (Fo)ϕ1 (ξ) + f2 (Fo)ϕ2 (ξ) + f1 (Fo)ϕ1 (ξ) + f2 (Fo)ϕ2 (ξ) ϕ1 (ξ)dξ = 0; 0 (16) 1 f1 (Fo)ϕ1 (ξ) + f2 (Fo)ϕ2 (ξ) + f1 (Fo)ϕ1 (ξ) + f2 (Fo)ϕ2 (ξ) ϕ2 (ξ)dξ = 0, 0 где f1 (Fo) = df1 (Fo)/dFo, f2 (Fo) = df2 (Fo)/dFo; ϕ1 (ξ) = d2 ϕ1 (ξ)/dξ 2 , ϕ2 (ξ) = d2 ϕ2 (ξ)/dξ 2 . После определения интегралов соотношения (16) относительно неизвестных функций времени f1 (Fo) и f2 (Fo) в силу ортогональности косинусов получаем следующие раздельные обыкновенные дифференциальные уравнения относительно f1 (Fo) и f2 (Fo): π2 f1 (Fo) = 0, 4 которые легко интегрируются: f1 (Fo) + f1 (Fo) = C1 exp - π2 Fo , 4 f2 (Fo) + 9π 2 f2 (Fo) = 0, 4 f2 (Fo) = C2 exp - 9π 2 Fo , 4 (17) 3π ξ . 2 (18) где C1 , C2 - постоянные интегрирования. Представляя (17) в (7), получаем Θ(ξ, Fo) = C1 e- π2 Fo 4 cos 9π 2 π ξ + C2 e- 4 Fo cos 2 Для определения постоянных интегрирования в (18) составим невязку начального условия (5) и потребуем ортогональности невязки к координатным функциям ϕ1 (ξ) и ϕ2 (ξ). Ввиду ортогональности косинусов неизвестные постоянные C1 и C2 в соотношениях для невязки разделятся: 1 C1 cos2 0 1 C2 cos2 0 π π ξ - cos ξ dξ = 0; 2 2 3π 3π ξ - cos ξ dξ = 0. 2 2 (19) (20) Определяя интегралы в (19), (20), находим C1 = 4/π; C2 = -4/(3π), c учетом которых соотношение (18) принимает вид Θ(ξ, Fo) = 4 -µ1 Fo π 4 -µ2 Fo 3π e cos ξ - e cos ξ , π 2 3π 2 (21) где µ1 = π 2 /4; µ2 = 9π 2 /4. Соотношение (21) представляет решение задачи (4)-(6) во втором приближении. Оно точно удовлетворяет уравнению (4) и граничным условиям (6) и приближенно (во втором приближении) - начальному условию (5). Отметим, что µ1 и µ2 совпадают с двумя первыми собственными числами, получаемыми при решении задачи (4)-(6) классическими аналитическими методами, а 201 К у д и н о в В. А., К л е б л е е в Р. М., К у к л о в а Е. А. константы интегрирования C1 и C2 совпадают с двумя первыми коэффициентами, получаемыми из выполнения начального условия в точном аналитическом решении. Результаты расчётов по формуле (21) в сравнении с точным аналитическим решением [5] приведены на рис. 1. Анализ численных данных позволяет заключить, что во втором приближении, по сравнению с первым, происходит существенное уточнение получаемого решения (при Fo 0.05 оно практически совпадает с точным). Найдем решение задачи (4)-(6) для любого приближения, то есть при n → ∞. Составим невязку уравнения (4) и потребуем ортогональности невязки ко всем координатным функциям ϕk (ξ): 1 0 ∞ k=1 dfk (Fo) ϕk (ξ) + µk dFo 1 ∞ fk (Fo)ϕk (ξ) ϕj (ξ)dξ = 0, j = 1, 2, . . . , (22) 0 k=1 где µk = (2k-1)2 π 2 /4. Соотношения (22) с учетом ортогональности косинусов приводятся к виду dfk (Fo) + µk fk (Fo) = 0, dFo k = 1, 2, . . . . (23) k = 1, 2, . . . , (24) После интегрирования (23) получаем fk (Fo) = Ck exp(-µk Fo), где Ck - постоянные интегрирования. Подставим (24) в (7) (при n → ∞): ∞ Θ(ξ, Fo) = Ck exp(-µk Fo)ϕk (ξ). (25) k=1 Для определения постоянных интегрирования Ck составим невязку начального условия (5) и потребуем ортогональности невязки ко всем координатным функциям ϕk (ξ): ∞ 1 Ck ϕk (ξ) - 1 ϕj (ξ)dξ = 0, 0 j = 1, 2, . . . . k=1 Последние соотношения ввиду ортогональности координатных функций (8) принимают вид 1 Ck ϕ2k (ξ) - ϕk (ξ) dξ = 0, k = 1, 2, . . . . 0 Определяя интегралы в (26), находим Ck : 2 Ck = (-1)k+1 √ , µk 202 k = 1, 2, . . . . (26) Получение точных аналитических решений нестационарных задач теплопроводности. . . Подставляя найденные Ck в (25), получаем формулу ∞ Θ(ξ, Fo) = k=1 2 (-1)k+1 √ exp(-µk Fo)ϕk (ξ), µk которая представляет собой точное аналитическое решение задачи (4)-(6). На рис. 2 в качестве иллюстрации предложенного метода приведено точное решение (сплошные линии) задачи (4)-(6) и решение, полученное по формуле (25), когда используется 100 членов ряда, при различных значениях Fo. Рис. 2. Точное решение (сплошные линии) задачи (4)-(6) и решение, полученное по формуле (25), когда используется 100 членов ряда, при различных значения Fo: 1 - Fo = 0.0001, 2 - Fo = 0.001, 3 - Fo = 0.0005, 4 - Fo = 0.02, 5 - Fo = 0.6 [Figure 2. The exact solution (solid lines) of the problem (4)-(6) and the solution obtained by the formula (25), when 100 members of the series are used for different Fo: 1 - Fo = 0.0001, 2 - Fo = 0.001, 3 - Fo = 0.0005, 4 - Fo = 0.02, 5 - Fo = 0.6] Выводы. На основе ортогональных методов Л. В. Канторовича, Бубнова-Галеркина и интегрального метода теплового баланса при использовании тригонометрических координатных функций получено точное аналитическое решение задачи теплопроводности для бесконечной пластины при симметричных граничных условиях первого рода. Благодаря ортогональности тригонометрических функций неизвестные функции времени в системе обыкновенных дифференциальных уравнений разделяются - в каждое уравнение входит лишь одна неизвестная величина. В итоге удается получить решение общего вида для неизвестных функций времени при любом числе приближений. Благодаря этому же свойству координатных функций при определении из начального условия констант интегрирования удается избежать решения больших систем алгебраических линейных уравнений с плохо обусловленными матрицами коэффициентов и получить общую формулу для нахождения этих констант. 203 К у д и н о в В. А., К л е б л е е в Р. М., К у к л о в а Е. А. Преимущества такого метода получения аналитического решения в том, что при совместном использовании указанных методов на вид дифференциального уравнения краевой задачи практически не накладывается каких-либо ограничений. Причем собственные числа определяются из решения дифференциального уравнения относительно неизвестных функций времени (уравнение первого порядка), а не из краевой задачи Штурма-Лиувилля (уравнение второго порядка), решение которой требуется находить при использовании классических точных аналитических методов. Ортогональные методы с тригонометрическими координатными функциями могут быть использованы и при решении краевых задач, дифференциальные уравнения которых не допускают разделения переменных (нелинейных, с переменными физическими свойствами среды и др.). Однако применительно к ним могут быть получены лишь приближенные аналитические решения. Трудности получения точных решений в форме бесконечных рядов связаны с тем, что для нелинейных краевых задач принцип линейной суперпозиции частных решений неприменим. Конкурирующие интересы. Мы не имеем конкурирующих интересов. Авторский вклад и ответственность. Все авторы принимали участие в разработке концепции статьи и в написании рукописи. Авторы несут полную ответственность за предоставление окончательной рукописи в печать. Окончательная версия рукописи была одобрена всеми авторами.

About the authors

Vasiliy A Kudinov

Samara State Technical University

Email: totig@yandex.ru
Dr. Phys. & Math. Sci., Professor; Head of Dept.; Dept. of Theoretical Fundamentals of Heat-Engineering and Hydromechanics 244, Molodogvardeyskaya st., Samara, 443100, Russian Federation

Ruslan M Klebleev

Samara State Technical University

Email: uio1123@list.ru
Graduate Student; Dept. of Theoretical Fundamentals of Heat-Engineering and Hydromechanics 244, Molodogvardeyskaya st., Samara, 443100, Russian Federation

Ekaterina A Kuklova

Samara State Technical University

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

References

Supplementary files