#*-------------------------------------------------------------------
* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
*
* This LIBRARY is free software; you can distribute it and/or modify
* it under the therms of the ALSOC FREE LICENSE as available at
* http://www.enq.ufrgs.br/alsoc.
*
* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
* from http://www.rps.eng.br Copyright (C) 2002-2004.
* All rights reserved.
*
* EMSO is distributed under the therms of the ALSOC LICENSE as
* available at http://www.enq.ufrgs.br/alsoc.
*
*-------------------------------------------------------------------
* Model of a static PH flash
*-------------------------------------------------------------------- 
* 	- Streams
*		* a liquid outlet stream
*		* a vapour outlet stream
*		* a feed stream
*
*	- Assumptions
*		* both phases are perfectly mixed
*
*	- Specify:
*		* the feed stream;
*		* the heat duty
*       * the outlet pressure
*
*----------------------------------------------------------------------
* Author: Rafael de P. Soares and Paula B. Staudt
* $Id: flashPH.mso 101 2007-01-09 15:59:59Z rafael $
*--------------------------------------------------------------------*#

using "streams";

Model FlashPHSteady
	PARAMETERS
ext PP as CalcObject;
ext NComp as Integer;
	B as Real(Default=1000, Brief="Regularization Factor");
	
	VARIABLES
in	Inlet as stream; #(Brief="Feed Stream");
out	OutletL as stream_therm; #(Brief="Liquid outlet stream");
out	OutletV as stream_therm; #(Brief="Vapour outlet stream");
in	Q as heat_rate (Brief="Rate of heat supply"); 
	vfrac as fraction(Brief="Real vaporization fraction");
	vsat as Real(Lower=-5, Upper=5, Brief="Vaporization fraction if saturated");
	Tsat as temperature(Lower=173, Upper=473, Brief="Temperature if saturated");
	xsat(NComp) as Real(Lower=-5, Upper=5, Brief="Liquid composition if saturated");
	ysat(NComp) as Real(Lower=-5, Upper=5, Brief="Vapour composition if saturated");
	
	zero_one as fraction(Brief="Regularization Variable");
	one_zero as fraction(Brief="Regularization Variable");

	EQUATIONS
	"Chemical equilibrium"
	PP.LiquidFugacityCoefficient(Tsat, OutletL.P, xsat)*xsat = 
		PP.VapourFugacityCoefficient(Tsat, OutletV.P, ysat)*ysat;

	"Global Molar Balance"
	Inlet.F = OutletV.F + OutletL.F;
	OutletV.F = Inlet.F * vfrac;

	"Component Molar Balance"
	Inlet.F*Inlet.z = OutletL.F*xsat + OutletV.F*ysat;
	sum(xsat) = sum(ysat);

	"Energy Balance if saturated"
	Inlet.F*Inlet.h  + Q =
		Inlet.F*(1-vsat)*PP.LiquidEnthalpy(Tsat, OutletL.P, xsat) +
		Inlet.F*vsat*PP.VapourEnthalpy(Tsat, OutletV.P, ysat);

	"Real Energy Balance"
	Inlet.F*Inlet.h  + Q =
		Inlet.F*(1-vfrac)*OutletL.h + Inlet.F*vfrac*OutletV.h;

	"Thermal Equilibrium"
	OutletV.T = OutletL.T;
	
	"Mechanical Equilibrium"
	OutletV.P = OutletL.P;
	
	"vaporization fraction "
	OutletV.v = 1.0;
	"vaporization fraction "
	OutletL.v = 0.0;
	
	# regularization functions
	zero_one = (1 + tanh(B * vsat))/2;
	one_zero = (1 - tanh(B * (vsat - 1)))/2;
	
	vfrac = zero_one * one_zero * vsat + 1 - one_zero;
	OutletL.z = zero_one*one_zero*xsat + (1-zero_one*one_zero)*Inlet.z;
	OutletV.z = zero_one*one_zero*ysat + (1-zero_one*one_zero)*Inlet.z;
end