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)

Chemical Reaction Engineering: Non Ideal Reactor

In general we study about ideal reactors like batch reactor, continuous stirred tank reactor (CSTR), plug flow reactor (PFR). Batch reactor is a unsteady state uniform reactor, PFR is a steady state non uniform reactor and CSTR is a steady state and uniform reactor. but in actual conditions the reactors deviate from their ideality. There are some reason like channeling, dead zone, bypass, back mixing etc. Due to these factors the behavior of an ideal reactor changes. The real flow pattern is found from the RTD (residence time distribution) curve. RTD shows the concentration of a tracer in output which was injected into the reactor. from the shape of RTD curve we can analyse the behavior of the reactor.

There are different models which can be used to find the non ideal behavior of the reactor like compartment model, dispersion model and Tanks in Series model. Dispersion model and Tanks in Series model shows the deviation from plug flow reactor. For an ideal plug flow behavior, the dispersion number (D/ul) will be 0 and the number of tanks in series will be infinite (∞). In PFR every cross section shows a CSTR so a PFR represent an infinite number of CSTRs connected in series.       

Programming based Projects in Chemical Engineering

Now a days computer programming is playing an important role in chemical engineering. before going for a pilot plant we study the system through simulation software like HYSYS. So learning computer programming is an important curriculum in engineering. there are a lots of projects that can be done with the help of programming in B. Tech and M. Tech. One is design of Distillation column by Shortcut method and Rigorous calculations. The further modification in the project can be done by considering the non-ideality of the solution i.e. calculating the activity coefficients.

There are many equations which can be used to find the activity coefficients, of which Margules and van Laar equations are used when the nature of the species is similar. Wilson equation, Non-random two-liquid (NRTL) equation, are based on the local composition concept instead of overall mixture composition, Universal quasi-chemical (UNIQUAC) equation and Universal functional activity coefficient (UNIFAC) method based on extension of quasi-chemical theory of liquid mixtures to solutions containing molecules of different sizes.

Thermodynamics: Residual Properties or Departure functions

Thermodynamic properties like enthalpy and entropy for any substance ideal or real, we can use Maxwell's equations, Thermodynamic Diagrams, property table like saturated steam and superheated steam  or activity coefficient. in addition to these, we can get the approximate values by using Residual Properties which are defined as the difference between the value of thermodynamics property at specified conditions and the thermodynamic property for ideal substance at same specified condition.

Chemical Reaction Engineering: Thiele modulus and its significance

Thiele modulus shows the effect of diffusion in the reaction kinetics. The Thiele modulus is given by mL = L(k/D)^0.5 The ratio of rate constant to diffusivity. The large vale of Thiele modulus represents small value of diffusivity i.e. diffusion resistance affects the overall rate of reaction. If mL<0.4, the diffusion resistance can be ignored. By calculating Thiele modulus, we can see the effect of diffusion resistance. If diffusion resistance is significant for a particular catalyst, we can decrease the resistance either by increasing the diffusivity i.e. increasing the temperature and decreasing the pressure, or by reducing the length of diffusion i.e. size the catalyst particles. Shorter the particles, lesser will be the length for diffusion and lesser will be the resistance.

Internal Energy and Enthalpy

A system have two kinds of stored energies, internal energy (U) and enthalpy (H). The internal energy is the energy required to establish the system and enthalpy is the energy required to establish the system and make a space for it (U + PV). A flowing stream has the energy equivalent to enthalpy, while a system in non flow condition has energy equivalent to internal energy.

Thermodynamics: Vapor Liquid Equilibrium

When we say the vapors are in equilibrium with the liquid, what property of the two phases will be equal? Whether it is temperature, pressure, concentration or something else that will be same? 
Equilibrium stands for no net driving force, so what is the driving force for vapor liquid system?
It is chemical potential or fugacity that is same in two phases.

Suppose we have a liquid mixture of ethyl alcohol (NBP 78.5°C) and water (NBP 100°C). At any temperature, there will be some amount of alcohol & water vapors in contact with the liquid phase. Ethyl alcohol is more volatile as compared to water, so ethyl alcohol has more concentration in the vapor phase and but how it will be calculated?

The amount of alcohol & water in vapor phase can be calculated by their vapor pressure at that temperature. The vapor pressure of alcohol is always greater than that of water. For example at 80°C the vapor pressure of alcohol is 1 atm while that of water is less, which will be 1 atm at 100°C. Therefore, at equilibrium: for ideal system i.e. ideal gas phase and ideal liquid phase (Raoult’s law):
P_Ty_A = Pv_Ax_A;  

Pv_{A is the vapor pressure of component A at specified temperature. yA is the mole fraction of component A in vapor phase and} xA _{is the mole fraction of component A} in liquid.

If partial pressure of component A (P_Ty_A) is less in vapor phase than corresponding to liquid phase (Pv_Ax_A), evaporation continues (increases y) till both are equal. and if partial pressure is greater than _{in vapor phase than corresponding to liquid phase,} condensation will start and continues (decreases y) till both are equal. 

So condensation can be made to happen by either increasing total pressure (P_T) or decreasing temperature which decreases the vapor pressure (PA)