#*-------------------------------------------------------------------
* 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.
*
*----------------------------------------------------------------------
* Author: Rafael de P. Soares and Paula B. Staudt
* $Id: sample_flashPH.mso 289 2007-06-20 15:07:08Z rafael $
*--------------------------------------------------------------------*#

using "streams";

Model FlashPHSteady
	ATTRIBUTES
	Pallete 	= true;
	Icon 		= "Flash"; 
	Brief 		= "Model of a static PH flash.";
	Info 		= "
	This model shows how to model a pressure enthalpy flash
	directly with the EMSO modeling language.

	This model is for demonstration purposes only, the flashPH
	routine available on VRTherm is much more robust.

	Assumptions:
	 * perfect mixing of both phases;

	Specify: 
	 * the feed stream;
	 * the heat duty;
	 * the outlet pressure.
	";	
	
	PARAMETERS
outer PP as Plugin(Brief = "External Physical Properties", Type="PP");
outer NComp as Integer;
	B as Real(Default=1000, Brief="Regularization Factor");

	VARIABLES
in	Inlet as stream (Brief="Feed Stream");
out	OutletL as liquid_stream (Brief="Liquid outlet stream");
out	OutletV as vapour_stream (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=-0.1, Upper=1.1, Brief="Vaporization fraction if saturated");
	Tsat as temperature(Lower=173, Upper=1473, Brief="Temperature if saturated");
	xsat(NComp) as Real(Lower=0, Upper=1, Brief="Liquid composition if saturated");
	ysat(NComp) as Real(Lower=0, Upper=1, 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;
	
	# 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

FlowSheet FlashPHTest
	PARAMETERS
	PP	as Plugin(Brief="Physical Properties",
		Type="PP",
		Components = ["methane", "isobutene", "benzene"],
		LiquidModel = "IdealLiquid",
		VapourModel = "Ideal"
	);
	NComp	as Integer;

	VARIABLES
	Q	as heat_rate (Brief="Heat supplied");
	
	SET
	NComp = PP.NumberOfComponents;
	
	DEVICES
	fl as FlashPHSteady;
	s1 as source;
	
	CONNECTIONS
	s1.Outlet to fl.Inlet;
	Q to fl.Q;
	
	SPECIFY
	s1.Outlet.F = 496.3 * 'kmol/h';
	s1.Outlet.T = 338 * 'K';
	s1.Outlet.P = 507.1 * 'kPa';
	s1.Outlet.z = [0.1, 0.7,0.2];
	
	fl.OutletL.P = 2.5 * 'atm';

	Q = 0 * 'kJ/h';
	#fl.OutletL.T = 315.06 * 'K';
	
	OPTIONS
	Dynamic = false;
end