Содержание к диссертации
Введение
Глава 1. Теоретические аспекты моделирования пористых сред 10
1.1. Структура пористых сред 10
1.1.1. Общие представления о пористых материалах 10
1.1.2. Классификация пористых материалов 11
1.1.3. Моделирование структуры пористых материалов 13
1.2. Влияние внешних условий на процесс термического разложения древесины 14
1.2.1. Конечная (максимальная) температура обугливания древесины 14
1.2.2. Скорость обугливания древесины 17
1.2.3. Давление газов и паров при обугливании древесины 18
1.2.4. Влажность древесины 19
1.3. Пиролиз древесины в различных средах 20
1.3.1. Пиролиз древесины в токе парогазов - продуктов разложения древесины 20
1.3.2. Пиролиз древесины в среде углеводородов 21
1.3.3. Разложение древесины в токе паров воды и газов 22
1.4. Модели пористых сред 24
1.4.1. Регулярные упаковки 24
1.4.2. Случайные упаковки 28
Глава 2. Математические модели 31
2.1. Математическая модель конвективного тепломассопереноса в слое пористой среды 31
2.1.1. Приведение уравнений к безразмерному виду 32
2.2. Одноканальная модель процесса термической деструкции пористой частицы 33
2.3. Многоканальная модель процесса термической деструкции пористой частицы 36
2.3.1. Приведение системы уравнений к безразмерному виду 38
2.4. Модель пористой структуры 40
Глава 3. Численные методы 45
3.1. Дискретизация уравнений 46
3.1.1. Дискретный аналог модели конвективного тепломассопереноса в слое пористой среды 46
3.1.2. Дискретный аналог одноканальной модели процесса термической деструкции пористой частицы 50
3.2. Методы решения систем уравнений 52
3.2.1. Линейно - итерационная процедура решения дискретного аналога 52
3.2.2. Уравнение количества движения 54
3.2.3. Поправки скорости и давления 55
3.2.4 Уравнение для поправки давления 56
3.2.5. Граничные условия 57
3.2.6. Алгоритм решения 58
3.3. Программная реализация 64
3.3.1. Блок схема процедуры решения одноканальной модели пиролиза 64
3.3.2. Блок-схема процедуры решения многоканальной модели пиролиза 66
3.3.3. Описание Delphi-приложения 67
Глава 4. Результаты расчётов 70
4.1. Одноканальная модель разложения частицы 70
4.1.1. Результаты моделирования пористой структуры 76
4.1.2. Аналитическое решение нестационарного уравнения теплопроводности 80
4.2. Многоканальная модель разложения частицы 87
4.3. Математическая модель конвективного тепломассопереноса в слое пористой среды 93
Выводы 100
Список литературы
- Конечная (максимальная) температура обугливания древесины
- Одноканальная модель процесса термической деструкции пористой частицы
- Дискретный аналог одноканальной модели процесса термической деструкции пористой частицы
- Аналитическое решение нестационарного уравнения теплопроводности
Введение к работе
Технические, экологические, экономические и иные системы, изучаемые современной наукой, обычно не поддаются исследованию в нужной полноте и точности существующими теоретическими методами. Прямой натурный эксперимент над ними долог, либо попросту невозможен, так как многие из этих систем существуют в «единственном экземпляре». Поэтому математическое моделирование является неизбежной составляющей научно-технического прогресса [1].
В соответствии с современными представлениями процесс моделирования может быть представлен в виде следующей последовательности: исследуемое явление — математические модели — численные алгоритмы — программирование — ЭВМ — вычисления и их анализ — обработка и хранение результатов, дополняющей известную триаду математического моделирования: модель — алгоритм — программа [1].
Будучи методологией, математическое моделирование не подменяет
її собой математику, физику, биологию и другие научные дисциплиы, не
конкурирует с ними. Наоборот, трудно переоценить его синтезирующую роль. Создание и применение триады невозможно без опоры на самые разные методы и подходы — от качественного анализа нелинейных моделей до современных языков программирования. Оно дает новые дополнительные стимулы самым разным направлениям науки [1].
Во всем многообразии задач, решаемых с помощью математического моделирования, особый интерес составляют процессы пиролиза.
Процессы пиролиза (термической деструкции) различных органиче
ских веществ лежат в основе многих технологических процессов (получение
активированного угля и т.д.). Кроме того, процессы пиролиза сопутствуют
горению и газификации природных твердых топлив (уголь и т.д.).
j Пиролиз органических веществ, представляет собой сложный физико-
химический процесс, при котором в результате нагрева из твердого топлива происходит выделение так называемых летучих веществ в виде газообразных
и жидких продуктов: углеводородов, смол, кислот и т.д. Известно, что общее количество и состав вышедших летучих веществ в определенной степени зависят от термических и других условий, в которых осуществляется процесс деструкции — скорости нагрева частиц, максимальной температуры и т.д.
В настоящее время существует большое количество работ, посвященных математическому моделированию процессов термического разложения органических веществ, в первую очередь углей [2-5]. Наиболее простой подход основан на однокомпонентной схеме и реакции первого порядка [2-4]. Двухкомпонентная схема, основанная на предположении, что термическая деструкция и выход летучих веществ происходит по двум параллельным каналам, предложена в работе [5]. Есть и более сложные модели, основанные на многокомпонентном подходе, основанном на представлении органической массы угля в виде совокупности функциональных групп, объединенных в комплексы из ароматических колец. Предполагается, что выход функциональных групп происходит параллельно и независимо друг от друга [6-8].
Характерной особенностью моделей является их применимость к достаточно мелким частицам (50 - 100 мкм) без учета структурных изменений в процессе пиролиза.
Существует широкий класс задач, связанных с процессами пиролиза органического топлива, когда уже нельзя считать частицы мелкими, например задачи, относящиеся к комплексной переработке растительной биомассы, отходов легкой промышленности и др. Прямое использование древесины в различных областях народного хозяйства, лесных продуктов (например, кедрового ореха) в пищевой промышленности и т.д. приводит к появлению большого количества отходов (некондиционной древесины, стружки и т.д.), которые необходимо утилизировать с получением ценных продуктов (активированный уголь и древесный уголь) на основе существующих и новых технологий. Поэтому требуются новые подходы и модели для описания пиролиза крупных частиц, размеры которых составляют несколько сантиметров.
7 Исследование естественной конвекции в насыщенных пористых средах
также представляет интерес для решения научных и технических задач. Этим проблемам посвящено большое число экспериментальных и теоретических работ, наиболее полный обзор которых представлен в [10]. В данной работе теоретическое изучение проводилось в рамках модели Дарси пористой среды. Одной из часто исследуемых конфигураций является горизонтальный цилиндрический слой. Как правило, авторы ограничивались наблюдением симметричного серповидного движения, которое имеет форму двух симметричных вихрей, образованных потоками, поднимающимися вблизи нагретой внутренней поверхности и опускающимися возле холодной внешней границы, и установлением зависимостей теплопереноса от определяющих параметров задачи. В тонких слоях (отношение радиусов цилиндров R>0.2), наряду с серповидным течением, наблюдались симметричные движения, с несколькими дополнительными вихрями в верхней части слоя. С ростом числа Грасгофа над внутренним цилиндром формируется симметричный относительно вертикали конвективный факел, а вблизи границ полости пограничные слои. Систематического исследования конвективных режимов в широком диапазоне значений толщины слоя, числа Прандтля, а также других определяющих параметров задачи, по моим сведениям не проводилось.
В то же время из экспериментов [20, 21] и численных расчетов [22] известно, что в горизонтальном цилиндрическом слое обычной жидкости с большим значением числа Прандтля наблюдается движение в виде конвективного факела, отклоненного от вертикали в ту или иную сторону случайным образом.
Актуальность работы определяется важностью разработки и созданием достаточно простых и надежных алгоритмов и машинных программ для численного моделирования процессов пиролиза органического топлива с целью использования их при решении проблем комплексной переработки растительной биомассы в ценные продукты медицинского и технического назначения, а также при совершенствовании существующих и проектировании
8 новых реакторов для переработки органических топлив.
В связи с этим целью работы является разработка математических моделей, численных методов решения дифференциальных уравнений и создание вычислительных программ, моделирующих процессы тепломассопе-реноса при пиролизе крупных частиц органического топлива с учетом структурных изменений в ходе процесса.
Для достижения этой цели потребовались разработки следующих основных положений, вынесенных на защиту:
математическая модель конвективного течения в горизонтальном цилиндрическом слое пористой среды (модель Дарси);
математическая модель реагирования частицы и ее структуры;
численная модель и программа расчета пиролиза крупных частиц с учетом пористости.
Научная новизна работы состоит в следующем:
разработана математическая модель конвективного тепломассо-переноса в слое пористой среды;
разработана математическая модель (одноканальная и многоканальная) процесса термической деструкции пористой частицы, применимая к крупным частицам;
разработана математическая модель пористой структуры частицы;
получено аналитическое решение нестационарного уравнения теплопроводности в приближении крупной частицы;
созданы программы расчета с использованием разработанных компонентов для отображения процесса термической деструкции крупных частиц органического топлива.
Апробация работы.
Диссертация подготовлена, обсуждена и одобрена на кафедре прикладной математики и информатики Московского государственного индустриального университета филиала в г. Сергиев Посад.
Материалы диссертации были доложены на конференции, посвященной проблемам тепловой конвекции (International Conference ADVANCED PROBLEMS IN THERMAL CONVECTION, г. Пермь, 2003 г.).
По теме диссертации опубликовано 6 печатных работ.
Структура и объем работы.
Структура диссертации определяется логикой поставленных задач и состоит из введения, четырех глав, заключения, библиографического списка использованной литературы (83 источника). Диссертация изложена на ПО страницах печатного текста, содержит 5 таблиц, 44 рисунка, приложение на 25 страницах.
Конечная (максимальная) температура обугливания древесины
Конечная (максимальная) температура обугливания древесины оказывает весьма большое влияние на выход и состав продуктов пиро-генетического разложения древесины.
В таблице 1 показан выход продуктов, полученный при разных конечных температурах обугливания сосновой и березовой древесины [23] при лабораторных исследованиях. Эта таблица показывает, что выход угля из сосновой и березовой древесины с увеличением конечной температуры обугливания все время понижается, причем в пределах температур 270-400С наблюдается быстрое понижение выхода, а в пределах 400-700С -медленное.
Выход жидких продуктов и неконденсирующихся газов в пределах температур 270-400С увеличивается быстро, а с 400С до 700С - медленно. Отсюда следует, что в том случае, когда древесину обугливают, чтобы получить жидкие продукты, можно ограничиться конечной температурой обугливания 380-400С, так как при этой температуре практически заканчивается полностью выделение кислот и спиртов. Дальнейшее повышение температуры вызовет дополнительный расход топлива, что не может быть оправдано.
Изменения в составе жидких продуктов, получаемой при обугливании древесины сосны и березы при разных конечных температурах, показаны в таблице 2 [23]. Из этой таблицы видно, что в начале процесса появляются вода и уксусная кислота, а затем - другие органические вещества. Состав жидких продуктов при изменении конечной температуры обугливания меняется быстро в пределах 270-400С, а при дальнейшем повышении температуры - медленно. Данные о составе жидких продуктов также подтверждают, что конечная температура обугливания 380-400С является вполне достаточной.
Приведенные выше результаты исследований показывают, что при одной и той же конечной температуре обугливания береза дает меньший выход угля и больший выход жидких продуктов, чем хвойные породы; при этом элементарный состав угля, полученного из всех пород, одинаковый.
Скорость обугливания зависит от скорости повышения температуры в камере, от размеров обугливаемой древесины, качества и породы, способа обугливания и т.д. При обугливании древесины в промышленном аппарате непрерывного действия очень трудно определить действительную скорость обугливания отдельного куска древесины, а в периодически действующем реакторе - вообще невозможно.
Поэтому приходится делать выводы о влиянии скорости обугливания на выход продуктов обугливания на основании лабораторных исследований, при которых обугливанию подвергают небольшое количество древесины, и теплота, выделяющаяся при экзотермической реакции, не влияет заметно на повышение температуры обугливаемой древесины.
Влияние скорости нагрева древесины в разных стадиях обугливания неодинаково: при температурах ниже начала экзотермической реакции разложения древесины, т.е. ниже 275С, при нагреве в атмосфере инертных газов можно нагревать древесину долгое время без заметного разложения ее, но при 275-3 00С в пределах экзотермической реакции увеличение быстроты нагрева оказывает очень большое влияние на выходы и состав продуктов обугливания [24].
В таблице 3 [24] приведены результаты исследования, полученные при разной продолжительности обугливания березовой древесины (по Класону).
Изменение продолжительности обугливания древесины в таких широких пределах, как 3 часа и 14 суток, не влияет существенно на выход уксусной кислоты и метилового спирта, но сильно влияет на выход угля и смолы: выход угля увеличивается с 25.51 кг до 39.14 кг, а выход смолы уменьшается с 18 кг до 1.80 кг.
При бурной экзотермической реакции древесины, внутри нее выделяется большое количество газов и паров, которые не могут быстро уйти по порам, имеющимсяся в древесине, и вызывают образование трещин. 1.2.3. Давление газов и паров при обугливании древесины
При обугливании древесины может изменяться давление газов и паров как внутри куска древесины (вследствие большей или меньшей скорости нагрева), так и в реакторе.
Большое давление газов и паров внутри куска древесины приводит к образованию трещиноватого, непрочного угля. Влияние общего давления газов и паров при обугливании древесины в лабораторном аппарате на выход продуктов обугливания показывает таблица 4 [23].
Одноканальная модель процесса термической деструкции пористой частицы
Рассмотрим одиночную частицу органического топлива (например, частицу древесины) радиуса Rp, состоящую из углеродного остатка и летучих веществ (рис. 2.1.). Тогда масса частицы выражается равенством тр=тс+тл=тр(\-У )+трУ\ (2.7) где Шр — масса частицы, тс - масса углеродного остатка, тл - масса летучих веществ, V - начальная доля летучих веществ в частице [39, 91]. Схема выхода летучих веществ
Под действием нагрева летучие вещества переходят из твердой фазы в газообразную и далее диффундируют по порам частицы в окружающую сре ДУ Процесс пиролиза нестационарный и одномерный в сферической системе координат. Пиролиз можно представить в виде глобальной реакции [39, 40]: Органическое топливо -» Углеродный остаток + Летучие вещества
Для моделирования процесса термической деструкции используем одноканальную модель [2]: dV (2.8) dt k(V -V), к = к0 exp Jj (2.9) где V — количество летучих веществ, к - константа скорости, к0 - предэкс-пненциальный множитель, Е - энергия активации, Т — температура частицы, Rg - газовая постоянная. Уравнение (2.8) можно записать через Сл - концентрацию летучих веществ в частице в твердой фазе в виде: дС - = -к0 ехр с„ dt (2.10) Уравнение неразрывности для вышедших летучих веществ Спр в объеме частицы имеет вид в сферической системе координат: ее, "р _ dt mr2 дг Dr дг + S (2.11) где m - пористость, S=ko expf--——) Сл - источниковый член, D- коэффици RgT ент диффузии летучих веществ в порах. Уравнение энергии имеет вид дрсТ , 1 дрсиТг2 _ д_(Лг2 дТ) дг г2 дг dt (2.12) дг где р, Л - плотность и коэффициент теплопроводности частицы, с - теплоемкость летучих веществ.
Движение летучих веществ, происходит по закону Дарси: ЗдР и /л дг где 3 — коэффициент проницаемости, //-динамическая вязкость летучих веществ, Р — давление. C„RT Р = пр s , где Mg - молекулярная масса летучих веществ. Изменение пористости частицы за счет выхода летучих веществ описывается уравнением: m = l-(l-V)(\-m0)- -, (2.13) где т0 — пористость частицы в начальный момент времени. Система дифференциальных уравнений (2.10) - (2.12) дополняется начальными Cn(t = 0)=Cl= - ,Cnp(t = Q) = Cnp, T(t = 0) = T0, (2.14) и граничными условиями = 0) = 0, f (Г = 0) = 0, дС D- -\{r = Rp) = hm{COKp-Cnp) , (2.15) где Токр и Сокр - температура окружающей среды и концентрация летучих веществ в окружающей среде. A. = . 2Rp h = ШЯокр 2Rp 2/3D..1/3 Nu = 2 + Re2/3Pr где D0Kp, Я окр - коэффициенты диффузии и теплопроводности окружающей среды; Рг - число Прандтля; Re - число Рейнольдса; Nu - число Нуссельта.
Количество летучих веществ, перешедших из твердой фазы в газообразную, может быть вычислено по формуле: )C4dV XR=\-i . (2.16) \сч dv о
Тогда выражение для Ху - количество летучих веществ, вышедших из частицы и покинувших ее, имеет вид: \A7iR2hmCnpdt Xv = i (2.17) (4/3)яСлД Процесс термической деструкции может протекать по нескольким (например, двум) параллельным каналам [41]. Скорости термической деструкции угля и выхода летучих веществ в данной модели пропорциональны количеству нетрансформировавшегося древесного угля (У) и описываются дифференциальными уравнениями: dY dt dV_ dt = -{kx +к2)У , (2.18) = {Yxk, +У2к2)У, (2.19) Et E, Ґ С \ Ґ J7 Л kx =A:0exp -- = , k2 =k0exp --zrhz , (2.20) Vy j где кх ,к2 - константы скорости выхода летучих веществ; Ех ,Е2 - энергия активации конкретной древесины; Т— температура частицы; Yx ,Y2-характеристики конкретной древесины.
Важной особенностью настоящей модели является то, что Ех Е2 [42]. Значение Yx равно V , в то время как Y2 представляет дополнительную долю летучих, образовавшихся при более высоких температурах нагрева (часто близко к единице). В литературе для Y приводятся различные данные. Например, рекомендуют значения Yx = 0.3 и Y2 = 1.0.
Изменение массы М j - частицы угля в результате выхода летучих со временем описывается дифференциальным уравнением вида: dMJ (2.21) dt -{Yx Л/ + Г2 V)exp(-W+ )0 Кроме того, данная модель может быть описана с помощью следующей системы уравнений:
Дискретный аналог одноканальной модели процесса термической деструкции пористой частицы
Такое зигзагообразное поле нельзя считать физичным. Можно заметить, что для каждой узловой точки р соответствующий перепад давления Pw РЕ = 0, так как значения давления в соседних с р узловых точках равны между собой. Таким образом, ошибочным следствием данной аппроксимации будет то, что такое волнистое поле давления будет восприниматься в уравнении количества движения как однородное.
Эта трудность ещё более усугубляется в двухмерном случае (для примера рассмотрим прямоугольную систему координат). Так же как на количество движения в направлении оси влияет перепад давления pw-pE, на количество движения в направлении оси Г влияет перепад давления ps-pN, при этом значение давления в точке Р не играет никакой роли. Имея это в виду, R можно сделать вывод о том, что пока- 0 занное на рис. 3.3 поле давления, об- Рис 3. Шахматное поле давления разованное из расположенных в шахматном порядке четырёх произвольных значений давления, не даст силу давления в направлениях осей ф или R . Таким образом, при рассмотренном способе дискретизации уравнений количества движения сильно не однородное поле давления будет восприниматься как однородное. Если бы в процессе итерационного решения возникли бы такие поля давления, они бы не сохранились в процессе, так как уравнения количества движения "забудут" об этих полях. Если в процессе решения будет получено некоторое гладкое поле давления, можно составить любое количество дополнительных решений путём прибавления к этому решению так называемого шахматного поля давления. Так как это поле даёт нулевой градиент давления, уравнение количества движения не почувствует его прибавления.
При заданном поле давления решение уравнений количества движения можно получить с помощью дискретного аналога уравнения для обобщённой переменной Ф. Для уравнения количества движения Ф обозначает одну из составляющих скорости. При использовании шахматной сетки дискретные аналоги уравнений количества движения отличаются от дискретных аналогов уравнения для других Ф, рассчитываемых в узлах основной сетки. Однако это отличие относится к несущественным деталям.
Эти отличия связаны с использованием для аппроксимации уравнений количества движения контрольных объёмов на шахматной сетке.
Контрольный объём для компоненты скорости U показан на рис. 3.4. Он смещён по отношению к обычному контрольному объёму, расположенному вокруг узловой точки Р. Смещение объёма произошло только в направлении оси (р таким образом, что перпендикулярные этому направлению грани проходят через основные узловые точки Р и Е. Отсюда видно одно из главных достоинств шахматной сетки; разность рр - рЕ можно использовать для расчёта силы давления, действующей на контрольный объём для скорости U.
Для расчёта коэффициента диффузии и массового расхода на гранях контрольного объёма, показанного на рис. 3.4, потребуется соответствующая интерполяция. Результирующий дискретный аналог можно записать в виде aeUe = Y.anbUnb +b + (pP -рЕ)Ае. (3.18) Уравнения количества движения можно решить только в том случае, если поле давления задано или каким-то образом найдено. Если при решении использовать не верное поле давления, найденное поле скорости не будет удовлетворять уравнению неразрывности. Выразим такое поле скорости, полученное с использованием приближённого поля давления р , через if , V . Это поле скорости находится в результате решения следующих уравнений: cteU:=ZanbU;b+b + (p;-pl)Ae, (3.19) ay;=5 „л;+ь+(Р; -PN)A„ . (3.20) Найдём способ улучшения приближённого поляр таким образом, чтобы результирующее поле скорости с каждой итерацией лучше удовлетворяло уравнению неразрывности. Предположим, что истинное давление находится из соотношения: Р=Р+Р , (3.21) где р назовём поправкой давления.
Аналитическое решение нестационарного уравнения теплопроводности
Для расчётов использовался ПК следующей конфигурации: процессор AMD Duron(tm) 900 MHz, видеокарта AGP 32 Мб SDRAM GeForce2 MX 100/200 Ver. 4.0, память DDR 256 Мб.
Проведён анализ полученных результатов при варьировании количества точек разбиения и выявлено, что оптимальное разбиение сетки - 50x50.
Значение относительной погрешности, % Как было показано выше, на точность расчёта влияет количество точек разбиения расчётной области, чем больше точек, тем точнее полученные значения, но увеличение узлов расчётной сетки увеличивает процессорное время, затрачиваемое на расчёт. t,c onno На рис. 4.25 изображены структуры течений, рассчитанные в модели
Дарси при R=0.3, Gr=400, Pr=\, Da=l0 6. Из рис. 4.25 следует, что при данных значениях параметров возможны как симметричные относительно плоскости 6= 0 , так и несимметричные движения.
Симметричное течение, представленное на рис. 4.25 а, состоящее из двух серповидных вихрей (далее S1), возникает при сколь угодно малых значениях числа Грасгофа и существует в широком диапазоне значений параметров задачи. С увеличением числа Грасгофа над внутренним цилиндром образуется конвективный факел, а вблизи границ полости формируются пограничные слои. При значениях Gr Gr, (О, 103) это стационарное решение становится неустойчивым. Значение критического числа Грасгофа Gr, зависит от толщины слоя, числа Прандтля и от модели насыщенной пористой среды. Так, для модели Дарси эти значения равны соответственно 3300. Другое симметричное движение (S2), изображенное на рис. 4.25 помимо двух серповидных вихрей содержит дополнительные конвективные ячейки в приполярной области цилиндрического слоя. Причиной образования этих дополнительных вихрей является неустойчивость релеевского типа: в относительно тонких цилиндрических слоях жидкости тепловые условия в верхней части становятся похожими на таковые в плоском слое, подогреваемом снизу. В слоях умеренной толщины, например с R=03, существует только одно симметричное решение такого вида с двумя дополнительными вихрями. В более тонких слоях существует несколько такого рода решений с разным числом вихрей: с двумя, с четырьмя и т.д. Течение с четырьмя дополнительными вихрями реализуется в слое с R=0.S, а с шестью вихрями - для R=0.9. Симметричное движение S2 возникает «жестким» образом при конечном числе Грасгофа и существует в ограниченном диапазоне его значений. В толстых слоях (R 0.\6) многовихревых решений не существуете наблюдается только серповидное течение.
Кроме симметричных стационарных решений в рассматриваемой задаче существует и асимметричное стационарное решение (А1), структура которого показана на рис. 4.25 с. Это движение по характеру близко к решению 67, но отличается от него потерей симметрии относительно плоскости 0=0: интенсивности вихревых движений в левой и правой половине слоя становятся различными, причем более интенсивное движение реализуется в любой половине слоя случайным образом. Данное решение возникает в слоях умеренной толщины для R 0.3 и существует, как и симметричное решение S2, в ограниченном интервале значений числа Грасгофа.
Стационарные решения разной симметрии с ростом числа Грасгофа теряют устойчивость и становятся нестационарными. Особый интерес представляет потеря устойчивости режима с двумя серповидными вихрями (S1). В слоях с R 0.7 (достаточно толстые слои) в окрестности факела возникают тепловые волны, которые распространяются вдоль внешней границы.
