Diagramas T-xy mediante el calculo de punto de burbuja usando ley de Raoult modificada con MatLab

Se muestra como generar un diagrama de equilibrio liquido - vapor para mezclas binarias usando la ley de Raoult Modificada en Matlab.

% Este programa genera el diagrama T - xy usando el calculo de la temperatura
% del punto de burbuja empleando:

% - La ley de Raoult (Usando la función LeyRaoutl.m)
% - La ley de Raoutl modificada (Usando la función LeyRaoutlMod.m)

clear
clc

global P x
global A B C
global A12 A21 VL1 VL2 R
global y

% Presion y composicion liquida
P = 2; % atm
x(1) = 0.0;
x(2) = 1 - x(1);

% Numero de datos calculados
n = 20;

% Parametros de la ecuacion de Antoine
A = [8.08097, 6.87987];
B = [1582.271, 1196.76];
C = [239.726, 219.161];

% Parametros del modelo de Wilson
A12 = 1666.441; %cal/mol
A21 = 227.2126; %cal/mol
VL1 = 40.73; %cm3/mol
VL2 = 89.41; %cm3/mol
R = 1.98721; %cal/mol*K

% Temperatura de ebullicion de los componentes
Tb = [80.1, 110.6];

for i = 1:n+1
% Temperatua supuesta
T = Tb(2);

% Solucion del Punto de Burbuja
T = fsolve('LeyRaoultMod',T)

% Asignacion de vectores para graficar
x_(i) = x(1);
y_(i) = y(1);
T_(i) = T;

% Reasignaci?n de nuevas xi
x(1) = x(1) + 1.0/n;
x(2) = 1.0 - x(1);
clc

end

plot(x_,T_,y_,T_)

xlabel('x - y')
ylabel('Temperatura, C')
axis square
axis normal
axis ([0 1 0.98*min(T_) 1.02*max(T_)])
legend('Linea de punto de burbuja, T-x','Linea de punto de rocio, T-y')


function Sy = LeyRaoult(T)

global P x
global A B C
global y

for i = 1:2
% Calculo de Ps (mmHg) con Ecc. de Antoine
Ps(i) = 10^(A(i)-B(i)/(T + C(i)));
% Calculo de K
k(i) = Ps(i)/(P*760);
% Calculo de la fraccion molar de vapor
y(i) = k(i)*x(i);
end

Sy = 1.0 - sum(y);

end


function Sy = LeyRaoultMod(T)

global P x
global A B C
global A12 A21 VL1 VL2 R
global y

for i = 1:2
% Calculo de Ps (mmHg) con Ecc. de Antoine
Ps(i) = 10^(A(i)-B(i)/(T + C(i)));

% Calculo de Gamma con modelo de Wilson
A_12 = VL2/VL1*exp(-A12/(R*(T+273)));
A_21 = VL1/VL2*exp(-A21/(R*(T+273)));
gamma(1) = exp(-log(x(1)+A_12*x(2))+x(2)*((A_12/(x(1)+A_12*x(2)))-(A_21/(A_21*x(1)+x(2)))));
gamma(2) = exp(-log(x(2)+A_21*x(1))-x(1)*((A_12/(x(1)+A_12*x(2)))-(A_21/(A_21*x(1)+x(2)))));

% Calculo de K
k(i) = (Ps(i)*gamma(i))/(P*760);

% Calculo de la fracci?n molar de vapor
y(i) = k(i)*x(i);

end

Sy = 1.0 - sum(y);

end