% % OPTIMALIZACE TEPLOTY U SILNE EXOTERMICKE REAKCE % % Zadani a predpoklady: % - vsadkovy reaktor [doba reakce znama] % - izotermicky prubeh reakce [hledame optimalni teplotu] % - konstantni objem smesi (reaktanty i produkty rozpustene) % - vratna reakce aA <---> bB tenda = 10*60; % doba reakce A tendb = 60*60; % doba reakce B cA0 = 1; % mol/dm3 a = 1; b = 1; R = 8.314; % kinetika reakce A = 2e3; Ea = 40e3; % aktivacni energie (citlivost na teplote), J/mol % rovnovaha reakce dHr = -70e3; % reakcni entalpie (exotermicka reakce), J/mol K0 = 200; % rovnovazna konstanta pri T0 T0 = 298.15; % K Z = exp(log(K0) + dHr/R/T0); % FUNKCE PRO VYPOCET ZAVISLOSTI "K" a "kp" na teplote FK = @(T) Z*exp(-dHr/R/T); % rovnovazna konstanta Fkp = @(T) A*exp(-Ea/R/T); % rychlostni konstanta % FUNKCE PRO VYPOCET ZAVISLOSTI RYCHLOSTI REAKCE na konverzi a teplote FR = @(X,T) Fkp(T)*(-cA0^a*(1-X)^a + 1/FK(T)*cA0^b*(b/a*X)^b); % prohledavany rozsah teplot T = linspace(T0,450,70); Xenda = []; % dosazena konverze za dobu A Xendb = []; % dosazena konverze za dobu B Xeq = []; % konverze pri rovnovaze for Ti=T fun = @(t,y) -a*FR(y,Ti)/cA0; % derivace konverze (z molarni bilance) [tt,yy] = ode15s(fun,[0 tenda], 0); Xenda = [Xenda, yy(end)]; [tt,yy] = ode15s(fun,[0 tendb], 0); Xendb = [Xendb, yy(end)]; Xeqnow = fzero(@(X) FR(X,Ti), [0.0000000000001,0.999999999999]); Xeq = [Xeq, Xeqnow]; end plot(T,Xenda,'o',T,Xendb,'+',T,Xeq,'-') xlabel('temperature / K') ylabel('conversion / [-]') legend('time A','time B','equilibrium') title('Batch reactor') |