RE:Proyecto fin de carrera
Pues resulta que tengo una solución de las funciones pero me dan error. Los resultados (x1 x2 x3 x4 x5 x6 x7 x8 x9 x10) de esta que envío como ejemplo deberían salir todos números positivos y entre 0 y 1. Desconozco la razón por la que no es así. Tengo unos documentos en word que explican un poco lo que está implementado en los códigos, si los necesitas, te los paso por correo electrónico.
%PROGRAMA DE CALCULO DE COMPOSICIONES DE MEZCLAS LÍQUIDO-VAPOR%
function ractor3
%VARIABLES DE OPERACION A INTRODUCIR
% Temperatura (K)
global t p a1 b c d a
t1=190;
t=t1+273;
%Presion (MPa)
p=0.45;
%Fracciones de las especies (H2O NH3 CO2 UREA)
a1=0.087109359;
b=0.81960057;
c=0.09329006;
d=0.00000000;
%COMANDOS PARA LA EJECUCION DEL PROGRAMA
iter=1000;
tol=0.001.*[1 1 1];
m=1;
n=1;
lnp1=-5231.82/t-0.06167*log(t)-0.003291*t+0.000001222*t^2+13.183;
lnp2=-25.070/t+56.321*log(t)-0.2625*t+0.0001753*t^-2-258.139;
lnp3=-2370.26/t-0.591*log(t)-0.01178*t+0.00001598*t^2+15.272;
f=exp([lnp1 lnp2 lnp3]);
K=f/p;
disp(t)
ac=[1 1 1 1 1 1 1];
while(m<=iter)
x1=fsolve(@ztres,[0.2479 0.0013 0.0098 0.0015 0.0089 0.0061 0.4818 0.237 0.0036 0.0055]);
x2=[x1(1) x1(2) x1(3) x1(4) x1(5) x1(6) x1(7)]./sum([x1(1) x1(2) x1(3) x1(4) x1(5) x1(6) x1(7)]);
y2=[x1(8) x1(9) x1(10)]./sum([x1(8) x1(9) x1(10)]);
ac=actr(x2,t);
ac2=[ac(1) ac(2) ac(3)];
fugac=fuga(t,y2);
K1=(f.*ac2)./(p.*fugac);
if(abs(K1-K)<=tol)
x2;
K;
break;
else
K=K1;
n=n+1;
end
end
if(i==iter)
disp('no converge');
end
%VARIABLES RESUELTAS POR EL PROGRAMA
nl=[x1(1) x1(2) x1(3) x1(4) x1(5) x1(6) x1(7) x1(8) x1(9) x1(10)];
disp('cantidades molares de las diez especies=')
disp(nl);
nlt=sum(nl);
L=nlt-x2(4);
disp('fracciones molares de las especies líquidas');
disp(L);
nv=[y2(1) y2(2) y2(3)];
nlv=sum(nv);
V=nlv;
disp('fracciones molares de las especies vapor-gas');
disp(V);
%SUBPROGRAMA DE CALCULO
function m=ztres(x);
global t p a1 b c d a
lnp1=-5231.82/t-0.06167*log(t)-0.003291*t+0.000001222*t^2+13.183;
lnp2=-25.070/t+56.321*log(t)-0.2625*t+0.0001753*t^-2-258.139;
lnp3=-2370.26/t-0.591*log(t)-0.01178*t+0.00001598*t^2+15.272;
f=exp([lnp1 lnp2 lnp3]);
x1=x(1);x2=x(2);x3=x(3);x4=x(4);x5=x(5);x6=x(6);x7=x(7);
y1=x(8);y2=x(9);y3=x(10);
m=zeros(12,1);
L=x1+x2+x3+x4+x5+x6+x7;
V=y1+y2+y3;
xi=[x1/L x2/L x3/L x4/L x5/L x6/L x7/L];
yi=[y1/V y2/V y3/V];
ac=actr(t,xi);
activ=[ac(1) ac(2) ac(3)];
fugac=fuga(t,yi);
K=f.*activ./(p.*fugac);
%Constante de equilibrio de la reaccion 7
K7=exp(9.9068e3./t+7.4296e-2*log(t)-5.3985e-3*t-20.2220);
%Constante de equilibrio de la reaccion 8
K8=exp(8.8226e3./t+0.8404e-2*log(t)+1.8736e-3*t-21.6135);
%BALANCES
%Equilibrios químicos
m(1)=K7-(((ac(4)*ac(6)./(ac(3)*ac(2).^2)).^2)*x4*x6*L./(x3*x2.^2));
m(2)=K8-(((ac(4)*ac(5)./(ac(1)*ac(2)*ac(3))).^2)*x4*x5*L./(x1*x2*x3));
%Balance al agua
m(3)=x1+x5+y1-a1;
%Balance al amoniaco
m(4)=x2+x4+x6+y2-b;
%Balance al dióxido de carbono
m(5)=x3+x5+x6+y3-c;
%Balance a la urea
m(6)=x7-d;
%Balance de cargas
m(7)=x4-x5-x6;
%Equilibrios líquido-vapor
m(8)=(x1/L)*K(1)-(y1/V);
m(9)=(x2/L)*K(2)-(y2/V);
m(10)=(x3/L)*K(3)-(y3/V);
m(11)=(x1+x2+x3+x4+x5+x6+x7)./L-1;
m(12)=(y1+y2+y3)./V-1;
m=real([m(1) m(2) m(3) m(4) m(5) m(6) m(7) m(8) m(9) m(10) m(11) m(12)]);