ANEXOS. Metodo de McCabe-Thiele

ANEXOS

%%%% metodo de McCabe-Thiele %%%%

clear all  clc

%%DATOS

% declaracion de variables globales

global R W zf xw xd Lo q x y an

%Constantes de Antoine para el benceno

A1=15.9008;

B1=2788.51;

C1=52.36;

%constantes de Antoine para el etilbenceno

A2=16.0195;

B2=3279.47;

C2=59.95;

pt=1; % presion total

R=1; %relacion de reflujo;

F=620; %alimentacion (kmol/h)

D=311.765;  %destilado (kmol/h)

xd=0.95; %composicion del destilado

zf=0.5; %composicion de la alimentacion

q= 1; %factor de calor de solucion

%%CALCULOS :

% determinacion de la curva de equilibrio

Te1=B1/(A1-log(pt))-C1;

Te2=B2/(A2-log(pt))-C2;

T=linspace(Te1,Te2,10); %creamos un vector entre los valores Teb1 y Teb2

Pv1=exp(A1-B1./(T+C1));

Pv2=exp(A2-B2./(T+C2));

%composiciones de equilibrio

xe=(pt-Pv2)./(Pv1-Pv2);

ye=Pv1.*xe./pt;

%balance global de materia

W=F-D;

%determinacion de recta recta diagonal.

ejx=linspace(0,1,30);

ejy=ejx;

L=R+D; %caudal de liquido

Lo=L+q*F; %fraccion de liquido en la mezcla L-V

[an,desv]=polyfit(ye,xe,7); %ajustamos la curva de eq. a un polinomio

%plato de alimentacion

Xp=fzero('plato_aliment',0);

yc=R/(R+1)*Xp+xd/(R+1);

%seccion de enrriquecimiento

n=1; y(n)=xd; x=xd;

while x>Xp

    x(n)=feval('Curva_equil',y(n));

    n=n+1;

    y(n)=feval('Ecc_enrriq',x(n-1));

end

% seccion de agotamiento

while  x>xw

    y(n)=feval('Ecc_agot',x(n-1));

    x(n)=feval('Curva_equil',y(n));

    n=n+1;

end

%rectificacion o enrriquecedora

xre=linspace(Xp,xd,10);

yre=R/(R+1)*xre+xd/(R+1);

% agotamiento

xag=linspace(xw,Xp,10);

yag=Lo/(Lo-W).*xag-W/(Lo-W)*xw;

if q>1 %liquido debajo punto de burbuja

    %alimentacion

    xa1=linspace(zf,Xp,10);

    ya1=q/(q-1)*xa1-zf/(q-1);

elseif q==1 %liquido saturado

    q=1.0000001; %para evitar que m=infinito

    %alimentacion

    xa1=[zf Xp];

    ya1=[zf,yc];

else %vapor saturado o recalentado

    %alimentacion

    xa1=linspace(Xp,zf,10);

    ya1=q/(q-q1)*xa1-zf/(q-1);

end

h=length(x);

X=zeros(1.2*h);

Y=zeros(1.2*h);

X(1)=xd;

X(2*h)=x(h);

i=2;

j=1;

while j

    X(i:i+1)=x(j);

    i=i+2;

    j=j+1;

end

Y(1:2)=xd;

k=3;

r=2;

while r

    Y(k:k+1)=y(z);

    k=k+2;

    r=r+1;

end

plot(xe,ye,ejx,ejy)

hold on;

plot(xag,yag,'c')

plot(xre,yre,'g')

...