Logo GenDocs.ru

Поиск по сайту:  


Загрузка...

Курсовая работа - Математическое моделирование в системе конечноэлементностных расчетов Femlab - файл 1.doc


Курсовая работа - Математическое моделирование в системе конечноэлементностных расчетов Femlab
скачать (320 kb.)

Доступные файлы (1):

1.doc320kb.19.11.2011 12:02скачать

содержание
Загрузка...

1.doc

Реклама MarketGid:
Загрузка...
Содержание

Введение

1

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




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

2

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

10

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




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

17

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

20

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

27

Введение

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

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

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

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

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

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

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

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

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

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

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

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

, где

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

 - плотность

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

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

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

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

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

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



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



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



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

, где

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

, где

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

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



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

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

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

  3. Симметрия

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


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

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

, где

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

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

, где

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

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

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

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

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

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

Величины подсистемы:

Обозначение


Описание

ts

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



Плотность среды

Сp

Теплоемкость среды

k

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



Матрица удельной теплопроводности среды

kij

ij-й компонент матрицы k

Q

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

htrans

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

Text

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

Ctrans

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

Tambtrans

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

dA

Толщина или площадь

^ Время-масштабный коэффициент

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

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

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

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

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

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

Описание



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



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



Заданная температура



Заданный ноль температуры

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



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

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

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

Теплообмен с окружающей средой посредством излучения моделируется . Tamb - температура радиации окружающей среды, которая может отличаться от Tinf. Cconst - произведение поверхностного излучения  и постоянной Стефана-Больцмана  = 5,67 * 10-8 . Поверхностное излучение – свойство материала, которое описано ниже.

Cconst =

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



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

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

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

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

T = T0

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

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 оС) предположим близкими к параметрам смеси предельных и ароматических углеводородов С6 – С10. Учитывая, что время задаем в минутах для удобства, вводим следующие параметры:

Таблица 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 м3/час, а сырья 10 м3/час.

Далее зададим граничные условия (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
После фильтрации получаем следующий переходной процесс.



Рис. 5. Переходной процесс в системе при подаче сырья (фильтрованный)

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

Смоделируем переходной процесс, происходящий в системе при подаче большего количества греющего пара. Для этого изменим время решения на 250 мин. с шагом 0,1 мин., а также скорость движения в доменах I и III (расход пара). Зададим повышение расхода пара с 50 до 250 м3/час на 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);

Получаем следующую кривую.



Рис. 6. Переходной процесс при подаче дополнительного греющего пара

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



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



Рис. 7. Исходный и аппроксимированный переходной процесс при подаче дополнительного греющего пара

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

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

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

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

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

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

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

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

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

, где

y0 – молярная доля кислорода

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

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

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

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

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

, где

T0 – температура стенки реактора

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



Рис. 8. Модель трубки реактора
Определение уравнения в частных производных

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

, где

uS – поверхностная скорость

g – плотность газа

b – плотность катализатора

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

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

Учтя, что 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

Плотность газа g

rhog

1293

Плотность катализатора b

rhob

1300

Теплопроводность среды

lambda

7,8 * 10-4

Теплоемкость среды CP

cp

0,992

Общая концентрация ctot

ctot

44,85

Коэффициент теплоотдачи

alpha

0,156

Температура стенки T0

T0

627

Тепловой эффект реакции H1

deltaH1

-1285000

Тепловой эффект реакции H3

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)

Время-масштабный коэффициент ts

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)

Время-масштабный коэффициент ts

1

Теплопроводность k (анизотропная)

[lambda/scaler^2 0 0 lambda/scalez^2]

Плотность 

rhog

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

cp

Внутреннее тепловыделение Q

rhob*(-deltaH1*rb-deltaH3*rc)

Скорость движения среды u по оси r

0

Скорость движения среды v по оси z

us/scalez

Начальная температура T(0)

0

Параметр стабилизации sd

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. Внешний вид сетки показан на рисунке.



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

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



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

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

Модифицируем программу для снятия кривой разгона соответственно изменениям модели (иные геометрические координаты).
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, чтобы получить абсолютные значения температуры.



Рис. 12. Переходной процесс после начала реакции

Как видно, график имеет характерную форму с максимумом. Это хорошо известный и имеющий место быть и на практике в трубчатых реакторах с применением катализа эффект, изученный еще в 1967 г. Г. Ф. Фроментом.

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

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

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

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

График переходного процесса показан на рисунке.



Рис. 13. График переходного процесса при повышении температуры стенки на 5 К

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

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

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

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

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

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



Скачать файл (320 kb.)

Поиск по сайту:  

© gendocs.ru
При копировании укажите ссылку.
обратиться к администрации
Рейтинг@Mail.ru