(* An Example in ASCEND to do Flash Calculations

Remember that everything between (* and *) is a comment

This uses the flash.a4l library and associated libraries in the
models directory of the ASCEND distribution -

This example written by Krishnan Chittur (kchittur@che.uah.edu) -
Chemical Engineering Department,
University of Alabama in Huntsville - BUT -

I have borrowed liberally from the examples in the models directory -
so any real credit goes to a lot of other people! - Errors are ofcourse
mine - I'd appreciate feedback *)

(* Please keep in mind that there is a lot of room for improvement in
   these models, methods - again, I'd appreciate specific feedback *)

REQUIRE "flash.a4l";

(* the flash.a4l library will in turn provide models, methods and
such provided in stream_holdup.a4l, thermodynamics.a4l,
components.a4l, phases.a4l,atoms.a4l, measures.a4l,
system.a4l and basemodel.a4l *)

MODEL do_a_flash;

(* You first need to tell the program what components are going 
to be in the flash - and what kind of thermodynamic model you want
to use for modeling the Flash - i.e. Can you consider the liquid phase
to be ideal?  Vapor phase ideal?  This is a decision YOU must make *)  
    
    cd IS_A components_data(['benzene','toluene'],'benzene');

(* Notice we have ONLY two components - and we have, arbitrarily 
   assigned toluene as the reference component - i.e. it will have 
   a relative volatility of 1.0 - all other relative volatilities will
   be computed with respect to toluene *)

    pdV IS_A phases_data('V', 'ideal_vapor_mixture', 'none', 'none');

    pdL IS_A phases_data('L', 'none', 'UNIFAC_liquid_mixture', 'none');
    pdVL IS_A phases_data('VL', 'ideal_vapor_mixture',
				'UNIFAC_liquid_mixture', 'none');

(* In this release of ASCEND, you have the option of modeling the 
   liquid mixture as a UNIFAC_liquid_mixture - where the mixture 
   properties are calculated using the group contribution method
   UNIFAC - OR you can model liquid mixtures using the 
   Wilson_liquid_mixture - the vapor phase can be modeled as 
   ideal_vapor_mixture OR as Pitzer_vapor_mixture - Not all the
   Wilson parameters are available at the time of this writing -
   this may change ofcourse - You can take a peek at the
   thermodynamics.a4l and components.a4l file - Add the parameters
   yourself *)

   equilibrated IS_A start_false;

(* Qin is unit specific *)
    Qin IS_A energy_rate;

    Feed IS_A stream(cd, pdVL, equilibrated);
    liqout IS_A stream(cd, pdL, equilibrated);
    vapout IS_A stream(cd, pdV, equilibrated);
    Flash IS_A vapor_liquid_flash(Qin, equilibrated, Feed, vapout, liqout);
    boundwidth IS_A bound_width;

(*  Just for the heck of it, I wanted to calculate the activity coefficient
    of benzene and toluene from the composition of the vapor phase,
    liquid phase, total pressure and the vapor pressure - Look how I
    can access these different values *)

    vpbenzene, vptoluene IS_A pressure;
    vpbenzene = Flash.liqout.phase['liquid1'].pure['benzene'].VP;
    vptoluene = Flash.liqout.phase['liquid1'].pure['toluene'].VP;

    ybenzene, ytoluene, xbenzene, xtoluene IS_A fraction;
    ybenzene = vapout.y['benzene'];
    ytoluene = vapout.y['toluene'];
    xbenzene = liqout.y['benzene'];
    xtoluene = liqout.y['toluene'];

    TotalPressure IS_A pressure;
    TotalPressure = vapout.P;

(* notice that TotalPressure = liqout.P is also OK - you do know why right? *)
(* Also TotalPressure = Flash.liqout.P is also OK *)

    gammabenzene, gammatoluene IS_A factor;
    gammabenzene = vapout.y['benzene']*vapout.P/(vpbenzene*liqout.y['benzene']);
    gammatoluene = vapout.y['toluene']*vapout.P/(vptoluene*liqout.y['toluene']);

METHODS

METHOD default_self;
END default_self;

METHOD ClearAll;
	equilibrated := FALSE;
	Feed.P.fixed := FALSE;
	Feed.T.fixed := FALSE;
	Feed.flow.fixed := FALSE;
	Feed.y['benzene'].fixed := FALSE;
	Feed.phase['vapor'].alpha['toluene'].fixed := FALSE;
	Feed.phase['vapor'].alpha['benzene'].fixed := FALSE;
	Flash.alpha['toluene'].fixed := FALSE;
	Flash.alpha['benzene'].fixed := FALSE;
	Flash.state.phase_fraction['vapor'].fixed := FALSE;
	Flash.P.fixed := FALSE;
	Flash.T.fixed := FALSE;
	Flash.Qin.fixed := FALSE;
END ClearAll;

METHOD specify_P_Calculate_T;
	equilibrated := TRUE;
	Feed.P.fixed := TRUE;
	Feed.T.fixed := TRUE;
	Feed.flow.fixed := TRUE;
	Feed.y['benzene'].fixed := TRUE;
	Feed.phase['vapor'].alpha['toluene'].fixed := TRUE;
	Feed.phase['vapor'].alpha['benzene'].fixed := FALSE;
	Flash.alpha['toluene'].fixed := FALSE;
	Flash.alpha['benzene'].fixed := TRUE;
	Flash.state.phase_fraction['vapor'].fixed := TRUE;
	Flash.P.fixed := TRUE;
	Flash.T.fixed := FALSE;
	Flash.Qin.fixed := FALSE;
END specify_P_Calculate_T;

METHOD specify_T_Calculate_P;

	Feed.P.fixed := TRUE;
	Feed.T.fixed := TRUE;
	Feed.flow.fixed := TRUE;
	Feed.y['benzene'].fixed := TRUE;
	Feed.phase['vapor'].alpha['toluene'].fixed := TRUE;
	Feed.phase['vapor'].alpha['benzene'].fixed := FALSE;
	Flash.alpha['toluene'].fixed := FALSE;
	Flash.alpha['benzene'].fixed := TRUE;
	Flash.state.phase_fraction['vapor'].fixed := TRUE;
	Flash.P.fixed := FALSE;
	Flash.T.fixed := TRUE;
	Flash.Qin.fixed := FALSE;
	equilibrated := TRUE;

END specify_T_Calculate_P;

METHOD specify_PT_Flash;
	Feed.P.fixed := TRUE;
	Feed.T.fixed := TRUE;
	Feed.flow.fixed := TRUE;
	Feed.y['benzene'].fixed := TRUE;
	Feed.phase['vapor'].alpha['toluene'].fixed := TRUE;
	Feed.phase['vapor'].alpha['benzene'].fixed := FALSE;
	Flash.alpha['toluene'].fixed := TRUE;
	Flash.alpha['benzene'].fixed := FALSE;
	Flash.state.phase_fraction['vapor'].fixed := FALSE;
	Flash.P.fixed := TRUE;
	Flash.T.fixed := TRUE;
	Flash.Qin.fixed := FALSE;
	equilibrated := TRUE;
END specify_PT_Flash;

METHOD values;
	Feed.P := 1 {atm};
	Feed.T := 298 {K};
	Feed.flow := 100 {mole/sec};
	Feed.y['benzene'] := 0.2;
	Feed.phase['vapor'].alpha['toluene'] := 1;
	Feed.phase['vapor'].alpha['benzene'] := 2;
	Flash.alpha['toluene'] := 1;
	Flash.alpha['benzene'] := 2;
	Flash.state.phase_fraction['vapor'] := 0.25;
	Flash.P := 1.5 {atm};
	Flash.T := 387.15 {K};
	Flash.Qin := 100.00 {kW};
END values;

METHOD seqmod_P_Calculate_T;
RUN Flash.ClearAll;
RUN ClearAll;
RUN specify_P_Calculate_T;
RUN values;
END seqmod_P_Calculate_T;

METHOD seqmod_T_Calculate_P;
RUN Flash.ClearAll;
RUN ClearAll;
RUN specify_T_Calculate_P;
RUN values;
END seqmod_T_Calculate_P;

METHOD seqmod_PT_Flash;
RUN Flash.ClearAll;
RUN ClearAll;
RUN specify_PT_Flash;
RUN values;
END seqmod_PT_Flash;

END do_a_flash;