Теперь рассмотрим еще одну программу по той же задаче.

function DATA

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

global xa0 TR t_a t0 T0 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. Начальное приближение при расчетах для правой границы изменения времени пребывания в реакторев  реакторе (с)

t0=0.5;

%5.Начальное приближение при расчетах для температуры в реакторе (C)

T0=500;

end

 

function dy = difpravreactor(t,y)

%Программныйкодфайлаdifpravreactor.m - расчетправыхчастейсистемы%дифференциальныхуравненийматематическогоописанияпроцесса%Программавключаетследующие%файлы:GLAV_model3_batch_reactor.m+DATA.m+difpravreactor.m+polyno3.m+REPORT.m%ПрограммарасчетаоптимальноговыходапродуктаРв%изотермическом периодическом реакторе со стехиометрической схемой%реакцииA - P - S

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 phi_p=model_pfr(regr)

%Программныйкодфайла model+pfr.m - расчетвыходацелевогопродуктаСпри%известнойтемпературеТивременипребывания t_b вреактореидеального%вытеснения%Программавключаетследующие%файлы:GLAV_optim5_batch_reactor_T_t.m+DATA.m+model_pfr+difpravreactor.m+REPORT.m

%Программа расчета максимального выхода продукта C в%изотермическом периодическом реакторе со стехиометрической схемой%реакцийA - B (k1); A - D (k2); A - D (k3); B - C (k4); 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;

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);

tk=t(n);yr(1)=y(n,1);yr(2)=y(n,2);yr(3)=y(n,3); yr(4)=y(n,4);

phi_p=-yr(3)/xa0;

end

 

functionGLAV_optim5_batch_reactor_T_t

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

clc;

clear all;

close all;

global xa0 t_a t_b t y n T k1 k2 k3 k4 k5 T0 t0 TR T_opt phi_p_max A1 E1 A2 E2 A3 E3 A4 E4 A5 E5 yr1 yr2 yr3 yr4;

DATA;

regr=[T0,t0];

disp('Информация об итерационном процессе расчета');

options=optimset('Display','iter');

%Реализация алгоритма поиска минимума функции двух переменных: T и t -

%методом многомерной оптимизации

[regrcalc,phi_p]=fminsearch('model_pfr',regr,options);

%Реализация алгоритма поиска минимума функции двух переменных: T и t -

%методом нелинейного программирования

phi_p_max=-phi_p;

T_opt=regrcalc(1);

t_b=regrcalc(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);

tk=t(n);yak=y(n,1);ybk=y(n,2);yck=y(n,3);

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

REPORT;

end

 

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

global xa0 t n t_a t_b t0 T0 TR T_opt phi_p_max A1 E1 A2 E2 A3 E3 A4 E4...

A5 E5 yr1 yr2 yr3 yr4;

disp('ПРОГРАММА ОПРЕДЕЛЕНИЕМ МАКСИМАЛЬНОГО ВЫХОДА ЦЕЛЕВОГО ПРОДУКТА "С" ')

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

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

disp('Программавключаетследующиефайлы:GLAV_optim5_batch_reactor_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') ' час']);

disp(['5. Начальное приближение при расчетах для правой границы изменения времени пребывания в реакторе(t0) = ' num2str(t0,'%10.2f') ' час']);

disp(['6. Начальное приближение при расчетах для температуры в реакторе(T0) = ' num2str(T0,'%10.2f') ' час']);

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

disp ('1.Оптимальные параметры выхода целевого продукта С ');

disp(' T_opt   t_opt      phi_p_max ')

disp(sprintf('%10.3f\t %10.5f\t %10.5f',T_opt,t_b,phi_p_max));

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

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

for i=1:n

y11= yr1(i);

y22=yr2(i);

y33=yr3(i);

y44=yr4(i);

tt=t(i);

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

end

plot(t,yr1,'k-.',t,yr2,'k--',t,yr3,'-',t,yr4,'*');

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

xlabel('t - C');ylabel('Конц. А ("-."), B("- -"), C("-"), D("...") - мольн. дол.');

gridon;   

end