“UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ”
“FACULTAD DE INGENIERÍA QUÍMICA”
UNIVERSIDAD CIENTIFICA, INNOVADORA UN INGENIERO QUIMICO, UNA
INVESTIGADORA Y FORMADORA. EMPRESA.
PROGRAMA EVL (ISOBÁRICO) MÉTODO
UNIFAC
CÁTEDRA:
TERMODINÁMICA DE LOS PROCESOS QUÍMICOS II
CATEDRÁTICO:
Dr. ORÉ VIDALÓN Salvador
ALUMNO:
QUISPE INGA Joel Alexander
Huancayo, Perú
EVL EN UN PROCESO ISOBÁRICO MODELO UNIFAC
%PRESENTADO POR: QUISPE INGA JOEL ALEXANDER
%E.A.P: INGENIERÍA QUÍMICA
%FECHA: 9-06-2019
clear
clc
disp('Coeficientes de Actividad para el sistema acetone(1)-
water(2):MODELO UNIFAC')
disp('---------------------------------------------------------------')
%PRESIÓN A LA QUE EL SISTEMA ES ISOBÁRICO:
%Presión en mm-Hg:
P=760;
%CONSTANTES DE ANTOINE, (TABLA I):
%Acetone(1):
A1=7.02447;
B1=1161.0;
C1=224.0;
%Water(2):
A2=7.96681;
B2=1668.21;
C2=228.0;
%Valor de z:
z=10;
%FRACCIONES MOLARES EN LA FASE LIQUIDA DE ACETONE(1) Y WATER(2):
x1=0.1:0.05:0.9;
x2=1-x1; n=length(x1)
%Declarar a la temperatura una variable
syms T
%ECUACIÓN DE ANTOINE PARA HALLA LA Pvp:
Pvp1=10^(A1-B1/(C1+T-273));
Pvp2=10^(A2-B2/(C2+T-273));
%CÁLCULO DE LAS TEMPERATURAS INICIALES:
%Calculo de las temperaturas iniciales
T1=B1/(A1-log10(P))-C1+273;
T2=B2/(A2-log10(P))-C2+273;
%DATOS GRUPOS FUNCIONALES:
%VOLUME (R):
R1=0.9011;
R2=1.6724;
R3=0.9200;
%SURFACE AREA (Q):
Q1=0.848;
Q2=1.488;
Q3=1.400;
%vk1(acetona):
v11=1;
v21=1;
v31=0;
%vk2(metanol):
v12=0;
v22=0;
v32=1;
%PARÁMETRO DE INTERACCIÓN BINARIA ajk:
%aj1:
a11=0;
a21=26.76;
a31=300.0;
%aj2:
a12=476.42;
a22=0;
a32=-195.4;
%aj3:
a13=1318;
a23=472.5;
a33=0;
%PARTE COMBINATORIAL:
%ri:
r1=v11*R1+v21*R2+v31*R3;
r2=v12*R1+v22*R2+v32*R3;
%qi:
q1=v11*Q1+v21*Q2+v31*Q3;
q2=v12*Q1+v22*Q2+v32*Q3;
%PARTE RESIDUAL:
%Hallando Xmi:
X11=v11/(v11+v21+v31);
X21=v21/(v11+v21+v31);
X31=v31/(v11+v21+v31);
X12=v12/(v12+v22+v32);
X22=v22/(v12+v22+v32);
X32=v32/(v12+v22+v32);
%Hallando theta(mi):
theta11=X11*Q1/(X11*Q1+X21*Q2+X31*Q3);
theta21=X21*Q2/(X11*Q1+X21*Q2+X31*Q3);
theta31=X31*Q3/(X11*Q1+X21*Q2+X31*Q3);
theta12=X12*Q1/(X12*Q1+X22*Q2+X32*Q3);
theta22=X22*Q2/(X12*Q1+X22*Q2+X32*Q3);
theta32=X32*Q3/(X12*Q1+X22*Q2+X32*Q3);
%Hallando psi(mn):
psi11=exp(-a11/T);
psi21=exp(-a21/T);
psi31=exp(-a31/T);
psi12=exp(-a12/T);
psi22=exp(-a22/T);
psi32=exp(-a32/T);
psi13=exp(-a13/T);
psi23=exp(-a23/T);
psi33=exp(-a33/T);
%Cálculo del coeficiente de actividad rki:
lnr11=Q1*[1-log(theta11*psi11+theta21*psi21+theta31*psi31)-
[(theta11*psi11/(theta11*psi11+theta21*psi21+theta31*psi31))+(theta21*psi
12/(theta11*psi12+theta21*psi22+theta31*psi32))+(theta31*psi13/(theta11*p
si13+theta21*psi23+theta31*psi33))]];
lnr21=Q2*[1-log(theta11*psi12+theta21*psi22+theta31*psi32)-
[(theta11*psi21/(theta11*psi11+theta21*psi21+theta31*psi31))+(theta21*psi
22/(theta11*psi12+theta21*psi22+theta31*psi32))+(theta31*psi23/(theta11*p
si13+theta21*psi23+theta31*psi33))]];
lnr31=Q3*[1-log(theta11*psi13+theta21*psi23+theta31*psi33)-
[(theta11*psi31/(theta11*psi11+theta21*psi21+theta31*psi31))+(theta21*psi
32/(theta11*psi12+theta21*psi22+theta31*psi32))+(theta31*psi33/(theta11*p
si13+theta21*psi23+theta31*psi33))]];
lnr12=Q1*[1-log(theta12*psi11+theta22*psi21+theta32*psi31)-
[(theta12*psi11/(theta12*psi11+theta22*psi21+theta32*psi31))+(theta22*psi
12/(theta12*psi12+theta22*psi22+theta32*psi32))+(theta32*psi13/(theta12*p
si13+theta22*psi23+theta32*psi33))]];
lnr22=Q2*[1-log(theta12*psi12+theta22*psi22+theta32*psi32)-
[(theta12*psi21/(theta12*psi11+theta22*psi21+theta32*psi31))+(theta22*psi
22/(theta12*psi12+theta22*psi22+theta32*psi32))+(theta32*psi23/(theta12*p
si13+theta22*psi23+theta32*psi33))]];
lnr32=Q3*[1-log(theta12*psi13+theta22*psi23+theta32*psi33)-
[(theta12*psi31/(theta12*psi11+theta22*psi21+theta32*psi31))+(theta22*psi
32/(theta12*psi12+theta22*psi22+theta32*psi32))+(theta32*psi33/(theta12*p
si13+theta22*psi23+theta32*psi33))]];
%Hallando Li:
L1=(z/2)*(r1-q1)-(r1-1);
L2=(z/2)*(r2-q2)-(r2-1);
for i=1:n
%fii:
fi1(i)=(x1(i)*r1)/(x1(i)*r1+x2(i)*r2);
fi2(i)=(x2(i)*r2)/(x1(i)*r1+x2(i)*r2);
%thetai:
thetaA(i)=(x1(i)*q1)/(x1(i)*q1+x2(i)*q2);
thetaB(i)=(x2(i)*q2)/(x1(i)*q1+x2(i)*q2);
%Coeficiente de actividad que queremos hallar:
lnyC1(i)=log(fi1(i)/x1(i))+(z/2)*q1*log(thetaA(i)/fi1(i))+L1-
(fi1(i)/x1(i))*(x1(i)*L1+x2(i)*L2);
lnyC2(i)=log(fi2(i)/x2(i))+(z/2)*q2*log(thetaB(i)/fi2(i))+L2-
(fi2(i)/x2(i))*(x1(i)*L1+x2(i)*L2);
%CÁLCULO DEL COEFICIENTE DE ACTIVIDAD rk:
%Hallando Xm:
X1(i)=(x1(i)*v11+x2(i)*v12)/(x1(i)*v11+x2(i)*v12+x1(i)*v21+x2(i)*v22+x1(i
)*v31+x2(i)*v32);
X2(i)=(x1(i)*v21+x2(i)*v22)/(x1(i)*v11+x2(i)*v12+x1(i)*v21+x2(i)*v22+x1(i
)*v31+x2(i)*v32);
X3(i)=(x1(i)*v31+x2(i)*v32)/(x1(i)*v11+x2(i)*v12+x1(i)*v21+x2(i)*v22+x1(i
)*v31+x2(i)*v32);
%Hallando theta(m):
theta1(i)=X1(i)*Q1/(X1(i)*Q1+X2(i)*Q2+X3(i)*Q3);
theta2(i)=X2(i)*Q2/(X1(i)*Q1+X2(i)*Q2+X3(i)*Q3);
theta3(i)=X3(i)*Q3/(X1(i)*Q1+X2(i)*Q2+X3(i)*Q3);
%hallando lnrk:
lnr1(i)=Q1*[1-log(theta1(i)*psi11+theta2(i)*psi21+theta3(i)*psi31)-
[(theta1(i)*psi11/(theta1(i)*psi11+theta2(i)*psi21+theta3(i)*psi31))+(the
ta2(i)*psi12/(theta1(i)*psi12+theta2(i)*psi22+theta3(i)*psi32))+(theta3(i
)*psi13/(theta1(i)*psi13+theta2(i)*psi23+theta3(i)*psi33))]];
lnr2(i)=Q2*[1-log(theta1(i)*psi12+theta2(i)*psi22+theta3(i)*psi32)-
[(theta1(i)*psi21/(theta1(i)*psi11+theta2(i)*psi21+theta3(i)*psi31))+(the
ta2(i)*psi22/(theta1(i)*psi12+theta2(i)*psi22+theta3(i)*psi32))+(theta3(i
)*psi23/(theta1(i)*psi13+theta2(i)*psi23+theta3(i)*psi33))]];
lnr3(i)=Q3*[1-log(theta1(i)*psi13+theta2(i)*psi23+theta3(i)*psi33)-
[(theta1(i)*psi31/(theta1(i)*psi11+theta2(i)*psi21+theta3(i)*psi31))+(the
ta2(i)*psi32/(theta1(i)*psi12+theta2(i)*psi22+theta3(i)*psi32))+(theta3(i
)*psi33/(theta1(i)*psi13+theta2(i)*psi23+theta3(i)*psi33))]];
%Hallando lnyR:
lnyR1(i)=v11*(lnr1(i)-lnr11)+v21*(lnr2(i)-lnr21)+v31*(lnr3(i)-lnr31);
lnyR2(i)=v12*(lnr1(i)-lnr12)+v22*(lnr2(i)-lnr22)+v32*(lnr3(i)-lnr32);
%Hallando lnyi:
lny1(i)=lnyC1(i)+lnyR1(i);
lny2(i)=lnyC2(i)+lnyR2(i);
%COEFICIENTES DE ACTIVIDAD:
gamma1(i)=exp(lny1(i));
gamma2(i)=exp(lny2(i));
%HALLANDO y1 y y2:
y1(i)=x1(i)*gamma1(i)*Pvp1/P;
y2(i)=x2(i)*gamma2(i)*Pvp2/P;
end
%DETERMINACIÓN DE LAS TEMPERATURAS DE EQUILIBRIO:
disp(' T ºC x1 x2 y1 Gamma 1 Gamma 2');
disp('-------------------------------------------------------------------
-----')
f=y1+y2-1;
for i = 1:17
xo=(T1+T2)/2;%valor inicial para la temperatura
e=10^-4;%criterio de convergencia
distancia=1;
while distancia > e
fxi=subs(f(i),sym('T'),xo);
derivada=diff(f(i),sym('T'));
Dxfxi=subs(derivada,sym('T'),xo);
x_1=single(xo-(fxi/Dxfxi));
distancia=single(abs(subs(f(i),sym('T'),x_1)));
xo=x_1;
end
Tm(i)=x_1;
t(i)=Tm(i);
Y1(i)=single(subs(y1(i),sym('T'),Tm(i)));
G1(i)=single(subs(gamma1(i),sym('T'),Tm(i)));
G2(i)=single(subs(gamma2(i),sym('T'),Tm(i)));
disp([t(i) x1(i) x2(i) Y1(i) G1(i) G2(i) ])
end
disp(' Pvp1 Pvp2')
disp('-----------------------');
for i=1:17
pvp1(i)=10^(A1-B1/(C1+t(i)));
pvp2(i)=10^(A2-B2/(C2+t(i)));
fprintf('%10.4f %12.4f\n',pvp1(i)',pvp2(i)');
end
%GRÁFICA 1: DIAGRAMA x1 vs y1:
figure
plot(x1,x1,'-k',x1,Y1,'-b');
legend('x1,x1','x1,y1')
grid on
title('Equilibrio Acetone(1)-Water(2) a 760 mm-Hg; Modelo UNIFAC')
xlabel('x1 (Acetone líquido)')
ylabel('y1 (Acetone vapor)')
%GRÁFICA 2: DIAGRAMA x1,y1 vs T:
figure, plot(x1,t,'-k',Y1,t,'-b');
legend('x1,T','y1,T');
grid on
title('Equilibrio Acetone(1)-Water(2) a 760 mm-Hg; Modelo UNIFAC');
xlabel(' x1,y1');
ylabel(' T °C ');
grid on
