CC BY
22
УДК 536.2:517.958:66
А. В. Садыков, А. А. Ахметов
НЕКОТОРЫЕ РЕЗУЛЬТАТЫ ЧИСЛЕННОГО РЕШЕНИЯ УРАВНЕНИЯ ЭНЕРГИИ МЕТОДОМ СПЛАЙН-КОЛЛОКАЦИИ
Предложена сплайн-схема для расчета двумерного температурного поля в цилиндрической камере Учтена переменность теплофизических свойств в рассматриваемом объеме Проведено сравнение полученных расчетных данных с имеющимися экспериментальными ирасчетнымиданными
Тепловой расчет топочных камер трубчатых печей сводится к решению задачи радиационно-конвективного теплообмена в условиях горения топлива. Комплекс процессов радиационно-конвективного теплообмена описывается системой дифференциальных и интегродифференциальных уравнений [1]. Одним из основных среди них является уравнение энергии. Осредненное уравнение энергии, учитывающее все механизмы переноса тепла, в стационарном случае имеет вид [2]
СрIV)-^(Лз gradT ) = Г, (1)
где f - источниковый член; Л3 - коэффициент эффективной теплопроводности. Остальные обозначения имеют общепринятый смысл. С учетом уравнения неразрывности это уравнение для двумерной цилиндрической геометрии принимает вид
d(cpT) d(cpT) Э
P Vz^T~ + P v^Hr~ = 1Т Э z Э г Э z
Э T
d z
1
Г д Г
д T
д Г
+ f, (2)
где f = Яу - div Яр; уг, Уг - компоненты вектора скорости V (Яу - плотность источников тепловыделений в объеме топочной камеры; ЯР - вектор плотности радиационного потока энергии). div Яр определяется путем решения уравнения переноса
излучения. Для определения распределения Яу в объеме решаются дифференциальные
уравнения модели горения. Поле скоростей в общем случае находится решением уравнений турбулентного движения среды.
Целью данной работы является численное решение уравнения энергии. Поэтому скоростное поле считается известным, а излучением пренебрегаем.
Уравнение (2) дополняется рядом граничных условий. На входном участке в топку ставится граничное условие I рода (рис. 1а). На оси симметрии (01) задается условие симметрии Э Т / Э Г = 0, на выходном участке - нулевой градиент температуры. На жесткой границе ставятся граничные условия I или III рода.
Для решения уравнения (2) применяется метод сплайн-коллокации. Этот метод имеет преимущества по сравнению с методом конечных разностей [3]. Использование сплайна позволяет получить более подробную информацию о точном решении. Все это послужило причиной применения сплайнов при численном решении рассматриваемого уравнения.
98
Рис. 1 - Размеры камеры, границы (а); теплопередача через левый торец (граница Г) (б)
В области Q = [о, L]x [0,R ] введем сетку А = Аz хАr, где A z ■
0 = Z0 < Z\ < ... < ZN = L; Ar ■ 0 = Г0 < f\ < ... < M = R. Кубический сплайн двух переменных на такой сетке можно представить в виде [3]
N+1
S(z, r)= X Wi (r)• Bi (z), (3)
i=-1
где
M+1 _
wi(r) = Xbij Bj(r) i = -1,..., N +1. (4)
j=-1
Здесь z - составляющая, представленная в виде дважды непрерывно дифференцируемой функции - кубического нормализованного B -сплайна, а r -составляющая - в виде Wi (r).
B -сплайны определенной степени образуют базис в пространстве сплайнов [3]. Поэтому другие функции пространства могут быть представлены через B -сплайны единственным образом. Здесь B -сплайны нумерованы по среднему узлу их интервалов носителей. Для полного определения базисных функций Bi (z), входящих в (3), Az
дополним узлами z_ 3 < z- 2 < z-1 < z0, zN+3 > zN+2 > zN+1 > zN. Ar дополним
узлами r_ 1 < ^ rM < rM+1.
Приближенное решение уравнения (2) ищем в виде (3). Для простоты рассмотрим равномерную сетку h = zi+1 - zi, t = rj+1 - rj (h, t - const). Зафиксировав r, будем
считать, что изменение температуры по z вблизи узла (zi, rj ) описывается сплайном
N+1
S(^ rj )= X wi, jBi (z), i=-1
где wi j = wi (rj ). С учетом финитности сплайна [3] имеем
99
5 (2/ ,Гу )= Ъ у'в/,(2/). (5)
/ '=/-1
Для равномерной сетки эта формула принимает вид
Э(г/, Гу) = 6 _ 1, у + 3 ^у + ^+1,у. (6)
Согласно идее метода коллокации приближенное решение Э(2, Г) в узлах
коллокации удовлетворяет уравнению (2). Узлы коллокации примем совпадающими с
узлами сплайна. Подставим (3) в уравнение (2) и потребуем совпадения левой и правой частей в узлах (2/, Г у )
P Vz
a(cpS) Э z
+
ij
P Vr
a(cpS)
Э г
ij
' Э f. 3S4
Л- лз^т~ д z v д z
+
ij
1
Г д Г
д L 3S4
1зГ— д Г
+ fjj, (7)
ij
V ^ Л
у = 0,1,..., М; / = 0,1,..., N.
Здесь нижний двойной индекс / у означает, что соответствующий элемент рассматривается в точке (2/, Г у). Вблизи точки (2/, Г у) значения Ср, Л3 меняются незначительно. Поэтому в окрестности этой точки Ср, Л3 можно вынести за знак производной. Тогда уравнение (7) примет вид
ГЪ 5\ ( \ (д
(p cpvz )ij
ij
+(P cp^ )jj
=(1з)
Ґ 2 \ Э 2S
ij
ij
Э z
2
+(1з)
ij
ij
1 Э f Э SA г
Г д Г
v
д г
+ fjj. (8)
ij
л,
Выведем формулы для производных ^ S . На основе (З) при Г = Г j имеем
Э Z -jJ J
I+1 I ) dБj'
,,Z ) d- •
і =1 -1
Учитывая, что для равномерной сетки [З]
Э z*
Э2S —1 / 4 d^.-j
Э z
л
(9)
Б-_ 1(zi ) = - — 1 1V " 2h
Б- (zj ) = 0,
Б/11 (zj ) = -L, Б-^і ) =
h
З
h
З
Б-+І^і ) = —
1 +п и 2h
Б[+1 (zj ) = 4
h
З
при Г = Гj получим
dS д z
z=z,
1 1 Э2S
+^hWj+1, j,
= 4t (Wi-1,j -Wj + Wi+1,j). (10)
z=zi
h
Для аппроксимации производной во втором конвективном члене в (8) используем разности против потока
L
100
Э S Э г
ij
Sij ~ Si, j -1
t
Si, j+1_ Sij
t
если vr > o;
если Vr < 0.
Производная во втором диффузионном члене в (8) аппроксимируется в виде
1
1 А
Г д Г
/.Э5Л Э г
i j j
г і t
І, j+1/З
І, j - І/З.
(ІЗ)
Индекс (у + 1/2) указывает на то, что производные берутся в точках, лежащих посередине между соответствующими узлами сетки. Члены в квадратных скобках в правой части формулы (12) аппроксимируем следующим образом:
/.Э5Л Э г
= О'+1/З'
Si, j+1 _ Sij
A j+1/2
Тогда формула (ІЗ) принимает вид
t
І, j -1/^
= гІ -1/2 ••
Sij ~ Si, j -1
t
1 A
Г д Г
f 'dSЛ Э г
l
■■ г-12 Ij гі t
-V2 'Si, j-1 "Ь-ІЗ + гі+1З)• Si j + гі+1З • Si, j+1].(1З)
Выразим производные во втором конвективном и втором диффузионном членах через Wjу :
3S Э г
ij
= ir(wi-1,j+p ~Wi-lj“1+p)+ Ь^І,j+p "Wi,j“1+p)+ + -t^i+1,j+p -Wi+1,j-1+pI
(14)
где
p
[o, если Vr ^ 0;
[1, если Vr < o;
1 'A
Г д Г
jj гі
.Л
Г і t
гІ - ІЗ • |б Wi -1, j -1 + З Wi, j -1 + б Wi+1, j -1
П-12 + гі+l2)' I б Wi- 1,j + З Wi, j + б Wi+1, j
+
1
З
1
+ O'+ІЗ I б Wi-l,j+1 + ЗWi,j+1 + б Wi+l,j+1
(15)
ч6 ' 1 " 3 1" 6
Чтобы избежать деления на ноль в последней формуле при у = 0, применяется следующий прием. Левую часть формулы (15) запишем в виде
1 А
Г д Г
f dS} 1 dS Э2S
+
Э г
1
1
А
1
101
Вследствие осевой симметрии дЭ/д Г = 0. Поэтому в первом слагаемом получаем
неопределенность (0/0). Применение правила Лопиталя дает д2э/д Г2. Если для аппроксимации второй производной использовать формулу центральной разности, то при у = 0 получим
1 А
Г д Г
dS
д Г
2
Si, -1 _ 2S1, 0 + Si, 1
t
(1б)
Затем, используя формулу (6), правую часть выразим через Wj ].
После подстановки правых частей формул (10), (14), (15) и несложных
преобразований уравнение (8) принимает вид
/ +1 ] +1
Е Е А/ ]' wj у= ^] (У = 0,1,..., М; / = 0,1,..., N), (17)
/ '=/-1 ] '=] -1
где
гі - V2
Vij-l =_(Л®h' '—ї+(иг^/ja/,J; Ai,
4A
A
A
A
j,j -1 І- 1,j -1 І+1,j -1 І- 1,j - І
(p v1cp)// (1з)// (1з)//(г] -1/2 + г] +1/2)
б Г] t
І - 1,j
2h
h
A
ij
2(із);] + 2 (із )jj(г]-ІЗ + г]+1/2)
^ + З
б Г] t
h
.Л
г]є
it2
+ 4(иг )ijbij
+ (иг )ijbjj;
A
i+1, j
(p v1cp )jj (1э )ij (із )ij (гі-ІЗ + г]+1/2) ^ h^+ б/2
A
(із)/j •г]+1/2
+ (иГ) ijbij;
І -1, j+1 . з + V )ijcjj; Ai, / +1 4' Ai -1, j+1; Ai+1, j+1 Ai -1, j+1;
б Г] t
a/j =
(p cp)/ i
IJ
-, если Vr > o; ь;і = < бt І] <
0,
если Vr < o;
(p cp )j i
/J
Pep I
если vr > o;
cij = <
/J
6є
, если vr < o;
. 0 если Vr > o;
(P cp )jj
, если Vr < 0.
6є
Следует заметить, что коэффициенты А/ у при у = 0 вычисляются с учетом (16).
В системе (17) на каждом у-м слое имеется N+3) неизвестных коэффициента wjj . Поэтому необходимо еще по два уравнения на каждом слое. Недостающие уравнения получаем с помощью граничных условий. С помощью граничных условий при
102
2 = 0 и 2 = L исключаем из системы неизвестные W_ 1, у и WN+1, у. Неизвестные Wj _ 1 исключаются с помощью условия симметрии на оси, неизвестные Wj М+1 - с помощью граничных условий на боковой поверхности (г = ).
Рассмотрим вывод дополнительных уравнений относительно wij для некоторых
из границ (остальные рассматриваются аналогично).
1. Входное сечение. На входе температура известна: Э0,у = Твх; 0 < у < увх. На
входные точки (/, у)(/' = 0; 0 < у < увх ) распространим уравнение энергии. Считаем, что в
этих точках выполняется дискретный аналог этого уравнения. При этом полагаем, что ьг = 0.
С учетом (6) условие Э0 у = Твх принимает вид
1 2 1 .
6 w~ 1,у + 3 W0, у + 6 W1,у = Твх, 0 “ у ~ увх . (18)
Отсюда
W-1,у = -4W0,у - W1,у + 6Твх (0 < у < увх). (19)
Подставив (19) в уравнение (17), получим
у+1
где
Z Z A(j Wi'j"=jl) 1 - і - івх -1 І =0 j = і -1
A(1) = 0. A(1) = 0. A(1) =2 (p vicp )oі +б (1э )oі .
A0, j -1 = 0; \j -1 = 0; A0,j - h + h 2 .
a(i)> vicp )o ]. a(i) - o. a(i) = 0
Ai,j = h ’ A0,j+1 =0; Ai,j+1 =0.
/ \
(20)
^ = f0 і б Твх
A i,j-1 + A- 1,j + A- 1,j+1
V ^ ^ ^ -у
Таким образом, в узлах входного сечения имеет место уравнение
ъ а/1) wj I=^1 0 ^у ^ увх. (21)
/=0 у
2. Граница Г1. На жесткой стенке возможны разные варианты граничных условий. Самое простое из них - это граничное условие I рода. В таком случае считается известным распределение температуры вдоль внутренней стороны стенки. Для границы Г1 это выглядит так: Э0 у = Ту (увх — у — М). Проведя почти аналогичные рассуждения, что и
в предыдущем случае, выводя уравнения для Wjу .
Больший практический интерес представляют граничные условия, учитывающие отвод тепла через стенку за счет теплопроводности. Допустим, что внешнюю поверхность стенки окружает среда с заданной температурой Тср и коэффициентом теплоотдачи ос2 (рис 1.б).
103
Как известно, при установившемся тепловом режиме для определения теплопередачи плоской стенки используется формула [2]
q=k (т„6 - Т,р I (22>
где к = 1/ (1 СС[ +S/Aqj + 1/ a2) - коэффициент теплопередачи (а 1 - коэффициент
теплоотдачи от теплоносителя к стенке, 8 - толщина стенки, А,СТ - коэффициент
теплопроводности стенки); Тоб - температура в объеме теплоносителя.
Удельный тепловой поток, передаваемый горячим теплоносителем стенке, определяется формулой [2]
q = а\ (Тоб - Тст\ (23)
где Т ст - температура внутренней поверхности стенки.
Приравнивая правые части формул (22), (23), получим
Тст = ^1тоб + d 2тср- (24)
где d\ = 1 - к/ОС{ , d2 = к/ОС\ . Если в качестве тоб взять температуру на расстоянии
одного шага от стенки, то уравнение (24) принимает вид
S0, j = d1S1 j + d2тср’ Jbx +1 - j - M. (25)
С учетом (6) из этого уравнения следует w_ 1, j = (d1 - 4)w0j + (4d1 - 1)w1, j + d\w2 j + 6d2Tcp, jBx +1 ^ j ^ M. (26)
С помощью этого уравнения исключаются неизвестные W_ 1 j в условиях коллокации.
3. Граница Г2. Рассматривается полностью аналогично предыдущему случаю. Отличие состоит только в том, что при использовании граничного условия III рода для удельной теплопередачи стенки используется формула для цилиндрической стенки.
4. Граница Г3. Рассматривается аналогично границе Г1.
5. Выходное сечение. В узлах выходного сечения ставится условие Э Т/Э Z = 0 (i = N; 0 < j < jBux). С учетом (10) это условие принимает вид
1 1 " 2hWN - 1j + 2hWN+1,j' = °.
Отсюда Wn+i , j = Wn-i , j.
6. Ось симметрии. На оси симметрии ставится условие симметрии ЭТ/ЭГ = 0.
Используя формулу центральной разности, получим Sj _ 1 = Sj \\ i = 0,1, ..., N . С учетом (6) это условие принимает вид Wj _ 1 = Wj 1; i = 0,1, ... , N .
В результате исключения ряда неизвестных с учетом граничных условий в системе (17) останется (N +1)- (M +1) неизвестных Wj, j (i = 0,1,... ,N; j = 0,1,... ,M).
Систему (17) запишем в матричной форме. Для этого сначала упорядочим неизвестные. В соответствии с лексикографическим упорядочиванием [4] неизвестные Wjj
объединяются в вектор
W = (W00> ^..., WN0, W0L W11..., WN1,..., W0M, W1M,..., WNM) , где ' означает транспонирование. Тогда система (17) примет вид
A - W = g , (27)
104
где А - матрица из N + і)(М + 1) строк и N + і)- (М + і) столбцов; д - вектор из правых частей.
Матрица А является разреженной. Ее диагональные элементы отличны от нуля. Структура матрицы показана на рис. 2.
Рис. 2 - Структура матрицы А (+ - ненулевые элементы)
Для решения системы уравнений применяется метод Гаусса с частичным выбором главного элемента по столбцу.
Для решения всей задачи в целом применяется следующий итерационный алгоритм:
1. Задание размеров камеры, расчетной сетки. Ввод скоростного поля.
2. Задание начального приближения для температурного поля.
3. Расчет теплофизических свойств при известном поле температуры.
4. Решение системы линейных уравнений (27). Расчет нового температурного поля.
5. Проверка условия
где к - номер итерации; е - заданная точность. При выполнении условия (28) итерационный процесс прекращается, в противном случае осуществляется возврат к пункту 3.
Вычислительная схема реализована в виде программы на Фортране (БРБ 4.0).
Для проверки работоспособности и эффективности алгоритма произведены тестовые расчеты и сравнение с результатами [5], полученными методом конечных разностей. Поле скоростей было рассчитано с помощью программы из [6] при тех же исходных данных, что и в работе [5]. Проведено также сравнение полученных расчетных данных с экспериментальными данными других авторов [7]. Во всех расчетах теплофизические свойства считались переменными в рассматриваемой области. Для их вычисления использовались те же зависимости, что и в [5].
105
В [5] произведен расчет турбулентной неизотермической струи в цилиндрическом канале при тех же условиях, что и в [7]. В работе [7] в качестве рабочего тела использовался нагретый воздух, подаваемый через сопло сверху (по направлению к оси
2 ). Размеры камеры следующие: радиус канала 0.315 м; радиус входного сечения 0.025 м; высота канала 3.25 м.
В уравнении (8) источниковый член в этом случае равен нулю. Турбулентное число Прандтля было взято равным 0.7. Коэффициент теплоотдачи ос (от нагретого воздуха к стальной стене) определяется из уравнения подобия [8].
На рис. 3 показано изменение температуры на оси струи при значении скорости на входе = 4 м/с. В расчетах для достижения точности £ = 0.01 (в условии (28)) потребовалось девять итераций (при N = 20, М = 10).
Г, К
380 -
340 -
300 4 0
Рис. 3 - Изменение температуры на оси струи (• — эксперимент [7]; — — — расчет [5]; — — метод сплайн-коллокации)
При простоте вычислительной схемы полученные результаты качественно и
количественно согласуются с имеющимися расчетными и экспериментальными данными.
Литература
1. АдриановВ.Н. Основы радиационного и сложного теплообмена. - М.: Энегрия, 1972. - 464 с.
2. Шорин С.Н. Теплопередача. - М.: Высшая школа, 1964. - 490 с.
3. ЗавьяловЮ.С, КвасовБ.И, МирошниченкоВ.Л. Методы сплайн-функций. - М.: Наука, 1980. - 352 с.
4. Ортега Дж, Пул У. Введение в численные методы решения дифференциальных уравнений / Пер. с англ./ Под ред. А.А. Абрамова. - М.: Наука, 1986. - 288 с.
5. Вафин ДБ, Садыков А. В. Некоторые результаты численного исследования аэродинамики топочных устройств / Казан. хим.-технол. ин-т. - Казань, 1988. - 14 с. Деп. в ОНИИТЭХИМ Черкассы 20.07.88. № 722-ХП 88.
6. Садыков А.В. Разработка численного метода расчета топочных камер трубчатых печей: Дис. ... канд. техн. наук/Казань, 1989. - 169 с.
7. Глебов Г. А. Козлов А.П. Турбулентная струя в канале при воздействии архимедовых сил // ИФЖ. - 1988. - Т 55. - № 2. - С. 191-198.
8. ЮдаевБН Теплопередача. - М.: Высшая школа, 1973. - 360 с.
106
© А. В. Садыков — канд. техн. наук, доц., декан факультета управления и автоматизации НХТИ КГТУ; А. А. Ахметов - студ. НХТИ КГТУ.
107
