CC BY
50
ГИДРОДИНАМИКА, ТЕПЛО- И МАССООБМЕННЫЕ ПРОЦЕССЫ, ЭНЕРГЕТИКА
УДК 536.3:535.34 А.В. Садыков
К СОВМЕСТНОМУ РЕШЕНИЮ УРАВНЕНИЙ ЭНЕРГИИ И ПЕРЕНОСА ИЗЛУЧЕНИЯ
Ключевые слова: теплообмен, излучение, горение, трубчатые печи.
Рассматривается вычислительная схема для совместного численного решения уравнений энергии, переноса излучения и уравнений модели горения. Для получения дискретных аналогов уравнений энергии и модели горения используется метод сплайн-коллокации в сочетании с методом конечных разностей.Радиационный
теплообмен рассматривается в рамках р - приближения метода сферических гармоник.
Keywords: heat exchange, radiation, combustion, tube furnace.
A computational scheme for joint numerical solution of equations for energy, energy radiative transference and equations for combustion model are considered. In order to obtain the discrete analogs of equations for energy and a model of combustion, use is made of the spline-collocation method in combination with the finite-difference method. Radiation heat is considered in the frames of P- approximation of spherical harmonics method.
Тепловой расчет топочных камер технологических трубчатых печей нефтехимической и нефтеперерабатывающей промышленности сводится к решению задачи радиационно-
конвективного теплообмена в условиях горения топлива. Комплекс процессов радиационно-
конвективного теплообмена описывается системой нелинейных дифференциальных и интегро-
дифференциальных уравнений.
Метод, предложенный в [1] для решения уравнения энергии, является простым в
реализацзии. Поэтому представляет интерес возможность применения этого метода при решении задачи радиационно-конвективного теплообмена.
В данной работе основное внимание уделено совместному решению уравнений энергии, переноса излучения, модели горения. Решается стационарная задача для двумерной цилиндрической геометрии (рис.1). При этом поле скоростей считается известным.
Уравнение энергии для рассматриваемой геометрии с учетом уравнения неразрывности потока имеет вид
д(СрТ) д(СрТ)
Pvz—I------+ pvr
дz
дг
д_Ґ
дz
і z дТ
1 э я
дz ,
+
1 _д г дг
і эгг f,
+f
(1)
где \ = Чу - , Уг - компоненты вектора
скорости V ; с - изобарная теплоемкость;, Л|, ЛЭ
— коэффициенты эффективной теплопроводности по осевому и радиальному направлениям; Т = Т(ъ,г)
- температура в точке с координатами ъ, г; -
плотность источников тепловыделений в объеме
топочной камеры; Чр - вектор плотности
радиационного потока энергии. Для нахождения
распределения в объеме топочной камеры
решаются дифференциальные уравнения двухшаговой модели горения. Дивергенция лучистых потоков определяется в результате решения уравнения переноса излучения.
Уравнение (1) дополняется рядом граничных условий. На входном участке в топку ставится граничное условие I рода. На оси потока(Ог) задается условие симметрии. На выходном участке задается нулевой градиент температуры. На жесткой границе ставятся граничные условия I или III рода.
Стационарное уравнение переноса излучения для конечного цилиндра имеет вид
8ІП0
С08ф-
аМ_
г
+
co89
о+
8ІПФ
IЛ М
Ф
k* I
,(м,ф
(2)
kp
4п
2л п
0
0
+кпЛ ь (т)
Здесь г, г пространственные координаты точки М; П - направление луча с угловыми
координатами 9 и Ф (рис.1); ку , кд, к = ку + кд -
спектральные коэффициенты поглощения,
рассеяния, ослабления соответственно ; I ь (т) -
функция Планка; 1А (м , П) -спектральная
интенсивность излучения.
r
I
z
Уравнение переноса излучения дополним граничными условиями. На оси симметрии ставится условие симметрии решения относительно оси ъ
1х (м ■ А) r
при г = 0.
Граничные условия для интенсивности излучения 1Л (м, п) на тепловоспринимающей поверхности, на футеровке свода и пода имеют вид
1л(м,ф,е) = £, • 1Ь() + Пдп,
при Пп<0 на Г . Здесь £ э, г э - полусферические излучательная и отражательная способности стенок
Рис. 1- Упрощенная схема топочной камеры и система координат. Обозначения: 1 - боковая стенка;2- трубчатый экран; 3 - горелка.
соответственно; п - внешняя нормаль к границе; Т^ - температура поверхности стенки со стороны
топки; Г -граница области; Цп - поверхностная
плотность теплового излучения в направлении внешней нормали:
ЧП = ^ !Л (м, )°,пс10, .
о’П>о у 7
В работе использована двухшаговая модель горения
СН4 +1,5 О2 ^ СО +2 Н2О ,
СО +0,5 О2 ^ СО2 .
При этом скорости реакций определяются по модели распада вихрей. В качестве окислителя используется воздух, в качестве газообразного топлива - природный газ, 98% которого составляет метан.
Топочные камеры трубчатых печей
работают при относительно низких давлениях
(порядка 1 атм.) и высоких температурах (1000^2000 К). Поэтому состояние дымовых газов рассматривается в приближении идеального газа.
Для описания турбулентного характера потока использовалась гипотеза пути смешения Прандтля. Коэффициенты турбулентной диффузии, теплопроводности и вязкости выражаются через турбулентные критерии подобия. Зависимости
коэффициентов динамической вязкости и
молекулярной теплопроводности от температуры выражаются формулами Сутерленда [2]. Удельные теплоемкости компонентов смеси при постоянном давлении задаются в виде эмпирических
зависимостей от температуры. Теплоемкость
газовой смеси определяется по данным для
компонентов по правилу аддитивности. Эффективные коэффициенты переноса равны сумме турбулентной и молекулярной составляющих.
Для решения уравнения энергии и
уравнений модели горения применяется метод сплайн-коллокации в сочетании с методом
конечных разностей. Рассмотрим сначала уравнение энергии.
В области [0, ь] х [0, р введем сетку А = Аъ х Аг, где А2 : 0 = ъ0 < ъ1 < ... < = I. ;
Дг : 0 = Го <г <... <м . Кубический сплайн
двух переменных на такой сетке можно представить в виде [3]
8 (ъ,г )Л%(г) х Б, (ъ), (3)
где
M+1
W:
(r)= I bjj Bj (r); i = -1,..., N + 1. (4)
j=-1
Здесь z -составляющая представлена в виде дважды непрерывно дифференцируемой функции -кубического нормализованного В-сплайна, а Г -
составляющая - в виде Wj (г ) .
В-сплайны определенной степени образуют базис в пространстве сплайнов [3]. Поэтому другие функции пространства могут быть представлены
через В-сплайны единственным образом. Здесь B -сплайны нумерованы по среднему узлу их интервалов носителей. Для полного определения
базисных функций Bj (z), входящих в (3), Az дополним узлами z_g <Z-2 <Z-1 < Zg ,
ZN+3 > ZN+2 > ZN+1 > ZN,а Ar - узлами r-i < rg, rM < rM+1.
Приближенное решение уравнения (1) ищем в виде (3). Для простоты рассмотрим равномерную сетку
h = Zi+1 -Zi, т = rj+1 - rj (h, т - const). Зафиксировав г , будем считать, что изменение температуры по z вблизи узла (Zj,rj) описывается
r
сплайном
Б ^) %2Лі ві (г),
где Ші j = Ші (г|). С учетом финитности сплайна
имеем
8 (ъ^ (5)
Для равномерной сетки эта формула принимает вид
(6)
+7'*«+^м,г
Узлы коллокации возьмем совпадающими с узлами сплайна. Подставим (3) в уравнение (1) и потребуем совпадения левой и правой частей
уравнения в узлах ^>г^. Для аппроксимации
производной во втором конвективном члене в (1) используется разность «против потока». Производная во втором диффузионном члене в (1) аппроксимируется в виде
1 .А
г д г
(
дБ
г—
д г
Ч У
1
г/т
( дБ ^ г —
. д г ,
Ч У
, /'+1/2
( дБ ^ г —
д г Ч У
, ./-1/2.
Индекс (] ± 12) указывает, что производные
берутся в точках, лежащих посередине между соответствующими узлами сетки. Члены в квадратных скобках в правой части последней формулы аппроксимируем следующим образом:
( д8Л г—
д г ч у
8
= г;
і,І+1/2
'І+12
У+1 - 8ІІ
^ д8Л г—
дг
ч У
=г
І-12
В результате простых преобразований получаем
1+1 |+1
£ £ А1| ^7 =^ (7)
1'=1-1]'=1-1
(| = 0,1,..., М; 1 = 0,1,..., N.
Выражения для коэффициентов приведены в
[4].
В системе (7) на каждом у-м слое имеется
(ы + э) неизвестных коэффициента .
Недостающие уравнения получаем с помощью граничных условий для уравнения энергии. С помощью граничных условий при ъ = 0 и ъ = .
исключаем из системы неизвестные Ш_1 j и
ш^+1 j. Неизвестные wj-1 исключаются с помощью условия симметрии на оси, а неизвестные
м+1 - с помощью граничных условий на боковой
поверхности (г = Р) .
Вывод дополнительных уравнений
относительно Wjj для разных границ для
переменной Т подробно рассмотрен в [4].
В результате исключения ряда неизвестных с учетом граничных условий в системе (7) останется
(ы + 1)(м + 1) неизвестных ш иС = 0,1...,м;
, = 0,1,..., Ы). Систему (7) запишем в виде
В0Ф0 + С0Ф1 =Ро
AjФ j—і + BjФ j + CjФ j+1 =Р j
(8)
АМфМ-1 + Вмфм =р
м
0 = 1,..., м -1)
где Aj,Bj,Cj - трехдиагональные матрицы; Фj = (Wоj,W1j,...,WNj)'; Р. - вектор из правых
частей. Здесь ( ) - означает транспонирование.
Получение дискретных аналогов для уравнений модели горения на основе вышеизложенной методики полностью аналогично предыдущему. Вывод дополнительных уравнений
для Ч>1/- на основе граничных условий для
У
параметров модели горения аналогичен выводу таких уравнений для Т. Для параметров модели горения получаем аналогичные системы уравнений вида (8).
Для решения полученных систем применен блочный метод последовательной верхней релаксации [5]. Алгоритм метода реализуется следующим образом
Ф
(к+1) = 0
= (1-т).Ф0к) + т В0'(-С0Ф^ + ^)
-1
ь(к).
Ф(к+1) = (1-т). Ф(к) + т В-1^^-+1) - ^Ф(к + ^)
Ф
(к+1) = (
М = (1-т).ФМ + т Вм1(-АМФМ-Г + >ЙМ)
'П-1
ь(к+1).
(j = 1,...,м -1)
где Т - итерационный параметр, к - номер
итерации. Матрицы Bj имеют диагональное
преобладание. Поэтому при реализации этого алгоритма не возникает проблем с обращением матриц.
Дивергенция лучистых потоков находится путем решения уравнения переноса излучения (2). Для решения уравнения переноса излучения использовано Р1-приближение метода сферических гармоник (МСГ). Из системы трех
дифференциальных уравнений Р1-приближения МСГ для первых трех моментов сферических функций можно получить уравнение диффузионного типа для нулевого момента ф :
т
т
_З_|0 ^Ф00 Зъ I Зъ
+ "гI(Г°ЗфЗ0° ^кПФ00 =-кПIb(T),
где Б = 1/(3(к - g1kp)) - коэффициент радиационной диффузии; g1 - первый коэффициент разложения индикатрисы рассеяния по полиномам Лежандра. Это уравнение решается методом конечных разностей [6].
Граничные поверхности считаются диффузно излучающими и отражающими. При этом ряд труб заменяется сплошной непрозрачной для теплового излучения лучевоспринимающей поверхностью с эффективной степенью черноты. Вопросы практической реализации Р1 - приближения МСГ для рассматриваемой геометрии, а также стыковки с уравнениями энергии, движения, модели турбулентности, модели горения подробно рассмотрены в [6].
Для вычисления интегрального по спектру коэффициента поглощения газообразных продуктов сгорания использована модель взвешенной суммы серых газов [7].
Дискретные аналоги дифференциальных уравнений энергии и модели горения образуют системы взаимосвязанных алгебраических уравнений. Входящие в них коэффициенты и источниковые члены зависят от ряда переменных. Теплофизические параметры, входящие в уравнения, являются функциями от температуры. Температура в свою очередь является переменной величиной. Поэтому для решения задачи применяется итерационный алгоритм.
Итерационный процесс организуется следующим образом:
1. Задаются исходные данные: температура газов на входе, массовые концентрации горючего и окислителя на входе, температура поверхности нагрева, степени черноты ограждающих поверхностей, параметры футеровки и др.
2. Задаются поле скоростей, начальные приближения полей температуры, массовых концентраций горючего и окислителя в объеме топочной камеры.
3. Расчет поля теплофизических параметров при известном поле температуры, коэффициента турбулентной теплопроводности.
4.Вычисляются интегральный коэффициент поглощения продуктов сгорания и другие параметры, необходимые для решения уравнения переноса излучения.
5. Решается уравнение переноса излучения при известном поле температуры, определяются лучистые тепловые потоки к ограждающим топочную камеру поверхностям и дивергенция лучистых потоков во всех точках расчетной сетки.
6. Решаются уравнения модели горения.
7.Решается уравнение энергии при известных значениях дивергенции лучистых потоков, массовых концентраций горючего и
окислителя. В результате определяется поле осредненной температуры.
8.Полученное поле осредненной температуры принимается как начальное и действия, начиная с 3-го пункта, повторяются.
Итерации заканчиваются при выполнении критерия сходимости
-(к +1)
тах
У
1 -
1,1
Т.
(к)
1,|
итерации; £ -
< е,
заданная точность,
где к - номер задается в пределах 0,001-0,01.
В результате определяются поле осредненных значений температуры и распределения плотностей лучистых и
конвективных тепловых потоков к поверхности нагрева и футеровке.
В ходе итерационного процесса совместного решения уравнений энергии и переноса излучения в первых внешних итерациях поле температуры сильно осциллирует. Для подавления неустойчивости применялась нижняя релаксация
Чр(к) = а • Чр(к) + (1 - а)*Чр(к-1), где к - номер внешней итерации, а - параметр релаксации.
Вычислительная схема реализована в виде программы на Фортране. Для проверки работоспособности метода проведены тестовые расчеты. Расчеты проводились для печи следующих размеров: диаметр - й = 1,68м (диаметр сечения, проходящего по осям труб, расположенных вертикально вдоль футеровки); высота - Ь = 9,6м;
диаметр туннеля сводовой горелки - йТ = 0,5м;
диаметр канала для отвода дымовых газов йвых =
0,9м.
Эффективная степень черноты и отражательная способность трубчатого экрана вычисляются по методике, предложенной в [8].
Горючая смесь состоит из природного газа (содержание метана 98%) и воздуха. Температура частично сгоревшей в туннеле горелки газовоздушной смеси на входе в топку Т=1270К, а скорость на входе уъ,0=10м/с. Поле скоростей
заимствовано из [6].
Для уравнения энергии на жесткой границе использовано граничное условие I рода. На поверхности нагрева задавалось характерное для таких процессов распределение температуры (на рис.2-нижняя кривая). В расчетах полагалось, что радиационные свойства поверхности нагрева постоянны.
Некоторые из полученных результатов приведены на рис.2, 3.
Рис. 2 - Профиль температуры продуктов
сгорания на оси топки. Обозначения: ---
^эф ~ 1; — " -; ^эф ~ 0,8 ; — — — ^эф ~ 0,5 .
Нижняя кривая - температура поверхности нагрева.
чГ
80
60
40
20
0
кВт /л 42
> \
' \ ^ \ N X N Л N \ ч\
\\ч
NV
0 0,25 0,5 0,75 5 / Ь 1
Рис. 3 - Распределение поверхностной плотности
рез
радиационного потока тепла по
поверхности нагрева. Обозначения как на рис. 2.
На рис.2 показано влияние эффективной степени черноты поверхности нагрева на распределение температур. Уменьшение эффективной степени черноты поверхности нагрева приводит к повышению уровня температуры в
топке. Это связано с тем, что увеличивается доля отраженного излучения. В результате этого повышается внутренняя энергия продуктов сгорания. С увеличением эффективной степени черноты трубчатого экрана возрастают плотности
лучистых потоков к ним (рис.3). Увеличение
приводит также к неравномерности обогрева труб по их длине.
Полученные результаты согласуются с данными работы [6, 9] и правильно отражают топочные процессы. Проведенные исследования показывают, что предложенный метод можно использовать при тепловых расчетах топок трубчатых печей цилиндрического типа.
Литература
1. Садыков, А.В. Решение двумерного уравнения энергии методом сплайн-коллокации / А.В.Садыков // Сеточные методы для краевых задач и приложения. Материалы Шестого Всероссийского семинара.-Казань: Казан. гос. университет, 2005. - С.194-197.
2. Иссерлин, А. С. Основы сжигания газового топлива / А.С.Иссерлин.-Л.:Недра,1987.-336с.
3. Завьялов, Ю.С. Методы сплайн-функций / Ю.С. Завьялов, Б.И. Квасов, В.Л. Мирошниченко.- М.: Наука, 1980. - 352 с.
4. Садыков, А.В. Некоторые результаты численного решения уравнения энергии методом сплайн-коллокации / А.В. Садыков, А.А. Ахметов // Вестник Казанского технологического университета / КГТУ.-2005.- №2,ч.11 - С.98-106.
5. 5. Марчук, Г.И. Методы вычислительной математики / Г.И. Марчук - М.: Наука, 1980. - 534 с.
6. Садыков, А.В. Разработка численного метода расчета
топочных камер трубчатых печей: Дис....
канд.техн.наук.- Казань: КХТИ, 1989. -169 с.
7. Smith, T. F. Evaluation of Coefficients for the Weighted Sum of Gray Gases Model/ T.F. Smith, Z.F. Shen , J.N. Friedman // J. Heat Transfer.- 1982. -№104.- P. 602-608.
8. Гориславец, С.П. Пиролиз углеводородного сырья / С.П. Гориславец. - Киев: Наукова думка, 1977. - 305 с.
9. Садыков, А.В. Влияние степени черноты поверхности нагрева на теплообмен в трубчатых печах цилиндрического типа / А.В.Садыков, И.М.Валеев, Д.Б. Вафин // Вестник Казанского технологического университета. - 2012. Т.15, №7 - С. 50 - 54.
© А. В.Садыков - канд. техн. наук, доцент, декан факультета управления и автоматизации НХТИ ФГБОУ ВПО «КНИТУ», sadykov@land.ru.
