Заключение

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

Рассматриваемый химический реактор представляет собой динамическую систему с пятью входными величинами F, Сн, Tн, Fc, Tc, каждая из которых приводит к изменению двух его выходных величин С и T.

Математическое описание нестационарного режима работы химического реактора должна представлять собой систему из двух уравнений (по числу выходных величин объекта): материального и теплового балансов.

Рисунок 1

Схемы химического реактора идеального перемешивания (а) и его динамических каналов (б) продуктов реакции Q --расход F реагентов, подаваемых в реактор. При этом изменение расхода хладоагента Fc вызывает также изменение состава продуктов реакции Q, а колебание расхода исходных реагентов F приводит к изменению температуры Т реакционной массы в реакторе. Кроме этого, выходные величины реактора (Q и T) зависят от концентрации Qн и температуры Тн входного продукта, а также от температуры хладоагента Тс. Помимо прямых рассматриваемый объект имеет и перекрестные каналы прохождения сигналов. Выходные величины такого реактора находят из уравнений динамики

Q=f₁(F,F_c,Q_н,T_н,T_c,t)

T=f₂(F_c,F,Q_н,T_н,T_c,t)

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

Объекты с сосредоточенными параметрами. К ним относятся объекты, регулируемые величины которых имеют одно числовое значение в данный момент времени (уровень жидкости в аппарате, давление газа в газгольдере и др.). (Типичными представителями объектов с сосредоточенными параметрами являются резервуар для жидкостей, испаритель, химический реактор.

Составить математическое описание химического реактора непрерывного действия, представляющего собой аппарат с мешалкой и охлаждающей рубашкой, в котором проводится экзотермическая реакция 1-го порядка. Объемный расход и температура подаваемой в него реакционной смеси составляют F, м³/с и Тн, К, а концентрация реагирующего вещества равна Qн. Объект обладает сосредоточенными параметрами, т. е. во всем объеме реакционного пространства температура Т, К и концентрация реагирующего компонента Q имеют одинаковые значения. Поток из реактора выходит с теми же значениями величин. Тепло реакции снимается хладоагентом, температура которого равна Тс, К, а расход через охлаждающую рубашку реактора поверхностью А, м² составляет Fc, м³/с Плотности с, кг/м³ исходного и прореагировавшего веществ одинаковы. Удельные теплоемкости с, Дж/(кг-град) продуктов на входе и выходе из реактора равны между собой. Объем реакционной массы V, м³, постоянен.

При принятых допущениях уравнение материального баланса реактора по реагирующему компоненту в переходном режиме:

FQ_н + VS = FQ + VdQ/dt (1.1)

где S = K(Q_max--Q)--скорость реакции; Q_max -- максимально достижимая концентрация целевого компонента; K=a*ехр(--E/RT)--константа скорости реакции; а -- постоянный коэффициент, Е -- энергия активации; R -- универсальная газовая постоянная.

С учетом последних выражений уравнение материального баланса реактора:

FQ_н + Va(Q_max-Q)exp(-E/RT) = FQ + dQ/dt (1.2)

Уравнение теплового баланса реактора в неравновесном режиме может быть записано следующим образом:

F_cсT_н + q_p = F_cсT + V_cсdT/dt + q_c (1.3)

Первые слагаемые левой и правой частей уравнения (1.3) определяют величины тепловых потоков, связанных с поступлением в реактор исходного и отводом из него прореагировавшего веществ; слагаемое q_p определяет тепловой поток, обусловленный теплотой реакции; второе слагаемое в правой части отражает теплоаккумулирующие свойства реактора; слагаемое q_c определяет тепловой поток, отводимый системой охлаждения. При этом

q_p = rPс(Q - Q_н)(1.4)

q_c = K_стA(T - T_c)(1.5)

где r --теплота реакций; К_ст -- коэффициент теплопередачи через охлаждаемую стенку.

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

Nu = 0.33Re^0.5Pr^0.33(1.6)

(где Nu, Re и Pr --критерии Нуссельта, Рейнольдса и Прандтля), можно считать, что коэффициент теплопередачи Кст пропорционален корню квадратному из расхода хладоагента Fc, т. е.

K_ст = bFc^0.5(1.7)

где b -- коэффициент пропорциональности.

Из уравнений (1.5) и (1.7) найдем

q_c = bAFc^0.5(T - Tc)(1.8)

Подставляя зависимости для q_p и q_c из (1.4) и (1.8) в (1.3), получим уравнение теплового баланса реактора

F_cсT_н + rFс(Q - Q_н) = F_cсT + V_cсdT/dt + bAF_c^0.5(T - T_c)(1.9)

Составим уравнения материального и теплового балансов реактора для равновесных режимов его работы

F₀Q_н0 + Va(Q_max - Q₀)exp(-E/RT) = F₀Q₀(1.10)

F_0cсT_н0 + rF₀с(Q₀ - Q_н0) = F_0cсT₀ + bAF_c0^0.5(T₀ - T_c0)(1.11)

где F_o, T_о, Q_o, Q_Ho, F_c0,T_Н0, T_со -- базисные значения соответствующих величин реактора.

После линеаризации нелинейных слагаемых в (1.2) и (1.9), составления системы уравнений реактора в приращениях и введения относительных входных и выходных величин

y₁ = ДQ/Q₀

y₂ = ДT/T₀

x₁ = ДF/F₀

x₂ = ДF_c/F_c0

x₃ = ДQ_н/Q_н0

x₄ = ДT_н/T_н0

x₅ = ДT_c/T_c0

получим:

VQ₀dy₁/dt + Q₀(F₀ + Vaexp(-E/RT₀)y₁ = -F₀(Q₀-Q_н0)x₁ + F₀Q_н0x₃ + VT₀(Q_max - Q₀)aexp(E/RT₀)(E/RT₀²)y₂(1.12)

V_cсT₀dy₂/dt + T₀(F₀cс + bAF_c0^0.5)y₂ = F₀с(r(Q₀ - Q_н0) - с(T₀ - T_н0))x₁ - 1/2 bAF_c0^0.5(T_с - T_c0)x₂ - rсF₀Q_н0x₃ + FcсT_н0x₄ + bAF_c0^0.5T_с0x₅ + rсF₀Q₀y₁ (1.13)

Если ввести обозначения

k₀ = aexp(-E/RT₀), k₀^* = E/RT₀² , k^*=bA F_c0^0.5

поведение рассматриваемого реактора в неравновесном режиме описывается системой уравнений:

T₁dy₁/dt + y₁ = -k₁x₁ + k₂x₃ + k₃y₂ (1.14)

T₂dy₂/dt + y₂ =k₄x₁ - k₅x₂ - k₆x₃ + k₇x₄ + k₈x₅ + k₉y₁(1.15)

где : T₁ = V/(F₀ + Vk₀)

T₂ = V/(F₀cс + k^*)

k₁ = F₀(Q₀ - Q_н0)/Q₀(F₀ + Vk₀)

k₂ = F₀Q_н0/Q₀(F₀ + Vk₀)

k₃ = VT₀(Q_max - Q₀)k₀k₀^*/Q₀(F₀+Vk₀)

k₄ = k^*(T₀ - Tc₀)/(F₀cс + k^*)T₀

k₅ = k^*(T₀ - Tc₀)/2(F₀cс + k^*)T₀

k₆ = rсF₀Q_н0/(F₀cс + k^*)T₀

k₇ = F₀cсT_н0/(F₀cс + k^*)T₀

k₈ = k^*T_c0/(F₀cс + k^*)T₀

k₉ = rсF₀Q₀/(F₀cс + k^*)T₀

Далее рассмотрим определение некоторых передаточных функций реактора. Для нахождения, например, передаточной функции W11(p) по каналу «расход исходного продукта F -- концентрация Q вещества в потоке на выходе», т. е. по каналу x₁--y₁, запишем уравнение динамики в операторной форме, исключив предварительно все остальные входные величины.

(T₁p + 1)y₁ = -k₁x₁ + k₃y₂(1.16)

(T₂p + 1)y₂ = k₄x₁ + k₉y₁ (1.17)

Выражая из полученного у₂, подставляя его в уравнение (1.16) и находя отношение y₁/x₁, найдем

W₁₁(p) = (k₃k₄ - k₁(T₂p + 1))/((T₁p+1)(T₂p+1) - k₃k₉)(1.18)

Покажем также нахождение передаточных функций по структурной схеме реактора, приведенной на рисунок 1.1 и составленной в соответствии с его математическим описанием. Для нахождения, например, передаточной функции W₂₂(p) по каналу «расход хладоагента Fc -- температура потока на выходе из реактора T», т. е. по каналу x₂--y₂, составим структурную схему реактора по каналу x₂--y₂ и ее упрощенные варианты (рисунок 1.1). Передаточные функции динамических звеньев упрощенных схем определяются равенствами:

W^*(p) = 1/(T₂p+1)

W^**(p) = 1/(T₁p + 1)

W_oc(p)=k₃k₉/(T₁p + 1)

Искомая же передаточная функция имеет вид:

W₂₂(p) = - k5(T₂p + 1)/((T₁p + 1)(T₂p + 1) - k₂k₉)(1.19)

Аналогичным образом находят передаточные функции и по другим каналам. Укажем, что знаменатели всех передаточных функций имеют одинаковое выражение, которое обозначим через D(p) = (T₁p+l)(T₂p+l)--k3k9. Выражения передаточных функций химического реактора с учетом этого обозначения сведены в таблица 1, где также приведены обозначения соответствующих входных и выходных величин.

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

Таблица 1. Передаточные функции химического реактора

Газофазная реакция N₂ O₄=2NO₂

2. Получение передаточной функции объекта по заданным динамическим каналам

Перенесем все рассчитанные формулы в Matlab и получим передаточные функции.

F0=0.31;%m3/c

Fc0=0.18;%m3/c

Qn0=410;%kg/m3

Tn0=335;%K

Tc0=295;%K

Q0=260;%kg/m3

T0=315;%K

Qmax=630;%kg/m3

V=5.4;%m3

A=12.1;%const

b=0.26;%const

c=0.28;% Dzh/kg*gradus

ro=760;%kg/m3

R=8.315; %termod const

a=1*10^16; %const koef

E=54.4; %kDzh

r=10; %kDzh/kg

k0 = a*exp(-E/(R*T0));

k01 = E/(R*T0^2);

kz=b*A*(Fc0);

T1 = V/(F0 + V*k0);

T2 = V/(F0*c*ro + kz);

k1 = (F0*(Q0 - Qn0))/(Q0*(F0 + V*k0));

k2 = F0*Qn0/(Q0*(F0 + V*k0));

k3 = (V*T0*(Qmax - Q0)*k0*k01)/(Q0*(F0+V*k0));

k4 = (kz*(T0 - Tc0))/((F0*c*ro + kz)*T0);

k5 = (kz*(T0 - Tc0))/(2*(F0*c*ro + kz)*T0);

k6 = (r*ro*F0*Qn0)/((F0*c*ro + kz)*T0);

k7 = F0*c*ro*Tn0/((F0*c*ro + kz)*T0);

k8 = kz*Tc0/((F0*c*ro + kz)*T0);

k9 = r*ro*F0*Q0/((F0*c*ro + kz)*T0);

D = [T1*T2 T1+T2 1-k3*k9];

W21=-tf([k5*k3],D)

W22=-tf([k5*T1 k5],D)

8.286*10^(-18 )s^2 + 0.08116 s + 0.1361 1-ая передаточная функция

-2.759*10^(-20) s - 0.0002702

8.286*10^(-18 )s^2 + 0.08116 s + 0.1361 2-ая передаточная функция

k5 = 0.000 27019

k3 = 0.0296

k0 = 9.7944*10^15

где i=m+nи для всех i>n a_i=0, а для всех i>m b_i=0.

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

Для расчета S1, S2 … Si используют численные методы (метод прямоугольников, метод трапеций и др.):

Переход от нормированной передаточной функции к обычной осуществляется путем ее умножения на коэффициент передачи :

(1.5.4)

W(p)=K?W₀(p)(1.5.5)

Для построения переходной характеристики воспользуемся пакетом Matlab:

figure(1)

step(W21);%переходная характеристика

[h t]=step(W21);%число экспериментальных данных

N=length(h);%нахождение коэффициента усиления

k=h(N);%нахождение нормированной переходной характеристики

h0=h/k;%нормированная переходная характеристика

i1=1-h0;

%нахождение первой площади

s1=trapz(t,i1);%первая площадь полученная методом трапеций

a1=s1;%найденный методом площадей первый коэффициент

x=t/s1;%относительное время

s2=s1.^2*trapz(x,(i1.*(-x+1)));%вторая площадь полученная методом трапеций

a2=s2;%найденный методом площадей второй коэффициент

Wrasc1=tf([k],[a2 a1 1]);%передаточная функция с найденными коэффициентами

figure(2)

step(Wrasc1);%реакция на единичное ступенчатое воздействие

[hr tr]=step(Wrasc1);%извлекаем данные

figure(3);

plot(t,h,tr,hr);%сравниваем экспериментальную и найденную характеристики

grid on % нанесение сетки на систему координат

Рисунок 1.5 - Сравнение переходных характеристик аналитической и рассчитанной передаточных функций по первому динамическому каналу

Полученная передаточная функция для 1 - го динамического канала:

figure(4) matlab математический реактор

step(W2);%переходная характеристика

[h t]=step(W2);%моделирование экспериментальных данных;

n=length(h);%число экспериментальных данных

K=h(n)%коэффициент передачи(усиления)

h0=h/K;%нормированная переходная характеристика

S1=trapz(t,1-h0)%первая площадь полученная методом трапеций

x=t/S1;%относительное время

S2=S1^2*trapz(x,(1-h0).*(1-x))%вторая площадь полученная методом трапеций

S3=S1^3*trapz(x,(1-h0).*(1-2*x+(x.^2)/2))%третья площадь полученная методом трапеций

%находим коэффициенты передаточной функции из системы линейных уравнений

A=[1 0 -1; 0 1 -S1; 0 0 -S2];

b=[S1;S2;S3];

x=A\b%переходная характеристика найденной функции

a1=x(1)%коэффициенты передаточной функции

a2=x(2)

b1=x(3)

Wrasc2=tf([b1 1],[a2 a1 1]);%передаточная функция с найденными коэффициентами

Wraschetnoe2=K*Wrasc2

figure(5)

step(Wraschetnoe2);%реакция на единичное ступенчатое воздействие

[hrasc trasc]=step(Wraschetnoe2);%извлекаем данные

figure(6);

plot(t,h,trasc,hrasc);%сравниваем экспериментальную и найденную характеристики

grid on % нанесение сетки на систему координат

Рисунок 1.6 -Сравнение переходных характеристик аналитической и рассчитанной передаточных функций по второму динамическому каналу

Полученная передаточная функция для 2 - го динамического канала:

-0.00102 s - 0.001479

0.3087 s^2 + 1.131 s + 1

Как видно из графиков переходные характеристики аналитической и рассчитанной переходных функций практически полностью совпадают.

Заключение

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

Все расчеты и моделирование были выполнены в пакете Matlab 7.0