Оптимизация равновесных экзотермических реакций

Рассмотрим химическую реакцию с целевым продуктом P, проходящую по схеме:

Скорость реакции по целевому продукту

выражается по закону действующих масс:

Необходимо определить

температуру T проведения реакции, при которой:

Необходимое условие экстремума:

 

 

 

В состоянии равновесия скорость реакции W равна 0:

Получаем связь равновесной и оптимальной температур проведения реакции:
Откуда:
После логарифмирования получаем:
Откуда следует:

Теперь рассмотрим программу по следующей задаче:

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

Схема №1	Схема №2

Система нелинейных уравнений (СНУ).В общем случае систему нелинейных уравнений можно записать как: или Решением СНУ является такой вектор при подстановке которого в систему последняя обращается в тождество.

Метод Ньютона-Рафсона.Пусть известно некоторое приближение к решению Запишем исходную систему в виде где Разложим функцию в ряд Тейлора и ограничимся линейными членами.

Это система линейных уравнений относительно

Матрица Якоби Тогда а новое приближение к решению СНУ будетиметь вид: или Условием окончания итерационного процесса является выполнения неравенства

Методом Ньютона-Рафсона воспользуемся для решения поставленной задачи.

ФункцияDATA1 содержит исходные данные для процесса, происходящего по вышеописанной Схеме 1.

aA -(k 1)-> pP + bB -(k 2)-> S

function DATA1

%Программный код файла DATA.m - задание исходной информации для расчетов %Программа расчета оптимального времени пребывания (tau - час) в%изотермическом реакторе идеального смешения со стехиометрической схемой%реакции

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%схема №1

% aA -(k1)-> pP + bB -(k2)-> S

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global xa0 xb0 k1 k2 k3 tau_a tau_b ii aa bb pp;

%Концентрация реагента А в долях в реакции A-P-S (мольные доли)

xa0=0.31;

%Константы скоростей реакций (час^(-1))

k1=0.88;

k2=0.68;

%3. Левая граница поиска оптимального времени пребывания в реакторе (час)

tau_a=5;

%4. Правая граница поиска оптимального времени пребывания в реакторе (час)

tau_b=20;

ii=0;

%5. Коэффициенты в схеме реакции a, b, p

aa=2;

bb=1;

pp=3;

end

ФункцияDATA2 содержит исходные данные для процесса, происходящего по вышеописанной Схеме 2.

 aA =(<-k2,k1->)= pP

 P + bB -(k3)-> sS

function DATA2

%Программный код файла DATA.m - задание исходной информации для расчетов; %Программа расчета оптимального времени пребывания (tau - час) в %изотермическом реакторе идеального смешения

global xa0 k1 k2 k3 tau_a tau_b ii aa bb pp NS xb0;

%Схемареакции

NS=2;%1

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%схема №2

% aA =(<-k2,k1->)= pP

% P + bB -(k3)-> sS

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Концентрация реагента в долях в реакции (мольные доли)

xa0=0.213;

%Константы скоростей реакций (час^(-1))

k1=0.75;

k2=0.01;

   k3=0.333;

   xb0=0.92;

%Константы схемы реакции

aa=3;

bb=3;

pp=3;

%Левая граница поиска оптимального времени пребывания в реакторе (час)

tau_a=5;

%Правая граница поиска оптимального времени пребывания в реакторе (час)

tau_b=50;

ii=0;

end

Функцияmodel1_stat_tau содержит систему нелинейных уравнений (СНУ), матрицу Якоби из частных производных и решает СНУ методом Ньютона-Рафсона. При этом найденное решение соответствует определенной продолжительности (определенному моменту времени в течение продолжительного процесса) реакции по Схеме 1. Следовательно, для построения зависимости концентраций от времени требуется циклически (многократно) вызывать функциюmodel1_stat_tau, задавая для нее время (аргумент функции) в пределах от начального до конечного с некоторым шагом. Это позволит накопить данные (точки графика) для построения кривой зависимости для каждой концентрации и, в целом, их совокупности, в течение всего заданного периода времени осуществления процесса в реакторе.

Система нелинейных уравнений для Схемы 1.

Описывается в тексте функцииmodel1_stat_tau следующими функциями – по одной на каждое уравнение системы.

function f1=fun1(xVector, tau)

   xa=xVector(1);

   f1=xa0-xa-aa*tau*k1*(xa^aa);

end% fun1

function f2=fun2(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   f2=-xp+tau*pp*(k1*(xa^aa)-k2*(xp^pp)*(xb^bb));

end% fun2

function f3=fun3(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   f3=-xb+tau*bb*(k1*(xa^aa)-k2*(xp^pp)*(xb^bb));

end% fun3

Затем описывается матрица Якоби

function x= model1_stat_tau(tau)

%Программный код файла model1_stat_tau.m - расчет выходных концентраций %продуктов из реактора %Программа расчета оптимального времени пребывания (tau - час) в %изотермическом реакторе идеального смешения со стехиометрической схемой

%реакции

МЕТОД НЬЮТОНА-РАФСОНА – РЕШЕНИЕ СИСТЕМЫ НЕЛИНЕЙНЫХ УРАВНЕНИЙ

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%схема №1

% aA -(k1)-> pP + bB -(k2)-> S

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global xa0 k1 k2 aa bb pp xaNR xbNR xpNR ii;

 

eps=0.0001;

 

function f1=fun1(xVector, tau)

   xa=xVector(1);

   f1=xa0-xa-aa*tau*k1*(xa^aa);

end% fun1

 

function f2=fun2(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   f2=-xp+tau*pp*(k1*(xa^aa)-k2*(xp^pp)*(xb^bb));

end% fun2

 

function f3=fun3(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   f3=-xb+tau*bb*(k1*(xa^aa)-k2*(xp^pp)*(xb^bb));

end% fun3

 

function d11=dfun1_dxa(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d11=-1-tau*k1*(aa^2)*(xa^(aa-1));

end% dfun1_dxa

 

function d21=dfun2_dxa(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d21=(tau*pp*k1*aa*xa^(aa-1));

end% dfun2_dxa

 

function d22=dfun2_dxb(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d22=-(tau*pp*k2*(xp^pp)*bb*(xb^(bb-1)));

end% dfun2_dxb

 

function d23=dfun2_dxp(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d23=-1-(tau*(pp^2)*k2*(xb^bb)*(xp^(pp-1)));

end% dfun3_dxp

 

function d31=dfun3_dxa(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d31=(tau*bb*k1*aa*(xa^(aa-1)));

end% dfun3_dxa

 

function d32=dfun3_dxb(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

      d32=-1-(tau*(bb^2)*k2*(xp^pp)*(xb^(bb-1)));

end% dfun3_dxb

 

function d33=dfun3_dxp(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d33=-(tau*bb*k2*(xb^bb)*pp*xp^(pp-1));

end% dfun3_dxp

 

%%% NEWTON - RAFSON %%%

 

if ii==0

   xVector=[xa0-0.01;0.01;0.01];

else

   xVector=[xaNR;xbNR;xpNR];

end% if

 

for i=1:100

 

   f1=fun1(xVector, tau);

   f2=fun2(xVector, tau);

   f3=fun3(xVector, tau);

   fVector=[f1;f2;f3];

 

   d11=dfun1_dxa(xVector, tau);

   d21=dfun2_dxa(xVector, tau);

   d22=dfun2_dxb(xVector, tau);

   d23=dfun2_dxp(xVector, tau);

   d31=dfun3_dxa(xVector, tau);

   d32=dfun3_dxb(xVector, tau);

   d33=dfun3_dxp(xVector, tau);

   Jac=[d11, 0, 0; d21, d22, d23; d31, d32, d33];

 

   JObr=Jac^(-1);

 

   hVector=JObr*fVector;

 

if norm(hVector)<eps

break;

end% if

 

 

   xVector=xVector-hVector;

 

end% for i

 

xaNR=xVector(1);

xbNR=xVector(2);

xpNR=xVector(3);

ii=ii+1;  

 

x=xVector;

 

end

Функцияmodel2_stat_tau содержит систему нелинейных уравнений (СНУ), матрицу Якоби из частных производных и решает СНУ методом Ньютона-Рафсона. При этом найденное решение соответствует определенной продолжительности (определенному моменту времени в течение продолжительного процесса) реакции по Схеме 2. Следовательно, для построения зависимости концентраций от времени требуется циклически (многократно) вызывать функциюmodel2_stat_tau, задавая для нее время (аргумент функции) в пределах от начального до конечного с некоторым шагом. Это позволит накопить данные (точки графика) для построения кривой зависимости для каждой концентрации и, в целом, их совокупности, в течение всего заданного периода времени осуществления процесса в реакторе.

Система нелинейных уравнений для Схемы 2.

 

 

Описывается в тексте функцииmodel2_stat_tau следующими функциями – по одной на каждое уравнение системы.

 

function f1=fun1(xVector, tau)

   xa=xVector(1);

   xp=xVector(3);

   f1=xa0-xa+tau*aa*(k2*(xp^pp)-k1*(xa^(aa)));

end% fun1   

 

function f2=fun2(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

f2=-xp+tau*(pp*k1*(xa^aa)-k2*pp*(xp^pp)-k3*xp*(xb^bb));

end% fun2

 

function f3=fun3(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   f3=xb0-xb-tau*bb*k3*xp*(xb^bb);

end% fun3

 

Затем описывается матрица Якоби

function x= model2_stat_tau(tau)

%Программный код файла model2_stat_tau.m - расчет выходных концентраций %продуктов из реактора %Программа расчета оптимального времени пребывания (tau - час) в %изотермическом реакторе идеального смешения со стехиометрической схемой %реакции

МЕТОД НЬЮТОНА-РАФСОНА – РЕШЕНИЕ СИСТЕМЫ НЕЛИНЕЙНЫХ УРАВНЕНИЙ

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%схема №2

% aA =(<-k2,k1->)= pP

% P + bB -(k3)-> sS

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global xa0 xb0 k1 k2 k3 aa bb pp xaNR xbNR xpNR ii;

 

eps=0.0001;

 

function f1=fun1(xVector, tau)

   xa=xVector(1);

   xp=xVector(3);

   f1=xa0-xa+tau*aa*(k2*(xp^pp)-k1*(xa^(aa)));

end% fun1   

 

function f2=fun2(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

f2=-xp+tau*(pp*k1*(xa^aa)-k2*pp*(xp^pp)-k3*xp*(xb^bb));

end% fun2

 

function f3=fun3(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   f3=xb0-xb-tau*bb*k3*xp*(xb^bb);

end% fun3

 

function d11=dfun1_dxa(xVector, tau)

   xa=xVector(1);

   d11=-1-tau*k1*(aa^2)*(xa^(aa-1));

end% dfun1_dxa

 

function d13=dfun1_dxp(xVector, tau)

   xa=xVector(1);

   xp=xVector(3);

   d13=tau*pp*k2*aa*(xp^(pp-1));

end% dfun1_dxp

 

function d21=dfun2_dxa(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d21=tau*pp*k1*aa*(xa^(aa-1));

end% dfun2_dxa

 

function d22=dfun2_dxb(xVector, tau)

    xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d22=-tau*bb*k3*xp*(xb^(bb-1));

end% dfun2_dxb

 

function d23=dfun2_dxp(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d23=-1-tau*((pp^2)*k2*(xp^(pp-1))+k3*(xb^(bb)));

end% dfun2_dxp

 

function d32=dfun3_dxb(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d32=-1-tau*(bb^2)*k3*xp*(xb^(bb-1));

end% dfun3_dxb

 

function d33=dfun3_dxp(xVector, tau)

   xa=xVector(1);

   xb=xVector(2);

   xp=xVector(3);

   d33=-tau*bb*k3*(xb^bb);

end% dfun3_dxp

 

%%% NEWTON - RAFSON %%%

 

if ii==0

   xVector=[xa0-0.01;xb0-0.01;0.01];

else

   xVector=[xaNR;xbNR;xpNR];

end% if

 

for i=1:100

 

   f1=fun1(xVector, tau);

   f2=fun2(xVector, tau);

   f3=fun3(xVector, tau);

   fVector=[f1;f2;f3];

 

   d11=dfun1_dxa(xVector, tau);

   d13=dfun1_dxp(xVector, tau);

   d21=dfun2_dxa(xVector, tau);

   d22=dfun2_dxb(xVector, tau);

   d23=dfun2_dxp(xVector, tau);

   d32=dfun3_dxb(xVector, tau);

   d33=dfun3_dxp(xVector, tau);

   Jac=[d11, 0, d13; d21, d22, d23; 0, d32, d33];

 

   JObr=Jac^(-1);

 

   hVector=JObr*fVector;

 

if norm(hVector)<eps

break;

end% if

 

   xVector=xVector-hVector;

end% for i

 

xaNR=xVector(1);

xbNR=xVector(2);

xpNR=xVector(3);

ii=ii+1;  

 

x=xVector;

 

end

 

function GLAV_model1_grafik

%Программныйкодфайла GLAV_ model1_grafik.m - главнаяуправляющаяпрограмма %Программавключаетследующие %файлы: GLAV_model1_grafik.m+ DATA1.m+ model1_stat_tau.m+ REPORT.m;%Программарасчетаоптимальноговременипребывания (tau - час) в %изотермическомреактореидеальногосмешениясостехиометрическойсхемой %реакции

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%схема №1

% aA -(k1)-> pP + bB -(k2)-> S

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc;

clear all;

close all;

global t xa0 xb0 xa xb xp tau_a tau_b xaNR xbNR xpNR ii phi_p;

i=0;

DATA1;

for tau=tau_a:0.1:tau_b

 i=i+1;  

x= model1_stat_tau(tau);

t(i)=tau;

xa(i)=x(1);

xb(i)=x(2);

xp(i)=x(3);

phi_p(i)=xp(i)/xa0;

end

REPORT;

end

 

function GLAV_model2_grafik

%Программныйкодфайла GLAV_model2_grafik.m - главнаяуправляющаяпрограмма%Программавключаетследующие%файлы:GLAV_model2_grafik.m+DATA2.m+model2_stat_tau.m+REPORT.m;%Программарасчетаоптимальноговременипребывания (tau - час) в%изотермическомреактореидеальногосмешениясостехиометрическойсхемой%реакции

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%схема №2

% aA =(<-k2,k1->)= pP

% P + bB -(k3)-> sS

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc;

clear all;

close all;

global t xa0 xb0 xa xb xp tau_a tau_b xaNR xbNR xpNR ii phi_p;

i=0;

DATA2;

for tau=tau_a:0.1:tau_b

 i=i+1;  

x= model2_stat_tau(tau);

t(i)=tau;

xa(i)=x(1);

xb(i)=x(2);

xp(i)=x(3);

phi_p(i)=xp(i)/xa0;

end

REPORT;

end

 

function REPORT;

%Программный код файла REPORT.m - отчет о работе программы

%Программа расчета оптимального времени пребывания (tau - час) в%изотермическом реакторе идеального смешения со стехиометрической схемой

%реакции A - P - S

global xa0 k1 k2 k3 t xa xb xp tau_a tau_b phi_p;

disp('ПРОГРАММА МОДЕЛИРОВАНИЯ ИЗОТЕРМИЧЕСКОГО РЕАКТОРА ИДЕАЛЬНОГО ПЕРЕМЕШИВАНИЯ ');

disp('ИСХОДНЫЕ ДАННЫЕ ');

disp(['1.Концентрация реагента A на входе в реактор (xa0) = ' num2str(xa0,'%10.2f') ' мольные доли']);

disp(['2.Константа скорости первой элементарной реакции (k1) = ' num2str(k1,'%10.2f') ' час^(-1)']);

disp(['3.Константа скорости второй элементарной реакции (k2) = ' num2str(k2,'%10.2f') ' час^(-1)']);

disp(['4.Константа скорости третьей элементарной реакции (k3) = ' num2str(k3,'%10.2f') ' час^(-1)']);

disp(['5.Левая граница интервала поиска оптимального времени пребывания в реакторе (tau_a) = ' num2str(tau_a,'%10.2f') ' час']);

disp(['6. Правая граница интервала поиска оптимального времени пребывания в реакторе(tau_b) = ' num2str(tau_b,'%10.2f') ' час']);

disp ('РЕЗУЛЬТАТЫ РАСЧЕТОВ');

disp ('1.Концентрации продуктов на выходе из реактора');

disp(' tau(час) xa(мол.д.) xb(мол.д.) xp(мол.д.)'); 

i=0;

for i=1:length(t)

   tt=t(i);xaa=xa(i);xbb=xb(i); xpp=xp(i);

disp(sprintf('%10.3f\t %10.5f\t %10.5f\t %10.5f',tt,xaa,xbb,xpp));

end

%Построение графиков;

hfig(1)=figure;

hfig(2)=figure;

hfig(3)=figure;

hfig(4)=figure;

figure(hfig(1));

plot(t,xa,'r-');

title('График зависимости выхода продукта A от времени пребывания в реакторе');

xlabel('tau - час');ylabel('Выходные концентрации: xa - мольные доли');

grid on;

figure(hfig(2));

plot(t,xb,'k-');

title('График зависимости выхода продукта B от времени пребывания в реакторе');

xlabel('tau - час');ylabel('Выходные концентрации: xb - мольные доли');

grid on;

figure(hfig(3));

plot(t,xp,'g-');

title('График зависимости выхода продукта P от времени пребывания в реакторе');

xlabel('tau - час');ylabel('Выходные концентрации: xp - мольные доли');

grid on;

figure(hfig(4));

plot(t,phi_p,'b-');

title('График зависимости выхода продукта PHI_P от времени пребывания в реакторе');

xlabel('tau - час');ylabel('Выходные концентрации: phi_p - мольные доли');

gridon;

end