CC BY
30
ГИДРОДИНАМИКА, ТЕПЛО- И МАССООБМЕННЫЕ ПРОЦЕССЫ,
ЭНЕРГЕТИКА
УДК 536.2:517.958:66
А. В. Сады ков, И. М. Валеев ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЯ ЭНЕРГИИ С УЧЕТОМ ГОРЕНИЯ
ГАЗООБРАЗНОГО ТОПЛИВА
Ключевые слова: дифференциальные уравнения в частных производных второго порядка, сплайн-коллокация, стационарная задача, differential equations in the partial derivatives of the second order, spline- collocation, the stationary task.
Предложена вычислительная схема для совместного численного решения уравнения энергии и уравнений модели горения. Для получения дискретных аналогов этих уравнений используется метод сплайн-коллокации в сочетании с методом конечных разностей. n this paper is given a calculation scheme for the numerical solution of the energy equation as well as the equations for combustion models. The discrete identical values are calculated using a combination of spline collocation and terminal difference methods.
Тепловой расчет топочных камер трубчатых печей сводится к совместному решению уравнений переноса излучения, осредненных уравнений энергии, неразрывности, модели турбулентности, модели горения. Система дифференциальных уравнений нелинейна, поэтому может быть решена численными или приближенными методами.
В настоящей работе основное внимание уделено численному интегрированию уравнения энергии с учетом горения газообразного топлива. Решается стационарная задача для двумерной цилиндрической геометрии.
Горение рассматривается в рамках модели простой химической реакции. Согласно этой модели реакция протекает в одну стадию [1 ]
1 кг (горючего) + A кг (окислителя) ^ (1 + A) кг (продуктов сгорания) + Q, где Q - теплота сгорания, A - стехиометрическое число. В любой точке объема имеет место соотношение Шг + ток + Шпр = 1, где тг ,ток, тпр - массовые концентрации
горючего, окислителя, продуктов сгорания соответственно. В рамках этой модели горение предварительно перемешанной газовой смеси описывается уравнениями для Шг, (Рг = Шг - Шок/A и уравнением энергии. Уравнение для этих параметров имеет вид [1]
d {pvz0)+ 1 -d(pvrr0) = -d(D3z ^1 + 1 D3rr ^1 + Эф, (1)
дг г дг дг ^ дг) г дг ^ дг,
7 №
где Ф = (тг ,уг}; йэ ,йэ - коэффициенты эффективной диффузии по осевому и радиальному направлениям, v2 , vг - компоненты вектора скорости v; р- плотность.
Источниковый член в уравнении (1) для тг определяется осредненной по времени моделью Аррениуса
Эг = Сг р2тгток ехр(- Е/НГ) (2)
где Ог, Е - постоянные Аррениуса; р - осредненное значение давления; Т = Т(г, г) -температура в точке с координатами г, Г; Я - универсальная газовая постоянная. Для метано-воздушных смесей Е/Я = 20000К [2], а 03 = 0.5 [3]. В уравнении для (р5
источниковый член равен нулю.
С учетом уравнения неразрывности уравнение (1) принимает вид
дФ дФ д (_ г дФ^1 1 д (^ г дФ
Руг^ + — = -\йэ2 — 1 +
дг дг
Уравнение энергии для неразрывности имеет вид [4]
йэ‘ Г
дг ^ дг ) г дг ^ дг
рассматриваемой геометрии
| + . с учетом
д(СрТ) д(СрТ) д
"г-еТ- + Р'>гЧрГ = дг
г
г
дг
т г дТ
+— і г г—
г дг { э дг
+ і,
(3)
уравнения
(4)
где ^ = Цу - divqр; qv - плотность источников тепловыделений в объеме топочной камеры; div Цр - определяется путем решения уравнения переноса энергии излучения;
г г
Ср - изобарная теплоемкость; Лэ, Лэ — коэффициенты эффективной теплопроводности
по осевому и радиальному направлениям. Целью данной работы является совместное численное решение уравнений энергии и модели горения. Поэтому поле скоростей считается известным, а излучением пренебрегается.
Источниковый член (Цу) в уравнении (4) выражает выделение теплоты в результате экзотермической реакции. При сгорании природного газа Цу рассчитывается по формуле [2]
Цу = Цотгток1Т ехр(- Е/НТ), где Цо - определяется из общего теплового баланса топки.
Топочные камеры трубчатых печей работают при относительно низких давлениях (порядка 1 атм.) и высоких температурах (1000^1900 К), поэтому состояние дымовых газов может быть описано в приближении идеального газа р = рТ/М, где М -молекулярная масса дымовых газов.
Уравнения дополним граничными условиями. На входе в топку задаются граничные условия I рода (рис. 1). На оси симметрии (02) для всех параметров задается условие симметрии, а на выходе - условие нулевого градиента. Краевые условия на жесткой стенке отличаются для разных переменных с учетом их физического переноса. Для уравнения энергии ставятся граничные условия I или III рода. Граничные условия для тг, Рг основаны на предположении отсутствия потока для них в нормальном направлении.
Рис. 1 - Размеры топки, границы
Для описания турбулентного характера потока использовалась гипотеза пути смешения Прандтля. Коэффициенты турбулентной диффузии, теплопроводности и вязкости выражаются через турбулентные критерии подобия. Зависимости коэффициентов динамической вязкости и молекулярной теплопроводности от температуры выражаются формулами Сутерленда [2]. Удельные теплоемкости компонентов смеси при постоянном давлении задаются в виде эмпирических зависимостей от температуры. Теплоемкость газовой смеси определяется по данным для компонентов по правилу аддитивности. Эффективные коэффициенты переноса равны сумме турбулентной и молекулярной составляющих.
Для решения уравнений (3), (4) применяется метод сплайн-коллокации в сочетании с методом конечных разностей. Рассмотрим сначала уравнение энергии.
В области Q = [0, L] х [0, R] введем сетку Д = Д z хД г, где Д z : 0 = zo < z\ < ... < zn = L; Дr : 0 = /q < r\ < ... < Гм = R. Кубический сплайн двух переменных на такой сетке можно представить в виде [5]
N+1
S(z,r )= Z Wi(r )Bi(z), (5)
i=-1
где
M+1
wi(r)= Zbij Bj(r) i = -1,...,N +1. (6)
j=-1
Здесь z - составляющая представлена в виде дважды непрерывно дифференцируемой функции - кубического нормализованного В-сплайна, а Г -составляющая - в виде Wi (r).
В-сплайны определенной степени образуют базис в пространстве сплайнов [5]. Поэтому другие функции пространства могут быть представлены через В-сплайны единственным образом. Здесь В-сплайны нумерованы по среднему узлу их интервалов носителей. Для полного определения базисных функций Bi (z), входящих в (5), Дz
дополним узлами z-3 < z-2 < z-1 < ^, zN+3 > zN+2 > zN+1 > zN . Д r дополним узлами Г-1 < ro, Гм < Гм+1.
Приближенное решение уравнения (4) ищем в виде (5). Для простоты рассмотрим равномерную сетку h = zj+1 - zj, т = rj+1 - rj (h, т - const). Зафиксировав Г, будем
считать, что изменение температуры по z вблизи узла (z\, Г j ) описывается сплайном
N+1
Z Wi,j Bi(z), i=-1
где Wi j = Wi (rj ). С учетом финитности сплайна имеем
i +1
S(zi,rj)= Z Wi'jBi,(zi). (7)
i'=i-1
Для равномерной сетки эта формула принимает вид
s(z(, rj )= 16W'-1,j + 3Wij + lWi+1,j . (8)
S(z,rj ) =
Согласно идее метода коллокации, приближенное решение 8(г, г) в узлах коллокации удовлетворяет уравнению (4). Узлы коллокации возьмем совпадающими с узлами сплайна. Подставим (5) в уравнение (4) и потребуем совпадения левой и правой частей в узлах ^ Г-)
рУг
д(ср5) д г
+
/У
РУг
фр^) д г
/У
д
дг
г
дв д г
+
/У
1 д_
г дг
ггг
дв
дг
(9)
/У
у = 0,1, ...,М; / = 0,1, ...,Ы.
Здесь нижний двойной индекс / у означает, что соответствующий элемент рассматривается в точке (г/, гу). Вблизи точки (г/, гу) значения Ср, Лэ меняются незначительно. Тогда уравнение (9) примет вид
д г
)/У
-Р г >Ц
(д в ^
~дг
( 2Л д2в
.2
+ (%г)..
)/У
/У
1 д
г д г
+Цу . (10)
/у
= г ,
;/у V ід г
Для аппроксимации производной во втором конвективном члене в (10) используется разность против потока. Производная во втором диффузионном члене в (10) аппроксимируется в виде
1 _д_ г д г
(
дв
V
г
V д г )
1
. . г; т /У У
/,У+12
( дв ^ г
дг
) І,ї-12
Индекс ( ± 1/2) указывает на то, что производные берутся в точках, лежащих посередине между соответствующими узлами сетки. Члены в квадратных скобках в правой части формулы (12) аппроксимируем следующим образом:
( двл г—
дг )
в/
/,]+12
г]+12
і, У+1 - в/У
( двл г
дг )
/,у-12
г]-12
віу - ві,у-1
После простых преобразований уравнение (10) принимает вид
/ +1 У+1
І I А/’]’ *■,'}'= Гу (] = 0,1..........М; / = 0,1..........N1
/'=/-1 у'=]-1
(11)
Выражения для коэффициентов А/у приведены в [4].
В системе (11) на каждому-м слое имеется (Ы + 3) неизвестных коэффициента ш/у.
Поэтому необходимо еще по два уравнения на каждом слое. Недостающие уравнения получаем с помощью граничных условий. С помощью граничных условий при г = 0 и г = L исключаем из системы неизвестные ш_1 у и +1 у. Неизвестные -1
исключаются с помощью условия симметрии на оси. Неизвестные м+1 исключаются с
помощью граничных условий на боковой поверхности (г = Н).
Вывод дополнительных уравнений относительно ш/у для разных границ для Т
подробно рассмотрен в [4].
В результате исключения ряда неизвестных с учетом граничных условий в системе (11) останется (/V + 1)М +1) неизвестных у (у = 0,1,М; / = 0,1,..., N). Систему (11)
т
т
запишем в матричной форме. Используя лексикографическое упорядочивание, неизвестные ш/у объединим в вектор ш = (ш00,ш10,...,WN0;W01,W11,..,WN1;...;W0M,
ЩМ ,...,WNM )', где ' - означает транспонирование. Тогда система примет вид
А • ш = д, (12)
где А - матрица из (V + 1)М +1) строк и (V + 1 )(М +1) столбцов; д - вектор из правых
частей. Матрица А является блочно-трехдиагональной.
Получение дискретных аналогов для уравнений (3), (4) на основе вышеизложенной методики полностью аналогично предыдущему. Решение также ищется в виде (5). В
7 Г
предположении, что вблизи точки (г/, гу) значения Оэ ,йэ меняются незначительно, имеем
(РУг)у •
(д2в ) дг 2
+
(°э)./у •
1 _^( дв г дг і дг
+
У
(вФ ) у
Это уравнение по виду аналогично уравнению (10). Дальнейший ход рассуждений такой же, что и для уравнения энергии.
Вывод дополнительных уравнений для Wij на основе граничных условий для
параметров тг, <р5 аналогичен выводу таких уравнений для Т. С учетом всех граничных условий для тг, срг окончательно получаем две системы линейных уравнений вида (12). Система (12) является блочно-трехдиагональной и имеет вид В0Ф 0 + С0Ф1 =Г 0
А у Ф у-1 + В у Ф у + С у Ф у+1 =Г у ( у = 1,---, М -1) АМ ФМ-1 + ВМ ФМ =гм
где Ау ,В у ,С у - известные трехдиагональные матрицы; Ф у = (W 0 у ,Wl у ,...^^У; Гі -
вектор из правых частей. Для решения таких систем имеются как прямые, так и итерационные методы. Как показывает вычислительная практика, итерационные методы более предпочтительны. Поэтому она решается блочным методом последовательной верхней релаксации [6]. Алгоритм метода реализуется следующим образом
Ф(0к+1) = (1 - т) • Ф(0к) + т • В-1(-С0Ф(к) + Г) )
Ф(ук+1) = (1 -т) Ф(ук) + т- В -1( - А у Ф(к+1) - Су Ф(к+1 + Гу) (у = 1,.., М -1)
Ф/М +1) =(1 -т) • Фм) +т • ВМ(-Ам Фм-1) + ГМ )
где т - итерационный параметр (те (1; 2)); к - номер итерации. Матрицы В у имеют
диагональное преобладание. Поэтому при реализации этого алгоритма не возникает проблем с обращением матриц.
Дискретные аналоги дифференциальных уравнений энергии и модели горения образуют системы взаимосвязанных совместных алгебраических уравнений. Входящие в них коэффициенты и источниковые члены зависят от ряда переменных. Теплофизические параметры, входящие в уравнения, являются функциями от температуры. Поэтому для решения задачи применяется итерационный алгоритм. Численное решение состоит в циклическом повторении следующих шагов:
1. задание начальных приближений температуры, массовых концентраций горючего и окислителя;
2. расчет теплофизических свойств при известном поле температуры;
3. нахождение следующего приближения поля температуры;
4. нахождение следующего приближения полей массовых концентраций горючего, окислителя;
5. использование новых значений полей температуры и концентраций как исходных и повторение всей процедуры, начиная с шага 2.
Вся процедура выполняется до тех пор, пока не выполнится критерий сходимости.
Вычислительная схема реализована в виде программы на Фортране (CVF V6.5). Такой выбор языка обусловлен тем, что при программировании сложных вычислительных задач для современных ЭВМ лидирующее положение занимает язык Фортран [7].
Для проверки работоспособности алгоритма проведены тестовые расчеты. Расчеты проводились для печи следующих размеров: диаметр D=1.68 м; высота L=9.6 м; диаметр туннеля горелки Dex=0.5 м; диаметр канала для отвода дымовых газов Dewx=0.9 м. Температура частично сгоревшей в туннеле горелки газовоздушной смеси на входе в топку Т3о=1270 К. На входном сечении задается скорость истечения газов, определенная исходя из расхода топлива, коэффициента избытка воздуха и других исходных данных.
Горючая смесь состоит из природного газа (содержание метана 98%) и воздуха. Скорость потока на входе vz,o=10 м/с, а поле скоростей такое же, что и в [8].
Температура боковой стенки (r = R) считается известной. Задавалось распределение температуры, характерное для таких процессов.
На рис. 2 показано изменение температуры потока по высоте топки для различных продольных сечений, а на рис. 3 - изменение температуры поперек потока в различных сечениях топки.
Рис. 2 - Профили температуры продуктов сгорания в различных продольных сечениях топки
Рис. 3 - Профили температуры продуктов сгорания в различных поперечных сечениях топки
Из рис. 2 видно, что повышение температуры на оси потока характеризует нагрев и воспламенение смеси, вызванное перемешиванием в камере высокотемпературных продуктов сгорания со свежей смесью. Как видно из рис. 3, распределение температур в области факела неравномерно. Наличие холодного ядра в факеле искривляет температурное поле, и оно выравнивается только на расстоянии Z/L=0.25. Таким образом, максимальная температура обнаруживается в глубине факела, а не на его оси. Наивысшая температура в факеле наблюдается в начальных сечениях на расстоянии Z/L=0.05, а затем по центру факела: температура холодного ядра повышается и максимальная температура продвигается к осевой линии, так как именно там расположена основная зона реакции.
Полученные результаты расчетов правильно отражают реальную картину топочных процессов.
Литература
1. Госмен, А.Д. Численные методы исследования течений вязкой жидкости / А.Д. Госмен [и др.] -М.: Мир, 1972. - 323 с.
2. Иссерлин, А.С. Основы сжигания газового топлива / А.С. Иссерлин. - Л.: Недра, 1987. - 336 с.
3. Лилли, Д.Г. Простой метод расчета скоростей и давления в сильно завихренных течениях / Д.Г. Лилли // Ракетная техника и космонавтика. - Т.14. - №6. - С.57-67.
4. Садыков, А.В. Некоторые результаты численного решения уравнения энергии методом сплайн-коллокации / А.В. Садыков, А.А. Ахметов // Вестник Казанского технологического университета / КГТУ. - Казань, 2005. - С.98-106.
5. Завьялов, Ю.С. Методы сплайн-функций / Ю.С. Завьялов, Б.И. Квасов, В.Л. Мирошниченко.-М.: Наука, 1980. - 352 с.
6. Марчук, Г.И. Методы вычислительной математики / Г.И.Марчук. - М.: Наука, 1980. - 534с.
7. Горелик, А.М. Современный Фортран для компьютеров традиционной архитектуры и для параллельных вычислительных систем / А.М. Горелик // Вычислительные методы и программирование. - 2004. - Т.5. - С.1-12.
8. Вафин, Д.Б. Влияние особенностей выгорания газообразного топлива на радиционно-конвективный теплообмен в цилиндрических печах / Д.Б. Вафин, А.В. Садыков // Тепло- и массообмен в химической технологии: межвуз. сб. науч. тр. - Казань: Изд-во Казан. хим.-технол ин-та, 1989. - С.31-37.
© А. В. Садыков - канд. техн. наук, доц., декан факультета управления и автоматизации НХТИ КГТУ; И. М. Валеев - асп. каф. теоретических основ теплотехники КГТУ.
