! Design and operation of a distillation column(distcolm.lng).
Key words: Distillation column, Petroleum Refining, Refining,
Chemical Equilibrium, Oil Refining, Process Industry;
!A distillation column is a vertical array of N trays in
an enclosed tube. A mixture of 2 or more components,
e.g., crude oil, is fed in at some tray i,  1 < i < N. 
Liquid, consisting of mainly the heavier, less volatile
components is removed from the bottom.
Gas, consisting of the lighter components is
removed from the top. Heat is applied at the bottom 
tray.  Heat is removed, e.g., via a heat exchanger, at
the top. Thus, temperature and pressure decrease from
bottom to top. At each tray, 3 things happen: 
a) per molecule, the lighter components in liquid form
 in the tray tend to evaporate more rapidly, and move
 to the tray above,
b) per molecule, heavier components in gaseous form
 tend to condense onto the tray more rapidly, and
c) liquid from the tray overflows onto the next
  lower tray.
The simplest design problem is to decide which tray is
the input tray, given the number and location of the trays.
The operational problem is to decide upon the rate of:
a) adding input to the feed tray, b) adding heat at the bottom,
c) removing heat at the top, d) removing liquid from the bottom,
   and e) removing gas at the top, so as to achieve a desired
 split of removing almost pure heavy components from the bottom
 and removing almost pure light components from the top.
  In a fractional distillation column, there are outlets at
 each tray to drain the liquid from each tray, e.g., gasoline
 towards the top, kerosene below it, diesel fuel below it, etc.;

!   Reference: L.T. Biegler, I.E. Grossmann and W. Westerberg,
               Systematic Methods of Chemical Process Design,
               Prentice-Hall, New Jersey, 1997;

! Phase equilibrium system of two components benzene-toluene,
  assuming thermodynamic model for both liquid and vapor as ideal,
  with the source of thermodynamic data from 
  R.C. Reid, J.M. Prausnitz and B.E. Poling,
  The properties of Gases and Liquids, McGraw-Hill, New York, 1987.
  Select a condenser type that totally reduces all vapor into liquid;

  Component: ! Each component of the input has..;
    tb, !normal boiling point; 
    tc, !critical temperature;
    pc, !critical pressure;
    vpa, vpb, vpc, vpd, !constants for the vapor pressure equation;
    cpa, cpb, cpc, cpd, !constants for the isobatic heat capacity equation;
    h0, !integration constant for enthalpy;
    xf; !mole fraction in feed stream;
  Tray : !the trays are numbered bottom upwards 
                     with reboiler as tray no. 1 and the condenser is the last tray;
    p, !pressure;
    t, !temperature;
    L, !molar flow rate of liquid;
    V, !molar flow rate of vapor;
    hl,!molar enthalpy of liquid ;
    hv,!molar enthalpy of vapor ;
    ft, !molar flow rate of feed entering each tray;
    zt, !binary variable to select the feed tray;
    hlq,!molar enthalpy of liquid; 
    hvp;!molar enthalpy of vapor;
  Mole_Frac (Tray, Component):
    x, !mole fraction of each component in liquid in each tray;
    y; !mole fraction of each component in vapor in each tray; 
DATA: Component = benzene, toluene; tray = 1..9; tb tc pc = !benzene; 353.2 562.2 48.9 !toluene; 383.8 591.8 41.0 vpa vpb vpc vpd = !benzene; -6.98273 1.33213 -2.62863 -3.33399 !toluene; -7.28607 1.38091 -2.83433 -2.79168 cpa cpb cpc cpd = !benzene; -3.392e+1 4.739e-1 -3.017e-4 7.130e-8 !toluene; -2.435e+1 5.125e-1 -2.765e-4 4.911e-8 xf = !benzene; 0.70 !toluene; 0.30 R = 8.314; !universal gas constant; treb = 380; !initial reboiler temperature; tbot = 375; !initial bottom temperature; ttop = 360; !initial top temperature; tcon = 355; !initial condenser temperature; preb = 1.2; !reboiler pressure; pbot = 1.12;!bottom-most tray pressure; ptop = 1.08;!top-most tray pressure; pcon = 1.05;!condenser pressure; f = 100; !total number of moles in the feed stream; pf = 1.12;!pressure of the feed stream; vf = 0.0; !vaport fraction in the feed stream (before expansion); ENDDATA SUBMODEL BUBBLE_POINT: ! assume that feed enters at bubble point this submodel solves for the bubble temperature, tbub; @SUM(Component(I): xf(I)*pc(I)* @EXP(tc(I)/tbub* (vpa(I)*(1-tbub/tc(I))+ vpb(I)*(1-tbub/tc(I))^1.5+ vpc(I)*(1-tbub/tc(I))^3+ vpd(I)*(1-tbub/tc(I))^6))/pf ) = 1 ; @BND(0,tbub,tb(@INDEX(toluene))); ENDSUBMODEL
SUBMODEL DIST_COL: !begin of the main model; !The objective chosen is indicative of the trade-offs between inceasing the throughput(P1) and the corresponding increase in reboiler duty - measure roughly by the value of the relux ratio; [obj] MAX = P1 - 50* reflux; !A. Phase Equilibrium; @FOR(Mole_Frac(I,J): y(I,J)*p(I) = x(I,J)*pc(J)* @EXP(tc(J)/t(I)* (vpa(J)*(1-t(I)/tc(J))+ vpb(J)*(1-t(I)/tc(J))^1.5+ vpc(J)*(1-t(I)/tc(J))^3+ vpd(J)*(1-t(I)/tc(J))^6));); !B. Phase equilibrium error; @FOR(Tray(I): @SUM(Component(J):x(I,J)) - @SUM(Component(J):y(I,J))=0;); !C. Total Material Balanaces (i.e, input equals output); !condensor, P1 is the top product rate; V(ncon-1)=(L(ncon)+P1); !other trays with possible feed; @FOR(Tray(I)|I#GT#nreb #AND# I#LT#ncon: V(I-1)+L(I+1)+ft(I)=V(I)+L(I)); !reboiler; L(nreb+1)=V(nreb)+L(nreb); !Do not feed in reboiler and condensor; @FOR(Tray(I)|I#EQ#nreb #OR# I#EQ#ncon: ft(I)=0;zt(I)=0); !reflux; L(ncon) = reflux * P1; !P2 is the bottom product rate; P2 = L(nreb); !D. Component Material Balances; !condensor; @FOR(Component(J): V(ncon-1)*y(ncon-1,J)=(L(ncon)+P1)*x(ncon,J)); !other trays with possible feed; @FOR(Tray(I)|I#GT#nreb #AND# I#LT#ncon: @FOR(Component(J): V(I-1)*y(I-1,J)+L(I+1)*x(I+1,J)+ft(I)*xf(J)= V(I)*y(I,J)+L(I)*x(I,J))); !reboiler; @FOR(Component(J): L(nreb+1)*x(nreb+1,J)=V(nreb)*y(nreb,J)+L(nreb)*x(nreb,J)); !E: Enthalpy Balances; !liquid enthalpy; @FOR(Tray(I): @FREE(hlq(I)); hlq(I) = @SUM(Component(J): x(I,J)* (cpa(J)*t(I) + cpb(J)/2*t(I)^2 + cpc(J)/3*t(I)^3 + cpd(J)/4*t(I)^4 - h0(J) + R*tc(J)*(vpa(J)*(1-t(I)/tc(J))+ vpb(J)*(1-t(I)/tc(J))^1.5+ vpc(J)*(1-t(I)/tc(J))^3+ vpd(J)*(1-t(I)/tc(J))^6) + R*t(I)*(vpa(J)+ 1.5*vpb(J)*(1-t(I)/tc(J))^0.5+ 3*vpc(J)*(1-t(I)/tc(J))^2+ 6*vpd(J)*(1-t(I)/tc(J))^5)))/hscal; ); !vapor enthalpy; @FOR(Tray(I): @FREE(hvp(I)); hvp(I) = @SUM(Component(J): y(I,J)* (cpa(J)*t(I) + cpb(J)/2*t(I)^2 + cpc(J)/3*t(I)^3 + cpd(J)/4*t(I)^4 -h0(J)))/hscal; ); !balance; !trays with possible feed; @FOR(Tray(I)|I#GT#nreb #AND# I#LT#ncon: V(I-1)*hvp(I-1)+L(I+1)*hlq(I+1)+ft(I)*hf/hscal=V(I)*hvp(I)+L(I)*hlq(I)); !F: Constraints on Feed Location; @FOR(Tray(I): @BIN(zt(I));); @SUM(Tray(I): zt(I))=1; @SUM(Tray(I): ft(I))=f; @FOR(Tray(I)|I#GT#nreb #AND# I#LT#ncon: ft(I)<=f*zt(I)); !G: Design Constraints; !reflux ratio; reflux <= 0.95; !purity of benzene in the distillate; x(ncon,@INDEX(benzene))>=0.95; !H: Variable Bounds; @FOR(Mole_Frac(I,J): @BND(0,X(I,J),1); @BND(0,Y(I,J),1);); @FOR(Tray(I): @BND(tcon-10,t(I),treb+10);); @FOR(Tray(I): @BND(0,ft(I),100);); @BND(0.01,reflux,0.95); @BND(50,P1,80); @BND(20,P2,50); @FOR(Tray(I): @BND((1-0.5*@SIGN(hlqlb))*hlqlb/hscal, hlq(I), (1+0.5*@SIGN(hlqub))*hlqub/hscal);); @FOR(Tray(I): @BND((1-0.5*@SIGN(hvplb))*hvplb/hscal, hvp(I), (1+0.5*@SIGN(hvpub))*hvpub/hscal);); ENDSUBMODEL
CALC: ntray = @SIZE(tray);!total number of trays; @FOR(Component(I): h0(I) = cpa(I)*tcon + cpb(I)/2*tcon^2 + cpc(I)/3*tcon^3 + cpd(I)/4*tcon^4); hlqlb = cpa(@INDEX(benzene))*tcon + cpb(@INDEX(benzene))/2*tcon^2 + cpc(@INDEX(benzene))/3*tcon^3 + cpd(@INDEX(benzene))/4*tcon^4 - h0(@INDEX(benzene)) + R*tc(@INDEX(benzene))*(vpa(@INDEX(benzene))*(1-tcon/tc(@INDEX(benzene)))+ vpb(@INDEX(benzene))*(1-tcon/tc(@INDEX(benzene)))^1.5+ vpc(@INDEX(benzene))*(1-tcon/tc(@INDEX(benzene)))^3+ vpd(@INDEX(benzene))*(1-tcon/tc(@INDEX(benzene)))^6) + R*tcon*(vpa(@INDEX(benzene))+ 1.5*vpb(@INDEX(benzene))*(1-tcon/tc(@INDEX(benzene)))^0.5+ 3*vpc(@INDEX(benzene))*(1-tcon/tc(@INDEX(benzene)))^2+ 6*vpd(@INDEX(benzene))*(1-tcon/tc(@INDEX(benzene)))^5); hlqub = cpa(@INDEX(toluene))*treb + cpb(@INDEX(toluene))/2*treb^2 + cpc(@INDEX(toluene))/3*treb^3 + cpd(@INDEX(toluene))/4*treb^4 - h0(@INDEX(toluene)) + R*tc(@INDEX(toluene))*(vpa(@INDEX(toluene))*(1-treb/tc(@INDEX(toluene)))+ vpb(@INDEX(toluene))*(1-treb/tc(@INDEX(toluene)))^1.5+ vpc(@INDEX(toluene))*(1-treb/tc(@INDEX(toluene)))^3+ vpd(@INDEX(toluene))*(1-treb/tc(@INDEX(toluene)))^6) + R*treb*(vpa(@INDEX(toluene))+ 1.5*vpb(@INDEX(toluene))*(1-treb/tc(@INDEX(toluene)))^0.5+ 3*vpc(@INDEX(toluene))*(1-treb/tc(@INDEX(toluene)))^2+ 6*vpd(@INDEX(toluene))*(1-treb/tc(@INDEX(toluene)))^5); hvplb = cpa(@INDEX(benzene))*tcon + cpb(@INDEX(benzene))/2*tcon^2 + cpc(@INDEX(benzene))/3*tcon^3 + cpd(@INDEX(benzene))/4*tcon^4 -h0(@INDEX(benzene)); hvpub = cpa(@INDEX(toluene))*treb + cpb(@INDEX(toluene))/2*treb^2 + cpc(@INDEX(toluene))/3*treb^3 + cpd(@INDEX(toluene))/4*treb^4 -h0(@INDEX(toluene)); hscal = @SMAX(@ABS(hlqlb),@ABS(hlqub),@ABS(hvplb),@ABS(hvpub)); nreb = 1;!tray index of the reboiler; ncon = ntray;!tray index of the condenser; p(nreb)=preb; p(ncon)=pcon; !assume a constant pressure gradient from the bottom to the top tray; @FOR(Tray(I)|I#GT#nreb #AND# I#LT#ncon: p(I) = pbot-(I-nreb-1)*(pbot-ptop)/((ncon-1)-(nreb+1));); !Use the global solver to solve for the bubble point; ! Make sure Global option is turned on; @SET( 'GLOBAL', 1); @SET( 'TERSEO', 1); @SOLVE(BUBBLE_POINT); @DIVERT(); @IFC( @STATUS() #EQ# 4: !set the feed temperature at bubble point; tf = tbub; @ELSE !Did not find a solution, call it quits; @STOP( ' Could not find the bubble temperature! ); !evaluate the enthalpy of feed; hf = @SUM(Component(I):xf(I)*(cpa(I)*tf + cpb(I)/2*tf^2 + cpc(I)/3*tf^3 + cpd(I)/4*tf^4 - h0(I) + R*tc(I)*(vpa(I)*(1-tf/tc(I))+ vpb(I)*(1-tf/tc(I))^1.5+ vpc(I)*(1-tf/tc(I))^3+ vpd(I)*(1-tf/tc(I))^6) + R*tf*(vpa(I)+ 1.5*vpb(I)*(1-tf/tc(I))^0.5+ 3*vpc(I)*(1-tf/tc(I))^2+ 6*vpd(I)*(1-tf/tc(I))^5))); !Solve for the best feed location; @SET( 'TERSEO', 0); @SET( 'GLOBAL', 1); @SOLVE(DIST_COL); ENDCALC END