Аналитический и численный методы решения уравнения теплопроводности Текст научной статьи по специальности «Математика»

ИНФОРМАТИКА И ТЕХНИЧЕСКИЕ НАУКИ

МОДЕЛИРОВАНИЕ ПРОЦЕССОВ И УПРАВЛЕНИЕ В ТЕХНИЧЕСКИХ СИСТЕМАХ

УДК 681.53

Д. С. Карпович, О. Н. Суша, Н. П. Коровкина, В. П. Кобринец

Белорусский государственный технологический университет

АНАЛИТИЧЕСКИЙ И ЧИСЛЕННЫЙ МЕТОДЫ РЕШЕНИЯ УРАВНЕНИЯ ТЕПЛОПРОВОДНОСТИ

В данной статье приводится сравнение методов решения уравнения теплопроводности. Подробно рассмотрено решение аналитическим методом. Заданы условия однозначности, а также начальные и граничные условия и приведены способы их задания с учетом физических особенностей моделирования теплопроводности режущего инструмента. Составлено и решено уравнение методом разделения переменных, в виде произведения двух функций, и разложено в ряд Фурье с заданными параметрами, определяемыми характеристическим уравнением. Получено окончательное выражение распределения температуры в инструменте. Также приведен пример решения численным методом тепловых балансов, выведено уравнение в конечно-разностной форме для расчета распределения температурного поля и получено приближенное решение для температур в узловых точках. Проанализированы характерные особенности каждого метода решения одномерной нестационарной задачи теплопроводности. Представлены графики распределения температуры в инструменте для интервала времени с разным количеством элементов ряда. Сделан вывод о точности и вычислительной сложности при решении каждой рассмотренной задачи.

В заключение раскрываются достоинства и недостатки аналитического и численного методов и приводится обоснование использования модифицированного численного метода в одномерной нестационарной задаче теплопроводности.

Ключевые слова: метод разделения переменных, дифференциальное уравнение, метод конечных разностей, закон Фурье, температурное поле, граничные условия.

D. S. Karpovich, O. N. Susha, N. P. Korovkina, V. P. Kobrinets

Belarusian State Technological University

THE ANALYTICAL AND NUMERICAL METHODS

OF SOLVING THE THERMAL CONDUCTIVITY EQUATION

In this article we present a comparison of the methods of solving the thermal conductivity equation. An analytical method of solving is considered in details. The conditions of unambiguity are specified, as well as the initial and boundary conditions and ways of their definition in consideration with physical features of modelling the thermal conductivity of the cutting tool. The equation is constituted and solved by the method of division of variables, as the product of two functions, and it is separated in Fourier series with the specified parameters defined by the characteristic equation. The final expression of temperature division in the tool is obtained as well. The example of solving by a numerical method of thermal balances is given too, the equation is deduced into finite-difference form to calculate

the distribution of the temperature field and the approximate solution for temperatures in central points is obtained as well. There are analyzed characteristic features of each method of solving of one-dimensional non-stationary problem of thermal conductivity. Diagrams of the temperature distibution in the tool are presented for a time interval with a different quantity of elements of series row. We made a conclusion on accuracy and computing complexity while solving each considered problem.

In conclusion we revealed the advantages and disadvantages of analytical and numerical methods and we give proof of using the modified numerical method in solving one-dimensional non-stationary problem of thermal conductivity.

Key words: the method of division of variables, the differential equation, the method of ultimate differences, Fourier' law, a temperature field, boundary conditions.

Введение. Существует несколько методов решения задач теплопроводности: аналитический, аналоговый, численный, графический и экспериментальный. Четыре из них исходят непосредственно из различных форм уравнений. Экспериментальным методом пользуются, когда остальные методы не дают результатов. Его применяют для определения теплофизических свойств, таких как теплопроводность и удельная теплоемкость.

Для решения задач теплопроводности в твердых телах сложной формы используются аналитические и численные методы. Решения возможны при известных краевых условиях, включающих начальное распределение температур в теле и граничные условия на поверхности тела, которые могут быть заданы одним из трех способов: температурой поверхности, тепловым потоком и коэффициентом теплоотдачи [1].

Аналитическое решение задач нестационарной теплопроводности методом разделения переменных Фурье предпочтительнее других методов при неравномерном начальном распределении температуры в теле и в тех случаях, когда нет необходимости в расчетах для очень малых времен от начала процесса.

Из численных методов решения задач теплопроводности в настоящее время наиболее ценным и широко используемым является метод конечных разностей.

Основная часть. Дифференциальное уравнение теплопроводности при отсутствии внутренних источников теплоты имеет вид

(

d т

• = a

t

dx2 dy

i2t ^

(1)

Данное уравнение известно как закон (или постулат) Фурье.

Условия однозначности задаются в виде: физических параметров X, с, р; формы и геометрических размеров объекта 10, А, /2, ..., 4; температуры тела в начальный момент времени т = 0:

* = *о = / (х, У,2).

Граничные условия могут быть заданы в следующем виде:

t (т = 0) = tHa4,

t (Х = Хнач ) = ^ t (Х = Хкон ) = t2 *

Решить задачу теплопроводности - значит установить зависимость между температурой t, временем т и координатами тела x, y, z.

Рассмотрим режущий инструмент толщиной 25. Если толщина инструмента мала по сравнению с длиной и шириной, то такую полоску обычно считают неограниченной.

При заданных граничных условиях коэффициент теплоотдачи а одинаков для всех точек поверхности инструмента. Изменение температуры происходит только в одном направлении x, в двух других направлениях температура не изменяется, следовательно, в пространстве задача является одномерной. Начальное распределение температуры задано некоторой функцией t(x, 0) = fx). Охлаждение происходит в среде с постоянной температурой tw = const. Отсчет температуры инструмента для любого момента времени будем вести от температуры окружающей среды. Поскольку задача в пространстве одномерная, то дифференциальное уравнение (1) принимает вид

д 2-Э — = a-дт

д2 x'

(2)

Начальные условия: при т = 0

3 = -Эо = / (х) -*ж = Р (х). (3)

Решение дифференциального уравнения (2) ищем в виде произведения двух функций, из которых одна является функцией только т, а другая - только х (метод разделения переменных):

S = &(т, x) = ф(т)у( x).

(4)

После подстановки последнего выражения в дифференциальное уравнение (2) имеем:

ф'(т) = a V''(x) ф(т) V(x)

(5)

Левая часть уравнения (5) есть функция только т, а правая - функция только х.

Решая уравнение (5), находим частное решение:

3 = Ле-

С08(&х).

(6)

Подчинив частное решение граничному условию, после простейших преобразований получим значения корней характеристического уравнения:

31 = Л1 С08| ц — I е

2 ат 82 .

(7)

ч _ц2 ат

32 = Л.2С08|Ц28|е 2 82 ;

2 ат

3п = Л Н Цп 8| е ^82

Полученные частные решения будут удовлетворять дифференциальному уравнению при любых значениях постоянных, но ни одно из этих решений не будет соответствовать действительному распределению температуры в начальный момент времени [3]. Однако путем наложения бесконечного числа таких распределений при соответствующем выборе величин Ап можно воспроизвести любую действительную температурную зависимость в начальный момент времени.

На основании сказанного общее решение можно представить суммой бесконечного ряда:

Л

С08

ч 2 ат X Ч ЦП "Г

Цп8)е 8 .

(8)

Известно, что если отдельные распределения (7) удовлетворяют дифференциальному уравнению (2) и граничным условиям, то и сумма их также удовлетворяет тем же условиям. Постоянная Ап в уравнении найдется из начальных условий.

Окончательное выражение для случая равномерного распределения температуры в инструменте в начальный момент времени:

3=1

3о 2зт Цп

Ц п

| Цп7 |е 81И Цп С08 Цп I 8

2 ат 82

(9)

Следует отметить, что если получение точного аналитического решения связано с трудностью удовлетворения граничных условий, которые не всегда осуществимы, то при помощи численного метода всегда возможно, по крайней мере приближенно, удовлетворить граничным условиям конкретной задачи.

Получим расчетную формулу для численного интегрирования одномерной нестацио-

нарной задачи методом тепловых балансов. Пусть в этом случае процесс теплопроводности описывается уравнением

=

дт " дх2

(10)

Рассмотрим применение метода к расчету температурного поля в режущем инструменте (уравнение (10)).

Разбиваем на элементарные объемы. Полагаем, что удельная теплоемкость с и коэффициент теплопроводности X в пределах элементарного участка постоянны. Очевидно, количество теплоты, подводимое стержнем к узловой точке, определится по закону Фурье. Если расстояние 5 достаточно мало, то можно выразить q через конечные:

X

( = qAтF = Ы Ат^.

(11)

Изменение внутренней энергии в рассматриваемой узловой точке за время Дт:

и = ерУА!' = ерУ(!' - t),

(12)

где с - удельная теплоемкость; р - плотность вещества; У - элементарный объем; ! - температура в данной узловой точке в момент времени т; !' - температура в момент времени т + Дт.

С учетом сказанного уравнение теплового баланса для узловой точки будет иметь вид

2

21

(31 = ^Дг(4 - о.

А!

Решая последнее уравнение относительно неизвестной температуры получаем:

( Ч

4 =

АтХ

ерУ

!2 + !3 + ДГ1

- 2и

ерУ

(13)

Если учесть, что X / ср = а - коэффициент температуропроводности вещества, У = 52 и Дта / 82 = F0 - число Фурье, то уравнение (13) принимает вид

= F '1 го

(!2 + !3 + !1)

(1 ч

V Fо)

- 2

(14)

Уравнение (14) является основой численного метода расчета нестационарной теплопроводности.

Используем полученные формулы для преобразования дифференциального уравнения к конечно-разностной форме. Преобразование проведем на примере одномерной нестацио-

п=1

п=1

нарной задачи теплопроводности безграничной стенки (уравнение (10)).

Поскольку температура ¿(х, т) является функцией двух переменных, удобно выбрать прямоугольную сетку. Весь интервал изменения х от 0 до 1 по оси абсцисс разобьем на одинаковые интервалы 5х, а отрезок времени от т = 0 до т = & разделим на равномерные интервалы 5т (рис. 1).

к5т

т 4-► §х

2 » 1 3 к \ 5т

1 -►

х = 5х х = 25х х = т5х х = I х

Рис. 1. Сетка для составления уравнений

Восстановленные перпендикуляры к координатным осям в точках деления при пересечении образуют расчетные узловые точки. Тогда температура для узловой точки 1 с координатами х = тЪх и т = к5т запишется так:

¿1 = {1(т5х, £5Т ) =

т,к'

Для точки 2 с координатами х = тЪх и т + 5т = (к + 1)5т имеем:

¿2 = ¿2(т5 х (к + 1)5Т) = 1п

,к+1'

для точки 3 с координатами х + 5х = (т + 1)5х и т + 5т = (к + 1)5т получим:

¿3 = ¿з[( т +1) 5 х,(к +1) 5Т) = г,

т+1,к+1

и т. д.

Заменим в точке 1 (т5х, £5т) частные производные в уравнении теплопроводности разностными отношениями:

т,к

= -Т({т,к+1 - ¿т,к ) + (15)

ах2

= ТТ(^+1,к - 2(т,к + ¿т-1,к ) + 8 2, (16)

о„

В этих выражениях е1 и е2 - остаточные члены, учитывающие переход от производных функций к разностным отношениям. Можно показать, что эти члены стремятся к нулю при стремлении к нулю интегралов разбиения 5х и 5т. Дифференциальное уравнение в конечно-разностной форме запишется следующим образом:

1

(^т,к+1 ¿т,к ) + 81 =

а

Ц2

t — 2t — t

т+1,к т, к 'т-1,к

(17)

Решая уравнение (17) относительно будущей температуры ¿тк+1 в рассматриваемой точке, получим:

а5

t,

т ,к+1

2а5т

>2 I- т+1, к т-1,к

+ с

-

-1

¿т,к + (а82 -81)5х .

(18)

Очевидно, остаточный член в уравнении (18) будет стремиться к нулю при стремлении к нулю 5т. Следовательно, чем более мелкие интервалы выбраны для сетки узловых точек, тем меньше ошибка перехода от дифференциального уравнения к уравнению в конечных разностях. Ошибки в1 и е2 можно оценить, воспользовавшись разложением функции t в ряд Тейлора.

Отбрасывая остаточный член в уравнении (18) и обозначая приближенное значение величины ¿тк через Ттк, получим приближенное решение для будущей температуры в узловой точке (т5х, к5т):

Т

, к+1

т+1,к

"Тт-1,к )"

(

2а5т

Л

-1

Т

(19)

Проанализируем формулы (9) и (19). И одна, и вторая формулы являются окончательными формулами для расчетов распределения теплового поля [3]. В данной программе мы производим вычисление температуры в ста точках инструмента для пятидесяти моментов времени. Для каждой из этих точек нахождение температуры требует вычисления суммы нескольких элементов. На рис. 2 показано поле распределения температуры при количестве элементов ряда п = 10.

1401 120 -10080 -6040 20 60

00

Рис. 2. Поле распределения температур при п = 10

д

т

2

Как видно из графика, при таком значении п решение по оси х имеет колебательные свойства. Это объясняется тем, что в формуле (9) присутствует тригонометрические функции синуса и косинуса.

При увеличении числа п данные колебательные свойства будут проявляться слабее. Действительно, на рис. 3 можно наблюдать то же поле распределения температур, но при величине п = 100.

103! 102 101! 100 99 98

977

9(5 60

120

100

10 00 20

Рис. 3. Поле распределения температур при п = 100

Сравнивая эти графики, можно прийти к выводу о том, что для получения распределения температуры в инструменте приемлемой точности необходимо выбирать достаточно большие значения п. Однако при этом вычислительная сложность задачи многократно возрастает. Для случая п = 100 компьютеру необходимо провести расчет 50^100^100 раз. Вычисление формулы (9) в количестве 500 000 раз является весьма длительным процессом даже на современных специализированных быстродействующих ЭВМ. Это объясняется также тем, что, кроме большого числа вычислений формулы, необходимо каждый раз рассчитывать значения синусов и косинусов. Для вычисления синуса или косинуса ЭВМ раскладывает

его в достаточно большой ряд, что еще более усиливает сложность расчетов.

В случае же реализации формулы (19) для такого же количества точек по осям времени t и длины х вычислительная сложность задачи значительно уменьшается (рис. 4).

100

00

Рис. 4. График решения уравнения теплопроводности численным методом

Действительно, для нахождения значения температур в ста точках инструмента при пятидесяти различных моментах времени необходимо будет вычислить 5000 раз формулу (19), что в сто раз меньше по сравнению с расчетом формулы (9). Кроме того, при вычислении формулы (19) отсутствует необходимость рассчитывать значения синусов и косинусов или любых других сложных математических функций.

Заключение. Проанализировав аналитическое решение уравнения теплопроводности, можно сделать вывод, что решение удается получить только для простейших условий, для тел простой формы, при этом отмечается высокая вычислительная сложность. При помощи численного метода всегда возможно удовлетворить граничным условиям задачи и решать сложные задачи, недоступные для аналитических методов. Вычислительная сложность задачи значительно меньше.

Литература

1. Суша О. Н., Карпович Д. С. Моделирование температурного поля в дереворежущем инструменте с использованием программы ANSYS // Энерго- и ресурсосберегающие технологии и оборудование, экологически безопасные технологии. Минск, 2014. С. 44-48.

2. Кудинов В. А., Кудинов И. В. Методы решения уравнений теплопроводности. Самара, 2012. 280 с.

3. Карпович Д. С. Моделирование и численное решение уравнения теплопроводности // Энерго-и ресурсосберегающие технологии и оборудование, экологически безопасные технологии. Минск, 2014. С. 311-313.

References

1. Susha O. N., Karpovich D. S. Simulation of temperature field in the wood-cutting tools with the ANSYS. Energo- i resursosberegayushchiye tekhnologii i oborudovanie, ekologicheski bezopasnyye

tekhnologii [Energy-saving technologies and equipment, environmentally friendly technologies], Minsk, 2014, pp. 44-48 (In Russian).

2. Kudinov V. A., Kudinov I. V. Metody reshenija uravnenij teploprovodnosti [Methods for solving the heat equation], Samara, 2012. 280 p.

3. Karpovich D. S. Modeling and numerical solution of the heat equation. Energo- i resurso-sberegayushchiye tekhnologii i oborudovanie, ekologicheski bezopasnyye tekhnologii [Energy-saving technologies and equipment, environmentally friendly technologies], Minsk, 2014, pp. 311-313 (In Russian).

Информация об авторах

Карпович Дмитрий Семенович - кандидат технических наук, доцент, заведующий кафедрой автоматизации производственных процессов и электротехники. Белорусский государственный технологический университет (220006, г. Минск, ул. Свердлова, 13а, Республика Беларусь). E-mail: karpovich@tut.by

Суша Оксана Николаевна - ассистент кафедры автоматизации производственных процессов и электротехники. Белорусский государственный технологический университет (220006, г. Минск, ул. Свердлова, 13а, Республика Беларусь). E-mail: s-oxana@inbox.ru

Коровкина Наталья Павловна - кандидат педагогических наук, доцент кафедры автоматизации производственных процессов и электротехники. Белорусский государственный технологический университет (220006, г. Минск, ул. Свердлова, 13а, Республика Беларусь). E-mail: korovkina@open.by

Кобринец Виктор Павлович - кандидат технических наук, доцент кафедры автоматизации производственных процессов и электротехники. Белорусский государственный технологический университет (220006, г. Минск, ул. Свердлова, 13а, Республика Беларусь). E-mail: kobrinec@open.by

Information about the authors

Karpovich Dmitriy Semenovich - Ph. D. (Engineering), Assistant Professor, Head of the Department of Automation of Production Processes and Electrical Engineering. Belarusian State Technological University (13a, Sverdlova str., 220006, Minsk, Republic of Belarus). E-mail: karpovich@tut.by

Susha Oksana Nikolaevna - assistant, the Department of Automation of Pproduction Processes and Electrical Engineering. Belarusian State Technological University (13a, Sverdlova str., 220006, Minsk, Republic of Belarus). E-mail: s-oxana@inbox.ru

Korovkina Natal'ya Pavlovna - Ph. D. (Pedagogics), Assistant Professor, the Department of Automation of Production Processes and Electrical Engineering. Belarusian State Technological University (13a, Sverdlova str., 220006, Minsk, Republic of Belarus). E-mail: korovkina@open.by

Kobrinets Viktor Pavlovich - Ph. D. (Engineering), Assistant Professor, the Department of Automation of Production Processes and Electrical Engineering. Belarusian State Technological University (13a, Sverdlova str., 220006, Minsk, Republic of Belarus). E-mail: kobrinec@open.by

Поступила 12.03.2015