%Binary distillation column (Luyben's column)
% This program simulates the dynamic mass balance of a binary
% distillation column.
% Assumptions:constant relative volatility
%            theoretical trays, 
%            liquid tray hydraulics
%            perfect level control in reflux drum & bottom
% 
%initialization
 trand_dat=[]; liq_o=[];liq=[];x=[];mole=[];mx=[];
 ttime=0.; errd=0;errb=0;errintd=0;errintb=0;
% Design parameters of the column;
% # of tray/# of feedstage/mole at reflux drum/mole at base/
% mole at each tray/reflux flw/vapor flw/parameters/feed flw/
  ntray=20;nfeed=10;mdo=100.;mbo=100.;mole_o=10.;rflx_o=128.01;
  vap_o=178.01;b_ta=0.1;alpha=2.0;ffeed=100.;
% initial state condition(steady state condition when feed composition =.5)
  xb=.02;xd=.98;
  x=[.035 .05719 .08885 .1318 .18622 .24951 .31618 .37948 .43391 .47688 .51526 .56295 .61896 .68052 .74345 .80319 .85603 .89995 .93458 .96079];
% PI controller setting; 
  kcd=1000.;kcb=1000.;taud=5.;taub=1.25;
% parameter setting for eular method;
  delta_t=.005;
  endtime=20; samp_t=.1;% total simulation time
% initial mole condition on each tray
 mole=mole_o*ones(1,ntray);
 mx=mole.*x;
% initial liquid flow at each tray
 liq_o=(rflx_o+ffeed)*ones(1,nfeed);
 liq_o=[liq_o rflx_o*ones(1,ntray-nfeed)];
 % simulation start
 zfeed=.55;
 trand_dat=[ttime,xb,xd,rflx_o,vap_o]
 tsamp_t=samp_t;
while ttime<endtime,
 % vapor phase concentration
 y=alpha*x./(1+(alpha-1)*x);
 liq=liq_o+(mole - mole_o)/b_ta;
  ybtm=alpha*xb/(1+(alpha-1)*xb);
% manipulated variables
% vap=vap_o; rflx=rflx_o; %non control
  errb=.02-xb;vap=vap_o-kcb*(errb+errintb/taub);
  errd=.98-xd;rflx=rflx_o+kcd*(errd+errintd/taud);
% Material balance at reflux drum
  distl=vap-rflx;
% Material balance at colum base
  btm=liq(1)-vap;
%
 if rflx<0|vap<0|distl<0|btm<0 ,
    disp('unrealistic stituation occurs');
 return;
 end;
% evaluate derivatives
xb_dot=(liq(1)*x(1)-vap*ybtm-btm*xb)/mbo;
mole_dot(1)=liq(2)-liq(1);
mx_dot(1)=vap*(ybtm-y(1))+liq(2)*x(2)-liq(1)*x(1);
% mass balance at the rectifying section
for i=2:nfeed-1,
   %mole balance at the i-th tray
   mole_dot(i)=liq(i+1)-liq(i);
   %component balance at the i-th tray
   mx_dot(i)=vap*(y(i-1)-y(i))+liq(i+1)*x(i+1)-liq(i)*x(i);
 end;
 % mass balance at the feed tray;
  mole_dot(nfeed)=liq(nfeed+1)-liq(nfeed)+ffeed;
  mx_dot(nfeed)=vap*(y(nfeed-1)-y(nfeed))+liq(nfeed+1)*x(nfeed+1)-liq(nfeed)*x(nfeed)+ffeed*zfeed;
 %mass balance at the enriching section
 for i=nfeed+1:ntray-1;
   %mole balance at the i-th tray
   mole_dot(i)=liq(i+1)-liq(i);
   %component balance at the i-th tray
   mx_dot(i)=vap*(y(i-1)-y(i))+liq(i+1)*x(i+1)-liq(i)*x(i);
 end;
 % mass balance at the reflux drum
 mole_dot(ntray)=rflx-liq(ntray);
 mx_dot(ntray)=vap*(y(ntray-1)-y(ntray))+rflx*xd-liq(ntray)*x(ntray);
 xd_dot=vap*(y(ntray)-xd)/mdo;
 
 %Euler integration
  ttime=ttime+delta_t;
   xb=xb+xb_dot*delta_t;
   mole=mole+mole_dot*delta_t;
   mx=mx+mx_dot*delta_t;
   x=mx./mole;
   if any(x<0)|any(x>1) ,
    disp('unrealistic stituation occurs');
 return;
 end;
  xd=xd+xd_dot*delta_t;
 errintd=errintd+errd*delta_t;
 errintb=errintb+errb*delta_t;
  % save the variables
  
  if ttime>=tsamp_t,
   xtemp=[ttime,xb,xd,rflx,vap]
   trand_dat=[trand_dat;xtemp];
   tsamp_t=tsamp_t+samp_t;
   end;
  end;
  % this is the end of simulation
  subplot(211);plot(trand_dat(:,1),trand_dat(:,2));
  title('concentration at the bottom');
  subplot(212);plot(trand_dat(:,1),trand_dat(:,3));
  title('concentration at the top tray');
  end;