% programa 5.7, W.L.Luyben, pag. 133
% "Multicomponent distillation dynamics" 
% modif. E.Quintero.M 5nov'13
%-------------------------------------------------------------------------
% assumptions: vapor and liquid feed, dynamics of the vapor is negligible,
%     vapor and liquid distillate, Murphree vapor-phase tray efficiency 
%     pressure is constant and varies linearly up de column,
%     liquid hydraulics from Francis formula,
%     perfect control in reflux drum and column base,
%     energy equation is algebraic,
%     reflux R ans heat input Qr are held constant (open loop response)
clear
format short e
global nc mw dens c1 c2 c3 bpt avp bvp
%---------------------------------------------------------
% column data, physical property data, feed and conditions
time=0; tprint=0; delta=0.0001; tf=0.001+eps;
nt=15; nf=5; nc=5; whs=0.75; whr=1.25; ds=72; dr=72; 
wls=48; wlr=48; mvb=10; mvd=10;
mw=[30,50,90,130,300]; dens=[40,40,60,70,90]; 
hvap=[100,90,70,80,80]; bpt=[10,90,150,210,360]; 
hcapv=[0.2,0.4,0.3,0.3,0.3]; hcapl=[0.6,0.6,0.5,0.4,0.4]; 
vp1=[14.7,14.7,14.7,14.7,14.7]; t1=[10,90,150,210,360];
vp2=[50,500,150,150,150]; t2=[30,200,200,300,420];
for j=1:nc
    avp(j)=(t1(j)+460)*(t2(j)+460)*log(vp2(j)/vp1(j))/(t1(j)-t2(j));
    bvp(j)=log(vp2(j))-avp(j)/(t2(j)+460);
    c2(j)=hcapv(j)*mw(j);
    c3(j)=hcapl(j)*mw(j);
    c1(j)=hvap(j)*mw(j)+(c3(j)-c2(j))*bpt(j);
end
fl=800; tfl=120; xf=[0.05,0.6,0.01,0.3,0.04];
fv=200; tfv=120; yf=[0.4,0.53,0.02,0.05,0];
[hlf,hvf]=enth(tfl,xf,yf);
pd=19.7; pb=21.2; qr=5; r=400; dv=200; eff=0.5;
%-------------------
% initial conditions
tb=201.58; xb=[0,.00725,.0488,.836,.108];
t=[154.9,132.6,120.2,114,108.4,101.2,98.2,96.9,96.2,95.8,...
    95.5,95.3,95.1,94.9,94.2];
lo=[740.1,814.4,892,960.1,986,320,381.9,409.6,423.7,431.2,...
    435.2,437.5,438.7,439.5,438.6];
x=[.999e-11,.11,.240,.607,.433e-1;        % tray 1
   .156e-8,.286,.202,.473,.393e-1;
   .182e-6,.457,.131,.376,.359e-1;
   .133e-4,.572,.803e-1,.314,.333e-1;
   .760e-3,.634,.496e-1,.284,.325e-1;     % tray 5
   .112e-2,.818,.866e-1,.942e-1,.174e-5;
   .129e-2,.910,.446e-1,.440e-1,.776e-6;
   .136e-2,.953,.238e-1,.218e-1,.371e-6;
   .140e-2,.975,.128e-1,.110e-1,.181e-6;
   .142e-2,.986,.694e-2,.563e-2,.893e-7;  % tray 10
   .143e-2,.992,.374e-2,.286e-2,.440e-7;
   .144e-2,.995,.199e-2,.145e-2,.216e-7;
   .144e-2,.997,.104e-2,.718e-3,.104e-7;
   .145e-2,.998,.519e-3,.342e-3,.484e-8;
   .175e-2,.998,.236e-3,.149e-3,.205e-8]; % tray 15
td=77.26; xd=[.174e-1,.982,.824e-4,.493e-4,.659e-9];
%--------------------------
% calculate initial holdups
[mwa,densa]=mwdens(xb);
mb=mvb*densa/mwa;
for n=1:nf
    for j=1:nc; xx(j)=x(n,j); end
    [mwa,densa]=mwdens(xx);
    lv(n)=lo(n)*mwa/densa;
    l(n)=lo(n);
    hfow=(lv(n)/(999*wls))^.66667;
    mv(n)=(hfow+whs/12)*3.1416*ds*ds/(4*144);
    m(n)=mv(n)*densa/mwa;
end
for n=nf+1:nt
    for j=1:nc; xx(j)=x(n,j); end
    [mwa,densa]=mwdens(xx);
    lv(n)=lo(n)*mwa/densa;
    l(n)=lo(n);
    hfow=(lv(n)/(999*wlr))^.66667;
    mv(n)=(hfow+whr/12)*3.1416*dr*dr/(4*144);
    m(n)=mv(n)*densa/mwa;
end
for n=1:nt
    for j=1:nc
        xm(n,j)=m(n)*x(n,j);
    end
end
[mwa,densa]=mwdens(xd);
md=mvd*densa/mwa;
%---------------------------
% calculate pressure profile
for n=1:nt
    p(n)=(pb-(n*(pb-pd))/nt);
end
v(5)=822; % initial guess of v(5) for first efficiency calc.
%-----------------------------
% main loop for each time step
while(time<=tf)
    [tb,yb]=bubpt(tb,xb,pb);
    [hlb,hvb]=enth(tb,xb,yb);
    for j=1:nc; xx(j)=x(1,j); end
    [t(1),yy]=bubpt(t(1),xx,p(1));
    for j=1:nc
        y(1,j)=yb(j)+eff*(yy(j)-yb(j));
        yy(j)=y(1,j);
    end
    [hl(1),hv(1)]=enth(t(1),xx,yy);
    for n=2:nf  % D0 110 WLL
        for j=1:nc; xx(j)=x(n,j); end
        [t(n),yy]=bubpt(t(n),xx,p(n));
        for j=1:nc
            y(n,j)=(yy(j)-y(n-1,j))*eff+y(n-1,j);
            yy(j)=y(n,j);
        end
        [hl(n),hv(n)]=enth(t(n),xx,yy);
    end
    for j=1:nc; xx(j)=x(nf+1,j); end
    [t(nf+1),yy]=bubpt(t(nf+1),xx,p(nf+1));
    for j=1:nc
        yav(j)=(yf(j)*fv+y(nf,j)*v(nf))/(v(nf)+fv);
        y(nf+1,j)=(yy(j)-yav(j))*eff+yav(j);
        yy(j)=y(nf+1,j);
    end
    [hl(nf+1),hv(nf+1)]=enth(t(nf+1),xx,yy);
    for n=nf+2:nt   % D0 115 WLL
        for j=1:nc; xx(j)=x(n,j); end
        [t(n),yy]=bubpt(t(n),xx,p(n));
        for j=1:nc
            y(n,j)=(yy(j)-y(n-1,j))*eff+y(n-1,j);
            yy(j)=y(n,j);
        end
        [hl(n),hv(n)]=enth(t(n),xx,yy);
    end
    [td,yd]=bubpt(td,xd,pd);
    [hld,hvd]=enth(td,xd,yd);
    %----------------------
    % calculate vapor rates
    vb=(qr*1000000-l(1)*(hlb-hl(1)))/(hvb-hlb);
    b=l(1)-vb;
    if (b<0); error('reboiler empty'); end
    v(1)=(hl(2)*l(2)+hvb*vb-hl(1)*l(1))/hv(1);
    for n=2:nf-1
        v(n)=(hl(n+1)*l(n+1)+hv(n-1)*v(n-1)-hl(n)*l(n))/hv(n);
    end
    v(nf)=(hl(nf+1)*l(nf+1)+hv(nf-1)*v(nf-1)-hl(nf)*l(nf)+hlf*fl)/hv(nf);
    v(nf+1)=(hl(nf+2)*l(nf+2)+hv(nf)*v(nf)-hl(nf+1)*l(nf+1)+hvf*fv)/hv(nf+1);    
    for n=nf+2:nt-1
        v(n)=(hl(n+1)*l(n+1)+hv(n-1)*v(n-1)-hl(n)*l(n))/hv(n);
    end
    v(nt)=(hld*r+hv(nt-1)*v(nt-1)-hl(nt)*l(nt))/hv(nt);
    dl=v(nt)-dv-r;
    %--------------------
    %evaluate derivatives
    dm(1)=l(2)+vb-v(1)-l(1);
    for n=2:nf-1; dm(n)=l(n+1)+v(n-1)-l(n)-v(n); end
    dm(nf)=l(nf+1)+fl+v(nf-1)-l(nf)-v(nf);
    dm(nf+1)=l(nf+2)+fv+v(nf)-l(nf+1)-v(nf+1);
    for n=nf+2:nt-1; dm(n)=l(n+1)+v(n-1)-l(n)-v(n); end
    dm(nt)=r+v(nt-1)-l(nt)-v(nt);
    for j=1:nc % DO 160 WLL
        dxb(j)=(x(1,j)*l(1)-yb(j)*vb-xb(j)*b)/mb;
        dxm(1,j)=x(2,j)*l(2)+yb(j)*vb-x(1,j)*l(1)-y(1,j)*v(1);
        for n=2:nf-1
            dxm(n,j)=x(n+1,j)*l(n+1)+y(n-1,j)*v(n-1)-x(n,j)*l(n)-y(n,j)*v(n);
        end
        dxm(nf,j)=x(nf+1,j)*l(nf+1)+y(nf-1,j)*v(nf-1)-x(nf,j)*l(nf)-y(nf,j)*v(nf)+fl*xf(j);
        dxm(nf+1,j)=x(nf+2,j)*l(nf+2)+y(nf,j)*v(nf)-x(nf+1,j)*l(nf+1)-y(nf+1,j)*v(nf+1)+fv*yf(j);
        for n=nf+2:nt-1
            dxm(n,j)=x(n+1,j)*l(n+1)+y(n-1,j)*v(n-1)-x(n,j)*l(n)-y(n,j)*v(n);
        end
        dxm(nt,j)=xd(j)*r+y(nt-1,j)*v(nt-1)-x(nt,j)*l(nt)-y(nt,j)*v(nt);
        dxd(j)=(v(nt)*y(nt,j)-dv*yd(j)-(r+dl)*xd(j))/md;
    end
    if(time>=tprint)
      xx11=[time,tb,xb(1:nc),b];
      xx12=[(1:nt)',t(1:nt)',x,l(1:nt)'];
      xx13=[nt+1,td,xd(1:nc),r];
      xx14=[0,0,yd(1:nc),dl];
      resultados=[xx11;xx12;xx13;xx14];
      tprint=tprint+0.0005;
    end
    %------------------
    % euler integration
    time=time+delta;
    for n=1:nt; m(n)=m(n)+dm(n)*delta; end
    for j=1:nc   % DO 220 WLL
        xb(j)=xb(j)+dxb(j)*delta;
        if (xb(j)<0); xb(j)=0; end
        if (xb(j)>1); xb(j)=1; end
        for n=1:nt
            xm(n,j)=xm(n,j)+dxm(n,j)*delta;
            x(n,j)=xm(n,j)/m(n);
            if (x(n,j)<0); x(n,j)=0; end
            if (x(n,j)>1); x(n,j)=1; end
        end
        xd(j)=xd(j)+dxd(j)*delta;
        if (xd(j)<0); xd(j)=0; end
        if (xd(j)>1); xd(j)=1; end
    end
    %----------------------
    % calc new liquid rates
    for n=1:nf
        for j=1:nc xx(j)=x(n,j); end
        [l(n)]=hydrau(m(n),xx,whs,wls,ds);
    end
    for n=nf+1:nt
        for j=1:nc xx(j)=x(n,j); end
        [l(n)]=hydrau(m(n),xx,whr,wlr,dr);
    end
end
