comment
                              Columbia University
                                  ChE4510Y81
                    Evaporator Optimisation Program
                Vasos-Peter John Panagiotopoulos II
                  Shien Chi,Chris Lin, Vincent Ng
   Process Design  and Engineering, Prof. Amir N. Nahavandi

          Central algorithm derived from:
                    Carnahan,Seider & Katz,Chemical Engineering,v.2.,
                                                ________ ___________
                              "Compuers in Engineering Design".

                    M a j o r   V a r i a b l e s :
F          Flow rate for effect                    X          Mass ratio water to glycol
W          Weight fraxns. of glycol          T,P          Temperatures and pressures
A          Heat exchange area                    U          Heat transfer coefficients
Q          Heat flux                              QQ          Initial heat flux
WG          Glycol weight                              LOSS          Heat loss
TS          Saturation temperatures                    AQ!PRESSURE!BUBBLE Aqueous Bubble Pressure
I!F          Fixed investment                    L!MOL!FRAXN!GLYCOL Liquid mole fraxn.
BP!ELEVATION Boiling point elvtn.          HLPH,TPF!..          Steam tables
Remaining major variables are self descriptive.;

begin "EVAPRT:Evaporator Project" comment preload steam tables in outer block;
          preload!with 505.6,327.6,290.3,321.6,356.0,419.1,504.0,639.4,717.5,
                    838.3,957.3;real array hlph[0:10];
          preload!with .363,.415,.474,.541,.615,.698,.790,.893,1.01,1.13,
                    1.27,1.69,2.69,4.74,7.51,11.53,17.2,25,35.4,49.2;
                    real array tpf!pressure[0:19];
          preload!with 70,74,78,82,86,90,94,98,102,106,110,120,140,160,180,200,
                    220,240,260,280;real array tpf!temperature[0:19];
          integer n!effects;print("How many effects? > ");
          n!effects:=cvd(inchwl);
begin"evaporator"          boolean skip; integer i,j,k;
real array u[0:n!effects]; 
real accuracy,a!accuracy,opt!steam,opt!p,opt!venture!profit,year;
require "{}{}" delimiters;
define mol!fraxn!glycol(g)={g/62.07/(g/62.07+(1-g)/18)},
          mass!ratio(g)={(1-g)/g};
Require "sai:abbrev.sai"  source!file;comment sail abbreviations package.;
Require "omni:disprm.sai" source!file;comment sail omnigraph     package.;
Require "sai:routin.hdr"  source!file;comment sail routine       package.;
Require "sai:ttyio.hdr"   source!file;comment sail i/o           package.;
Real procedure realrd(string s); begin "realrd"          
          ! input real constants and assign to real variables;
          real f;string a;print (13&10&"Please input ",s,">  ");
          a:=inchwl;f:=realscan(a,32);return (f);end "realrd" ;
Procedure read!arai(reference real array a; integer n;string name);
          begin"araird" string ary;
              ! similarly input arrays;
             print(13&10&"Type in array,"&name&", elements sep by space,",
                    "<cr> when finished;",n,"elements."&13&10&" > ");
             ary:=inchwl;
             for i:=1 step 1 until n do a[i]:=realscan(ary,32);
          end"araird";
Real procedure enthalpy(real WG;string PHASE;real T);
          Begin "ENTHALPY" real m!mol!fraxn!glycol, h60, hlval, MWAV, HL, LAT;
                         real array mol!fraxn!glycol!vector[0:10];
                         label HOT,LINEAR,VAPOR,INTERP;
                    m!mol!fraxn!glycol:=mol!fraxn!glycol(wg); 
                    for j:=0 step 1 until 10 do mol!fraxn!glycol!vector[j]:=j/10;
                    if PHASE="heat" or phase="HEAT" then go HOT;
                    for j:=0 step 1 until 10 do 
                              if mol!fraxn!glycol!vector[j]>m!mol!fraxn!glycol 
                                        then go INTERP;
                              return(950);
       INTERP:           if j=0 then return(505);
                    h60:=hlph[j-1]+(hlph[j]-hlph[j-1])/(mol!fraxn!glycol!vector[j]
                              -mol!fraxn!glycol!vector[j-1])*(m!mol!fraxn!glycol
                                        -mol!fraxn!glycol!vector[j-1]);
                    hlval:=h60+(.5175*(t-60)+3.25@-4*(t^2-3600))*62.07
                              *m!mol!fraxn!glycol
                               +(.9375*(t-60)+25@-5*(t^2-3600))*18*
                                        (1-m!mol!fraxn!glycol);
                    mwav:=m!mol!fraxn!glycol*44.+(1.-m!mol!fraxn!glycol)*18;  
                    hl := hlval/mwav;
          hot:          lat:=1002 *((1 -(460 +T)/1165 )/(1 -620 /1165 ))^.38;
                          if  phase="LIQUID" or phase="liquid"  
                              then return (hl) 
                              else if  phase="VAPOR" or phase="vapor"  
                                      then return (hl+lat) else return(lat);
          END"ENTHALPY";

Real Procedure Hotpf (string wanted;real given);
          Begin"hotpf"
               real psi; label vapor;
                    if wanted="PRESSURE" then go vapor else psi := given*14.7;
                    If Tpf!pressure[0]>psi then return(70) 
                              else if  Tpf!pressure[19]<psi then return(280);
                    for j:=0 step 1 until 19 do  if tpf!pressure[j]>psi then 
                          return(tpf!temperature[j-1]+(tpf!temperature[j]
                                        -tpf!temperature[j-1])
                              /(tpf!pressure[j]-tpf!pressure[j-1])
                              *(psi-tpf!pressure[j-1]));
          vapor:          if tpf!temperature[0]>given then return(.025) else 
                              if tpf!temperature[19]<given then return(2.8);
                    for j:=0 step 1 until 19 do 
                              if tpf!temperature[j]>given then 
                                   RETURN((TPF!PRESSURE[J-1]+(TPF!PRESSURE[J] 
                                                  -TPF!PRESSURE[J-1])
                                        /(TPF!TEMPERATURE[J]-TPF!TEMPERATURE[J-1])
                                              *(GIVEN-TPF!TEMPERATURE[J-1]))/14.7);
          END"hotpf";
real procedure f(real p!test );
begin"central procedure"
label calculate,loop,bypass,each!effect;
real array enthalpy!vapor,f,x,w,t,p,liquid!enthalpy,yh,aq!pressure!bubble,ts,
          steam!enthalpy!latent,temp!difference,bp!elevation,vapor!rate,
          new!vapor!rate,l!mol!fraxn!glycol,a,q[0:n!effects];
real guess!bp!elevation,totvap,sum!bp!elevation,
          sum!temp!difference,sum!vap,det,sum!a,qq,average!a,wg,loss,
        steam!econ,condenser!load,venture!profit;
!                           **************
                            *MAIN PROGRAM*
                              **************;
guess!bp!elevation:=10;arrclr(bp!elevation);
  loss:=sum!bp!elevation:=sum!temp!difference:=sum!vap:=0;
  f[0]:=1@4; w[0]:=.02; w[n!effects]:=.9; t[0]:=100;
  aq!pressure!bubble[0]:=2;
!                    OPTIMISING P;          p[n!effects]:=p!test;
totvap:=f[0]*w[0]*((x[0]:=mass!ratio(w[0]))
          -(x[n!effects]:=mass!ratio(w[n!effects])));
  accuracy:=totvap*(a!accuracy);skip:=0;
aq!pressure!bubble[n!effects]:=p[n!effects]
          /(1-(l!mol!fraxn!glycol[n!effects]:=mol!fraxn!glycol(w[n!effects])));
t[n!effects]:=hotpf("temperature",aq!pressure!bubble[n!effects]);          
ts[0]:=hotpf("temperature",aq!pressure!bubble[0]);
steam!enthalpy!latent[0]:=enthalpy(0,"heat",ts[0]);
          for i:=1 step 1 until n!effects do 
                    begin"averaging conditions"
                              bp!elevation[i]:=guess!bp!elevation;
                              sum!bp!elevation:=(n!effects-1)*guess!bp!elevation;
                              temp!difference[i]:=(ts[0]-sum!bp!elevation
                                        -t[n!effects])/n!effects;
                              vapor!rate[i]:=totvap/n!effects;
                    end"averaging conditions";
sum!bp!elevation:=0;
loop:                              liquid!enthalpy[0]:=enthalpy(w[0],"liquid",t[0]);
          for i:=1 step 1 until n!effects do 
                    begin "effect conditions"
                              sum!bp!elevation:=sum!bp!elevation+bp!elevation[i-1];
                              sum!temp!difference:=sum!temp!difference
                                        +temp!difference[i];
                              sum!vap:=sum!vap+vapor!rate[i];
                              t[i]:=ts[0]-sum!bp!elevation-sum!temp!difference;
                              x[i]:=x[0]-sum!vap/(f[0]*w[0]);
                              w[i]:=1/(1+x[i]);
                              l!mol!fraxn!glycol[i]:=mol!fraxn!glycol(w[i]);
                              liquid!enthalpy[i]:=enthalpy(w[i],"liquid",t[i]);
                              enthalpy!vapor[i]:=enthalpy(w[i],"vapor" ,t[i]);
                              aq!pressure!bubble[i]:=hotpf("PRESSURE",t[i]);
                              p[i]:=aq!pressure!bubble[i]*(1-l!mol!fraxn!glycol[i]);
                              ts[i]:=hotpf("temperature",p[i]);
                              steam!enthalpy!latent[i]:=enthalpy(0,"heat",ts[i]);
                              f[i]:=f[0]-f[0]*w[0]*(x[0]-x[i]);                              
                    end   "effect conditions";
sum!bp!elevation:=0;
                    for i:=1 step 1 until n!effects-1 do
                              begin "BPE summation"
                                    bp!elevation[i]:=t[i]-ts[i];
                                    sum!bp!elevation:=sum!bp!elevation
                                        +bp!elevation[i];
                              end "BPE summation";
                    if skip then go bypass;
                         for k:=1 step 1 until n!effects do
                          temp!difference[k]:=(ts[0]-t[n!effects]-sum!bp!elevation)
                                        /n!effects;
                    bypass: sum!bp!elevation:=sum!vap:=sum!temp!difference:=0;
                                     if n!effects=1 then 
                                                  begin "1 effect" 
                                                            vapor!rate[1]:=
                                                               new!vapor!rate[1]:=totvap;
                                                            go calculate;
                                                  end "1 effect";
                                       if n!effects=2 then 
                                                  begin "2 effects"
                                                            new!vapor!rate[1]:=
                                                                 (totvap*enthalpy!vapor[2]
                                                                   +f[2]*liquid!enthalpy[2]
                                                                         +loss-f[1]
                                                                   *liquid!enthalpy[1])
                                                                 /(steam!enthalpy!latent[1]
                                                                   +enthalpy!vapor[2]);
                                                            new!vapor!rate[2]:=totvap
                                                                      -new!vapor!rate[1];
                                                            if abs(vapor!rate[1]
                                                                 -new!vapor!rate[1])
                                                                <accuracy
                                                                           then go calculate;
                                                            for k:=1 step 1 
                                                                until n!effects do
                                                                    vapor!rate[k]:=
                                                                           new!vapor!rate[k];
                                                            go loop;          
                                                  end   "2 effects";
                                     if n!effects=3 then
                                                  begin "3 effect" ! load forCramer soln;
                                                  real array d!matr[1:2,0:2];
                                                            d!matr[1,0]:=liquid!enthalpy[2]
                                                                     *f[2]
                                                                   -liquid!enthalpy[1]*f[1]
                                                                   +loss;
                                                            d!matr[1,1]:=
                                                                 steam!enthalpy!latent[1];
                                                            d!matr[1,2]:=
                                                                      -enthalpy!vapor[2];
                                                            d!matr[2,0]:=liquid!enthalpy[3]
                                                                      *f[3]
                                                                   -liquid!enthalpy[2]*f[2]
                                                                   +totvap
                                                                         *enthalpy!vapor[3]
                                                                   +loss;
                                                            d!matr[2,1]:=enthalpy!vapor[3];
                                                            d!matr[2,2]:=
                                                                   steam!enthalpy!latent[2]
                                                                 +enthalpy!vapor[3];
                                                            det:=d!matr[1,1]*d!matr[2,2]
                                                                      -d!matr[2,1]*d!matr[1,2];
                                                            new!vapor!rate[1]:=
                                                                 (d!matr[1,0]*d!matr[2,2]
                                                                  -d!matr[2,0]*d!matr[1,2])
                                                                 /det;
                                                            new!vapor!rate[2]:=
                                                                 (d!matr[1,1]*d!matr[2,0]
                                                                  -d!matr[2,1]*d!matr[1,0])
                                                                       /det;
                                                                   new!vapor!rate[3]:=totvap
                                                                -new!vapor!rate[1]
                                                                -new!vapor!rate[2];
                                                            if abs(vapor!rate[1]
                                                                      -new!vapor!rate[1])
                                                                 <accuracy
                                                               and abs(vapor!rate[2]
                                                                      -new!vapor!rate[2])
                                                                 <accuracy 
                                                               then go calculate 
                                                               else 
                                                                      begin"newvap"
                                                                        for k:=1 step 1 
                                                                             until n!effects do
                                                                        vapor!rate[k]:=
                                                                           new!vapor!rate[k];
                                                                        go loop;          
                                                                      end"newvap";
                                                     end   "3 effect";
Calculate:                    sum!a:=0;
                    q[1]:=qq:=(vapor!rate[1]:=new!vapor!rate[1])*enthalpy!vapor[1]
                              +f[1]*liquid!enthalpy[1]-f[0]*liquid!enthalpy[0]+loss;
                    sum!a:=a[1]:=qq/u[1]/temp!difference[1];
                    vapor!rate[0]:=qq/steam!enthalpy!latent[0];
                    if n!effects=1 then go each!effect;
                    for i:=2 step 1 until n!effects do
                              begin "multi"
                                        vapor!rate[i-1]:=new!vapor!rate[i-1];
                                        q[i]:=vapor!rate[i-1]
                                                  *steam!enthalpy!latent[i-1];
                                        a[i]:=q[i]/u[i]/temp!difference[i];
                                      sum!a:=sum!a+a[i];
                              end"multi";
each!effect:                    average!a:=sum!a/n!effects;
                    for i:=1 step 1 until n!effects do 
                           if abs(a[i]-average!a)>a!accuracy*average!a then 
                                  begin" each effect "
                                        for k:=1 step 1 until n!effects do 
                                                  temp!difference[k]:=temp!difference[k]
                                                            *a[k]/average!a;
                                        skip:=1;go loop;
                                  end" each effect ";
                    steam!econ:=totvap/vapor!rate[0];
                    condenser!load:=vapor!rate[n!effects]*
                              steam!enthalpy!latent[n!effects];
                    p[0]:=aq!pressure!bubble[0];
                    begin "economics"
                                        ! Correlation for Inflation:
                                                  WPI($1961) from 1947 to 1980:
                                                  .810,.876,.833,.866,.964,.938,.925,
                                                  .927,.929,.960,.987,1.001,1.003,1.004,
                                                  1.000,1.003,1.000,1.002,1.022,1.056,
                                                  1.058,1.085,1.127,1.168,1.206,1.259,
                                                  1.426,1.693,1.850,1.939,2.083,2.246,
                                                2.448,2.747.
                                                  Sources:Organisation for Economic
                                                    Co-operation and Development,
                                                    March,1981,Main Economic Indicators
                                                                 ________________________
                                                   U.S.Dept. Commerce, Bur. Econ. Anal.,
                                                    Business Statistics,1977(v.21)
                                                    ___________________ 
                                            Fitting this data to an exponential curve:
                                                  110@-27*exp(year*29@-3)
                                                    sum of squares=780@-3
                                                     use of program derived from listing
                                                       HP41C Standard Applications Manual
                                                       __________________________________
                                                          Corvallis,OR:HP,1979;
                     real density,i!f,cost;
                      density:=p[n!effects]/(hotpf("PRESSURE",p[n!effects])+460)/
                              .7302; ! Ideal gas law ;
                      i!f:=2.9*n!effects*(.24@6+.22@3*a[1]^.61);
                    ! Rudd & Watson, Strategy of Process Engineering,NY:Wiley,1968.
                                 _______________________________ 
                        Bauman, Fundamentals of Cost Engineering,NY:Reinhold, 1964.
                                ________________________________
                       ;    
                      cost:=(.2*i!f +1633+6.57*vapor!rate[n!effects]
                           +.087*2.9*abs(p[n!effects]-p[0])/density
                               *f[n!effects]*8000)
                                    /.95;
                      venture!profit:=(-.5*cost-.2*i!f)*110@-27*exp(year*29@-3);
                    end   "economics";
          return(venture!profit);
end"central procedure";
real procedure fibonacci(real aa,bb);
          begin"fibonacci" 
          !  MAXIMISATION
             Algorithm due to:
                    Process Design and Engineering, Prof. A.N.Nahavandi,
                    ______________________________
                              class notes,ChE4510Y81,1981,Columbia Univ.
                    Beveridge and Schechter,Optimisation,Mc Graw Hill
                                                  ____________
                    Van Nostrand Scientific Encyclopedia,VNR
                    ____________________________________
                    ;
          real new!fibonacci,fibonacci!old,temporary,tau,l,r,a,b;
              fibonacci!old:=0;new!fibonacci:=1;l:=a:=aa;r:=b:=bb;
               while abs((l-r)/r)>a!accuracy do
                    begin"fibonacci section"
                              new!fibonacci:=fibonacci!old+new!fibonacci
                                        +0*(fibonacci!old:=new!fibonacci);
                              tau:=fibonacci!old/new!fibonacci;
                              l:=b-tau*abs(b-a);r:=a+tau*abs(b-a);
                              if f(l)>f(r) then b:=r else a:=l;
                              if a<aa or b>bb then done;
                    end"fibonacci section";
           return((a+b)/2);
          end"fibonacci";
  read!arai(u,n!effects,"Heat Transfer Coefficient");
  a!accuracy:=realrd("Area accuracy ");
  year:=realrd("year for which evaluation is to be performed");
  Print(13&10&"Pressure guess should be between .025 and 2.8 atm ");
  opt!p:=fibonacci(realrd("Lower Boundary"),
                       realrd("Upper Boundary"));
             opt!venture!profit:=f(opt!p);
            print(13&10&"Optimum(max): P,Venture Profit:",opt!p,
                                                                  opt!venture!profit);
end"evaporator";          end "EVAPRT:Evaporator Project";
comment

          Mathematics is music,
                    FORTRAN is a drum,
                    Array-oriented languages are pianos,
                              MLAB is a Yamaha
                              APLSF is a Steinway
                    Structured languages are violins
                              SAIL is a Stradivarius.


                    Evaporator Steam Economy Determination Program 
                              Panagiotopoulos,Chi,Lin, & Ng
                                  FLOW CHART
_________________________________________________________________________
                    I/O # EFFECTS& THUS SET DIMENSIONS
                                        |
                                     \|/
                    I/O ENDPOINT:FLOW,T,P,VAPOR RATES,ALL:U AS PRECONDITIONS
                                        |
                                     \|/
                    COMPUTE TEMPERATURES,FLOW RATES,ENTHALPIES
                                        |
                                     \|/
                    FOR EACH EFFECT
                              DETERMINE BOILING POINT ELEVATION,VAPOR RATE,
                                        TEMPERATURE CHANGE
                                        |
                                     \|/
LOOP:                    FOR EACH EFFECT
                              EVALUATE TEMPERATURE,COMPOSITION,ENTHALPY,PRESSURE,
                                        LATENT HEATS,AND FLOWS
                                        |
                                     \|/
                    FOR EACH EFFECT EVALUATE BOILING POINT ELEVATION
                                        |
                                     \|/
                    FOR EACH EFFECT DETERMINE THE TEMPERATURE DROP
                                        |
                                     \|/
                    HOW MANY EFFECTS ARE THERE(VALUE INPUT EARLIER)
                                        |
          _______________________\|/________________________
          |                              |                               |
       \|/1                           \|/2                            \|/3
          VAPOR RATE IS SAME          VAPOR RATE FROM                     VAPOR RATE FROM
          AS TOTAL                    ENTHALPIES                     ENTHALPIES
          |_______________________|________________________|
                                     \|/
                                 TOLERANCE?-----(NO)---APPLY NEW VAPOR RATES,GO LOOP
                                        |
                                     \|/
                               DETERMINE AREAS
                               ARE TOLERABLE?----(NO)---MODIFY T DROPS VIA AREAS
                                        |                              THEN GO TO LOOP
                                     \|/
                    DETERMINE STEAM ECONOMY AND CONDENSER LOAD----->PRINTOUT
                                        |
                                     \|/
                                     END

;