Определение кинетических коэффициентов гомогенной реакции

Определить кинетические коэффициенты гомогенной реакции в периодическом изотермическом реакторе с мешалкой, использовав в качестве квадратичного критерия оптимальности рассогласование между расчетными и экспериментальными значениями концентрации компонентов А и Р в диапазоне времени от 0 до 10 ч.

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

Порядок обеих стадий реакции - первый. Начальные приближения для констант скоростей равны:k₁ = 0,01 час⁻¹; k₂ = 1 час⁻¹.

 

function DATA

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

global xa0 k10 k20 t_a t_b aaexp;

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

xa0=1;

%Экспериментальные данные: для каждой опытной точки построчно - время (t-%час),концентрация исходного компонента (А - мольн.доли), концентрация%целевого компонента (P - мольн.доли);

aaexp=[0.000,1.000,0.000; 0.104,0.969,0.037;…

0.523,0.855,0.141;1.023,0.736,0.251;…

 3.773,0.322,0.545; 4.023,0.299,0.554;…

5.523,0.191,0.577; 6.273,0.153,0.573;…

 9.023,0.067,0.50833;10.000,0.0498,0.477];

%Начальные приближения констант скоростей реакций (час^(-1))

k10=0.01;k20=1;

%3. Левая граница изменения времени в реакторе при определении констант скоростей реакций(час)

t_a=0;

%4. Правая граница изменения времени в реакторе при определении констант скоростей реакций(час)

t_b=10;

end

 

functionCr = Criterium(aregr)

%Программныйкодфайла Criterium.m - расчетквадратичногокритериярассогласованиямежду %опытнымиирасчетнымиконцентрациямикомпонентовреакции;%Программавключаетследующие %файлы: GLAV_ident4_batch_reactor.m+ DATA.m+ difpravreactor.m+Criterium.m+REPORT.m

%Программа определения кинетических констант по экспериментальным данным в %изотермическом периодическом реакторе со стехиометрической схемой%реакции A - P - S

global tt t_a t_b xa0 k1 k2 yar ypr xae n xpe;

k1=aregr(1);

k2=aregr(2);

yar(1)=xa0;ypr(1)=0;

for i=2:n

t_b=tt(i); 

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

nn=length(t);

yar(i)=y(nn,1);

ypr(i)=y(nn,2);

end

Cr=0;

for i=1:10

Cr=Cr+(yar(i)-xae(i))^2+(ypr(i)-xpe(i))^2;

end

end

 

function dy = difpravreactor(t,y)

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

%Программа определения кинетических констант по экспериментальным данным в%изотермическом периодическом реакторе со стехиометрической схемой%реакции A - P - S

global k1 k2

dy=zeros(2,1);

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

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

end

 

functionGLAV_ident4_batch_reactor

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

clc;

clear all;

close all;

global n aaexp tt xae xpe k10 k20 crit;

DATA;

%pause

n=length(aaexp);

for i=1:n

tt(i)=aaexp(i,1);xae(i)=aaexp(i,2);xpe(i)=aaexp(i,3);

end

aregr=[k10 k20];

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

[aregr,crit]=fminsearch('Criterium',aregr,options);

k1=aregr(1);k2=aregr(2);

REPORT;

end

 

function REPORT

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

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

global xa0 k1 k2 t n y t_a t_b tt yar ypr xae xpe k10 k20 crit;

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

disp(' РЕАКЦИЕЙ: A - P - S ');

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

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

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

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

disp(['3.Начальное приближение константы скорости второй элементарной реакции (k20) = ' num2str(k20,'%10.2f') ' час^(-1)']);

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

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

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

disp(['1.Минимальное значения критерия рассогласования опытных и расчетных значения концентраций (crit) = ' num2str(crit,'%10.6f') ' мольные доли']);

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

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

disp ('4.Сравнение изменения расчетных и опытных концентраций компонентов во времени в периодическом реакторе');

disp(' t(час)        xa(мол.д.) xaexp(мол.д.) xp(мол.д.) xpexp(мол.д.)'); 

for i=1:n

y1(i)=yar(i); 

y2(i)=ypr(i);

y11= y1(i);

y22=y2(i);

y11e=xae(i);

y22e=xpe(i);

ttt=tt(i);

disp(sprintf('%10.3f\t %10.5f\t %10.5f\t %10.5f\t %10.5f',ttt,y11,y11e,y22,y22e));

end

hfig(1)=figure;

hfig(2)=figure;

figure(hfig(1));

   plot(tt,y1,'k-.',tt,xae,'k*',tt,y2,'k-',tt,xpe,'k*');

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

xlabel('t - часы');ylabel('концентрации А ("-"), P(расч. - "***"), P(аппр. - "-.-")- мольн. дол.');

grid on;  

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

nn=length(t);

disp ('5.Изменение концентрации компонентов во времени в реакторе при найденных значениях констант скоростей реакций');

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

for i=1:nn

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

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

y11= yy1(i);

y22=yy2(i);

tt=t(i);

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

end

figure(hfig(2));

plot(t,yy1,'k-.',t,yy2,'k-');

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

xlabel('t - часы');ylabel('концентрации А ("-"), P(расч. - "***"), P(аппр. - "-.-")- мольн. дол.');

grid on;

end