Následující kód v Matlabu ukazuje příklad výpočtu reakční entalpie, reakční Gibbsovy energie, rovnovážné konstanty a rovnovážného složení pro oxidaci oxidu siřičitého. Vycházíme ze znalosti standardních slučovacích entalpií a standardních Gibbsových energií. Závislost těchto veličin na teplotě vypočítáme ze znalosti tepelných kapacit pomocí tzv. Kirchofovy rovnice a van’t Hoffovy rovnice. Rovnovážné složení potom počítáme pro vsádkový izotermický reaktor, který na začátku obsahuje směs oxidu siřičitého a kyslíku ve stechiometrickém složení při tlaku 1 bar a uvažované teplotě. Je důležité poznamenat, že po dosažení rovnováhy bude (vzhledem ke konstantnímu objemu reaktoru) tlak nižší než 1 bar.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % PRIKLAD - OXIDACE SO2 NA SO3 % % Vypocet dHr, dGr, K a rovnovazneho slozeni % % SO2 + 1/2 O2 --> SO3, SLOZKY=[SO2(g), O2(g), SO3(g)] % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% nu = [-1 -1/2 +1]; % stechiometricke koeficienty %%%%%%%%%%%%%%%%%%%%%%% % TERMODYNAMICKA DATA % %%%%%%%%%%%%%%%%%%%%%%% % koeficienty pro tepelne kapacity % Cp = a + bT + c/T^2 + dT^2 [J/mol.K] a = [43.43 29.96 57.145]; b = [10.63 4.18 27.347]*1e-3; c = [-5.94 -1.67 -12.912]*1e5; d = [0 0 -7.728]*1e-6; Tr = 298.15; % referencni teplota / K R = 8.314; % slucovaci enthalpie pri ref. teplote Hf = [-296.842 0 -395.765]*1e3; % J/mol % td. funkce "-(G-H298)/T" pri ref. teplote Tr TDF298 = [248.212, 205.147 256.769]; % T = 298.15 K % td. funkce "-(G_H298)/T" pri zvysene teplote (pro kontrolu) % TDF2= [271.339, 220.875 287.768]; % T = 1000 K TDF2= [252.979, 208.524 262.992]; % T = 500 K % T2 = 1000; T2 = 500; %%%%%%%%%%%%%%%%%%%% % REAKCNI ENTALPIE % %%%%%%%%%%%%%%%%%%%% dHr298 = sum(nu.*Hf) % J/mol ar = sum(nu.*a); br = sum(nu.*b); cr = sum(nu.*c); dr = sum(nu.*d); % zavislost reakcni entalpie na teplote FdHr = @(T) dHr298 + ar*(T-Tr) ... +br/2*(T.^2-Tr^2) ... -cr*(1./T-1/Tr) ... +dr/3*(T.^3-Tr^3); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % REAKCNI GIBBSOVA ENERGIE / ROVNOVAZNA KONSTANTA % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% dGr298 = sum(nu.*(-TDF298*Tr))+dHr298 % J/mol lnK298 = -dGr298/R/Tr % zavislost rovnovazne konstanty na teplote HH = dHr298 - ar*Tr - br/2*Tr^2 + cr/Tr - dr/3*Tr^3; FlnK21 = @(T1,T2) 1/R*(-HH*(1./T2-1./T1) ... +ar*log(T2./T1) ... +br/2*(T2-T1) ... +cr/2*(1./T2.^2-1./T1.^2) ... +dr/6*(T2.^2-T1.^2)); FlnK = @(T) lnK298 + FlnK21(Tr,T); % van't Hoffova rovnice FdGr = @(T) -FlnK(T) * R .* T; % gibsova slucovaci enthalpie pri zvysene teplote (pro kontrolu) dGr2 = sum(nu.*(-TDF2*T2))+dHr298; lnK2 = -dGr2/R/T2; %%%%%%%%%%%%%%%%%%%%%% % VYSLEDKY A GRAFY I % %%%%%%%%%%%%%%%%%%%%%% T = linspace(273.15,1000,40); % sledovany rozsah yH = FdHr(T)/1000; % kJ/mol yG = FdGr(T)/1000; subplot(2,2,1) plot(T,yH) xlabel('temp / K') ylabel('dHr / kJ/mol') title('SO2 + 1/2 O2 --> SO3') subplot(2,2,2) plot(T,yG) xlabel('temp / K') ylabel('dGr / kJ/mol') title('SO2 + 1/2 O2 --> SO3') T2 = T2 % hodnoty dGr a lnK pri teplote T2 dGr2 = [ FdGr(T2) dGr2] % - srovnani tabulky a van't Hoff integral dlnK2 = [ FlnK(T2) lnK2] %%%%%%%%%%%%%%%%%%%%%% % ROVNOVAZNE SLOZENI % %%%%%%%%%%%%%%%%%%%%%% n0 = [1, 0.5, 0]; % pocatecni latkova mnozstvi (odpovida stech.) % - zavislost kvocientu Q na rozsahu reakce "xi" % - odvozeno pro [V,T]=konst. % - pocatecni tlak (pri xi=0) je "p0" FQ = @(xi,p0) (n0(3)+xi) / ... ((n0(1)-xi)*(n0(2)-0.5*xi)^0.5) ... /(p0/sum(n0))^0.5; % zavislost rovnovazneho slozeni na teplote ... p0 = 1; % ... pro tlak 1 bar xi = []; for Ti=T; xii = fzero(@(xi) log(FQ(xi,p0))-FlnK(Ti), [0.000000000001,0.999999999999]) xi = [xi, xii]; % rozsah reakce v rovnovaze end xSO3 = (n0(3)+xi)./(sum(n0)-0.5*xi); % slozeni smesi pri rovnovaze xSO2 = (n0(1)-xi)./(sum(n0)-0.5*xi); xO2 = (n0(2)-0.5*xi)./(sum(n0)-0.5*xi); p = p0/sum(n0)*(sum(n0)-0.5*xi); % tlak pri rovnovaze subplot(2,2,3) plot(T,xi,T,xSO3,T,xSO2,T,xO2) xlabel('temp / K') ylabel('[-]') legend('xi','xS03','xSO2','xO2','Location','northwest') title('izothermic batch reactor [T,V]=const.') subplot(2,2,4) plot(T,p) xlabel('temp / K') ylabel('p / bar') title('izothermic batch reactor [T,V]=const.') |
