МИНИСТЕРСТВО ВЫСШЕГО ОБРАЗОВАНИЯ РОССИИ
АРХАНГЕЛЬСКИЙ ЛЕСОТЕХНИЧЕСКИЙ ИНСТИТУТ
К а ф е д р а т е п л о т е х н и к и
РАЗРАБОТКА ПРОГРАММЫ
ДЛЯ РЕШЕНИЯ НЕОДНОМЕРНЫХ СТАЦИОНАРНЫХ ЗАДАЧ
ТЕПЛОПРОВОДНОСТИ ЧИСЛЕННЫМ МЕТОДОМ С
ИСПОЛЬЗОВАНИЕМ КОНСЕРВАТИВНО-РАЗНОСТНОЙ СХЕМЫ
А Р Х А Н Г Е Л Ь С К
1 9 9 3
…………………………………………………………………………………………………………………………………………………………………………………
О Г Л А В Л Е Н И Е
Введение ................................………………………………….......
1.Основные положения методики построения консервативно- разностной схемы при решении неодномерных задач стационарной теплопроводности ...........…………………...........
2. Методика подготовки и решения задачи на ЭВМ ....…………...
2.1. Постановка задачи, разработка математической модели ...................................………………………………….....
2.2. Выбор метода численного решения .......…………………......
2.3. Разработка алгоритма и структуры .........…………………......
2.4. Написание программы и подготовка ее к вводу в ЭВМ .....................………………………………...............
2.5. Тестирование, отладка программы и решение на ЭВМ
Литература .......................…………………………………................
В В Е Д Е Н И Е
Базовый уровень подготовки инженера-энергетика в области информатики
и вычислительной техники определяется необходимым набором знаний, умений
и навыков в применении ЭВМ для решения различных технических задач.
Специалисты этой категории, помимо умения использовать прикладное
программное обеспечение, должны быть программирующими пользователями, т.к.
их профессиональная деятельность связана с выполнением большого количества
теплотехнических расчетов.
Для соблюдения принципа фундаментальности высшего образования работа
построена на базе рассмотрения вопросов применения ЭВМ для решения
основных задач теории теплообмена. К одной из таких задач относится задача,
связанная с определением температурного поля не одномерных тел численными методами.
Рассмотрим методику подготовки и решения указанной задачи на
персональном компьютере.
1. О С Н О В Н Ы Е П О Л О Ж Е Н И Я М Е Т О Д И К И
П О С Т Р О Е Н И Я К О Н С Е Р В А Т И В Н О-Р А З Н О С Т Н О Й С Х Е
М Ы ПРИ Р Е Ш Е Н И И Н Е О Д Н О М Е Р Н Ы Х З А Д А Ч С Т А Ц И О Н
А Р Н О Й Т Е П Л О П Р О В О Д Н О С Т И
Определение температурного поля в любой момент времени является
основной задачей теории теплопроводности. Для изотропного тела {с
постоянным по различным направлениям коэффициентом теплопроводности (} она
может быть описана дифференциальным уравнением теплопроводности
? T + Qv/( = 1/a*( dT/d(()),
(1)
где Т - температура; а - коэффициент температуропроводности, а=(/((*c);
( - плотность материала, с - удельная теплоемкость при постоянном
давлении, ? -обозначение оператора Лапласа {?= d /dx + d /dy + d /dz - в
декартовых координатах x, y, z }; ( - время, Qv - объемная плотность
теплового потока.
Уравнение теплопроводности является математическим выражением закона
сохранения энергии в твердом теле.
При решении задачи к дифференциальному уравнению теплопроводности
необходимо добавить краевые условия. В описание краевых условий входят:
поле температур для какого-нибудь предшествующего момента времени
{начальные условия}, геометрия тела {геометрические условия},
теплофизические характеристики тела {физические условия} и закон
теплообмена между поверхностью тела и окружающей средой {граничные
условия}.
Если процесс теплопроводности не только стационарный
{dT/d(tay)=0}, но и происходит без тепловыделения внутри материала (Qv =
0), то уравнение принимает вид
?(Т) = 0 .
(2)
Ввиду сложности и трудоемкости решения неодномерных задач
теплопроводности аналитическими методами в инженерной практике наиболее
часто используют приближенные. Один из них – метод конечных разностей,
непосредственно базирующийся на дифференциальном уравнении
теплопроводности и граничных условиях, представляет наибольший интерес.
В настоящее время значительное распространение получили конечно-
разностные методы, построенные с использованием известных законов
сохранения. В этом случае разностные схемы получили название
консервативные. Такой подход к построению схемы, сохраняющий физическую
сущность задачи, предпочтительнее чисто аналитического подхода,
заключающегося в непосредственной записи дифференциальных уравнений
конечно-разностными аналогами.
Следует заметить, что теория конечно-разностных численных методов
является самостоятельным разделом вычислительной математики и широко
представлена в специальной литературе[1,2,]. С основными методами
построения конечно-разностных схем, алгоритмами расчета, программным
обеспечением применительно к задачам теплообмена можно ознакомиться в
учебной литературе [3,4,5].
При изложении указанного метода особое внимание уделено физическому
смыслу построения консервативной разностной схемы и ее реализации на
ПЭВМ в задачах теплопроводности.
При использовании численного метода с консервативной разностной
схемой твердое тело разбивают на элементарные объемы. Предполагается,
что масса такого элементарного объема сосредотачивается в его центре,
называемом узлом. Для каждого узла на основе закона сохранения энергии
составляется уравнение теплового баланса, которое включает значения всех тепловых потоков на границах объемов (ячеек). Если ячейка прилегает к
поверхности тела, то выражения для определения тепловых потоков должны
описывать теплообмен между телом и окружающей средой, то есть учитывать
граничные условия. После выполнения преобразований с уравнениями
теплового баланса получают алгебраические уравнения для температуры в
каждом узле. Поскольку число узлов и число ячеек совпадают, то
образованная система алгебраических уравнений является конечно-
разностным аналогом дифференциального уравнения теплопроводности и
заменяет его с соответствующими граничными условиями. Такой подход к
составлению конечно-разностного аналога, увязанного с тепловым балансом,
позволяет получать правдоподобные решения даже при грубом выборе
расстояния между узлами (размера ячейки сетки).
Рассмотрим некоторые конкретные примеры составления конечно-
разностных схем для узлов двумерной задачи теплопроводности. В этом
случае уравнение (2) принимает вид dT/dx + dT/dy = 0 .
(3)
Внутренняя область типичного двумерного тела показана на рис.1.
Рис.1. Расположение узла внутри двумерного тела толщиной б.
Каждый элементарный прямоугольник (ячейка сетки) имеет длину -х и
высоту -у в направлениях осей х и у. Внутренний узел, обозначенный
символом 0, окружен четырьмя соседними узлами: 1,2,3,4. Кондуктивный
перенос теплоты, который в действительности происходит в твердом теле
через поверхности y*б и x*б (б -толщина тела) будем считать как перенос
теплоты от соответствующих узлов к центральному. В установившихся
условиях уравнение баланса тепловых потоков для узла 0 при отсутствии
внутреннего тепловыделения будет иметь вид
Q(1-0) + Q(2-0) + Q(3-0) + Q(4-0)
= 0 , (4)
где Q(I-0) - тепловой поток; индекс (I-0) указывает направление переноса в
узлах.
Для определения кондуктивного теплового потока может быть применен
закон Фурье
Q = - lamda * F *
dT/dn, (5)
где Т - температура, n - направление переноса теплового потока, F -
поверхность, через которую переносится тепловой поток.
Для построения расчетной схемы градиент температуры в выражении (5)
заменим разностью температур в соседних узлах. В этом случае первый член
выражения (4) примет вид
Q(1-0) = y*б*(T[1] -
T[0])/x. (6)
Здесь градиент температуры определяется на границе двух узлов 1 и 0,
имеющих температуры соответственно Т[1] и Т[0].
Аналогичные уравнения могут быть получены и для остальных трех
членов уравнения (1):
Q(2-0) = x*б*(T[2] - T[0])/y,
(7)
Q(3-0) = y*б*(T[3] - T[0])/x,
(8)
Q(4-0) = x*б*(T[4] - T[0])/y .
(9)
Точность аппроксимации градиента зависит от размера ячейки. Если
ячейка имеет квадратную форму, то уравнение теплового потока становится
независимым от формы тела.
Подставляя зависимости (6)...(9) в выражение (4), можно увидеть,
что при постоянном коэффициенте теплопроводности для квадратной сетки
(x = y) оно сводится к соотношению между температурами в рассматриваемом
узле и близлежащих:
T[1]+ T[2] + T[3] + T[4] - 4*T[0] = 0.
(10)
Выражение (10) применимо ко всем внутренним узлам.
Рассмотрим узел, расположенный на поверхности твердого тела, толщиной
б в двухмерной задаче (рис.2).
Рис.2.Расположение узлов на поверхности двумерного тела, омываемого жидкостью
Пусть узел 0, расположенный на границе твердого тела,
контактирует с окружающей средой, имеющей температуру Тc. Интенсивность
теплообмена с окружающей средой характеризуется коэффициентом теплоотдачи
alfa. Узел 0 может также обмениваться кондуктивным потоком теплоты с
тремя соседними узлами: 1,2,3. В этом случае тепловой баланс для узла 0 запишется следующим образом:
Q(1-0) + Q(2-0) + Q(3-0) + Q(c-0) = 0,
(11)
где Q(c 0)-тепловой поток, передаваемый от среды узлу 0 конвекцией.
По закону Ньютона - Рихмана
Q(c-0) = alfa*F*(T[c] - T[0]) .
(12)
В результате преобразований выражения (11), по аналогии с ранее
выполненными, для внутреннего узла, получим
y*б*(T[1] -T[0])/ x + (x/2)*б*(T[2] -T[0])/ y
+ ( x/2)*
*б*(T[3] -T[0])/ y + alfa*
y*б*(Tc -T[0]) = 0 . (13)
Соотношение (13) значительно упрощается при выборе квадратной
сетки. В этом случае при постоянном коэффициенте теплопроводности оно
приводится к виду
T[1] + 0,5*(T[2] + T[3]) + Bi*Tc - (2+Bi)*T[0] =
0, (14)
где Bi =alfa* x/lamda - число Био.
Ниже приведены уравнения теплового баланса при других граничных
условиях для двухмерных тел (x=y):
Узел Схема
Расчетное
уравнение
.....|/ Т
. 2 */ Е
. |/ П
Плоская поверх- -+--.---- + - |/ Л ность с тепло- | . |/ О изолированной x . * --?- *|/ И границей | . 1 0 |/ З
0,5(T[2] + T[3]) +
-+--.---- +- -?/ О
+ T[1] -2*T[0] = 0
. |/ Л
. ->+ x++---? x +