Расчет профилей изменения концентраций во времени в изотермическом периодическом реакторе.

Программа расчета профилей изменения концентраций во времени визотермическом периодическом реакторе со стехиометрической схемойреакций (5 реакций): 1.) A -(k1)-> B; 2.) A -(k2)-> D; 3.) A -(k3)-> D;

 4.) B -(k4)-> C; 5.) C -(k5)-> D.

Целью является определение максимального выхода целевого продукта "С" в изотермическом реакторе периодического действия

functionDATA

%Программный код файла DATA.m - задание исходной информации для расчетов;%Программа включает следующие%файлы:GLAV_optim5_batch_reactor_MODEL_T_t.m+DATA.m+model_pfr+difpravreactor.m+REPORT.m%Программа расчета профилей изменения концентраций во времени в%изотермическом периодическом реакторе со стехиометрической схемой%реакций 5 реакций: 1.)A - B (k1);2.) A - D (k2);3.) A - D (k3); 4.)B - C (k4); 5.) C - D (k5);

global xa0 TR t_a t_b T_a T_b A1 E1 A2 E2 A3 E3 A4 E4 A5 E5;

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

xa0=1;

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

%а) Первая реакция - k1:

%a.1) предэкспоненциальный множитель (c^(-1))

A1=1.02;

%a.2) экспоненциальный множитель (С)

E1=808.08;

%b) Вторая реакция - k2:

%b.1) предэкспоненциальный множитель (c^(-1))

A2=0.93;

%b.2) экспоненциальный множитель (С)

E2=707.07;

%с) Третья реакция - k3:

%с.1) предэкспоненциальный множитель (c^(-1))

A3=0.386;

%с.2) экспоненциальный множитель (С)

E3=757.57;

%d) Четвертая реакция - k4:

%d.1) предэкспоненциальный множитель (c^(-1))

A4=3.28;

%d.2) экспоненциальный множитель (С)

E4=505.05;

%e) Пятая реакция - k5:

%e.1) предэкспоненциальный множитель (c^(-1))

A5=1.084;

%e.2) экспоненциальный множитель (С)

E5=757.57;

%f.Базовая температура (С)

TR=658;

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

t_a=0;

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

t_b=4;

%5.Нижняя граница изменения температуры в реакторе (C)

T_a=500;

%6.Верхняя граница изменения температуры в реакторе (C)

T_b=1000;

end

 

function dy = difpravreactor(t,y)

%Программный код файла difpravreactor.m - расчет правых частей системы%дифференциальных уравнений математического описания процесса%Программа включает следующие%файлы:GLAV_optim5_batch_reactor_MODEL_T_t.m+DATA.m+model_pfr+difpravreactor.m+REPORT.m%Программа расчета профилей изменения концентраций во времени в%изотермическом периодическом реакторе со стехиометрической схемой%реакций 5 реакций: 1.)A - B (k1);2.) A - D (k2);3.) A - D (k3); 4.)B - C (k4); 5.) C - D (k5);

global k1 k2 k3 k4 k5;

dy=zeros(4,1);

dy(1)=-(k1+k2+k3)*y(1);

dy(2)=k1*y(1)-k4*y(2);

dy(3)=k4*y(2)-k5*y(3);

dy(4)=(k2+k3)*y(1)+k5*y(3);

end

 

function yr=model_pfr(regr)

%Программный код файла model_pfr.m - расчет выходных концентраций продуктов реакции при%известной температуре Т и времени пребывания t_b в изотермическом реакторе%периодического действия%Программа включает следующие%файлы:GLAV_optim5_batch_reactor_MODEL_T_t.m+DATA.m+model_pfr+difpravreactor.m+REPORT.m%Программа расчета профилей изменения концентраций во времени в%изотермическом периодическом реакторе со стехиометрической схемой%реакций 5 реакций: 1.)A - B (k1);2.) A - D (k2);3.) A - D (k3); 4.)B - C (k4); 5.) C - D (k5);

global xa0 t_a t_b t y n T k1 k2 k3 k4 k5 TR A1 E1 A2 E2 A3 E3 A4 E4 A5 E5;

global yr1 yr2 yr3 yr4 ii xa xb xc xd yfa yfb yfc yfd TT nn tt;

yr=zeros(4,1);

T=regr(1); 

t_b=regr(2);

k1=A1*exp(-E1*(1/T-1/TR));

k2=A2*exp(-E2*(1/T-1/TR));

k3=A3*exp(-E3*(1/T-1/TR));

k4=A4*exp(-E4*(1/T-1/TR));

k5=A5*exp(-E5*(1/T-1/TR));

[t,y]=ode45(@difpravreactor,[t_a t_b],[xa0 0 0 0]);

n=length(t);

for i=1:n

yr1(i)=y(i,1);

yr2(i)=y(i,2);

yr3(i)=y(i,3);

yr4(i)=y(i,4);

end

if ii==1

xa=zeros(n,6);xb=xa;xc=xa;xd=xa;tt=xa;

end

for j=1:n

tt(j,ii)=t(j);

xa(j,ii)=yr1(j);

xb(j,ii)=yr2(j);

xc(j,ii)=yr3(j);

xd(j,ii)=yr4(j);

end

nn(ii)=n;

TT(ii)=T;

yfa(ii)=xa(n,ii);

yfb(ii)=xb(n,ii);

yfc(ii)=xc(n,ii);

yfd(ii)=xd(n,ii);

yr(1,1)=yfa(ii);

yr(2,1)=yfb(ii);

yr(3,1)=yfc(ii);

yr(4,1)=yfd(ii);

end

 

functionGLAV_optim5_batch_reactor_MODEL_T_t

%Программныйкодфайла GLAV_optim5_batch_reactor_MODEL_T_t.m - основнаяуправляющаяпрограмма%Программавключаетследующие%файлы:GLAV_optim5_batch_reactor_MODEL_T_t.m+DATA.m+model_pfr+difpravreactor.m+REPORT.m%Программарасчетапрофилейизмененияконцентрацийвовременив%изотермическомпериодическомреакторесостехиометрическойсхемой%реакций 5 реакций: 1.)A - B (k1);2.) A - D (k2);3.) A - D (k3); 4.)B - C (k4); 5.) C - D (k5);

clc;

clear all;

close all;

global phi_c t_b t T T_a T_b xa0 yr3 ii;

DATA;

ii=0;

for T=T_a:100:T_b

ii=ii+1; 

regr=[T,t_b];

x= model_pfr(regr);

switch ii

case 1

      ssr='+';

case 2

      ssr='*';

case 3

      ssr=':';

case 4

      ssr='-';

case 5

  ssr='-.';

otherwise

     ssr='--';

end

phi_c=yr3/xa0;

hold on;    

plot(t,phi_c,ssr);

title('Выход целевого продукта С в периодическом реакторе при 6 температурах');

xlabel('t - с');ylabel('Т (С): 500-"х", 600-"*", 700-"..", 800-"-",900-"-.-",1000-"- -"-мол.дол.');

% axis([0 t_b 0 0.21]);

grid on;   

end

REPORT

end

 

functionREPORT

%ПрограммныйкодфайлаREPORT.m - отчетоработепрограммы%Программавключаетследующие%файлы:GLAV_optim5_batch_reactor_MODEL_T_t.m+DATA.m+model_pfr+difpravreactor.m+REPORT.m%Программарасчетапрофилейизмененияконцентрацийвовременив%изотермическомпериодическомреакторесостехиометрическойсхемой%реакций 5 реакций: 1.)A - B (k1);2.) A - D (k2);3.) A - D (k3); 4.)B - C (k4); 5.) C - D (k5);

global xa0 t n t_a t_b T_a T_b T0 TR A1 E1 A2 E2 A3 E3 A4 E4 A5 E5 TT;

global yfa yfb yfc yfd ii xa xb xc xd nn tt;

disp('ПРОГРАММА МОДЕЛИРОВАНИЯ ПРОЦЕССА ')

disp('В ИЗОТЕРМИЧЕСКОМ РЕАКТОРЕ ПЕРИОДИЧЕСКОГО ДЕЙСТВИЯ')

disp('С 5 РЕАКЦИЯМИ: A - B (k1); A - D (k2); A - D (k3); B - C (k4); C - D (k5) ');

disp('Программавключаетследующиефайлы:GLAV_optim5_batch_reactor_MODEL_T_t.m+DATA.m+model_pfr+difpravreactor.m+REPORT.m ');

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

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

disp(['2.1.Предэкспоненциальный множитель первой элементарной реакции (A1) = ' num2str(A1,'%10.4f') ' c^(-1)']);

disp(['2.2.Экспоненциальный множитель первой элементарной реакции (E1) = ' num2str(E1,'%10.4f') ' C']);

disp(['3.1.Предэкспоненциальный множитель второй элементарной реакции (A2) = ' num2str(A2,'%10.4f') ' c^(-1)']);

disp(['3.2.Экспоненциальный множитель второй элементарной реакции (E2) = ' num2str(E2,'%10.4f') ' C']);

disp(['4.1.Предэкспоненциальный множитель третьей элементарной реакции (A3) = ' num2str(A3,'%10.4f') ' c^(-1)']);

disp(['2.2.Экспоненциальный множитель третьей элементарной реакции (E3) = ' num2str(E3,'%10.4f') ' C']);

disp(['5.1.Предэкспоненциальный множитель четвертой элементарной реакции (A4) = ' num2str(A4,'%10.4f') ' c^(-1)']);

disp(['5.2.Экспоненциальный множитель четвертой элементарной реакции (E4) = ' num2str(E4,'%10.4f') ' C']);

disp(['2.1.Предэкспоненциальный множитель пятой элементарной реакции (A5) = ' num2str(A5,'%10.4f') ' c^(-1)']);

disp(['2.2.Экспоненциальный множитель пятой элементарной реакции (E5) = ' num2str(E5,'%10.4f') ' C']);

disp(['3.Базовая температура для выражения констант скоростей реакций (TR) - k=A*exp(-E*(1/T-1/TR)) = ' num2str(TR,'%10.2f') ' C']);

disp(['4.Левая граница изменения времени пребывания в реакторе (t_a) = ' num2str(t_a,'%10.2f') ' c ']);

disp(['5.Правая граница изменения времени пребыванияв реакторе (t_b) = ' num2str(t_b,'%10.2f') ' c ']);

disp(['5. Нижняя граница изменения температуры в реакторе(T_a) = ' num2str(T_a,'%10.2f') ' C']);

disp(['6. Верхняя граница изменения температуры в реакторе(T_b) = ' num2str(T_b,'%10.2f') ' C ']);

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

forj=1:6

hfig(j)=figure;

end

for i=1:ii

TV=TT(i);

TVS=num2str(TV);

disp(['1. Температура в реакторе(T) = ' num2str(TV,'%10.2f') ' С']);

disp ('2.Параметры выхода продуктов ');

disp(' A     B      C           D ')

yfaa=yfa(i);yfbb=yfb(i);yfcc=yfc(i);yfdd=yfd(i);

disp(sprintf('%10.3f\t %10.5f\t %10.5f\t %10.5f',yfaa,yfbb,yfcc,yfdd))

disp ('2.Изменение концентрации продуктов во времени в периодическом реакторе');

disp(' t(c)        xa(мол.д.)    xb(мол.д.) xc(мол.д.) xd(мол.д.) '); 

n=nn(i);

for j=1:n

t=tt(j,i);

y11= xa(j,i);

y22=xb(j,i);

y33=xc(j,i);

y44=xd(j,i);

tp(j)=t;y1p(j)=y11;y2p(j)=y22;y3p(j)=y33;y4p(j)=y44;

disp(sprintf(' %10.5f\t %10.5f\t %10.5f\t %10.5f\t %10.5f',t,y11,y22,y33,y44));

end

figure(hfig(i));

 plot(tp,y1p,'-.',tp,y2p,'--',tp,y3p,'-',tp,y4p,':');

title(['Результат моделирования периодического реактора с 5 реакциями при Т =' TVS ' C']);

xlabel('t - с');ylabel('Составы: А ("-."), B("- -"), C(- "-"), D("...")-мол.дол.');

grid on;   

end

end