Program to design distillation column by FUG (Short cut) Method

Program to design distillation column by FUG (Short cut) Method
there are two methods to analyze a distillation column. one is graphical ( McCabe–Thiele method & Ponchon-Savarit ) and other is analytical (FUG: Fenske-Underwood-Gilliland). In FUG method minimum trays are calculated by Fenske equation, minimum reflux by undrewood equation and theoretical stages by Gilliland equation.

Take composition of pentane, hexane, heptane, octane as an input from the user:
    xipent=input('enter the composition of pentane');
    xihex=input('enter the composition of hexane');
    xihept=input('enter the composition of heptane');
    xioct=input('enter the composition of octane');
    pt=input('enter the total pressure');

Program to calculate bubble point and dew point at the Feed composition
To calculate bubble point and dew point we require the values of saturation pressure and distributive coefficient K for each component at different temperatures. So we  have prepared a separate function called psat() for the claculation of pressure in which the value of assumed temperature is passed as an argument. This has also served to reduce the program length because it is required to calculate pressure in every iteration untill the condition of ∑Ki*xi=1.0 is satisfied. The definition of the function is:
            function[prpent,prhex,prhept,proct]=psat(t)
            penta=13.7667;pentb=2451.88;pentc=-41.136;
            hexa=13.8193;hexb=2696.04;hexc=-48.833;
            hepta=13.8622;heptb=2910.26;heptc=-56.718;
          octa=13.9346;octb=3123.13;octc=-63.515;
            prpent=exp(penta-((pentb/(t+pentc))));
            prhex=exp(hexa-((hexb/(t+hexc))));
            prhept=exp(hepta-((heptb/(t+heptc))));
            proct=exp(octa-((octb/(t+octc))));
         end

Now the main body of the programme for calculation of bubble point where the function psat() has been called is:
disp('now firstly bubble point and dew point claculations will be done')
    for t=323:0.1:423
            [p,q,r,s]=psat(t);
            kipent1=p/pt;
            kihex1=q/pt;
            kihept1=r/pt;
            kioct1=s/pt;
            sumkixi=(kipent1*xipent)+(kihex1*xihex)+(kihept1*xihept)+(kioct1*xioct);
            if(sumkixi>0.99&sumkixi<1)
                disp('the bubble point of the mixture is')
                disp(t)
                break
            end
    end
    disp('the value of equilibrium constant of light key at bubble point is:')
    disp(kihex1)
    disp('the value of equilibrium constant of heavy key at bubble point is:')
    disp(kihept1)
    for temp=343:0.1:573
        [a,b,c,d]=psat(temp);
            kipent=a/pt;
            kihex=b/pt;
            kihept=c/pt;
            kioct=d/pt;
            sumy=(xipent/kipent)+(xihex/kihex)+(xihept/kihept)+(xioct/kioct);
        if(sumy>0.99&sumy<1)
            disp('the dew point of the mixture is')
            disp(temp)
            break
        end
    end

program to calculate minimum number of trays by Fenske equation
alpha=kihex1/kihept1;
disp('the value of relative volatility is')
disp(alpha)
disp('the composition of light key that is of hexane in distillate is:')
mlkd=(0.98*xihex);
disp(mlkd)
disp('the composition of heavy key that is of heptane in distillate is:')
mhkd=(0.01*xihept);
disp(mhkd)
disp('the composition of light key that is of hexane in bottom is:')
mlkb=(xihex-mlkd);
disp(mlkb)
disp('the composition of heavy key that is of heptane in bottom is:')
mhkb=(xihept-mhkd);
disp(mhkb)
disp('now the minimum number of stages by Fenske equation is calculated')
disp('the fenske equation is:')
disp('Nmin=ln[(xdi/xbi)/(xdj/xbj)]/((ln alphaij)-1)')
disp('where, Nmin= minimum number of trays')
disp('alpha= relative volatility')
nmin=(log(((mlkd/mhkd)/(mlkb/mhkb)))/log(alpha))-1;
disp('the minimum number of trays required are')
disp(nmin)
disp('minimum number of stages including reboiler are')
nminn=nmin+1;

program to calculate minimum reflux ratio is:
disp('now the calculation of minimum reflux by undrewood method is done:')
disp('Rdm+1= summation[(alphai*xdi)/(alphai-phi)]')
disp('where phi= root of the equation whose value lie between the relative volativity of heavy key and light key')
disp('Rdm=minimum reflux ratio')
disp('alphai=relative volativity of respective components')
disp('xdi=compositions in distillate of respective components')
alpha1=kipent1/kihept1;
alpha2=kihex1/kihept1;
alpha3=kihept1/kihept1;
alpha4=kioct1/kihept1;
for phi=alpha3:0.1:alpha2
  sumf=((alpha1*xipent)/(alpha1-phi))+((alpha2*xihex)/(alpha2-    phi))+((alpha3*xihept)/(alpha3-phi))+((alpha4*xioct)/(alpha4-phi));
if(sumf>-0.9999&sumf<1)
 disp('the value of phi is')
       disp(phi)
       break;
   end;
end;
disp('we need 98% recovery of hexane and 1% heptane in distillate')
 disp('the moles of all the four components in distillate are calculated below by simple material balnce calculations:')
mpent=xipent*100;
mhex=mlkd*100;
    mhept=mhkd*100;
moct=0;
disp('the moles of pentane in distillate is:');
disp(mpent)
disp('the moles of hexane in distillate is:');
disp(mhex)
disp('the moles of heptane in distillate is:');
disp(mhept)
disp('the moles of octane in distillate is:');
disp(moct)
summoles=mpent+mhex+mhept+moct;
disp('the total moles in distillate are')
disp(summoles)
sumrd=(((alpha1*(mpent/summoles))/(alpha1-phi))+((alpha2*(mhex/summoles))/(alpha2-phi))+((alpha3*(mhept/summoles))/(alpha3-phi))+((alpha4*(mpent/summoles))/(alpha4-phi)))-1;
    disp('minimum reflux ratio is')
disp(sumrd)

5. program to calculate theoretical stages by Gilliland method is:

disp('now the calculations will be done on the basis of Gilliland corelations for claculating actual number of stages')
disp('actual reflux ration is taken as 1.5 times the minimum reflux ratio')
rd=1.5*sumrd;
x=(rd-sumrd)/(rd+1);
y=1-exp(((1+(54.4*x))/(11+(117.2*x)))*((x-1)/(x.^0.5)));
n=(y+nminn)/(1-y);
disp('required number of stages are')
disp(n)

Coded by: Ravisha Goswami and group
(B.Tech 2010-2014, Chemical Engg., BTKIT Dwarahat)