(* This is written in the ASCEND Language ... ASCEND is an
environment to pose modeling problems and solve the resulting
equations in an object oriented fashion ... ASCEND is
copyrighted by Carnegie Mellon University ... the project
has been directed by Prof. Art Westerberg in the Chemical
Engineering Department ...*)

(* ASCEND is free software, released under GNU ... and you
can get it from http://www.cs.cmu.edu/~ascend ... The
software has been released WITH source code ... it compiles
under Solaris, Linux and IRIX (SGI) ... thanks to the excellent
job in putting the software together by the folks at CMU! ...
and autoconf ... I have no idea of how people compiled things
before ./configure ! *)

(* This example written by Krishnan Chittur (chitturk@uah.edu)
Chemical Engineering Department, University of Alabama in Huntsville,
Huntsville, AL 35899 (256) 890 6850 (V), (256) 890 6839 (FAX) *)

(* Models for Bubble Point, Dew Point and Flash Calculations
   Bubble Point T Calculations ...
Specify P, Antoine Coefficients, Feed Composition ...
Calculate Temperature and Composition of the First few drops
of Vapor formed ...

   Dew Point T Calculations ...
Specify P, Antoine Coefficients, Feed Composition ...
Calculate Temperature and Composition of the First few drops
of liquid formed ...

   Flash
Specify P,T, Antoine Coefficients, Feed Composition ...
Calculate the amount of Vapor and Liquid Formed and their
compositions

... we assume ideal gas, ideal liquid
... obviously it will be more interesting
to include activity coefficient models, fugacity calculations
from equation of state ...

   We also have only two components here - you can easily extend it
   systems with more components

Look at VanNess and Smith - Thermodynamics, 5th Edition -
page 444 *)

REQUIRE "atoms.a4l";

MODEL antoine_coefficients;

	components IS_A set OF symbol_constant;
	components :== ['acetonitrile','nitromethane'];

	Antoine_A[components],
	Antoine_B[components],
	Antoine_C[components] IS_A factor;

METHODS
	METHOD specify;
	Antoine_A[components].fixed := TRUE;
	Antoine_B[components].fixed := TRUE;
	Antoine_C[components].fixed := TRUE;
	END specify;

	METHOD values;
	Antoine_A['acetonitrile'] := 14.2724;
	Antoine_A['nitromethane'] := 14.2043;
	Antoine_B['acetonitrile'] := 2945.47;
	Antoine_B['nitromethane'] := 2972.64;
	Antoine_C['acetonitrile'] := 224.00;
	Antoine_C['nitromethane'] := 209.00;
	END values;

END antoine_coefficients;

MODEL bubble_point;

	components IS_A set OF symbol_constant;
	components :== ['acetonitrile','nitromethane'];
	T, T_degC IS_A temperature;
	P IS_A pressure;
	Kvalue[components] IS_A positive_factor;

	antoine IS_A antoine_coefficients;

	Vapor_Pressure[components] IS_A pressure;

	z[components] IS_A fraction;
	y[components] IS_A fraction;

	FOR i IN components CREATE 	
	eq1[i]: ln(Vapor_Pressure[i]/1.0{kPa})
	= antoine.Antoine_A[i]
	- antoine.Antoine_B[i]/(T_degC + antoine.Antoine_C[i]);
	END FOR;

	T_degC = T - 273.15{K};

	FOR i IN components CREATE 	
	eq2[i]:Kvalue[i] = Vapor_Pressure[i]/P;
	END FOR;

	eq3: SUM[Kvalue[i]*z[i] | i IN components] = 1.0;

	eq4: SUM[z[i] | i IN components] = 1.0;

	FOR i IN components CREATE
	eq5[i]: y[i] = Kvalue[i]*z[i];
	END FOR;

METHODS
	METHOD clear;
	z[components].fixed := FALSE;
	y[components].fixed := FALSE;
	T.fixed := FALSE;
	P.fixed := FALSE;
	Kvalue[components].fixed := FALSE;
	Vapor_Pressure[components].fixed := FALSE;
	END clear;

	METHOD specify;
	z[components].fixed := FALSE;
	z[CHOICE[components]].fixed := TRUE;
	y[components].fixed := FALSE;
	T.fixed := FALSE;
	P.fixed := TRUE;
	Kvalue[components].fixed := FALSE;
	Vapor_Pressure[components].fixed := FALSE;
	END specify;

	METHOD values;
	P := 70.0 {kPa};
	z['acetonitrile'] := 0.7378;
	END values;

	METHOD seqmod;
	RUN antoine.clear;
	RUN antoine.specify;
	RUN antoine.values;
	RUN clear;
	RUN specify;
	RUN values;
	END seqmod;

END bubble_point;

MODEL dew_point;

	components IS_A set OF symbol_constant;
	components :== ['acetonitrile','nitromethane'];
	T, T_degC IS_A temperature;
	P IS_A pressure;
	Kvalue[components] IS_A positive_factor;

	Vapor_Pressure[components] IS_A pressure;
	antoine IS_A antoine_coefficients;
	z[components] IS_A fraction;
	x[components] IS_A fraction;

	FOR i IN components CREATE 	
	eq1[i]: ln(Vapor_Pressure[i]/1.0{kPa})
	= antoine.Antoine_A[i]
	- antoine.Antoine_B[i]/(T_degC + antoine.Antoine_C[i]);
	END FOR;

	T_degC = T - 273.15{K};

	FOR i IN components CREATE 	
	eq2[i]:Kvalue[i] = Vapor_Pressure[i]/P;
	END FOR;

	eq3: SUM[z[i]/Kvalue[i] | i IN components] = 1.0;

	eq4: SUM[z[i] | i IN components] = 1.0;

	FOR i IN components CREATE
	eq5[i]: x[i] = z[i]/Kvalue[i];
	END FOR;

METHODS
	METHOD clear;
	z[components].fixed := FALSE;
	x[components].fixed := FALSE;
	T.fixed := FALSE;
	P.fixed := FALSE;
	Kvalue[components].fixed := FALSE;
	Vapor_Pressure[components].fixed := FALSE;
	END clear;

	METHOD specify;
	z[components].fixed := FALSE;
	z[CHOICE[components]].fixed := TRUE;
	x[components].fixed := FALSE;
	T.fixed := FALSE;
	P.fixed := TRUE;
	Kvalue[components].fixed := FALSE;
	Vapor_Pressure[components].fixed := FALSE;
	END specify;

	METHOD values;
	P := 70.0 {kPa};
	z['acetonitrile'] := 0.7378;
	END values;

	METHOD seqmod;
	RUN antoine.clear;
	RUN antoine.specify;
	RUN antoine.values;
	RUN clear;
	RUN specify;
	RUN values;
	END seqmod;

END dew_point;

MODEL pt_flash;

	components IS_A set OF symbol_constant;
	components :== ['acetonitrile','nitromethane'];
	T, T_degC IS_A temperature;
	P IS_A pressure;
	Kvalue[components] IS_A positive_factor;
	antoine IS_A antoine_coefficients;
	Vapor_Pressure[components] IS_A pressure;

	z[components] IS_A fraction;
	x[components] IS_A fraction;
	y[components] IS_A fraction;

	F, L, V IS_A mole;

	FOR i IN components CREATE
	eq1[i]: z[i]*F = y[i]*V + x[i]*L;
	END FOR;

	FOR i IN components CREATE 	
	eq2[i]: ln(Vapor_Pressure[i]/1.0{kPa})
	= antoine.Antoine_A[i]
	- antoine.Antoine_B[i]/(T_degC + antoine.Antoine_C[i]);
	END FOR;

	T_degC = T - 273.15{K};

	FOR i IN components CREATE 	
	eq3[i]:Kvalue[i] = Vapor_Pressure[i]/P;
	END FOR;

	eq4: F = L + V;

	eq5: SUM[z[i] | i IN components] = 1.0;

	eq6: SUM[y[i] | i IN components] - SUM[x[i] | i IN components] = 0.0;

	FOR i IN components CREATE 	
	eq7[i]: y[i] = Kvalue[i]*x[i];
	END FOR;

METHODS
	METHOD clear;
	z[components].fixed := FALSE;
	x[components].fixed := FALSE;
	y[components].fixed := FALSE;
	T.fixed := FALSE;
	P.fixed := FALSE;
	Kvalue[components].fixed := FALSE;
	Vapor_Pressure[components].fixed := FALSE;
	F.fixed := FALSE;
	L.fixed := FALSE;
	V.fixed := FALSE;
	END clear;

	METHOD specify;
	z[components].fixed := FALSE;
	z[CHOICE[components]].fixed := TRUE;
	x[components].fixed := FALSE;
	y[components].fixed := FALSE;
	T.fixed := TRUE;
	P.fixed := TRUE;
	Kvalue[components].fixed := FALSE;
	Vapor_Pressure[components].fixed := FALSE;
	F.fixed := TRUE;
	L.fixed := FALSE;
	V.fixed := FALSE;
	END specify;

	METHOD values;
	P := 70.0 {kPa};
	T := 352 {K};
	z['acetonitrile'] := 0.50;
	F := 100 {mole};
	END values;

	METHOD seqmod;
	RUN antoine.clear;
	RUN antoine.specify;
	RUN antoine.values;
	RUN clear;
	RUN specify;
	RUN values;
	END seqmod;

END pt_flash;