Процессы и аппараты химической технологии

Размещено на http://www.allbest.ru/

Введение

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

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

FEMLAB имеет большую встроенную библиотеку готовых моделей, которая является важной частью этого программного пакета. Библиотека содержит готовые модели, описывающие явления из различных областей инженерии. На каждую модель предоставляется подробнейшая документация, включающая технические выкладки, обсуждение результатов и пошаговое описание всех установок модели. Поскольку модели уже идут со определенными сеткой и решением, вы можете «поиграть» с различными опциями результатов моделирования.

Начиная с версии 3.0 FEMLAB полностью интегрирован с MATLAB. Фактически, необходимость в интеграции была и раньше, но появилась только в 3-ем выпуске FEMLABа. Теперь же, стало возможным запускать FEMLAB вместе с MATLAB, что позволило программировать в FEMLAB используя язык программирования MATLAB.

1. Моделирование теплообменного процесса

1.1 Теоретическая часть

Перенос тепла определяется как движение энергии вследствие разности температур. Она осуществляется по следующим механизмам:

Проводимость - перенос тепла посредством диффузии в стационарной среде по градиенту температуры. Среда может быть твердой или жидкой.

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

Радиация (излучение) - перенос тепла между поверхностью A с температурой T₁ и поверхностью B с температурой T₂ посредством электромагнитных волн, причем предполагается, что T₁ T₂ и поверхность A является бесконечно малой для наблюдателя на поверхности B.

Тепловое уравнение

Математическая модель посредством проводимости описывается тепловым уравнением:

, где

T - температура

r - плотность

C - теплоемкость

C_P - теплоемкость при постоянном давлении

C_V - теплоемкость при постоянном объеме

k - удельная теплопроводность

Q - количество теплоты (переданной или поглощенной).

Оператор «набла» есть взятие градиента:

В установившемся режиме температура постоянна и первый член теплового уравнения стремятся к нулю:

Если среда анизотропна, то k задается в виде матрицы:

Для модели тепловой проводимости и конвекции в жидкости или газе тепловое уравнение содержит так же и конвективную составляющую. В FEMLAB такая модель описывается уравнением:

, где

U - поле скоростей среды.

Вектор теплового потока (одновременно для проводимости и конвекции) определяется из вышестоящего выражения в скобках:

, где

q - вектор теплового потока.

Для одной проводимости, без учета конвекции, вектор теплового потока принимает вид:

Доступными краевыми условиями в данном режиме являются:

1) Тепловой поток

2) Изоляционный материал

3) Симметрия

4) Диапазон температур

Определение уравнения в частных производных

Математическая модель переноса тепла посредством проводимости описывается следующей разновидностью теплового уравнения:

, где

d_ts - время-масштабный коэффициент.

В некоторых случаях, исследователь может решить добавить локально-проявляющуюся поперечную конвекцию или излучение в одно- или двумерном пространстве. Это осуществляется добавлением двух дополнительных членов справа:

, где

h_trans - коэффициент теплообмена посредством поперечной конвекции

T_ext - поперечная внешняя температура

C_trans - константа, определяемая пользователем

T_ambtrans - поперечная температура окружающей среды

dA - толщина в одномерном пространстве и площадь в двумерном.

Примечание: Когда используется поперечная конвекция или излучение, тепловое уравнение необходимо модифицировать для использования области dA. Для всех практических целей, когда использование этого расширения исключается, члены h_trans и C_trans устанавливаются равными нулю.

Обычно этот коэффициент равен 1, но вы можете изменить временной масштаб (установив, например, значение данного коэффициента равным 1/60, что будет означать, что за секунду реального времени FEMLAB промоделирует 60 секунд).

Плотность

Введите плотность материала. Вводимая величина должна представлять отношение единицы массы к единице объема.

Теплоемкость

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

Удельная теплопроводность

Удельная теплопроводность характеризует отношение между вектором теплового потока q и температурным градиентом T: - что является законом теплопроводности Фурье. Вводимая величина должна представлять единицу мощности на единицу длины и на единицу температуры.

Источник тепла

Источник тепла характеризует выработку тепла внутри области. Нагревание и охлаждение выражается соответственно положительным и отрицательным значениями. Значение Q должно вводится как единица мощности на единицу длины (в одномерном пространстве), единица мощности на единицу площади (в двумерном пространстве), единица мощности на единицу объема (в трехмерном пространстве).

Виды граничных условий

Граничные условия	Описание
	Тепловой поток
	Обособленность или симметрия
	Заданная температура
	Заданный ноль температуры

Тепловой поток

Диалоговое окно граничных условий включает начальные условия теплового потока, конвекции и излучения под видом граничных условий для теплового потока. Тепловой поток q₀ предполагается направленным нормально внутрь. Составляющие конвекции и излучения предполагаются направленными нормально наружу.

Значение q₀ определяется как тепловой поток, входящий в данную область (домен). Например, любой электрический нагреватель хорошо описывается этим условием и при этом можно пренебречь геометрией самого нагревателя. Величина q₀ вводится как единица мощности в одномерном пространстве, как единица мощности на единицу длины в двумерном пространстве, как единица мощности на единицу площади в трехмерном пространстве.

Слагаемое h(T_inf - T) моделирует конвективный обмен теплом с окружающей средой, где h - коэффициент теплообмена и T_inf - окружающая температура. Значение h зависит от геометрии и общих условий потока. Исчерпывающее представление о том, как считаются коэффициенты теплообмена, будет дано ниже.

Теплообмен с окружающей средой посредством излучения моделируется . T_amb - температура радиации окружающей среды, которая может отличаться от T_inf. C_const - произведение поверхностного излучения ? и постоянной Стефана-Больцмана ? = 5,67 * 10^-8 . Поверхностное излучение - свойство материала, которое описано ниже.

C_const =?es

Примечание: Задача с границами радиации нелинейна, поэтому для стационарных решений используется нелинейный алгоритм решения.

Обособленность или симметрия

Это условие используется, чтобы задать, где домен каким-либо образом изолирован или для уменьшения размера модели, если она симметрична. Глядя на выражение, можно предположить, что градиент поперек границы должен быть равен нулю. Для этого температура по одну сторону границы должна быть равна температуре по другую сторону границы. Поскольку если нет разности температур на границе, тепло через нее переноситься не будет.

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

Осевая симметрия

Это граничное условие доступно только для моделей теплообмена с осевой симметрией. Используйте его для описания оси симметрии, для которой r = 0.

Задание температуры

T = T₀

Это условие используется для задания температуры T₀. Температура поверхности, для которой задано такое граничное условие, считается постоянной.

Задание нуля температуры

T = 0

Это ограничение используется для ввода нуля раздела температур.

Граничные предустановки

Граничные условия вводятся в диалоговом окне граничных предустановок.

Переменная	Описание
q0	Внутренний тепловой поток
h	Коэффициент теплообмена посредством конвекции
Tinf	Температура окружающей среды
Cconst	Константа излучения: произведение излучательной способности и константы Стефана-Больцмана.
Tamb	Температура излучения окружающей среды
T0	Заданная температура

Режим конвекции и проводимости

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

Переменные и размерность пространства

Режим конвекции и проводимости может работать в одно-, двух-, и трехмерном пространстве, также осесимметричном одно- и двумерном пространстве.

Определение уравнения в частных производных

В этом прикладном режиме вы можете выбрать одну из двух предложенных формулировок:

(нераспространенная)

(распространенная)

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

Вы можете переключаться между распространенной и нераспространенной формами в диалоговом окне application mode properties посредством из выбора в закладке equation form.

Так же можно добавить искусственную стабилизацию (можно выбрать «isotropic diffusion», «streamline diffusion» или «crosswind diffusion») для подгонки численного решения вручную.

Настройки подсистемы

Таблица ниже содержит величины, содержащиеся в уравнениях:

Обозначение	Описание
?ts	Время-масштабный коэффициент
?	Плотность среды
CP	Теплоемкость среды
k, kij	Удельная теплопроводность среды
Q	Тепловой поток
u, v, w	Компоненты вектора скорости движения среды
?id	Параметр настройки для «isotropic diffusion»
?sd	Параметр настройки для «streamline diffusion»
?cd	Параметр настройки для «crosswind diffusion»

Граничные условия

Единственное отличие от режима теплопереноса заключается в том, что тепловой поток содержит конвективную составляющую. В случаях, когда конвекция через границу превышает диффузию, используют граничные условия потока конвекции. Они устанавливают, что поток диффузии на границе равен нулю, однако поток конвекции может проникать через границу. Это задается граничным условием «конвективный поток».

1.2 Практическая часть

Для моделирования выберем теплообменник типа «труба в трубе». Смоделируем его в системе FEMLAB и получим кривые переходных процессов. В качестве реального прототипа выберем процесс нагрева сырья колонны К-101 блока ВТ установки У-731. В качестве теплоносителя выберем перегретый водяной пар.

При моделировании приняты следующие допущения:

1) теплопроводность стенок труб велика по сравнению с теплопроводностью среды внутри теплообменника

2) пар достаточно перегрет для того, чтобы вклад процесса конденсации греющего пара в теплообмен был мал по сравнению с процессом охлаждения

3) конвективные потоки, вызванные изменением плотности жидкости при нагреве, слабы по сравнению с ее движением вдоль теплообменника

Создадим модель в системе FEMLAB, для чего выберем из предложенного списка моделей heat transfer, convection and conduction, transient analysis. Укажем, что пространство двумерное (2D).

Примем следующие геометрические размеры теплообменника: длина 12 м, радиус внутренней трубы 0,5 м, внешней трубы 1,5 м. После построения получаем следующее:

Рис. 1. Модель теплообменника типа «труба в трубе»

Арабскими цифрами помечены границы, римскими - области с одинаковыми параметрами среды, или домены.

Переходим к заданию параметров областей (physics, subdomain settings). Во внешней рубашке (области I и III) идет перегретый пар, начальная температура которого 250^оС. Смоделируем процесс, когда в разогретый теплообменник подают на нагрев холодное сырье с температурой 20^оС. Физические параметры сырья (стабильного гидрогенизата фракция НК - 350^оС) предположим близкими к параметрам смеси предельных и ароматических углеводородов С₆ - С₁₀. Учитывая, что время задаем в минутах для удобства, вводим следующие параметры:

Таблица 1

Параметр	Размерность	Область I	Область II	Область III
Время-масштабный коэффициент ?ts	безразмерный	1	1	1
Теплопроводность k (изотропная)		1,75	7	1,75
Плотность ?		8	800	8
Теплоемкость Cp		2000	1880	2000
Внутреннее тепловыделение Q		0	0	0
Скорость в горизонтальном направлении u		0,27	0,42	0,27
Скорость в вертикальном направлении v		0	0	0
Начальная температура T(0)	К	523	523	523

Скорости потоков рассчитываются исходя из расхода пара 50 м³/час, а сырья 10 м³/час. Далее зададим граничные условия (physics, boundary settings). Учтем потери тепла в окружающую среду через стенки внешней трубы (границы 4 и 5). Принимая температуру окружающей среды 20^оС, потери с достаточной точностью могут быть описаны формулой:



Таблица 2

Граница	Вид условия	T0, K	q0,
1, 3	Температура	523	-
2	Температура	293	-
4, 5	Тепловой поток	-
6, 7, 8	Конвективный теплообмен	-	-
9, 10	Расчетная граница	-	-

Зададим параметры решения задачи (solve, solver parameters). Выберем решение, зависимое от времени (time dependent) и время до 50 мин с шагом 0,1 мин. Относительную и абсолютную точность зададим 0,0001. Все остальные параметры оставим без изменения.

Инициализируем сетку (mesh, initialize mesh). Предложенная сетка весьма грубая, поэтому сделаем ее более мелкой (mesh, refine mesh). Чем мельче элементы сетки, тем точнее решение, но тем больше ресурсов оно требует. Я остановился на сетке из 3456 областей. Статистику разбиения можно увидеть в меню mesh, mesh statistics. Внешний вид сетки показан на рисунке.

Рис. 2. Разбиение модели на расчетные области

Запускаем решение (solve, solve problem). Через некоторое время мы можем увидеть результат решения - распределение температур в теплообменнике.

Рис. 3. Распределение температур внутри теплообменника

Как видно, FEMLAB предоставляет возможность наглядной демонстрации результатов расчета. Однако, необходимо представить результаты моделирования в виде кривой разгона. Для этого воспользуемся сопряжением со средой MATLAB 6.5. Запустив FEMLAB в связке c MATLAB, экспортируем структуру с параметрами модели в MATLAB (file, export, FEM structure).

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

filename: get_h.m

function T = get_h (fem, t, step_t, x, y, r, step_r)

T = zeros (length (t(1):step_t:t(2)), 1);

for i = 0:step_r:r

if (i == 0) | (i == r)

k = 0.5;

else

k = 1;

end

T = T + postinterp (fem, 'T', [x y+i]', 'T', t(1):step_t:t(2))*i*k*step_r;

end

T = 2 * T / r^2;

plot (t(1):step_t:t(2), T);

Запустим ее командой

T = get_h (fem, [0 50], 0.1, 6, 0, 0.25, 0.01);

Получим массив T, содержащий в себе переходной процесс в системе, заданный с шагом 0,1 мин. Построим переходной процесс.

Рис. 4. Переходной процесс в системе при подаче сырья

Как видно, на графике есть всплески температуры, не могущие иметь места в реальности и обусловленные ошибками моделирования. Для их сглаживания воспользуемся фильтром скользящего среднего.

filename: mean_filter.m

function Signal_filtered = mean_filter (Signal_noised, m)

Signal_filtered = [];

for i = 1: (length (Signal_noised) - m)

s = [];

for j = 0:m - 1

s = [s Signal_noised (i+j)];

end

Signal_filtered = [Signal_filtered mean(s)];

end

for i = (length (Signal_noised) - m + 1):length (Signal_noised)

s = [];

for j = 0: (length (Signal_noised) - i)

s = [s Signal_noised (i+j)];

end

Signal_filtered = [Signal_filtered mean(s)];

end

После фильтрации получаем следующий переходной процесс.

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

Смоделируем переходной процесс, происходящий в системе при подаче большего количества греющего пара. Для этого изменим время решения на 250 мин. с шагом 0,1 мин., а также скорость движения в доменах I и III (расход пара). Зададим повышение расхода пара с 50 до 250 м³/час на 150-й минуте моделирования выражением:

u = 0.27 + 1.08 * flsmhs (t - 150, 1)

После завершения моделирования и экспорта в MATLAB, выполним команду

T = get_h (fem, [150 250], 0.1, 6, 0, 0.25, 0.01);

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

Построим исходную и аппроксимированную кривые разгона. На рисунке исходная кривая показана красным.

2. Моделирование химической реакции

2.1 Теоретическая часть

Рассмотрим моделирование химического реактора, в котором происходит процесс частичного окисления о-ксилола до фталевого ангидрида. Процесс происходит в многотрубном реакторе при температуре порядка 400^оС на катализаторе из оксида ванадия и сульфата калия, нанесенных на основу из диоксида кремния. Время протекания реакции колеблется от 0,5 до 5 с. Наиболее важным параметром в процессе является температура в реакторе, поскольку реакции, проходящие в реакторе, сильно экзотермичны, и при нарушении охлаждения реактора он теряет устойчивость. Распределение температуры внутри реактора влияет на его КПД и выход фталевого ангидрида, который, естественно, необходимо повысить. Температуру можно контролировать путем варьирования диаметра трубок реактора, скорости потока в трубках, температуры стенок трубки (степени охлаждения).

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

Химизм процесса, протекающего в реакторе, достаточно сложен, и не может быть здесь детально рассмотрен, но упрощенно его можно представить как совокупность 3 реакций:

1) окисление о-ксилола до фталевого ангидрида.

2) горение о-ксилола с образованием смеси оксидов углерода.

3) горение фталевого ангидрида с образованием смеси оксидов углерода.

Как видно, целевая реакция - первая. Поскольку процесс ведется при значительном избытке кислорода, можно считать, что реакции первого порядка. Обозначив константы скорости реакций k1, k2 и k3 соотвественно, запишем уравнения:

, где

y₀ - молярная доля кислорода

yA0 - молярная доля о-ксилола

xA - общая скорость расхода о-ксилола

xB - скорость превращения о-ксилола во фталевый ангидрид

xС - скорость образования оксидов углерода

Коэффициенты скорости реакций в соответствии с законом Аррениуса определяются из выражения:

, где

T₀ - температура стенки реактора

- перегрев среды относительно стенки

Рис. 8. Модель трубки реактора

Определение уравнения в частных производных

Уравнения, описывающие данный процесс, выглядят следующим образом.

, где

u_S - поверхностная скорость

rg - плотность газа

rb - плотность катализатора

сtot - общая концентрация

leff - эффективная теплопроводность среды

Учтя, что xA = xB + xC, решаем уравнения только для xB и xC. После подстановок и упрощений получаем:

Уравнения скорости реакций:

Граничные условия



Полагаем, что на выходе реактора определяющее влияние имеют конвекция и теплопередача.

Задание модели

Ввиду того, что длина трубки реактора значительно превосходит ее диаметр, масштабируем уравнения для большей наглядности моделирования. В данном случае, длина превышает диаметр более чем в 200 раз - трубки имеют в длину 3 метра при диаметре 1,27 см. Масштабирование также избавит нас от проблем при формировании сетки.

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

2.2 Практическая часть

Создадим модель химических процессов в реакторе в системе FEMLAB, для чего выберем из предложенного списка моделей chemical engineering module, mass balance, convection and diffusion. Укажем, что пространство двумерное осесимметричное (axial symmetry 2D). Введем через пробел зависимые переменные (dependent variables): xb и xc. Добавим модель, нажав add. Теперь создадим модель тепловых процессов в реакторе, для чего выберем из предложенного списка chemical engineering module, energy balance, convection and conduction, и выставим распространенную форму уравнений теплового баланса (application mode properties, equation form, conservative).

Создаем прямоугольник, как показано на рисунке.

Рис. 9. Модель трубки реактора

Зададим параметры среды в полученной области. Для удобства и наглядности определим несколько констант (options, constants), которыми будем пользоваться в дальнейшем. Константы, имена переменных и значения изложены в таблице 3.

Таблица 3

Константа	Имя переменной	Значение
Коэффициент диффузии Deff	Deff	3,19 * 10-7
Поверхностная скорость us	us	1,064 * 10-3
Плотность газа rg	rhog	1293
Плотность катализатора rb	rhob	1300
Теплопроводность среды l	lambda	7,8 * 10-4
Теплоемкость среды CP	cp	0,992
Общая концентрация ctot	ctot	44,85
Коэффициент теплоотдачи a	alpha	0,156
Температура стенки T0	T0	627
Тепловой эффект реакции DH1	deltaH1	-1285000
Тепловой эффект реакции DH3	deltaH3	-4564000
Молярная доля о-ксилола yA0	ya0	0,00924
Молярная доля кислорода y0	y0	0,208
Коэффициент уравнения Аррениуса B1	B1	13588
Коэффициент уравнения Аррениуса B2	B2	15803
Коэффициент уравнения Аррениуса B3	B3	14394
Коэффициент уравнения Аррениуса A1	A1	114500
Коэффициент уравнения Аррениуса A2	A2	318480
Коэффициент уравнения Аррениуса A3	A3	48113
Масштабирование по радиусу scaler	scaler	0,0127
Масштабирование по высоте scalez	scalez	0,6

Определим несколько выражений (options, expressions, scalar expressions).

Таблица 4

Переменная	Имя переменной	Значение
Константа скорости реакции k1	k1	A1*exp (-B1/(T+T0))
Константа скорости реакции k2	k2	A2*exp (-B2/(T+T0))
Константа скорости реакции k3	k3	A3*exp (-B3/(T+T0))
Скорость реакции rB	rb	ya0*y0*(k1*(1-xb-xc) - k2*xb)
Скорость реакции rC	rc	ya0*y0*(k3*(1-xb-xc)+k2*xb)

Зададим параметры среды внутри трубки. Для энергетического баланса зададим стабилизацию вида «streamline diffusion».

Таблица 5

Параметр	Значение
Уравнения массового баланса (convection and diffusion)
Время-масштабный коэффициент dts	1
Коэффициент диффузии D (анизотропный)	[Deff/scaler^2 0 0 Deff/scalez^2]
Скорость реакции R	rb*rhob/ctot/ya0 rc*rhob/ctot/ya0
Скорость движения среды u по оси r	0
Скорость движения среды v по оси z	us/scalez
Начальная концентрация xB(0)	0
Начальная концентрация xC(0)	0
Уравнения энергетического баланса (convection and conduction)
Время-масштабный коэффициент dts	1
Теплопроводность k (анизотропная)	[lambda/scaler^2 0 0 lambda/scalez^2]
Плотность r	rhog
Теплоемкость Cp	cp
Внутреннее тепловыделение Q	rhob*(-deltaH1*rb-deltaH3*rc)
Скорость движения среды u по оси r	0
Скорость движения среды v по оси z	us/scalez
Начальная температура T(0)	0
Параметр стабилизации dsd	0.25

Таблица 6

Граница	Вид условия	Значение параметра
Уравнения массового баланса (convection and diffusion)
1, 4	Изоляция или симметрия	-
2	Концентрация	0
3	Конвективный поток	-
Уравнения энергетического баланса (convection and conduction)
1	Теплоизоляция	-
2	Температура	0
3	Конвективный поток	-
4	Тепловой поток	-alpha*T/scaler

Зададим параметры решения задачи (solve, solver parameters). Выберем решение, зависимое от времени (time dependent) и время до 3000 сек с шагом 10 сек. Относительную и абсолютную точность зададим 0,01. Все остальные параметры оставим без изменения.

Инициализируем сетку (mesh, initialize mesh). Предложенная сетка весьма грубая, поэтому сделаем ее более мелкой (mesh, refine mesh). Я остановился на сетке из 744 областей. Статистику разбиения можно увидеть в меню mesh, mesh statistics. Внешний вид сетки показан на рисунке.

Запускаем решение (solve, solve problem). Через некоторое время мы можем увидеть результат решения - распределение температур в трубке реактора. Шкала показывает превышение температуры относительно стенки реактора (627 К).

Как видно, в начале трубки имеется зона сильно повышенной температуры. Это так называемая «горячая точка», явление, достаточно распространенное в реакторах, в которых проводятся экзотермические реакции с охлаждением.

Модифицируем программу для снятия кривой разгона соответственно изменениям модели (иные геометрические координаты).

filename: get_h.m

function T = get_h (fem, t, step_t, z, r, step_r)

T = zeros (length (t(1):step_t:t(2)), 1);

for i = 0:step_r:r

if (i == 0) | (i == r)

k = 0.5;

else

k = 1;

end

T = T + postinterp (fem, 'T', [i z]', 'T', t(1):step_t:t(2))*i*k*step_r;

end

T = 2 * T / r^2;

plot (t(1):step_t:t(2), T);

Запустим ее командой

T = get_h (fem, [0 3000], 10, 5, 1, 0.05);

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

Для повышения выработки фталевого ангидрида, целевого продукта реакции, температуру необходимо повышать. Однако, это следует делать крайне осторожно. Смоделируем повышение температуры стенки на 5 К и построим переходной процесс, принимая за выходной сигнал среднюю температуру по сечению трубки на уровне «горячей точки». Для этого зададим:

T0 = 627 + 5 * flsmhs (t - 3000, 1);

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

T = get_h (fem, [3000 3350], 10, 1, 1, 0.5);

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

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

теплообменный конвекция химический реактор

1. А.М. Плановский, В.М. Рамм, С.З. Каган. Процессы и аппараты химической технологии. М.: «Химия», 1968 г.

2. Варгафтик Н.Б. Справочник по теплофизическим свойствам газов и жидкостей. М.Ж. «Физматгиз», 1963 г.

3. G.F. Froment. Fixed bed catalytic reactors. Ind. Eng. Chem. №59 (2), 1967 г.

4. Справочное руководство по системе конечноэлементностных расчетов Femlab.

Размещено на Allbest.ru