function yEqfuel = functionEqfim0(Pcalculo,Temperatura,PHIcalculo,percentalcool); warning off % Programa Combustão % Departamento de pós-graduação em Engenharia Mecânica % Autor: Juan Canellas Bosch Neto Versão em MATLAB (Aluno doutorado) % Revisor Prof.Mautone (Orientador) % Referencias: % FERGUSON, Internal Combustion Engines % RUGIERO, Cálculo Numérico % clear; % clc; % disp('Programa de equilíbrio dos produtos de combustão') % disp('') % disp('Combustível considerado gasolina e alcool') % disp('Entre com a percentagem de etanol (sugestão 28.739 % para resultar em 50 % de fração molar dos componentes alcool e gasolina') % disp('ou para comparações com os resultados do GASEQ entre com a percentagem de etanol igual a zero para resultar em 100 % (gasolina pura) e compare com o GASEQ') % percentalcool=input('Percentual de etanol na mistura gasolina/alcool: ') % disp('Entre com a pressão em atm : '); % Pcalculo=input('P: ') % disp('Entre com a temperatura (700 a 3000) em K : '); % Temperatura=input('T: ') % disp('Entre com a relação combustível/ar Faixa: 0.6 a 2 , podendo ultrapassar o PHI=2 (molar ou volumétrica: '); % PHIcalculo=input('PHI: ') % formula minima artigo Tadeu equifuel=[1 2.2 0.13;1 2.5 0.26;1 2.7 0.34;1 2.9 0.46;1 3 0.53]; if percentalcool==0 ALFA=equifuel(1,1); BETA=equifuel(1,2); GAMA=equifuel(1,3); DELTA=0; end if percentalcool==30 ALFA=equifuel(2,1); BETA=equifuel(2,2); GAMA=equifuel(2,3); DELTA=0; end if percentalcool==50 ALFA=equifuel(3,1); BETA=equifuel(3,2); GAMA=equifuel(3,3); DELTA=0; end if percentalcool==80 ALFA=equifuel(4,1); BETA=equifuel(4,2); GAMA=equifuel(4,3); DELTA=0; end if percentalcool==100 ALFA=equifuel(5,1); BETA=equifuel(5,2); GAMA=equifuel(5,3); DELTA=0; end e=0.210/(ALFA+BETA/4-GAMA/2); format short e % Pressão, temperatura e PHI ponto de partida P=1; T=1800; % Passos PHI=0.6; passot=0.5; passop=0.5; if Temperatura<800 passoPHI=0.001; else passoPHI=0.001; end passoPHI2=0.0001; %%%%% % Resolucao de sistemas nao-lineares Metodo de Newton % 23 equacoes e 23incognitas %y1y2y3y4y5y6y7y8y9y10 N y12 y13 y14 y15 y16y17y18y19y20y21y22y23y24 % CO2;H2O;N2;O2;C0;H2;H;O;0H;N0;N;CH4;NH3;HCN;CH2O;;N2O;NO2;HO2;N;NH2; % HCO;CN,NH % Estimativas iniciais para o Metodo de Newton y3=0.75425; y2=0.08649; y1=0.07707; y5=5.56E-05; y4=0.0789; y9=5.13E-04; y7=1.45E-06; y8=3.38E-05; y6=1.65E-05; y10=2.66E-03; y21=4.59E-14; y15=6.41E-16; y12=5.99E-25; y18=6.01E-07; y17=3.86E-06; y13=1.75E-12; y20=4.09E-13; y19=3.16E-11; y14=6.57E-17; y22=1.15E-19; y16=1.34E-07; K1=[0.432168 -.112464e5 0.267269e1 -.745744e-4 0.242484e-08]; K2=[0.310805 -.129540e5 0.321779e1 -0.738336e-4 0.344645e-8]; K3=[-0.141784 -0.213308e04 0.853461 0.355015e-4 -0.310227e-8]; K4=[0.150879e-01 -0.470959e04 0.646096 0.272805e-5 -0.154444e-8]; K5=[-0.752364 0.124210e05 -.260286e1 0.259656e-3 -0.162687e-07]; K6=[-0.415302e-02 0.148627e05 -.475746e1 0.124699e-3 -.900227e-8]; K7=[2.4456145 12197.386 -11.040409 -0.0031378133 5.5687468e-7]; K8=[-3.7548739 1051.1775 -6.2837125 0.0022003969 -2.1301305e-7]; K9=[-7.0150422 -16498.583 9.8642759 0.0069034474 -1.0590558e-6]; K10=[-1.1813411 -27320.566 0.083588227 0.0013032983 -2.1193592e-007]; % K11=[0.72911664 18616.818 -0.68104229 -0.00044548766 4.119092e-008]; CH3 K12=[0.23083205 5900.9553 0.74020448 -0.00027580823 2.5583576e-008]; K13=[0.94696264 2057.1509 -0.80652256 -0.00046615688 3.7696605e-008]; K14=[6.1493112 14860.941 0.11347102 -0.0037818991 3.6258144e-007]; K15=[-0.37240054 16214.984 -0.52570649 0.0001565831 -1.1439649e-008]; K16=[-0.078081888 4686.4491 0.25416542 -7.0632838e-5 8.5401494e-9]; K17=[-0.26174852 19127.31 -0.29007952 -5.0007562e-005 7.2647514e-009]; K18=[56.603446 41515.128 -15.440239 -0.025765384 1.9183826e-006]; K19=[-0.22005882 15382.953 -0.92406951 0.00020533884 -2.1757966e-008]; N=1; if Temperatura>1800 for T=1800:passot:Temperatura % Inicio do programa KP1=10^((K1(1)*log(T/1000)+(K1(2)/T)+K1(3)+K1(4)*T+K1(5)*T^2)); KP2=10^((K2(1)*log(T/1000)+(K2(2)/T)+K2(3)+K2(4)*T+K2(5)*T^2)); KP3=10^((K3(1)*log(T/1000)+(K3(2)/T)+K3(3)+K3(4)*T+K3(5)*T^2)); KP4=10^((K4(1)*log(T/1000)+(K4(2)/T)+K4(3)+K4(4)*T+K4(5)*T^2)); KP5=10^((K5(1)*log(T/1000)+(K5(2)/T)+K5(3)+K5(4)*T+K5(5)*T^2)); KP6=10^((K6(1)*log(T/1000)+(K6(2)/T)+K6(3)+K6(4)*T+K6(5)*T^2)); KP7=10^((K7(1)*log(T/1000)+(K7(2)/T)+K7(3)+K7(4)*T+K7(5)*T^2)); KP8=10^((K8(1)*log(T/1000)+(K8(2)/T)+K8(3)+K8(4)*T+K8(5)*T^2)); KP9=10^((K9(1)*log(T/1000)+(K9(2)/T)+K9(3)+K9(4)*T+K9(5)*T^2)); KP10=10^((K10(1)*log(T/1000)+(K10(2)/T)+K10(3)+K10(4)*T+K10(5)*T^2)); %KP11=10^((K11(1)*log(T/1000)+(K11(2)/T)+K11(3)+K11(4)*T+K11(5)*T^2)); %Ch3 KP12=10^((K12(1)*log(T/1000)+(K12(2)/T)+K12(3)+K12(4)*T+K12(5)*T^2)); KP13=10^((K13(1)*log(T/1000)+(K13(2)/T)+K13(3)+K13(4)*T+K13(5)*T^2)); KP14=10^((K14(1)*log(T/1000)+(K14(2)/T)+K14(3)+K14(4)*T+K14(5)*T^2)); KP15=10^((K15(1)*log(T/1000)+(K15(2)/T)+K15(3)+K15(4)*T+K15(5)*T^2)); KP16=10^((K16(1)*log(T/1000)+(K16(2)/T)+K16(3)+K16(4)*T+K16(5)*T^2)); KP17=10^((K17(1)*log(T/1000)+(K17(2)/T)+K17(3)+K17(4)*T+K17(5)*T^2)); KP18=10^((K18(1)*log(T/1000)+(K18(2)/T)+K18(3)+K18(4)*T+K18(5)*T^2)); KP19=10^((K19(1)*log(T/1000)+(K19(2)/T)+K19(3)+K19(4)*T+K19(5)*T^2)); c1=KP1/sqrt(P); c2=KP2/sqrt(P); c3=KP3; c4=KP4; c5=KP5*sqrt(P); c6=KP6*sqrt(P); c7=KP7*(P^2); c8=KP8*P; c9=KP9/(P^2); c10=KP10; %c11=1/KP11; % CH3 c12=1/KP12; c13=KP13; c14=1/KP14; c15=1/KP15; c16=KP16; c17=1/KP17; c18=1/KP18; c19=1/KP19; d1=1/(c1^2); d2=1/(c2^2); d3=1/(c3^2); d4=1/(c4^2); d5=1/(c5^2); d6=1/(c6^2); d8=1/(c8^2); for i=1:1 F5Y7=2*d1*y7; F6Y8=2*d2*y8; F7Y4=(-(y9^2)*d3)/(y4^2); F7Y9=2*y9*d3/y4; F8Y4=(-d4*(y10)^2)/(y4^2); F8Y10=(2*y10*d4)/y4; F9Y2=2*y2*d5/(y6^2); F9Y6=-2*(y2^2)*d5/(y6^3); F10Y1=(2*y1*d6)/(y5^2); F10Y5=-2*(y1^2)*d6/(y5^3); Fn1=(y1+y5+y12+y14+y15+y21+y22); Fn2=(2*y2+2*y6+y7+y9+4*y12+3*y13+y14+2*y15+y18+2*y20+y21); Fn3=(2*y1+y2+2*y4+y5+y8+y9+y10+y15+y16+2*y17+2*y18+y21); Fn4=(2*y3+y10+y13+y14+2*y16+y17+y19+y20+y22); F12Y2=-c7*y5*(y6^3)/(y2^2); F12Y5=c7*(y6^3)/y2; F12Y6=c7*y5*3*(y6^2)/y2; F13Y13=2*d8*y13/(y6^3); F13Y6=-(d8*3*(y13^2)/(y6^4)); F14Y6=-3*c9*y13*y12/(y6^4); F14Y12=c9*y13/(y6^3); F14Y13=c9*y12/(y6^3); F15Y1=c10*y2/y4; F15Y2=c10*y1/y4; F15Y4=-c10*y1*y2/(y4^2); F16Y3=c12*y18/y9; F16Y9=-c12*y3*y18/(y9^2); F16Y18=c12*y3/y9; F17Y9=-c13*y10*y18/(y9^2); F17Y10=c13*y18/y9; F17Y18=c13*y10/y9; F18Y4=c14*y9/y8; F18Y8=-c14*y9*y4/(y8^2); F18Y9=c14*y4/y8; F19Y3=c15*y8/y10; F19Y8=c15*y3/y10; F19Y10=-c15*y3*y8/(y10^2); F20Y5=-c16*y14*y9/(y5^2); F20Y9=c16*y14/y5; F20Y14=c16*y9/y5; F21Y5=c17*y6/y7; F21Y6=c17*y5/y7; F21Y7=-c17*y5*y6/(y7^2); F22Y2=-c18*y14*y9/(y2^2); F22Y9=c18*y14/y2; F22Y14=c18*y9/y2; % Matriz Jacobiana de derivadas parciais a=[N 0 0 0 N 0 0 0 0 0 Fn1 N 0 N N 0 0 0 0 0 N N; 0 2*N 0 0 0 2*N N 0 N 0 Fn2 4*N 3*N N 2*N 0 0 N 0 2*N N 0; 2*N N 0 2*N N 0 0 N N N Fn3 0 0 0 0 N 2*N 2*N 0 0 N 1; 0 0 2*N 0 0 0 0 0 0 N Fn4 0 N N 0 2*N N 0 N N 0 N; 0 0 0 0 0 -1 F5Y7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 -1 0 0 0 F6Y8 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 F7Y4 0 -1 0 0 F7Y9 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 -1 F8Y4 0 0 0 0 0 F8Y10 0 0 0 0 0 0 0 0 0 0 0 0; 0 F9Y2 0 -1 0 F9Y6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; F10Y1 0 0 -1 F10Y5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1; 0 F12Y2 0 0 F12Y5 F12Y6 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0; 0 0 -1 0 0 F13Y6 0 0 0 0 0 0 F13Y13 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 F14Y6 0 0 0 0 0 F14Y12 F14Y13 -1 0 0 0 0 0 0 0 0; F15Y1 F15Y2 0 F15Y4 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0; 0 0 F16Y3 0 0 0 0 0 F16Y9 0 0 0 0 0 0 -1 0 F16Y18 0 0 0 0; 0 0 0 0 0 0 0 0 F17Y9 F17Y10 0 0 0 0 0 0 -1 F17Y18 0 0 0 0; 0 0 0 F18Y4 0 0 0 F18Y8 F18Y9 0 0 0 0 0 0 0 0 -1 0 0 0 0; 0 0 F19Y3 0 0 0 0 F19Y8 0 F19Y10 0 0 0 0 0 0 0 0 -1 0 0 0; 0 0 0 0 F20Y5 0 0 0 F20Y9 0 0 0 0 F20Y14 0 0 0 0 0 -1 0 0; 0 0 0 0 F21Y5 F21Y6 F21Y7 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0; 0 F22Y2 0 0 0 0 0 0 F22Y9 0 0 0 0 F22Y14 0 0 0 0 0 0 0 -1]; % Eq. 3.56 b1=((y1+y5+y12+y14+y15+y21+y22)*N)-(e*PHI*ALFA); % Eq. 3.57 b2=((2*y2+2*y6+y7+y9+4*y12+3*y13+y14+2*y15+y18+2*y20+y21)*N)-(e*PHI*BETA); % Eq. 3.58 b3=((2*y1+y2+2*y4+y5+y8+y9+y10+y15+y16+2*y17+2*y18+y21)*N)-((e*PHI*GAMA)+0.42); % Eq. 3.59 b4=((2*y3+y10+y13+y14+2*y16+y17+y19+y20+y22)*N)-((e*PHI*DELTA)+1.58); % Eq. 3.68 b5=(d1*y7^2)-y6; % Eq. 3.69 b6=(d2*y8^2)-y4; % Eq. 3.70 b7=((d3*y9^2)/y4)-y6; % Eq. 3.71 b8=((d4*(y10^2))/y4)-y3; % Eq. 3.72 b9=((d5*y2^2)/(y6^2))-y4; % Eq. 3.73 b10=((d6*y1^2)/(y5^2))-y4; % Eq.3.60 b11=y1+y2+y3+y4+y5+y6+y7+y8+y9+y10+y12+y13+y14+y15+y16+y17+y18+y19+y20+y21+y22-1; % Equacao do metano CO + 3H2 === CH4 + H20 b12=(c7*y5*(y6^3)/y2)-y12; % Equação da amônia 1/2 N2 + 3/2 H2 === NH3 b13=(((y13^2)*d8)/y6^3)-y3; % Equação do ácido cianídrico NH3 + CH4 === HCN+ 3H2 b14=(c9*y13*y12/(y6^3))-y14; % Equação do Formol C02 + H20 ===CH20+ O2 b15=(c10*y1*y2/y4)-y15; % Equação do CH3 CH3 + NO == HCN + H2O %b16=(c11*y14*y2/y10)-y16; % Equação do N2O N2O + OH == N2 + HO2 b16=(c12*y3*y18/y9)-y16; % Equação do NO2 HO2 + NO === NO2 + OH b17=(c13*y10*y18/y9)-y17; % Equação do HO2 0 + H02 === OH + 02 (substituta) b18=(c14*y9*y4/y8)-y18; % Equação do N N + NO ==== N2 + 0 b19=(c15*y3*y8/y10)-y19; % Equação do NH2 HCN + 0H == NH2 + CO b20=(c16*y14*y9/y5)-y20; % Equacao do HCO H + HCO == H2 + CO b21=(c17*y5*y6/y7)-y21; % Equacao do CN CN + H20 == HCN + OH b22=(c18*y14*y9/y2)-y22; % Equacao NH NH + 0 === NO + H %b24=(c19*y10*y7/y8)-y24; %\Montando o vetor B com as equacoes b=[-b1;-b2;-b3;-b4;-b5;-b6;-b7;-b8;-b9;-b10;-b11;-b12;-b13;-b14;-b15;-b16;-b17;-b18;-b19;-b20;-b21;-b22]; %Resolucao de sistemas lineares pelo MATLAB ao inves de x=A\B, tem-se x=inv(a)*b y=a\b; y1=y(1)+y1; y2=y(2)+y2; y3=y(3)+y3; y4=y(4)+y4; y5=y(5)+y5; y6=y(6)+y6; y7=y(7)+y7; y8=y(8)+y8; y9=y(9)+y9; y10=y(10)+y10; N=y(11)+N; y12=y(12)+y12; y13=y(13)+y13; y14=y(14)+y14; y15=y(15)+y15; y16=y(16)+y16; y17=y(17)+y17; y18=y(18)+y18; y19=y(19)+y19; y20=y(20)+y20; y21=y(21)+y21; y22=y(22)+y22; %y23=y(23)+y23; %y24=y(24)+y24; end ytotal=y1+y2+y3+y4+y5+y6+y7+y8+y9+y10+y12+y13+y14+y15+y16+y17+y18+y19+y20+y21+y22; end y1=abs(y1); y2=abs(y2); y3=abs(y3); y4=abs(y4); y5=abs(y5); y6=abs(y6); y7=abs(y7); y8=abs(y8); y9=abs(y9); y10=abs(y10); y12=abs(y12); y13=abs(y13); y14=abs(y14); y15=abs(y15); y16=abs(y16); y17=abs(y17); y18=abs(y18); y19=abs(y19); y20=abs(y20); y21=abs(y21); y22=abs(y22); else for T=1800:-passot:Temperatura % Inicio do programa KP1=10^((K1(1)*log(T/1000)+(K1(2)/T)+K1(3)+K1(4)*T+K1(5)*T^2)); KP2=10^((K2(1)*log(T/1000)+(K2(2)/T)+K2(3)+K2(4)*T+K2(5)*T^2)); KP3=10^((K3(1)*log(T/1000)+(K3(2)/T)+K3(3)+K3(4)*T+K3(5)*T^2)); KP4=10^((K4(1)*log(T/1000)+(K4(2)/T)+K4(3)+K4(4)*T+K4(5)*T^2)); KP5=10^((K5(1)*log(T/1000)+(K5(2)/T)+K5(3)+K5(4)*T+K5(5)*T^2)); KP6=10^((K6(1)*log(T/1000)+(K6(2)/T)+K6(3)+K6(4)*T+K6(5)*T^2)); KP7=10^((K7(1)*log(T/1000)+(K7(2)/T)+K7(3)+K7(4)*T+K7(5)*T^2)); KP8=10^((K8(1)*log(T/1000)+(K8(2)/T)+K8(3)+K8(4)*T+K8(5)*T^2)); KP9=10^((K9(1)*log(T/1000)+(K9(2)/T)+K9(3)+K9(4)*T+K9(5)*T^2)); KP10=10^((K10(1)*log(T/1000)+(K10(2)/T)+K10(3)+K10(4)*T+K10(5)*T^2)); %KP11=10^((K11(1)*log(T/1000)+(K11(2)/T)+K11(3)+K11(4)*T+K11(5)*T^2)); %Ch3 KP12=10^((K12(1)*log(T/1000)+(K12(2)/T)+K12(3)+K12(4)*T+K12(5)*T^2)); KP13=10^((K13(1)*log(T/1000)+(K13(2)/T)+K13(3)+K13(4)*T+K13(5)*T^2)); KP14=10^((K14(1)*log(T/1000)+(K14(2)/T)+K14(3)+K14(4)*T+K14(5)*T^2)); KP15=10^((K15(1)*log(T/1000)+(K15(2)/T)+K15(3)+K15(4)*T+K15(5)*T^2)); KP16=10^((K16(1)*log(T/1000)+(K16(2)/T)+K16(3)+K16(4)*T+K16(5)*T^2)); KP17=10^((K17(1)*log(T/1000)+(K17(2)/T)+K17(3)+K17(4)*T+K17(5)*T^2)); KP18=10^((K18(1)*log(T/1000)+(K18(2)/T)+K18(3)+K18(4)*T+K18(5)*T^2)); KP19=10^((K19(1)*log(T/1000)+(K19(2)/T)+K19(3)+K19(4)*T+K19(5)*T^2)); c1=KP1/sqrt(P); c2=KP2/sqrt(P); c3=KP3; c4=KP4; c5=KP5*sqrt(P); c6=KP6*sqrt(P); c7=KP7*(P^2); c8=KP8*P; c9=KP9/(P^2); c10=KP10; %c11=1/KP11; % CH3 c12=1/KP12; c13=KP13; c14=1/KP14; c15=1/KP15; c16=KP16; c17=1/KP17; c18=1/KP18; c19=1/KP19; d1=1/(c1^2); d2=1/(c2^2); d3=1/(c3^2); d4=1/(c4^2); d5=1/(c5^2); d6=1/(c6^2); d8=1/(c8^2); for i=1:50 F5Y7=2*d1*y7; F6Y8=2*d2*y8; F7Y4=(-(y9^2)*d3)/(y4^2); F7Y9=2*y9*d3/y4; F8Y4=(-d4*(y10)^2)/(y4^2); F8Y10=(2*y10*d4)/y4; F9Y2=2*y2*d5/(y6^2); F9Y6=-2*(y2^2)*d5/(y6^3); F10Y1=(2*y1*d6)/(y5^2); F10Y5=-2*(y1^2)*d6/(y5^3); Fn1=(y1+y5+y12+y14+y15+y21+y22); Fn2=(2*y2+2*y6+y7+y9+4*y12+3*y13+y14+2*y15+y18+2*y20+y21); Fn3=(2*y1+y2+2*y4+y5+y8+y9+y10+y15+y16+2*y17+2*y18+y21); Fn4=(2*y3+y10+y13+y14+2*y16+y17+y19+y20+y22); F12Y2=-c7*y5*(y6^3)/(y2^2); F12Y5=c7*(y6^3)/y2; F12Y6=c7*y5*3*(y6^2)/y2; F13Y13=2*d8*y13/(y6^3); F13Y6=-(d8*3*(y13^2)/(y6^4)); F14Y6=-3*c9*y13*y12/(y6^4); F14Y12=c9*y13/(y6^3); F14Y13=c9*y12/(y6^3); F15Y1=c10*y2/y4; F15Y2=c10*y1/y4; F15Y4=-c10*y1*y2/(y4^2); F16Y3=c12*y18/y9; F16Y9=-c12*y3*y18/(y9^2); F16Y18=c12*y3/y9; F17Y9=-c13*y10*y18/(y9^2); F17Y10=c13*y18/y9; F17Y18=c13*y10/y9; F18Y4=c14*y9/y8; F18Y8=-c14*y9*y4/(y8^2); F18Y9=c14*y4/y8; F19Y3=c15*y8/y10; F19Y8=c15*y3/y10; F19Y10=-c15*y3*y8/(y10^2); F20Y5=-c16*y14*y9/(y5^2); F20Y9=c16*y14/y5; F20Y14=c16*y9/y5; F21Y5=c17*y6/y7; F21Y6=c17*y5/y7; F21Y7=-c17*y5*y6/(y7^2); F22Y2=-c18*y14*y9/(y2^2); F22Y9=c18*y14/y2; F22Y14=c18*y9/y2; % Matriz Jacobiana de derivadas parciais a=[N 0 0 0 N 0 0 0 0 0 Fn1 N 0 N N 0 0 0 0 0 N N; 0 2*N 0 0 0 2*N N 0 N 0 Fn2 4*N 3*N N 2*N 0 0 N 0 2*N N 0; 2*N N 0 2*N N 0 0 N N N Fn3 0 0 0 0 N 2*N 2*N 0 0 N 1; 0 0 2*N 0 0 0 0 0 0 N Fn4 0 N N 0 2*N N 0 N N 0 N; 0 0 0 0 0 -1 F5Y7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 -1 0 0 0 F6Y8 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 F7Y4 0 -1 0 0 F7Y9 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 -1 F8Y4 0 0 0 0 0 F8Y10 0 0 0 0 0 0 0 0 0 0 0 0; 0 F9Y2 0 -1 0 F9Y6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; F10Y1 0 0 -1 F10Y5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1; 0 F12Y2 0 0 F12Y5 F12Y6 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0; 0 0 -1 0 0 F13Y6 0 0 0 0 0 0 F13Y13 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 F14Y6 0 0 0 0 0 F14Y12 F14Y13 -1 0 0 0 0 0 0 0 0; F15Y1 F15Y2 0 F15Y4 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0; 0 0 F16Y3 0 0 0 0 0 F16Y9 0 0 0 0 0 0 -1 0 F16Y18 0 0 0 0; 0 0 0 0 0 0 0 0 F17Y9 F17Y10 0 0 0 0 0 0 -1 F17Y18 0 0 0 0; 0 0 0 F18Y4 0 0 0 F18Y8 F18Y9 0 0 0 0 0 0 0 0 -1 0 0 0 0; 0 0 F19Y3 0 0 0 0 F19Y8 0 F19Y10 0 0 0 0 0 0 0 0 -1 0 0 0; 0 0 0 0 F20Y5 0 0 0 F20Y9 0 0 0 0 F20Y14 0 0 0 0 0 -1 0 0; 0 0 0 0 F21Y5 F21Y6 F21Y7 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0; 0 F22Y2 0 0 0 0 0 0 F22Y9 0 0 0 0 F22Y14 0 0 0 0 0 0 0 -1]; % Termos independentes % Sistema de eq. não lineares 10x10 % Eqs. 3.56 3.57 3.58 3.59 3.68 3.69 3.70 3.71 3.72 3.73 pag 122 e 123 % Eq. 3.56 b1=((y1+y5+y12+y14+y15+y21+y22)*N)-(e*PHI*ALFA); % Eq. 3.57 b2=((2*y2+2*y6+y7+y9+4*y12+3*y13+y14+2*y15+y18+2*y20+y21)*N)-(e*PHI*BETA); % Eq. 3.58 b3=((2*y1+y2+2*y4+y5+y8+y9+y10+y15+y16+2*y17+2*y18+y21)*N)-((e*PHI*GAMA)+0.42); % Eq. 3.59 b4=((2*y3+y10+y13+y14+2*y16+y17+y19+y20+y22)*N)-((e*PHI*DELTA)+1.58); % Eq. 3.68 b5=(d1*y7^2)-y6; % Eq. 3.69 b6=(d2*y8^2)-y4; % Eq. 3.70 b7=((d3*y9^2)/y4)-y6; % Eq. 3.71 b8=((d4*(y10^2))/y4)-y3; % Eq. 3.72 b9=((d5*y2^2)/(y6^2))-y4; % Eq. 3.73 b10=((d6*y1^2)/(y5^2))-y4; % Eq.3.60 b11=y1+y2+y3+y4+y5+y6+y7+y8+y9+y10+y12+y13+y14+y15+y16+y17+y18+y19+y20+y21+y22-1; % Equacao do metano CO + 3H2 === CH4 + H20 b12=(c7*y5*(y6^3)/y2)-y12; % Equação da amônia 1/2 N2 + 3/2 H2 === NH3 b13=(((y13^2)*d8)/y6^3)-y3; % Equação do ácido cianídrico NH3 + CH4 === HCN+ 3H2 b14=(c9*y13*y12/(y6^3))-y14; % Equação do Formol C02 + H20 ===CH20+ O2 b15=(c10*y1*y2/y4)-y15; % Equação do CH3 CH3 + NO == HCN + H2O %b16=(c11*y14*y2/y10)-y16; % Equação do N2O N2O + OH == N2 + HO2 b16=(c12*y3*y18/y9)-y16; % Equação do NO2 HO2 + NO === NO2 + OH b17=(c13*y10*y18/y9)-y17; % Equação do HO2 0 + H02 === OH + 02 (substituta) b18=(c14*y9*y4/y8)-y18; % Equação do N N + NO ==== N2 + 0 b19=(c15*y3*y8/y10)-y19; % Equação do NH2 HCN + 0H == NH2 + CO b20=(c16*y14*y9/y5)-y20; % Equacao do HCO H + HCO == H2 + CO b21=(c17*y5*y6/y7)-y21; % Equacao do CN CN + H20 == HCN + OH b22=(c18*y14*y9/y2)-y22; % Equacao NH NH + 0 === NO + H %b24=(c19*y10*y7/y8)-y24; %\Montando o vetor B com as equacoes b=[-b1;-b2;-b3;-b4;-b5;-b6;-b7;-b8;-b9;-b10;-b11;-b12;-b13;-b14;-b15;-b16;-b17;-b18;-b19;-b20;-b21;-b22]; %Resolucao de sistemas lineares pelo MATLAB ao inves de x=A\B, tem-se x=inv(a)*b y=a\b; y1=y(1)+y1; y2=y(2)+y2; y3=y(3)+y3; y4=y(4)+y4; y5=y(5)+y5; y6=y(6)+y6; y7=y(7)+y7; y8=y(8)+y8; y9=y(9)+y9; y10=y(10)+y10; N=y(11)+N; y12=y(12)+y12; y13=y(13)+y13; y14=y(14)+y14; y15=y(15)+y15; y16=y(16)+y16; y17=y(17)+y17; y18=y(18)+y18; y19=y(19)+y19; y20=y(20)+y20; y21=y(21)+y21; y22=y(22)+y22; end ytotal=y1+y2+y3+y4+y5+y6+y7+y8+y9+y10+y12+y13+y14+y15+y16+y17+y18+y19+y20+y21+y22; end y1=abs(y1); y2=abs(y2); y3=abs(y3); y4=abs(y4); y5=abs(y5); y6=abs(y6); y7=abs(y7); y8=abs(y8); y9=abs(y9); y10=abs(y10); y12=abs(y12); y13=abs(y13); y14=abs(y14); y15=abs(y15); y16=abs(y16); y17=abs(y17); y18=abs(y18); y19=abs(y19); y20=abs(y20); y21=abs(y21); y22=abs(y22); end for P=1:passop:Pcalculo KP1=10^((K1(1)*log(T/1000)+(K1(2)/T)+K1(3)+K1(4)*T+K1(5)*T^2)); KP2=10^((K2(1)*log(T/1000)+(K2(2)/T)+K2(3)+K2(4)*T+K2(5)*T^2)); KP3=10^((K3(1)*log(T/1000)+(K3(2)/T)+K3(3)+K3(4)*T+K3(5)*T^2)); KP4=10^((K4(1)*log(T/1000)+(K4(2)/T)+K4(3)+K4(4)*T+K4(5)*T^2)); KP5=10^((K5(1)*log(T/1000)+(K5(2)/T)+K5(3)+K5(4)*T+K5(5)*T^2)); KP6=10^((K6(1)*log(T/1000)+(K6(2)/T)+K6(3)+K6(4)*T+K6(5)*T^2)); KP7=10^((K7(1)*log(T/1000)+(K7(2)/T)+K7(3)+K7(4)*T+K7(5)*T^2)); KP8=10^((K8(1)*log(T/1000)+(K8(2)/T)+K8(3)+K8(4)*T+K8(5)*T^2)); KP9=10^((K9(1)*log(T/1000)+(K9(2)/T)+K9(3)+K9(4)*T+K9(5)*T^2)); KP10=10^((K10(1)*log(T/1000)+(K10(2)/T)+K10(3)+K10(4)*T+K10(5)*T^2)); %KP11=10^((K11(1)*log(T/1000)+(K11(2)/T)+K11(3)+K11(4)*T+K11(5)*T^2)); %Ch3 KP12=10^((K12(1)*log(T/1000)+(K12(2)/T)+K12(3)+K12(4)*T+K12(5)*T^2)); KP13=10^((K13(1)*log(T/1000)+(K13(2)/T)+K13(3)+K13(4)*T+K13(5)*T^2)); KP14=10^((K14(1)*log(T/1000)+(K14(2)/T)+K14(3)+K14(4)*T+K14(5)*T^2)); KP15=10^((K15(1)*log(T/1000)+(K15(2)/T)+K15(3)+K15(4)*T+K15(5)*T^2)); KP16=10^((K16(1)*log(T/1000)+(K16(2)/T)+K16(3)+K16(4)*T+K16(5)*T^2)); KP17=10^((K17(1)*log(T/1000)+(K17(2)/T)+K17(3)+K17(4)*T+K17(5)*T^2)); KP18=10^((K18(1)*log(T/1000)+(K18(2)/T)+K18(3)+K18(4)*T+K18(5)*T^2)); KP19=10^((K19(1)*log(T/1000)+(K19(2)/T)+K19(3)+K19(4)*T+K19(5)*T^2)); c1=KP1/sqrt(P); c2=KP2/sqrt(P); c3=KP3; c4=KP4; c5=KP5*sqrt(P); c6=KP6*sqrt(P); c7=KP7*(P^2); c8=KP8*P; c9=KP9/(P^2); c10=KP10; %c11=1/KP11; % CH3 c12=1/KP12; c13=KP13; c14=1/KP14; c15=1/KP15; c16=KP16; c17=1/KP17; c18=1/KP18; c19=1/KP19; d1=1/(c1^2); d2=1/(c2^2); d3=1/(c3^2); d4=1/(c4^2); d5=1/(c5^2); d6=1/(c6^2); d8=1/(c8^2); % Inicio do programa for i=1:50 F5Y7=2*d1*y7; F6Y8=2*d2*y8; F7Y4=(-(y9^2)*d3)/(y4^2); F7Y9=2*y9*d3/y4; F8Y4=(-d4*(y10)^2)/(y4^2); F8Y10=(2*y10*d4)/y4; F9Y2=2*y2*d5/(y6^2); F9Y6=-2*(y2^2)*d5/(y6^3); F10Y1=(2*y1*d6)/(y5^2); F10Y5=-2*(y1^2)*d6/(y5^3); Fn1=(y1+y5+y12+y14+y15+y21+y22); Fn2=(2*y2+2*y6+y7+y9+4*y12+3*y13+y14+2*y15+y18+2*y20+y21); Fn3=(2*y1+y2+2*y4+y5+y8+y9+y10+y15+y16+2*y17+2*y18+y21); Fn4=(2*y3+y10+y13+y14+2*y16+y17+y19+y20+y22); F12Y2=-c7*y5*(y6^3)/(y2^2); F12Y5=c7*(y6^3)/y2; F12Y6=c7*y5*3*(y6^2)/y2; F13Y13=2*d8*y13/(y6^3); F13Y6=-(d8*3*(y13^2)/(y6^4)); F14Y6=-3*c9*y13*y12/(y6^4); F14Y12=c9*y13/(y6^3); F14Y13=c9*y12/(y6^3); F15Y1=c10*y2/y4; F15Y2=c10*y1/y4; F15Y4=-c10*y1*y2/(y4^2); F16Y3=c12*y18/y9; F16Y9=-c12*y3*y18/(y9^2); F16Y18=c12*y3/y9; F17Y9=-c13*y10*y18/(y9^2); F17Y10=c13*y18/y9; F17Y18=c13*y10/y9; F18Y4=c14*y9/y8; F18Y8=-c14*y9*y4/(y8^2); F18Y9=c14*y4/y8; F19Y3=c15*y8/y10; F19Y8=c15*y3/y10; F19Y10=-c15*y3*y8/(y10^2); F20Y5=-c16*y14*y9/(y5^2); F20Y9=c16*y14/y5; F20Y14=c16*y9/y5; F21Y5=c17*y6/y7; F21Y6=c17*y5/y7; F21Y7=-c17*y5*y6/(y7^2); F22Y2=-c18*y14*y9/(y2^2); F22Y9=c18*y14/y2; F22Y14=c18*y9/y2; % Matriz Jacobiana de derivadas parciais a=[N 0 0 0 N 0 0 0 0 0 Fn1 N 0 N N 0 0 0 0 0 N N; 0 2*N 0 0 0 2*N N 0 N 0 Fn2 4*N 3*N N 2*N 0 0 N 0 2*N N 0; 2*N N 0 2*N N 0 0 N N N Fn3 0 0 0 0 N 2*N 2*N 0 0 N 1; 0 0 2*N 0 0 0 0 0 0 N Fn4 0 N N 0 2*N N 0 N N 0 N; 0 0 0 0 0 -1 F5Y7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 -1 0 0 0 F6Y8 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 F7Y4 0 -1 0 0 F7Y9 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 -1 F8Y4 0 0 0 0 0 F8Y10 0 0 0 0 0 0 0 0 0 0 0 0; 0 F9Y2 0 -1 0 F9Y6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; F10Y1 0 0 -1 F10Y5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1; 0 F12Y2 0 0 F12Y5 F12Y6 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0; 0 0 -1 0 0 F13Y6 0 0 0 0 0 0 F13Y13 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 F14Y6 0 0 0 0 0 F14Y12 F14Y13 -1 0 0 0 0 0 0 0 0; F15Y1 F15Y2 0 F15Y4 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0; 0 0 F16Y3 0 0 0 0 0 F16Y9 0 0 0 0 0 0 -1 0 F16Y18 0 0 0 0; 0 0 0 0 0 0 0 0 F17Y9 F17Y10 0 0 0 0 0 0 -1 F17Y18 0 0 0 0; 0 0 0 F18Y4 0 0 0 F18Y8 F18Y9 0 0 0 0 0 0 0 0 -1 0 0 0 0; 0 0 F19Y3 0 0 0 0 F19Y8 0 F19Y10 0 0 0 0 0 0 0 0 -1 0 0 0; 0 0 0 0 F20Y5 0 0 0 F20Y9 0 0 0 0 F20Y14 0 0 0 0 0 -1 0 0; 0 0 0 0 F21Y5 F21Y6 F21Y7 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0; 0 F22Y2 0 0 0 0 0 0 F22Y9 0 0 0 0 F22Y14 0 0 0 0 0 0 0 -1]; % Termos independentes % Sistema de eq. não lineares 10x10 % Eqs. 3.56 3.57 3.58 3.59 3.68 3.69 3.70 3.71 3.72 3.73 pag 122 e 123 % Eq. 3.56 b1=((y1+y5+y12+y14+y15+y21+y22)*N)-(e*PHI*ALFA); % Eq. 3.57 b2=((2*y2+2*y6+y7+y9+4*y12+3*y13+y14+2*y15+y18+2*y20+y21)*N)-(e*PHI*BETA); % Eq. 3.58 b3=((2*y1+y2+2*y4+y5+y8+y9+y10+y15+y16+2*y17+2*y18+y21)*N)-((e*PHI*GAMA)+0.42); % Eq. 3.59 b4=((2*y3+y10+y13+y14+2*y16+y17+y19+y20+y22)*N)-((e*PHI*DELTA)+1.58); % Eq. 3.68 b5=(d1*y7^2)-y6; % Eq. 3.69 b6=(d2*y8^2)-y4; % Eq. 3.70 b7=((d3*y9^2)/y4)-y6; % Eq. 3.71 b8=((d4*(y10^2))/y4)-y3; % Eq. 3.72 b9=((d5*y2^2)/(y6^2))-y4; % Eq. 3.73 b10=((d6*y1^2)/(y5^2))-y4; % Eq.3.60 b11=y1+y2+y3+y4+y5+y6+y7+y8+y9+y10+y12+y13+y14+y15+y16+y17+y18+y19+y20+y21+y22-1; % Equacao do metano CO + 3H2 === CH4 + H20 b12=(c7*y5*(y6^3)/y2)-y12; % Equação da amônia 1/2 N2 + 3/2 H2 === NH3 b13=(((y13^2)*d8)/y6^3)-y3; % Equação do ácido cianídrico NH3 + CH4 === HCN+ 3H2 b14=(c9*y13*y12/(y6^3))-y14; % Equação do Formol C02 + H20 ===CH20+ O2 b15=(c10*y1*y2/y4)-y15; % Equação do CH3 CH3 + NO == HCN + H2O %b16=(c11*y14*y2/y10)-y16; % Equação do N2O N2O + OH == N2 + HO2 b16=(c12*y3*y18/y9)-y16; % Equação do NO2 HO2 + NO === NO2 + OH b17=(c13*y10*y18/y9)-y17; % Equação do HO2 0 + H02 === OH + 02 (substituta) b18=(c14*y9*y4/y8)-y18; % Equação do N N + NO ==== N2 + 0 b19=(c15*y3*y8/y10)-y19; % Equação do NH2 HCN + 0H == NH2 + CO b20=(c16*y14*y9/y5)-y20; % Equacao do HCO H + HCO == H2 + CO b21=(c17*y5*y6/y7)-y21; % Equacao do CN CN + H20 == HCN + OH b22=(c18*y14*y9/y2)-y22; % Equacao NH NH + 0 === NO + H %b24=(c19*y10*y7/y8)-y24; %\Montando o vetor B com as equacoes b=[-b1;-b2;-b3;-b4;-b5;-b6;-b7;-b8;-b9;-b10;-b11;-b12;-b13;-b14;-b15;-b16;-b17;-b18;-b19;-b20;-b21;-b22]; %Resolucao de sistemas lineares pelo MATLAB ao inves de x=A\B, tem-se x=inv(a)*b y=a\b; y1=y(1)+y1; y2=y(2)+y2; y3=y(3)+y3; y4=y(4)+y4; y5=y(5)+y5; y6=y(6)+y6; y7=y(7)+y7; y8=y(8)+y8; y9=y(9)+y9; y10=y(10)+y10; N=y(11)+N; y12=y(12)+y12; y13=y(13)+y13; y14=y(14)+y14; y15=y(15)+y15; y16=y(16)+y16; y17=y(17)+y17; y18=y(18)+y18; y19=y(19)+y19; y20=y(20)+y20; y21=y(21)+y21; y22=y(22)+y22; end end y1=abs(y1); y2=abs(y2); y3=abs(y3); y4=abs(y4); y5=abs(y5); y6=abs(y6); y7=abs(y7); y8=abs(y8); y9=abs(y9); y10=abs(y10); y12=abs(y12); y13=abs(y13); y14=abs(y14); y15=abs(y15); y16=abs(y16); y17=abs(y17); y18=abs(y18); y19=abs(y19); y20=abs(y20); y21=abs(y21); y22=abs(y22); if PHIcalculo<0.99 for PHI=0.6:passoPHI:PHIcalculo KP1=10^((K1(1)*log(T/1000)+(K1(2)/T)+K1(3)+K1(4)*T+K1(5)*T^2)); KP2=10^((K2(1)*log(T/1000)+(K2(2)/T)+K2(3)+K2(4)*T+K2(5)*T^2)); KP3=10^((K3(1)*log(T/1000)+(K3(2)/T)+K3(3)+K3(4)*T+K3(5)*T^2)); KP4=10^((K4(1)*log(T/1000)+(K4(2)/T)+K4(3)+K4(4)*T+K4(5)*T^2)); KP5=10^((K5(1)*log(T/1000)+(K5(2)/T)+K5(3)+K5(4)*T+K5(5)*T^2)); KP6=10^((K6(1)*log(T/1000)+(K6(2)/T)+K6(3)+K6(4)*T+K6(5)*T^2)); KP7=10^((K7(1)*log(T/1000)+(K7(2)/T)+K7(3)+K7(4)*T+K7(5)*T^2)); KP8=10^((K8(1)*log(T/1000)+(K8(2)/T)+K8(3)+K8(4)*T+K8(5)*T^2)); KP9=10^((K9(1)*log(T/1000)+(K9(2)/T)+K9(3)+K9(4)*T+K9(5)*T^2)); KP10=10^((K10(1)*log(T/1000)+(K10(2)/T)+K10(3)+K10(4)*T+K10(5)*T^2)); %KP11=10^((K11(1)*log(T/1000)+(K11(2)/T)+K11(3)+K11(4)*T+K11(5)*T^2)); %Ch3 KP12=10^((K12(1)*log(T/1000)+(K12(2)/T)+K12(3)+K12(4)*T+K12(5)*T^2)); KP13=10^((K13(1)*log(T/1000)+(K13(2)/T)+K13(3)+K13(4)*T+K13(5)*T^2)); KP14=10^((K14(1)*log(T/1000)+(K14(2)/T)+K14(3)+K14(4)*T+K14(5)*T^2)); KP15=10^((K15(1)*log(T/1000)+(K15(2)/T)+K15(3)+K15(4)*T+K15(5)*T^2)); KP16=10^((K16(1)*log(T/1000)+(K16(2)/T)+K16(3)+K16(4)*T+K16(5)*T^2)); KP17=10^((K17(1)*log(T/1000)+(K17(2)/T)+K17(3)+K17(4)*T+K17(5)*T^2)); KP18=10^((K18(1)*log(T/1000)+(K18(2)/T)+K18(3)+K18(4)*T+K18(5)*T^2)); KP19=10^((K19(1)*log(T/1000)+(K19(2)/T)+K19(3)+K19(4)*T+K19(5)*T^2)); c1=KP1/sqrt(P); c2=KP2/sqrt(P); c3=KP3; c4=KP4; c5=KP5*sqrt(P); c6=KP6*sqrt(P); c7=KP7*(P^2); c8=KP8*P; c9=KP9/(P^2); c10=KP10; %c11=1/KP11; % CH3 c12=1/KP12; c13=KP13; c14=1/KP14; c15=1/KP15; c16=KP16; c17=1/KP17; c18=1/KP18; c19=1/KP19; d1=1/(c1^2); d2=1/(c2^2); d3=1/(c3^2); d4=1/(c4^2); d5=1/(c5^2); d6=1/(c6^2); d8=1/(c8^2); % Inicio do programa for i=1:50 F5Y7=2*d1*y7; F6Y8=2*d2*y8; F7Y4=(-(y9^2)*d3)/(y4^2); F7Y9=2*y9*d3/y4; F8Y4=(-d4*(y10)^2)/(y4^2); F8Y10=(2*y10*d4)/y4; F9Y2=2*y2*d5/(y6^2); F9Y6=-2*(y2^2)*d5/(y6^3); F10Y1=(2*y1*d6)/(y5^2); F10Y5=-2*(y1^2)*d6/(y5^3); Fn1=(y1+y5+y12+y14+y15+y21+y22); Fn2=(2*y2+2*y6+y7+y9+4*y12+3*y13+y14+2*y15+y18+2*y20+y21); Fn3=(2*y1+y2+2*y4+y5+y8+y9+y10+y15+y16+2*y17+2*y18+y21); Fn4=(2*y3+y10+y13+y14+2*y16+y17+y19+y20+y22); F12Y2=-c7*y5*(y6^3)/(y2^2); F12Y5=c7*(y6^3)/y2; F12Y6=c7*y5*3*(y6^2)/y2; F13Y13=2*d8*y13/(y6^3); F13Y6=-(d8*3*(y13^2)/(y6^4)); F14Y6=-3*c9*y13*y12/(y6^4); F14Y12=c9*y13/(y6^3); F14Y13=c9*y12/(y6^3); F15Y1=c10*y2/y4; F15Y2=c10*y1/y4; F15Y4=-c10*y1*y2/(y4^2); F16Y3=c12*y18/y9; F16Y9=-c12*y3*y18/(y9^2); F16Y18=c12*y3/y9; F17Y9=-c13*y10*y18/(y9^2); F17Y10=c13*y18/y9; F17Y18=c13*y10/y9; F18Y4=c14*y9/y8; F18Y8=-c14*y9*y4/(y8^2); F18Y9=c14*y4/y8; F19Y3=c15*y8/y10; F19Y8=c15*y3/y10; F19Y10=-c15*y3*y8/(y10^2); F20Y5=-c16*y14*y9/(y5^2); F20Y9=c16*y14/y5; F20Y14=c16*y9/y5; F21Y5=c17*y6/y7; F21Y6=c17*y5/y7; F21Y7=-c17*y5*y6/(y7^2); F22Y2=-c18*y14*y9/(y2^2); F22Y9=c18*y14/y2; F22Y14=c18*y9/y2; % Matriz Jacobiana de derivadas parciais a=[N 0 0 0 N 0 0 0 0 0 Fn1 N 0 N N 0 0 0 0 0 N N; 0 2*N 0 0 0 2*N N 0 N 0 Fn2 4*N 3*N N 2*N 0 0 N 0 2*N N 0; 2*N N 0 2*N N 0 0 N N N Fn3 0 0 0 0 N 2*N 2*N 0 0 N 1; 0 0 2*N 0 0 0 0 0 0 N Fn4 0 N N 0 2*N N 0 N N 0 N; 0 0 0 0 0 -1 F5Y7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 -1 0 0 0 F6Y8 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 F7Y4 0 -1 0 0 F7Y9 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 -1 F8Y4 0 0 0 0 0 F8Y10 0 0 0 0 0 0 0 0 0 0 0 0; 0 F9Y2 0 -1 0 F9Y6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; F10Y1 0 0 -1 F10Y5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1; 0 F12Y2 0 0 F12Y5 F12Y6 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0; 0 0 -1 0 0 F13Y6 0 0 0 0 0 0 F13Y13 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 F14Y6 0 0 0 0 0 F14Y12 F14Y13 -1 0 0 0 0 0 0 0 0; F15Y1 F15Y2 0 F15Y4 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0; 0 0 F16Y3 0 0 0 0 0 F16Y9 0 0 0 0 0 0 -1 0 F16Y18 0 0 0 0; 0 0 0 0 0 0 0 0 F17Y9 F17Y10 0 0 0 0 0 0 -1 F17Y18 0 0 0 0; 0 0 0 F18Y4 0 0 0 F18Y8 F18Y9 0 0 0 0 0 0 0 0 -1 0 0 0 0; 0 0 F19Y3 0 0 0 0 F19Y8 0 F19Y10 0 0 0 0 0 0 0 0 -1 0 0 0; 0 0 0 0 F20Y5 0 0 0 F20Y9 0 0 0 0 F20Y14 0 0 0 0 0 -1 0 0; 0 0 0 0 F21Y5 F21Y6 F21Y7 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0; 0 F22Y2 0 0 0 0 0 0 F22Y9 0 0 0 0 F22Y14 0 0 0 0 0 0 0 -1]; % Termos independentes % Sistema de eq. não lineares 10x10 % Eqs. 3.56 3.57 3.58 3.59 3.68 3.69 3.70 3.71 3.72 3.73 pag 122 e 123 % Eq. 3.56 b1=((y1+y5+y12+y14+y15+y21+y22)*N)-(e*PHI*ALFA); % Eq. 3.57 b2=((2*y2+2*y6+y7+y9+4*y12+3*y13+y14+2*y15+y18+2*y20+y21)*N)-(e*PHI*BETA); % Eq. 3.58 b3=((2*y1+y2+2*y4+y5+y8+y9+y10+y15+y16+2*y17+2*y18+y21)*N)-((e*PHI*GAMA)+0.42); % Eq. 3.59 b4=((2*y3+y10+y13+y14+2*y16+y17+y19+y20+y22)*N)-((e*PHI*DELTA)+1.58); % Eq. 3.68 b5=(d1*y7^2)-y6; % Eq. 3.69 b6=(d2*y8^2)-y4; % Eq. 3.70 b7=((d3*y9^2)/y4)-y6; % Eq. 3.71 b8=((d4*(y10^2))/y4)-y3; % Eq. 3.72 b9=((d5*y2^2)/(y6^2))-y4; % Eq. 3.73 b10=((d6*y1^2)/(y5^2))-y4; % Eq.3.60 b11=y1+y2+y3+y4+y5+y6+y7+y8+y9+y10+y12+y13+y14+y15+y16+y17+y18+y19+y20+y21+y22-1; % Equacao do metano CO + 3H2 === CH4 + H20 b12=(c7*y5*(y6^3)/y2)-y12; % Equação da amônia 1/2 N2 + 3/2 H2 === NH3 b13=(((y13^2)*d8)/y6^3)-y3; % Equação do ácido cianídrico NH3 + CH4 === HCN+ 3H2 b14=(c9*y13*y12/(y6^3))-y14; % Equação do Formol C02 + H20 ===CH20+ O2 b15=(c10*y1*y2/y4)-y15; % Equação do CH3 CH3 + NO == HCN + H2O %b16=(c11*y14*y2/y10)-y16; % Equação do N2O N2O + OH == N2 + HO2 b16=(c12*y3*y18/y9)-y16; % Equação do NO2 HO2 + NO === NO2 + OH b17=(c13*y10*y18/y9)-y17; % Equação do HO2 0 + H02 === OH + 02 (substituta) b18=(c14*y9*y4/y8)-y18; % Equação do N N + NO ==== N2 + 0 b19=(c15*y3*y8/y10)-y19; % Equação do NH2 HCN + 0H == NH2 + CO b20=(c16*y14*y9/y5)-y20; % Equacao do HCO H + HCO == H2 + CO b21=(c17*y5*y6/y7)-y21; % Equacao do CN CN + H20 == HCN + OH b22=(c18*y14*y9/y2)-y22; % Equacao NH NH + 0 === NO + H %b24=(c19*y10*y7/y8)-y24; %\Montando o vetor B com as equacoes b=[-b1;-b2;-b3;-b4;-b5;-b6;-b7;-b8;-b9;-b10;-b11;-b12;-b13;-b14;-b15;-b16;-b17;-b18;-b19;-b20;-b21;-b22]; %Resolucao de sistemas lineares pelo MATLAB ao inves de x=A\B, tem-se x=inv(a)*b y=a\b; y1=y(1)+y1; y2=y(2)+y2; y3=y(3)+y3; y4=y(4)+y4; y5=y(5)+y5; y6=y(6)+y6; y7=y(7)+y7; y8=y(8)+y8; y9=y(9)+y9; y10=y(10)+y10; N=y(11)+N; y12=y(12)+y12; y13=y(13)+y13; y14=y(14)+y14; y15=y(15)+y15; y16=y(16)+y16; y17=y(17)+y17; y18=y(18)+y18; y19=y(19)+y19; y20=y(20)+y20; y21=y(21)+y21; y22=y(22)+y22; end end y1=abs(y1); y2=abs(y2); y3=abs(y3); y4=abs(y4); y5=abs(y5); y6=abs(y6); y7=abs(y7); y8=abs(y8); y9=abs(y9); y10=abs(y10); y12=abs(y12); y13=abs(y13); y14=abs(y14); y15=abs(y15); y16=abs(y16); y17=abs(y17); y18=abs(y18); y19=abs(y19); y20=abs(y20); y21=abs(y21); y22=abs(y22); end if PHIcalculo>=0.99 && PHIcalculo<=1.01 for PHI=0.6:passoPHI:0.99 KP1=10^((K1(1)*log(T/1000)+(K1(2)/T)+K1(3)+K1(4)*T+K1(5)*T^2)); KP2=10^((K2(1)*log(T/1000)+(K2(2)/T)+K2(3)+K2(4)*T+K2(5)*T^2)); KP3=10^((K3(1)*log(T/1000)+(K3(2)/T)+K3(3)+K3(4)*T+K3(5)*T^2)); KP4=10^((K4(1)*log(T/1000)+(K4(2)/T)+K4(3)+K4(4)*T+K4(5)*T^2)); KP5=10^((K5(1)*log(T/1000)+(K5(2)/T)+K5(3)+K5(4)*T+K5(5)*T^2)); KP6=10^((K6(1)*log(T/1000)+(K6(2)/T)+K6(3)+K6(4)*T+K6(5)*T^2)); KP7=10^((K7(1)*log(T/1000)+(K7(2)/T)+K7(3)+K7(4)*T+K7(5)*T^2)); KP8=10^((K8(1)*log(T/1000)+(K8(2)/T)+K8(3)+K8(4)*T+K8(5)*T^2)); KP9=10^((K9(1)*log(T/1000)+(K9(2)/T)+K9(3)+K9(4)*T+K9(5)*T^2)); KP10=10^((K10(1)*log(T/1000)+(K10(2)/T)+K10(3)+K10(4)*T+K10(5)*T^2)); %KP11=10^((K11(1)*log(T/1000)+(K11(2)/T)+K11(3)+K11(4)*T+K11(5)*T^2)); %Ch3 KP12=10^((K12(1)*log(T/1000)+(K12(2)/T)+K12(3)+K12(4)*T+K12(5)*T^2)); KP13=10^((K13(1)*log(T/1000)+(K13(2)/T)+K13(3)+K13(4)*T+K13(5)*T^2)); KP14=10^((K14(1)*log(T/1000)+(K14(2)/T)+K14(3)+K14(4)*T+K14(5)*T^2)); KP15=10^((K15(1)*log(T/1000)+(K15(2)/T)+K15(3)+K15(4)*T+K15(5)*T^2)); KP16=10^((K16(1)*log(T/1000)+(K16(2)/T)+K16(3)+K16(4)*T+K16(5)*T^2)); KP17=10^((K17(1)*log(T/1000)+(K17(2)/T)+K17(3)+K17(4)*T+K17(5)*T^2)); KP18=10^((K18(1)*log(T/1000)+(K18(2)/T)+K18(3)+K18(4)*T+K18(5)*T^2)); KP19=10^((K19(1)*log(T/1000)+(K19(2)/T)+K19(3)+K19(4)*T+K19(5)*T^2)); c1=KP1/sqrt(P); c2=KP2/sqrt(P); c3=KP3; c4=KP4; c5=KP5*sqrt(P); c6=KP6*sqrt(P); c7=KP7*(P^2); c8=KP8*P; c9=KP9/(P^2); c10=KP10; %c11=1/KP11; % CH3 c12=1/KP12; c13=KP13; c14=1/KP14; c15=1/KP15; c16=KP16; c17=1/KP17; c18=1/KP18; c19=1/KP19; d1=1/(c1^2); d2=1/(c2^2); d3=1/(c3^2); d4=1/(c4^2); d5=1/(c5^2); d6=1/(c6^2); d8=1/(c8^2); % Inicio do programa for i=1:50 F5Y7=2*d1*y7; F6Y8=2*d2*y8; F7Y4=(-(y9^2)*d3)/(y4^2); F7Y9=2*y9*d3/y4; F8Y4=(-d4*(y10)^2)/(y4^2); F8Y10=(2*y10*d4)/y4; F9Y2=2*y2*d5/(y6^2); F9Y6=-2*(y2^2)*d5/(y6^3); F10Y1=(2*y1*d6)/(y5^2); F10Y5=-2*(y1^2)*d6/(y5^3); Fn1=(y1+y5+y12+y14+y15+y21+y22); Fn2=(2*y2+2*y6+y7+y9+4*y12+3*y13+y14+2*y15+y18+2*y20+y21); Fn3=(2*y1+y2+2*y4+y5+y8+y9+y10+y15+y16+2*y17+2*y18+y21); Fn4=(2*y3+y10+y13+y14+2*y16+y17+y19+y20+y22); F12Y2=-c7*y5*(y6^3)/(y2^2); F12Y5=c7*(y6^3)/y2; F12Y6=c7*y5*3*(y6^2)/y2; F13Y13=2*d8*y13/(y6^3); F13Y6=-(d8*3*(y13^2)/(y6^4)); F14Y6=-3*c9*y13*y12/(y6^4); F14Y12=c9*y13/(y6^3); F14Y13=c9*y12/(y6^3); F15Y1=c10*y2/y4; F15Y2=c10*y1/y4; F15Y4=-c10*y1*y2/(y4^2); F16Y3=c12*y18/y9; F16Y9=-c12*y3*y18/(y9^2); F16Y18=c12*y3/y9; F17Y9=-c13*y10*y18/(y9^2); F17Y10=c13*y18/y9; F17Y18=c13*y10/y9; F18Y4=c14*y9/y8; F18Y8=-c14*y9*y4/(y8^2); F18Y9=c14*y4/y8; F19Y3=c15*y8/y10; F19Y8=c15*y3/y10; F19Y10=-c15*y3*y8/(y10^2); F20Y5=-c16*y14*y9/(y5^2); F20Y9=c16*y14/y5; F20Y14=c16*y9/y5; F21Y5=c17*y6/y7; F21Y6=c17*y5/y7; F21Y7=-c17*y5*y6/(y7^2); F22Y2=-c18*y14*y9/(y2^2); F22Y9=c18*y14/y2; F22Y14=c18*y9/y2; % Matriz Jacobiana de derivadas parciais a=[N 0 0 0 N 0 0 0 0 0 Fn1 N 0 N N 0 0 0 0 0 N N; 0 2*N 0 0 0 2*N N 0 N 0 Fn2 4*N 3*N N 2*N 0 0 N 0 2*N N 0; 2*N N 0 2*N N 0 0 N N N Fn3 0 0 0 0 N 2*N 2*N 0 0 N 1; 0 0 2*N 0 0 0 0 0 0 N Fn4 0 N N 0 2*N N 0 N N 0 N; 0 0 0 0 0 -1 F5Y7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 -1 0 0 0 F6Y8 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 F7Y4 0 -1 0 0 F7Y9 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 -1 F8Y4 0 0 0 0 0 F8Y10 0 0 0 0 0 0 0 0 0 0 0 0; 0 F9Y2 0 -1 0 F9Y6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; F10Y1 0 0 -1 F10Y5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1; 0 F12Y2 0 0 F12Y5 F12Y6 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0; 0 0 -1 0 0 F13Y6 0 0 0 0 0 0 F13Y13 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 F14Y6 0 0 0 0 0 F14Y12 F14Y13 -1 0 0 0 0 0 0 0 0; F15Y1 F15Y2 0 F15Y4 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0; 0 0 F16Y3 0 0 0 0 0 F16Y9 0 0 0 0 0 0 -1 0 F16Y18 0 0 0 0; 0 0 0 0 0 0 0 0 F17Y9 F17Y10 0 0 0 0 0 0 -1 F17Y18 0 0 0 0; 0 0 0 F18Y4 0 0 0 F18Y8 F18Y9 0 0 0 0 0 0 0 0 -1 0 0 0 0; 0 0 F19Y3 0 0 0 0 F19Y8 0 F19Y10 0 0 0 0 0 0 0 0 -1 0 0 0; 0 0 0 0 F20Y5 0 0 0 F20Y9 0 0 0 0 F20Y14 0 0 0 0 0 -1 0 0; 0 0 0 0 F21Y5 F21Y6 F21Y7 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0; 0 F22Y2 0 0 0 0 0 0 F22Y9 0 0 0 0 F22Y14 0 0 0 0 0 0 0 -1]; % Termos independentes % Sistema de eq. não lineares 10x10 % Eqs. 3.56 3.57 3.58 3.59 3.68 3.69 3.70 3.71 3.72 3.73 pag 122 e 123 % Eq. 3.56 b1=((y1+y5+y12+y14+y15+y21+y22)*N)-(e*PHI*ALFA); % Eq. 3.57 b2=((2*y2+2*y6+y7+y9+4*y12+3*y13+y14+2*y15+y18+2*y20+y21)*N)-(e*PHI*BETA); % Eq. 3.58 b3=((2*y1+y2+2*y4+y5+y8+y9+y10+y15+y16+2*y17+2*y18+y21)*N)-((e*PHI*GAMA)+0.42); % Eq. 3.59 b4=((2*y3+y10+y13+y14+2*y16+y17+y19+y20+y22)*N)-((e*PHI*DELTA)+1.58); % Eq. 3.68 b5=(d1*y7^2)-y6; % Eq. 3.69 b6=(d2*y8^2)-y4; % Eq. 3.70 b7=((d3*y9^2)/y4)-y6; % Eq. 3.71 b8=((d4*(y10^2))/y4)-y3; % Eq. 3.72 b9=((d5*y2^2)/(y6^2))-y4; % Eq. 3.73 b10=((d6*y1^2)/(y5^2))-y4; % Eq.3.60 b11=y1+y2+y3+y4+y5+y6+y7+y8+y9+y10+y12+y13+y14+y15+y16+y17+y18+y19+y20+y21+y22-1; % Equacao do metano CO + 3H2 === CH4 + H20 b12=(c7*y5*(y6^3)/y2)-y12; % Equação da amônia 1/2 N2 + 3/2 H2 === NH3 b13=(((y13^2)*d8)/y6^3)-y3; % Equação do ácido cianídrico NH3 + CH4 === HCN+ 3H2 b14=(c9*y13*y12/(y6^3))-y14; % Equação do Formol C02 + H20 ===CH20+ O2 b15=(c10*y1*y2/y4)-y15; % Equação do CH3 CH3 + NO == HCN + H2O %b16=(c11*y14*y2/y10)-y16; % Equação do N2O N2O + OH == N2 + HO2 b16=(c12*y3*y18/y9)-y16; % Equação do NO2 HO2 + NO === NO2 + OH b17=(c13*y10*y18/y9)-y17; % Equação do HO2 0 + H02 === OH + 02 (substituta) b18=(c14*y9*y4/y8)-y18; % Equação do N N + NO ==== N2 + 0 b19=(c15*y3*y8/y10)-y19; % Equação do NH2 HCN + 0H == NH2 + CO b20=(c16*y14*y9/y5)-y20; % Equacao do HCO H + HCO == H2 + CO b21=(c17*y5*y6/y7)-y21; % Equacao do CN CN + H20 == HCN + OH b22=(c18*y14*y9/y2)-y22; % Equacao NH NH + 0 === NO + H %b24=(c19*y10*y7/y8)-y24; %\Montando o vetor B com as equacoes b=[-b1;-b2;-b3;-b4;-b5;-b6;-b7;-b8;-b9;-b10;-b11;-b12;-b13;-b14;-b15;-b16;-b17;-b18;-b19;-b20;-b21;-b22]; %Resolucao de sistemas lineares pelo MATLAB ao inves de x=A\B, tem-se x=inv(a)*b y=a\b; y1=y(1)+y1; y2=y(2)+y2; y3=y(3)+y3; y4=y(4)+y4; y5=y(5)+y5; y6=y(6)+y6; y7=y(7)+y7; y8=y(8)+y8; y9=y(9)+y9; y10=y(10)+y10; N=y(11)+N; y12=y(12)+y12; y13=y(13)+y13; y14=y(14)+y14; y15=y(15)+y15; y16=y(16)+y16; y17=y(17)+y17; y18=y(18)+y18; y19=y(19)+y19; y20=y(20)+y20; y21=y(21)+y21; y22=y(22)+y22; end end for PHI=0.99:passoPHI2:PHIcalculo KP1=10^((K1(1)*log(T/1000)+(K1(2)/T)+K1(3)+K1(4)*T+K1(5)*T^2)); KP2=10^((K2(1)*log(T/1000)+(K2(2)/T)+K2(3)+K2(4)*T+K2(5)*T^2)); KP3=10^((K3(1)*log(T/1000)+(K3(2)/T)+K3(3)+K3(4)*T+K3(5)*T^2)); KP4=10^((K4(1)*log(T/1000)+(K4(2)/T)+K4(3)+K4(4)*T+K4(5)*T^2)); KP5=10^((K5(1)*log(T/1000)+(K5(2)/T)+K5(3)+K5(4)*T+K5(5)*T^2)); KP6=10^((K6(1)*log(T/1000)+(K6(2)/T)+K6(3)+K6(4)*T+K6(5)*T^2)); KP7=10^((K7(1)*log(T/1000)+(K7(2)/T)+K7(3)+K7(4)*T+K7(5)*T^2)); KP8=10^((K8(1)*log(T/1000)+(K8(2)/T)+K8(3)+K8(4)*T+K8(5)*T^2)); KP9=10^((K9(1)*log(T/1000)+(K9(2)/T)+K9(3)+K9(4)*T+K9(5)*T^2)); KP10=10^((K10(1)*log(T/1000)+(K10(2)/T)+K10(3)+K10(4)*T+K10(5)*T^2)); %KP11=10^((K11(1)*log(T/1000)+(K11(2)/T)+K11(3)+K11(4)*T+K11(5)*T^2)); %Ch3 KP12=10^((K12(1)*log(T/1000)+(K12(2)/T)+K12(3)+K12(4)*T+K12(5)*T^2)); KP13=10^((K13(1)*log(T/1000)+(K13(2)/T)+K13(3)+K13(4)*T+K13(5)*T^2)); KP14=10^((K14(1)*log(T/1000)+(K14(2)/T)+K14(3)+K14(4)*T+K14(5)*T^2)); KP15=10^((K15(1)*log(T/1000)+(K15(2)/T)+K15(3)+K15(4)*T+K15(5)*T^2)); KP16=10^((K16(1)*log(T/1000)+(K16(2)/T)+K16(3)+K16(4)*T+K16(5)*T^2)); KP17=10^((K17(1)*log(T/1000)+(K17(2)/T)+K17(3)+K17(4)*T+K17(5)*T^2)); KP18=10^((K18(1)*log(T/1000)+(K18(2)/T)+K18(3)+K18(4)*T+K18(5)*T^2)); KP19=10^((K19(1)*log(T/1000)+(K19(2)/T)+K19(3)+K19(4)*T+K19(5)*T^2)); c1=KP1/sqrt(P); c2=KP2/sqrt(P); c3=KP3; c4=KP4; c5=KP5*sqrt(P); c6=KP6*sqrt(P); c7=KP7*(P^2); c8=KP8*P; c9=KP9/(P^2); c10=KP10; %c11=1/KP11; % CH3 c12=1/KP12; c13=KP13; c14=1/KP14; c15=1/KP15; c16=KP16; c17=1/KP17; c18=1/KP18; c19=1/KP19; d1=1/(c1^2); d2=1/(c2^2); d3=1/(c3^2); d4=1/(c4^2); d5=1/(c5^2); d6=1/(c6^2); d8=1/(c8^2); % Inicio do programa for i=1:50 F5Y7=2*d1*y7; F6Y8=2*d2*y8; F7Y4=(-(y9^2)*d3)/(y4^2); F7Y9=2*y9*d3/y4; F8Y4=(-d4*(y10)^2)/(y4^2); F8Y10=(2*y10*d4)/y4; F9Y2=2*y2*d5/(y6^2); F9Y6=-2*(y2^2)*d5/(y6^3); F10Y1=(2*y1*d6)/(y5^2); F10Y5=-2*(y1^2)*d6/(y5^3); Fn1=(y1+y5+y12+y14+y15+y21+y22); Fn2=(2*y2+2*y6+y7+y9+4*y12+3*y13+y14+2*y15+y18+2*y20+y21); Fn3=(2*y1+y2+2*y4+y5+y8+y9+y10+y15+y16+2*y17+2*y18+y21); Fn4=(2*y3+y10+y13+y14+2*y16+y17+y19+y20+y22); F12Y2=-c7*y5*(y6^3)/(y2^2); F12Y5=c7*(y6^3)/y2; F12Y6=c7*y5*3*(y6^2)/y2; F13Y13=2*d8*y13/(y6^3); F13Y6=-(d8*3*(y13^2)/(y6^4)); F14Y6=-3*c9*y13*y12/(y6^4); F14Y12=c9*y13/(y6^3); F14Y13=c9*y12/(y6^3); F15Y1=c10*y2/y4; F15Y2=c10*y1/y4; F15Y4=-c10*y1*y2/(y4^2); F16Y3=c12*y18/y9; F16Y9=-c12*y3*y18/(y9^2); F16Y18=c12*y3/y9; F17Y9=-c13*y10*y18/(y9^2); F17Y10=c13*y18/y9; F17Y18=c13*y10/y9; F18Y4=c14*y9/y8; F18Y8=-c14*y9*y4/(y8^2); F18Y9=c14*y4/y8; F19Y3=c15*y8/y10; F19Y8=c15*y3/y10; F19Y10=-c15*y3*y8/(y10^2); F20Y5=-c16*y14*y9/(y5^2); F20Y9=c16*y14/y5; F20Y14=c16*y9/y5; F21Y5=c17*y6/y7; F21Y6=c17*y5/y7; F21Y7=-c17*y5*y6/(y7^2); F22Y2=-c18*y14*y9/(y2^2); F22Y9=c18*y14/y2; F22Y14=c18*y9/y2; % Matriz Jacobiana de derivadas parciais a=[N 0 0 0 N 0 0 0 0 0 Fn1 N 0 N N 0 0 0 0 0 N N; 0 2*N 0 0 0 2*N N 0 N 0 Fn2 4*N 3*N N 2*N 0 0 N 0 2*N N 0; 2*N N 0 2*N N 0 0 N N N Fn3 0 0 0 0 N 2*N 2*N 0 0 N 1; 0 0 2*N 0 0 0 0 0 0 N Fn4 0 N N 0 2*N N 0 N N 0 N; 0 0 0 0 0 -1 F5Y7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 -1 0 0 0 F6Y8 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 F7Y4 0 -1 0 0 F7Y9 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 -1 F8Y4 0 0 0 0 0 F8Y10 0 0 0 0 0 0 0 0 0 0 0 0; 0 F9Y2 0 -1 0 F9Y6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; F10Y1 0 0 -1 F10Y5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1; 0 F12Y2 0 0 F12Y5 F12Y6 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0; 0 0 -1 0 0 F13Y6 0 0 0 0 0 0 F13Y13 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 F14Y6 0 0 0 0 0 F14Y12 F14Y13 -1 0 0 0 0 0 0 0 0; F15Y1 F15Y2 0 F15Y4 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0; 0 0 F16Y3 0 0 0 0 0 F16Y9 0 0 0 0 0 0 -1 0 F16Y18 0 0 0 0; 0 0 0 0 0 0 0 0 F17Y9 F17Y10 0 0 0 0 0 0 -1 F17Y18 0 0 0 0; 0 0 0 F18Y4 0 0 0 F18Y8 F18Y9 0 0 0 0 0 0 0 0 -1 0 0 0 0; 0 0 F19Y3 0 0 0 0 F19Y8 0 F19Y10 0 0 0 0 0 0 0 0 -1 0 0 0; 0 0 0 0 F20Y5 0 0 0 F20Y9 0 0 0 0 F20Y14 0 0 0 0 0 -1 0 0; 0 0 0 0 F21Y5 F21Y6 F21Y7 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0; 0 F22Y2 0 0 0 0 0 0 F22Y9 0 0 0 0 F22Y14 0 0 0 0 0 0 0 -1]; % Termos independentes % Sistema de eq. não lineares 10x10 % Eqs. 3.56 3.57 3.58 3.59 3.68 3.69 3.70 3.71 3.72 3.73 pag 122 e 123 % Eq. 3.56 b1=((y1+y5+y12+y14+y15+y21+y22)*N)-(e*PHI*ALFA); % Eq. 3.57 b2=((2*y2+2*y6+y7+y9+4*y12+3*y13+y14+2*y15+y18+2*y20+y21)*N)-(e*PHI*BETA); % Eq. 3.58 b3=((2*y1+y2+2*y4+y5+y8+y9+y10+y15+y16+2*y17+2*y18+y21)*N)-((e*PHI*GAMA)+0.42); % Eq. 3.59 b4=((2*y3+y10+y13+y14+2*y16+y17+y19+y20+y22)*N)-((e*PHI*DELTA)+1.58); % Eq. 3.68 b5=(d1*y7^2)-y6; % Eq. 3.69 b6=(d2*y8^2)-y4; % Eq. 3.70 b7=((d3*y9^2)/y4)-y6; % Eq. 3.71 b8=((d4*(y10^2))/y4)-y3; % Eq. 3.72 b9=((d5*y2^2)/(y6^2))-y4; % Eq. 3.73 b10=((d6*y1^2)/(y5^2))-y4; % Eq.3.60 b11=y1+y2+y3+y4+y5+y6+y7+y8+y9+y10+y12+y13+y14+y15+y16+y17+y18+y19+y20+y21+y22-1; % Equacao do metano CO + 3H2 === CH4 + H20 b12=(c7*y5*(y6^3)/y2)-y12; % Equação da amônia 1/2 N2 + 3/2 H2 === NH3 b13=(((y13^2)*d8)/y6^3)-y3; % Equação do ácido cianídrico NH3 + CH4 === HCN+ 3H2 b14=(c9*y13*y12/(y6^3))-y14; % Equação do Formol C02 + H20 ===CH20+ O2 b15=(c10*y1*y2/y4)-y15; % Equação do CH3 CH3 + NO == HCN + H2O %b16=(c11*y14*y2/y10)-y16; % Equação do N2O N2O + OH == N2 + HO2 b16=(c12*y3*y18/y9)-y16; % Equação do NO2 HO2 + NO === NO2 + OH b17=(c13*y10*y18/y9)-y17; % Equação do HO2 0 + H02 === OH + 02 (substituta) b18=(c14*y9*y4/y8)-y18; % Equação do N N + NO ==== N2 + 0 b19=(c15*y3*y8/y10)-y19; % Equação do NH2 HCN + 0H == NH2 + CO b20=(c16*y14*y9/y5)-y20; % Equacao do HCO H + HCO == H2 + CO b21=(c17*y5*y6/y7)-y21; % Equacao do CN CN + H20 == HCN + OH b22=(c18*y14*y9/y2)-y22; % Equacao NH NH + 0 === NO + H %b24=(c19*y10*y7/y8)-y24; %\Montando o vetor B com as equacoes b=[-b1;-b2;-b3;-b4;-b5;-b6;-b7;-b8;-b9;-b10;-b11;-b12;-b13;-b14;-b15;-b16;-b17;-b18;-b19;-b20;-b21;-b22]; %Resolucao de sistemas lineares pelo MATLAB ao inves de x=A\B, tem-se x=inv(a)*b y=a\b; y1=y(1)+y1; y2=y(2)+y2; y3=y(3)+y3; y4=y(4)+y4; y5=y(5)+y5; y6=y(6)+y6; y7=y(7)+y7; y8=y(8)+y8; y9=y(9)+y9; y10=y(10)+y10; N=y(11)+N; y12=y(12)+y12; y13=y(13)+y13; y14=y(14)+y14; y15=y(15)+y15; y16=y(16)+y16; y17=y(17)+y17; y18=y(18)+y18; y19=y(19)+y19; y20=y(20)+y20; y21=y(21)+y21; y22=y(22)+y22; end end y1=abs(y1); y2=abs(y2); y3=abs(y3); y4=abs(y4); y5=abs(y5); y6=abs(y6); y7=abs(y7); y8=abs(y8); y9=abs(y9); y10=abs(y10); y12=abs(y12); y13=abs(y13); y14=abs(y14); y15=abs(y15); y16=abs(y16); y17=abs(y17); y18=abs(y18); y19=abs(y19); y20=abs(y20); y21=abs(y21); y22=abs(y22); end if PHIcalculo>=1.01 for PHI=0.6:passoPHI:0.99 KP1=10^((K1(1)*log(T/1000)+(K1(2)/T)+K1(3)+K1(4)*T+K1(5)*T^2)); KP2=10^((K2(1)*log(T/1000)+(K2(2)/T)+K2(3)+K2(4)*T+K2(5)*T^2)); KP3=10^((K3(1)*log(T/1000)+(K3(2)/T)+K3(3)+K3(4)*T+K3(5)*T^2)); KP4=10^((K4(1)*log(T/1000)+(K4(2)/T)+K4(3)+K4(4)*T+K4(5)*T^2)); KP5=10^((K5(1)*log(T/1000)+(K5(2)/T)+K5(3)+K5(4)*T+K5(5)*T^2)); KP6=10^((K6(1)*log(T/1000)+(K6(2)/T)+K6(3)+K6(4)*T+K6(5)*T^2)); KP7=10^((K7(1)*log(T/1000)+(K7(2)/T)+K7(3)+K7(4)*T+K7(5)*T^2)); KP8=10^((K8(1)*log(T/1000)+(K8(2)/T)+K8(3)+K8(4)*T+K8(5)*T^2)); KP9=10^((K9(1)*log(T/1000)+(K9(2)/T)+K9(3)+K9(4)*T+K9(5)*T^2)); KP10=10^((K10(1)*log(T/1000)+(K10(2)/T)+K10(3)+K10(4)*T+K10(5)*T^2)); %KP11=10^((K11(1)*log(T/1000)+(K11(2)/T)+K11(3)+K11(4)*T+K11(5)*T^2)); %Ch3 KP12=10^((K12(1)*log(T/1000)+(K12(2)/T)+K12(3)+K12(4)*T+K12(5)*T^2)); KP13=10^((K13(1)*log(T/1000)+(K13(2)/T)+K13(3)+K13(4)*T+K13(5)*T^2)); KP14=10^((K14(1)*log(T/1000)+(K14(2)/T)+K14(3)+K14(4)*T+K14(5)*T^2)); KP15=10^((K15(1)*log(T/1000)+(K15(2)/T)+K15(3)+K15(4)*T+K15(5)*T^2)); KP16=10^((K16(1)*log(T/1000)+(K16(2)/T)+K16(3)+K16(4)*T+K16(5)*T^2)); KP17=10^((K17(1)*log(T/1000)+(K17(2)/T)+K17(3)+K17(4)*T+K17(5)*T^2)); KP18=10^((K18(1)*log(T/1000)+(K18(2)/T)+K18(3)+K18(4)*T+K18(5)*T^2)); KP19=10^((K19(1)*log(T/1000)+(K19(2)/T)+K19(3)+K19(4)*T+K19(5)*T^2)); c1=KP1/sqrt(P); c2=KP2/sqrt(P); c3=KP3; c4=KP4; c5=KP5*sqrt(P); c6=KP6*sqrt(P); c7=KP7*(P^2); c8=KP8*P; c9=KP9/(P^2); c10=KP10; %c11=1/KP11; % CH3 c12=1/KP12; c13=KP13; c14=1/KP14; c15=1/KP15; c16=KP16; c17=1/KP17; c18=1/KP18; c19=1/KP19; d1=1/(c1^2); d2=1/(c2^2); d3=1/(c3^2); d4=1/(c4^2); d5=1/(c5^2); d6=1/(c6^2); d8=1/(c8^2); % Inicio do programa for i=1:50 F5Y7=2*d1*y7; F6Y8=2*d2*y8; F7Y4=(-(y9^2)*d3)/(y4^2); F7Y9=2*y9*d3/y4; F8Y4=(-d4*(y10)^2)/(y4^2); F8Y10=(2*y10*d4)/y4; F9Y2=2*y2*d5/(y6^2); F9Y6=-2*(y2^2)*d5/(y6^3); F10Y1=(2*y1*d6)/(y5^2); F10Y5=-2*(y1^2)*d6/(y5^3); Fn1=(y1+y5+y12+y14+y15+y21+y22); Fn2=(2*y2+2*y6+y7+y9+4*y12+3*y13+y14+2*y15+y18+2*y20+y21); Fn3=(2*y1+y2+2*y4+y5+y8+y9+y10+y15+y16+2*y17+2*y18+y21); Fn4=(2*y3+y10+y13+y14+2*y16+y17+y19+y20+y22); F12Y2=-c7*y5*(y6^3)/(y2^2); F12Y5=c7*(y6^3)/y2; F12Y6=c7*y5*3*(y6^2)/y2; F13Y13=2*d8*y13/(y6^3); F13Y6=-(d8*3*(y13^2)/(y6^4)); F14Y6=-3*c9*y13*y12/(y6^4); F14Y12=c9*y13/(y6^3); F14Y13=c9*y12/(y6^3); F15Y1=c10*y2/y4; F15Y2=c10*y1/y4; F15Y4=-c10*y1*y2/(y4^2); F16Y3=c12*y18/y9; F16Y9=-c12*y3*y18/(y9^2); F16Y18=c12*y3/y9; F17Y9=-c13*y10*y18/(y9^2); F17Y10=c13*y18/y9; F17Y18=c13*y10/y9; F18Y4=c14*y9/y8; F18Y8=-c14*y9*y4/(y8^2); F18Y9=c14*y4/y8; F19Y3=c15*y8/y10; F19Y8=c15*y3/y10; F19Y10=-c15*y3*y8/(y10^2); F20Y5=-c16*y14*y9/(y5^2); F20Y9=c16*y14/y5; F20Y14=c16*y9/y5; F21Y5=c17*y6/y7; F21Y6=c17*y5/y7; F21Y7=-c17*y5*y6/(y7^2); F22Y2=-c18*y14*y9/(y2^2); F22Y9=c18*y14/y2; F22Y14=c18*y9/y2; % Matriz Jacobiana de derivadas parciais a=[N 0 0 0 N 0 0 0 0 0 Fn1 N 0 N N 0 0 0 0 0 N N; 0 2*N 0 0 0 2*N N 0 N 0 Fn2 4*N 3*N N 2*N 0 0 N 0 2*N N 0; 2*N N 0 2*N N 0 0 N N N Fn3 0 0 0 0 N 2*N 2*N 0 0 N 1; 0 0 2*N 0 0 0 0 0 0 N Fn4 0 N N 0 2*N N 0 N N 0 N; 0 0 0 0 0 -1 F5Y7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 -1 0 0 0 F6Y8 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 F7Y4 0 -1 0 0 F7Y9 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 -1 F8Y4 0 0 0 0 0 F8Y10 0 0 0 0 0 0 0 0 0 0 0 0; 0 F9Y2 0 -1 0 F9Y6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; F10Y1 0 0 -1 F10Y5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1; 0 F12Y2 0 0 F12Y5 F12Y6 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0; 0 0 -1 0 0 F13Y6 0 0 0 0 0 0 F13Y13 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 F14Y6 0 0 0 0 0 F14Y12 F14Y13 -1 0 0 0 0 0 0 0 0; F15Y1 F15Y2 0 F15Y4 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0; 0 0 F16Y3 0 0 0 0 0 F16Y9 0 0 0 0 0 0 -1 0 F16Y18 0 0 0 0; 0 0 0 0 0 0 0 0 F17Y9 F17Y10 0 0 0 0 0 0 -1 F17Y18 0 0 0 0; 0 0 0 F18Y4 0 0 0 F18Y8 F18Y9 0 0 0 0 0 0 0 0 -1 0 0 0 0; 0 0 F19Y3 0 0 0 0 F19Y8 0 F19Y10 0 0 0 0 0 0 0 0 -1 0 0 0; 0 0 0 0 F20Y5 0 0 0 F20Y9 0 0 0 0 F20Y14 0 0 0 0 0 -1 0 0; 0 0 0 0 F21Y5 F21Y6 F21Y7 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0; 0 F22Y2 0 0 0 0 0 0 F22Y9 0 0 0 0 F22Y14 0 0 0 0 0 0 0 -1]; % Termos independentes % Sistema de eq. não lineares 10x10 % Eqs. 3.56 3.57 3.58 3.59 3.68 3.69 3.70 3.71 3.72 3.73 pag 122 e 123 % Eq. 3.56 b1=((y1+y5+y12+y14+y15+y21+y22)*N)-(e*PHI*ALFA); % Eq. 3.57 b2=((2*y2+2*y6+y7+y9+4*y12+3*y13+y14+2*y15+y18+2*y20+y21)*N)-(e*PHI*BETA); % Eq. 3.58 b3=((2*y1+y2+2*y4+y5+y8+y9+y10+y15+y16+2*y17+2*y18+y21)*N)-((e*PHI*GAMA)+0.42); % Eq. 3.59 b4=((2*y3+y10+y13+y14+2*y16+y17+y19+y20+y22)*N)-((e*PHI*DELTA)+1.58); % Eq. 3.68 b5=(d1*y7^2)-y6; % Eq. 3.69 b6=(d2*y8^2)-y4; % Eq. 3.70 b7=((d3*y9^2)/y4)-y6; % Eq. 3.71 b8=((d4*(y10^2))/y4)-y3; % Eq. 3.72 b9=((d5*y2^2)/(y6^2))-y4; % Eq. 3.73 b10=((d6*y1^2)/(y5^2))-y4; % Eq.3.60 b11=y1+y2+y3+y4+y5+y6+y7+y8+y9+y10+y12+y13+y14+y15+y16+y17+y18+y19+y20+y21+y22-1; % Equacao do metano CO + 3H2 === CH4 + H20 b12=(c7*y5*(y6^3)/y2)-y12; % Equação da amônia 1/2 N2 + 3/2 H2 === NH3 b13=(((y13^2)*d8)/y6^3)-y3; % Equação do ácido cianídrico NH3 + CH4 === HCN+ 3H2 b14=(c9*y13*y12/(y6^3))-y14; % Equação do Formol C02 + H20 ===CH20+ O2 b15=(c10*y1*y2/y4)-y15; % Equação do CH3 CH3 + NO == HCN + H2O %b16=(c11*y14*y2/y10)-y16; % Equação do N2O N2O + OH == N2 + HO2 b16=(c12*y3*y18/y9)-y16; % Equação do NO2 HO2 + NO === NO2 + OH b17=(c13*y10*y18/y9)-y17; % Equação do HO2 0 + H02 === OH + 02 (substituta) b18=(c14*y9*y4/y8)-y18; % Equação do N N + NO ==== N2 + 0 b19=(c15*y3*y8/y10)-y19; % Equação do NH2 HCN + 0H == NH2 + CO b20=(c16*y14*y9/y5)-y20; % Equacao do HCO H + HCO == H2 + CO b21=(c17*y5*y6/y7)-y21; % Equacao do CN CN + H20 == HCN + OH b22=(c18*y14*y9/y2)-y22; % Equacao NH NH + 0 === NO + H %b24=(c19*y10*y7/y8)-y24; %\Montando o vetor B com as equacoes b=[-b1;-b2;-b3;-b4;-b5;-b6;-b7;-b8;-b9;-b10;-b11;-b12;-b13;-b14;-b15;-b16;-b17;-b18;-b19;-b20;-b21;-b22]; %Resolucao de sistemas lineares pelo MATLAB ao inves de x=A\B, tem-se x=inv(a)*b y=a\b; y1=y(1)+y1; y2=y(2)+y2; y3=y(3)+y3; y4=y(4)+y4; y5=y(5)+y5; y6=y(6)+y6; y7=y(7)+y7; y8=y(8)+y8; y9=y(9)+y9; y10=y(10)+y10; N=y(11)+N; y12=y(12)+y12; y13=y(13)+y13; y14=y(14)+y14; y15=y(15)+y15; y16=y(16)+y16; y17=y(17)+y17; y18=y(18)+y18; y19=y(19)+y19; y20=y(20)+y20; y21=y(21)+y21; y22=y(22)+y22; end end for PHI=0.99:passoPHI2:1.01 KP1=10^((K1(1)*log(T/1000)+(K1(2)/T)+K1(3)+K1(4)*T+K1(5)*T^2)); KP2=10^((K2(1)*log(T/1000)+(K2(2)/T)+K2(3)+K2(4)*T+K2(5)*T^2)); KP3=10^((K3(1)*log(T/1000)+(K3(2)/T)+K3(3)+K3(4)*T+K3(5)*T^2)); KP4=10^((K4(1)*log(T/1000)+(K4(2)/T)+K4(3)+K4(4)*T+K4(5)*T^2)); KP5=10^((K5(1)*log(T/1000)+(K5(2)/T)+K5(3)+K5(4)*T+K5(5)*T^2)); KP6=10^((K6(1)*log(T/1000)+(K6(2)/T)+K6(3)+K6(4)*T+K6(5)*T^2)); KP7=10^((K7(1)*log(T/1000)+(K7(2)/T)+K7(3)+K7(4)*T+K7(5)*T^2)); KP8=10^((K8(1)*log(T/1000)+(K8(2)/T)+K8(3)+K8(4)*T+K8(5)*T^2)); KP9=10^((K9(1)*log(T/1000)+(K9(2)/T)+K9(3)+K9(4)*T+K9(5)*T^2)); KP10=10^((K10(1)*log(T/1000)+(K10(2)/T)+K10(3)+K10(4)*T+K10(5)*T^2)); %KP11=10^((K11(1)*log(T/1000)+(K11(2)/T)+K11(3)+K11(4)*T+K11(5)*T^2)); %Ch3 KP12=10^((K12(1)*log(T/1000)+(K12(2)/T)+K12(3)+K12(4)*T+K12(5)*T^2)); KP13=10^((K13(1)*log(T/1000)+(K13(2)/T)+K13(3)+K13(4)*T+K13(5)*T^2)); KP14=10^((K14(1)*log(T/1000)+(K14(2)/T)+K14(3)+K14(4)*T+K14(5)*T^2)); KP15=10^((K15(1)*log(T/1000)+(K15(2)/T)+K15(3)+K15(4)*T+K15(5)*T^2)); KP16=10^((K16(1)*log(T/1000)+(K16(2)/T)+K16(3)+K16(4)*T+K16(5)*T^2)); KP17=10^((K17(1)*log(T/1000)+(K17(2)/T)+K17(3)+K17(4)*T+K17(5)*T^2)); KP18=10^((K18(1)*log(T/1000)+(K18(2)/T)+K18(3)+K18(4)*T+K18(5)*T^2)); KP19=10^((K19(1)*log(T/1000)+(K19(2)/T)+K19(3)+K19(4)*T+K19(5)*T^2)); c1=KP1/sqrt(P); c2=KP2/sqrt(P); c3=KP3; c4=KP4; c5=KP5*sqrt(P); c6=KP6*sqrt(P); c7=KP7*(P^2); c8=KP8*P; c9=KP9/(P^2); c10=KP10; %c11=1/KP11; % CH3 c12=1/KP12; c13=KP13; c14=1/KP14; c15=1/KP15; c16=KP16; c17=1/KP17; c18=1/KP18; c19=1/KP19; d1=1/(c1^2); d2=1/(c2^2); d3=1/(c3^2); d4=1/(c4^2); d5=1/(c5^2); d6=1/(c6^2); d8=1/(c8^2); % Inicio do programa for i=1:50 F5Y7=2*d1*y7; F6Y8=2*d2*y8; F7Y4=(-(y9^2)*d3)/(y4^2); F7Y9=2*y9*d3/y4; F8Y4=(-d4*(y10)^2)/(y4^2); F8Y10=(2*y10*d4)/y4; F9Y2=2*y2*d5/(y6^2); F9Y6=-2*(y2^2)*d5/(y6^3); F10Y1=(2*y1*d6)/(y5^2); F10Y5=-2*(y1^2)*d6/(y5^3); Fn1=(y1+y5+y12+y14+y15+y21+y22); Fn2=(2*y2+2*y6+y7+y9+4*y12+3*y13+y14+2*y15+y18+2*y20+y21); Fn3=(2*y1+y2+2*y4+y5+y8+y9+y10+y15+y16+2*y17+2*y18+y21); Fn4=(2*y3+y10+y13+y14+2*y16+y17+y19+y20+y22); F12Y2=-c7*y5*(y6^3)/(y2^2); F12Y5=c7*(y6^3)/y2; F12Y6=c7*y5*3*(y6^2)/y2; F13Y13=2*d8*y13/(y6^3); F13Y6=-(d8*3*(y13^2)/(y6^4)); F14Y6=-3*c9*y13*y12/(y6^4); F14Y12=c9*y13/(y6^3); F14Y13=c9*y12/(y6^3); F15Y1=c10*y2/y4; F15Y2=c10*y1/y4; F15Y4=-c10*y1*y2/(y4^2); F16Y3=c12*y18/y9; F16Y9=-c12*y3*y18/(y9^2); F16Y18=c12*y3/y9; F17Y9=-c13*y10*y18/(y9^2); F17Y10=c13*y18/y9; F17Y18=c13*y10/y9; F18Y4=c14*y9/y8; F18Y8=-c14*y9*y4/(y8^2); F18Y9=c14*y4/y8; F19Y3=c15*y8/y10; F19Y8=c15*y3/y10; F19Y10=-c15*y3*y8/(y10^2); F20Y5=-c16*y14*y9/(y5^2); F20Y9=c16*y14/y5; F20Y14=c16*y9/y5; F21Y5=c17*y6/y7; F21Y6=c17*y5/y7; F21Y7=-c17*y5*y6/(y7^2); F22Y2=-c18*y14*y9/(y2^2); F22Y9=c18*y14/y2; F22Y14=c18*y9/y2; % Matriz Jacobiana de derivadas parciais a=[N 0 0 0 N 0 0 0 0 0 Fn1 N 0 N N 0 0 0 0 0 N N; 0 2*N 0 0 0 2*N N 0 N 0 Fn2 4*N 3*N N 2*N 0 0 N 0 2*N N 0; 2*N N 0 2*N N 0 0 N N N Fn3 0 0 0 0 N 2*N 2*N 0 0 N 1; 0 0 2*N 0 0 0 0 0 0 N Fn4 0 N N 0 2*N N 0 N N 0 N; 0 0 0 0 0 -1 F5Y7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 -1 0 0 0 F6Y8 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 F7Y4 0 -1 0 0 F7Y9 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 -1 F8Y4 0 0 0 0 0 F8Y10 0 0 0 0 0 0 0 0 0 0 0 0; 0 F9Y2 0 -1 0 F9Y6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; F10Y1 0 0 -1 F10Y5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1; 0 F12Y2 0 0 F12Y5 F12Y6 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0; 0 0 -1 0 0 F13Y6 0 0 0 0 0 0 F13Y13 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 F14Y6 0 0 0 0 0 F14Y12 F14Y13 -1 0 0 0 0 0 0 0 0; F15Y1 F15Y2 0 F15Y4 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0; 0 0 F16Y3 0 0 0 0 0 F16Y9 0 0 0 0 0 0 -1 0 F16Y18 0 0 0 0; 0 0 0 0 0 0 0 0 F17Y9 F17Y10 0 0 0 0 0 0 -1 F17Y18 0 0 0 0; 0 0 0 F18Y4 0 0 0 F18Y8 F18Y9 0 0 0 0 0 0 0 0 -1 0 0 0 0; 0 0 F19Y3 0 0 0 0 F19Y8 0 F19Y10 0 0 0 0 0 0 0 0 -1 0 0 0; 0 0 0 0 F20Y5 0 0 0 F20Y9 0 0 0 0 F20Y14 0 0 0 0 0 -1 0 0; 0 0 0 0 F21Y5 F21Y6 F21Y7 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0; 0 F22Y2 0 0 0 0 0 0 F22Y9 0 0 0 0 F22Y14 0 0 0 0 0 0 0 -1]; % Termos independentes % Sistema de eq. não lineares 10x10 % Eqs. 3.56 3.57 3.58 3.59 3.68 3.69 3.70 3.71 3.72 3.73 pag 122 e 123 % Eq. 3.56 b1=((y1+y5+y12+y14+y15+y21+y22)*N)-(e*PHI*ALFA); % Eq. 3.57 b2=((2*y2+2*y6+y7+y9+4*y12+3*y13+y14+2*y15+y18+2*y20+y21)*N)-(e*PHI*BETA); % Eq. 3.58 b3=((2*y1+y2+2*y4+y5+y8+y9+y10+y15+y16+2*y17+2*y18+y21)*N)-((e*PHI*GAMA)+0.42); % Eq. 3.59 b4=((2*y3+y10+y13+y14+2*y16+y17+y19+y20+y22)*N)-((e*PHI*DELTA)+1.58); % Eq. 3.68 b5=(d1*y7^2)-y6; % Eq. 3.69 b6=(d2*y8^2)-y4; % Eq. 3.70 b7=((d3*y9^2)/y4)-y6; % Eq. 3.71 b8=((d4*(y10^2))/y4)-y3; % Eq. 3.72 b9=((d5*y2^2)/(y6^2))-y4; % Eq. 3.73 b10=((d6*y1^2)/(y5^2))-y4; % Eq.3.60 b11=y1+y2+y3+y4+y5+y6+y7+y8+y9+y10+y12+y13+y14+y15+y16+y17+y18+y19+y20+y21+y22-1; % Equacao do metano CO + 3H2 === CH4 + H20 b12=(c7*y5*(y6^3)/y2)-y12; % Equação da amônia 1/2 N2 + 3/2 H2 === NH3 b13=(((y13^2)*d8)/y6^3)-y3; % Equação do ácido cianídrico NH3 + CH4 === HCN+ 3H2 b14=(c9*y13*y12/(y6^3))-y14; % Equação do Formol C02 + H20 ===CH20+ O2 b15=(c10*y1*y2/y4)-y15; % Equação do CH3 CH3 + NO == HCN + H2O %b16=(c11*y14*y2/y10)-y16; % Equação do N2O N2O + OH == N2 + HO2 b16=(c12*y3*y18/y9)-y16; % Equação do NO2 HO2 + NO === NO2 + OH b17=(c13*y10*y18/y9)-y17; % Equação do HO2 0 + H02 === OH + 02 (substituta) b18=(c14*y9*y4/y8)-y18; % Equação do N N + NO ==== N2 + 0 b19=(c15*y3*y8/y10)-y19; % Equação do NH2 HCN + 0H == NH2 + CO b20=(c16*y14*y9/y5)-y20; % Equacao do HCO H + HCO == H2 + CO b21=(c17*y5*y6/y7)-y21; % Equacao do CN CN + H20 == HCN + OH b22=(c18*y14*y9/y2)-y22; % Equacao NH NH + 0 === NO + H %b24=(c19*y10*y7/y8)-y24; %\Montando o vetor B com as equacoes b=[-b1;-b2;-b3;-b4;-b5;-b6;-b7;-b8;-b9;-b10;-b11;-b12;-b13;-b14;-b15;-b16;-b17;-b18;-b19;-b20;-b21;-b22]; %Resolucao de sistemas lineares pelo MATLAB ao inves de x=A\B, tem-se x=inv(a)*b y=a\b; y1=y(1)+y1; y2=y(2)+y2; y3=y(3)+y3; y4=y(4)+y4; y5=y(5)+y5; y6=y(6)+y6; y7=y(7)+y7; y8=y(8)+y8; y9=y(9)+y9; y10=y(10)+y10; N=y(11)+N; y12=y(12)+y12; y13=y(13)+y13; y14=y(14)+y14; y15=y(15)+y15; y16=y(16)+y16; y17=y(17)+y17; y18=y(18)+y18; y19=y(19)+y19; y20=y(20)+y20; y21=y(21)+y21; y22=y(22)+y22; end end for PHI=1.01:passoPHI:PHIcalculo KP1=10^((K1(1)*log(T/1000)+(K1(2)/T)+K1(3)+K1(4)*T+K1(5)*T^2)); KP2=10^((K2(1)*log(T/1000)+(K2(2)/T)+K2(3)+K2(4)*T+K2(5)*T^2)); KP3=10^((K3(1)*log(T/1000)+(K3(2)/T)+K3(3)+K3(4)*T+K3(5)*T^2)); KP4=10^((K4(1)*log(T/1000)+(K4(2)/T)+K4(3)+K4(4)*T+K4(5)*T^2)); KP5=10^((K5(1)*log(T/1000)+(K5(2)/T)+K5(3)+K5(4)*T+K5(5)*T^2)); KP6=10^((K6(1)*log(T/1000)+(K6(2)/T)+K6(3)+K6(4)*T+K6(5)*T^2)); KP7=10^((K7(1)*log(T/1000)+(K7(2)/T)+K7(3)+K7(4)*T+K7(5)*T^2)); KP8=10^((K8(1)*log(T/1000)+(K8(2)/T)+K8(3)+K8(4)*T+K8(5)*T^2)); KP9=10^((K9(1)*log(T/1000)+(K9(2)/T)+K9(3)+K9(4)*T+K9(5)*T^2)); KP10=10^((K10(1)*log(T/1000)+(K10(2)/T)+K10(3)+K10(4)*T+K10(5)*T^2)); %KP11=10^((K11(1)*log(T/1000)+(K11(2)/T)+K11(3)+K11(4)*T+K11(5)*T^2)); %Ch3 KP12=10^((K12(1)*log(T/1000)+(K12(2)/T)+K12(3)+K12(4)*T+K12(5)*T^2)); KP13=10^((K13(1)*log(T/1000)+(K13(2)/T)+K13(3)+K13(4)*T+K13(5)*T^2)); KP14=10^((K14(1)*log(T/1000)+(K14(2)/T)+K14(3)+K14(4)*T+K14(5)*T^2)); KP15=10^((K15(1)*log(T/1000)+(K15(2)/T)+K15(3)+K15(4)*T+K15(5)*T^2)); KP16=10^((K16(1)*log(T/1000)+(K16(2)/T)+K16(3)+K16(4)*T+K16(5)*T^2)); KP17=10^((K17(1)*log(T/1000)+(K17(2)/T)+K17(3)+K17(4)*T+K17(5)*T^2)); KP18=10^((K18(1)*log(T/1000)+(K18(2)/T)+K18(3)+K18(4)*T+K18(5)*T^2)); KP19=10^((K19(1)*log(T/1000)+(K19(2)/T)+K19(3)+K19(4)*T+K19(5)*T^2)); c1=KP1/sqrt(P); c2=KP2/sqrt(P); c3=KP3; c4=KP4; c5=KP5*sqrt(P); c6=KP6*sqrt(P); c7=KP7*(P^2); c8=KP8*P; c9=KP9/(P^2); c10=KP10; %c11=1/KP11; % CH3 c12=1/KP12; c13=KP13; c14=1/KP14; c15=1/KP15; c16=KP16; c17=1/KP17; c18=1/KP18; c19=1/KP19; d1=1/(c1^2); d2=1/(c2^2); d3=1/(c3^2); d4=1/(c4^2); d5=1/(c5^2); d6=1/(c6^2); d8=1/(c8^2); % Inicio do programa for i=1:50 F5Y7=2*d1*y7; F6Y8=2*d2*y8; F7Y4=(-(y9^2)*d3)/(y4^2); F7Y9=2*y9*d3/y4; F8Y4=(-d4*(y10)^2)/(y4^2); F8Y10=(2*y10*d4)/y4; F9Y2=2*y2*d5/(y6^2); F9Y6=-2*(y2^2)*d5/(y6^3); F10Y1=(2*y1*d6)/(y5^2); F10Y5=-2*(y1^2)*d6/(y5^3); Fn1=(y1+y5+y12+y14+y15+y21+y22); Fn2=(2*y2+2*y6+y7+y9+4*y12+3*y13+y14+2*y15+y18+2*y20+y21); Fn3=(2*y1+y2+2*y4+y5+y8+y9+y10+y15+y16+2*y17+2*y18+y21); Fn4=(2*y3+y10+y13+y14+2*y16+y17+y19+y20+y22); F12Y2=-c7*y5*(y6^3)/(y2^2); F12Y5=c7*(y6^3)/y2; F12Y6=c7*y5*3*(y6^2)/y2; F13Y13=2*d8*y13/(y6^3); F13Y6=-(d8*3*(y13^2)/(y6^4)); F14Y6=-3*c9*y13*y12/(y6^4); F14Y12=c9*y13/(y6^3); F14Y13=c9*y12/(y6^3); F15Y1=c10*y2/y4; F15Y2=c10*y1/y4; F15Y4=-c10*y1*y2/(y4^2); F16Y3=c12*y18/y9; F16Y9=-c12*y3*y18/(y9^2); F16Y18=c12*y3/y9; F17Y9=-c13*y10*y18/(y9^2); F17Y10=c13*y18/y9; F17Y18=c13*y10/y9; F18Y4=c14*y9/y8; F18Y8=-c14*y9*y4/(y8^2); F18Y9=c14*y4/y8; F19Y3=c15*y8/y10; F19Y8=c15*y3/y10; F19Y10=-c15*y3*y8/(y10^2); F20Y5=-c16*y14*y9/(y5^2); F20Y9=c16*y14/y5; F20Y14=c16*y9/y5; F21Y5=c17*y6/y7; F21Y6=c17*y5/y7; F21Y7=-c17*y5*y6/(y7^2); F22Y2=-c18*y14*y9/(y2^2); F22Y9=c18*y14/y2; F22Y14=c18*y9/y2; % Matriz Jacobiana de derivadas parciais a=[N 0 0 0 N 0 0 0 0 0 Fn1 N 0 N N 0 0 0 0 0 N N; 0 2*N 0 0 0 2*N N 0 N 0 Fn2 4*N 3*N N 2*N 0 0 N 0 2*N N 0; 2*N N 0 2*N N 0 0 N N N Fn3 0 0 0 0 N 2*N 2*N 0 0 N 1; 0 0 2*N 0 0 0 0 0 0 N Fn4 0 N N 0 2*N N 0 N N 0 N; 0 0 0 0 0 -1 F5Y7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 -1 0 0 0 F6Y8 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 F7Y4 0 -1 0 0 F7Y9 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 -1 F8Y4 0 0 0 0 0 F8Y10 0 0 0 0 0 0 0 0 0 0 0 0; 0 F9Y2 0 -1 0 F9Y6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; F10Y1 0 0 -1 F10Y5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1; 0 F12Y2 0 0 F12Y5 F12Y6 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0; 0 0 -1 0 0 F13Y6 0 0 0 0 0 0 F13Y13 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 F14Y6 0 0 0 0 0 F14Y12 F14Y13 -1 0 0 0 0 0 0 0 0; F15Y1 F15Y2 0 F15Y4 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0; 0 0 F16Y3 0 0 0 0 0 F16Y9 0 0 0 0 0 0 -1 0 F16Y18 0 0 0 0; 0 0 0 0 0 0 0 0 F17Y9 F17Y10 0 0 0 0 0 0 -1 F17Y18 0 0 0 0; 0 0 0 F18Y4 0 0 0 F18Y8 F18Y9 0 0 0 0 0 0 0 0 -1 0 0 0 0; 0 0 F19Y3 0 0 0 0 F19Y8 0 F19Y10 0 0 0 0 0 0 0 0 -1 0 0 0; 0 0 0 0 F20Y5 0 0 0 F20Y9 0 0 0 0 F20Y14 0 0 0 0 0 -1 0 0; 0 0 0 0 F21Y5 F21Y6 F21Y7 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0; 0 F22Y2 0 0 0 0 0 0 F22Y9 0 0 0 0 F22Y14 0 0 0 0 0 0 0 -1]; % Termos independentes % Sistema de eq. não lineares 10x10 % Eqs. 3.56 3.57 3.58 3.59 3.68 3.69 3.70 3.71 3.72 3.73 pag 122 e 123 % Eq. 3.56 b1=((y1+y5+y12+y14+y15+y21+y22)*N)-(e*PHI*ALFA); % Eq. 3.57 b2=((2*y2+2*y6+y7+y9+4*y12+3*y13+y14+2*y15+y18+2*y20+y21)*N)-(e*PHI*BETA); % Eq. 3.58 b3=((2*y1+y2+2*y4+y5+y8+y9+y10+y15+y16+2*y17+2*y18+y21)*N)-((e*PHI*GAMA)+0.42); % Eq. 3.59 b4=((2*y3+y10+y13+y14+2*y16+y17+y19+y20+y22)*N)-((e*PHI*DELTA)+1.58); % Eq. 3.68 b5=(d1*y7^2)-y6; % Eq. 3.69 b6=(d2*y8^2)-y4; % Eq. 3.70 b7=((d3*y9^2)/y4)-y6; % Eq. 3.71 b8=((d4*(y10^2))/y4)-y3; % Eq. 3.72 b9=((d5*y2^2)/(y6^2))-y4; % Eq. 3.73 b10=((d6*y1^2)/(y5^2))-y4; % Eq.3.60 b11=y1+y2+y3+y4+y5+y6+y7+y8+y9+y10+y12+y13+y14+y15+y16+y17+y18+y19+y20+y21+y22-1; % Equacao do metano CO + 3H2 === CH4 + H20 b12=(c7*y5*(y6^3)/y2)-y12; % Equação da amônia 1/2 N2 + 3/2 H2 === NH3 b13=(((y13^2)*d8)/y6^3)-y3; % Equação do ácido cianídrico NH3 + CH4 === HCN+ 3H2 b14=(c9*y13*y12/(y6^3))-y14; % Equação do Formol C02 + H20 ===CH20+ O2 b15=(c10*y1*y2/y4)-y15; % Equação do CH3 CH3 + NO == HCN + H2O %b16=(c11*y14*y2/y10)-y16; % Equação do N2O N2O + OH == N2 + HO2 b16=(c12*y3*y18/y9)-y16; % Equação do NO2 HO2 + NO === NO2 + OH b17=(c13*y10*y18/y9)-y17; % Equação do HO2 0 + H02 === OH + 02 (substituta) b18=(c14*y9*y4/y8)-y18; % Equação do N N + NO ==== N2 + 0 b19=(c15*y3*y8/y10)-y19; % Equação do NH2 HCN + 0H == NH2 + CO b20=(c16*y14*y9/y5)-y20; % Equacao do HCO H + HCO == H2 + CO b21=(c17*y5*y6/y7)-y21; % Equacao do CN CN + H20 == HCN + OH b22=(c18*y14*y9/y2)-y22; % Equacao NH NH + 0 === NO + H %b24=(c19*y10*y7/y8)-y24; %\Montando o vetor B com as equacoes b=[-b1;-b2;-b3;-b4;-b5;-b6;-b7;-b8;-b9;-b10;-b11;-b12;-b13;-b14;-b15;-b16;-b17;-b18;-b19;-b20;-b21;-b22]; %Resolucao de sistemas lineares pelo MATLAB ao inves de x=A\B, tem-se x=inv(a)*b y=a\b; y1=y(1)+y1; y2=y(2)+y2; y3=y(3)+y3; y4=y(4)+y4; y5=y(5)+y5; y6=y(6)+y6; y7=y(7)+y7; y8=y(8)+y8; y9=y(9)+y9; y10=y(10)+y10; N=y(11)+N; y12=y(12)+y12; y13=y(13)+y13; y14=y(14)+y14; y15=y(15)+y15; y16=y(16)+y16; y17=y(17)+y17; y18=y(18)+y18; y19=y(19)+y19; y20=y(20)+y20; y21=y(21)+y21; y22=y(22)+y22; end end end y1=abs(y1); y2=abs(y2); y3=abs(y3); y4=abs(y4); y5=abs(y5); y6=abs(y6); y7=abs(y7); y8=abs(y8); y9=abs(y9); y10=abs(y10); y12=abs(y12); y13=abs(y13); y14=abs(y14); y15=abs(y15); y16=abs(y16); y17=abs(y17); y18=abs(y18); y19=abs(y19); y20=abs(y20); y21=abs(y21); y22=abs(y22); N2 =1; H2O=2; CO2=3; CO=4; O2=5; OH=6; H= 7; O= 8; H2=9; NO=10; HCO=12; CH2O=13; CH4= 14; HO2= 15; NO2= 16; NH3= 17; NH2= 18; Nat = 19; HCN = 20; CN =21; N2O = 22; yEqfuel(N2)=y3; yEqfuel(H2O)=y2; yEqfuel(CO2)=y1; yEqfuel(CO)=y5; yEqfuel(O2)=y4; yEqfuel(OH)=y9; yEqfuel(H)=y7; yEqfuel(O)=y8; yEqfuel(H2)=y6; yEqfuel(NO)=y10; yEqfuel(HCO)=y21; yEqfuel(CH2O)=y15; yEqfuel(CH4)=y12; yEqfuel(HO2)=y18; yEqfuel(NO2)=y17; yEqfuel(NH3)=y13; yEqfuel(NH2)=y20; yEqfuel(Nat)=y19; yEqfuel(HCN)=y14; yEqfuel(CN)=y22; yEqfuel(N2O)=y16; end