METHOD OF ADDITIONAL REQUIRED FUNCTIONS IN TASKS OF HEAT CONDUCTION WITH VARIABLE PHYSICAL PROPERTIES OF THE ENVIRONMENT

Full Text

Известные точные аналитические решения нестационарных задач теплопроводности с переменными по пространственной координате физическими свойствами среды выражаются сложными бесконечными рядами, включающими функции Бесселя первого и второго рода. Такие ряды плохо сходятся при малых значениях временной и пространственной переменных и поэтому малопригодны не только для инженерных приложений, но и для выполнения научных исследований [1, 2]. Применительно к решению краевых задач большое распространение получили ортогональные методы взвешенных невязок (Л.В. Канторовича, Бубнова - Галеркина, интегральный метод теплового баланса и др.) [3-16]. Они могут быть эффективно применены к задачам, не поддающимся решению с помощью классических аналитических методов, что связано с возможностью сведения исходного дифференциального уравнения в частных производных к интегрированию обыкновенного дифференциального уравнения, минуя процесс разделения переменных в исходном уравнении. Такая возможность основана на использовании дополнительных искомых функций, зависящих лишь от одной переменной и дополнительных граничных условий, определяемых так, чтобы искомое решение удовлетворяло дифференциальному уравнению краевой задачи в граничных точках. Основной проблемой этих методов является низкая точность, связанная с необходимостью решения алгебраических уравнений высоких степеней, а также больших систем алгебраических линейных уравнений с заполненными квадратными матрицами, которые, как правило, плохо обусловлены. Для решения указанных проблем в настоящей работе используются дополнительные искомые функции и дополнительные граничные условия, позволяющие получить сильно разреженные цепочного типа системы алгебраических линейных уравнений. Основную идею метода рассмотрим на примере решения краевой задачи для бесконечной пластины с экспоненциальным изменением коэффициента теплопроводности от координаты ; ; ; (1) ; (2) (3) (4) где - температура; - координата; - время; - коэффициент теплопроводности; - теплоемкость; - плотность; - начальная температура; - температура стенки при ; - половина толщины пластины. Найдем решение задачи (1)-(4) в случае, когда с = const и = const (в общем случае они могут быть и переменными), а коэффициент теплопроводности является экспоненциальной функцией координаты : , (5) где m > 0 - коэффициент, характеризующий интенсивность изменения коэффициента теплопроводности по координате ; = const - коэффициент теплопроводности пластины при = 0. Из формулы (5) следует, что теплопроводность пластины с увеличением координаты уменьшается, достигая минимума при . Введем следующие безразмерные переменные и параметры: (6) где - безразмерная температура; - число Фурье; - коэффициент температуропроводности пластины при = 0; - безразмерная координата. С учетом обозначений (6) задача (1)-(4) принимает вид (рис. 1) ; (7) ; (8) (9) (10) где . Рис. 1. Схема теплообмена в пластине Введем дополнительную искомую функцию (11) представляющую изменение температуры во времени в центре пластины. Температура в точке будет изменяться сразу после приложения граничного условия первого рода на поверхности пластины ( ), что связано с бесконечной скоростью распространения теплоты. Так как температура в центре пластины является искомой величиной задачи (7)-(10), то введение функции не изменяет ее, а лишь упрощает процесс получения аналитического решения. Решение задачи (7)-(10) с учетом (11) разыскивается в виде (12) где - координатные функции. Для получения решения задачи (7)-(10) в первом приближении подставим (12) (ограничиваясь тремя членами ряда) в (9)-(11). Для определения неизвестных коэффициентов , (k = 0, 1, 2) будем иметь систему трех алгебраических уравнений. После нахождения соотношение (12) будет . (13) Интеграл теплового баланса для уравнения (7) записывается в виде . (14) Подставляя (13) в (14), находим , (15) где . Интегрируя уравнение (15), получаем , (16) где - постоянная интегрирования. Подставляя (16) в (13), находим , (17) где . Для определения постоянной интегрирования подставим (17) в начальное условие (8) и потребуем ортогональности полученной невязки к функции , являющейся, по сути, собственной функцией: . (18) Определяя интеграл в (18), находим = . Соотношение (17) с учетом найденного значения постоянной интегрирования представляет решение задачи (7)-(10) в первом приближении. Оно точно удовлетворяет граничным условиям (9)-(10), интегралу теплового баланса (14) (осредненному уравнению (7)) и приближенно (в первом приближении) - уравнению (7) и начальному условию (8). Для повышения точности необходимо увеличивать число членов ряда (12), неизвестные коэффициенты которого будем находить из основных граничных условий (9), (10), соотношения (11) и некоторых дополнительных граничных условий. Для их получения применительно к точке продифференцируем соотношение (11) и граничное условие (9) по переменной : (19) (20) Соотношения (19), (20) с учетом уравнения (7) приводятся к следующим дополнительным граничным условиям: (21) (22) Для получения дополнительного граничного условия в точке продифференцируем соотношение (10) по переменной . Сравнивая полученное соотношение с уравнением (7), находим (23) Основные граничные условия (9), (10), соотношение (11) и дополнительные условия (21)-(23) используются для получения решения во втором приближении. При получении дополнительных граничных условий третьего приближения соотношения (21)-(23) дифференцируются по переменной и выполняется сравнение полученных соотношений с уравнением (7). Для получения решения задачи (7)-(10) во втором приближении, подставляя (12) (ограничиваясь шестью членами ряда) в соотношения (9)-(11), (21)-(23), относительно неизвестных коэффициентов , будем иметь цепочную систему шести алгебраических линейных уравнений. Соотношение (12) после определения принимает вид (24) где ; ; ; ; . Подставляя (24) в интеграл теплового баланса (14), после определения интегралов относительно неизвестной функции получаем обыкновенное дифференциальное уравнение вида (25) где . (26) Интегрируя уравнение (25), находим (27) где , - постоянные интегрирования. Подставляя (27) в (24), получаем (28) где . Подставляя (28) в начальное условие (8) и требуя ортогональности полученной невязки к собственным функциям и , относительно констант интегрирования и получаем систему двух алгебраических линейных уравнений: (29) Из решения системы уравнений (29) находим (30) где ; ; . Если положить , то соотношения (27), (28) приводятся к виду (31) (32) Отметим, что коэффициенты под знаком экспонент в соотношениях (31), (32) незначительно отличаются от точных значений собственных чисел краевой задачи Штурма - Лиувилля, имеющих вид = 2,4674; = 22,2066. Коэффициенты , перед экспонентами соотношения (31) незначительно отличаются от коэффициентов классического точного аналитического решения, получаемых из выполнения начального условия краевой задачи, точные значения которых = 1,2732; = 0,4241. Анализ результатов расчетов при по формуле (32) в сравнении с классическим точным аналитическим решением задачи (7)-(10) (при ) [17] позволяет заключить, что в диапазоне 0,1≤ < ∞ их расхождение не превышает 1 %. Аналогично могут быть найдены решения и для других приближений. При их получении не возникают какие-либо принципиальные трудности - возрастает лишь объем вычислительной работы. Результаты расчетов температуры по формуле (28) при в сравнении с расчетом по методу конечных разностей и по методу [4] даны на рис. 2. Их анализ позволяет заключить, что в диапазоне 0,05≤ < ∞ результаты расчетов практически совпадают. На рис. 2 представлены также результаты расчетов температуры по формуле (12) в третьем и в четвертом приближениях . Из их анализа следует, что при увеличении числа членов ряда (12) точность получаемых решений существенно возрастает. Результаты расчетов во втором приближении можно существенно уточнить, если при нахождении констант интегрирования и соотношения (28) ввести поправку на время достижения фронтом температурного возмущения координаты при его движении от поверхности пластины ( ) к ее центру ( ). Несмотря на то, что в задаче (7)-(10) заложена бесконечная скорость распространения теплоты (так как решению подлежит параболическое уравнение (7)), сколько-нибудь заметное изменение температуры в центре пластины будет происходить лишь по истечении некоторого времени , при котором определяемый ниже (см. решение задачи (37)-(40)) фронт температурного возмущения достигает центра пластины, то есть = 1. Положив в формуле (44) = 1, находим = 0,05 (при ). Вычитая это время из величины числа соотношения (27), получаем . (33) Рис. 2. Распределение температуры в пластине при : - по формуле (3.278) из [4]; 2, 3, 4 - по формуле (12) соответственно во втором, третьем и четвертом приближении; - численное решение Следовательно, отсчет времени в соотношении (33) осуществляется от , то есть является начальным временем задачи (7)-(10) при получении ее решения во втором приближении. Начальные условия для функции имеют вид ; . (34) Подставляя (33) в (34), относительно и получаем систему двух алгебраических уравнений вида (35) Из решения системы уравнений (35) при и = 0,05 находим = 1,2741; = 0,38. Как видно, произошло существенное уточнение коэффициента по сравнению с его определением из начального условия (8). Подставляя (33) в (24), с учетом найденных значений и находим (36) Анализ результатов расчетов по формуле (36) в сравнении с формулой (32) позволяет сделать вывод о возрастании точности получаемых результатов. К тому же отметим, что процесс нахождения постоянных интегрирования и из начальных условий (34) значительно упрощается по сравнению с их определением из начального условия (8). Рис. 3. Расчетная схема теплообмена Для получения решения задачи (7)-(10) при малых и сверхмалых значениях времени будем использовать метод, основанный на определении фронта температурного возмущения и дополнительных граничных условий (модель с допущением конечной скорости распространения теплоты) (рис. 3). С введением понятия движущейся границы (фронта температурного возмущения) исходная область разделяется на две подобласти и , где - функция, определяющая продвижение границы раздела во времени (см. рис. 3). Первая стадия заканчивается при достижении движущейся границей центра пластины ( ), то есть когда . Во второй стадии изменение температуры происходит по всему объему тела . Отметим, что для этой стадии процесса решение было получено выше в виде (12). Следовательно, при использовании допущения о конечной скорости распространения теплоты решение находится лишь для первой стадии процесса. Применительно к расчетной схеме, представленной на рис. 3, математическая постановка задачи будет ( ; ) (37) ; (38) ; (39) . (40) Так как на первой стадии процесса задача (37)-(40) за пределами фронта температурного возмущения не определена, то нет необходимости выполнения начального условия вида (38) по всей ширине пластины. В данном случае вполне достаточно выполнения условия . Также нет необходимости выполнения граничного условия (40), так как оно не влияет на процесс теплообмена в первой его стадии. Таким образом, при получении решения задачи для первой стадии используется граничное условие (29) и следующие условия, выполняющиеся на фронте температурного возмущения, то есть в точке : ; (41) . (42) Решение задачи (37), (39), (41), (42) принимается в виде (43) где - неизвестные коэффициенты. В первом приближении подставим (43) (ограничиваясь тремя членами ряда) в (39), (41), (42). Для определения неизвестных коэффициентов (k = 0, 1, 2) будем иметь систему трех алгебраических линейных уравнений. После их определения соотношение (43) принимает вид (44) Интеграл теплового баланса в данном случае имеет вид (45) Подставляя (44) в (45), после вычисления интегралов находим (46) Интегрируя уравнение (46), при начальном условии получаем . (47) Положив в (47) , находим время достижения фронтом температурного возмущения координаты (при ). Соотношение (44) с учетом (47) представляет решение задачи (37), (39), (41), (42) в первом приближении. Результаты расчетов по формуле (44) (при ) в сравнении с точным решением [17] приведены на рис. 4. Их анализ позволяет заключить, что в диапазоне расхождение с точным решением составляет 8 %. Для повышения точности необходимо увеличить число членов ряда (43), неизвестные коэффициенты которого находятся из основных (39), (41), (42) и дополнительных граничных условий, определяемых с использованием заданных граничных условий (39), (41), (42) и уравнения (37). Отметим, что в каждом последующем приближении вводятся три новых дополнительных граничных условия [11-15]. Формулы для их получения во втором приближении имеют вид (48) Дополнительные граничные условия находятся в таком виде, чтобы их выполнение искомым решением было эквивалентно выполнению уравнения (37) в граничной точке и на фронте температурного возмущения . Используя основные (39), (41), (42) и дополнительные (48) граничные условия, можно найти шесть неизвестных коэффициентов , ряда (43), для определения которых будем иметь цепочную систему шести алгебраических линейных уравнений. После определения из решения этой системы , соотношение (43) принимает вид (49) Подставляя (49) в интеграл теплового баланса (45), относительно неизвестной функции получаем следующее нелинейное обыкновенное дифференциальное уравнение: (50) Интегрируя уравнение (50), находим (51) где С - постоянная интегрирования, определяемая из условия . Отсюда Положив в (51) , находим время окончания первой стадии процесса во втором приближении = 0,05 (при ). Конкретные значения , найденные из соотношения (51), для приведены в таблице: 0,05 0,0247 0,0124 0,0045 0,0019 0,0007 0,00011 5∙10-6 1 0,7 0,5 0,3 0,2 0,1 0,05 0,01 Результаты расчетов перемещения фронта температурного возмущения по координате во времени позволяют заключить, что с увеличением числа приближений величина времени , при которой фронт температурного возмущения достигает координаты , уменьшается. И в пределе при n→∞ → 0. Этот результат находится в полном соответствии с гипотезой о бесконечной скорости распространения теплоты, лежащей в основе вывода параболического уравнения (37). И, в частности, для задачи с постоянными физическими свойствами среды было показано, что величина уменьшается от значения = 0,0833 в первом приближении до значения = 0,00784 в четырнадцатом приближении. Соотношения (49), (51) представляют решение задачи (37), (39), (41), (42) во втором приближении. Результаты расчетов безразмерной температуры по формуле (49) отличаются от результатов, полученных численным методом прогонки (в диапазоне чисел Фурье ), не более чем на 1,5 % (при ). Рис. 4. Распределение температуры в пластине : - первое приближение первой стадии; - 14 - е приближение первой стадии; - точное решение Если положить , то будем иметь задачу теплопроводности с постоянными физическими свойствами, решение которой можно сравнить с известным точным аналитическим решением. И, в частности, при было получено решение задачи в первом, втором, пятом, седьмом и четырнадцатом приближениях. Из их анализа следует, что в диапазоне числа Фурье значения полученных по формуле (41) температур отличаются от точных их значений в зависимости от числа приближений. И, в частности, в первом приближении - на 8 %, во втором - на 3 %, в седьмом - на 1,5 % и в четырнадцатом - на 0,1 % (см. рис. 4). Отметим, что при необходимости число приближений можно увеличить и, следовательно, имеется возможность получения решений практически с заданной степенью точности, причем без каких-либо ограничений на величину числа Фурье в области его малых и сверхмалых значений. Отметим, что при использовании точных аналитических решений в форме бесконечных рядов, полученных классическими методами, для обеспечения их сходимости при малых значениях числа Фурье необходимо использовать большое число приближений. Например, в диапазоне чисел необходимо использовать от 5000 до 500000 членов ряда точного решения. Выводы На основе использования дополнительных искомых функций и дополнительных граничных условий в интегральном методе теплового баланса получены приближенные аналитические решения задачи теплопроводности для бесконечной пластины с переменными физическими свойствами среды. Сочетание двух моделей теплопроводности - с конечной и бесконечной скоростью распространения теплоты - позволило в достаточно простом виде получить аналитические решения во всем диапазоне времени нестационарного процесса практически с заданной степенью точности, включая малые и сверхмалые его значения. 1. Дополнительной искомой функцией в модели c конечной скоростью распространения теплоты является фронт температурного возмущения, а в модели с бесконечной скоростью - функция, характеризующая изменение температуры в центре пластины. Использование дополнительных искомых функций позволяет в обеих моделях сводить решение исходного уравнения в частных производных к интегрированию обыкновенного дифференциального уравнения, что приводит к существенному упрощению процессов получения решений и окончательных формул для них. 2. Дополнительные граничные условия в обеих моделях находятся в таком виде, чтобы их выполнение искомым решением было эквивалентно выполнению исходного дифференциального уравнения в граничных точках. Показано, что выполнение уравнения в граничных точках приводит к его выполнению и внутри области с точностью, зависящей от числа приближений.

About the authors

Marina P Skvortsova

Samara State Technical University

Email: totig@yandex.ru
Assistant. 244, Molodogvardeyskaya st., Samara, 443100, Russian Federation

Igor V Kudinov

Samara State Technical University

Email: totig@yandex.ru
(Ph.D. (Techn.)), Associate Professor. 244, Molodogvardeyskaya st., Samara, 443100, Russian Federation

Ruslan M Klebleev

Samara State Technical University

Email: totig@yandex.ru
Graduate Student. 244, Molodogvardeyskaya st., Samara, 443100, Russian Federation

Tatiana E Gavrilova

Samara State Technical University

Email: totig@yandex.ru
Postgraduate Student. 244, Molodogvardeyskaya st., Samara, 443100, Russian Federation

Vladimir A Rodionov

Samara State Technical University

Email: totig@yandex.ru
Associate Professor. 244, Molodogvardeyskaya st., Samara, 443100, Russian Federation

References

Supplementary files