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

Подождите немного. Документ загружается.

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.

Константа Имя переменной Значение

Коэффициент диффузии D

eff

Deff

3,19 * 10

-7

Поверхностная скорость u

1,064 * 10

-3

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

rhog

1293

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

rhob

1300

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

lambda

7,8 * 10

-4

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

0,992

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

tot

ctot

44,85

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

alpha

0,156

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

627

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

deltaH1

-1285000

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

deltaH3

-4564000

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

ya0

0,00924

Молярная доля кислорода y

0,208

Коэффициент уравнения Аррениуса B

13588

Коэффициент уравнения Аррениуса B

15803

Коэффициент уравнения Аррениуса B

14394

Коэффициент уравнения Аррениуса A

114500

Коэффициент уравнения Аррениуса A

318480

Коэффициент уравнения Аррениуса A

48113

Масштабирование по радиусу scale

scaler

0,0127

Масштабирование по высоте scale

scalez

0,6

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

Таблица 4.

Переменная Имя переменной Значение

Константа скорости реакции k

A1*exp(-B1/(T+T0))

Константа скорости реакции k

A2*exp(-B2/(T+T0))

Константа скорости реакции k

A3*exp(-B3/(T+T0))

Скорость реакции r

ya0*y0*(k1*(1-xb-xc)-k2*xb)

Скорость реакции r

ya0*y0*(k3*(1-xb-xc)+k2*xb)

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

стабилизацию вида «streamline diffusion».

Таблица 5.

Параметр Значение

Уравнения массового баланса (convection and diffusion)

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

Коэффициент диффузии D (анизотропный)

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

Скорость реакции R

rb*rhob/ctot/ya0

rc*rhob/ctot/ya0

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

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

us/scalez

Начальная концентрация x

(0)

Начальная концентрация x

(0)

Уравнения энергетического баланса (convection and conduction)

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

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

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

Плотность 

rhog

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

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

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

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

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

us/scalez

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

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

0.25

Зададим граничные условия.

Таблица 6.

Граница Вид условия Значение параметра

Уравнения массового баланса (convection and diffusion)

1, 4 Изоляция или симметрия

2 Концентрация

3 Конвективный поток

Уравнения энергетического баланса (convection and conduction)

1 Теплоизоляция

2 Температура

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 сек. Построим переходной процесс, предварительно прибавив к каждому

элементу массива значение T

, чтобы получить абсолютные значения температуры.

Рис. 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.