Численное моделирование физико-химических процессов в слоевом газификаторе Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 10, № 5, 2005

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ФИЗИКО-ХИМИЧЕСКИХ ПРОЦЕССОВ В СЛОЕВОМ ГАЗИФИКАТОРЕ

И. А. КузовАтов Красноярский государственный технический университет, Россия

А. А. Гроо, С. Г. Степанов Карбоника-Ф, Красноярск, Россия e-mail: grooaa@yandex.ru

The system of equations of gas flow balance in the fixed-bed gasifier is solved by a decomposition of physical-chemical processes. The decomposed system describes heat-and-mass transfer with geterogeneous reactions. An algorithm of numerical solution based on the exponential fitting has been developed.

Введение

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

Математическое моделирование процессов тепло- и массообмена при газификации угля осложнено гетерогенными и гомогенными химическими реакциями. Скорость гомогенных реакций при характерных условиях (P = 101.3 кПа, T = 1000 K) намного превосходит скорость гетерогенных реакций на поверхности угольных частиц. Разномасштабность процессов затрудняет построение эффективного численного алгоритма.

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

1. Постановка задачи и формулировка проблемы

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

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2005.

Рис. 1. Схема слоевого газификатора.

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

Задача математического моделирования процессов тепло- и массопереноса в слоевом газификаторе заключается в расчете концентрационных и температурных профилей и определении их динамики при заданных расходах и составах газифицирующего агента и твердого топлива.

Примем следующие допущения: реактор одномерный, адиабатический. Газовый поток предполагается идеально перемешанным, а химический состав — термодинамически равновесным. Так как скорость фильтрации газа достаточно велика (И,е > 50), диффузионным переносом вдоль основного направления потока пренебрегаем [3]. Взаимодействие углерода с газовыми реагентами описывается следующими брутто-реакциями:

с + 02 ^ СО2 + дс2, с + СО2 ^ 2со - дСо2, с + Н2О ^ со + Н2 - дн2о-

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

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

Определение теплопроводности засыпки угольных частиц, претерпевающих физико-химические превращения — очень сложная задача. Существуют методики расчета этого

параметра с учетом кондуктивной передачи тепла в газе, угольных частицах, теплообмен излучением и др. (см., например [4]). В настоящей работе для простоты теплопроводность твердой фазы принята постоянной и определяется из условий соответствия интегральных расчетных параметров процесса (скорость "тепловой волны", максимальная температура на фронте) экспериментальным данным.

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

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

Уравнение энергии твердой фазы с краевыми условиями отсутствия теплового потока через границу:

с дТс _ д2Тс + + Рс_ Лс дх2 + УЬе + ЯсЪет, (1)

3T(

C

дх

О, дТс

x=0 дх

= О,

x=L

где Тс и рс — температура и кажущаяся плотность твердой фазы, К и кг/м3; — тепловой поток от газа к угольным частицам, Вт/м3; дсьет — источник (сток) тепловой энергии при химическом реагировании, Вт/м3.

Начальные условия имитируют зажигание в верхней части слоя, они имеют вид

4 4

Тс (х, 0) _ Тт[п, К при 0 < х < -Ь; Тс (х, 0) _ Ттах, К при -Ь < х < Ь,

5 5

где Ь — высота газификатора.

Скорость движения "тепловой волны" в диапазоне режимных параметров, представляющих интерес для практического использования, достигает 0.2 м/ч и на три порядка ниже скорости газового потока [1, 5, 6]. Так как перемещение температурных и концентрационных профилей практически изоморфно, можно допустить, что газовая фаза в каждый момент времени т0 достигает установившегося состояния для данного мгновенного распределения Тс (х,т0).

Тепло- и массообмен с химическими реакциями формально описывается следующей системой обыкновенных дифференциальных уравнений первого порядка:

± (ир) _ г°2 + + гН20 + тГАР + тГЬ + г?; (2)

d

— (UpCg Tg) = qhe + qg, (3)

где pi и p — концентрация i-й (i = O2, CO2, CO, H2O, H2, N2, VOL) компоненты и плотность газа, кг/м3; U — скорость газового потока, м/с; Cg и Tg — теплоемкость и температура газа, Дж/(кг-К) и K; r°2, rC°2 и rH2° — скорости изменения концентрации i-го компонента при химическом реагировании углерода с O2, CO2 и H2O соответственно, кг/(м3-е); rVAP = О (i = H2O) и rVOL — скорость сушки и выхода летучих веществ, кг/(м3-е); rh — скорости изменения концентрации i-й компоненты в гомогенных реакциях (включая окисление и разложение летучих веществ), кг/(м3-е); qg — источник, описывающий перенос тепла с продуктами газификации, Вт/(м3).

Расчет химического состава, температуры и плотности газа (термохимического поля) с использованием уравнений переноса газовых компонентов (2), (3) достаточно трудоемок и нецелесообразен, так как гомогенные реакции протекают быстро и образуют равновесную систему, которая определяется исходным набором реагентов и энтальпией веществ.

Для построения эффективного алгоритма расчета термохимического поля воспользуемся известным методом расщепления [7], который заключается в следующем: на первом этапе из правых частей уравнений (2) отбросим сильно нелинейные источники, описывающие гомогенные превращения, а на втором — выход летучих веществ и сушку. Полученная

система включает только гетерогенное реагирование. Запишем ее в виде

&

ах (до2 ) = - ь02до2; (4)

(дсо2) = - ьсо2 дсо2 + ^^ ь02 до2; (5)

ах «о2

а

ах (зно) = - ьн2одн2о; (6)

(дсо) = «со ьсо2дсо2 + «со Ьн2одн2о; (7) ах «со2 «н2 о

(зн) = «н2 ьн2одн2о; (8) ах «н2о

а

ах

12 (Зя) = /я; (9)

а

ах (зт) = -Ьт зт +/т' (10)

где д. (г = О2, СО2, Н2О, СО, Н2) = р%и; Зд = Ряи; дт = РяиСдТ; р., рд — концентрации ре-

8¿к;

и

агирующих компонентов и плотность газа, кг/м3; Тя — температура газа, К; Ьг =

Ьт = — (аБ + /я) и /т = Тс (аБ + /я) — коэффициенты и источник, м-1, м-1 и Вт/(м3-е); дя

к.[ — коэффициенты реакционного газообмена для частиц коксового остатка, м/с; Б а и 8 — удельная поверхность "сжимающегося углеродного ядра" и контурная поверхность угольных частиц, м2/м3; а — коэффициент теплоотдачи, Вт/(м2-К);

( _ «с , ас «с ,

/я =-Ьо2 Зо2 +--Ьсо2 Зсо2 +--Ьн2оЗн2о

«о2 «со2 «н2о

— вспомогательная функция.

Система уравнений (4)-(10) дополняется начальными условиями

Ро2 (0) = Р0о2, Р° (0) = Р°°, Ря (0) = £Рг (0), и (0) = и0, Тя (0) = Тт1п. (11)

.

Для определения плотности и скорости газового потока используются уравнение состояния идеального газа и допущение о постоянстве давления в реакторе [2, 8].

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

газовой фазы производится методом минимизации функции Гиббса [7]. Окончательный результат определяется итерационно. Затем производится решение системы (4)-(10) в следующем узле и т. д.

Полная система уравнений, описывающая нестационарные процессы в слоевом газификаторе, содержит параболическое уравнение (1), систему линейных дифференциальных уравнений первого порядка (4)-(10) и замыкающие уравнения, которые здесь не приводятся.

Аппроксимация нестационарного уравнения теплопроводности (1) производится инте-гроинтерполяционным методом с использованием чисто неявной схемы по времени. Решение разностного аналога находится методом прогонки [9].

Наличие быстро исчезающих (до2) и возникающих (дсо2, дн2) компонентов обусловливает жесткость системы обыкновенных дифференциальных уравнений (4)-(10), описывающей тепло- и массообмен с химическими реакциями в газовом потоке. В математическом отношении проявление свойств жесткости в задачах с начальными данными связано с присутствием в уравнении малого параметра при старшей производной. В данном случае использование классических методов численного интегрирования затруднительно.

2. Алгоритм решения системы уравнений тепло- и массообмена с гетерогенным реагированием

Вычислительный алгоритм основан на применении одношаговой схемы расчета с использованием экспоненциальной подгонки.

На отрезке интегрирования [0, Ь] зададим равномерную сетку

Ь

и =| хг | XI = Нхц г = 0,Жж; Нх = — | •

Рассмотрим линейное дифференциальное уравнение

еп'(х) + а(х)и(х) = f (х), и (0) = и0. (12)

Схема расчета вытекает из общей теории линейных дифференциальных уравнений и опирается на формулу точного решения задачи (12):

х

А(х) 1 / А(£)

и (х) е= и (хо) +— f (С) е~~ ^, (13)

хо

где А (х) = / а (С) ^ - первообразная для функции а (х).

хо

Различные способы дискретизации точного соотношения (13) для двух соседних узлов сетки хг и х+\ приводят к различным вариантам разностной схемы. Эти способы отли-

Х А(£)

чаются аппроксимацией для функции А (х) и последующим вычислением J f (С) е ~~ ^

хо

(см. [10- 14]. Переменные коэффициенты системы дифференциальных уравнений (4)-(10) будем считать кусочно-постоянными, т.е. положим а (х) = а (х), х £ [х^х^]. Значение

а (х) вычисляется как среднее на интервале, а интеграл — по формуле трапеций:

Жг+1

- 1 [ /х- - а (х.+1) + а (х.)

а = т а (х) ах, а ^--—1. (14)

к X 2

хг

После упрощения схему можно представить в виде

гЦ'+1ек - "' = ¿/ / (х. + {) е¥¿С. (15)

0

В данной работе, если не удается вычислить интеграл в формуле (15) точно, используется линейная интерполяция функции / (х) для приближенного вычисления. При указанном способе построения схема (15) обладает вторым порядком аппроксимации и первым порядком сходимости по малому параметру е. Достаточно полное изложение доказательства абсолютной устойчивости схемы и обоснование равномерной сходимости по малому параметру и порядка аппроксимации приведены в обзорной монографии Дулана с соавторами [10].

Система дифференциальных уравнений (4)-(10) может быть записана в следующем виде:

а

-г (У1) = ЬиУ1; (16)

dx

dx

dx

dx

dx

d

-jz Ы = b2iyi + 622^2; (17) d

-JZ (Уз) = 6ззУз; (18) d

Ы = Ь42У2 + Ь43У3; (19) d

-г (У5) = Ь5зУз; (20)

d

— (Уб) = 6б1 У1 + Ьб2У2 + ЬбзУз; (21)

d

— (У7) = 671У1 + 672У2 + ЬтзУз + 67. (22)

dx

Нетрудно заметить, что система уравнений (16)-(22) обладает нижнетреугольной матрицей, что позволяет построить алгоритм численного решения, основанный на последовательном интегрировании отдельных уравнений.

Рассмотрим алгоритм построения разностных схем для уравнений (4)-(10) в порядке их решения более подробно. Для уравнения (4) в соответствии с общей формулой (15) и с учетом нулевой правой части получим

gO+ = gO 2 (23)

где 6o2 ( x) — , г , а

6о2 определяется по формуле (14).

U (x)

Величину до2 (х) = дО2 е Ьо2(х можно найти не только в узлах сетки, но и внутри интервала. Это позволяет построить схему для уравнения (5) по формуле (15), используя

подстановку 'с°2 Ьо2д°2е Ьо2^ вместо f (С), и вычислить интеграл точно. В итоге получим '°2

схему

дС° = (■■д° 2 Ь°2 е('С0д - "с" }Ь' -1 + дС 0) е-ьсо2\ (24)

\ '°2 Ьс°2 - Ь°2 )

Бкс°2 (х) — , ,

где Ьс°2(х) = —и ( )— , а Ьс°2 определяется по формуле (14).

Схему для уравнения (6) получим аналогично (4):

дН+° = дн 2°е-Ьн2^ (25)

Уравнение (7) с учетом постоянства молярных масс может быть проинтегрировано аналитически. Действительно, из уравнений (4)-(6) следует

дС°,Ьс°2 = ^ д°2- I (дс°2) ■ д°2Ь°2 = ^ ^ • = ^ <**>)•

Сделав соответствующие замены, запишем уравнение (7) в виде

d ( 2-co . 2-co . -co

~T~ gCO +--go 2 +--gCO2 +--

ax \ -02 -CO2 -H20

Тогда

. 2-co . 2-co . -co _ г1

gCO +--go 2 +--gC02 +--gH20 — cco,

-O2 -CO2 -H2O

где Cco — const определяется из начальных условий

C — d + 2-co gi + 2-co g , -co

eC0 — gC0 + --g02 + --gC02 + --gH20-

-O2 -CO2 -H2O

В результате с учетом уже найденных значений g^1, gC+)2 и gi+0 получаем расчетную формулу

gi+i — C 2-co gi+i 2-co gi+i -co gi+i (96)

gco — cco - —— g02 - --gC02 - "-gH20- (96)

-O2 -CO2 -H2O

Для уравнения (8) после аналогичных преобразований имеем

gj+1 — CH2 - gH+O, (97)

-H2O

ГУ i I — H2 i

где — gH2 + - gH20 •

-H20

Уравнение (9) также допускает аналитическое интегрирование:

gg+1 — C - ^ gOV - gCO - gH20, (98)

-O2 -CO2 -H2 0

i , 2-C i . -C i . -C i где C9 — gg + -— go2 + --gc02 + --gH20-

-02 -C02 -H20

d

1 d

В уравнении (10) сделаем подстановку fg = — (gg) из формулы (9), замену — ~¡~9g =

d ln gg

dx

dx

и примем следующие обозначения для коэффициентов:

gg dx

d ln (gg (x))

¡т (1)=Гс ^. (ж) 5+(ж) ^, 6т (ж) ^—^ Аппроксимируя правую часть кусочно-линейной функцией, получим схему расчета:

a (x) S d ln (gg (x))

gT+1 = gT e-bT h +

л fTi+1 fTi t

M + —h—*

3bT * d£.

(29)

Таким образом, формулы (23)-(29) позволяют найти решение жесткой системы обыкновенных дифференциальных уравнений (4)-(10) в г + 1-м узле и перейти к дальнейшим вычислениям согласно алгоритму расщепления.

h

3. Результаты численных расчетов

Расчет производился при следующих условиях: Ттах = 1000 К, Ттп = 300 К. Температура, скорость газа и мольные доли газовых компонентов на входе: Тд (0) = 300 К, ид = 0.065 м/с, объемное содержание 02 и N2 на входе соответственно о!рЬло2 = 0.2, о!рка^2 = 0.8.

На рис. 2 представлено мгновенное распределение концентраций газовых компонентов по высоте реактора. По оси ординат отложены объемные доли газовых компонентов, по оси абсцисс — высота реактора. Проанализируем изменения концентраций реагирующих компонентов в зонах, указанных на рис. 1, 2. Первой по потоку начинает возрастать концентрация паров воды вследствие нагрева и сушки слоя. Затем следует зона пиролиза. Так как кислорода для полного сгорания всего количества выделившихся летучих веществ достаточно, в данной зоне имеются продукты полного окисления — Н20, С02 и 02. Далее свободный кислород вступает во взаимодействие с углеродом и полностью расходуется. Следующая зона — восстановления, кислород отсутствует, растет концентрация Н2 и СО вследствие гетерогенного реагирования.

X, % об

щ

10

Í 0.05 0 1 015 0 2 l h.u

Рис. 2. Мгновенное распределение концентраций газовых реагентов по высоте газификатора.

т,кл

Рис. 3. Динамика температурных профилей газа и слоя по высоте газификатора.

На рис. 3 представлена динамика профилей температуры газового потока и слоя. Результаты расчетов показывают, что перемещение концентрационных и температурных профилей происходит изоморфно. Скорость "тепловой волны" 15 см/ч « 4.17 10-5 м/с и состав газа качественно соответствуют данным натурных измерений в промышленном слоевом газификаторе [5, 6].

Заключение

Выполнено моделирование процессов тепло- и массопереноса в слоевом угольном газификаторе. Предложен вычислительный алгоритм для решения системы обыкновенных дифференциальных уравнений первого порядка, описывающих тепло- и массоперенос с гетерогенным реагированием. Данная система получена расщеплением по физико-химическим процессам исходной системы, учитывающей процессы сушки, выхода летучих веществ и химическое реагирование на поверхности угольной частицы и в объеме газа. Алгоритм основан на методе экспоненциальной подгонки, который дает решение с погрешностью О (И2) и обладает равномерной сходимостью по малому параметру. Используемые сильные упрощения (модель "сжимающегося углеродного ядра", постоянный коэффициент теплопроводности твердой фазы и др.) все же позволяют адекватно отразить физико-химические процессы, происходящие в слоевом газификаторе. Для получения достоверных количественных зависимостей необходимо провести идентификацию параметров модели по экспериментальным данным.

Список литературы

[1] СТЕПАНОВ С.Г. Промышленные технологии переработки угля: перспективы использования в Канско-Ачинском угольном бассейне. Красноярск: Краснояр. гос. ун-т, 2002. 85 с.

[2] Гроо А.А., Кузоватов И.А., Степанов С.Г. Численное решение системы обыкновенных дифференциальных уравнений, описывающих тепло- и массообмен с химическими реакциями в слоевом газификаторе // Математические методы и моделирование: Вест. КГТУ. 2003. Т. 33. С. 72-79.

[3] Канторович Б.В. Введение в теорию горения и газификации твердого топлива. М.: Ме-таллургиздат, 1961.

[4] Аэров М.Э., Тодес О.М., НАринский Д.А. Аппараты со стационарным зернистым слоем: гидравлические и тепловые основы работы. Л.: Химия, 1979.

[5] Патент 2014882 РФ. МКИ B 01 J 26/20, C 01 B 31/08. Способ получения адсорбента / С.Р. Исламов, С.Г. Степанов, А.Б. Морозов, В.С. Славин (РФ). — 92004035/26; Заявлено 11.11.92; Опубл. 30.06.94, Бюл. 12. Приоритет 11.11.92.

[6] Патент 2014883 РФ. МКИ B 01 J 26/20. Способ получения углеродного адсорбента / С.Р. Исламов, С.Г. Степанов, А.Б. Морозов (РФ). — 93039409/26; Заявлено 16.08.93; Опубл. 30.06.94, Бюл. 12. Приоритет 16.08.93.

[7] Pratt P.T., Wormeck J.J. CREK. A Computer Program for Calculation of Combustion Reaction Equilibrium and Kinetics in Laminar or Turbulent Flow. Report WSU-ME-TEL-76-1. Washington State Univ., 1976.

[8] Степанов С.Г., Исламов С.Р. Математическая модель газификации угля в слоевом реакторе // Химия твердого топлива. 1991. № 2. С. 52-58.

[9] Самарский А.А. Теория разностных схем. М.: Наука, 1983.

[10] Дулан Э., Миллер Дж., Шилдерс У. Равномерные численные методы решения задач с пограничным слоем. М.: Мир, 1983.

[11] Титов В.А., Шишкин Г.И. Численное решение задачи Коши для обыкновенного дифференциального уравнения с малым параметром при производной // Численные методы механики сплошной среды: Сб. науч. тр. / АН СССР. Сиб. отд-ние. ВЦ; ИТПМ. 1978. Т. 9, № 7. С. 112-121.

[12] Боглаев И.П. О численном интегрировании сингулярно-возмущенной задачи Коши для обыкновенного дифференциального уравнения // Журн. вычисл. математики и мат. физики. 1985. Т. 25, № 7. С. 1009-1022.

[13] Doolan E.P., Schilders W.H.A. Uniformly convergent methods for initial value problems // Boundary and Interior Layers — Computational and Asymptotic Methods. Dublin: Boole Press, 1980. P. 256-259.

[14] Зверев В.Г. О численном решении сингулярно-возмущенной задачи Коши для дифференциального уравнения первого порядка // Вычисл. гидродинамика. Томск: Том. гос. ун-т. 1999. С. 73-82.

Поступила в редакцию 6 декабря 2004 г., в переработанном виде — 21 февраля 2005 г.