Numerical simulation оf coal combustion in the combustion chamber of a boiler unit

全文:

Введение

С целью снижения угрозы изменения климата уголь должен постепенно исключаться из использования в качестве источника электроэнергии [1]. Однако замена угольных электростанций устойчивыми и возобновляемыми источниками энергии при одновременном удовлетворении постоянно растущего спроса на энергию представляет собой серьезную проблему даже для самых богатых и сознательных стран. Потребность в огромных инвестициях и реструктуризации энергетической инфраструктуры, а также неизбежный длительный переходный процесс, вероятно, обусловят сохранение угля в качестве преобладающего топлива для производства электроэнергии в обозримом будущем во многих странах. Передовые чистые угольные технологии и методы снижения выбросов как варианты решения экологических проблем могут продлить использование угля на некоторое время. Тем не менее подавляющее большинство существующих станций будут стремиться продлить свою текущую эксплуатацию как можно дольше, обеспечивая, по крайней мере, первоначальную проектную эффективность и надежность [2, 3].

Широкое использование численного моделирования для анализа, проектирования и оптимизации процессов во всех отраслях промышленности считается бесценным инструментом для решения текущих и будущих задач в энергетике [4, 5]. Несмотря на некоторые неопределенности, его потенциал в предоставлении исчерпывающей и полной информации о распределении полей всех интересующих параметров пока не может быть оспорен никакими измерениями или другими методами. Таким образом, численное моделирование особенно удобно для диагностики, устранения и предотвращения неисправностей и неоптимальной работы стареющих существующих энергетических установок, а также для оптимизации их работы в непроектных условиях.

Поскольку большинство отечественных тепловых электрических станций введены в эксплуатацию более 40 лет назад, то большая часть энергетического оборудования нуждается в модернизации и ремонте [6]. Одним из наиболее эффективных методов выбора оптимальных вариантов модернизации энергетических установок является применение численного моделирования. Но при работе с пакетами прикладных программ можно использовать различные численные алгоритмы, методы и схемы дискретизации [7, 8]. Для решения задач горения в энергетических установках доступно большое количество методов, в целом их можно разделить на эйлеровы, лагранжевы и смешанные лагранжево-эйлеровы методы. В эйлеровых методах уравнение переноса дискретизируется с помощью метода конечных разностей или метода конечных элементов на фиксированной сетке, тогда как в лагранжевом подходе сетка движется вместе с потоком или остается фиксированной в деформирующей системе координат. В лагранжево-эйлеровом подходе используется двухшаговый метод. Первым шагом является оценка конвективного переноса с использованием траекторий частиц в лагранжевом подходе, а все остальные процессы моделируются в эйлеровом подходе на следующем этапе [9, 10].

Поэтому актуальным является проведение численных экспериментов с целью установления наиболее оптимальных к применению схем аппроксимации при численном решении процессов горения в объеме топочной камеры. Также актуальным является отработка методики проведения численных исследований тепломассообменных процессов с учетом применения различных численных алгоритмов.

Целью работы является численное моделирование процессов окисления полифракционного твердого топлива, формирование аэродинамических потоков и тепломассобмена с учетом применения различных схем аппроксимации.

Физико-техническая постановка задачи

Объектом исследования является барабанный паровой котел (рис. 1) паропроизводительностью 220 т/час с параметрами пара (давление 10.8 МПа, температура 540°C) [11]. Котельный агрегат с естественной циркуляцией, испарительные поверхности, экранирующие стены топочной камеры, выполнены вертикально и свободно удлиняются вниз.

 

Рис. 1. Эскиз топочной камеры: 1 – горелочные устройства.

 

Камера сгорания котельного агрегата имеет форму призмы, размер в плане по осям труб 6 656 × 9 536 мм. Фронтовой и задний экраны в нижней их части образуют скаты “холодной” воронки. Верхняя часть топочного объема ограничена радиационным пароперегревателем. Перед горизонтальным газоходом трубы заднего экрана образуют аэродинамический выступ.

Вихревые горелочные устройства скомпонованы встречно в два яруса треугольником вниз на боковых стенах.

В качестве топлива используется каменный кузнецкий уголь марки Д [12]. При проведении численного моделирования приняты следующие характеристики топлива. Теплота сгорания топлива на рабочую массу составляет 4 949 ккал/кг. Также учтено, что влажность топлива составляет 17.6%, а зольность – 10.2%. Элементный состав кузнецкого угля, следующий: углерод – 56%, водород – 4%, сера – 0.3%, кислород – 10%, азот – 1.9%. Таким образом, при численном расчете учитывается забалластированность топлива, определяющая эффективность его выгорания.

Математическая постановка задачи

В качестве расчетной программы применен пакет прикладных программ FIRE 3D, имеющий высокие качественные показатели при решении подобных задач [13]. Поля скоростей и уровень температуры газовой фазы моделируется с использованием реализуемой k-ε модели турбулентности с высоким числом Re, хорошо учитывающей отрыв, завихрение и другие неравновесные эффекты. Течение и теплообмен пылевидных частиц решались в лагранжевой системе отсчета с применением метода дискретного случайного блуждания для учета влияния турбулентных флуктуаций на их движение, а несущая фаза – на основе эйлерова подхода. Распределение размеров частиц и связанный с ним массовый расход, полученные с помощью ситового анализа, описываются регрессией Розина–Раммлера. Метод источника частиц в ячейке (PSIC) использовался для соединения дискретной и непрерывной фаз с помощью дополнительных условий источника массы, импульса и энергии частиц в соответствующих уравнениях сохранения.

Горение угольных частиц моделируется как пять параллельных и последовательных процессов: нагрев, испарение влаги, выход летучих, воспламенение и горение летучих веществ, воспламенение и выгорание коксового остатка; температура частицы определяется из уравнения энергии с учетом конвективного и радиационного теплообмена [13].

Уравнения изменения массы и температуры угольной частицы [14]:

$\frac{d{m}_i}{dt}=-{\dot{M}}_i^{evap}-{\dot{M}}_i^{vap}-{\dot{M}}_i^{char},$

$m_i{c}_p\frac{d{T}_i}{dt}=Q_{char}{\dot{M}}_i^{char}+\pi {d^2}_{pi}\left[{\alpha }_{conv}\left(T-T_i\right)+{\epsilon }_p\left(H-4\sigma {T^4}_i\right)\right]-$

$- Q_{vap}{\dot{M}}_i^{vap}-L_{evap}\cdot {\dot{M}}_i^{evap}.$

Для расчета массовой скорости выхода влаги из угольной частицы применен следующий закон [15]:

${\dot{M}}_i^{evap}=\pi {d_i}^2{k}_c\left({\left(\rho C_{H_{2}O}\right)}_{surf}-\left(\rho C_{H_{2}O}\right)\right),$

где k_c = (D_H₂O /d_i) (2 + 0.6 Re_i^0,5 Sc^0,33) – коэффициент массопереноса; D_H₂O – коэффициент диффузии паров воды; (ρC_H₂O)_surf – концентрация водяного пара вокруг частицы.

При численном расчете необходимо учитывать массовую скорость выхода летучих Ṁ_i^vap, определяется по закону Аррениуса [16]:

${\dot{M}}_i^{vap}=k_{vap} \text{exp}\left(-\frac{E_{vap}}{R_0{T}_i}\right)\cdot f_{vo}{m}_{oi}.$

Скорость горения углеродного остатка является одной из ключевых величин, которая зависит от нескольких параметров. Одним из них является диффузия окислителя к поверхности угольной частицы. При достижении окислителем поверхности частицы начинается процесс горения. Вторым фактором, влияющим на скорость горения углерода, является скорость кинетической реакции горения [15]:

$M_i^{char}=\pi d_{pi}^2{C}_{O_2}\frac{K_K{K}_D}{K_K+K_D}.$

В данном уравнении $K_k=k_{char}\text{exp}\left(-\frac{E_{char}}{R_0{T}_i}\right)$, kchar, Echar – константы, определяемые по виду топлива согласно [17]. А $K_d=\frac{D_{O2}}{d_{pi}}\left(2+0.16\text{Re}\right)$.

Метод QUICK использует квадратичную интерполяцию вверх по потоку для получения номинальных значений концентраций. Метод можно интерпретировать как линейную интерполяцию, представленную первым членом в правой части уравнения, скорректированную пропорциональным членом по кривизне потока, в правой части. Также стоит отметить, что для расчета граничных значений в схеме QUICK используются концентрации в начале временного шага [18].

Запись данной схемы для грани W в случае неравномерной сетки имеет вид:

${\left(\rho u\Phi A\right)}_W=\left\{\begin{array}{c}{F}_W\left[a^{+}{\Phi }_{WW}+b^{+}{\Phi }_W+c^{+}{\Phi }_P\right], \text{åñëè }{u}_W>\mathrm{0,}\\ F_W\left[a^{-}{\Phi }_{WW}+b^{-}{\Phi }_W+c^{-}{\Phi }_P\right],\text{ åñëè }{u}_W\le \mathrm{0,}\end{array}\right.$

где Φ_WW, Φ_W, Φ_P, Φ_E – результаты значения функции в расчетных узлах; u_W – значение скорости потока на поверхности; a±, b±, c± – коэффициенты интерполяции в многочлене Лагранжа. При этом значение функции на поверхности рассчитывается как интерполяционное значение с учетом трех точек. Аналогично определяются значения и на других гранях.

При рассмотрении монотонизированной противопотоковой схемы (MLU) [19, 20] расчетная траектория s функции между новым значением на грани W и предыдущим определяется с учетом соотношения:

$S_W=\min \text{mod}\left(\frac{\left(a+b\right)}{2.2·\text{min} \text{mod}\left(a,b\right)}\right).$

Функция min mod рассчитывается:

$\min \text{mod}\left(a,b\right)=\left\{\begin{array}{c}a, \text{\hspace{0.33em}если} \left|à\right|\le \left|b\right| и a\cdot b>\mathrm{0,}\\ b, \text{\hspace{0.33em}если} \left|à\right|\ge \left|b\right| и a\cdot b>\mathrm{0,} \\ \mathrm{0,} \text{\hspace{0.33em}если} a\cdot b\le \mathrm{0,}\end{array}\right.$

где a и b демонстрируют кривую функции Φ, с учетом соотношений:

$a=\frac{{\Phi }_P-{\Phi }_W}{\delta {\chi }_W},$

$b=\left\{\begin{array}{c}\frac{{\Phi }_{WW}-{\Phi }_W}{\delta {\chi }_{WW}}, \text{\hspace{0.33em}если} u_W>\mathrm{0,}\\ \frac{{\Phi }_E-{\Phi }_P}{\delta {\chi }_{â}}, \text{\hspace{0.33em}если} u_W\le 0.\end{array}\right.$

В общем случае на немонотонных траекториях функции схема MLU имеет пониженный порядок точности, но остается более предпочтительной к применению в сравнении со схемой QUICK, т.к. она не допускает нефизичных значений функции Φ.

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

Результаты и обсуждение

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

На рис. 2 и 3 представлены данные о скоростях потоков и векторах скоростей в продольных и поперечных сечениях топки. Эти результаты показывают, что внутри топки наблюдается подъемное движение с достижением 8…12 м/с. Это подтверждает, что происходит эффективное перемешивание воздушных потоков внутри горелочных устройств. Также на рис. 2 и 3 видно распространение потоков, которые формируются при подаче пылевоздушной смеси и вторичного воздуха на выходе из горелочных устройств. Горелочные струи формируют форму факела в соответствие с направлением и параметром крутки в вихревых горелочных устройствах. При этом стоит отметить, что при организации расчета по схеме MLU формирование закрученных потоков происходит более активно, нежели чем при расчете с применением схемы аппроксимации QUICK (рис. 2, 3). В частности, при схеме QUICK наблюдается существенная раскрутка факела на выходе из горелочных устройств, и это приводит к малому проникновению горелочных струй по ширине топочной камеры и к значительным зонам обратных токов вдоль оси горелок (рис. 3). Из-за этого на рис. 2б не наблюдается сформированный закрученный аэродинамический поток от нижней горелки в центральном сечении по продольной оси топки.

 

Рис. 2. Векторное поле скоростей в продольном сечении по оси топочной камеры (м/с): (а) MLU; (б) QUICK.

 

Рис. 3. Векторное поле скоростей в поперечном сечении по оси топочной камеры (м/с): (а) MLU; (б) QUICK.

 

Данный факт оказывает влияние и на изменение температуры в топочной камере (рис. 4, 5). Результаты численного моделирования с использованием схемы MLU (рис. 4а) имеют высокотемпературные зоны, которые представляются более равномерными, а область активного горения формируется во всем объеме топочной камеры до пода аэродинамического выступа. Во втором случае по схеме QUICK в центральной части сечения наблюдается зона пониженных температур, которая соответствует положению факела от горелочного устройства нижнего яруса.

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

 

Рис. 4. Распределение температуры в вертикальном сечении по продольной оси топки (К): (a) MLU; (б) QUICK.

 

Рис. 5. Распределение температуры в горизонтальном сечении по оси верхнего яруса горелок (К): (a) MLU; (б) QUICK.

 

На рис. 6 и 7 представлено распределение концентрации кислорода в вертикальных сечениях по оси топки и в горизонтальных сечениях по оси горелочных устройств верхнего яруса. Аналогично, как и с результатами формирования температурных полей в топочной камере, вследствие малого взаимопроникающего учета потоков воздуха и топлива при формировании факела в горелках нижнего яруса при расчете по схеме QUICK, наблюдается повышенная концентрация О₂ (рис. 6б).

 

Рис. 6. Концентрация О₂ в вертикальном сечении по продольной оси топки (%): (a) MLU; (б) QUICK.

 

Рис. 7. Концентрация О₂ в горизонтальном сечении по оси верхнего яруса горелочных устройств (%): (a) MLU; (б) QUICK.

 

Несмотря на качественное различие результатов моделирования в близи фронтовой и тыльной стен топочной камеры наблюдаются высокотемпературные области (рис. 5) с пониженной концентрацией кислорода (рис. 7). Соответственно данная среда совместно с высоким содержанием CO (рис. 8) создает благоприятные условия для возникновения высокотемпературной коррозии. Таким образом, при значительном содержании CO и отсутствии O₂, а также в случае увеличения содержания серы в поставляемом угле, будут образовывать сульфиды различных элементов, а не оксиды сульфидов, которые преобладают в угольных котлах, особенно H₂S, который агрессивно действует на стальные трубы экранных поверхностей. Считается, если значение CO превышает 0.01, то это является признаком неполного горения и при непосредственном расположении данных зон у стен топочной камеры угрожает возникновением высокотемпературной коррозии. Подобные локальные зоны наблюдаются у фронтовой и тыльной стен (рис. 8).

 

Рис. 8. Концентрация CО₂ в горизонтальном сечении по оси верхнего яруса горелочных устройств: (a) MLU; (б) QUICK.

 

Температурный уровень по высоте топочной камеры представлен на рис. 9. Также представлены значения натурных измерений температуры через лючки топки. В области холодной воронки при расчете с учетом схемы MLU температурный уровень примерно на 200°C ниже, чем при расчете с применением схемы QUICK. Соответственно во втором варианте наблюдаются нехарактерно высокие температуры в данной зоне. В дальнейшем при достижении высоты расположения горелочных устройств температурный уровень снижается из-за ввода через горелочные устройства топливно-воздушной смеси и вторичного воздуха с относительно низкими значениями температур.

 

Рис. 9. Распределение температуры по высоте топки, °C.

 

При этом при аппроксимации по схеме QUICK снижение температурного уровня на нижнем ярусе горелок выражено более крутым падением. Далее температуры возрастают и достигают максимума 1500°C, после зоны активного горения восходящий поток продуктов сгорания охлаждается, отдавая тепло посредством конвективного и лучистого теплообмена ограждающим поверхностям топочной камеры. При этом, согласно рис. 8, значения температур, полученных при численном моделировании и натурных измерениях на уровне ядра горения, отличаются незначительно (не более 100°C). В дальнейшем расхождение значений натурных измерений и результатов моделирования представляется крайне низкими значениями.

В работе предварительно также оценивалось влияние дискретизации структурированной сетки в трехмерной постановке на результаты численного моделирования. К расчету принимались три варианта с количеством ячеек 104 760, 271 760 и 516 344. При рассмотрении изменения температуры по высоте топочной камеры было выявлено, что разность максимальных температур между первым и вторым вариантами в зоне активного горения составила 170°C, а между вторым и третьим 50°C. Таким образом к последующим расчетам принят вариант с 271 760 ячеек.

Выводы

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

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

При этом в обоих вариантах моделирования вблизи фронтовой и тыльной стен топочной камеры наблюдаются локальные зоны с пониженным содержанием O₂ и повышенным CO, что можно рассматривать как индикаторы высокотемпературной коррозии.

Работа выполнена при финансовой поддержке Российского научного фонда (проект № 23-29-00274, https://rscf.ru/project/23-29-00274/).

作者简介

A. Gil

National Research Tomsk Polytechnic University

编辑信件的主要联系方式.
Email: andgil@tpu.ru
俄罗斯联邦, Tomsk

K. Maltsev

National Research Tomsk Polytechnic University

Email: andgil@tpu.ru
俄罗斯联邦, Tomsk

N. Abramov

National Research Tomsk Polytechnic University

Email: andgil@tpu.ru
俄罗斯联邦, Tomsk

S. Puzyrev

National Research Tomsk Polytechnic University

Email: andgil@tpu.ru
俄罗斯联邦, Tomsk

参考

补充文件